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(107.11,107.57,107.81,108.75,109.43,109.62,109.54,109.53,109.84,109.67,109.79,109.56,110.22,110.40,110.69,110.72,110.89,110.58,110.94,110.91,111.22,111.09,111.00,111.06,111.55,112.32,112.64,112.36,112.04,112.37,112.59,112.89,113.22,112.85,113.06,112.99,113.32,113.74,113.91,114.52,114.96,114.91,115.30,115.44,115.52,116.08,115.94,115.56,115.88,116.66,117.41,117.68,117.85,118.21,118.92,119.03,119.17,118.95,118.92,118.90) > 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.1044402986 0.0007045387 0.0000000000 0.0000000000 > m$fitted level slope sea Jan 1 107.1100 0.00000000 0.000000000 Feb 1 107.5461 0.02487104 0.023867837 Mar 1 107.7854 0.02842449 0.024643931 Apr 1 108.7221 0.04914888 0.027857758 May 1 109.4000 0.06721414 0.030018180 Jun 1 109.5896 0.07141011 0.030424305 Jul 1 109.5101 0.06545812 0.029943189 Aug 1 109.5003 0.06213582 0.029713973 Sep 1 109.8096 0.07411282 0.030430518 Oct 1 109.6402 0.06137193 0.029761650 Nov 1 109.7601 0.06463506 0.029913350 Dec 1 109.5308 0.04732866 0.029195764 Jan 2 110.2406 0.06955050 -0.020559240 Feb 2 110.3961 0.07602144 0.003891358 Mar 2 110.6859 0.08998910 0.004065154 Apr 2 110.7160 0.08595884 0.004019705 May 2 110.8859 0.09174422 0.004078994 Jun 2 110.5762 0.06351207 0.003815490 Jul 2 110.9360 0.08471450 0.003996049 Aug 2 110.9061 0.07638875 0.003931259 Sep 2 111.2159 0.09355778 0.004053505 Oct 2 111.0861 0.07695198 0.003945210 Nov 2 110.9961 0.06443852 0.003870397 Dec 2 111.0561 0.06410328 0.003868559 Jan 3 111.5084 0.09141112 0.041552169 Feb 3 112.3132 0.14793828 0.006777406 Mar 3 112.6332 0.16114444 0.006799876 Apr 3 112.3533 0.12714824 0.006746697 May 3 112.0333 0.09257088 0.006696963 Jun 3 112.3633 0.11098449 0.006721323 Jul 3 112.5833 0.11946001 0.006731639 Aug 3 112.8833 0.13352574 0.006747391 Sep 3 113.2132 0.14886014 0.006763195 Oct 3 112.8433 0.10830304 0.006724720 Nov 3 113.0533 0.11626244 0.006731671 Dec 3 112.9833 0.10166859 0.006719938 Jan 4 113.3550 0.12240231 -0.035048343 Feb 4 113.7340 0.14267407 0.005962595 Mar 4 113.9040 0.14482013 0.005961911 Apr 4 114.5140 0.18137392 0.005951268 May 4 114.9541 0.20170650 0.005945817 Jun 4 114.9040 0.18190986 0.005950705 Jul 4 115.2941 0.19828174 0.005946982 Aug 4 115.4341 0.19369497 0.005947942 Sep 4 115.5141 0.18474495 0.005949669 Oct 4 116.0741 0.21429111 0.005944419 Nov 4 115.9341 0.18639057 0.005948986 Dec 4 115.5540 0.14178036 0.005955712 Jan 5 115.9212 0.15938818 -0.041218307 Feb 5 116.6514 0.20431920 0.008627679 Mar 5 117.4014 0.24731489 0.008596959 Apr 5 117.6714 0.24910244 0.008595783 May 5 117.8414 0.24286893 0.008599560 Jun 5 118.2014 0.25209968 0.008594407 Jul 5 118.9114 0.28818699 0.008575851 Aug 5 119.0214 0.27414349 0.008582503 Sep 5 119.1614 0.26357087 0.008587116 Oct 5 118.9414 0.22545688 0.008602435 Nov 5 118.9114 0.20532187 0.008609889 Dec 5 118.8914 0.18756176 0.008615946 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.80932358 0.66015830 2.78846270 1.92390559 0.37339110 2 2.37722696 0.21972005 0.64014432 -0.17928192 0.25094016 -1.19865192 3 1.26395209 1.93598588 0.51157497 -1.31139232 -1.32904900 0.70561918 4 0.85279371 0.71485211 0.08116714 1.38170618 0.76817153 -0.74760598 5 0.70105560 1.61318169 1.62062785 0.06737272 -0.23492655 0.34786711 Jul Aug Sep Oct Nov Dec 1 -0.45921109 -0.22829395 0.74821258 -0.73538807 0.17629344 -0.88451569 2 0.88401803 -0.34184007 0.69579901 -0.66557291 -0.49686630 -0.01320537 3 0.32395113 0.53644623 0.58374853 -1.54149695 0.30211811 -0.55331700 4 0.61804975 -0.17310078 -0.33767828 1.11451135 -1.05224017 -1.68216200 5 1.35991387 -0.52919506 -0.39839035 -1.43614684 -0.75867521 -0.66917746 > 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/13k0i1259954634.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/2dj2g1259954634.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/37jt31259954634.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/47x4g1259954634.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/5sowl1259954634.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/6g3ub1259954634.tab") > system("convert tmp/13k0i1259954634.ps tmp/13k0i1259954634.png") > system("convert tmp/2dj2g1259954634.ps tmp/2dj2g1259954634.png") > system("convert tmp/37jt31259954634.ps tmp/37jt31259954634.png") > system("convert tmp/47x4g1259954634.ps tmp/47x4g1259954634.png") > system("convert tmp/5sowl1259954634.ps tmp/5sowl1259954634.png") > > > proc.time() user system elapsed 1.536 0.806 1.762