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(126.51,131.02,136.51,138.04,132.92,129.61,122.96,124.04,121.29,124.56,118.53,113.14,114.15,122.17,129.23,131.19,129.12,128.28,126.83,138.13,140.52,146.83,135.14,131.84,125.7,128.98,133.25,136.76,133.24,128.54,121.08,120.23,119.08,125.75,126.89,126.6,121.89,123.44,126.46,129.49,127.78,125.29,119.02,119.96,122.86,131.89,132.73,135.01,136.71,142.73,144.43,144.93,138.75,130.22,122.19,128.4,140.43,153.5,149.33,142.97) > 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.9586113 21.0457459 0.2686388 0.0000000 > m$fitted level slope sea Jan 1 126.5100 0.00000000 0.00000000 Feb 1 130.9172 3.92246240 0.10279744 Mar 1 136.4751 5.42119633 0.03494923 Apr 1 138.1733 2.03974118 -0.13330136 May 1 133.1300 -4.40470662 -0.20997412 Jun 1 129.4802 -3.71771194 0.12977445 Jul 1 123.0452 -6.19027273 -0.08518392 Aug 1 123.7168 0.05352586 0.32317230 Sep 1 121.4229 -2.08245193 -0.13293068 Oct 1 124.3141 2.44319451 0.24587286 Nov 1 118.8599 -4.74289200 -0.32986962 Dec 1 113.0721 -5.69368359 0.06793481 Jan 2 113.8201 0.13729392 0.32994304 Feb 2 121.9486 6.93900047 0.22137293 Mar 2 129.2511 7.26134584 -0.02113368 Apr 2 131.2271 2.59179533 -0.03714027 May 2 129.6101 -1.13022731 -0.49013523 Jun 2 127.7620 -1.76540834 0.51796204 Jul 2 127.2216 -0.68158133 -0.39164831 Aug 2 136.8915 8.47639143 1.23846731 Sep 2 141.4961 5.05100181 -0.97612579 Oct 2 146.1820 4.72798641 0.64797096 Nov 2 136.0553 -8.41550491 -0.91534365 Dec 2 131.4997 -5.00325928 0.34031165 Jan 3 126.3610 -5.12301347 -0.66099773 Feb 3 128.6209 1.19443232 0.35906713 Mar 3 132.8519 3.84495003 0.39808179 Apr 3 136.5475 3.71454246 0.21249725 May 3 133.9654 -1.78161479 -0.72540995 Jun 3 127.9117 -5.51375707 0.62828248 Jul 3 122.3107 -5.59002397 -1.23065990 Aug 3 118.5839 -3.96237669 1.64608938 Sep 3 120.2041 0.91422982 -1.12413179 Oct 3 123.6345 3.11237491 2.11547870 Nov 3 128.0255 4.22934776 -1.13552432 Dec 3 126.3985 -0.88035703 0.20149494 Jan 4 123.4092 -2.72377114 -1.51915706 Feb 4 123.2442 -0.52913109 0.19582175 Mar 4 126.0859 2.39064071 0.37410155 Apr 4 128.8503 2.71512018 0.63966718 May 4 128.3411 -0.08106188 -0.56111501 Jun 4 124.7515 -3.12645945 0.53847054 Jul 4 120.2159 -4.34972999 -1.19594433 Aug 4 118.4721 -2.08770906 1.48789077 Sep 4 123.7085 4.26988405 -0.84848127 Oct 4 130.0478 6.06638992 1.84217731 Nov 4 133.5840 3.87038929 -0.85396790 Dec 4 134.6876 1.47156316 0.32237101 Jan 5 138.2255 3.26604773 -1.51545829 Feb 5 142.8027 4.39140545 -0.07273961 Mar 5 144.4563 2.02957335 -0.02625677 Apr 5 144.2974 0.13636580 0.63259491 May 5 139.3353 -4.26872744 -0.58533954 Jun 5 129.6752 -8.93054894 0.54479657 Jul 5 123.1174 -6.87835320 -0.92744095 Aug 5 126.8472 2.29615983 1.55275343 Sep 5 140.8544 12.42450911 -0.42442506 Oct 5 151.6596 11.02398526 1.84039670 Nov 5 150.6606 0.62962919 -1.33059752 Dec 5 143.3603 -6.22189288 -0.39033951 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.91384762 0.32322854 -0.73530100 -1.40484470 0.14975122 2 1.30617669 1.47143578 0.07029562 -1.01410309 -0.81138113 -0.13845531 3 -0.02644075 1.36993223 0.57956955 -0.02832662 -1.19793769 -0.81353046 4 -0.40421114 0.47676819 0.63857087 0.07053901 -0.60930773 -0.66385220 5 0.39223653 0.24475163 -0.51629860 -0.41190698 -0.95970999 -1.01622860 Jul Aug Sep Oct Nov Dec 1 -0.53897032 1.36102693 -0.46560172 0.98650312 -1.56642743 -0.20725412 2 0.23625334 1.99626005 -0.74666857 -0.07041111 -2.86507464 0.74403147 3 -0.01662463 0.35479558 1.06300566 0.47915650 0.24347986 -1.11480801 4 -0.26664566 0.49307669 1.38583281 0.39161036 -0.47867938 -0.52354422 5 0.44733094 1.99986097 2.20779437 -0.30529262 -2.26577834 -1.49545917 > 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/1arjq1259956683.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/2xdcr1259956683.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/3gznf1259956683.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/4ye6o1259956683.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/5hs8d1259956683.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/62jbw1259956683.tab") > system("convert tmp/1arjq1259956683.ps tmp/1arjq1259956683.png") > system("convert tmp/2xdcr1259956683.ps tmp/2xdcr1259956683.png") > system("convert tmp/3gznf1259956683.ps tmp/3gznf1259956683.png") > system("convert tmp/4ye6o1259956683.ps tmp/4ye6o1259956683.png") > system("convert tmp/5hs8d1259956683.ps tmp/5hs8d1259956683.png") > > > proc.time() user system elapsed 1.350 0.806 1.568