par1 <- as.numeric(par1)
nx <- length(x)
x <- ts(x,frequency=par1)
m <- StructTS(x,type='BSM')
m$coef
m$fitted
m$resid
mylevel <- as.numeric(m$fitted[,'level'])
myslope <- as.numeric(m$fitted[,'slope'])
myseas <- as.numeric(m$fitted[,'sea'])
myresid <- as.numeric(m$resid)
myfit <- mylevel+myseas
bitmap(file='test1.png')
op <- par(mfrow = c(2,1))
plot(as.numeric(x),main='Interpolation Fit',ylab='Observed (black) / Fitted (red)',xlab='time',type='b')
lines(myfit, col = 'red')
grid()
plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time')
grid()
par(op)
dev.off()
mylagmax <- nx/2
bitmap(file='test2.png')
op <- par(mfrow = c(3,2))
acf(as.numeric(x),lag.max = mylagmax,main='Observed')
acf(myfit,na.action=na.pass,lag.max = mylagmax,main='Fitted')
acf(mylevel,na.action=na.pass,lag.max = mylagmax,main='Level')
acf(myslope,na.action=na.pass,lag.max = mylagmax,main='Slope')
acf(myseas,na.action=na.pass,lag.max = mylagmax,main='Seasonal')
acf(myresid,na.action=na.pass,lag.max = mylagmax,main='Standardized Residals')
par(op)
dev.off()
bitmap(file='test3.png')
op <- par(mfrow = c(3,2))
spectrum(as.numeric(x),main='Observed')
spectrum(myfit,main='Fitted')
spectrum(mylevel,main='Level')
spectrum(myslope,main='Slope')
spectrum(myseas,main='Seasonal')
spectrum(myresid,main='Standardized Residals')
par(op)
dev.off()
bitmap(file='test4.png')
op <- par(mfrow = c(2,2))
cpgram(as.numeric(x),main='Observed')
cpgram(myfit,main='Fitted')
cpgram(myseas,main='Seasonal')
cpgram(myresid,main='Standardized Residals')
par(op)
dev.off()
bitmap(file='test5.png')
op <- par(mfrow = c(2,2))
hist(m$resid,main='Residual Histogram')
plot(density(m$resid),main='Residual Kernel Density')
qqnorm(m$resid,main='Residual Normal QQ Plot')
qqline(m$resid)
plot(m$resid^2, myfit^2,main='Sq.Resid vs. Sq.Fit',xlab='Squared residuals',ylab='Squared Fit')
par(op)
dev.off()
load(file='createtable')
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Structural Time Series Model',7,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'t',header=TRUE)
a<-table.element(a,'Observed',header=TRUE)
a<-table.element(a,'Fitted',header=TRUE)
a<-table.element(a,'Level',header=TRUE)
a<-table.element(a,'Slope',header=TRUE)
a<-table.element(a,'Seasonal',header=TRUE)
a<-table.element(a,'Stand. Residuals',header=TRUE)
a<-table.row.end(a)
for (i in 1:nx) {
a<-table.row.start(a)
a<-table.element(a,i,header=TRUE)
a<-table.element(a,x[i])
a<-table.element(a,myfit[i])
a<-table.element(a,mylevel[i])
a<-table.element(a,myslope[i])
a<-table.element(a,myseas[i])
a<-table.element(a,myresid[i])
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file='mytable.tab')
|
