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(161,149,139,135,130,127,122,117,112,113,149,157,157,147,137,132,125,123,117,114,111,112,144,150,149,134,123,116,117,111,105,102,95,93,124,130,124,115,106,105,105,101,95,93,84,87,116,120,117,109,105,107,109,109,108,107,99,103,131,137) > 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 4.728582 33.281372 4.442201 0.000000 > m$fitted level slope sea Jan 1 161.00000 0.0000000 0.00000000 Feb 1 150.65280 -7.5470350 -1.65280495 Mar 1 138.98266 -10.5000706 0.01734429 Apr 1 133.86416 -6.5805255 1.13583866 May 1 129.95414 -4.6558722 0.04586033 Jun 1 126.98887 -3.4335796 0.01113100 Jul 1 122.48080 -4.2114182 -0.48080246 Aug 1 117.29904 -4.9137827 -0.29903963 Sep 1 112.15806 -5.0782044 -0.15805972 Oct 1 112.15298 -1.4069536 0.84702006 Nov 1 143.07098 21.9860350 5.92901919 Dec 1 160.04781 18.3609133 -3.04781350 Jan 2 159.90360 5.1029640 -2.90359981 Feb 2 149.56081 -5.5817287 -2.56080666 Mar 2 138.24535 -9.4932981 -1.24535241 Apr 2 130.54711 -8.2385509 1.45289185 May 2 124.92995 -6.4150610 0.07004873 Jun 2 122.31113 -3.7771323 0.68886845 Jul 2 117.44246 -4.5372067 -0.44246482 Aug 2 112.64461 -4.7188050 1.35538681 Sep 2 110.20499 -3.1308828 0.79501378 Oct 2 116.88362 3.7031908 -4.88362097 Nov 2 135.55933 14.1277381 8.44066512 Dec 2 150.24913 14.5182840 -0.24913285 Jan 3 151.24142 5.0979810 -2.24141694 Feb 3 138.21639 -7.3716296 -4.21639005 Mar 3 125.00078 -11.3489070 -2.00077875 Apr 3 114.70623 -10.6226326 1.29376645 May 3 115.50666 -2.7521551 1.49334246 Jun 3 110.70887 -4.1579063 0.29112579 Jul 3 104.90382 -5.2912328 0.09618369 Aug 3 99.95160 -5.0577121 2.04839874 Sep 3 95.99942 -4.2960379 -0.99941958 Oct 3 100.57595 1.8151441 -7.57594736 Nov 3 114.71673 10.2949290 9.28326913 Dec 3 127.45892 11.9773661 2.54107702 Jan 4 124.77855 1.8823059 -0.77854501 Feb 4 118.49747 -3.7034384 -3.49747238 Mar 4 109.17403 -7.5246303 -3.17403236 Apr 4 105.75553 -4.7157422 -0.75552627 May 4 102.74295 -3.5485604 2.25705210 Jun 4 99.86328 -3.0911935 1.13672055 Jul 4 94.65886 -4.5365574 0.34114435 Aug 4 90.23964 -4.4562061 2.76035513 Sep 4 86.90851 -3.6855482 -2.90850818 Oct 4 95.85870 4.9662572 -8.85870249 Nov 4 107.16189 9.3009290 8.83810945 Dec 4 114.58960 8.0194499 5.41040296 Jan 5 116.77775 4.0267761 0.22224875 Feb 5 112.57208 -1.5861985 -3.57207736 Mar 5 109.46706 -2.6178130 -4.46705955 Apr 5 108.08167 -1.7778688 -1.08167361 May 5 106.78939 -1.4462795 2.21061247 Jun 5 106.60158 -0.5881911 2.39841814 Jul 5 106.71783 -0.1080146 1.28217253 Aug 5 104.47596 -1.5637824 2.52403708 Sep 5 104.86123 -0.2335596 -5.86122737 Oct 5 112.41802 5.0813649 -9.41802440 Nov 5 121.52837 7.8284185 9.47163125 Dec 5 130.13688 8.3605269 6.86311600 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.63771583 -0.52693971 0.66855028 0.33232247 0.21195430 2 -2.38883808 -1.81943303 -0.68705583 0.21768047 0.31426668 0.45708604 3 -1.64794338 -2.14798335 -0.69240953 0.12637748 1.35936481 -0.24335361 4 -1.75505818 -0.96512570 -0.66353860 0.48851457 0.20197230 0.07915837 5 -0.69282149 -0.97101271 -0.17897944 0.14595564 0.05744035 0.14853997 Jul Aug Sep Oct Nov Dec 1 -0.13482803 -0.12174772 -0.02850092 0.63637523 4.05494630 -0.62837947 2 -0.13177289 -0.03147664 0.27525158 1.18467504 1.80694467 0.06783432 3 -0.19648851 0.04047779 0.13202939 1.05934968 1.47009362 0.29234737 4 -0.25053228 0.01392859 0.13358574 1.49970820 0.75166242 -0.22256430 5 0.08321189 -0.25235087 0.23058081 0.92131739 0.47642336 0.09236487 > 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/1oexb1259697887.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/2kc7f1259697887.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/32dvd1259697887.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/4y5011259697887.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/5nwl61259697887.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/6sw9y1259697887.tab") > system("convert tmp/1oexb1259697887.ps tmp/1oexb1259697887.png") > system("convert tmp/2kc7f1259697887.ps tmp/2kc7f1259697887.png") > system("convert tmp/32dvd1259697887.ps tmp/32dvd1259697887.png") > system("convert tmp/4y5011259697887.ps tmp/4y5011259697887.png") > system("convert tmp/5nwl61259697887.ps tmp/5nwl61259697887.png") > > > proc.time() user system elapsed 1.299 0.798 1.513