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(216234,213587,209465,204045,200237,203666,241476,260307,243324,244460,233575,237217,235243,230354,227184,221678,217142,219452,256446,265845,248624,241114,229245,231805,219277,219313,212610,214771,211142,211457,240048,240636,230580,208795,197922,194596,194581,185686,178106,172608,167302,168053,202300,202388,182516,173476,166444,171297,169701,164182,161914,159612,151001,158114,186530,187069,174330,169362,166827,178037,186412,189226,191563,188906,186005,195309,223532,226899,214126) > 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 31709268 10119994 0 > m$fitted level slope sea Jan 1 216234.0 0.0000 0.0000 Feb 1 214477.2 -1830.3871 -890.1816 Mar 1 209827.8 -3680.1839 -362.7763 Apr 1 204418.5 -4877.8123 -373.5261 May 1 199912.7 -4623.0269 324.2752 Jun 1 201577.9 -375.5222 2088.0792 Jul 1 231683.7 20286.0819 9792.2518 Aug 1 260697.1 26211.6814 -390.0843 Sep 1 255302.9 4750.9617 -11978.8590 Oct 1 245969.5 -4810.3439 -1509.5073 Nov 1 234448.1 -9365.4069 -873.1219 Dec 1 233621.7 -3569.8250 3595.2709 Jan 2 234449.0 -594.9376 794.0007 Feb 2 231487.7 -2202.7724 -1133.7337 Mar 2 226389.0 -4111.7988 794.9508 Apr 2 220112.2 -5523.9872 1565.8346 May 2 217319.2 -3722.6174 -177.2076 Jun 2 225574.1 4118.6324 -6122.0737 Jul 2 246288.9 14948.5735 10157.1440 Aug 2 258828.0 13373.1862 7017.0450 Sep 2 259837.8 5281.7206 -11213.7629 Oct 2 246096.8 -7168.5229 -4982.8474 Nov 2 233605.5 -10650.5840 -4360.4987 Dec 2 227827.2 -7467.0254 3977.8071 Jan 3 218543.6 -8654.2315 733.3755 Feb 3 215992.2 -4660.4496 3320.8317 Mar 3 210180.0 -5407.5408 2430.0493 Apr 3 210854.2 -1481.9664 3916.8190 May 3 214841.0 2066.7506 -3699.0252 Jun 3 223122.9 6096.0388 -11665.9407 Jul 3 230442.0 6886.2516 9606.0432 Aug 3 231564.4 3161.0779 9071.6109 Sep 3 234650.4 3112.4973 -4070.3841 Oct 3 218734.2 -9202.4661 -9939.1872 Nov 3 204202.9 -12649.5277 -6280.9052 Dec 3 190122.4 -13574.7295 4473.5765 Jan 4 190251.1 -4707.7443 4329.8981 Feb 4 182988.4 -6360.5418 2697.5704 Mar 4 176859.1 -6211.6132 1246.8787 Apr 4 170764.9 -6136.1897 1843.1154 May 4 171497.3 -1710.9663 -4195.3391 Jun 4 178477.3 3891.6514 -10424.2915 Jul 4 190018.1 8812.9157 12281.8932 Aug 4 195001.1 6350.4349 7386.9100 Sep 4 184852.7 -4263.7633 -2336.7243 Oct 4 180114.0 -4569.3869 -6638.0367 Nov 4 173617.8 -5809.1287 -7173.7994 Dec 4 170483.2 -4088.2522 813.8327 Jan 5 164676.6 -5194.5454 5024.4208 Feb 5 160119.7 -4784.2723 4062.3078 Mar 5 158822.5 -2546.5907 3091.4810 Apr 5 159426.4 -527.4639 185.6016 May 5 159375.0 -221.8611 -8374.0495 Jun 5 168869.0 6019.4993 -10754.9993 Jul 5 173678.9 5243.3242 12851.0789 Aug 5 175021.3 2742.3721 12047.7442 Sep 5 176604.2 1999.0350 -2274.2025 Oct 5 176282.1 510.6004 -6920.1253 Nov 5 175188.7 -518.0857 -8361.7003 Dec 5 175800.6 206.7833 2236.4411 Jan 6 180014.1 2778.0245 6397.9271 Feb 6 185377.1 4435.7082 3848.8809 Mar 6 189375.2 4155.5724 2187.8052 Apr 6 190385.5 2143.7915 -1479.5144 May 6 196606.2 4754.0918 -10601.2229 Jun 6 204886.3 7013.3561 -9577.2727 Jul 6 210211.3 5932.2150 13320.7070 Aug 6 214691.9 5003.5365 12207.0808 Sep 6 216343.4 2859.7123 -2217.3880 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.394798933 -0.339948722 -0.224145044 0.044785362 0.754017126 2 0.530187415 -0.290158594 -0.336336347 -0.254681917 0.321211514 1.387083663 3 -0.211501835 0.710392026 -0.132288677 0.700184461 0.632709486 0.714353568 4 1.577339009 -0.293491191 0.026410818 0.013419154 0.787870194 0.994589206 5 -0.196622435 0.072828188 0.397062001 0.358942433 0.054364249 1.108551633 6 0.456780488 0.294238190 -0.049721925 -0.357512666 0.464137381 0.401351579 Jul Aug Sep Oct Nov Dec 1 3.674339895 1.052243326 -3.811087217 -1.697983737 -0.808911182 1.029209687 2 1.923749107 -0.279940690 -1.436814428 -2.211089482 -0.618383957 0.565603425 3 0.140188518 -0.661749837 -0.008627875 -2.186945497 -0.612189103 -0.164426695 4 0.873031646 -0.437232024 -1.885098760 -0.054277662 -0.220209076 0.305824936 5 -0.137731494 -0.443955593 -0.132007691 -0.264371079 -0.182749364 0.128806520 6 -0.191894446 -0.164841281 -0.380687410 > 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/1uetd1259790489.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/2tyjb1259790489.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/3hskk1259790489.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/4sb8i1259790489.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/5kkth1259790489.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/6tgzd1259790489.tab") > system("convert tmp/1uetd1259790489.ps tmp/1uetd1259790489.png") > system("convert tmp/2tyjb1259790489.ps tmp/2tyjb1259790489.png") > system("convert tmp/3hskk1259790489.ps tmp/3hskk1259790489.png") > system("convert tmp/4sb8i1259790489.ps tmp/4sb8i1259790489.png") > system("convert tmp/5kkth1259790489.ps tmp/5kkth1259790489.png") > > > proc.time() user system elapsed 1.296 0.815 3.291