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(274412,272433,268361,268586,264768,269974,304744,309365,308347,298427,289231,291975,294912,293488,290555,284736,281818,287854,316263,325412,326011,328282,317480,317539,313737,312276,309391,302950,300316,304035,333476,337698,335932,323931,313927,314485,313218,309664,302963,298989,298423,301631,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) > 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 26862702 5572579 0 > m$fitted level slope sea Jan 1 274412.0 0.0000 0.000000 Feb 1 272915.3 -1536.5308 -482.331243 Mar 1 268799.1 -3330.9233 -438.122134 Apr 1 267702.2 -1681.8796 883.775392 May 1 265200.8 -2274.2944 -432.843149 Jun 1 268402.9 1663.7241 1571.104955 Jul 1 297595.2 21522.7182 7148.831830 Aug 1 313331.8 17346.1506 -3966.834411 Sep 1 312843.3 4473.9524 -4496.304041 Oct 1 301377.7 -7028.2588 -2950.693011 Nov 1 289315.8 -10660.3342 -84.814876 Dec 1 288777.3 -3356.6839 3197.719565 Jan 2 293436.0 2412.3431 1476.012590 Feb 2 294349.8 1329.7439 -861.808137 Mar 2 291817.3 -1360.2616 -1262.282139 Apr 2 283538.8 -6168.4666 1197.206358 May 2 281539.2 -3253.0152 278.838744 Jun 2 291768.2 6101.2716 -3914.189043 Jul 2 307372.7 12693.2656 8890.254727 Aug 2 324635.2 15867.7251 776.837684 Sep 2 329293.5 8077.0912 -3282.453699 Oct 2 330488.2 3293.9928 -2206.168014 Nov 2 323033.8 -4173.1577 -5553.844223 Dec 2 316907.9 -5528.4457 631.133940 Jan 3 311841.3 -5207.8368 1895.661436 Feb 3 310242.2 -2699.9613 2033.823076 Mar 3 307433.8 -2774.5482 1957.228521 Apr 3 302579.6 -4200.9939 370.417624 May 3 302623.1 -1275.1884 -2307.141364 Jun 3 309533.6 4350.9579 -5498.626230 Jul 3 324302.1 11496.3892 9173.935552 Aug 3 335079.8 11002.9818 2618.174078 Sep 3 339343.4 6372.3739 -3411.367580 Oct 3 327750.3 -5971.7510 -3819.312711 Nov 3 318827.3 -7998.5715 -4900.285108 Dec 3 313287.5 -6311.0202 1197.473828 Jan 4 311163.1 -3435.5949 2054.942265 Feb 4 307129.7 -3846.0021 2534.278842 Mar 4 300096.5 -6023.4304 2866.519292 Apr 4 297944.4 -3382.5409 1044.604557 May 4 301737.4 1526.5434 -3314.421643 Jun 4 309852.9 6031.3261 -8221.917814 Jul 4 320537.7 9207.0652 9227.313424 Aug 4 330547.0 9754.6592 4536.015312 Sep 4 328591.4 1755.9973 -975.391099 Oct 4 315616.7 -8306.8081 -6497.708272 Nov 4 302611.8 -11515.0378 -6695.768810 Dec 4 291415.2 -11297.5892 -2.209161 Jan 5 286805.1 -6728.7865 4736.901515 Feb 5 280371.3 -6527.4181 4306.734601 Mar 5 273925.6 -6471.7823 2549.420408 Apr 5 272134.8 -3288.0907 431.249115 May 5 270572.3 -2112.3122 -5591.314117 Jun 5 273017.5 993.2397 -9727.474822 Jul 5 285819.2 9029.5910 10986.836898 Aug 5 295852.4 9712.4350 7745.613354 Sep 5 288017.9 -2230.7695 -1023.905979 Oct 5 281589.0 -5088.5750 -5162.037516 Nov 5 274020.9 -6776.2525 -7596.912884 Dec 5 268704.9 -5782.2181 -1551.926201 Jan 6 263247.6 -5560.9839 5133.419510 Feb 6 257671.8 -5571.0209 4850.153154 Mar 6 253324.8 -4739.7566 2217.239328 Apr 6 251823.0 -2541.7159 1334.978358 May 6 251048.1 -1340.9257 -7245.122719 Jun 6 262013.2 7025.4766 -11272.223162 Jul 6 270261.7 7856.2139 10183.347173 Aug 6 273298.8 4584.2790 11958.230593 Sep 6 272059.5 629.9607 -1083.522663 Oct 6 266969.3 -3254.6432 -5893.333068 Nov 6 263581.1 -3345.3753 -7978.089840 Dec 6 261541.6 -2458.3175 -1165.563047 Jan 7 258529.4 -2834.6371 5373.619493 Feb 7 258408.7 -991.9839 5882.336557 Mar 7 260583.6 1155.2963 2692.406892 Apr 7 261992.6 1327.2302 579.447714 May 7 266871.5 3736.6587 -10704.473009 Jun 7 275268.6 6899.5264 -11047.620019 Jul 7 282527.0 7142.9190 11333.008420 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.338623655 -0.366553965 0.327333274 -0.113348645 0.760817481 2 1.116438853 -0.212154990 -0.516326195 -0.942283389 0.562316102 1.800731415 3 0.062044593 0.484576269 -0.014357627 -0.276731280 0.565653289 1.083446999 4 0.555785713 -0.079174399 -0.419638236 0.510806463 0.948859165 0.868319498 5 0.882345692 0.038836932 0.010727518 0.615117380 0.227162823 0.598958672 6 0.042705311 -0.001935726 0.160319978 0.424470423 0.231922252 1.614095158 7 -0.072622797 0.355379864 0.414183125 0.033194225 0.465263807 0.610293303 Jul Aug Sep Oct Nov Dec 1 3.833383059 -0.805779869 -2.483633933 -2.219257592 -0.700776746 1.409176054 2 1.272830353 0.612594340 -1.503034406 -0.922919277 -1.440735346 -0.261574901 3 1.378447497 -0.095228514 -0.893408863 -2.381771871 -0.391063486 0.325809075 4 0.612377862 0.105666368 -1.543322409 -1.941537391 -0.619061200 0.041982445 5 1.549594892 0.131737935 -2.304438598 -0.551414776 -0.325687246 0.191897923 6 0.160200632 -0.631159447 -0.762968725 -0.749588610 -0.017510878 0.171228690 7 0.046941292 > 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/1ski51259939888.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/2lhio1259939888.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/3vhjg1259939888.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/48ed71259939888.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/5dn1x1259939888.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/6nfpa1259939888.tab") > system("convert tmp/1ski51259939888.ps tmp/1ski51259939888.png") > system("convert tmp/2lhio1259939888.ps tmp/2lhio1259939888.png") > system("convert tmp/3vhjg1259939888.ps tmp/3vhjg1259939888.png") > system("convert tmp/48ed71259939888.ps tmp/48ed71259939888.png") > system("convert tmp/5dn1x1259939888.ps tmp/5dn1x1259939888.png") > > > proc.time() user system elapsed 1.338 0.817 1.542