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(6802.96,7132.68,7073.29,7264.5,7105.33,7218.71,7225.72,7354.25,7745.46,8070.26,8366.33,8667.51,8854.34,9218.1,9332.9,9358.31,9248.66,9401.2,9652.04,9957.38,10110.63,10169.26,10343.78,10750.21,11337.5,11786.96,12083.04,12007.74,11745.93,11051.51,11445.9,11924.88,12247.63,12690.91,12910.7,13202.12,13654.67,13862.82,13523.93,14211.17,14510.35,14289.23,14111.82,13086.59,13351.54,13747.69,12855.61,12926.93,12121.95,11731.65,11639.51,12163.78,12029.53,11234.18,9852.13,9709.04,9332.75,7108.6,6691.49,6143.05) > 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 172017.897 7571.409 0.000 0.000 > m$fitted level slope sea Jan 1 6802.960 0.00000 0.000000 Feb 1 7115.776 21.54037 16.904368 Mar 1 7056.524 15.06189 16.765832 Apr 1 7247.466 34.50371 17.034079 May 1 7088.551 8.60187 16.778588 Jun 1 7201.814 24.41570 16.895948 Jul 1 7208.840 21.57615 16.879633 Aug 1 7337.287 39.93003 16.962683 Sep 1 7728.273 102.23836 17.187069 Oct 1 8052.957 142.58010 17.303465 Nov 1 8348.961 170.80036 17.368979 Dec 1 8650.096 194.99665 17.414301 Jan 2 8986.284 219.97303 -131.943909 Feb 2 9206.108 219.94582 11.991542 Mar 2 9320.848 200.16300 12.052172 Apr 2 9346.176 167.23635 12.133955 May 2 9236.421 115.01637 12.239096 Jun 2 9388.972 122.09764 12.227536 Jul 2 9639.845 146.40323 12.195363 Aug 2 9945.217 176.41741 12.163144 Sep 2 10098.463 172.04162 12.166953 Oct 2 10157.078 150.61826 12.182079 Nov 2 10331.601 155.13362 12.179493 Dec 2 10738.053 202.60935 12.157446 Jan 3 11382.614 284.88534 -45.113995 Feb 3 11781.210 305.92423 5.749675 Mar 3 12077.286 304.06262 5.754205 Apr 3 12001.844 232.34408 5.895892 May 3 11739.884 138.94255 6.045569 Jun 3 11045.260 -18.55343 6.250285 Jul 3 11439.732 59.48151 6.168015 Aug 3 11918.780 138.75075 6.100231 Sep 3 12241.554 173.51832 6.076118 Oct 3 12684.863 224.48990 6.047445 Nov 3 12904.652 223.60186 6.047850 Dec 3 13196.077 236.41579 6.043109 Jan 4 13687.617 284.18321 -32.947239 Feb 4 13861.010 263.57796 1.810396 Mar 4 13521.911 149.66794 2.019041 Apr 4 14209.302 251.28389 1.868058 May 4 14508.493 260.33617 1.857149 Jun 4 14287.284 169.35116 1.946090 Jul 4 14109.822 103.82506 1.998043 Aug 4 13084.455 -109.51858 2.135239 Sep 4 13349.442 -38.76192 2.098334 Oct 4 13745.626 43.41406 2.063571 Nov 4 12853.486 -133.34397 2.124218 Dec 4 12924.817 -94.67382 2.113457 Jan 5 12206.165 -211.77973 -84.215063 Feb 5 11726.266 -261.80881 5.383752 Mar 5 11634.173 -229.73365 5.337021 Apr 5 12158.611 -87.21936 5.168620 May 5 12024.353 -96.10752 5.177138 Jun 5 11228.900 -228.24423 5.279861 Jul 5 9846.713 -446.26863 5.417331 Aug 5 9703.652 -388.98169 5.388034 Sep 5 9327.363 -386.58359 5.387040 Oct 5 7103.096 -733.78652 5.503845 Nov 5 6686.002 -673.95200 5.487519 Dec 5 6137.568 -650.23726 5.482271 > m$resid Jan Feb Mar Apr May 1 0.0000000000 0.4493203984 -0.1871573103 0.4005856645 -0.4346183776 2 0.3596350518 -0.0002727981 -0.2284252020 -0.3795805730 -0.6013578692 3 1.0414088557 0.2235802836 -0.0213735504 -0.8236917802 -1.0729575015 4 0.5855231857 -0.2241044450 -1.3080917121 1.1672205143 0.1039978480 5 -1.4128708870 -0.5509298714 0.3683947167 1.6371721877 -0.1021193617 Jun Jul Aug Sep Oct 1 0.2327857478 -0.0383911627 0.2347081544 0.7681879357 0.4855525550 2 0.0814900739 0.2795764713 0.3451357006 -0.0503076431 -0.2462680958 3 -1.8095086162 0.8966470809 0.9108863730 0.3995321381 0.5857565418 4 -1.0454085286 -0.7529447970 -2.4516025197 0.8131146728 0.9443624206 5 -1.5183053320 -2.5053403557 0.6583159790 0.0275585458 -3.9900726620 Nov Dec 1 0.3343305056 0.2836914752 2 0.0519010639 0.5456703927 3 -0.0102053043 0.1472597583 4 -2.0313233189 0.4444059193 5 0.6876289233 0.2725359801 > 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/1bf481259956440.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/22kyn1259956440.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/3gcj21259956440.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/4ekd81259956440.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/5byp71259956440.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/62q1x1259956440.tab") > system("convert tmp/1bf481259956440.ps tmp/1bf481259956440.png") > system("convert tmp/22kyn1259956440.ps tmp/22kyn1259956440.png") > system("convert tmp/3gcj21259956440.ps tmp/3gcj21259956440.png") > system("convert tmp/4ekd81259956440.ps tmp/4ekd81259956440.png") > system("convert tmp/5byp71259956440.ps tmp/5byp71259956440.png") > > > proc.time() user system elapsed 1.371 0.835 1.624