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(105.7,105.7,111.1,82.4,60,107.3,99.3,113.5,108.9,100.2,103.9,138.7,120.2,100.2,143.2,70.9,85.2,133,136.6,117.9,106.3,122.3,125.5,148.4,126.3,99.6,140.4,80.3,92.6,138.5,110.9,119.6,105,109,129.4,148.6,101.4,134.8,143.7,81.6,90.3,141.5,140.7,140.2,100.2,125.7,119.6,134.7,109,116.3,146.9,97.4,89.4,132.1,139.8,129,112.5,121.9,121.7,123.1) > 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 10.16448 0.00000 63.68676 53.97844 > m$fitted level slope sea Jan 1 105.70000 0.000000000 0.00000000 Feb 1 105.70000 0.000000000 0.00000000 Mar 1 107.34464 0.192645853 2.61208369 Apr 1 100.32690 -0.600075382 -13.37284982 May 1 87.51743 -1.750536641 -18.69116218 Jun 1 89.05057 -1.490247603 15.56053539 Jul 1 91.55106 -1.221516121 4.19590289 Aug 1 97.08037 -0.826493598 10.08425502 Sep 1 101.10989 -0.573229050 3.07978723 Oct 1 101.86393 -0.510117977 -2.97895842 Nov 1 102.40103 -0.463902980 0.44666367 Dec 1 110.63200 -0.103094588 19.24495240 Jan 2 114.05954 -0.190180561 1.92359980 Feb 2 111.89767 -0.207307663 -9.64748942 Mar 2 115.86934 -0.067403366 23.63028327 Apr 2 107.83305 -0.408816135 -30.27364121 May 2 105.32484 -0.501441986 -18.35007022 Jun 2 107.97319 -0.369766808 22.27070695 Jul 2 114.96668 -0.087905668 14.96599364 Aug 2 116.09664 -0.045663223 0.67023297 Sep 2 114.16867 -0.104755625 -6.08193441 Oct 2 116.67665 -0.030951197 3.10592738 Nov 2 120.69218 0.068516431 0.85846433 Dec 2 124.14275 0.132392438 20.90296799 Jan 3 125.06103 0.138413571 0.43181112 Feb 3 122.45253 0.105562095 -20.10754199 Mar 3 119.34456 0.041453000 24.10722603 Apr 3 116.76470 -0.025361998 -34.05590297 May 3 115.88938 -0.049086167 -22.51428692 Jun 3 116.97628 -0.017180040 20.47913334 Jul 3 113.04812 -0.122538945 1.50083198 Aug 3 113.27521 -0.113758542 5.99364190 Sep 3 113.68365 -0.101778183 -9.18479299 Oct 3 113.09102 -0.111845558 -3.61473435 Nov 3 116.49050 -0.049563073 9.46964458 Dec 3 119.33368 -0.007401578 26.40636849 Jan 4 115.55059 -0.051324948 -10.38797560 Feb 4 121.68390 0.026279953 7.00925803 Mar 4 122.33291 0.035917884 20.76389462 Apr 4 121.01015 0.011101003 -38.11509489 May 4 118.92750 -0.030747810 -26.64550869 Jun 4 118.21959 -0.044669813 23.92234173 Jul 4 122.44277 0.041623218 14.18564310 Aug 4 126.11693 0.111547949 10.58897320 Sep 4 124.02750 0.072211893 -21.69291231 Oct 4 125.43427 0.093882278 -1.03816030 Nov 4 123.20599 0.060267091 -1.32286004 Dec 4 119.54275 0.012600880 18.83617649 Jan 5 119.60085 0.013132937 -10.64590836 Feb 5 117.46792 -0.012476650 0.94841017 Mar 5 118.32576 -0.001048403 27.72304770 Apr 5 121.81054 0.049544935 -27.79298744 May 5 121.86791 0.049666614 -32.47547218 Jun 5 120.06696 0.020029236 13.82064734 Jul 5 121.16196 0.037155975 17.59665877 Aug 5 120.45676 0.025723074 9.26568050 Sep 5 122.93494 0.061407085 -12.83357260 Oct 5 122.88047 0.059841215 -0.86659187 Nov 5 122.28289 0.051667344 0.06574196 Dec 5 118.51542 0.007761493 8.36315542 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.000000000 0.338190942 -1.422402444 -2.697476213 0.800690552 2 1.238239944 -0.605839568 1.125174541 -2.063844497 -0.552217588 0.854582903 3 0.243563898 -0.830383279 -0.931303001 -0.741379347 -0.239630628 0.323196170 4 -1.146895593 1.862918490 0.184605385 -0.397874961 -0.610582233 -0.197972368 5 0.013784512 -0.647659949 0.260848333 1.038437644 0.002325309 -0.550063540 Jul Aug Sep Oct Nov Dec 1 1.039067822 1.831629861 1.352011872 0.375707425 0.299715813 2.507556963 2 2.056481918 0.347822953 -0.546295481 0.767090685 1.199182747 1.014258152 3 -1.127443570 0.102086436 0.154110378 -0.146100934 1.052584956 0.873169679 4 1.255374096 1.076166209 -0.656471340 0.400356465 -0.699949345 -1.126413863 5 0.320454294 -0.222191855 0.737088651 -0.034958080 -0.198935307 -1.158002376 > 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/1ka3f1259698791.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/2xaer1259698791.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/3a86c1259698791.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/404bl1259698791.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/5cd5v1259698791.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/69ngy1259698791.tab") > system("convert tmp/1ka3f1259698791.ps tmp/1ka3f1259698791.png") > system("convert tmp/2xaer1259698791.ps tmp/2xaer1259698791.png") > system("convert tmp/3a86c1259698791.ps tmp/3a86c1259698791.png") > system("convert tmp/404bl1259698791.ps tmp/404bl1259698791.png") > system("convert tmp/5cd5v1259698791.ps tmp/5cd5v1259698791.png") > > > proc.time() user system elapsed 1.885 0.815 2.124