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(14.3,14.2,15.9,15.3,15.5,15.1,15,12.1,15.8,16.9,15.1,13.7,14.8,14.7,16,15.4,15,15.5,15.1,11.7,16.3,16.7,15,14.9,14.6,15.3,17.9,16.4,15.4,17.9,15.9,13.9,17.8,17.9,17.4,16.7,16,16.6,19.1,17.8,17.2,18.6,16.3,15.1,19.2,17.7,19.1,18,17.5,17.8,21.1,17.2,19.4,19.8,17.6,16.2,19.5,19.9,20,17.3) > 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 2.236437e-02 8.498127e-06 2.802491e-01 9.289072e-02 > m$fitted level slope sea Jan 1 14.30000 0.000000000 0.000000000 Feb 1 14.29501 0.001666172 -0.078026503 Mar 1 14.71161 0.049058153 0.991909141 Apr 1 15.07550 0.086860039 0.118197521 May 1 15.31003 0.102913073 0.133530815 Jun 1 15.35633 0.097576087 -0.230424390 Jul 1 15.31967 0.086751539 -0.249454492 Aug 1 14.67506 0.036003433 -2.157691399 Sep 1 14.67034 0.033529125 1.154184900 Oct 1 15.10719 0.055483068 1.541570030 Nov 1 15.28485 0.061573077 -0.262571162 Dec 1 15.06704 0.048584192 -1.186935301 Jan 2 15.04189 0.048302266 -0.188765020 Feb 2 15.10303 0.048460374 -0.411868938 Mar 2 15.17668 0.049208336 0.807872597 Apr 2 15.24825 0.050122613 0.139121678 May 2 15.17833 0.044635738 -0.112223283 Jun 2 15.20702 0.043896803 0.301871068 Jul 2 15.16062 0.039870969 -0.008923493 Aug 2 14.95225 0.029524020 -3.106067464 Sep 2 14.99465 0.030017302 1.297587230 Oct 2 15.04772 0.030818766 1.638104656 Nov 2 15.06905 0.030524925 -0.063119639 Dec 2 15.25849 0.034749726 -0.459089932 Jan 3 15.25908 0.034028327 -0.637080809 Feb 3 15.35059 0.035232163 -0.087393659 Mar 3 15.69527 0.042658330 2.011183067 Apr 3 15.90909 0.047366036 0.386473660 May 3 15.91333 0.046071217 -0.487399426 Jun 3 16.20762 0.053819217 1.543489629 Jul 3 16.24327 0.053251914 -0.332304975 Aug 3 16.42490 0.057147050 -2.603119324 Sep 3 16.53959 0.058808944 1.224986051 Oct 3 16.57567 0.058193493 1.338451753 Nov 3 16.80878 0.062591828 0.481606915 Dec 3 17.00310 0.065658231 -0.386172679 Jan 4 17.09021 0.066126867 -1.103799138 Feb 4 17.15638 0.066127803 -0.556409569 Mar 4 17.23981 0.066511607 1.849315553 Apr 4 17.35248 0.067585160 0.418738356 May 4 17.53169 0.070292218 -0.400892840 Jun 4 17.55814 0.069201650 1.068987210 Jul 4 17.46950 0.065252021 -1.071747457 Aug 4 17.53285 0.065205139 -2.431677219 Sep 4 17.68382 0.067274561 1.462640822 Oct 4 17.57878 0.063258011 0.229264658 Nov 4 17.78008 0.066353145 1.233022379 Dec 4 17.99287 0.069518180 -0.085310080 Jan 5 18.21114 0.072648558 -0.805199507 Feb 5 18.35597 0.074152314 -0.601592523 Mar 5 18.60926 0.077907725 2.377705517 Apr 5 18.42137 0.072248660 -1.054089611 May 5 18.61436 0.074862619 0.709814612 Jun 5 18.71544 0.075437007 1.068107392 Jul 5 18.78614 0.075332676 -1.183160114 Aug 5 18.83322 0.074714335 -2.615471453 Sep 5 18.74607 0.071213701 0.855825511 Oct 5 18.95909 0.074228417 0.851469733 Nov 5 19.07732 0.075146222 0.894873405 Dec 5 18.89415 0.069853314 -1.430694375 > m$resid Jan Feb Mar Apr May 1 0.0000000000 -0.1352068145 1.8109067720 1.1133151501 0.5893142490 2 -0.5060420891 0.0844505345 0.1500911078 0.1251015674 -0.6620659174 3 -0.2156398170 0.3608379016 1.9038111975 1.0325654051 -0.2575175085 4 0.1341257949 0.0002713284 0.1073958712 0.2847176132 0.6856514081 5 0.9299803443 0.4510509718 1.1176810535 -1.6549018990 0.7505469951 Jun Jul Aug Sep Oct 1 -0.2616060188 -0.6898258340 -4.0246113708 -0.2337448020 2.3763742145 2 -0.0892618040 -0.5180058376 -1.4588869482 0.0771677660 0.1402913295 3 1.4821597514 -0.1091966238 0.7782667819 0.3519121174 -0.1400422973 4 -0.2690516218 -0.9700423046 -0.0116829236 0.5308215758 -1.0704131594 5 0.1628880333 -0.0294863541 -0.1757963584 -1.0089586160 0.8853857816 Nov Dec 1 0.7316248447 -1.6901631187 2 -0.0584266714 0.9885354997 3 1.0847823523 0.8209275962 4 0.8601460878 0.9144362456 5 0.2751551176 -1.6169162778 > 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/1pldg1259957917.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/2dw9l1259957917.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/358ih1259957917.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/4aen01259957917.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/544k21259957917.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/6cnnj1259957917.tab") > system("convert tmp/1pldg1259957917.ps tmp/1pldg1259957917.png") > system("convert tmp/2dw9l1259957917.ps tmp/2dw9l1259957917.png") > system("convert tmp/358ih1259957917.ps tmp/358ih1259957917.png") > system("convert tmp/4aen01259957917.ps tmp/4aen01259957917.png") > system("convert tmp/544k21259957917.ps tmp/544k21259957917.png") > > > proc.time() user system elapsed 1.758 0.835 3.582