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,564) > 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 57.31257 24.11027 0.00000 > m$fitted level slope sea Jan 1 611.0000 0.0000000 0.0000000 Feb 1 600.8901 -10.6793508 -6.8900770 Mar 1 590.1329 -10.7287539 4.8670728 Apr 1 587.1586 -5.6019979 3.8414146 May 1 587.4580 -1.6994334 1.5420242 Jun 1 584.8033 -2.3186712 -0.8033140 Jul 1 575.3506 -6.9540753 -2.3505988 Aug 1 566.4891 -8.1967650 0.5108532 Sep 1 565.5347 -3.4768164 3.4652831 Oct 1 604.6916 24.3028962 16.3084276 Nov 1 632.6436 26.6800781 -3.6436110 Dec 1 636.9463 12.1041838 -8.9463475 Jan 2 621.0974 -6.0368129 -9.0974341 Feb 2 603.6069 -13.5044740 -8.6068895 Mar 2 593.0170 -11.6530326 3.9829949 Apr 2 589.1569 -6.7801455 3.8431083 May 2 586.6714 -4.0597327 3.3285535 Jun 2 577.9869 -6.9731107 2.0130799 Jul 2 571.3757 -6.7463964 2.6243246 Aug 2 572.1989 -1.9965182 0.8011423 Sep 2 582.4876 5.7243033 -9.4875947 Oct 2 602.5313 14.7253286 17.4686532 Nov 2 621.8787 17.6290794 4.1213462 Dec 2 625.0184 8.5371646 -5.0184124 Jan 3 603.7374 -10.1792332 -15.7374121 Feb 3 579.4528 -19.0451486 -13.4527716 Mar 3 556.8242 -21.2798320 0.1758315 Apr 3 551.3686 -11.4706763 9.6313879 May 3 542.3934 -9.9162278 6.6065742 Jun 3 529.2883 -11.9039195 2.7117456 Jul 3 522.6186 -8.6533177 3.3813872 Aug 3 515.5568 -7.6656034 -4.5567652 Sep 3 514.1037 -3.8054357 -15.1036921 Oct 3 533.0719 10.3514803 21.9280946 Nov 3 553.5187 16.6247066 11.4813479 Dec 3 545.5982 1.3772243 -3.5981712 Jan 4 538.1219 -4.1269943 -11.1218876 Feb 4 525.8319 -9.2008946 -15.8319155 Mar 4 518.7236 -7.9054367 -4.7235810 Apr 4 508.6346 -9.2535850 8.3654112 May 4 499.0489 -9.4590583 8.9510593 Jun 4 490.1032 -9.1410088 2.8968112 Jul 4 484.4349 -6.9930131 5.5650691 Aug 4 477.5446 -6.9295963 -8.5445758 Sep 4 494.6621 7.9291414 -16.6621298 Oct 4 509.3104 12.0824502 18.6895737 Nov 4 515.8611 8.6630798 18.1389140 Dec 4 519.5459 5.5851741 -1.5458686 Jan 5 517.3087 0.7461167 -11.3086781 Feb 5 518.1632 0.8131116 -16.1631604 Mar 5 519.5609 1.1736617 -3.5609226 Apr 5 519.3099 0.2963425 8.6901116 May 5 521.8910 1.7050595 11.1089972 Jun 5 530.2626 5.8202861 5.7374112 Jul 5 533.7440 4.3776770 3.2559811 Aug 5 541.2105 6.2804838 -17.2105424 Sep 5 553.6194 10.0542088 -17.6193504 Oct 5 566.0709 11.5306555 20.9291312 Nov 5 577.0476 11.1894186 19.9523958 Dec 5 581.9699 7.3268804 -0.9698893 Jan 6 578.8903 0.9111039 -14.8902990 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.803372008 -0.006627903 0.726000184 0.511064435 -0.081564256 2 -2.405711700 -1.002686750 0.242380254 0.652622414 0.362119643 -0.383179201 3 -2.480238721 -1.173161841 -0.294269187 1.300115716 0.206362143 -0.262304976 4 -0.728236535 -0.670173352 0.170862477 -0.178331663 -0.027219996 0.042022746 5 -0.639668806 0.008845378 0.047583244 -0.115981796 0.186430479 0.543911553 6 -0.847732020 Jul Aug Sep Oct Nov Dec 1 -0.613579480 -0.164174374 0.623424007 3.669595218 0.314007345 -1.925349279 2 0.029927935 0.627993866 1.019877518 1.188959415 0.383580524 -1.201599624 3 0.428682138 0.130480696 0.509989560 1.870035092 0.828747655 -2.015748397 4 0.283393626 0.008372785 1.962849998 0.548708236 -0.451836215 -0.406891413 5 -0.190418594 0.251181020 0.498433498 0.195082947 -0.045100159 -0.510572478 6 > 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/1qssk1259953789.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/2e0y71259953789.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/37nqw1259953790.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/4neqc1259953790.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/5buv41259953790.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/65bhu1259953790.tab") > system("convert tmp/1qssk1259953789.ps tmp/1qssk1259953789.png") > system("convert tmp/2e0y71259953789.ps tmp/2e0y71259953789.png") > system("convert tmp/37nqw1259953790.ps tmp/37nqw1259953790.png") > system("convert tmp/4neqc1259953790.ps tmp/4neqc1259953790.png") > system("convert tmp/5buv41259953790.ps tmp/5buv41259953790.png") > > > proc.time() user system elapsed 1.266 0.814 1.462