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(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.187072 2.836062 30.129542 0.000000 > m$fitted level slope sea Jan 1 274.0000 0.00000000 0.0000000 Feb 1 287.4328 0.80743789 3.5671837 Mar 1 283.7453 0.43963794 -3.7453423 Apr 1 266.9070 -1.19694932 -8.9069913 May 1 254.3923 -2.44605809 -2.3922645 Jun 1 249.7035 -2.72809905 1.2965211 Jul 1 231.6477 -4.84624695 -7.6477438 Aug 1 223.5402 -5.32626835 1.4598222 Sep 1 228.1078 -3.80850306 5.8922457 Oct 1 231.2774 -2.70834483 1.7225687 Nov 1 229.7200 -2.52354400 -0.7200474 Dec 1 214.6631 -4.56080120 -6.6630546 Jan 2 219.4092 -3.11720808 4.5908085 Feb 2 219.0583 -2.66465515 6.9416542 Mar 2 219.5812 -2.14545480 3.4187665 Apr 2 213.6092 -2.76877686 -8.6092339 May 2 205.0964 -3.70934852 -4.0963943 Jun 2 196.6574 -4.48536450 5.3425976 Jul 2 190.1606 -4.81551122 -7.1605729 Aug 2 188.4610 -4.30388641 -0.4609564 Sep 2 192.6109 -2.91496239 7.3890791 Oct 2 199.7381 -1.26450245 6.2618870 Nov 2 205.5581 -0.10109698 5.4418691 Dec 2 210.1190 0.66212049 -9.1190156 Jan 3 263.0930 9.21227755 35.9069891 Feb 3 255.7177 6.48888516 -11.7176984 Mar 3 250.3994 4.56139787 0.6006372 Apr 3 249.3774 3.65241728 -8.3773802 May 3 248.4148 2.89875341 -4.4147616 Jun 3 247.6053 2.29168111 4.3947008 Jul 3 245.0763 1.50221075 -11.0762970 Aug 3 248.2654 1.77837713 -2.2653875 Sep 3 257.3102 2.96771635 7.6898077 Oct 3 269.6296 4.49774005 7.3704078 Nov 3 284.0101 6.11282539 2.9899336 Dec 3 299.4558 7.63612991 -24.4557836 Jan 4 289.1416 4.70341540 30.8584178 Feb 4 323.8546 9.61322929 14.1453509 Mar 4 339.6681 10.62467236 2.3319011 Apr 4 336.9402 8.45176113 -14.9402028 May 4 331.6971 6.22065260 -8.6970760 Jun 4 335.8273 5.87942321 7.1727404 Jul 4 332.3251 4.34694876 -17.3251089 Aug 4 337.6213 4.50198239 -3.6213214 Sep 4 350.6161 5.88838436 8.3839427 Oct 4 358.6892 6.24479230 3.3108241 Nov 4 372.7143 7.51293966 5.2856731 Dec 4 371.8239 6.14364425 -26.8239276 Jan 5 393.4168 8.66371685 28.5831912 Feb 5 412.8621 10.42355723 17.1379013 Mar 5 430.9692 11.67564091 12.0308119 Apr 5 441.7980 11.53789199 -10.7980194 May 5 439.9113 9.35358306 -14.9112553 Jun 5 430.6211 6.31614018 1.3788627 Jul 5 415.6688 2.84837784 -28.6688499 Aug 5 406.9145 0.95668417 -10.9145242 Sep 5 404.0788 0.33870041 6.9212052 Oct 5 413.4171 1.80413186 7.5829305 Nov 5 416.8710 2.07260954 7.1290176 Dec 5 434.9152 4.67203282 -24.9152306 Jan 6 442.0148 5.06746456 21.9851934 Feb 6 463.2368 7.69989952 22.7632435 Mar 6 476.1366 8.54635986 13.8634462 Apr 6 471.4657 6.39752640 -12.4657477 May 6 465.1662 4.33338766 -11.1662317 Jun 6 446.7876 0.63797630 -0.7875568 Jul 6 434.2435 -1.50831502 -28.2434999 Aug 6 424.8427 -2.79339933 -12.8426863 Sep 6 422.9654 -2.64430767 5.0345884 Oct 6 422.3453 -2.31510234 6.6546587 Nov 6 423.1876 -1.80182984 1.8124182 Dec 6 423.3474 -1.48286214 -27.3474028 Jan 7 417.5086 -2.19162758 11.4914299 Feb 7 416.2936 -2.03267577 22.7063883 Mar 7 407.7340 -3.09417911 16.2659765 Apr 7 391.2014 -5.27776253 -12.2014188 May 7 375.2168 -7.01722188 -5.2168195 Jun 7 354.5254 -9.24001022 -1.5254373 Jul 7 345.3545 -9.22877210 -23.3545228 Aug 7 335.8528 -9.27316898 -13.8528000 Sep 7 331.4321 -8.48422038 6.5678611 Oct 7 335.9720 -6.36789644 12.0279644 Nov 7 342.7929 -4.22543535 7.2070739 Dec 7 339.3008 -4.10628308 -27.3008069 Jan 8 342.1397 -2.97722100 15.8602935 Feb 8 347.2470 -1.66273955 30.7530325 Mar 8 335.7826 -3.25563186 16.2173546 Apr 8 322.6734 -4.85593713 -10.6733808 May 8 311.7229 -5.84556552 -1.7229462 Jun 8 297.5257 -7.20223080 -5.5256956 Jul 8 295.4864 -6.36318428 -19.4863917 Aug 8 287.8280 -6.57369553 -18.8279523 Sep 8 284.1745 -6.09924253 1.8254828 Oct 8 283.0586 -5.28997711 8.9413953 Nov 8 280.9509 -4.77330005 7.0490898 Dec 8 282.6010 -3.73017255 -27.6010475 Jan 9 286.5578 -2.48138034 17.4422143 Feb 9 271.6011 -4.50830532 27.3988835 Mar 9 268.9958 -4.19922204 24.0041886 Apr 9 275.9126 -2.39463022 -0.9126007 May 9 272.6819 -2.53033519 -0.6819257 Jun 9 270.7621 -2.43121395 -6.7620938 Jul 9 258.9310 -3.95789815 -24.9309525 Aug 9 252.0021 -4.44046020 -21.0020729 Sep 9 257.3779 -2.84645921 5.6220529 Oct 9 257.0079 -2.44448546 6.9921469 Nov 9 258.1932 -1.85537383 5.8068403 Dec 9 267.8467 0.01282222 -22.8467161 Jan 10 273.1528 0.87231025 23.8471670 Feb 10 285.3567 2.71239481 31.6433413 Mar 10 295.4391 3.90889434 22.5608981 Apr 10 309.8869 5.61925432 5.1130775 May 10 313.9019 5.35892593 -1.9018874 Jun 10 314.8123 4.63693952 -4.8123274 Jul 10 325.9462 5.69166115 -19.9462442 Aug 10 336.1346 6.42169870 -23.1346042 Sep 10 344.7811 6.78282606 5.2188792 Oct 10 350.4066 6.59501788 3.5934116 Nov 10 365.2259 7.92945303 5.7740522 Dec 10 380.0930 9.05524944 -23.0929909 Jan 11 396.8363 10.30309275 22.1637354 Feb 11 401.6139 9.40619745 23.3861133 Mar 11 406.9325 8.74282000 17.0674884 Apr 11 400.9826 6.35892644 -1.9826128 May 11 397.6390 4.78483877 -4.6389949 Jun 11 391.3806 2.99302296 -13.3806189 Jul 11 392.6438 2.71229295 -21.6437883 Aug 11 390.7668 1.96748444 -26.7668020 Sep 11 384.1259 0.57060910 -0.1259243 Oct 11 379.8159 -0.22120814 -2.8158834 Nov 11 380.2609 -0.11312834 2.7390780 Dec 11 378.8923 -0.31682380 -26.8922524 > m$resid Jan Feb Mar Apr May 1 0.000000000 1.032224647 -0.438523977 -1.857323073 -1.194472244 2 0.992244843 0.268319744 0.296267431 -0.372203872 -0.564076893 3 5.215005876 -1.621469903 -1.124757252 -0.538856064 -0.450204397 4 -1.761120991 2.919571680 0.595437848 -1.286902390 -1.330123563 5 1.505594549 1.045494110 0.739745281 -0.081601123 -1.300243755 6 0.235697509 1.563278230 0.500982912 -1.273553076 -1.227570994 7 -0.421928825 0.094381139 -0.628898811 -1.294715324 -1.033908277 8 0.671619779 0.780477779 -0.944344573 -0.949206906 -0.588021296 9 0.742470660 -1.203502158 0.183321972 1.070671544 -0.080616453 10 0.510834012 1.092584160 0.709887851 1.014976430 -0.154628247 11 0.741470574 -0.532557942 -0.393674772 -1.414895656 -0.934879955 Jun Jul Aug Sep Oct 1 -0.232140171 -1.566274366 -0.330231751 0.995210477 0.698585263 2 -0.461944097 -0.196056808 0.303864705 0.825138348 0.980427182 3 -0.361183986 -0.468151635 0.163700496 0.705310653 0.907454706 4 -0.203262815 -0.910102868 0.091926766 0.821773917 0.211305947 5 -1.809348748 -2.061225078 -1.122146153 -0.366265383 0.868849988 6 -2.200287035 -1.276125527 -0.762570622 0.088377473 0.195223022 7 -1.322749905 0.006681688 -0.026351639 0.467796625 1.255302880 8 -0.806931129 0.498775644 -0.124968751 0.281403798 0.480125269 9 0.058932207 -0.907360433 -0.286500140 0.945673008 0.238531213 10 -0.429116978 0.626736468 0.433456485 0.214293228 -0.111463029 11 -1.064716086 -0.166788101 -0.442244776 -0.829050582 -0.469997868 Nov Dec 1 0.114832378 -1.247499724 2 0.690751244 0.454936412 3 0.958565479 0.908648658 4 0.753097174 -0.816386687 5 0.159487879 1.548566233 6 0.304930311 0.189880670 7 1.272762212 0.070888360 8 0.306914776 0.620303133 9 0.349912417 1.110529236 10 0.792551626 0.669031302 11 0.064187382 -0.121025519 > 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/rcomp/tmp/1y1j31322754877.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/266o71322754877.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/35jip1322754877.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/4v41s1322754877.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/www/rcomp/tmp/5ansj1322754877.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/6rcad1322754877.tab") > > try(system("convert tmp/1y1j31322754877.ps tmp/1y1j31322754877.png",intern=TRUE)) character(0) > try(system("convert tmp/266o71322754877.ps tmp/266o71322754877.png",intern=TRUE)) character(0) > try(system("convert tmp/35jip1322754877.ps tmp/35jip1322754877.png",intern=TRUE)) character(0) > try(system("convert tmp/4v41s1322754877.ps tmp/4v41s1322754877.png",intern=TRUE)) character(0) > try(system("convert tmp/5ansj1322754877.ps tmp/5ansj1322754877.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.996 0.364 3.350