R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-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(61,65,55,56,91,80,135,129,129,130,109,126,73,68,74,95,105,108,127,108,126,154,127,103,95,59,68,82,92,124,139,167,138,146,128,145,91,66,89,98,113,130,127,157,157,136,145,112,71,95,95,105,116,104,128,181,130,124,123,152) > par1 = '12' > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 413.470584 0.000000 1.633987 54.658269 > m$fitted level slope sea Jan 1 61.00000 0.0000000 0.000000000 Feb 1 64.62446 0.2052716 0.212837861 Mar 1 56.01052 0.1035777 0.007210297 Apr 1 55.89514 0.1022914 0.130496053 May 1 86.81001 0.2668353 0.579509734 Jun 1 80.68808 0.2331597 0.060563464 Jul 1 128.44541 0.4821988 0.985576070 Aug 1 128.71276 0.4810788 0.312420939 Sep 1 128.59721 0.4779849 0.472699233 Oct 1 129.46791 0.4800109 0.486075461 Nov 1 111.02718 0.3829017 0.190032092 Dec 1 123.89736 0.4466644 0.639328065 Jan 2 96.47994 1.8090456 -20.192843156 Feb 2 69.54219 1.0740659 0.896410736 Mar 2 72.51284 1.0853265 1.268507373 Apr 2 91.67526 1.1354501 1.222585469 May 2 102.19950 1.1591488 1.708370495 Jun 2 107.01677 1.1683407 0.557701458 Jul 2 123.26418 1.2061568 1.981771779 Aug 2 109.18195 1.1679072 0.596444836 Sep 2 123.02485 1.1995461 1.500753788 Oct 2 149.22022 1.2619191 1.872173607 Nov 2 129.24622 1.2082755 0.224027658 Dec 2 104.74398 1.1524178 1.247015844 Jan 3 107.11761 1.1246551 -12.260469893 Feb 3 63.76743 0.4003080 -0.486926570 Mar 3 66.96989 0.4127523 0.708866650 Apr 3 79.85167 0.4362281 0.707952561 May 3 89.76560 0.4522180 1.139334871 Jun 3 120.21459 0.5030730 0.319599918 Jul 3 135.50826 0.5281777 1.782844282 Aug 3 163.35083 0.5744794 0.493293139 Sep 3 141.25252 0.5360823 -0.632925616 Oct 3 144.05183 0.5399378 1.686679191 Nov 3 130.16134 0.5145671 -0.494053862 Dec 3 142.53458 0.5277632 1.095091632 Jan 4 106.89098 1.0511181 -11.684976901 Feb 4 70.86000 0.6157214 -1.092017963 Mar 4 85.70072 0.6698054 1.677824663 Apr 4 95.81049 0.6839289 1.105337572 May 4 110.30479 0.7017792 1.108756122 Jun 4 127.40565 0.7232966 0.710578969 Jul 4 126.32911 0.7209220 0.877635572 Aug 4 151.98007 0.7538005 2.156242203 Sep 4 157.16354 0.7596464 -0.672369090 Oct 4 136.95349 0.7315386 1.455281130 Nov 4 144.55236 0.7411271 -0.341216260 Dec 4 114.62422 0.7253662 0.899267414 Jan 5 85.52974 1.0340084 -11.079775188 Feb 5 94.84364 1.1087049 -0.710288826 Mar 5 93.47609 1.1001365 1.804745651 Apr 5 102.51854 1.1105218 1.573672799 May 5 113.52146 1.1210876 1.347577791 Jun 5 104.40922 1.1099037 0.760669113 Jul 5 124.44360 1.1308249 1.392970127 Aug 5 171.72262 1.1818919 4.001770939 Sep 5 136.01393 1.1409291 -1.796639260 Oct 5 124.20335 1.1261424 1.277249502 Nov 5 122.37653 1.1226634 0.961057383 Dec 5 146.34320 1.1253756 3.044906067 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.10911964 -0.43215788 -0.01078797 1.51832321 -0.31481266 2 -1.67017072 -1.19853607 0.09305288 0.88980153 0.46216075 0.18006982 3 0.06640213 -1.98304531 0.13749121 0.61358953 0.46639586 1.47612393 4 -1.89862985 -1.70426583 0.69782626 0.46447447 0.67951239 0.80687424 5 -1.53641389 0.38724491 -0.12144073 0.39075850 0.48670021 -0.50346439 Jul Aug Sep Oct Nov Dec 1 2.34177880 -0.01058706 -0.02939820 0.01935021 -0.93229185 0.61528624 2 0.74225900 -0.75256064 0.62391674 1.23040481 -1.04536334 -1.26531561 3 0.72783558 1.34412087 -1.11571155 0.11137321 -0.71016479 0.58335696 4 -0.08855603 1.22661289 0.21794986 -1.03177946 0.33792215 -1.50839401 5 0.93105142 2.27040949 -1.81495852 -0.63721497 -0.14529671 1.12345460 > 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/1ocbc1324659306.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/2apqy1324659306.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/322cv1324659306.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/4mno41324659306.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/5elj91324659306.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/6wkzp1324659306.tab") > > try(system("convert tmp/1ocbc1324659306.ps tmp/1ocbc1324659306.png",intern=TRUE)) character(0) > try(system("convert tmp/2apqy1324659306.ps tmp/2apqy1324659306.png",intern=TRUE)) character(0) > try(system("convert tmp/322cv1324659306.ps tmp/322cv1324659306.png",intern=TRUE)) character(0) > try(system("convert tmp/4mno41324659306.ps tmp/4mno41324659306.png",intern=TRUE)) character(0) > try(system("convert tmp/5elj91324659306.ps tmp/5elj91324659306.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.976 0.341 2.318