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(539,548,563,581,572,519,521,531,540,548,556,551,549,564,586,604,601,545,537,552,563,575,580,575,558,564,581,597,587,536,524,537,536,533,528,516,502,506,518,534,528,478,469,490,493,508,517,514,510,527,542,565,555,499,511,526,532,549,561,557,566,588,620,626,620,573,573,574,580,590) > 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 120.84853 0.00000 24.57873 0.00000 > m$fitted level slope sea Jan 1 539.0000 0.0000000 0.0000000 Feb 1 546.7535 0.3971431 1.2464552 Mar 1 558.5657 0.7819797 4.4342749 Apr 1 575.2757 0.9770008 5.7242755 May 1 574.1638 0.9641524 -2.1637519 Jun 1 533.8626 0.7160273 -14.8626224 Jul 1 519.9485 0.6200354 1.0514633 Aug 1 526.0772 0.6572274 4.9228333 Sep 1 536.3763 0.7221568 3.6236782 Oct 1 545.4307 0.7777761 2.5692583 Nov 1 553.5653 0.8264980 2.4347129 Dec 1 552.1355 0.8116574 -1.1355208 Jan 2 554.5425 0.7546352 -5.5425440 Feb 2 563.8763 0.8497917 0.1236596 Mar 2 581.1427 1.1805743 4.8572548 Apr 2 594.2691 1.3348835 9.7309004 May 2 593.8788 1.3235656 7.1211624 Jun 2 568.0874 1.2064934 -23.0873875 Jul 2 544.6071 1.1040635 -7.6071020 Aug 2 545.1718 1.1016509 6.8281966 Sep 2 555.9141 1.1480788 7.0858770 Oct 2 569.2007 1.2090100 5.7992725 Nov 2 575.5936 1.2299838 4.4064086 Dec 2 577.3885 1.2294632 -2.3885088 Jan 3 572.4349 1.2914525 -14.4349364 Feb 3 571.6161 1.2825138 -7.6161299 Mar 3 578.2264 1.3502141 2.7736194 Apr 3 584.3223 1.4041643 12.6777467 May 3 575.5454 1.3301534 11.4546278 Jun 3 558.5056 1.2446902 -22.5055929 Jul 3 539.1737 1.1660725 -15.1736883 Aug 3 533.7816 1.1404305 3.2184012 Sep 3 532.7480 1.1313348 3.2520435 Oct 3 530.2165 1.1168871 2.7835071 Nov 3 525.6720 1.1042873 2.3280137 Dec 3 519.1086 1.1146258 -3.1086165 Jan 4 516.6389 1.1270030 -14.6389237 Feb 4 516.1546 1.1221629 -10.1546247 Mar 4 516.7609 1.1176185 1.2391194 Apr 4 517.5512 1.1145237 16.4487714 May 4 512.7024 1.0714647 15.2976451 Jun 4 500.0386 1.0022469 -22.0386233 Jul 4 487.5478 0.9489178 -18.5478133 Aug 4 485.8128 0.9390092 4.1872410 Sep 4 487.4687 0.9415949 5.5313381 Oct 4 497.9671 0.9702626 10.0329367 Nov 4 508.2425 0.9831217 8.7575364 Dec 4 514.4191 0.9793687 -0.4191341 Jan 5 521.7526 0.9734284 -11.7526253 Feb 5 533.0532 1.0022250 -6.0532161 Mar 5 539.4145 1.0382023 2.5855449 Apr 5 544.6107 1.0710015 20.3892945 May 5 538.5027 1.0221324 16.4972946 Jun 5 523.6199 0.9394355 -24.6198642 Jul 5 525.8547 0.9447373 -14.8547065 Aug 5 523.7558 0.9339343 2.2442168 Sep 5 527.7444 0.9435196 4.2555623 Oct 5 537.6607 0.9647216 11.3393162 Nov 5 549.5819 0.9763933 11.4180999 Dec 5 557.9014 0.9755821 -0.9014186 Jan 6 574.7944 0.9795077 -8.7943957 Feb 6 590.8702 1.0208483 -2.8702345 Mar 6 610.9428 1.1257302 9.0572234 Apr 6 607.1355 1.0925999 18.8645451 May 6 601.4209 1.0498327 18.5791121 Jun 6 598.8644 1.0312448 -25.8644234 Jul 6 590.8581 0.9938875 -17.8580913 Aug 6 578.6126 0.9485756 -4.6126095 Sep 6 578.0486 0.9443391 1.9513714 Oct 6 580.9214 0.9481691 9.0785985 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.44712859 0.97254615 1.45208178 -0.19057292 -3.76175760 2 0.16345290 0.73504927 1.40649244 1.07800202 -0.15721646 -2.47089418 3 -0.58080113 -0.18859073 0.46625968 0.42600684 -0.92579551 -1.67364002 4 -0.32956831 -0.14505191 -0.04575143 -0.02936139 -0.54109306 -1.25068923 5 0.57988971 0.93148971 0.47874536 0.37348851 -0.65048508 -1.44725766 6 1.44805094 1.36297857 1.70923227 -0.44378596 -0.61638042 -0.32793774 Jul Aug Sep Oct Nov Dec 1 -1.33324719 0.50196493 0.87861499 0.75927347 0.67037643 -0.20559690 2 -2.24851901 -0.04911249 0.87788245 1.10570275 0.47203614 0.05145944 3 -1.87434777 -0.59727359 -0.19800771 -0.33357059 -0.51512235 -0.69913781 4 -1.22921883 -0.24450030 0.06529815 0.87001005 0.84648694 0.47327005 5 0.11800633 -0.27729359 0.27821447 0.81660733 0.99664529 0.66855804 6 -0.82311983 -1.20608528 -0.13773462 0.17547330 > 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/1hn7l1324645736.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/2ai1a1324645736.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/3m9961324645736.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/4e9zi1324645736.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/5zlet1324645736.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/68hg71324645736.tab") > > try(system("convert tmp/1hn7l1324645736.ps tmp/1hn7l1324645736.png",intern=TRUE)) character(0) > try(system("convert tmp/2ai1a1324645736.ps tmp/2ai1a1324645736.png",intern=TRUE)) character(0) > try(system("convert tmp/3m9961324645736.ps tmp/3m9961324645736.png",intern=TRUE)) character(0) > try(system("convert tmp/4e9zi1324645736.ps tmp/4e9zi1324645736.png",intern=TRUE)) character(0) > try(system("convert tmp/5zlet1324645736.ps tmp/5zlet1324645736.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.646 0.306 1.971