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(89.1,82.6,102.7,91.8,94.1,103.1,93.2,91,94.3,99.4,115.7,116.8,99.8,96,115.9,109.1,117.3,109.8,112.8,110.7,100,113.3,122.4,112.5,104.2,92.5,117.2,109.3,106.1,118.8,105.3,106,102,112.9,116.5,114.8,100.5,85.4,114.6,109.9,100.7,115.5,100.7,99,102.3,108.8,105.9,113.2,95.7,80.9,113.9,98.1,102.8,104.7,95.9,94.6,101.6,103.9,110.3,114.1) > 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.74149331 0.01665299 12.84007369 0.97948701 > m$fitted level slope sea Jan 1 89.10000 0.00000000 0.00000000 Feb 1 87.67357 0.07066134 -4.83289900 Mar 1 92.53742 0.44140065 9.88193886 Apr 1 93.72323 0.50674100 -1.96641681 May 1 94.21271 0.50535296 -0.11146086 Jun 1 97.03698 0.66522230 5.86933569 Jul 1 96.95596 0.62048229 -3.68996230 Aug 1 95.25425 0.49227136 -4.04539168 Sep 1 94.37995 0.41862221 0.04303385 Oct 1 95.53010 0.45857839 3.80440430 Nov 1 101.88303 0.79018593 13.29292052 Dec 1 108.68721 1.13961687 7.58180918 Jan 2 109.70691 1.13509571 -9.89617924 Feb 2 109.34826 1.06916975 -13.21452269 Mar 2 108.90754 0.98186748 7.11765741 Apr 2 109.86647 0.98033426 -0.76469517 May 2 112.95038 1.12844121 4.18933416 Jun 2 111.01892 0.91243732 -0.98041636 Jul 2 111.72656 0.89825327 1.08980897 Aug 2 113.23194 0.93929978 -2.58118068 Sep 2 111.21623 0.74377259 -10.97552680 Oct 2 112.16805 0.75729908 1.11502308 Nov 2 113.01151 0.76281737 9.38149827 Dec 2 110.95092 0.58439230 1.77759737 Jan 3 110.98711 0.55003645 -6.74264463 Feb 3 109.42409 0.41568089 -16.75262968 Mar 3 109.36012 0.38431364 7.87838008 Apr 3 109.78621 0.38711677 -0.48950466 May 3 107.05485 0.17459549 -0.71078662 Jun 3 109.94442 0.36063743 8.64266671 Jul 3 109.10941 0.27883122 -3.71495641 Aug 3 108.17161 0.19603842 -2.07493289 Sep 3 109.26352 0.25658202 -7.33483299 Oct 3 110.07705 0.29396300 2.77867847 Nov 3 108.79354 0.18868608 7.83157295 Dec 3 109.31000 0.21048617 5.46401428 Jan 4 108.40811 0.13652119 -7.81982455 Feb 4 106.13078 -0.02449014 -20.53907097 Mar 4 105.36875 -0.07394244 9.28970420 Apr 4 106.19397 -0.01334443 3.63505966 May 4 105.43220 -0.06396280 -4.67331364 Jun 4 105.26180 -0.07117286 10.24657046 Jul 4 104.79228 -0.09815825 -4.06087983 Aug 4 103.68114 -0.16670601 -4.60113640 Sep 4 104.93578 -0.07068671 -2.74810185 Oct 4 105.20383 -0.04784575 3.56942454 Nov 4 102.94097 -0.19696235 3.13376480 Dec 4 103.03390 -0.17746596 10.14323945 Jan 5 102.59463 -0.19507099 -6.87397400 Feb 5 101.95000 -0.22532245 -21.01450678 Mar 5 102.52992 -0.17107652 11.30653352 Apr 5 100.06655 -0.32567063 -1.78585664 May 5 101.27246 -0.22229741 1.40693776 Jun 5 99.43368 -0.33144310 5.39357434 Jul 5 98.81697 -0.35070535 -2.89449844 Aug 5 98.82745 -0.32632198 -4.25591749 Sep 5 100.15075 -0.21499843 1.31919846 Oct 5 100.06828 -0.20605849 3.82127035 Nov 5 102.03858 -0.05931617 8.08996161 Dec 5 103.05453 0.01316357 10.96077116 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.263767108 1.832958204 0.302144510 -0.008391694 1.247925038 2 -0.068571540 -0.854500864 -0.806237121 -0.011627677 1.064588046 -1.584381498 3 -0.292175817 -1.126620065 -0.253127173 0.021760560 -1.615950755 1.412280380 4 -0.584310445 -1.268705617 -0.386825814 0.469928404 -0.390316174 -0.055537964 5 -0.137002837 -0.235361753 0.421389309 -1.198329246 0.800027353 -0.844394940 Jul Aug Sep Oct Nov Dec 1 -0.416890858 -1.308379285 -0.768405531 0.409135792 3.276806207 3.324360517 2 -0.108246170 0.324104625 -1.579883652 0.110936266 0.045790404 -1.499062643 3 -0.626267074 -0.640119773 0.471591316 0.292616068 -0.827295294 0.171907489 4 -0.208301318 -0.530579051 0.744609937 0.177298497 -1.158383765 0.151602761 5 -0.149108534 0.188901504 0.862837597 0.069298464 1.137557470 0.562005064 > 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/1kzss1259930014.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/2bmpy1259930014.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/35rwp1259930014.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/4t2bi1259930014.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/5mz971259930014.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/6kaxn1259930014.tab") > > system("convert tmp/1kzss1259930014.ps tmp/1kzss1259930014.png") > system("convert tmp/2bmpy1259930014.ps tmp/2bmpy1259930014.png") > system("convert tmp/35rwp1259930014.ps tmp/35rwp1259930014.png") > system("convert tmp/4t2bi1259930014.ps tmp/4t2bi1259930014.png") > system("convert tmp/5mz971259930014.ps tmp/5mz971259930014.png") > > > proc.time() user system elapsed 1.624 0.815 3.171