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(1.2999,1.3074,1.3242,1.3516,1.3511,1.3419,1.3716,1.3622,1.3896,1.4227,1.4684,1.457,1.4718,1.4748,1.5527,1.5751,1.5557,1.5553,1.577,1.4975,1.437,1.3322,1.2732,1.3449,1.3239,1.2785,1.305,1.319,1.365,1.4016,1.4088,1.4268,1.4562,1.4816,1.4914,1.4614,1.4272,1.3686,1.3569,1.3406,1.2565,1.2209,1.277,1.2894,1.3067,1.3898,1.3661,1.322,1.336,1.3649,1.3999,1.4442,1.4349,1.4388,1.4264,1.4343,1.377,1.3706,1.3556,1.3179,1.2905,1.3224,1.3201,1.3162,1.2789,1.2526,1.2288,1.24,1.2856) > 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 0.00140668 0.00000000 0.00000000 0.00000000 > m$fitted level slope sea Jan 1 1.299900 0.000000e+00 0.000000e+00 Feb 1 1.307010 3.900000e-04 3.900000e-04 Mar 1 1.323745 4.553785e-04 4.553785e-04 Apr 1 1.351038 5.623015e-04 5.623015e-04 May 1 1.350542 5.581027e-04 5.581027e-04 Jun 1 1.341380 5.196850e-04 5.196850e-04 Jul 1 1.370966 6.341176e-04 6.341176e-04 Aug 1 1.361605 5.949218e-04 5.949218e-04 Sep 1 1.388901 6.992217e-04 6.992217e-04 Oct 1 1.421875 8.248061e-04 8.248061e-04 Nov 1 1.467402 9.980694e-04 9.980694e-04 Dec 1 1.456050 9.503845e-04 9.503845e-04 Jan 2 1.474234 2.212598e-04 -2.433858e-03 Feb 2 1.474576 2.241818e-04 2.241818e-04 Mar 2 1.552335 3.651542e-04 3.651542e-04 Apr 2 1.574695 4.050724e-04 4.050724e-04 May 2 1.555331 3.692586e-04 3.692586e-04 Jun 2 1.554932 3.678700e-04 3.678700e-04 Jul 2 1.576594 4.063063e-04 4.063063e-04 Aug 2 1.497237 2.625899e-04 2.625899e-04 Sep 2 1.436846 1.535009e-04 1.535009e-04 Oct 2 1.332235 -3.458782e-05 -3.458782e-05 Nov 2 1.273340 -1.400716e-04 -1.400716e-04 Dec 2 1.344912 -1.178572e-05 -1.178572e-05 Jan 3 1.327525 3.295154e-04 -3.624669e-03 Feb 3 1.278930 -4.303529e-04 -4.303529e-04 Mar 3 1.305399 -3.987074e-04 -3.987074e-04 Apr 3 1.319382 -3.818075e-04 -3.818075e-04 May 3 1.365327 -3.274326e-04 -3.274326e-04 Jun 3 1.401884 -2.841920e-04 -2.841920e-04 Jul 3 1.409075 -2.754386e-04 -2.754386e-04 Aug 3 1.427054 -2.540888e-04 -2.540888e-04 Sep 3 1.456419 -2.194866e-04 -2.194866e-04 Oct 3 1.481790 -1.896270e-04 -1.896270e-04 Nov 3 1.491578 -1.779977e-04 -1.779977e-04 Dec 3 1.461613 -2.126744e-04 -2.126744e-04 Jan 4 1.429367 1.970438e-04 -2.167482e-03 Feb 4 1.369094 -4.943478e-04 -4.943477e-04 Mar 4 1.357404 -5.040834e-04 -5.040834e-04 Apr 4 1.341118 -5.177951e-04 -5.177951e-04 May 4 1.257090 -5.902862e-04 -5.902862e-04 Jun 4 1.221521 -6.206239e-04 -6.206239e-04 Jul 4 1.277572 -5.715152e-04 -5.715152e-04 Aug 4 1.289960 -5.602941e-04 -5.602941e-04 Sep 4 1.307245 -5.448574e-04 -5.448574e-04 Oct 4 1.390273 -4.726252e-04 -4.726252e-04 Nov 4 1.366593 -4.926661e-04 -4.926661e-04 Dec 4 1.322530 -5.302586e-04 -5.302586e-04 Jan 5 1.329396 -6.003903e-04 6.604293e-03 Feb 5 1.365171 -2.713104e-04 -2.713103e-04 Mar 5 1.400147 -2.470021e-04 -2.470021e-04 Apr 5 1.444416 -2.163223e-04 -2.163223e-04 May 5 1.435123 -2.225740e-04 -2.225740e-04 Jun 5 1.439020 -2.197387e-04 -2.197387e-04 Jul 5 1.426628 -2.281100e-04 -2.281100e-04 Aug 5 1.434523 -2.225275e-04 -2.225275e-04 Sep 5 1.377262 -2.617021e-04 -2.617021e-04 Oct 5 1.370866 -2.659122e-04 -2.659122e-04 Nov 5 1.355876 -2.760110e-04 -2.760110e-04 Dec 5 1.318202 -3.016438e-04 -3.016438e-04 Jan 6 1.289533 -8.791905e-05 9.671094e-04 Feb 6 1.322242 1.575428e-04 1.575429e-04 Mar 6 1.319944 1.561393e-04 1.561393e-04 Apr 6 1.316046 1.538242e-04 1.538242e-04 May 6 1.278768 1.324586e-04 1.324586e-04 Jun 6 1.252483 1.173888e-04 1.173888e-04 Jul 6 1.228696 1.037607e-04 1.037607e-04 Aug 6 1.239890 1.100797e-04 1.100797e-04 Sep 6 1.285464 1.359704e-04 1.359704e-04 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.113824588 0.436660604 0.716986476 -0.028267676 -0.259663600 2 0.556237940 0.002743341 2.069158061 0.586973589 -0.527576805 -0.020491904 3 -0.510875294 -1.185714998 0.717610880 0.383681425 1.235933876 0.984004262 4 -0.912175997 -1.506143990 -0.298641857 -0.420977165 -2.227551198 -0.933045010 5 0.207218023 0.919854110 0.940100266 1.187329246 -0.242110989 0.109880683 6 -0.787058251 0.837390396 -0.065505734 -0.108116300 -0.998330974 -0.704556798 Jul Aug Sep Oct Nov Dec 1 0.776495654 -0.267012492 0.713300920 0.862212741 1.194176702 -0.329928291 2 0.568257135 -2.128592054 -1.618632887 -2.795823256 -1.570763567 1.913731176 3 0.199431192 0.486985767 0.790193295 0.682683607 0.266194220 -0.794669310 4 1.511663744 0.345704528 0.475995809 2.229225129 -0.619034672 -1.162183068 5 -0.324645658 0.216642005 -1.521310737 -0.163606765 -0.392714244 -0.997478585 6 -0.637517753 0.295770193 1.212533751 > 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/1853l1353792696.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/20g8c1353792696.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/3aqnv1353792696.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/4u4bs1353792696.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/5lpr41353792696.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/6zkh61353792696.tab") > > try(system("convert tmp/1853l1353792696.ps tmp/1853l1353792696.png",intern=TRUE)) character(0) > try(system("convert tmp/20g8c1353792696.ps tmp/20g8c1353792696.png",intern=TRUE)) character(0) > try(system("convert tmp/3aqnv1353792696.ps tmp/3aqnv1353792696.png",intern=TRUE)) character(0) > try(system("convert tmp/4u4bs1353792696.ps tmp/4u4bs1353792696.png",intern=TRUE)) character(0) > try(system("convert tmp/5lpr41353792696.ps tmp/5lpr41353792696.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.989 0.449 3.446