R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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(286602,283042,276687,277915,277128,277103,275037,270150,267140,264993,287259,291186,292300,288186,281477,282656,280190,280408,276836,275216,274352,271311,289802,290726,292300,278506,269826,265861,269034,264176,255198,253353,246057,235372,258556,260993,254663,250643,243422,247105,248541,245039,237080,237085,225554,226839,247934,248333,246969,245098,246263,255765,264319,268347,273046,273963,267430,271993,292710,295881,294563) > par1 = '12' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 13624333.8 695753.7 4062962.7 0.0 > m$fitted level slope sea Jan 1 286602.0 0.00000 0.00000 Feb 1 283653.5 -194.91878 -611.51468 Mar 1 278738.1 -633.58180 -2051.11211 Apr 1 277701.9 -680.25227 213.09009 May 1 277271.5 -645.89077 -143.53024 Jun 1 277213.7 -554.77982 -110.66623 Jul 1 275810.0 -696.86167 -772.96735 Aug 1 271748.6 -1287.72186 -1598.63110 Sep 1 267979.8 -1736.67874 -839.81954 Oct 1 265323.0 -1906.33467 -330.00305 Nov 1 280417.1 1265.66684 6841.86909 Dec 1 290454.5 2914.69562 731.49441 Jan 2 291827.6 2635.81463 472.39733 Feb 2 289105.3 1630.43337 -919.30617 Mar 2 285177.4 601.23235 -3700.44790 Apr 2 283202.1 121.71928 -546.09991 May 2 281322.9 -253.22498 -1132.89691 Jun 2 280363.3 -385.76641 44.74147 Jul 2 277419.9 -865.53841 -583.91582 Aug 2 275491.2 -1065.00033 -275.15770 Sep 2 274612.0 -1030.13292 -259.95875 Oct 2 276026.2 -571.32676 -4715.19891 Nov 2 282777.3 802.02631 7024.74212 Dec 2 288224.9 1670.53225 2501.08415 Jan 3 289964.1 1683.34735 2335.88242 Feb 3 281998.1 -124.04551 -3492.08965 Mar 3 274827.8 -1432.84780 -5001.81802 Apr 3 267747.9 -2481.02532 -1886.94835 May 3 267801.9 -2008.27966 1232.11529 Jun 3 264104.8 -2323.88077 71.17128 Jul 3 257203.5 -3179.21790 -2005.49419 Aug 3 253381.6 -3299.27138 -28.57048 Sep 3 248384.5 -3616.38993 -2327.47665 Oct 3 243648.9 -3825.38359 -8276.86511 Nov 3 248979.8 -2117.60317 9576.23309 Dec 3 255068.4 -589.40855 5924.57689 Jan 4 251037.5 -1231.03727 3625.50182 Feb 4 250527.0 -1096.54168 115.96335 Mar 4 247730.7 -1412.43091 -4308.69934 Apr 4 248642.9 -981.20942 -1537.86615 May 4 247027.3 -1099.16017 1513.66853 Jun 4 244050.1 -1449.09027 988.92431 Jul 4 239692.6 -1991.14702 -2612.62593 Aug 4 236430.5 -2227.95475 654.50653 Sep 4 229940.0 -3021.87532 -4386.02908 Oct 4 234615.3 -1589.01480 -7776.25174 Nov 4 238931.8 -490.69056 9002.15808 Dec 4 240720.3 -66.99031 7612.66580 Jan 5 242754.7 324.15135 4214.26419 Feb 5 244328.8 556.90417 769.23370 Mar 5 249469.1 1408.28018 -3206.09365 Apr 5 255350.1 2237.72037 414.90006 May 5 260783.6 2831.00444 3535.44242 Jun 5 265445.6 3171.47295 2901.41765 Jul 5 272785.4 3946.95829 260.61853 Aug 5 273978.8 3434.81371 -15.81915 Sep 5 275520.9 3082.95113 -8090.91993 Oct 5 280665.5 3465.97178 -8672.54459 Nov 5 284488.0 3532.16319 8221.96763 Dec 5 288639.7 3647.18967 7241.31795 Jan 6 291496.4 3500.26913 3066.56364 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.51350417 -1.10580595 -0.09970652 0.06023274 0.13914372 2 -0.38176891 -1.17241416 -1.19525698 -0.58052520 -0.45168185 -0.15866973 3 0.01582875 -2.15946616 -1.54256079 -1.25987165 0.56969686 -0.37826232 4 -0.77886818 0.16108781 -0.37513901 0.51687803 -0.14203243 -0.42016687 5 0.47200839 0.27881746 1.01476387 0.99337097 0.71358123 0.40910235 6 -0.17683248 Jul Aug Sep Oct Nov Dec 1 -0.19855381 -0.78068784 -0.57258213 -0.21157137 3.89954991 2.00887184 2 -0.57424836 -0.23896916 0.04179673 0.55005764 1.64583404 1.04459986 3 -1.02322730 -0.14370409 -0.37986353 -0.25037365 2.04668896 1.84078304 4 -0.64916543 -0.28346434 -0.95057001 1.71599627 1.31711125 0.51018690 5 0.92957939 -0.61316040 -0.42117504 0.45868569 0.07940800 0.13839909 6 > 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/html/rcomp/tmp/1efgx1292521936.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/html/rcomp/tmp/2efgx1292521936.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/html/rcomp/tmp/36og01292521936.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/html/rcomp/tmp/46og01292521936.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/html/rcomp/tmp/5hyf31292521936.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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/html/rcomp/tmp/6d7db1292521936.tab") > > try(system("convert tmp/1efgx1292521936.ps tmp/1efgx1292521936.png",intern=TRUE)) character(0) > try(system("convert tmp/2efgx1292521936.ps tmp/2efgx1292521936.png",intern=TRUE)) character(0) > try(system("convert tmp/36og01292521936.ps tmp/36og01292521936.png",intern=TRUE)) character(0) > try(system("convert tmp/46og01292521936.ps tmp/46og01292521936.png",intern=TRUE)) character(0) > try(system("convert tmp/5hyf31292521936.ps tmp/5hyf31292521936.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.319 0.897 4.889