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(102.1,102.86,102.99,103.73,105.02,104.43,104.63,104.93,105.87,105.66,106.76,106,107.22,107.33,107.11,108.86,107.72,107.88,108.38,107.72,108.41,109.9,111.45,112.18,113.34,113.46,114.06,115.54,116.39,115.94,116.97,115.94,115.91,116.43,116.26,116.35,117.9,117.7,117.53,117.86,117.65,116.51,115.93,115.31,115) > 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.249261459 0.031903323 0.007957765 0.054894428 > m$fitted level slope sea Jan 1 102.1000 0.00000000 0.0000000000 Feb 1 102.7652 0.06421066 0.0472981722 Mar 1 102.9373 0.08308504 0.0375618300 Apr 1 103.5880 0.21195164 0.0677148641 May 1 104.8127 0.47019880 0.0812581721 Jun 1 104.5406 0.27020506 -0.0205258726 Jul 1 104.6134 0.21560744 0.0402172398 Aug 1 104.8801 0.22993291 0.0437790001 Sep 1 105.7275 0.40405037 0.0692493145 Oct 1 105.7060 0.28366325 0.0044316007 Nov 1 106.6092 0.45917421 0.0775061607 Dec 1 106.1331 0.19403653 -0.0224910439 Jan 2 107.1251 0.41119476 0.0045184293 Feb 2 107.3428 0.35872162 0.0038212067 Mar 2 107.2307 0.22769268 -0.0680454006 Apr 2 108.6425 0.56135759 0.0835838856 May 2 107.9338 0.20350971 -0.0704465515 Jun 2 107.9342 0.14630821 -0.0312448269 Jul 2 108.3198 0.21368630 0.0330781990 Aug 2 107.8890 0.03225139 -0.0961080205 Sep 2 108.2654 0.12912405 0.1056624496 Oct 2 109.7143 0.50071545 0.0362588532 Nov 2 111.1936 0.77625484 0.1455878915 Dec 2 112.2523 0.85567194 -0.1042037934 Jan 3 113.2248 0.88808997 0.1020490692 Feb 3 113.5202 0.72514629 -0.0027609636 Mar 3 114.2448 0.72499699 -0.1847621137 Apr 3 115.2857 0.81351990 0.2193995386 May 3 116.3980 0.89729213 -0.0409189243 Jun 3 116.2516 0.60480141 -0.1965828917 Jul 3 116.8445 0.60148601 0.1267798210 Aug 3 116.3004 0.28052278 -0.2340656989 Sep 3 116.0036 0.11878067 -0.0299736904 Oct 3 116.3995 0.19642493 -0.0000953093 Nov 3 116.2379 0.09611013 0.0615700667 Dec 3 116.4358 0.12457156 -0.0970117213 Jan 4 117.4178 0.36319875 0.3871266084 Feb 4 117.7276 0.34841902 -0.0221249495 Mar 4 117.8269 0.27947210 -0.2702318671 Apr 4 117.7733 0.18648787 0.1229697933 May 4 117.6634 0.10363854 0.0188137011 Jun 4 116.9230 -0.13216370 -0.3213564496 Jul 4 115.8821 -0.38600464 0.1466623013 Aug 4 115.4668 -0.39418284 -0.1536213097 Sep 4 115.0654 -0.39620421 -0.0645984849 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.811371258 0.189765895 0.954909891 1.654878582 -1.195593663 2 1.456673280 -0.260008449 -0.739681293 1.865019835 -1.997471269 -0.319764662 3 0.196644176 -0.857876759 -0.000839676 0.495729607 0.467859691 -1.635287535 4 1.396611504 -0.079848440 -0.386690011 -0.521370824 -0.462937653 -1.318368019 Jul Aug Sep Oct Nov Dec 1 -0.315576160 0.081455709 0.982188463 -0.676473794 0.984369587 -1.485690970 2 0.376964779 -1.015458587 0.542271328 2.080269909 1.542610262 0.444689781 3 -0.018550613 -1.796485186 -0.905431962 0.434686592 -0.561600289 0.159495425 4 -1.420320711 -0.045775821 -0.011315908 > 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/1ye721260548568.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/25o8q1260548568.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/3bxdp1260548568.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/4cq0p1260548568.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/51e4a1260548568.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/63rjw1260548568.tab") > > system("convert tmp/1ye721260548568.ps tmp/1ye721260548568.png") > system("convert tmp/25o8q1260548568.ps tmp/25o8q1260548568.png") > system("convert tmp/3bxdp1260548568.ps tmp/3bxdp1260548568.png") > system("convert tmp/4cq0p1260548568.ps tmp/4cq0p1260548568.png") > system("convert tmp/51e4a1260548568.ps tmp/51e4a1260548568.png") > > > proc.time() user system elapsed 1.306 0.770 1.547