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(519,517,510,509,501,507,569,580,578,565,547,555,562,561,555,544,537,543,594,611,613,611,594,595,591,589,584,573,567,569,621,629,628,612,595,597,593,590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516) > 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.00000 79.79558 18.65712 0.00000 > m$fitted level slope sea Jan 1 519.0000 0.0000000 0.0000000 Feb 1 517.5349 -1.5093197 -0.5348858 Mar 1 511.1608 -4.8364044 -1.1608179 Apr 1 508.0371 -3.5929768 0.9628622 May 1 501.9312 -5.3826198 -0.9312109 Jun 1 504.4815 0.2252912 2.5185456 Jul 1 554.9516 35.8673319 14.0483577 Aug 1 586.0620 32.4901132 -6.0620277 Sep 1 586.9886 10.0833547 -8.9886245 Oct 1 570.1722 -9.0087067 -5.1722310 Nov 1 548.4150 -18.0560585 -1.4150167 Dec 1 548.5595 -5.1393779 6.4404879 Jan 2 558.8378 5.7727463 3.1621589 Feb 2 562.3085 4.1370203 -1.3084589 Mar 2 557.6871 -1.8726752 -2.6870507 Apr 2 542.4039 -11.0334912 1.5961136 May 2 536.6845 -7.3750244 0.3155427 Jun 2 551.2683 7.6204892 -8.2682558 Jul 2 577.4631 20.2926890 16.5368522 Aug 2 607.7603 27.1302151 3.2396539 Sep 2 619.0736 16.3148821 -6.0735859 Oct 2 616.1609 3.1678682 -5.1609070 Nov 2 604.1560 -7.2024828 -10.1560105 Dec 2 593.8949 -9.2907212 1.1050551 Jan 3 586.7861 -7.8008454 4.2138937 Feb 3 584.6794 -3.9076774 4.3205752 Mar 3 580.4915 -4.0974197 3.5085462 Apr 3 572.3764 -6.8079245 0.6236184 May 3 571.2755 -2.9376714 -4.2755279 Jun 3 581.1398 5.7227690 -12.1398075 Jul 3 604.7835 17.8162552 16.2165446 Aug 3 622.9143 18.0287501 6.0856782 Sep 3 631.2732 11.4919827 -3.2731859 Oct 3 619.0344 -4.5508446 -7.0343727 Nov 3 605.6026 -10.5515908 -10.6026424 Dec 3 596.1180 -9.8311042 0.8820151 Jan 4 589.5056 -7.6559416 3.4944099 Feb 4 583.8266 -6.3203680 6.1734198 Mar 4 574.3702 -8.4287426 5.6297692 Apr 4 571.9363 -4.4053037 2.0637058 May 4 578.7753 3.1630382 -5.7752700 Jun 4 590.4236 8.8725566 -17.4236162 Jul 4 604.4790 12.3531944 15.5210458 Aug 4 617.4134 12.7435338 8.5866224 Sep 4 618.9484 5.2109569 1.0515906 Oct 4 599.1281 -11.6137170 -11.1280708 Nov 4 578.6353 -17.5798043 -12.6352955 Dec 4 557.8685 -19.7210447 -0.8684780 Jan 5 552.6738 -9.9554152 8.3261825 Feb 5 541.4878 -10.7821494 7.5121978 Mar 5 527.8698 -12.6818517 4.1301643 Apr 5 524.3538 -6.5476709 1.6462395 May 5 520.0266 -5.0590382 -9.0265770 Jun 5 518.7608 -2.5152268 -19.7607848 Jul 5 535.8846 10.6385501 19.1153553 Aug 5 552.4114 14.5805360 12.5886215 Sep 5 540.4827 -3.1726379 1.5173276 Oct 5 534.1924 -5.2608748 -7.1923829 Nov 5 523.7013 -8.7637958 -13.7013220 Dec 5 518.3767 -6.4600258 -4.3766936 Jan 6 508.6638 -8.6398233 8.3361961 Feb 6 498.9714 -9.3446543 9.0285924 Mar 6 490.4645 -8.7847179 2.5354766 Apr 6 486.5635 -5.5225796 3.4365251 May 6 480.4165 -5.9401229 -11.4165193 Jun 6 498.6902 10.2600401 -20.6901510 Jul 6 511.4183 11.9100071 16.5816895 Aug 6 515.1410 6.4401180 18.8589786 Sep 6 516.5340 3.0681251 1.4660018 Oct 6 513.5876 -0.9507128 -7.5875801 Nov 6 515.8532 1.1986329 -13.8532485 Dec 6 518.8181 2.3792441 -2.8181084 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.19591743 -0.39223786 0.14409339 -0.19850345 0.62847050 2 1.22544629 -0.18601712 -0.66868505 -1.04188732 0.40983641 1.67404699 3 0.16729824 0.43647860 -0.02118779 -0.30502529 0.43438059 0.96763169 4 0.24393591 0.14949651 -0.23573780 0.45146219 0.84896933 0.63862003 5 1.09422817 -0.09251342 -0.21251827 0.68759018 0.16689237 0.28469801 6 -0.24412821 -0.07886892 0.06265542 0.36549103 -0.04679433 1.81360381 Jul Aug Sep Oct Nov Dec 1 3.99269022 -0.37803760 -2.50841294 -2.13730329 -1.01281832 1.44597794 2 1.41962649 0.76565402 -1.21064162 -1.47186502 -1.16094207 -0.23385382 3 1.35334777 0.02379634 -0.73176046 -1.79597672 -0.67178099 0.08071070 4 0.38936473 0.04370039 -0.84329004 -1.88348049 -0.66796634 -0.23986638 5 1.47152213 0.44122328 -1.98749946 -0.23378832 -0.39223334 0.25804694 6 0.18460897 -0.61216518 -0.37748692 -0.44996427 0.24069086 0.13222780 > 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/1w4w91259787556.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/2pmcc1259787556.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/3msnf1259787556.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/46r8f1259787556.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/5bvd21259787556.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/6qj8g1259787556.tab") > > system("convert tmp/1w4w91259787556.ps tmp/1w4w91259787556.png") > system("convert tmp/2pmcc1259787556.ps tmp/2pmcc1259787556.png") > system("convert tmp/3msnf1259787556.ps tmp/3msnf1259787556.png") > system("convert tmp/46r8f1259787556.ps tmp/46r8f1259787556.png") > system("convert tmp/5bvd21259787556.ps tmp/5bvd21259787556.png") > > > proc.time() user system elapsed 1.297 0.816 2.323