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(19876,45335,48674,156392,100837,101605,532850,294189,80763,105995,25045,90474,48481,50730,68694,207716,99132,104012,422632,364974,82687,66834,28408,97073,40284,24421,116346,72120,108751,91738,402216,390070,106045,110070,70668,167841,28607,95371,30605,131063,81214,85451,455196,454570,63114,74287,42350,113375) > 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 1161302758 2582741299 0 > m$fitted level slope sea Jan 1 19876.000 0.0000 0.000 Feb 1 23705.163 5616.7522 21629.837 Mar 1 49924.222 17299.1503 -1250.222 Apr 1 118643.277 41742.0878 37748.723 May 1 137795.319 30205.5178 -36958.319 Jun 1 126812.870 9132.4693 -25207.870 Jul 1 350459.558 115664.4211 182390.442 Aug 1 406741.200 86369.1059 -112552.200 Sep 1 255397.690 -31541.5340 -174634.690 Oct 1 118683.916 -83889.6420 -12688.916 Nov 1 11682.218 -95400.3098 13362.782 Dec 1 8578.678 -49443.8837 81895.322 Jan 2 23850.841 -17460.9106 24630.159 Feb 2 33836.855 -3814.1412 16893.145 Mar 2 63100.136 12555.2461 5593.864 Apr 2 129052.841 38285.1891 78663.159 May 2 177105.579 42972.6438 -77973.579 Jun 2 234895.571 50175.1171 -130883.571 Jul 2 227588.397 22243.2365 195043.603 Aug 2 288810.071 41023.6983 76163.929 Sep 2 258982.194 7003.7746 -176295.194 Oct 2 148865.315 -49297.1089 -82031.315 Nov 2 63317.181 -66742.2167 -34909.181 Dec 2 19313.624 -55796.3401 77759.376 Jan 3 4584.494 -36002.3266 35699.506 Feb 3 9726.434 -16144.4713 14694.566 Mar 3 82780.067 26748.6075 33565.933 Apr 3 64702.067 5351.5332 7417.933 May 3 134511.652 36051.6728 -25760.652 Jun 3 202486.610 51322.9654 -110748.610 Jul 3 238150.196 43811.5560 164065.804 Aug 3 263383.394 34925.1283 126686.606 Sep 3 245406.392 9731.4342 -139361.392 Oct 3 199958.699 -16491.9172 -89888.699 Nov 3 139516.722 -37387.9313 -68848.722 Dec 3 103895.282 -36546.6039 63945.718 Jan 4 43906.692 -47729.6669 -15299.692 Feb 4 61223.728 -16706.2918 34147.272 Mar 4 22335.559 -27253.3710 8269.441 Apr 4 89285.720 17368.7703 41777.280 May 4 132342.072 29524.6515 -51128.072 Jun 4 181030.376 38614.2063 -95579.376 Jul 4 252992.541 54464.0991 202203.459 Aug 4 302167.561 51952.9795 152402.439 Sep 4 245619.144 606.4161 -182505.144 Oct 4 168819.181 -35933.3516 -94532.181 Nov 4 106500.038 -48387.3760 -64150.038 Dec 4 45956.171 -54134.7015 67418.829 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.37779541 0.32349254 0.78922864 -0.36782051 -0.61368441 2 0.94418152 0.40715496 0.47827859 0.74938135 0.13975848 0.21413182 3 0.58279439 0.58354973 1.25454303 -0.62632275 0.90495428 0.45125178 4 -0.32886695 0.90989727 -0.30872058 1.30787451 0.35760775 0.26786147 Jul Aug Sep Oct Nov Dec 1 3.08504736 -0.86135391 -3.47143212 -1.53680028 -0.33769229 1.34858582 2 -0.81744738 0.54820571 -0.99739825 -1.65389196 -0.51260285 0.32190919 3 -0.22083746 -0.26004049 -0.73732098 -0.76951128 -0.61461037 0.02476725 4 0.46604631 -0.07358855 -1.50325370 -1.07172167 -0.36619019 -0.16914995 > 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/1w1xh1292764750.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/www/html/rcomp/tmp/2w1xh1292764750.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/www/html/rcomp/tmp/3psw11292764750.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/www/html/rcomp/tmp/4psw11292764750.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/www/html/rcomp/tmp/502wm1292764750.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/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/6etcv1292764750.tab") > > try(system("convert tmp/1w1xh1292764750.ps tmp/1w1xh1292764750.png",intern=TRUE)) character(0) > try(system("convert tmp/2w1xh1292764750.ps tmp/2w1xh1292764750.png",intern=TRUE)) character(0) > try(system("convert tmp/3psw11292764750.ps tmp/3psw11292764750.png",intern=TRUE)) character(0) > try(system("convert tmp/4psw11292764750.ps tmp/4psw11292764750.png",intern=TRUE)) character(0) > try(system("convert tmp/502wm1292764750.ps tmp/502wm1292764750.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.204 0.849 2.896