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(611,594,595,591,589,584,573,567,569,621,629,628,612,595,597,593,590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516,528,533,536,537,524,536,587,597,581) > 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.00000 58.48163 23.91491 0.00000 > m$fitted level slope sea Jan 1 611.0000 0.0000000 0.0000000 Feb 1 600.7647 -10.7943517 -6.7647157 Mar 1 590.1307 -10.6922035 4.8692848 Apr 1 587.2285 -5.5162060 3.7714827 May 1 587.5184 -1.6627503 1.4816374 Jun 1 584.8107 -2.3429813 -0.8107208 Jul 1 575.3146 -7.0114014 -2.3145631 Aug 1 566.4669 -8.2129283 0.5331141 Sep 1 565.5634 -3.4288573 3.4365894 Oct 1 604.8745 24.5391522 16.1254680 Nov 1 632.7257 26.7059717 -3.7257324 Dec 1 636.8585 11.9394896 -8.8585444 Jan 2 620.9271 -6.2267084 -8.9270669 Feb 2 603.4536 -13.5910086 -8.4536023 Mar 2 592.9638 -11.6133764 4.0361713 Apr 2 589.2043 -6.6808679 3.7956509 May 2 586.7800 -3.9731719 3.2200459 Jun 2 578.0738 -6.9666548 1.9262390 Jul 2 571.4024 -6.7808841 2.5976270 Aug 2 572.1490 -2.0370971 0.8509789 Sep 2 582.4314 5.7380097 -9.4313692 Oct 2 602.6073 14.8517376 17.3927163 Nov 2 622.0237 17.7315891 3.9762912 Dec 2 625.0515 8.4662158 -5.0515422 Jan 3 603.5719 -10.4094329 -15.5718929 Feb 3 579.2420 -19.1958760 -13.2420063 Mar 3 556.7036 -21.2888492 0.2963683 Apr 3 551.4416 -11.3123835 9.5583729 May 3 542.5338 -9.8081939 6.4662176 Jun 3 529.3814 -11.9012179 2.6186166 Jul 3 522.6325 -8.6883145 3.3674648 Aug 3 515.4665 -7.7395824 -4.4664599 Sep 3 514.0191 -3.8133723 -15.0191047 Oct 3 533.1602 10.5158987 21.8398123 Nov 3 553.7048 16.7738909 11.2951646 Dec 3 545.6488 1.2852169 -3.6488292 Jan 4 538.0580 -4.2560035 -11.0580150 Feb 4 525.6734 -9.3295998 -15.6733715 Mar 4 518.6106 -7.9204370 -4.6106472 Apr 4 508.6331 -9.1957935 8.3668872 May 4 499.1394 -9.3809009 8.8605768 Jun 4 490.1860 -9.1149426 2.8140116 Jul 4 484.4424 -7.0211135 5.5576052 Aug 4 477.4447 -7.0066157 -8.4446528 Sep 4 494.6336 8.0067533 -16.6336023 Oct 4 509.3733 12.1860687 18.6267303 Nov 4 515.9848 8.7258507 18.0152102 Dec 4 519.6303 5.5717839 -1.6303232 Jan 5 517.2777 0.6491838 -11.2777192 Feb 5 518.0752 0.7412213 -16.0751577 Mar 5 519.4942 1.1609682 -3.4942056 Apr 5 519.3071 0.3273849 8.6928753 May 5 521.9542 1.7636126 11.0458411 Jun 5 530.3307 5.8627615 5.6692993 Jul 5 533.7050 4.3215377 3.2949931 Aug 5 541.1090 6.2283074 -17.1089923 Sep 5 553.5778 10.0872894 -17.5777964 Oct 5 566.1260 11.6093156 20.8739888 Nov 5 577.1543 11.2498657 19.8457151 Dec 5 582.0352 7.3081306 -1.0351880 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.794174666 0.013593824 0.724337657 0.499410966 -0.088725081 2 -2.384767528 -0.978848735 0.256319207 0.654112631 0.356670497 -0.389762016 3 -2.476213365 -1.150953156 -0.272843819 1.309148740 0.197665336 -0.273405424 4 -0.725774794 -0.663402146 0.183994200 -0.167015905 -0.024275025 0.034784689 5 -0.644180662 0.012029775 0.054839726 -0.109095347 0.188159321 0.536318992 Jul Aug Sep Oct Nov Dec 1 -0.611703208 -0.157136434 0.625550366 3.657352210 0.283344523 -1.930930388 2 0.024279517 0.620873163 1.016714986 1.191760132 0.376601358 -1.212211759 3 0.419478619 0.124076269 0.513498787 1.873778486 0.818423416 -2.027046589 4 0.273474029 0.001894957 1.963366428 0.546590437 -0.452632450 -0.412767216 5 -0.201389601 0.249182208 0.504581333 0.199082831 -0.047029054 -0.515801063 > 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/1q03b1259938548.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/2tu861259938548.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/3egtz1259938548.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/4jso51259938548.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/58ipg1259938548.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/6uzde1259938548.tab") > system("convert tmp/1q03b1259938548.ps tmp/1q03b1259938548.png") > system("convert tmp/2tu861259938548.ps tmp/2tu861259938548.png") > system("convert tmp/3egtz1259938548.ps tmp/3egtz1259938548.png") > system("convert tmp/4jso51259938548.ps tmp/4jso51259938548.png") > system("convert tmp/58ipg1259938548.ps tmp/58ipg1259938548.png") > > > proc.time() user system elapsed 1.361 0.805 1.574