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(95.5,76.7,79.4,55.2,60,64.8,82.3,210.5,106,80.8,97.3,189.5,90,69.3,87.3,57.4,56.2,61.6,77.7,177.2,97.6,81.6,96.8,191.3,106,75.1,72,63.5,57.4,62.3,79.4,178.1,109.3,85.2,102.7,193.7,108.4,73.4,85.9,58.5,58.6,62.7,77.5,180.5,102.2,82.6,97.8,197.8,93.8,72.4,77.7,58.7,53.1,64.3,76.4,188.4,105.5,79.8,96.1,202.5) > 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 94.37676 0.00000 333.70614 0.00000 > m$fitted level slope sea Jan 1 95.50000 0.000000000 0.000000 Feb 1 89.96471 0.140802314 -13.264705 Mar 1 81.40203 -0.459989954 -2.002031 Apr 1 69.27608 -1.341430298 -14.076082 May 1 62.99779 -1.653136719 -2.997786 Jun 1 61.56256 -1.642996413 3.237437 Jul 1 67.51762 -1.387489130 14.782376 Aug 1 114.29496 -0.136275047 96.205042 Sep 1 128.00610 0.173220851 -22.006099 Oct 1 117.98865 -0.041389198 -37.188646 Nov 1 107.98023 -0.248780871 -10.680233 Dec 1 130.11518 0.217622488 59.384823 Jan 2 123.30042 0.282788229 -33.300416 Feb 2 105.80923 0.395814878 -36.509229 Mar 2 91.26584 0.214734450 -3.965841 Apr 2 79.00719 -0.114384805 -21.607191 May 2 71.22236 -0.361627021 -15.022356 Jun 2 71.69323 -0.335202581 -10.093231 Jul 2 79.88027 -0.096702134 -2.180270 Aug 2 84.41333 0.011732415 92.786673 Sep 2 94.54906 0.206414978 3.050938 Oct 2 105.34426 0.370649313 -23.744256 Nov 2 112.60729 0.452520566 -15.807289 Dec 2 116.74461 0.482648565 74.555387 Jan 3 118.99308 0.491994336 -12.993083 Feb 3 112.69962 0.450329664 -37.599624 Mar 3 96.13738 0.271229062 -24.137376 Apr 3 86.59990 0.115866702 -23.099896 May 3 80.92824 0.001628779 -23.528242 Jun 3 79.70867 -0.024391302 -17.408675 Jul 3 82.86797 0.041600859 -3.467965 Aug 3 88.51295 0.146339633 89.587053 Sep 3 98.99958 0.311092371 10.300422 Oct 3 107.11332 0.412662587 -21.913323 Nov 3 113.38246 0.473388736 -10.682464 Dec 3 115.56431 0.487741403 78.135694 Jan 4 114.02198 0.472383849 -5.621985 Feb 4 107.71363 0.416998348 -34.313634 Mar 4 103.98498 0.375531961 -18.084978 Apr 4 95.17027 0.261982488 -36.670272 May 4 89.41898 0.175409249 -30.818976 Jun 4 87.20971 0.138330828 -24.509714 Jul 4 87.85508 0.146247785 -10.355078 Aug 4 92.61493 0.214175157 87.885071 Sep 4 96.28811 0.259832157 5.911886 Oct 4 101.08823 0.311755017 -18.488231 Nov 4 104.38791 0.341079032 -6.587906 Dec 4 108.31794 0.372148233 89.482064 Jan 5 104.23811 0.335710036 -10.438108 Feb 5 101.48565 0.309594627 -29.085652 Mar 5 96.32765 0.258590000 -18.627648 Apr 5 93.28421 0.223960193 -34.584210 May 5 89.65271 0.179384305 -36.552711 Jun 5 89.49645 0.175273917 -25.196448 Jul 5 90.52332 0.185817810 -14.123320 Aug 5 95.47926 0.242924314 92.920743 Sep 5 99.65493 0.286797050 5.845069 Oct 5 101.18896 0.299457044 -21.388964 Nov 5 102.01258 0.304275260 -5.912585 Dec 5 103.92108 0.317886205 98.578920 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.69678003 -0.57045023 -0.87300355 -0.44332713 0.02133846 2 -0.74533209 -1.89677975 -1.47875155 -1.16712315 -0.71992142 0.08052992 3 0.18107230 -0.69186622 -1.70341721 -0.96385048 -0.56521188 -0.12010057 4 -0.20670884 -0.68729250 -0.41694189 -0.91722377 -0.59799876 -0.23756157 5 -0.45225487 -0.31287126 -0.55181544 -0.33204706 -0.38695443 -0.03370214 Jul Aug Sep Oct Nov Dec 1 0.76688109 4.90587188 1.41329664 -1.03956184 -1.01551536 2.27822802 2 0.84703975 0.46763406 1.03042337 1.08077871 0.70407199 0.37691978 3 0.31700066 0.56403407 1.04807739 0.79367078 0.59669015 0.17416761 4 0.05077527 0.46473008 0.35004452 0.46082432 0.30373648 0.36496615 5 0.08570531 0.48152764 0.39810656 0.12651439 0.05323251 0.16297317 > 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/14tgk1259948670.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/2kimm1259948670.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/3o4eg1259948670.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/4a34y1259948670.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/5vgvn1259948670.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/6g2jz1259948670.tab") > system("convert tmp/14tgk1259948670.ps tmp/14tgk1259948670.png") > system("convert tmp/2kimm1259948670.ps tmp/2kimm1259948670.png") > system("convert tmp/3o4eg1259948670.ps tmp/3o4eg1259948670.png") > system("convert tmp/4a34y1259948670.ps tmp/4a34y1259948670.png") > system("convert tmp/5vgvn1259948670.ps tmp/5vgvn1259948670.png") > > > proc.time() user system elapsed 1.349 0.824 1.576