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(9.3,9.3,8.7,8.2,8.3,8.5,8.6,8.5,8.2,8.1,7.9,8.6,8.7,8.7,8.5,8.4,8.5,8.7,8.7,8.6,8.5,8.3,8,8.2,8.1,8.1,8,7.9,7.9,8,8,7.9,8,7.7,7.2,7.5,7.3,7,7,7,7.2,7.3,7.1,6.8,6.4,6.1,6.5,7.7,7.9,7.5,6.9,6.6,6.9,7.7,8,8,7.7,7.3,7.4,8.1,8.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 0.070212130 0.000000000 0.003468003 0.000000000 > m$fitted level slope sea Jan 1 9.300000 0.000000000 0.000000000 Feb 1 9.300000 0.000000000 0.000000000 Mar 1 8.775252 -0.007674555 -0.075252306 Apr 1 8.252280 -0.009937357 -0.052279812 May 1 8.282142 -0.009761783 0.017858368 Jun 1 8.483991 -0.008749233 0.016008704 Jul 1 8.601599 -0.008140276 -0.001598616 Aug 1 8.522866 -0.008478384 -0.022866336 Sep 1 8.240938 -0.009780985 -0.040937930 Oct 1 8.113126 -0.010340543 -0.013125993 Nov 1 7.927894 -0.011165767 -0.027894235 Dec 1 8.524503 -0.008311441 0.075497232 Jan 2 8.612335 -0.012167485 0.087664673 Feb 2 8.646979 -0.011236464 0.053021094 Mar 2 8.541500 -0.012389540 -0.041499832 Apr 2 8.461709 -0.012654554 -0.061709410 May 2 8.493203 -0.012550050 0.006797437 Jun 2 8.670755 -0.012054283 0.029244549 Jul 2 8.707139 -0.011918915 -0.007139140 Aug 2 8.618649 -0.012134912 -0.018649357 Sep 2 8.538344 -0.012327570 -0.038344118 Oct 2 8.317420 -0.012943740 -0.017419552 Nov 2 8.118659 -0.013506129 -0.118658810 Dec 2 8.123058 -0.013510309 0.076942092 Jan 3 8.035906 -0.012396943 0.064094382 Feb 3 8.024531 -0.012386633 0.075468734 Mar 3 8.023357 -0.012272315 -0.023356565 Apr 3 7.973567 -0.012442197 -0.073566686 May 3 7.929405 -0.012512526 -0.029405470 Jun 3 7.968498 -0.012408486 0.031501525 Jul 3 7.994314 -0.012324143 0.005685947 Aug 3 7.928570 -0.012448048 -0.028569700 Sep 3 7.989191 -0.012270963 0.010808970 Oct 3 7.742827 -0.012873988 -0.042826602 Nov 3 7.385458 -0.013579951 -0.185458406 Dec 3 7.385070 -0.013601290 0.114930048 Jan 4 7.257006 -0.012745341 0.042994420 Feb 4 6.973931 -0.014365745 0.026069444 Mar 4 6.993921 -0.014079171 0.006078615 Apr 4 7.050091 -0.013741697 -0.050091401 May 4 7.204169 -0.013334641 -0.004168977 Jun 4 7.268780 -0.013191572 0.031219841 Jul 4 7.127779 -0.013436622 -0.027778650 Aug 4 6.896137 -0.013893388 -0.096136547 Sep 4 6.441597 -0.014891700 -0.041597263 Oct 4 6.130192 -0.015565855 -0.030191696 Nov 4 6.567842 -0.014979443 -0.067842043 Dec 4 7.383888 -0.016470880 0.316112377 Jan 5 7.739018 -0.017967873 0.160982472 Feb 5 7.539857 -0.018706126 -0.039857058 Mar 5 7.025310 -0.022100513 -0.125309620 Apr 5 6.739333 -0.023360445 -0.139332885 May 5 6.881261 -0.022922251 0.018739063 Jun 5 7.510254 -0.021727306 0.189746353 Jul 5 7.920125 -0.020958759 0.079874698 Aug 5 8.008045 -0.020745546 -0.008044811 Sep 5 7.749225 -0.021250210 -0.049225432 Oct 5 7.474393 -0.021745448 -0.174392774 Nov 5 7.589220 -0.021636756 -0.189219524 Dec 5 7.799473 -0.021998279 0.300526654 Jan 6 8.033824 -0.022566307 0.266176146 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.000000000 -1.970884653 -1.949535859 0.150511606 0.800213213 2 0.426479469 0.156250978 -0.350195989 -0.254804525 0.166823769 0.718216131 3 -0.295275808 0.003660628 0.041523274 -0.141689227 -0.119854609 0.194940091 4 -0.444270438 -0.992453666 0.127358033 0.264930561 0.634115095 0.294440075 5 1.422363761 -0.672434601 -1.842282858 -0.993893188 0.624450051 2.462670217 6 0.974579543 Jul Aug Sep Oct Nov Dec 1 0.477817965 -0.266942510 -1.034037355 -0.446323798 -0.661329766 2.298203042 2 0.182991857 -0.289273523 -0.257541197 -0.788169999 -0.702266478 0.067586906 3 0.144373967 -0.201766120 0.276004804 -0.884472879 -1.301295624 0.049845095 4 -0.482720137 -0.824109094 -1.664462543 -1.120246689 1.711188254 3.142493053 5 1.630123393 0.411219468 -0.899308835 -0.957896052 0.515519068 0.877003018 6 > 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/1nd6x1259941088.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/2oc3c1259941088.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/3s7v91259941088.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/4jz9p1259941088.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/57r311259941088.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/67cqu1259941088.tab") > system("convert tmp/1nd6x1259941088.ps tmp/1nd6x1259941088.png") > system("convert tmp/2oc3c1259941088.ps tmp/2oc3c1259941088.png") > system("convert tmp/3s7v91259941088.ps tmp/3s7v91259941088.png") > system("convert tmp/4jz9p1259941088.ps tmp/4jz9p1259941088.png") > system("convert tmp/57r311259941088.ps tmp/57r311259941088.png") > > > proc.time() user system elapsed 1.380 0.819 2.849