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(10284.5,12792,12823.61538,13845.66667,15335.63636,11188.5,13633.25,12298.46667,15353.63636,12696.15385,12213.93333,13683.72727,11214.14286,13950.23077,11179.13333,11801.875,11188.82353,16456.27273,11110.0625,16530.69231,10038.41176,11681.25,11148.88235,8631,9386.444444,9764.736842,12043.75,12948.06667,10987.125,11648.3125,10633.35294,10219.3,9037.6,10296.31579,11705.41176,10681.94444,9362.947368,11306.35294,10984.45,10062.61905,8118.583333,8867.48,8346.72,8529.307692,10697.18182,8591.84,8695.607143,8125.571429,7009.758621,7883.466667,7527.645161,6763.758621,6682.333333,7855.681818,6738.88,7895.434783,6361.884615,6935.956522,8344.454545,9107.944444) > 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 267319.7715 987.5048 0.0000 1754672.6020 > m$fitted level slope sea Jan 1 10284.500 0.000000 0.0000 Feb 1 11178.001 111.819531 110.9737 Mar 1 11805.021 160.528437 158.4458 Apr 1 12584.113 208.035118 203.3247 May 1 13622.770 262.895644 252.8787 Jun 1 12849.417 201.275522 200.2919 Jul 1 13183.993 208.670161 206.2016 Aug 1 12941.190 184.690067 188.3518 Sep 1 13832.403 221.277947 213.6515 Oct 1 13509.667 193.479256 195.8063 Nov 1 13119.416 163.821104 178.1219 Dec 1 13360.361 167.735473 180.2924 Jan 2 13455.789 171.974990 -2075.0780 Feb 2 13681.019 175.302355 190.8096 Mar 2 12786.365 114.535444 156.8254 Apr 2 12446.837 90.348588 145.3504 May 2 12005.267 63.108201 134.3896 Jun 2 13572.308 138.462524 160.3138 Jul 2 12737.476 90.194705 145.9384 Aug 2 14080.490 152.136107 162.1298 Sep 2 12697.955 76.162884 144.4687 Oct 2 12337.881 54.487599 139.9355 Nov 2 11904.118 30.109050 135.3040 Dec 2 10719.276 -30.879645 124.6965 Jan 3 10708.156 -30.455625 -1358.8309 Feb 3 10285.089 -52.803222 114.4034 Mar 3 10853.570 -18.859856 125.1041 Apr 3 11551.294 19.309052 134.3777 May 3 11313.689 5.826106 131.8273 Jun 3 11389.802 9.485150 132.3758 Jul 3 11079.063 -7.121394 130.3603 Aug 3 10721.574 -25.267841 128.5352 Sep 3 10059.118 -58.276751 125.7212 Oct 3 10061.388 -55.136878 125.9528 Nov 3 10567.643 -25.980254 127.8475 Dec 3 10546.102 -25.749303 127.8609 Jan 4 10596.924 -22.540364 -1369.6492 Feb 4 10800.099 -10.163871 128.1217 Mar 4 10814.066 -8.864101 128.3950 Apr 4 10491.656 -25.574582 125.8657 May 4 9579.165 -72.552063 120.7135 Jun 4 9234.528 -86.913839 119.5577 Jul 4 8818.134 -104.273497 118.5143 Aug 4 8605.368 -109.985043 118.2524 Sep 4 9241.426 -70.719522 119.6590 Oct 4 8920.525 -83.888329 119.2809 Nov 4 8743.381 -88.798390 119.1651 Dec 4 8422.326 -101.031880 118.9219 Jan 5 8320.210 -101.084565 -1308.5637 Feb 5 8051.801 -110.088258 116.9013 Mar 5 7749.456 -120.351615 115.2729 Apr 5 7276.418 -139.094534 113.2227 May 5 6933.481 -149.899682 112.4095 Jun 5 7127.599 -131.690821 113.3529 Jul 5 6863.149 -138.712084 113.1019 Aug 5 7103.692 -118.663031 113.5983 Sep 5 6720.850 -132.625841 113.3576 Oct 5 6672.282 -128.183348 113.4113 Nov 5 7149.222 -96.202269 113.6852 Dec 5 7749.452 -59.393360 113.9112 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 1.465566840 0.835148916 1.018982895 1.395403175 -1.768942591 2 -0.177190418 0.077326903 -1.698038918 -0.752090679 -0.899678740 2.569515755 3 0.038401264 -0.612334328 1.014941198 1.196011614 -0.433050642 0.118992462 4 0.138772585 0.361810251 0.039859506 -0.524824269 -1.492974621 -0.459042860 5 -0.001911175 -0.271848965 -0.319378039 -0.591423223 -0.343118223 0.579842259 Jul Aug Sep Oct Nov Dec 1 0.230111018 -0.784802391 1.232786997 -0.950668440 -1.020092915 0.134653735 2 -1.670453347 2.153693487 -2.638474587 -0.749474527 -0.837926576 -2.082459019 3 -0.543066277 -0.594467743 -1.080987510 0.102677158 0.951565876 0.007520941 4 -0.556291184 -0.183198761 1.259570498 -0.422294456 -0.157371667 -0.391855630 5 -0.223841078 0.639417549 -0.445327217 0.141671873 1.019706089 1.173418999 > 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/1ah0t1259888307.ps",horizontal=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/2qcf71259888307.ps",horizontal=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/32tll1259888307.ps",horizontal=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/48hv21259888307.ps",horizontal=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/5vrnl1259888307.ps",horizontal=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/6hroz1259888307.tab") > system("convert tmp/1ah0t1259888307.ps tmp/1ah0t1259888307.png") > system("convert tmp/2qcf71259888307.ps tmp/2qcf71259888307.png") > system("convert tmp/32tll1259888307.ps tmp/32tll1259888307.png") > system("convert tmp/48hv21259888307.ps tmp/48hv21259888307.png") > system("convert tmp/5vrnl1259888307.ps tmp/5vrnl1259888307.png") > > > proc.time() user system elapsed 1.445 0.796 1.702