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(37,30,47,35,30,43,82,40,47,19,52,136,80,42,54,66,81,63,137,72,107,58,36,52,79,77,54,84,48,96,83,66,61,53,30,74,69,59,42,65,70,100,63,105,82,81,75,102,121,98,76,77,63,37,35,23,40,29,37,51,20,28,13,22,25,13,16,13,16,17,9,17,25,14,8,7,10,7,10,3) > 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 66.9018 0.0000 0.0000 308.9231 > m$fitted level slope sea Jan 1 37.000000 0.0000000000 0.0000000000 Feb 1 34.067108 -0.3156477376 -0.3156477364 Mar 1 39.275426 0.1130180117 0.1130180117 Apr 1 37.598825 0.0060909456 0.0060909456 May 1 34.631992 -0.1387018410 -0.1387018410 Jun 1 37.829331 0.0002856515 0.0002856515 Jul 1 54.734562 0.6263475592 0.6263475592 Aug 1 49.275539 0.4201592630 0.4201592630 Sep 1 48.514012 0.3827646169 0.3827646169 Oct 1 37.441837 0.0393364740 0.0393364740 Nov 1 42.949318 0.1961707162 0.1961707162 Dec 1 78.093425 1.1614093447 1.1614093447 Jan 2 82.405129 0.8036404206 -8.8400446305 Feb 2 65.051615 0.0542732849 0.0542732963 Mar 2 60.629143 -0.0882433309 -0.0882433309 Apr 2 62.690874 -0.0332820723 -0.0332820723 May 2 69.659982 0.1173069798 0.1173069798 Jun 2 67.173087 0.0681274579 0.0681274579 Jul 2 93.436458 0.5179943263 0.5179943263 Aug 2 85.531368 0.3829471369 0.3829471369 Sep 2 93.660838 0.5010775950 0.5010775950 Oct 2 80.458127 0.2997695641 0.2997695641 Nov 2 63.928195 0.0595839554 0.0595839554 Dec 2 59.490069 -0.0031260372 -0.0031260372 Jan 3 65.169006 -0.2969582114 3.2665403118 Feb 3 69.932718 -0.1633159521 -0.1633159407 Mar 3 63.696620 -0.2869526065 -0.2869526065 Apr 3 71.361935 -0.1564286817 -0.1564286817 May 3 62.518547 -0.2770175998 -0.2770175998 Jun 3 74.997921 -0.1206463388 -0.1206463388 Jul 3 77.957356 -0.0861220409 -0.0861220409 Aug 3 73.475655 -0.1323633291 -0.1323633291 Sep 3 68.793778 -0.1781178077 -0.1781178077 Oct 3 62.867823 -0.2340777401 -0.2340777401 Nov 3 50.575605 -0.3486754103 -0.3486754103 Dec 3 59.201718 -0.2649593900 -0.2649593900 Jan 4 61.189328 -0.3403712425 3.7440836556 Feb 4 60.250705 -0.3519874725 -0.3519874609 Mar 4 53.165463 -0.4527542024 -0.4527542024 Apr 4 57.526350 -0.3945691737 -0.3945691737 May 4 62.104415 -0.3436056393 -0.3436056393 Jun 4 76.168024 -0.2128682335 -0.2128682335 Jul 4 71.207759 -0.2523643977 -0.2523643977 Aug 4 83.716784 -0.1524647526 -0.1524647526 Sep 4 83.039400 -0.1564029468 -0.1564029468 Oct 4 82.241369 -0.1610743407 -0.1610743407 Nov 4 79.509801 -0.1793862102 -0.1793862102 Dec 4 87.817101 -0.1199088442 -0.1199088442 Jan 5 98.177682 -0.3793096880 4.1724065508 Feb 5 98.025275 -0.3758261392 -0.3758261225 Mar 5 89.551325 -0.4716364865 -0.4716364865 Apr 5 84.714260 -0.5133973338 -0.5133973338 May 5 76.467998 -0.5762115213 -0.5762115213 Jun 5 61.615390 -0.6790513572 -0.6790513572 Jul 5 51.542148 -0.7411892469 -0.7411892469 Aug 5 40.746674 -0.8038613744 -0.8038613744 Sep 5 40.262546 -0.8019484705 -0.8019484705 Oct 5 35.875405 -0.8227936279 -0.8227936279 Nov 5 36.080629 -0.8169368008 -0.8169368008 Dec 5 41.405660 -0.7824643147 -0.7824643147 Jan 6 30.815377 -0.5899092288 6.4890015092 Feb 6 29.592893 -0.5979399288 -0.5979399206 Mar 6 23.164318 -0.6549773984 -0.6549773984 Apr 6 22.563943 -0.6545452075 -0.6545452075 May 6 23.306534 -0.6451458442 -0.6451458442 Jun 6 19.306722 -0.6651804371 -0.6651804371 Jul 6 17.907274 -0.6692111759 -0.6692111759 Aug 6 15.913266 -0.6760714474 -0.6760714474 Sep 6 15.771051 -0.6734154697 -0.6734154697 Oct 6 16.053013 -0.6687918733 -0.6687918733 Nov 6 13.264738 -0.6788522571 -0.6788522571 Dec 6 14.474299 -0.6700135525 -0.6700135525 Jan 7 15.108449 -0.6912306510 7.6035371538 Feb 7 14.522471 -0.6900915932 -0.6900915862 Mar 7 11.892175 -0.7062698249 -0.7062698249 Apr 7 9.884120 -0.7150572274 -0.7150572274 May 7 9.744991 -0.7117506403 -0.7117506403 Jun 7 8.542051 -0.7142558711 -0.7142558711 Jul 7 8.899340 -0.7092287084 -0.7092287084 Aug 7 6.527864 -0.7165904301 -0.7165904301 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.291557374 0.578107230 -0.194826690 -0.332719357 0.380389932 2 0.530695411 -1.772657942 -0.486103541 0.244989300 0.817044166 -0.307662325 3 0.816165035 0.539272848 -0.689459783 0.929872417 -1.030813041 1.525512492 4 0.306763104 -0.066247760 -0.780050188 0.569637434 0.594750080 1.732936482 5 1.389563176 0.025675290 -0.949061241 -0.520171727 -0.929148217 -1.723024047 6 -1.279585967 -0.072578703 -0.688443433 0.006535501 0.168379160 -0.405776124 7 0.168326296 0.012193350 -0.230287014 -0.156306744 0.069556556 -0.059506994 Jul Aug Sep Oct Nov Dec 1 1.951635376 -0.708355609 -0.138315763 -1.345904435 0.644219298 4.125546633 2 3.116206370 -1.006076266 0.927553151 -1.643443696 -2.020427971 -0.540347391 3 0.369942032 -0.529281928 -0.548659557 -0.693843430 -1.456498568 1.084522091 4 -0.572834758 1.542663918 -0.063525685 -0.077703577 -0.311435531 1.028534302 5 -1.136605812 -1.218230871 0.038774415 -0.435015167 0.124779244 0.745672599 6 -0.088997805 -0.160764798 0.064830216 0.116070036 -0.257571299 0.229532516 7 0.130042001 -0.201933370 > 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/wessaorg/rcomp/tmp/1h0l91355145043.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/25e5j1355145043.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/3yms01355145043.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/41voq1355145043.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/wessaorg/rcomp/tmp/5ic9y1355145043.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/6lr801355145043.tab") > > try(system("convert tmp/1h0l91355145043.ps tmp/1h0l91355145043.png",intern=TRUE)) character(0) > try(system("convert tmp/25e5j1355145043.ps tmp/25e5j1355145043.png",intern=TRUE)) character(0) > try(system("convert tmp/3yms01355145043.ps tmp/3yms01355145043.png",intern=TRUE)) character(0) > try(system("convert tmp/41voq1355145043.ps tmp/41voq1355145043.png",intern=TRUE)) character(0) > try(system("convert tmp/5ic9y1355145043.ps tmp/5ic9y1355145043.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.135 0.589 3.717