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(12.6,15.7,13.2,20.3,12.8,8,0.9,3.6,14.1,21.7,24.5,18.9,13.9,11,5.8,15.5,22.4,31.7,30.3,31.4,20.2,19.7,10.8,13.2,15.1,15.6,15.5,12.7,10.9,10,9.1,10.3,16.9,22,27.6,28.9,31,32.9,38.1,28.8,29,21.8,28.8,25.6,28.2,20.2,17.9,16.3,13.2,8.1,4.5,-0.1,0,2.3,2.8,2.9,0.1,3.5,8.6,13.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 25.29345 0.00000 0.00000 0.00000 > m$fitted level slope sea Jan 1 12.600000000 0.000000000 0.000000000 Feb 1 15.538800017 0.161199983 0.161199983 Mar 1 13.049402406 0.150597594 0.150597594 Apr 1 20.121825416 0.178174584 0.178174584 May 1 12.652173928 0.147826072 0.147826072 Jun 1 7.871653556 0.128346444 0.128346444 Jul 1 0.800000010 0.099999990 0.099999990 Aug 1 3.489843761 0.110156239 0.110156239 Sep 1 13.949416358 0.150583642 0.150583642 Oct 1 21.520542654 0.179457346 0.179457346 Nov 1 24.310424730 0.189575270 0.189575270 Dec 1 18.732692325 0.167307675 0.167307675 Jan 2 16.743831821 0.258530166 -2.843831821 Feb 2 10.889454551 0.110545449 0.110545449 Mar 2 5.699092563 0.100907437 0.100907437 Apr 2 15.381702904 0.118297096 0.118297096 May 2 22.269439427 0.130560573 0.130560573 Jun 2 31.552888094 0.147111906 0.147111906 Jul 2 30.155675682 0.144324318 0.144324318 Aug 2 31.253956841 0.146043159 0.146043159 Sep 2 20.074326756 0.125673244 0.125673244 Oct 2 19.575448034 0.124551966 0.124551966 Nov 2 10.691592134 0.108407866 0.108407866 Dec 2 13.087500005 0.112499995 0.112499995 Jan 3 15.780029344 0.061820850 -0.680029344 Feb 3 15.542823534 0.057176469 0.057176466 Mar 3 15.443008227 0.056991773 0.056991773 Apr 3 12.646361504 0.053638496 0.053638496 May 3 10.848534585 0.051465415 0.051465415 Jun 3 9.949648713 0.050351287 0.050351287 Jul 3 9.050760235 0.049239765 0.049239765 Aug 3 10.249415889 0.050584111 0.050584111 Sep 3 16.841773630 0.058226370 0.058226370 Oct 3 21.935897437 0.064102563 0.064102563 Nov 3 27.529452854 0.070547146 0.070547146 Dec 3 28.828023258 0.071976742 0.071976742 Jan 4 31.435066248 0.039551477 -0.435066248 Feb 4 32.844782613 0.055217390 0.055217387 Mar 4 38.040312773 0.059687227 0.059687227 Apr 4 28.748437501 0.051562499 0.051562499 May 4 28.948308761 0.051691239 0.051691239 Jun 4 21.754592722 0.045407278 0.045407278 Jul 4 28.748571430 0.051428570 0.051428570 Aug 4 25.551384084 0.048615916 0.048615916 Sep 4 28.149178912 0.050821088 0.050821088 Oct 4 20.156131262 0.043868738 0.043868738 Nov 4 17.858153581 0.041846419 0.041846419 Dec 4 16.259568966 0.040431034 0.040431034 Jan 5 13.895550337 0.063231849 -0.695550337 Feb 5 8.089862070 0.010137931 0.010137930 Mar 5 4.492350103 0.007649897 0.007649897 Apr 5 -0.104476584 0.004476584 0.004476584 May 5 -0.004542327 0.004542327 0.004542327 Jun 5 2.293878954 0.006121046 0.006121046 Jul 5 2.793539519 0.006460481 0.006460481 Aug 5 2.893475274 0.006524726 0.006524726 Sep 5 0.095401510 0.004598490 0.004598490 Oct 5 3.493072702 0.006927298 0.006927298 Nov 5 8.589581905 0.010418095 0.010418095 Dec 5 13.786027397 0.013972603 0.013972603 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.35085703 -0.52808836 1.37904993 -1.52368058 -0.98186949 2 -0.51898449 -1.03612776 -1.05497123 1.90691959 1.34722977 1.82157228 3 0.56571817 -0.05404598 -0.03123404 -0.56774055 -0.36835474 -0.18907518 4 0.53835720 0.25450231 1.02252569 -1.86023880 0.02950198 -1.44127546 5 -0.50240875 -1.10676632 -0.71757951 -0.91585308 0.01898700 0.45626361 Jul Aug Sep Oct Nov Dec 1 -1.43443762 0.51596399 2.06185622 1.47834193 0.52005245 -1.14896249 2 -0.30734493 0.18985217 -2.25398057 -0.12429510 -1.79280389 0.45524497 3 -0.18885391 0.22867936 1.30150242 1.00190386 1.10009714 0.24431784 4 1.38222758 -0.64622273 0.50708880 -1.60010512 -0.46584537 -0.32631812 5 0.09816737 0.01859268 -0.55784780 0.67489792 1.01234124 1.03152445 > 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/1u9qy1259943964.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/2o6051259943964.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/3no8u1259943964.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/4xdcw1259943964.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/5n9tt1259943964.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/6ofyg1259943964.tab") > system("convert tmp/1u9qy1259943964.ps tmp/1u9qy1259943964.png") > system("convert tmp/2o6051259943964.ps tmp/2o6051259943964.png") > system("convert tmp/3no8u1259943964.ps tmp/3no8u1259943964.png") > system("convert tmp/4xdcw1259943964.ps tmp/4xdcw1259943964.png") > system("convert tmp/5n9tt1259943964.ps tmp/5n9tt1259943964.png") > > > proc.time() user system elapsed 1.294 0.782 1.531