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(12.610,10.862,52.929,56.902,81.776,87.876,82.103,72.846,60.632,33.521,15.342,7.758,8.668,13.082,38.157,58.263,81.153,88.476,72.329,75.845,61.108,37.665,12.755,2.793,12.935,19.533,33.404,52.074,70.735,69.702,61.656,82.993,53.990,32.283,15.686,2.713,12.842,19.244,48.488,54.464,84.192,84.458,85.793,75.163,68.212,49.233,24.302,5.402,15.058,33.559,70.358,85.934,94.452,129.305,113.882,107.256,94.274,57.842,26.611,14.521) > 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.0000000 164.4306600 0.9378151 21.7473735 > m$fitted level slope sea Jan 1 12.610000 0.0000000 0.00000000 Feb 1 11.082867 -1.5086032 -0.02256684 Mar 1 48.996294 31.2540590 0.49556542 Apr 1 59.177071 13.5654350 -0.47005561 May 1 80.930887 20.4465498 0.14866385 Jun 1 89.150610 10.1692318 -0.23455173 Jul 1 83.643185 -3.0079780 -0.20668549 Aug 1 73.529083 -8.9809535 -0.07863657 Sep 1 60.996312 -11.9663703 -0.06219743 Oct 1 34.936946 -23.8119868 -0.21720588 Nov 1 14.953496 -20.5939759 0.06285099 Dec 1 6.605839 -10.3005421 0.11049722 Jan 2 5.817539 -2.4032949 2.05128338 Feb 2 13.192733 4.9144924 -0.72958515 Mar 2 35.435929 19.4053781 1.29951169 Apr 2 58.708928 22.6174030 -0.76207766 May 2 81.006892 22.3517924 0.17249188 Jun 2 89.961560 11.2032456 -0.37780460 Jul 2 75.193759 -10.4089161 -0.71730093 Aug 2 74.766705 -2.1028622 0.25297158 Sep 2 62.127772 -10.8699443 -0.14863809 Oct 2 39.128579 -20.9628440 -0.46070668 Nov 2 13.307075 -25.0058677 -0.15034284 Dec 2 1.382822 -14.1239582 0.32892355 Jan 3 7.989202 3.0286793 3.23726767 Feb 3 20.046037 10.0015744 -1.13510616 Mar 3 32.156482 11.7469119 1.07716399 Apr 3 52.231202 18.6136292 -0.82476760 May 3 70.630589 18.4368745 0.12172715 Jun 3 71.255292 3.7294212 -0.11180767 Jul 3 64.048771 -5.3006278 -1.50773575 Aug 3 80.152069 12.3719912 1.10881922 Sep 3 58.019859 -16.1167299 -1.23764957 Oct 3 33.209494 -23.2949206 -0.22295191 Nov 3 15.347350 -18.8091537 -0.10100514 Dec 3 2.388042 -13.9821879 -0.14810992 Jan 4 7.433685 1.6961147 3.86547678 Feb 4 19.366648 9.7330808 -0.85829455 Mar 4 45.975029 23.5956491 1.16731101 Apr 4 57.221713 13.4655843 -1.78106470 May 4 81.962227 22.7148558 1.33231820 Jun 4 86.149216 7.5027559 -0.21429381 Jul 4 89.217309 3.8614624 -3.07078582 Aug 4 74.861272 -11.0959683 1.75392238 Sep 4 68.324682 -7.3525548 -0.47612508 Oct 4 50.786180 -15.7159179 -0.74119193 Nov 4 25.387558 -23.6657058 -0.31372495 Dec 4 5.614009 -20.4728881 -0.52197067 Jan 5 8.616414 -1.2026789 4.56234400 Feb 5 32.794989 18.9140962 -1.10358251 Mar 5 66.765661 31.2178664 2.40348466 Apr 5 89.696537 24.4442224 -3.11283650 May 5 94.603209 8.4812608 1.38772781 Jun 5 126.500437 27.6296193 0.95704555 Jul 5 120.186806 -0.1308739 -3.62642714 Aug 5 107.242543 -10.6097214 1.02449180 Sep 5 94.504069 -12.3506039 -0.06210227 Oct 5 60.444621 -30.1046834 -0.88963964 Nov 5 26.729104 -33.0573603 0.16677785 Dec 5 13.317364 -17.0074905 -0.34485727 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.12624935 2.61768315 -1.38565519 0.53664081 -0.80148497 2 0.65593947 0.55055322 1.11850227 0.24912712 -0.02071130 -0.86940760 3 1.38744459 0.53230507 0.13559497 0.53291949 -0.01378198 -1.14694751 4 1.25103834 0.61802439 1.08064952 -0.78657303 0.72112908 -1.18630836 5 1.52525780 1.55360080 0.96078217 -0.52622661 -1.24444063 1.49328823 Jul Aug Sep Oct Nov Dec 1 -1.02762168 -0.46579984 -0.23281640 -0.92377514 0.25095514 0.80272886 2 -1.68540782 0.64774389 -0.68369699 -0.78709038 -0.31529372 0.84865588 3 -0.70420021 1.37819121 -2.22168026 -0.55978850 0.34982139 0.37650647 4 -0.28396217 -1.16644826 0.29192847 -0.65221513 -0.61996062 0.24910681 5 -2.16486342 -0.81718759 -0.13576202 -1.38455256 -0.23026187 1.25257362 > 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/14t4g1259922893.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/2ht021259922893.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/37lqg1259922893.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/4ygft1259922893.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/5oae51259922893.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/6h4ta1259922893.tab") > system("convert tmp/14t4g1259922893.ps tmp/14t4g1259922893.png") > system("convert tmp/2ht021259922893.ps tmp/2ht021259922893.png") > system("convert tmp/37lqg1259922893.ps tmp/37lqg1259922893.png") > system("convert tmp/4ygft1259922893.ps tmp/4ygft1259922893.png") > system("convert tmp/5oae51259922893.ps tmp/5oae51259922893.png") > > > proc.time() user system elapsed 1.444 0.789 1.861