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(31,36,24,22,17,8,12,5,6,5,8,15,16,17,23,24,27,31,40,47,43,60,64,65,65,55,57,57,57,65,69,70,71,71,73,68,65,57,41,21,21,17,9,11,6,-2,0,5,3,7,4,8,9,14,12,12,7,15,14,19) > 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 2.723508 6.753002 0.000000 8.109394 > m$fitted level slope sea Jan 1 31.0000000 0.0000000 0.00000000 Feb 1 34.2012639 1.5260413 0.16838869 Mar 1 26.8485585 -3.5133126 0.11316279 Apr 1 22.2069200 -4.1431841 0.11912547 May 1 17.1397873 -4.6541243 0.12416545 Jun 1 8.8948574 -6.6420792 0.13548952 Jul 1 9.7252169 -2.4961110 0.12593439 Aug 1 5.3999322 -3.5121160 0.12649567 Sep 1 4.9824967 -1.7925727 0.12662891 Oct 1 4.4969591 -1.0662848 0.12676316 Nov 1 6.8800044 0.8503633 0.12701060 Dec 1 13.2763091 3.9318938 0.12719838 Jan 2 17.3180027 3.9903331 -1.34828131 Feb 2 17.6638208 2.1767860 0.03150150 Mar 2 22.2644411 3.5070914 0.04103159 Apr 2 24.3588955 2.7198492 0.04491652 May 2 26.9823203 2.6663164 0.04520922 Jun 2 30.6650744 3.2301748 0.04338616 Jul 2 38.6056212 5.8452611 0.03988816 Aug 2 46.3996310 6.9278022 0.03952378 Sep 2 45.2778047 2.4551970 0.03937898 Oct 2 57.2270202 7.7305491 0.03991114 Nov 2 64.1830731 7.3002074 0.03987987 Dec 2 66.4182512 4.4859298 0.03978134 Jan 3 67.2936390 2.5221405 -1.27618606 Feb 3 57.7685263 -3.7020818 -0.08620090 Mar 3 56.4080977 -2.4114861 -0.08012502 Apr 3 56.3967701 -1.0752651 -0.08449456 May 3 56.6942455 -0.3129449 -0.08725721 Jun 3 63.1511804 3.4445820 -0.09530559 Jul 3 68.5376607 4.5229673 -0.09626109 Aug 3 70.7585537 3.2440644 -0.09597595 Sep 3 71.7456847 1.9900553 -0.09600285 Oct 3 71.6861027 0.8511962 -0.09607895 Nov 3 72.9711709 1.0922699 -0.09606735 Dec 3 69.4299322 -1.4822253 -0.09612705 Jan 4 65.1639987 -3.0071505 0.62606399 Feb 4 58.0814746 -5.1519430 -0.10841551 Mar 4 43.7597989 -10.2162200 -0.12617042 Apr 4 23.8840364 -15.5912062 -0.11302784 May 4 18.2682003 -10.0508800 -0.12804273 Jun 4 15.1438850 -6.2053703 -0.13420099 Jul 4 9.0905230 -6.1209450 -0.13425692 Aug 4 9.3099188 -2.5984891 -0.13484404 Sep 4 6.2637243 -2.8472447 -0.13484803 Oct 4 -0.6845081 -5.1258930 -0.13496186 Nov 4 -1.1940156 -2.5608933 -0.13486956 Dec 4 3.1477927 1.2744108 -0.13480306 Jan 5 2.3487737 0.1345321 1.24179229 Feb 5 6.1615882 2.0917271 -0.07561891 Mar 5 5.0119879 0.2994119 -0.08062464 Apr 5 7.4647948 1.4972956 -0.08296359 May 5 9.0561022 1.5495139 -0.08307661 Jun 5 13.3081551 3.0501541 -0.08499552 Jul 5 13.0390856 1.2065133 -0.08402039 Aug 5 12.5670415 0.2739508 -0.08389628 Sep 5 8.3707329 -2.2098256 -0.08392808 Oct 5 13.0893349 1.6398225 -0.08377452 Nov 5 14.2280361 1.3613859 -0.08378252 Dec 5 18.3027151 2.8689527 -0.08376165 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 1.00260733 -2.08436027 -0.24131841 -0.19621664 -0.76531448 2 0.02670341 -0.61637525 0.51757860 -0.30042917 -0.02048907 0.21669638 3 -0.82817301 -2.22237612 0.50024784 0.51137086 0.29230308 1.44469497 4 -0.62465970 -0.78227177 -1.95935121 -2.05987193 2.12629838 1.47884709 5 -0.45995557 0.72236582 -0.69268181 0.45945021 0.02005149 0.57716875 Jul Aug Sep Oct Nov Dec 1 1.59591424 -0.39099788 0.66170696 0.27948583 0.73755416 1.18581809 2 1.00618771 0.41657915 -1.72112207 2.03002798 -0.16560164 -1.08297500 3 0.41494137 -0.49214174 -0.48256070 -0.43824896 0.09276864 -0.99070330 4 0.03248592 1.35549493 -0.09572474 -0.87685617 0.98704900 1.47588095 5 -0.70942207 -0.35886420 -0.95579305 1.48139951 -0.10714646 0.58013372 > 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/1xaou1324540528.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/2ik0f1324540528.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/3l6mh1324540528.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/4fdre1324540528.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/5ugqn1324540529.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/69dtx1324540529.tab") > > try(system("convert tmp/1xaou1324540528.ps tmp/1xaou1324540528.png",intern=TRUE)) character(0) > try(system("convert tmp/2ik0f1324540528.ps tmp/2ik0f1324540528.png",intern=TRUE)) character(0) > try(system("convert tmp/3l6mh1324540528.ps tmp/3l6mh1324540528.png",intern=TRUE)) character(0) > try(system("convert tmp/4fdre1324540528.ps tmp/4fdre1324540528.png",intern=TRUE)) character(0) > try(system("convert tmp/5ugqn1324540529.ps tmp/5ugqn1324540529.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.611 0.252 1.889