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(13328,12873,14000,13477,14237,13674,13529,14058,12975,14326,14008,16193,14483,14011,15057,14884,15414,14440,14900,15074,14442,15307,14938,17193,15528,14765,15838,15723,16150,15486,15986,15983,15692,16490,15686,18897,16316,15636,17163,16534,16518,16375,16290,16352,15943,16362,16393,19051,16747,16320,17910,16961,17480,17049,16879,17473,16998,17307,17418,20169,17871,17226,19062,17804,19100,18522,18060,18869,18127,18871,18890,21263,19547,18450,20254,19240,20216,19420,19415,20018,18652,19978,19509,21971) > 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 15123.89433 51.98838 56600.92406 0.00000 > m$fitted level slope sea Jan 1 13328.00 0.000000 0.00000 Feb 1 13200.61 3.821531 -327.60962 Mar 1 13460.25 21.998825 539.74732 Apr 1 13535.99 26.277874 -58.98692 May 1 13789.38 42.162266 447.61933 Jun 1 13834.08 42.304893 -160.07859 Jul 1 13756.19 36.769596 -227.19425 Aug 1 13824.11 38.041400 233.89185 Sep 1 13582.50 27.030708 -607.49839 Oct 1 13729.75 31.841451 596.25301 Nov 1 13869.99 36.341648 138.00504 Dec 1 14657.75 68.755205 1535.24856 Jan 2 14940.97 73.377369 -457.96872 Feb 2 14973.31 72.310956 -962.30546 Mar 2 14910.20 66.740634 146.79578 Apr 2 14966.77 66.209601 -82.76658 May 2 14983.39 63.397112 430.61484 Jun 2 14841.90 51.759857 -401.90172 Jul 2 14877.57 50.878188 22.42749 Aug 2 14855.28 47.043456 218.71694 Sep 2 14996.14 51.771207 -554.13900 Oct 2 15086.87 53.673258 220.12551 Nov 2 15227.87 57.817220 -289.86819 Dec 2 15465.19 66.119214 1727.81145 Jan 3 15699.36 73.785741 -171.35755 Feb 3 15768.90 73.587164 -1003.89750 Mar 3 15770.23 70.023404 67.77441 Apr 3 15767.79 66.259050 -44.79442 May 3 15725.55 60.428797 424.44986 Jun 3 15780.16 60.111802 -294.15554 Jul 3 15877.23 62.117935 108.77082 Aug 3 15928.55 61.539010 54.45160 Sep 3 16100.19 67.347619 -408.19224 Oct 3 16281.84 73.282929 208.16354 Nov 3 16325.10 71.743770 -639.10177 Dec 3 16660.86 85.183806 2236.13709 Jan 4 16722.82 84.001348 -406.82119 Feb 4 16713.58 79.216834 -1077.57710 Mar 4 16831.78 81.242856 331.22430 Apr 4 16796.04 75.081041 -262.03614 May 4 16595.76 60.429532 -77.76334 Jun 4 16599.69 57.407494 -224.69355 Jul 4 16513.57 49.727151 -223.57404 Aug 4 16484.46 45.521400 -132.45876 Sep 4 16491.13 43.459685 -548.13228 Oct 4 16429.67 37.921581 -67.66523 Nov 4 16668.18 48.461550 -275.18041 Dec 4 16787.46 52.173455 2263.53751 Jan 5 16927.77 56.792743 -180.77222 Feb 5 17113.10 63.544379 -793.10420 Mar 5 17270.62 68.498012 639.37914 Apr 5 17247.87 63.669406 -286.86552 May 5 17325.07 64.387687 154.93180 Jun 5 17312.03 60.271036 -263.02794 Jul 5 17242.32 53.358112 -363.32340 Aug 5 17351.46 56.322624 121.53643 Sep 5 17478.67 60.084164 -480.67220 Oct 5 17568.38 61.653652 -261.38452 Nov 5 17693.24 64.996615 -275.23926 Dec 5 17860.16 70.382992 2308.84047 Jan 6 18029.55 75.614393 -158.54567 Feb 6 18110.44 75.893576 -884.44050 Mar 6 18234.47 78.441185 827.52957 Apr 6 18242.10 74.689289 -438.09946 May 6 18466.82 82.646018 633.17641 Jun 6 18650.02 87.981592 -128.02300 Jul 6 18689.53 85.408826 -629.52613 Aug 6 18776.60 85.496977 92.40371 Sep 6 18804.52 82.443510 -677.52118 Oct 6 18980.91 87.423691 -109.91368 Nov 6 19161.23 92.345370 -271.22699 Dec 6 19208.34 89.949539 2054.65943 Jan 7 19420.03 96.397176 126.97324 Feb 7 19505.23 95.804417 -1055.23376 Mar 7 19548.79 93.036167 705.20617 Apr 7 19679.91 95.054365 -439.90631 May 7 19740.47 93.225724 475.52972 Jun 7 19722.55 87.332085 -302.55446 Jul 7 19854.63 89.704696 -439.62632 Aug 7 19957.11 90.382535 60.88654 Sep 7 19856.40 80.250384 -1204.40466 Oct 7 19952.88 81.110615 25.11854 Nov 7 19964.53 77.428600 -455.52839 Dec 7 20006.57 75.553237 1964.42666 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.30910987 1.31766722 0.31190209 1.58143730 0.01924002 2 1.70441395 -0.32902091 -1.01366957 -0.07218915 -0.35218503 -1.49386935 3 1.26382724 -0.03188330 -0.53599545 -0.52967142 -0.78914862 -0.04257440 4 -0.17202310 -0.69027653 0.28759255 -0.85920712 -2.01772610 -0.41444463 5 0.64952498 0.94730095 0.69187122 -0.67085876 0.09939843 -0.56879139 6 0.72843195 0.03885563 0.35411268 -0.52066655 1.10276755 0.73907937 7 0.89524656 -0.08229761 -0.38420452 0.27996650 -0.25356375 -0.81706654 Jul Aug Sep Oct Nov Dec 1 -0.93687537 0.24393249 -2.18563592 0.93563574 0.83989279 5.79778367 2 -0.11989599 -0.55070544 0.70731937 0.29317482 0.65522058 1.34613445 3 0.27249806 -0.08004500 0.81698324 0.84680104 -0.22207404 1.95349655 4 -1.05580414 -0.58126280 -0.28653898 -0.77335593 1.47741619 0.52164513 5 -0.95569742 0.41050326 0.52176265 0.21801433 0.46490503 0.74974634 6 -0.35643002 0.01221849 -0.42345897 0.69094265 0.68305593 -0.33260226 7 0.32893993 0.09399092 -1.40517485 0.11931483 -0.51074560 -0.26016161 > 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/wessaorg/rcomp/tmp/10xkb1355399739.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/wessaorg/rcomp/tmp/2mr3g1355399739.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/39tpr1355399739.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/4ijpd1355399739.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/5d5mr1355399739.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/6v8st1355399739.tab") > > try(system("convert tmp/10xkb1355399739.ps tmp/10xkb1355399739.png",intern=TRUE)) character(0) > try(system("convert tmp/2mr3g1355399739.ps tmp/2mr3g1355399739.png",intern=TRUE)) character(0) > try(system("convert tmp/39tpr1355399739.ps tmp/39tpr1355399739.png",intern=TRUE)) character(0) > try(system("convert tmp/4ijpd1355399739.ps tmp/4ijpd1355399739.png",intern=TRUE)) character(0) > try(system("convert tmp/5d5mr1355399739.ps tmp/5d5mr1355399739.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.529 0.598 4.113