R version 2.12.1 (2010-12-16) Copyright (C) 2010 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(68897,38683,44720,39525,45315,50380,40600,36279,42438,38064,31879,11379,70249,39253,47060,41697,38708,49267,39018,32228,40870,39383,34571,12066,70938,34077,45409,40809,37013,44953,37848,32745,39401,34931,33008,8620,68906,39556,50669,36432,40891,48428,36222,33425,39401,37967,34801,12657,69116,41519,51321,38529,41547,52073,38401,40898,40439,41888,37898,8771,68184,50530,47221,41756,45633,48138,39486,39341,41117,41629,29722,7054,56676,34870,35117,30169,30936,35699,33228,27733,33666,35429,27438,8170) > 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 15864827.87 27437.62 8288058.17 0.00 > m$fitted level slope sea Jan 1 68897.00 0.00000 0.0000 Feb 1 46220.02 -1041.16529 -7537.0247 Mar 1 42974.65 -1154.48649 1745.3509 Apr 1 40880.68 -1184.14883 -1355.6823 May 1 43436.64 -1112.99065 1878.3574 Jun 1 48260.32 -1018.01800 2119.6829 Jul 1 45041.50 -1055.69635 -4441.4958 Aug 1 39770.42 -1136.00065 -3491.4181 Sep 1 40758.55 -1091.99185 1679.4474 Oct 1 39738.63 -1090.40167 -1674.6321 Nov 1 35371.15 -1166.53681 -3492.1509 Dec 1 21050.27 -1486.44291 -9671.2743 Jan 2 35363.05 -1624.07398 34885.9468 Feb 2 42236.18 -1447.55171 -2983.1815 Mar 2 45047.44 -1287.79845 2012.5627 Apr 2 45707.22 -1213.41443 -4010.2237 May 2 42346.35 -1285.51003 -3638.3543 Jun 2 43824.02 -1202.20094 5442.9840 Jul 2 43008.31 -1190.94399 -3990.3095 Aug 2 39432.03 -1261.55176 -7204.0338 Sep 2 38432.88 -1253.52757 2437.1172 Oct 2 37892.71 -1231.22625 1490.2942 Nov 2 35216.12 -1275.65883 -645.1245 Dec 2 31249.69 -1351.58898 -19183.6913 Jan 3 35296.80 -1215.79509 35641.2039 Feb 3 37673.98 -1105.04792 -3596.9830 Mar 3 41317.09 -928.33649 4091.9115 Apr 3 43248.71 -815.51957 -2439.7079 May 3 42624.95 -808.11344 -5611.9507 Jun 3 40613.56 -852.72343 4339.4372 Jul 3 40230.32 -835.77685 -2382.3156 Aug 3 39809.43 -820.91318 -7064.4325 Sep 3 37823.45 -862.70248 1577.5476 Oct 3 34176.97 -962.22111 754.0307 Nov 3 32666.07 -981.51794 341.9329 Dec 3 31032.09 -1003.93296 -22412.0922 Jan 4 32984.10 -902.39924 35921.8958 Feb 4 39983.86 -616.06846 -427.8623 Mar 4 45157.93 -392.72505 5511.0747 Apr 4 42311.73 -490.47729 -5879.7251 May 4 44006.66 -403.28391 -3115.6559 Jun 4 44400.00 -371.93358 4028.0001 Jul 4 41161.64 -483.10780 -4939.6356 Aug 4 39642.93 -522.88285 -6217.9286 Sep 4 37375.98 -589.43551 2025.0237 Oct 4 36463.28 -601.68613 1503.7150 Nov 4 35068.54 -631.49031 -267.5385 Dec 4 36021.82 -572.24845 -23364.8225 Jan 5 36770.32 -522.56495 32345.6777 Feb 5 40859.22 -345.47425 659.7758 Mar 5 43622.30 -223.19962 7698.7008 Apr 5 44983.39 -159.91861 -6454.3904 May 5 44967.25 -154.15520 -3420.2467 Jun 5 45831.51 -113.51088 6241.4933 Jul 5 44139.15 -176.06565 -5738.1462 Aug 5 44978.85 -136.09852 -4080.8527 Sep 5 41233.44 -277.27346 -794.4445 Oct 5 40055.35 -312.32399 1832.6454 Nov 5 39270.76 -330.62460 -1372.7631 Dec 5 35841.29 -450.58688 -27070.2878 Jan 6 36672.11 -400.76893 31511.8942 Feb 6 44815.72 -65.52916 5714.2762 Mar 6 43337.86 -121.51117 3883.1355 Apr 6 45917.02 -13.69395 -4161.0213 May 6 48172.09 77.10274 -2539.0934 Jun 6 44700.78 -64.70257 3437.2198 Jul 6 44509.82 -69.73045 -5023.8202 Aug 6 42869.99 -131.99883 -3528.9921 Sep 6 41722.61 -172.11098 -605.6097 Oct 6 40069.89 -230.41092 1559.1149 Nov 6 34269.63 -449.27212 -4547.6274 Dec 6 34282.77 -431.10364 -27228.7659 Jan 7 31240.33 -533.95173 25435.6698 Feb 7 29449.57 -583.66533 5420.4252 Mar 7 30850.70 -504.79076 4266.3001 Apr 7 33086.30 -395.52158 -2917.3045 May 7 32536.33 -401.68869 -1600.3263 Jun 7 32036.17 -405.61834 3662.8286 Jul 7 34772.35 -280.44592 -1544.3506 Aug 7 32637.86 -354.14744 -4904.8584 Sep 7 32335.34 -352.09974 1330.6595 Oct 7 31615.39 -366.66175 3813.6080 Nov 7 31487.87 -357.20505 -4049.8742 Dec 7 32772.50 -292.27651 -24602.5042 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -3.90585274 -0.46046973 -0.23035056 0.93621753 1.48395394 2 4.23373201 2.08248194 0.97174803 0.46253895 -0.52466781 0.67840998 3 1.33923181 0.87494033 1.12109339 0.67892018 0.04625895 -0.29220628 4 0.71973099 1.90998488 1.38083715 -0.58464739 0.52493748 0.19241862 5 0.31943842 1.11108776 0.74391578 0.37859179 0.03449799 0.24535207 6 0.30898496 2.05574244 -0.33849303 0.64648432 0.54430847 -0.85369748 7 -0.62864488 -0.30214869 0.47603693 0.65667055 -0.03705952 -0.02367048 Jul Aug Sep Oct Nov Dec 1 -0.54879341 -1.04897289 0.52771382 0.01788063 -0.81199606 -3.25536398 2 0.09478771 -0.58420445 0.06420114 0.17436529 -0.35262754 -0.65670801 3 0.11405057 0.10069784 -0.28249844 -0.67409357 -0.13269919 -0.15812512 4 -0.69310399 -0.25017328 -0.42073470 -0.07787826 -0.19102512 0.38254862 5 -0.38081885 0.24477785 -0.86839456 -0.21649620 -0.11354102 -0.74639857 6 -0.03040679 -0.37782701 -0.24398565 -0.35546869 -1.33788775 0.11123659 7 0.75586433 -0.44577294 0.01239856 -0.08827912 0.05741827 0.39464865 > 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 > postscript(file="/var/www/rcomp/tmp/15r651322232262.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') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/rcomp/tmp/2of9v1322232262.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/rcomp/tmp/3iaaf1322232262.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/rcomp/tmp/444f51322232262.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/rcomp/tmp/5r4vv1322232262.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/61sum1322232262.tab") > > try(system("convert tmp/15r651322232262.ps tmp/15r651322232262.png",intern=TRUE)) character(0) > try(system("convert tmp/2of9v1322232262.ps tmp/2of9v1322232262.png",intern=TRUE)) character(0) > try(system("convert tmp/3iaaf1322232262.ps tmp/3iaaf1322232262.png",intern=TRUE)) character(0) > try(system("convert tmp/444f51322232262.ps tmp/444f51322232262.png",intern=TRUE)) character(0) > try(system("convert tmp/5r4vv1322232262.ps tmp/5r4vv1322232262.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.708 0.340 3.051