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(16111,15554,15220,14807,14291,14653,17006,18032,16558,16102,15055,15484,14596,14609,13923,14226,14056,14278,16142,16509,15680,14086,13129,13086,13096,12280,11534,11135,10903,10926,13220,13581,11788,11088,10434,11061,10828,10270,10360,9899,9395,9944,12117,12474,11106,10643,10227,11273,11516,11583,11605,11414,11181,12000,14007,14582,13251,12806,12645,13869,13342,13079,12513,12331,11882,12388,14394,14635,13218,12554,12031,12429) > 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 36658.75 133808.76 54692.67 0.00 > m$fitted level slope sea Jan 1 16111.000 0.000000 0.000000 Feb 1 15707.158 -238.373199 -153.157870 Mar 1 15218.651 -383.884415 1.348696 Apr 1 14798.938 -405.991721 8.061998 May 1 14319.564 -450.653491 -28.563558 Jun 1 14452.308 -96.978687 200.691556 Jul 1 16376.496 1133.226648 629.503631 Aug 1 18073.612 1476.957460 -41.611589 Sep 1 17389.624 159.784953 -831.623562 Oct 1 16306.209 -597.775326 -204.209108 Nov 1 15129.998 -950.203300 -74.998262 Dec 1 15107.357 -385.063817 376.642693 Jan 2 14631.300 -439.776966 -35.299988 Feb 2 14538.314 -233.792240 70.686132 Mar 2 13964.307 -428.502934 -41.306802 Apr 2 13881.568 -226.249234 344.431553 May 2 14121.365 48.604352 -65.365467 Jun 2 14718.030 369.803651 -440.030156 Jul 2 15652.426 700.682982 489.573725 Aug 2 16062.112 529.725556 446.888352 Sep 2 16202.033 300.478961 -522.032979 Oct 2 14729.778 -741.969901 -643.777543 Nov 2 13516.138 -1018.960611 -387.137571 Dec 2 12600.346 -958.494814 485.654218 Jan 3 12752.819 -305.873286 343.181433 Feb 3 12183.448 -459.826903 96.551853 Mar 3 11633.596 -511.649988 -99.595747 Apr 3 10983.817 -591.535875 151.183013 May 3 10996.579 -239.444085 -93.579485 Jun 3 11459.010 168.544810 -533.010326 Jul 3 12495.516 672.103419 724.484157 Aug 3 13090.213 627.136947 490.787206 Sep 3 12210.147 -249.519361 -422.147050 Oct 3 11573.908 -474.426694 -485.907991 Nov 3 10946.244 -563.444472 -512.243551 Dec 3 10788.181 -327.996115 272.819150 Jan 4 10391.373 -368.022010 436.627073 Feb 4 10025.512 -366.769541 244.488170 Mar 4 10161.212 -77.562177 198.787606 Apr 4 10013.482 -117.998198 -114.481565 May 4 9842.806 -148.497620 -447.805637 Jun 4 10568.474 357.519265 -624.473897 Jul 4 11288.107 566.738868 828.893496 Aug 4 11566.196 399.958416 907.804080 Sep 4 11515.584 139.421981 -409.583617 Oct 4 11221.001 -111.538778 -578.001495 Nov 4 10914.722 -224.102576 -687.722125 Dec 4 10874.729 -117.628299 398.271326 Jan 5 11001.521 23.812024 514.478650 Feb 5 11359.547 216.708587 223.453295 Mar 5 11399.889 115.315061 205.110536 Apr 5 11519.588 117.836220 -105.588048 May 5 11902.315 270.646729 -721.315293 Jun 5 12643.572 542.275914 -643.571553 Jul 5 13138.705 515.104288 868.294983 Aug 5 13535.550 446.990197 1046.450424 Sep 5 13607.395 230.883637 -356.394960 Oct 5 13426.404 -6.412505 -620.403865 Nov 5 13385.646 -26.204168 -740.645584 Dec 5 13514.596 63.266141 354.403702 Jan 6 13098.846 -212.996686 243.153974 Feb 6 12763.637 -283.337262 315.363123 Mar 6 12340.825 -363.410741 172.174886 Apr 6 12399.460 -121.063401 -68.459563 May 6 12679.098 109.507113 -797.097546 Jun 6 12990.744 225.906344 -602.743941 Jul 6 13424.623 345.542883 969.377227 Aug 6 13541.516 214.146641 1093.484169 Sep 6 13515.093 75.939872 -297.092898 Oct 6 13295.210 -94.039991 -741.210332 Nov 6 12906.785 -263.283787 -875.784717 Dec 6 12071.405 -592.418610 357.594663 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.986921521 -0.423177336 -0.060115583 -0.120062173 0.965095210 2 -0.155811572 0.555417887 -0.532763421 0.560948790 0.747557502 0.873590714 3 1.799225476 -0.419185053 -0.141447310 -0.220080146 0.963125443 1.111207934 4 -0.109704268 0.003415742 0.789807501 -0.110990360 -0.083527113 1.380703462 5 0.386963452 0.526427258 -0.277007849 0.006908948 0.418493496 0.741989084 6 -0.755284975 -0.192036319 -0.218807543 0.663611388 0.631319806 0.318144827 Jul Aug Sep Oct Nov Dec 1 3.364632199 0.939643055 -3.600754316 -2.070984730 -0.963447022 1.544945666 2 0.904340365 -0.467457509 -0.626665144 -2.849919360 -0.757261137 0.165760647 3 1.374773752 -0.122933432 -2.396618368 -0.614844357 -0.243479849 0.645335125 4 0.571044856 -0.455783470 -0.712272523 -0.686195345 -0.307984460 0.291600706 5 -0.074176620 -0.186091585 -0.590786033 -0.648982760 -0.054159081 0.244887437 6 0.326688808 -0.358932873 -0.377807990 -0.464951582 -0.463142900 -0.900535266 > 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/wessaorg/rcomp/tmp/128y91322508023.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() null device 1 > mylagmax <- nx/2 > postscript(file="/var/wessaorg/rcomp/tmp/22zpj1322508023.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/356bq1322508023.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/47wab1322508023.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/5y9m11322508023.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/6xvbj1322508023.tab") > > try(system("convert tmp/128y91322508023.ps tmp/128y91322508023.png",intern=TRUE)) character(0) > try(system("convert tmp/22zpj1322508023.ps tmp/22zpj1322508023.png",intern=TRUE)) character(0) > try(system("convert tmp/356bq1322508023.ps tmp/356bq1322508023.png",intern=TRUE)) character(0) > try(system("convert tmp/47wab1322508023.ps tmp/47wab1322508023.png",intern=TRUE)) character(0) > try(system("convert tmp/5y9m11322508023.ps tmp/5y9m11322508023.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.748 0.368 2.138