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(8,8.1,7.7,7.5,7.6,7.8,7.8,7.8,7.5,7.5,7.1,7.5,7.5,7.6,7.7,7.7,7.9,8.1,8.2,8.2,8.2,7.9,7.3,6.9,6.6,6.7,6.9,7,7.1,7.2,7.1,6.9,7,6.8,6.4,6.7,6.6,6.4,6.3,6.2,6.5,6.8,6.8,6.4,6.1,5.8,6.1,7.2,7.3,6.9,6.1,5.8,6.2,7.1,7.7,7.9,7.7,7.4,7.5,8) > 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 0.0000000 0.1044248 0.0004892 0.0000000 > m$fitted level slope sea Jan 1 8.000000 0.000000000 0.000000e+00 Feb 1 8.099304 0.099361933 6.955508e-04 Mar 1 7.709670 -0.379823097 -9.670232e-03 Apr 1 7.495451 -0.218449602 4.549188e-03 May 1 7.595752 0.092296634 4.247627e-03 Jun 1 7.799987 0.201429194 1.348273e-05 Jul 1 7.804406 0.009358138 -4.405785e-03 Aug 1 7.799897 -0.004161079 1.029332e-04 Sep 1 7.505462 -0.287157700 -5.462421e-03 Oct 1 7.494490 -0.017896203 5.510020e-03 Nov 1 7.108080 -0.377171140 -8.080477e-03 Dec 1 7.485896 0.358888262 1.410430e-02 Jan 2 7.508172 0.030766745 -8.171527e-03 Feb 2 7.585776 0.076443187 1.422416e-02 Mar 2 7.708085 0.120444403 -8.084903e-03 Apr 2 7.705034 0.002242907 -5.033961e-03 May 2 7.894944 0.181877838 5.056015e-03 Jun 2 8.095222 0.199492864 4.777677e-03 Jul 2 8.211569 0.119896450 -1.156864e-02 Aug 2 8.191363 -0.014222644 8.636829e-03 Sep 2 8.217621 0.024529150 -1.762094e-02 Oct 2 7.883348 -0.318951137 1.665188e-02 Nov 2 7.336641 -0.536982620 -3.664055e-02 Dec 2 6.872748 -0.467013724 2.725169e-02 Jan 3 6.603229 -0.277984677 -3.228808e-03 Feb 3 6.680778 0.062488485 1.922172e-02 Mar 3 6.895539 0.206713095 4.461066e-03 Apr 3 7.016708 0.125703229 -1.670799e-02 May 3 7.098777 0.084397081 1.222903e-03 Jun 3 7.193595 0.094263466 6.405136e-03 Jul 3 7.109702 -0.074417065 -9.701635e-03 Aug 3 6.905005 -0.197765369 -5.005252e-03 Sep 3 6.996835 0.076422524 3.165106e-03 Oct 3 6.784868 -0.196624016 1.513152e-02 Nov 3 6.439987 -0.336995347 -3.998676e-02 Dec 3 6.651503 0.182325881 4.849669e-02 Jan 4 6.625609 -0.014782824 -2.560915e-02 Feb 4 6.399111 -0.215285325 8.890629e-04 Mar 4 6.285692 -0.119487006 1.430812e-02 Apr 4 6.211853 -0.076526870 -1.185284e-02 May 4 6.491816 0.258765321 8.183830e-03 Jun 4 6.787135 0.293151257 1.286550e-02 Jul 4 6.806088 0.035205207 -6.088276e-03 Aug 4 6.444438 -0.338123842 -4.443838e-02 Sep 4 6.083219 -0.359850099 1.678095e-02 Oct 4 5.758319 -0.326972000 4.168107e-02 Nov 4 6.147582 0.346809114 -4.758226e-02 Dec 4 7.125677 0.940648154 7.432300e-02 Jan 5 7.345755 0.262895393 -4.575539e-02 Feb 5 6.928880 -0.376615232 -2.888049e-02 Mar 5 6.090372 -0.809210616 9.628403e-03 Apr 5 5.806305 -0.316903598 -6.305136e-03 May 5 6.175410 0.325773629 2.458969e-02 Jun 5 7.064527 0.853623163 3.547304e-02 Jul 5 7.699526 0.648766878 4.736443e-04 Aug 5 7.959281 0.284256019 -5.928111e-02 Sep 5 7.680737 -0.243093421 1.926295e-02 Oct 5 7.380316 -0.296809910 1.968356e-02 Nov 5 7.584801 0.172913430 -8.480115e-02 Dec 5 7.901772 0.307889335 9.822788e-02 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.30851125 -1.49769066 0.49919367 0.96163576 0.33771673 2 -1.01562862 0.14172132 0.13715506 -0.36553550 0.55592200 0.05451059 3 0.58511811 1.05523539 0.44834626 -0.25053751 -0.12783229 0.03053203 4 -0.61016753 -0.62089781 0.29726893 0.13289122 1.03760437 0.10640968 5 -2.09813128 -1.97945587 -1.34095569 1.52319813 1.98873886 1.63348204 Jul Aug Sep Oct Nov Dec 1 -0.59437442 -0.04183596 -0.87574857 0.83324447 -1.11179600 2.27777620 2 -0.24631552 -0.41503889 0.11991955 -1.06291859 -0.67471044 0.21652276 3 -0.52199143 -0.38170809 0.84848949 -0.84495750 -0.43438707 1.60707073 4 -0.79822817 -1.15528735 -0.06723310 0.10174308 2.08505651 1.83767299 5 -0.63393818 -1.12799918 -1.63191183 -0.16622864 1.45358847 0.41769384 > 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/1pgt61260545539.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/2z4qc1260545539.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/3mzqn1260545539.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/4s15z1260545539.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/5eyd41260545539.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/6v3m11260545540.tab") > system("convert tmp/1pgt61260545539.ps tmp/1pgt61260545539.png") > system("convert tmp/2z4qc1260545539.ps tmp/2z4qc1260545539.png") > system("convert tmp/3mzqn1260545539.ps tmp/3mzqn1260545539.png") > system("convert tmp/4s15z1260545539.ps tmp/4s15z1260545539.png") > system("convert tmp/5eyd41260545539.ps tmp/5eyd41260545539.png") > > > proc.time() user system elapsed 1.349 0.799 2.333