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(267413,267366,264777,258863,254844,254868,277267,285351,286602,283042,276687,277915,277128,277103,275037,270150,267140,264993,287259,291186,292300,288186,281477,282656,280190,280408,276836,275216,274352,271311,289802,290726,292300,278506,269826,265861,269034,264176,255198,253353,246057,235372,258556,260993,254663,250643,243422,247105,248541,245039,237080,237085,225554,226839,247934,248333,246969,245098,246263,255765,264319,268347,273046,273963,267430,271993,292710,295881,293299) > 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 0 9995583 3972943 0 > m$fitted level slope sea Jan 1 267413.0 0.00000 0.00000 Feb 1 267384.4 -30.15951 -18.35765 Mar 1 265475.8 -1230.64730 -698.76307 Apr 1 260376.2 -3814.08699 -1513.17386 May 1 254967.3 -4876.39017 -123.29687 Jun 1 253450.9 -2679.70951 1417.08451 Jul 1 270271.0 10102.88450 6996.00853 Aug 1 285788.0 13660.77704 -437.00352 Sep 1 290736.3 7934.00113 -4134.27637 Oct 1 286687.5 59.23771 -3645.50077 Nov 1 278491.4 -5364.70454 -1804.42028 Dec 1 275926.5 -3525.34751 1988.51161 Jan 2 275888.6 -1242.35391 1239.38108 Feb 2 276119.6 -273.44201 983.35995 Mar 2 274310.5 -1256.50466 726.45338 Apr 2 270425.8 -2914.36198 -275.80322 May 2 267870.7 -2684.87711 -730.66080 Jun 2 269620.8 131.33077 -4627.75893 Jul 2 279812.7 6487.89198 7446.27653 Aug 2 289067.1 8238.78671 2118.91402 Sep 2 293757.2 5989.62412 -1457.16536 Oct 2 291683.1 878.09063 -3497.11604 Nov 2 286133.3 -3194.37723 -4656.28454 Dec 2 281944.9 -3823.35446 711.06035 Jan 3 279026.6 -3250.51792 1163.43005 Feb 3 277630.5 -2075.11210 2777.53838 Mar 3 274670.0 -2631.76033 2165.98942 Apr 3 273993.9 -1409.30215 1222.13894 May 3 275843.9 638.32005 -1491.89232 Jun 3 279467.6 2514.41173 -8156.55359 Jul 3 283709.1 3595.88142 6092.85372 Aug 3 287633.9 3801.70085 3092.12341 Sep 3 290947.5 3495.88292 1352.47031 Oct 3 283703.8 -3236.44053 -5197.84538 Nov 3 275263.7 -6497.21045 -5437.65478 Dec 3 266174.6 -8120.73640 -313.58248 Jan 4 265211.7 -3633.50684 3822.26704 Feb 4 261286.6 -3816.25619 2889.35050 Mar 4 254427.3 -5715.73511 770.68141 Apr 4 251794.8 -3796.22282 1558.15703 May 4 248859.1 -3259.28167 -2802.10509 Jun 4 245091.6 -3576.77771 -9719.62502 Jul 4 250185.4 1830.32170 8370.56668 Aug 4 256661.5 4724.16035 4331.52624 Sep 4 252630.0 -731.56120 2032.96051 Oct 4 252529.1 -338.48683 -1886.07339 Nov 4 249524.3 -2000.44578 -6102.25902 Dec 4 249320.3 -880.46467 -2215.30510 Jan 5 245548.3 -2684.14785 2992.67449 Feb 5 241221.4 -3708.24407 3817.59932 Mar 5 237068.3 -3984.86210 11.69206 Apr 5 234614.6 -3034.09109 2470.37837 May 5 229476.3 -4342.42688 -3922.27453 Jun 5 236099.6 2483.23188 -9260.63748 Jul 5 240991.2 3980.98270 6942.81438 Aug 5 242009.1 2140.45409 6323.85962 Sep 5 244702.5 2483.81189 2266.47770 Oct 5 246686.4 2173.28268 -1588.35775 Nov 5 251949.6 4092.86270 -5686.60854 Dec 5 256802.8 4565.41443 -1037.82634 Jan 6 260744.5 4177.64728 3574.48715 Feb 6 264204.2 3731.62971 4142.79504 Mar 6 271478.9 5928.85175 1567.10035 Apr 6 273042.4 3224.03589 920.64154 May 6 275607.5 2815.42548 -8177.50102 Jun 6 281017.3 4425.81673 -9024.30066 Jul 6 284951.7 4120.90337 7758.31555 Aug 6 289357.1 4297.23225 6523.91262 Sep 6 291689.9 3080.48100 1609.08628 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.01205733 -0.38720778 -0.87293996 -0.33292301 0.69325099 2 0.72489753 0.31150299 -0.30821847 -0.53189606 0.07309066 0.88695114 3 0.18176994 0.37241994 -0.17552705 0.38805948 0.65079018 0.59272695 4 1.42162440 -0.05779934 -0.59991217 0.60805545 0.17031602 -0.10043587 5 -0.57092977 -0.32377375 -0.08741744 0.30098773 -0.41459229 2.16002895 6 -0.12268889 -0.14099950 0.69456626 -0.85604334 -0.12941736 0.50966443 Jul Aug Sep Oct Nov Dec 1 4.05104096 1.12545652 -1.81127216 -2.49085227 -1.71558561 0.58178346 2 2.00973341 0.55428377 -0.71139881 -1.61678483 -1.28817390 -0.19904542 3 0.34155314 0.06511052 -0.09674531 -2.12942824 -1.03149128 -0.51393788 4 1.70823975 0.91494108 -1.72577205 0.12434602 -0.52584855 0.35452496 5 0.47338302 -0.58180650 0.10859677 -0.09824597 0.60747794 0.14957021 6 -0.09639971 0.05573808 -0.38479430 > 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/1v3491259945854.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/2utzm1259945854.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/3vgfx1259945854.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/4945l1259945854.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/5neg91259945854.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/61rew1259945854.tab") > system("convert tmp/1v3491259945854.ps tmp/1v3491259945854.png") > system("convert tmp/2utzm1259945854.ps tmp/2utzm1259945854.png") > system("convert tmp/3vgfx1259945854.ps tmp/3vgfx1259945854.png") > system("convert tmp/4945l1259945854.ps tmp/4945l1259945854.png") > system("convert tmp/5neg91259945854.ps tmp/5neg91259945854.png") > > > proc.time() user system elapsed 1.305 0.812 1.523