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(115.6,111.9,107,107.1,100.6,99.2,108.4,103,99.8,115,90.8,95.9,114.4,108.2,112.6,109.1,105,105,118.5,103.7,112.5,116.6,96.6,101.9,116.5,119.3,115.4,108.5,111.5,108.8,121.8,109.6,112.2,119.6,104.1,105.3,115,124.1,116.8,107.5,115.6,116.2,116.3,119,111.9,118.6,106.9,103.2,118.6,118.7,102.8,100.6,94.9,94.5,102.9,95.3,92.5,102.7,91.5,89.5) > 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 5.554876 0.000000 8.580506 0.000000 > m$fitted level slope sea Jan 1 115.60000 0.00000000 0.000000000 Feb 1 113.68774 -0.04718809 -1.787741698 Mar 1 109.65065 -0.30042027 -2.650650696 Apr 1 107.81449 -0.38559516 -0.714489771 May 1 104.41920 -0.49698261 -3.819200031 Jun 1 101.41984 -0.55526911 -2.219839936 Jul 1 103.56624 -0.51038253 4.833758983 Aug 1 103.96558 -0.49721243 -0.965582785 Sep 1 102.39022 -0.51265654 -2.590220127 Oct 1 107.27461 -0.43354734 7.725390358 Nov 1 101.77290 -0.50896689 -10.972904482 Dec 1 97.88724 -0.55923678 -1.987241232 Jan 2 101.58600 -0.64449499 12.813996370 Feb 2 104.01411 -0.66373284 4.185887593 Mar 2 108.47397 -0.57972380 4.126030219 Apr 2 109.09196 -0.54780293 0.008037241 May 2 108.69617 -0.54376372 -3.696170420 Jun 2 108.61043 -0.53376396 -3.610425413 Jul 2 110.70853 -0.48934279 7.791469072 Aug 2 108.72351 -0.50955074 -5.023510442 Sep 2 111.24077 -0.47454980 1.259226740 Oct 2 109.20800 -0.49008361 7.391998563 Nov 2 107.51422 -0.49912716 -10.914220213 Dec 2 106.90678 -0.49952062 -5.006780163 Jan 3 106.61602 -0.49954203 9.883976417 Feb 3 111.16744 -0.48596159 8.132562046 Mar 3 112.46017 -0.46852601 2.939830831 Apr 3 111.14918 -0.48186070 -2.649181882 May 3 112.73540 -0.44400712 -1.235400796 Jun 3 113.36596 -0.42504979 -4.565961532 Jul 3 113.46820 -0.41695729 8.331800905 Aug 3 114.22879 -0.40201039 -4.628788659 Sep 3 112.35694 -0.41702042 -0.156940499 Oct 3 111.39449 -0.42132131 8.205507133 Nov 3 112.60461 -0.41212828 -8.504607149 Dec 3 112.43082 -0.41122263 -7.130824121 Jan 4 110.70651 -0.41554748 4.293485712 Feb 4 112.71636 -0.40346227 11.383637410 Mar 4 113.57362 -0.39301545 3.226377404 Apr 4 112.88329 -0.39646318 -5.383286112 May 4 114.52494 -0.36871405 1.075062246 Jun 4 117.41743 -0.32319078 -1.217426363 Jul 4 114.21968 -0.36042445 2.080317923 Aug 4 116.73200 -0.32831643 2.268000141 Sep 4 115.31590 -0.33824443 -3.415896545 Oct 4 113.26903 -0.35043228 5.330973454 Nov 4 113.54331 -0.34701645 -6.643308158 Dec 4 112.19596 -0.35147874 -8.995956955 Jan 5 113.34095 -0.34480762 5.259050555 Feb 5 111.44486 -0.35345910 7.255141217 Mar 5 106.39406 -0.38855206 -3.594057356 Apr 5 105.77669 -0.39071329 -5.176692232 May 5 101.30768 -0.43491133 -6.407677568 Jun 5 97.55239 -0.47242260 -3.052394679 Jul 5 98.44829 -0.45759839 4.451712023 Aug 5 95.95595 -0.47732768 -0.655953101 Sep 5 94.90264 -0.48206375 -2.402643990 Oct 5 95.59664 -0.47413209 7.103359523 Nov 5 96.40899 -0.46698332 -4.908992980 Dec 5 97.20591 -0.46077845 -7.705910033 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.69461689 -1.17986137 -0.56690834 -1.23743256 -1.05663906 2 1.89998441 1.34566829 2.03765012 0.46774186 0.06140922 0.19050085 3 0.08904262 2.13501456 0.73235660 -0.34176275 0.84462055 0.44560441 4 -0.55521031 1.01863608 0.52352510 -0.12252193 0.84094118 1.35571372 5 0.63097254 -0.65103652 -1.95933031 -0.09504001 -1.69421195 -1.38471498 Jul Aug Sep Oct Nov Dec 1 1.14565085 0.38565675 -0.45651788 2.28293494 -2.14266126 -1.42734437 2 1.10937592 -0.63333748 1.28309305 -0.66028082 -0.50958268 -0.04591203 3 0.22131362 0.49738291 -0.62223094 -0.23092266 0.69047279 0.10088996 4 -1.20471851 1.21057250 -0.45955430 -0.72232744 0.26408834 -0.42265051 5 0.57352587 -0.85629922 -0.24295161 0.49652024 0.54319170 0.53331266 > 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/10ku41259861988.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/2pii41259861988.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/32itg1259861988.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/40r8n1259861988.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/52pst1259861988.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/6oey11259861988.tab") > > system("convert tmp/10ku41259861988.ps tmp/10ku41259861988.png") > system("convert tmp/2pii41259861988.ps tmp/2pii41259861988.png") > system("convert tmp/32itg1259861988.ps tmp/32itg1259861988.png") > system("convert tmp/40r8n1259861988.ps tmp/40r8n1259861988.png") > system("convert tmp/52pst1259861988.ps tmp/52pst1259861988.png") > > > proc.time() user system elapsed 1.412 0.817 1.787