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(5560,3922,3759,4138,4634,3996,4308,4143,4429,5219,4929,5755,5592,4163,4962,5208,4755,4491,5732,5731,5040,6102,4904,5369,5578,4619,4731,5011,5299,4146,4625,4736,4219,5116,4205,4121,5103,4300,4578,3809,5526,4247,3830,4394,4826,4409,4569,4106,4794,3914,3793,4405,4022,4100,4788,3163,3585,3903,4178,3863,4187) > 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 59366.816 0.000 2263.590 164217.464 > m$fitted level slope sea Jan 1 5560.000 0.00000 0.0000000 Feb 1 4714.525 -74.86342 -89.9284610 Mar 1 4251.253 -98.59778 -102.8767432 Apr 1 4188.743 -97.00692 -90.2626406 May 1 4380.814 -87.00910 -76.2344692 Jun 1 4201.358 -89.71068 -97.9182549 Jul 1 4238.676 -86.39521 -79.8413625 Aug 1 4187.542 -85.54220 -86.1848467 Sep 1 4285.768 -81.32555 -74.5169003 Oct 1 4692.917 -70.54702 -53.8980841 Nov 1 4793.511 -66.88599 -67.9896016 Dec 1 5214.388 -56.72257 -39.8752156 Jan 2 5103.991 -51.40162 561.7679919 Feb 2 4639.175 -66.45553 -93.6842779 Mar 2 4780.642 -60.71215 -40.0207282 Apr 2 4968.636 -55.64366 -43.2199308 May 2 4864.981 -56.42269 -54.0020854 Jun 2 4697.800 -57.96165 -76.1980433 Jul 2 5140.713 -51.66599 -2.5062126 Aug 2 5402.432 -47.96984 -43.9732834 Sep 2 5244.574 -49.21386 -73.7482507 Oct 2 5607.408 -44.67574 3.6134465 Nov 2 5305.700 -47.44934 -95.2663634 Dec 2 5316.260 -46.87256 -16.5047300 Jan 3 5190.951 -43.37087 487.4889975 Feb 3 4940.272 -48.04046 -108.0265387 Mar 3 4833.819 -49.12281 -38.5477609 Apr 3 4897.780 -47.56521 -16.5996155 May 3 5064.961 -45.21092 -16.8562203 Jun 3 4682.332 -48.38029 -139.0805842 Jul 3 4630.156 -48.41260 -0.6713349 Aug 3 4665.292 -47.74328 -28.1963504 Sep 3 4487.007 -48.75061 -113.3200237 Oct 3 4724.327 -46.59478 52.5057738 Nov 3 4520.640 -47.75416 -129.3190747 Dec 3 4322.456 -48.66308 -22.7167990 Jan 4 4375.383 -51.45881 601.1161495 Feb 4 4363.209 -50.84105 -105.4525042 Mar 4 4452.502 -48.84612 -30.5041847 Apr 4 4162.496 -51.43451 -76.3579461 May 4 4716.954 -46.26684 103.0092406 Jun 4 4560.848 -47.06765 -185.0714711 Jul 4 4219.384 -49.01170 -43.2350545 Aug 4 4268.432 -48.40189 10.1047778 Sep 4 4540.355 -46.47990 -91.8186119 Oct 4 4437.288 -46.81203 38.4224060 Nov 4 4516.674 -46.09349 -96.5020133 Dec 4 4329.468 -46.65314 -56.6796655 Jan 5 4263.101 -46.27918 555.0315019 Feb 5 4133.084 -47.25788 -126.9046587 Mar 5 3959.753 -48.73244 -25.6541133 Apr 5 4175.791 -46.34690 -74.6660196 May 5 4040.998 -46.98003 83.7232869 Jun 5 4105.324 -46.30124 -135.2786297 Jul 5 4395.028 -44.44891 -0.2994514 Aug 5 3857.744 -47.00696 -117.2087037 Sep 5 3739.528 -47.36389 -71.0265917 Oct 5 3768.219 -46.99113 45.5643284 Nov 5 3954.635 -45.90234 -50.5506900 Dec 5 3918.677 -45.87564 -67.3752714 Jan 6 3774.590 -44.52052 531.1934173 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -2.64728315 -1.39053745 0.13751198 1.13290462 -0.36734508 2 -0.29104680 -1.39183124 0.77900677 0.97752437 -0.19230513 -0.44726682 3 -0.36949544 -0.75801131 -0.22552495 0.45062433 0.86684071 -1.36962201 4 0.45504115 0.14920227 0.54858693 -0.96694443 2.45436570 -0.44691610 5 -0.08598670 -0.32507400 -0.49757336 1.06523696 -0.35897535 0.45349353 6 -0.42135038 Jul Aug Sep Oct Nov Dec 1 0.50842625 0.14169042 0.74012362 1.97003267 0.69083860 1.97015701 2 2.03044740 1.27287580 -0.44679977 1.67636960 -1.04609019 0.23619741 3 -0.01544677 0.34046904 -0.53234443 1.16703873 -0.64102463 -0.61419131 4 -1.20033228 0.40022700 1.30810653 -0.23115726 0.51562213 -0.57700605 5 1.37142306 -2.01329649 -0.29103431 0.31092235 0.95438798 0.04069952 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/106d01259937230.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/2p1601259937230.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/3tmi31259937230.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/4loyq1259937230.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/5c5xc1259937230.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/6u1lt1259937230.tab") > system("convert tmp/106d01259937230.ps tmp/106d01259937230.png") > system("convert tmp/2p1601259937230.ps tmp/2p1601259937230.png") > system("convert tmp/3tmi31259937230.ps tmp/3tmi31259937230.png") > system("convert tmp/4loyq1259937230.ps tmp/4loyq1259937230.png") > system("convert tmp/5c5xc1259937230.ps tmp/5c5xc1259937230.png") > > > proc.time() user system elapsed 1.472 0.818 1.676