R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-pc-linux-gnu (32-bit) 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(15579,16348,15928,16171,15937,15713,15594,15683,16438,17032,17696,17745,19394,20148,20108,18584,18441,18391,19178,18079,18483,19644,19195,19650,20830,23595,22937,21814,21928,21777,21383,21467,22052,22680,24320,24977,25204,25739,26434,27525,30695,32436,30160,30236,31293,31077,32226,33865,32810,32242,32700,32819,33947,34148,35261,39506,41591,39148,41216,40225,41126,42362,40740,40256,39804,41002,41702,42254,43605,43271) > par1 = '12' > 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 1161150.195 3126.726 6068.639 0.000 > m$fitted level slope sea Jan 1 15579.00 0.00000 0.000000 Feb 1 16306.20 40.45360 41.804137 Mar 1 15898.80 35.84692 29.201096 Apr 1 16127.57 38.06299 43.428353 May 1 15904.77 34.36191 32.230496 Jun 1 15679.41 30.04543 33.585606 Jul 1 15559.31 27.19937 34.689400 Aug 1 15645.48 28.44939 37.519369 Sep 1 16389.07 45.14145 48.929802 Oct 1 16986.94 59.17020 45.061057 Nov 1 17647.37 75.59338 48.629787 Dec 1 17705.99 75.09891 39.005083 Jan 2 19196.54 66.65010 197.457780 Feb 2 20124.25 107.97279 23.747083 Mar 2 20128.23 104.44943 -20.230788 Apr 2 18619.96 49.59951 -35.957082 May 2 18433.33 41.20337 7.671152 Jun 2 18389.92 38.08554 1.078403 Jul 2 19153.01 65.66230 24.989618 Aug 2 18122.66 22.77032 -43.657243 Sep 2 18458.04 35.32081 24.957146 Oct 2 19616.71 81.46923 27.293328 Nov 2 19203.72 60.70848 -8.719248 Dec 2 19656.60 77.35327 -6.603968 Jan 3 20640.08 107.14936 189.918608 Feb 3 23499.31 245.82214 95.691672 Mar 3 22953.90 209.91604 -16.897859 Apr 3 21866.29 151.55110 -52.294206 May 3 21881.65 145.34633 46.353062 Jun 3 21792.03 134.52680 -15.032597 Jul 3 21325.78 106.60789 57.215571 Aug 3 21481.00 108.88554 -14.003270 Sep 3 22025.64 129.45430 26.358301 Oct 3 22603.51 150.76567 76.486788 Nov 3 24252.67 222.46640 67.330289 Dec 3 24984.96 246.82473 -7.960862 Jan 4 25451.07 256.58157 -247.067914 Feb 4 25632.54 252.75184 106.459459 Mar 4 26381.03 277.13648 52.972845 Apr 4 27541.91 320.26583 -16.909352 May 4 30532.94 451.06032 162.059865 Jun 4 32398.57 520.56990 37.427149 Jul 4 30251.04 389.10913 -91.043979 Aug 4 30244.82 369.58410 -8.821335 Sep 4 31230.12 400.05948 62.878912 Oct 4 31088.06 373.17004 -11.057078 Nov 4 32139.89 406.90286 86.110951 Dec 4 33829.86 470.46756 35.141458 Jan 5 33198.51 417.22867 -388.507294 Feb 5 32196.54 345.03831 45.456366 Mar 5 32629.63 349.46029 70.368566 Apr 5 32952.17 348.11500 -133.174166 May 5 33853.30 375.77363 93.698273 Jun 5 34011.81 364.89503 136.192210 Jul 5 35279.72 410.14664 -18.720658 Aug 5 39363.95 594.38831 142.049767 Sep 5 41467.01 670.09109 123.991522 Oct 5 39360.25 530.65606 -212.254441 Nov 5 41106.34 591.72173 109.664844 Dec 5 40238.79 518.59790 -13.788529 Jan 6 41301.73 545.61682 -175.730174 Feb 6 42300.70 568.59083 61.301700 Mar 6 40795.39 463.96957 -55.386497 Apr 6 40440.41 422.75993 -184.414130 May 6 39726.49 365.56151 77.512510 Jun 6 40810.11 401.70829 191.886000 Jul 6 41876.31 435.16728 -174.314502 Aug 6 42236.81 431.40682 17.185323 Sep 6 43297.46 463.10355 307.540020 Oct 6 43583.24 454.16910 -312.242901 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.40558794 -0.41537387 0.17866485 -0.24125885 -0.23988189 2 1.55207641 0.67420487 -0.09500761 -1.47178624 -0.21542391 -0.07710304 3 0.89673831 2.29338580 -0.71742523 -1.17570444 -0.12338377 -0.21280233 4 0.20936745 -0.06433148 0.44811147 0.79860931 2.41329121 1.27818024 5 -1.03331505 -1.23457281 0.07947853 -0.02430149 0.49919965 -0.19614011 6 0.50508229 0.39825852 -1.87042666 -0.73911543 -1.02561932 0.64796912 Jul Aug Sep Oct Nov Dec 1 -0.13850450 0.05432969 0.65806865 0.50803050 0.55202912 -0.01556284 2 0.66021763 -0.99743658 0.28432843 1.02113959 -0.44928268 0.35600042 3 -0.54397918 0.04400502 0.39437999 0.40577039 1.35582712 0.46087354 4 -2.41066317 -0.35715978 0.55622622 -0.48973449 0.61311546 1.15810976 5 0.81520323 3.31671108 1.36190332 -2.50682160 1.09726714 -1.31651471 6 0.59964314 -0.06737987 0.56783185 -0.16002790 > 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/wessaorg/rcomp/tmp/1aa701353683588.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2ip0u1353683588.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3hpwd1353683589.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4b2zi1353683589.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5skam1353683589.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/64g3u1353683589.tab") > > try(system("convert tmp/1aa701353683588.ps tmp/1aa701353683588.png",intern=TRUE)) character(0) > try(system("convert tmp/2ip0u1353683588.ps tmp/2ip0u1353683588.png",intern=TRUE)) character(0) > try(system("convert tmp/3hpwd1353683589.ps tmp/3hpwd1353683589.png",intern=TRUE)) character(0) > try(system("convert tmp/4b2zi1353683589.ps tmp/4b2zi1353683589.png",intern=TRUE)) character(0) > try(system("convert tmp/5skam1353683589.ps tmp/5skam1353683589.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.558 0.542 4.384