R version 2.12.0 (2010-10-15) 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(274,291,280,258,252,251,224,225,234,233,229,208,224,226,223,205,201,202,183,188,200,206,211,201,299,244,251,241,244,252,234,246,265,277,287,275,320,338,342,322,323,343,315,334,359,362,378,345,422,430,443,431,425,432,387,396,411,421,424,410,464,486,490,459,454,446,406,412,428,429,425,396,429,439,424,379,370,353,322,322,338,348,350,312,358,378,352,312,310,292,276,269,286,292,288,255,304,299,293,275,272,264,234,231,263,264,264,245,297,317,318,315,312,310,306,313,350,354,371,357,419,425,424,399,393,378,371,364,384,377,383,352) > 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 74.186873 2.836064 30.129639 0.000000 > m$fitted level slope sea Jan 1 274.0000 0.00000000 0.0000000 Feb 1 287.4328 0.80743710 3.5671971 Mar 1 283.7453 0.43963815 -3.7453483 Apr 1 266.9070 -1.19695134 -8.9070175 May 1 254.3923 -2.44606383 -2.3922822 Jun 1 249.7035 -2.72810673 1.2965197 Jul 1 231.6478 -4.84625595 -7.6477582 Aug 1 223.5402 -5.32627872 1.4598201 Sep 1 228.1077 -3.80851184 5.8922605 Oct 1 231.2774 -2.70834978 1.7225794 Nov 1 229.7200 -2.52354618 -0.7200473 Dec 1 214.6631 -4.56080315 -6.6630722 Jan 2 219.4092 -3.11721113 4.5908152 Feb 2 219.0583 -2.66465879 6.9416713 Mar 2 219.5812 -2.14545684 3.4187816 Apr 2 213.6092 -2.76877596 -8.6092396 May 2 205.0964 -3.70934796 -4.0964062 Jun 2 196.6574 -4.48536628 5.3425933 Jul 2 190.1606 -4.81551231 -7.1605822 Aug 2 188.4610 -4.30388720 -0.4609625 Sep 2 192.6109 -2.91496393 7.3890891 Oct 2 199.7381 -1.26450321 6.2619065 Nov 2 205.5581 -0.10109500 5.4418816 Dec 2 210.1190 0.66212604 -9.1190210 Jan 3 263.0930 9.21228491 35.9070339 Feb 3 255.7177 6.48889705 -11.7177095 Mar 3 250.3994 4.56140493 0.6006288 Apr 3 249.3774 3.65242134 -8.3773842 May 3 248.4148 2.89875635 -4.4147684 Jun 3 247.6053 2.29168312 4.3946920 Jul 3 245.0763 1.50221225 -11.0763107 Aug 3 248.2654 1.77837748 -2.2653936 Sep 3 257.3102 2.96771594 7.6898151 Oct 3 269.6296 4.49774083 7.3704199 Nov 3 284.0101 6.11282970 2.9899357 Dec 3 299.4558 7.63613678 -24.4557899 Jan 4 289.1416 4.70341637 30.8584182 Feb 4 323.8546 9.61322925 14.1453941 Mar 4 339.6681 10.62467782 2.3319184 Apr 4 336.9402 8.45176748 -14.9402131 May 4 331.6971 6.22065630 -8.6970932 Jun 4 335.8273 5.87942472 7.1727303 Jul 4 332.3251 4.34695008 -17.3251311 Aug 4 337.6213 4.50198218 -3.6213344 Sep 4 350.6161 5.88838358 8.3839441 Oct 4 358.6892 6.24479226 3.3108244 Nov 4 372.7143 7.51294015 5.2856798 Dec 4 371.8239 6.14364526 -26.8239350 Jan 5 393.4168 8.66371794 28.5832015 Feb 5 412.8621 10.42355689 17.1379327 Mar 5 430.9692 11.67564183 12.0308455 Apr 5 441.7980 11.53789623 -10.7980083 May 5 439.9113 9.35358841 -14.9112718 Jun 5 430.6212 6.31614293 1.3788339 Jul 5 415.6689 2.84837781 -28.6688908 Aug 5 406.9146 0.95668013 -10.9145550 Sep 5 404.0788 0.33869326 6.9211929 Oct 5 413.4171 1.80412465 7.5829368 Nov 5 416.8710 2.07260431 7.1290210 Dec 5 434.9152 4.67203086 -24.9152223 Jan 6 442.0148 5.06746385 21.9851990 Feb 6 463.2367 7.69989746 22.7632781 Mar 6 476.1365 8.54635893 13.8634824 Apr 6 471.4657 6.39752660 -12.4657340 May 6 465.1662 4.33338766 -11.1662326 Jun 6 446.7876 0.63797584 -0.7875829 Jul 6 434.2435 -1.50831680 -28.2435338 Aug 6 424.8427 -2.79340247 -12.8427198 Sep 6 422.9654 -2.64431215 5.0345674 Oct 6 422.3454 -2.31510746 6.6546486 Nov 6 423.1876 -1.80183391 1.8124107 Dec 6 423.3474 -1.48286613 -27.3474043 Jan 7 417.5086 -2.19163193 11.4914292 Feb 7 416.2936 -2.03268364 22.7064148 Mar 7 407.7340 -3.09418841 16.2660116 Apr 7 391.2014 -5.27777093 -12.2013966 May 7 375.2168 -7.01722915 -5.2168091 Jun 7 354.5255 -9.24001475 -1.5254533 Jul 7 345.3545 -9.22877503 -23.3545441 Aug 7 335.8528 -9.27317024 -13.8528289 Sep 7 331.4322 -8.48422172 6.5678400 Oct 7 335.9720 -6.36789757 12.0279593 Nov 7 342.7929 -4.22543511 7.2070779 Dec 7 339.3008 -4.10628222 -27.3008061 Jan 8 342.1397 -2.97722127 15.8603089 Feb 8 347.2469 -1.66274162 30.7530692 Mar 8 335.7826 -3.25563441 16.2173842 Apr 8 322.6734 -4.85593989 -10.6733595 May 8 311.7229 -5.84556709 -1.7229370 Jun 8 297.5257 -7.20222953 -5.5257126 Jul 8 295.4864 -6.36318177 -19.4864113 Aug 8 287.8280 -6.57369123 -18.8279869 Sep 8 284.1745 -6.09923888 1.8254515 Oct 8 283.0586 -5.28997540 8.9413784 Nov 8 280.9509 -4.77330003 7.0490856 Dec 8 282.6010 -3.73017300 -27.6010410 Jan 9 286.5578 -2.48138197 17.4422380 Feb 9 271.6011 -4.50830919 27.3989073 Mar 9 268.9958 -4.19922730 24.0042270 Apr 9 275.9126 -2.39463259 -0.9125609 May 9 272.6819 -2.53033272 -0.6819144 Jun 9 270.7621 -2.43120801 -6.7621057 Jul 9 258.9310 -3.95789160 -24.9309856 Aug 9 252.0021 -4.44045470 -21.0021098 Sep 9 257.3780 -2.84645492 5.6220301 Oct 9 257.0079 -2.44448190 6.9921274 Nov 9 258.1932 -1.85537206 5.8068327 Dec 9 267.8467 0.01282211 -22.8466996 Jan 10 273.1528 0.87230921 23.8471958 Feb 10 285.3566 2.71239512 31.6433768 Mar 10 295.4391 3.90889632 22.5609320 Apr 10 309.8869 5.61925756 5.1131148 May 10 313.9019 5.35893240 -1.9018753 Jun 10 314.8123 4.63694754 -4.8123368 Jul 10 325.9463 5.69167078 -19.9462634 Aug 10 336.1346 6.42170959 -23.1346354 Sep 10 344.7811 6.78283528 5.2188503 Oct 10 350.4066 6.59502494 3.5933853 Nov 10 365.2260 7.92945751 5.7740449 Dec 10 380.0930 9.05525186 -23.0929815 Jan 11 396.8362 10.30309435 22.1637569 Feb 11 401.6139 9.40619840 23.3861305 Mar 11 406.9325 8.74281925 17.0675100 Apr 11 400.9826 6.35892414 -1.9825990 May 11 397.6390 4.78483647 -4.6389908 Jun 11 391.3806 2.99302228 -13.3806362 Jul 11 392.6438 2.71229264 -21.6438092 Aug 11 390.7668 1.96748467 -26.7668317 Sep 11 384.1260 0.57060851 -0.1259597 Oct 11 379.8159 -0.22121041 -2.8159139 Nov 11 380.2609 -0.11313317 2.7390662 Dec 11 378.8922 -0.31683073 -26.8922478 > m$resid Jan Feb Mar Apr May 1 0.000000000 1.032225570 -0.438522197 -1.857322676 -1.194474071 2 0.992243088 0.268319342 0.296268193 -0.372201922 -0.564076859 3 5.215004343 -1.621466626 -1.124759626 -0.538857594 -0.450204887 4 -1.761123743 2.919569907 0.595440851 -1.286901295 -1.330124582 5 1.505593944 1.045492830 0.739745722 -0.081599118 -1.300242532 6 0.235698147 1.563276784 0.500983373 -1.273551867 -1.227570586 7 -0.421928862 0.094379005 -0.628899395 -1.294714248 -1.033907156 8 0.671618811 0.780476391 -0.944344467 -0.949206633 -0.588020341 9 0.742469641 -1.203502989 0.183321079 1.070672801 -0.080613544 10 0.510833244 1.092584505 0.709888545 1.014976752 -0.154626266 11 0.741469772 -0.532558107 -0.393675612 -1.414895988 -0.934879559 Jun Jul Aug Sep Oct 1 -0.232141474 -1.566273576 -0.330232371 0.995210649 0.698587150 2 -0.461945286 -0.196056309 0.303864749 0.825137535 0.980427219 3 -0.361184390 -0.468151748 0.163699748 0.705309889 0.907455030 4 -0.203264039 -0.910102604 0.091925819 0.821773221 0.211306295 5 -1.809349561 -2.061225885 -1.122148065 -0.366267065 0.868849581 6 -2.200286386 -1.276125783 -0.762571112 0.088376643 0.195222560 7 -1.322747714 0.006682629 -0.026350635 0.467796383 1.255302480 8 -0.806929091 0.498776173 -0.124967641 0.281403293 0.480123924 9 0.058934253 -0.907359691 -0.286500640 0.945671889 0.238530683 10 -0.429115875 0.626737156 0.433457062 0.214292141 -0.111464266 11 -1.064714667 -0.166787814 -0.442244267 -0.829050722 -0.469998669 Nov Dec 1 0.114834023 -1.247498805 2 0.690752575 0.454938346 3 0.958567168 0.908649806 4 0.753097178 -0.816386037 5 0.159488989 1.548567523 6 0.304930803 0.189880642 7 1.272762482 0.070888704 8 0.306913637 0.620302594 9 0.349911206 1.110527642 10 0.792549758 0.669029798 11 0.064185839 -0.121026723 > 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/1et121322754974.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/2ivcb1322754974.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/33m181322754974.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/4xiha1322754974.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/5o60u1322754974.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/69id91322754974.tab") > > try(system("convert tmp/1et121322754974.ps tmp/1et121322754974.png",intern=TRUE)) character(0) > try(system("convert tmp/2ivcb1322754974.ps tmp/2ivcb1322754974.png",intern=TRUE)) character(0) > try(system("convert tmp/33m181322754974.ps tmp/33m181322754974.png",intern=TRUE)) character(0) > try(system("convert tmp/4xiha1322754974.ps tmp/4xiha1322754974.png",intern=TRUE)) character(0) > try(system("convert tmp/5o60u1322754974.ps tmp/5o60u1322754974.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.790 0.160 1.959