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(101.3,106.3,94,102.8,102,105.1,92.4,81.4,105.8,120.3,100.7,88.8,94.3,99.9,103.4,103.3,98.8,104.2,91.2,74.7,108.5,114.5,96.9,89.6,97.1,100.3,122.6,115.4,109,129.1,102.8,96.2,127.7,128.9,126.5,119.8,113.2,114.1,134.1,130,121.8,132.1,105.3,103,117.1,126.3,138.1,119.5,138,135.5,178.6,162.2,176.9,204.9,132.2,142.5,164.3,174.9,175.4,143) > par20 = '' > par19 = '' > par18 = '' > par17 = '' > par16 = '' > par15 = '' > par14 = '' > par13 = '' > par12 = '' > par11 = '' > par10 = '' > par9 = '' > par8 = '' > par7 = '' > par6 = '' > par5 = '' > par4 = '' > par3 = '' > par2 = '' > par1 = '12' > ylab = '' > xlab = '' > main = '' > #'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 26.88559263 0.03210685 31.98526634 0.00000000 > m$fitted level slope sea Jan 1 101.30000 0.00000000 0.0000000 Feb 1 104.21401 0.09500628 2.0859935 Mar 1 99.39967 -0.20792223 -5.3996671 Apr 1 100.19656 -0.15734224 2.6034413 May 1 101.19975 -0.11920145 0.8002527 Jun 1 103.22354 -0.07112923 1.8764634 Jul 1 98.82673 -0.15280490 -6.4267284 Aug 1 89.82725 -0.32020953 -8.4272511 Sep 1 94.89651 -0.21177249 10.9034906 Oct 1 107.75525 0.06726744 12.5447531 Nov 1 107.73904 0.06540496 -7.0390443 Dec 1 99.11365 -0.13493099 -10.3136548 Jan 2 95.15931 -0.10082135 -0.8593102 Feb 2 94.29229 -0.10743062 5.6077083 Mar 2 101.25760 0.10180991 2.1424032 Apr 2 102.81274 0.15563304 0.4872642 May 2 100.84288 0.08006969 -2.0428814 Jun 2 99.23318 0.02727209 4.9668240 Jul 2 96.65690 -0.04519534 -5.4569042 Aug 2 92.54885 -0.15172334 -17.8488505 Sep 2 96.81610 -0.03768460 11.6839038 Oct 2 99.06817 0.02076947 15.4318334 Nov 2 99.19890 0.02344719 -2.2989048 Dec 2 98.72955 0.01269026 -9.1295493 Jan 3 98.94329 0.01659020 -1.8432942 Feb 3 99.04431 0.01851105 1.2556889 Mar 3 108.78226 0.29927202 13.8177359 Apr 3 112.81835 0.42325607 2.5816504 May 3 112.16100 0.38608549 -3.1610031 Jun 3 116.50253 0.51871644 12.5974686 Jul 3 113.39328 0.40267233 -10.5932772 Aug 3 114.42459 0.42191680 -18.2245878 Sep 3 115.89987 0.45306758 11.8001282 Oct 3 115.08153 0.41662603 13.8184660 Nov 3 120.89650 0.56629785 5.6035000 Dec 3 126.09489 0.69124434 -6.2948937 Jan 4 123.66694 0.60684762 -10.4669368 Feb 4 121.49417 0.52752566 -7.3941700 Mar 4 121.69358 0.51741352 12.4064189 Apr 4 124.48418 0.59189786 5.5158170 May 4 126.01780 0.62365598 -4.2178003 Jun 4 122.78296 0.49358656 9.3170434 Jul 4 119.54386 0.37002984 -14.2438650 Aug 4 120.14976 0.37764384 -17.1497561 Sep 4 114.28842 0.18128037 2.8115817 Oct 4 114.12966 0.17081840 12.1703439 Nov 4 122.99254 0.43360121 15.1074601 Dec 4 125.01463 0.48126509 -5.5146253 Jan 5 135.19708 0.77469201 2.8029195 Feb 5 141.00530 0.93063402 -5.5052966 Mar 5 154.10590 1.31980626 24.4941012 Apr 5 157.96264 1.40321127 4.2373559 May 5 167.99553 1.69152202 8.9044662 Jun 5 181.47668 2.08659464 23.4233163 Jul 5 168.67951 1.59142105 -36.4795082 Aug 5 161.39681 1.30000796 -18.8968137 Sep 5 161.43968 1.25933029 2.8603172 Oct 5 165.09643 1.33588073 9.8035727 Nov 5 165.71891 1.31330901 9.6810936 Dec 5 162.03659 1.15568239 -19.0365930 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.44744706 -0.68710802 0.17534576 0.21949869 0.41090127 2 -0.76814114 -0.14941446 1.25307160 0.25692468 -0.38973172 -0.31665534 3 0.03809933 0.01586214 1.77836717 0.67650663 -0.19759575 0.73288693 4 -0.58213922 -0.51634506 -0.06033471 0.41552487 0.17262976 -0.71215549 5 1.79917344 0.93117236 2.24070444 0.46562304 1.58528663 2.17334990 Jul Aug Sep Oct Nov Dec 1 -0.82928387 -1.69234534 1.02868501 2.49028569 -0.01588194 -1.65182397 2 -0.49100639 -0.76664146 0.83305658 0.43098761 0.02065756 -0.09262582 3 -0.67696425 0.11749955 0.19671328 -0.23702144 1.00525025 0.86337235 4 -0.69238805 0.04381539 -1.15798089 -0.06302408 1.61023620 0.29450143 5 -2.75211569 -1.64250237 -0.23252281 0.44298804 -0.13179343 -0.92343210 > 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/1bwy41259777101.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/271in1259777101.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/3fam11259777101.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/4wqn01259777101.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/5sgiq1259777101.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/6ubrc1259777101.tab") > system("convert tmp/1bwy41259777101.ps tmp/1bwy41259777101.png") > system("convert tmp/271in1259777101.ps tmp/271in1259777101.png") > system("convert tmp/3fam11259777101.ps tmp/3fam11259777101.png") > system("convert tmp/4wqn01259777101.ps tmp/4wqn01259777101.png") > system("convert tmp/5sgiq1259777101.ps tmp/5sgiq1259777101.png") > > > proc.time() user system elapsed 1.533 0.806 1.792