x <- c(235.1 ,280.7 ,264.6 ,240.7 ,201.4 ,240.8 ,241.1 ,223.8 ,206.1 ,174.7 ,203.3 ,220.5 ,299.5 ,347.4 ,338.3 ,327.7 ,351.6 ,396.6 ,438.8 ,395.6 ,363.5 ,378.8 ,357 ,369 ,464.8 ,479.1 ,431.3 ,366.5 ,326.3 ,355.1 ,331.6 ,261.3 ,249 ,205.5 ,235.6 ,240.9 ,264.9 ,253.8 ,232.3 ,193.8 ,177 ,213.2 ,207.2 ,180.6 ,188.6 ,175.4 ,199 ,179.6 ,225.8 ,234 ,200.2 ,183.6 ,178.2 ,203.2 ,208.5 ,191.8 ,172.8 ,148 ,159.4 ,154.5 ,213.2 ,196.4 ,182.8 ,176.4 ,153.6 ,173.2 ,171 ,151.2 ,161.9 ,157.2 ,201.7 ,236.4 ,356.1 ,398.3 ,403.7 ,384.6 ,365.8 ,368.1 ,367.9 ,347 ,343.3 ,292.9 ,311.5 ,300.9 ,366.9 ,356.9 ,329.7 ,316.2 ,269 ,289.3 ,266.2 ,253.6 ,233.8 ,228.4 ,253.6 ,260.1 ,306.6 ,309.2 ,309.5 ,271 ,279.9 ,317.9 ,298.4 ,246.7 ,227.3 ,209.1 ,259.9 ,266 ,320.6 ,308.5 ,282.2 ,262.7 ,263.5 ,313.1 ,284.3 ,252.6 ,250.3 ,246.5 ,312.7 ,333.2 ,446.4 ,511.6 ,515.5 ,506.4 ,483.2 ,522.3 ,509.8 ,460.7 ,405.8 ,375 ,378.5 ,406.8 ,467.8 ,469.8 ,429.8 ,355.8 ,332.7 ,378 ,360.5 ,334.7 ,319.5 ,323.1 ,363.6 ,352.1 ,411.9 ,388.6 ,416.4 ,360.7 ,338 ,417.2 ,388.4 ,371.1 ,331.5 ,353.7 ,396.7 ,447 ,533.5 ,565.4 ,542.3 ,488.7 ,467.1 ,531.3 ,496.1 ,444 ,403.4 ,386.3 ,394.1 ,404.1 ,462.1 ,448.1 ,432.3 ,386.3 ,395.2 ,421.9 ,382.9 ,384.2 ,345.5 ,323.4 ,372.6 ,376 ,462.7 ,487 ,444.2 ,399.3 ,394.9 ,455.4 ,414 ,375.5 ,347 ,339.4 ,385.8 ,378.8 ,451.8 ,446.1 ,422.5 ,383.1 ,352.8 ,445.3 ,367.5 ,355.1 ,326.2 ,319.8 ,331.8 ,340.9 ,394.1 ,417.2 ,369.9 ,349.2 ,321.4 ,405.7 ,342.9 ,316.5 ,284.2 ,270.9 ,288.8 ,278.8 ,324.4 ,310.9 ,299 ,273 ,279.3 ,359.2 ,305 ,282.1 ,250.3 ,246.5 ,257.9 ,266.5 ,315.9 ,318.4 ,295.4 ,266.4 ,245.8 ,362.8 ,324.9 ,294.2 ,289.5 ,295.2 ,290.3 ,272 ,307.4 ,328.7 ,292.9 ,249.1 ,230.4 ,361.5 ,321.7 ,277.2 ,260.7 ,251 ,257.6 ,241.8 ,287.5 ,292.3 ,274.7 ,254.2 ,230 ,339 ,318.2 ,287 ,295.8 ,284 ,271 ,262.7 ,340.6 ,379.4 ,373.3 ,355.2 ,338.4 ,466.9 ,451 ,422 ,429.2 ,425.9 ,460.7 ,463.6 ,541.4 ,544.2 ,517.5 ,469.4 ,439.4 ,549 ,533 ,506.1 ,484 ,457 ,481.5 ,469.5 ,544.7 ,541.2 ,521.5 ,469.7 ,434.4 ,542.6 ,517.3 ,485.7 ,465.8 ,447 ,426.6 ,411.6 ,467.5 ,484.5 ,451.2 ,417.4 ,379.9 ,484.7 ,455 ,420.8 ,416.5 ,376.3 ,405.6 ,405.8 ,500.8 ,514 ,475.5 ,430.1 ,414.4 ,538 ,526 ,488.5 ,520.2 ,504.4 ,568.5 ,610.6 ,818 ,830.9 ,835.9 ,782 ,762.3 ,856.9 ,820.9 ,769.6 ,752.2 ,724.4 ,723.1 ,719.5 ,817.4 ,803.3 ,752.5 ,689 ,630.4 ,765.5 ,757.7 ,732.2 ,702.6 ,683.3 ,709.5 ,702.2 ,784.8 ,810.9 ,755.6 ,656.8 ,615.1 ,745.3 ,694.1 ,675.7 ,643.7 ,622.1 ,634.6 ,588 ,689.7 ,673.9 ,647.9 ,568.8 ,545.7 ,632.6 ,643.8 ,593.1 ,579.7 ,546 ,562.9 ,572.5) par1 = '12' 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 postscript(file="/var/fisher/rcomp/tmp/1kwkb1352555245.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time') grid() dev.off() mylagmax <- nx/2 postscript(file="/var/fisher/rcomp/tmp/2g9az1352555245.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) op <- par(mfrow = c(2,2)) acf(as.numeric(x),lag.max = mylagmax,main='Observed') acf(mylevel,na.action=na.pass,lag.max = mylagmax,main='Level') 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() postscript(file="/var/fisher/rcomp/tmp/3jgy11352555245.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) op <- par(mfrow = c(2,2)) spectrum(as.numeric(x),main='Observed') spectrum(mylevel,main='Level') spectrum(myseas,main='Seasonal') spectrum(myresid,main='Standardized Residals') par(op) dev.off() postscript(file="/var/fisher/rcomp/tmp/4h1an1352555245.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) op <- par(mfrow = c(2,2)) cpgram(as.numeric(x),main='Observed') cpgram(mylevel,main='Level') cpgram(myseas,main='Seasonal') cpgram(myresid,main='Standardized Residals') par(op) dev.off() postscript(file="/var/fisher/rcomp/tmp/52xro1352555245.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 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() #Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/fisher/rcomp/createtable") a<-table.start() a<-table.row.start(a) a<-table.element(a,'Structural Time Series Model',6,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,'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,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="/var/fisher/rcomp/tmp/6trwn1352555245.tab") try(system("convert tmp/1kwkb1352555245.ps tmp/1kwkb1352555245.png",intern=TRUE)) try(system("convert tmp/2g9az1352555245.ps tmp/2g9az1352555245.png",intern=TRUE)) try(system("convert tmp/3jgy11352555245.ps tmp/3jgy11352555245.png",intern=TRUE)) try(system("convert tmp/4h1an1352555245.ps tmp/4h1an1352555245.png",intern=TRUE)) try(system("convert tmp/52xro1352555245.ps tmp/52xro1352555245.png",intern=TRUE))