R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- c(590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516,528,533,536,537,524,536,587,597,581,564,558,575,580,575,563,552,537,545,601,604,586,564,549,551,556,548,540,531,521,519,572,581,563,548,539) > par1 = '12' > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 18.27530 36.31871 25.85900 0.00000 > m$fitted level slope sea Jan 1 590.0000 0.0000000 0.0000000 Feb 1 583.2839 -2.9752149 -3.2838653 Mar 1 574.5932 -5.8085175 -0.5931573 Apr 1 571.2731 -4.4489715 1.7268831 May 1 571.6231 -1.8439185 1.3769400 Jun 1 605.7954 17.5296828 14.2046091 Jul 1 629.1533 20.6752954 -3.1533428 Aug 1 629.9655 9.9234477 -9.9655147 Sep 1 602.0409 -10.5759351 -14.0408945 Oct 1 570.2289 -22.0754965 -4.2288732 Nov 1 552.5475 -19.6966098 4.4524678 Dec 1 553.3208 -8.6158912 7.6792400 Jan 2 546.3188 -7.7565383 2.6812442 Feb 2 532.8342 -10.7952271 -0.8342000 Mar 2 523.3868 -10.1100426 2.6131994 Apr 2 512.6489 -10.4345520 -1.6489162 May 2 509.9793 -6.3580633 -10.9793127 Jun 2 533.5085 9.2546499 21.4914926 Jul 2 558.2252 17.3110459 6.7748476 Aug 2 550.5183 4.2464727 -8.5182876 Sep 2 537.6867 -4.6890040 -10.6866967 Oct 2 520.4847 -11.2375471 -10.4847086 Nov 2 513.8651 -8.8243020 0.1348613 Dec 2 508.3046 -7.1221049 8.6953831 Jan 3 501.6873 -6.8582463 6.3127270 Feb 3 492.5939 -8.0224080 0.4061364 Mar 3 484.9135 -7.8469521 5.0865349 Apr 3 476.1299 -8.3279389 -7.1299103 May 3 492.3764 4.4070642 -14.3763603 Jun 3 509.5612 11.0299996 18.4387870 Jul 3 519.1268 10.2732073 14.8731599 Aug 3 522.2294 6.5672607 -4.2294168 Sep 3 516.9429 0.4330654 -10.9428938 Oct 3 514.8179 -0.8909793 -12.8179376 Nov 3 516.1135 0.2398522 -0.1135043 Dec 3 518.1462 1.1671887 9.8537731 Jan 4 522.8375 2.9927381 10.1624525 Feb 4 531.5962 5.9709754 4.4038111 Mar 4 534.4448 4.3692999 2.5552304 Apr 4 540.8689 5.4227141 -16.8688707 May 4 553.3830 9.0756625 -17.3830270 Jun 4 566.9794 11.4079346 20.0205783 Jul 4 578.8594 11.6510137 18.1405954 Aug 4 583.3881 7.9876788 -2.3880578 Sep 4 578.7079 1.4710824 -14.7078936 Oct 4 573.3268 -2.0544544 -15.3268213 Nov 4 574.0048 -0.6484518 0.9952241 Dec 4 572.3550 -1.1641340 7.6450014 Jan 5 567.6619 -2.9827424 7.3381330 Feb 5 559.0436 -5.8802524 3.9563856 Mar 5 551.0751 -6.9497107 0.9249088 Apr 5 553.3367 -2.2350684 -16.3366795 May 5 561.4523 3.0777349 -16.4522920 Jun 5 576.8609 9.4170252 24.1390648 Jul 5 584.3784 8.4415276 19.6215609 Aug 5 585.8595 4.8725454 0.1404515 Sep 5 580.5295 -0.3555700 -16.5294692 Oct 5 569.6351 -5.7564983 -20.6351237 Nov 5 553.2900 -11.1864478 -2.2899683 Dec 5 544.1415 -10.1402674 11.8584542 Jan 6 536.3857 -8.9160027 11.6142796 Feb 6 532.7381 -6.2160829 7.2618985 Mar 6 532.1460 -3.3414689 -1.1459806 Apr 6 538.9386 1.8362646 -17.9386494 May 6 541.7611 2.3410781 -22.7611017 Jun 6 546.6605 3.6524627 25.3395101 Jul 6 556.3949 6.7681243 24.6051183 Aug 6 560.2574 5.2819097 2.7426453 Sep 6 560.9741 2.9492007 -12.9740798 Oct 6 557.4625 -0.3521959 -18.4624665 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.88888321 -0.52578661 0.22617693 0.42151059 3.19282578 2 0.14947402 -0.49812066 0.11282278 -0.05485493 0.67560841 2.57081723 3 0.04419121 -0.19253069 0.02899038 -0.08042681 2.12098604 1.09516323 4 0.30374512 0.49301682 -0.26519620 0.17543235 0.60827612 0.38666547 5 -0.30199503 -0.47988146 -0.17722783 0.78398087 0.88404400 1.05214565 6 0.20315099 0.44730075 0.47656618 0.86040911 0.08395282 0.21773935 Jul Aug Sep Oct Nov Dec 1 0.52206908 -1.78431483 -3.40139455 -1.90816122 0.39473845 1.83866398 2 1.33404670 -2.16823512 -1.48272449 -1.08661955 0.40050656 0.28338598 3 -0.12521611 -0.61456853 -1.01790677 -0.21973502 0.18781343 0.15433209 4 0.04024163 -0.60715486 -1.08118054 -0.58529288 0.23362062 -0.08574705 5 -0.16160728 -0.59140682 -0.86726543 -0.89684514 -0.90237256 0.17385412 6 0.51642286 -0.24628559 -0.38692673 -0.54824867 > 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 > mylagmax <- nx/2 > postscript(file="/var/wessaorg/rcomp/tmp/1ii001324657149.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/21o5e1324657149.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3mtf41324657149.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/44tmb1324657149.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',type='b') > grid() > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/5q7ak1324657149.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() null device 1 > > #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/6395a1324657149.tab") > > try(system("convert tmp/1ii001324657149.ps tmp/1ii001324657149.png",intern=TRUE)) character(0) > try(system("convert tmp/21o5e1324657149.ps tmp/21o5e1324657149.png",intern=TRUE)) character(0) > try(system("convert tmp/3mtf41324657149.ps tmp/3mtf41324657149.png",intern=TRUE)) character(0) > try(system("convert tmp/44tmb1324657149.ps tmp/44tmb1324657149.png",intern=TRUE)) character(0) > try(system("convert tmp/5q7ak1324657149.ps tmp/5q7ak1324657149.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.672 0.308 1.986