R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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(325412,326011,328282,317480,317539,313737,312276,309391,302950,300316,304035,333476,337698,335932,323931,313927,314485,313218,309664,302963,298989,298423,310631,329765,335083,327616,309119,295916,291413,291542,284678,276475,272566,264981,263290,296806,303598,286994,276427,266424,267153,268381,262522,255542,253158,243803,250741,280445,285257,270976,261076,255603,260376,263903,264291,263276,262572,256167,264221,293860,300713,287224) > par1 = '12' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 44693630 0 6390054 0 > m$fitted level slope sea Jan 1 325412.0 0.000000 0.00000 Feb 1 325942.6 27.752756 68.44916 Mar 1 327718.7 76.525533 563.27905 Apr 1 319986.3 7.480058 -2506.25852 May 1 317443.8 -5.787950 95.17930 Jun 1 314384.5 -22.929012 -647.46083 Jul 1 312519.7 -34.012704 -243.65743 Aug 1 309980.4 -49.246260 -589.39188 Sep 1 304397.9 -82.656697 -1447.93353 Oct 1 300840.6 -103.476526 -524.56271 Nov 1 303054.5 -89.680041 980.52130 Dec 1 326654.1 50.518956 6821.87692 Jan 2 335235.9 -266.862226 2462.08079 Feb 2 336715.2 -242.484768 -783.20007 Mar 2 326237.3 -433.235919 -2306.26316 Apr 2 318390.4 -507.018859 -4463.44202 May 2 314626.3 -522.396858 -141.29110 Jun 2 313748.8 -523.648132 -530.76125 Jul 2 310805.4 -532.588715 -1141.37644 Aug 2 304558.9 -555.255230 -1595.93356 Sep 2 300238.0 -570.987029 -1249.04143 Oct 2 299128.5 -573.350834 -705.51039 Nov 2 310613.6 -528.009733 17.43356 Dec 2 322025.2 -539.454475 7739.77707 Jan 3 329559.1 -631.820253 5523.88460 Feb 3 326666.6 -644.617825 949.41118 Mar 3 314118.8 -794.790094 -4999.81938 Apr 3 303114.8 -892.576432 -7198.82519 May 3 294177.5 -936.380996 -2764.45632 Jun 3 291407.6 -942.811519 134.41098 Jul 3 285795.1 -957.554687 -1117.12526 Aug 3 278771.4 -978.141342 -2296.40524 Sep 3 274508.1 -990.302103 -1942.13528 Oct 3 270422.7 -1001.498680 -5441.71378 Nov 3 267626.1 -1005.314073 -4336.08902 Dec 3 283793.1 -1034.228259 13012.92171 Jan 4 293129.3 -1080.628902 10468.67903 Feb 4 285471.1 -1103.640611 1522.87256 Mar 4 279858.1 -1144.262771 -3431.11244 Apr 4 273542.6 -1188.010107 -7118.62938 May 4 270130.2 -1200.759316 -2977.17666 Jun 4 267359.7 -1206.726870 1021.30096 Jul 4 262992.9 -1216.596316 -470.91954 Aug 4 257929.4 -1228.704377 -2387.39726 Sep 4 254454.9 -1236.030223 -1296.89960 Oct 4 250362.4 -1244.142800 -6559.36467 Nov 4 257031.4 -1234.213338 -6290.44044 Dec 4 266336.4 -1246.878892 14108.57793 Jan 5 270928.1 -1256.800040 14328.93922 Feb 5 268845.4 -1259.163113 2130.62944 Mar 5 264450.4 -1280.718863 -3374.42225 Apr 5 262408.4 -1286.285696 -6805.42362 May 5 262351.8 -1279.318115 -1975.83729 Jun 5 261637.4 -1277.035798 2265.64622 Jul 5 262485.1 -1270.196656 1805.93009 Aug 5 263710.7 -1262.717451 -434.67299 Sep 5 263126.1 -1260.784216 -554.06014 Oct 5 264490.0 -1254.916550 -8323.02117 Nov 5 270948.1 -1248.204199 -6727.08339 Dec 5 278477.4 -1253.606639 15382.61544 Jan 6 284007.2 -1256.087608 16705.82532 Feb 6 284283.5 -1252.056783 2940.48855 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.04952179 0.25136428 -1.17182664 -0.38238507 -0.45769140 2 1.45510684 0.24136407 -1.45975595 -1.10617166 -0.48818961 -0.05320207 3 1.25528892 -0.32964530 -1.71829631 -1.51571342 -1.20466977 -0.27467047 4 1.57351519 -0.97093567 -0.65780180 -0.76615534 -0.33262526 -0.23511781 5 0.87802710 -0.12237452 -0.46053668 -0.11276550 0.18364474 0.08457478 6 1.01617356 0.22743801 Jul Aug Sep Oct Nov Dec 1 -0.27603137 -0.37544378 -0.82922368 -0.52072741 0.34728514 3.55004923 2 -0.36236683 -0.85554488 -0.56387482 -0.08066881 1.80596069 1.78781620 3 -0.69931623 -0.90830551 -0.49195141 -0.46355703 -0.26863419 2.57470227 4 -0.47329944 -0.57612320 -0.33634029 -0.42767449 1.18382262 1.57987078 5 0.31826358 0.37383236 0.10155324 0.39285559 1.15377720 1.31483654 6 > 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/www/html/rcomp/tmp/1qq0d1259689828.ps",horizontal=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/www/html/rcomp/tmp/2ykzu1259689828.ps",horizontal=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/www/html/rcomp/tmp/3rm9b1259689828.ps",horizontal=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/www/html/rcomp/tmp/4hvwo1259689828.ps",horizontal=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/www/html/rcomp/tmp/5f5f61259689828.ps",horizontal=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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/www/html/rcomp/tmp/6hsof1259689828.tab") > system("convert tmp/1qq0d1259689828.ps tmp/1qq0d1259689828.png") > system("convert tmp/2ykzu1259689828.ps tmp/2ykzu1259689828.png") > system("convert tmp/3rm9b1259689828.ps tmp/3rm9b1259689828.png") > system("convert tmp/4hvwo1259689828.ps tmp/4hvwo1259689828.png") > system("convert tmp/5f5f61259689828.ps tmp/5f5f61259689828.png") > > > proc.time() user system elapsed 1.398 0.815 1.615