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(96.2,96.8,109.9,88,91.1,106.4,68.6,100.1,108,106,108.6,91.5,99.2,98,96.6,102.8,96.9,110,70.5,101.9,109.6,107.8,113,93.8,108,102.8,116.3,89.2,106.7,112.1,74.2,108.8,111.5,118.8,118.9,97.6,116.4,107.9,121.2,97.9,113.4,117.6,79.6,115.9,115.7,129.1,123.3,96.7,121.2,118.2,102.1,125.4,116.7,121.3,85.3,114.2,124.4,131,118.3,99.6) > 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.000000000 0.002781178 26.255266265 11.932067941 > m$fitted level slope sea Jan 1 96.20000 0.00000000 0.0000000 Feb 1 96.16187 -0.01268651 0.4909717 Mar 1 97.85650 0.41559317 9.0401322 Apr 1 97.41626 0.24387164 -8.7766611 May 1 95.91796 -0.04782976 -3.6495059 Jun 1 97.40005 0.17240903 7.8266554 Jul 1 93.07611 -0.39705911 -20.3647863 Aug 1 92.49152 -0.41835554 7.8103771 Sep 1 94.24363 -0.19369015 11.0653809 Oct 1 96.30972 0.02267565 6.5308608 Nov 1 98.55895 0.22254801 6.5934445 Dec 1 98.79600 0.22378370 -7.3204845 Jan 2 99.19466 0.23879406 -0.3613112 Feb 2 99.50639 0.24484112 -1.6566646 Mar 2 98.53118 0.14817655 0.4296296 Apr 2 99.17737 0.18616227 2.6836842 May 2 99.60747 0.20428212 -3.1710576 Jun 2 99.80588 0.20385309 10.2055444 Jul 2 99.38958 0.15872488 -27.6434549 Aug 2 99.21491 0.13441579 3.3745734 Sep 2 99.47484 0.14363602 9.8592392 Oct 2 100.01934 0.17341834 6.9156351 Nov 2 100.84288 0.22241447 10.7366335 Dec 2 101.33264 0.24291398 -8.1209634 Jan 3 102.16604 0.28917337 4.5293092 Feb 3 102.74602 0.31214077 -0.5889073 Mar 3 104.25177 0.40642592 9.4198819 Apr 3 103.74707 0.33460598 -12.5487334 May 3 104.15694 0.34052928 2.3783543 Jun 3 104.26255 0.32204412 8.3510411 Jul 3 104.29293 0.29906891 -29.4558214 Aug 3 104.63165 0.30219827 4.0818695 Sep 3 104.81966 0.29316649 6.9287696 Oct 3 105.74122 0.34297801 11.6964125 Nov 3 106.60205 0.38410653 11.1798386 Dec 3 107.22291 0.40293825 -10.1319977 Jan 4 108.15006 0.44467167 7.1262473 Feb 4 108.89072 0.46823300 -1.6244948 Mar 4 109.56041 0.48425793 11.2083084 Apr 4 110.16342 0.49369550 -12.5176636 May 4 110.71456 0.49825693 2.5624343 Jun 4 111.08430 0.48805706 6.7910002 Jul 4 111.37963 0.47276692 -31.3667181 Aug 4 111.84851 0.47245919 4.0597974 Sep 4 112.18595 0.46175170 3.8031994 Oct 4 113.08797 0.49665980 15.0697169 Nov 4 113.72749 0.50798492 9.2669313 Dec 4 113.80570 0.47392412 -16.1868956 Jan 5 114.15046 0.46368984 7.3256297 Feb 5 115.00373 0.49455246 2.3634591 Mar 5 113.70430 0.35244254 -7.7683818 Apr 5 115.23150 0.44550125 7.6561466 May 5 115.97128 0.46881352 0.0992920 Jun 5 116.38284 0.46427847 5.0396089 Jul 5 116.81600 0.46181293 -31.4494268 Aug 5 116.70353 0.41631273 -1.2748836 Sep 5 117.20753 0.42326040 7.0048579 Oct 5 117.49949 0.41285713 13.7814226 Nov 5 117.07865 0.34680751 3.0047833 Dec 5 116.92722 0.30733591 -16.2614276 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.086021577 1.883045736 -0.467456494 -0.896635600 0.887236301 2 0.251425388 0.103254253 -1.641363653 0.658413846 0.326128408 -0.008038344 3 0.910942360 0.448987937 1.836541063 -1.397461854 0.115256056 -0.359528535 4 0.788577828 0.444878282 0.302742569 0.178470243 0.086350324 -0.193258835 5 -0.193895447 0.584900466 -2.693984277 1.764411853 0.442045514 -0.085998538 Jul Aug Sep Oct Nov Dec 1 -3.025205311 -0.144907069 1.892894851 2.187181073 2.356569741 0.016569223 2 -0.874814264 -0.482935014 0.185851047 0.603463025 0.989712452 0.409803892 3 -0.446067606 0.060556966 -0.173977528 0.954345355 0.783587310 0.357005753 4 -0.289865667 -0.005835014 -0.202995190 0.661574160 0.214560501 -0.645210221 5 -0.046756366 -0.862879213 0.131756650 -0.197284711 -1.252531146 -0.748532767 > 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/1q3u61259957855.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/2eh9m1259957855.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/3zuf51259957855.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/4s9mk1259957855.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/5rp4k1259957855.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/66q4e1259957856.tab") > > system("convert tmp/1q3u61259957855.ps tmp/1q3u61259957855.png") > system("convert tmp/2eh9m1259957855.ps tmp/2eh9m1259957855.png") > system("convert tmp/3zuf51259957855.ps tmp/3zuf51259957855.png") > system("convert tmp/4s9mk1259957855.ps tmp/4s9mk1259957855.png") > system("convert tmp/5rp4k1259957855.ps tmp/5rp4k1259957855.png") > > > proc.time() user system elapsed 1.441 0.807 1.681