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(8.3,8.2,8,7.9,7.6,7.6,8.3,8.4,8.4,8.4,8.4,8.6,8.9,8.8,8.3,7.5,7.2,7.4,8.8,9.3,9.3,8.7,8.2,8.3,8.5,8.6,8.5,8.2,8.1,7.9,8.6,8.7,8.7,8.5,8.4,8.5,8.7,8.7,8.6,8.5,8.3,8,8.2,8.1,8.1,8,7.9,7.9,8,8,7.9,8,7.7,7.2,7.5,7.3,7,7,7,7.2,7.3,7.1,6.8,6.4,6.1,6.5,7.7,7.9,7.5,6.9,6.6,6.9) > 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.000000000 0.103322238 0.001056948 0.000000000 > m$fitted level slope sea Jan 1 8.300000 0.000000000 0.0000000000 Feb 1 8.201508 -0.098616680 -0.0015079432 Mar 1 8.002831 -0.194512079 -0.0028305775 Apr 1 7.896198 -0.110992222 0.0038018623 May 1 7.607103 -0.280349045 -0.0071029031 Jun 1 7.589435 -0.030515472 0.0105647217 Jul 1 8.277564 0.652962935 0.0224364997 Aug 1 8.424115 0.171335733 -0.0241146509 Sep 1 8.402298 -0.012363897 -0.0022979408 Oct 1 8.398119 -0.004579820 0.0018813300 Nov 1 8.399980 0.001545771 0.0000203504 Dec 1 8.593530 0.184154019 0.0064697472 Jan 2 8.897237 0.297819722 0.0027627888 Feb 2 8.812549 -0.066156765 -0.0125485789 Mar 2 8.332106 -0.449554429 -0.0321055174 Apr 2 7.496548 -0.806802769 0.0034524896 May 2 7.168255 -0.364379269 0.0317453182 Jun 2 7.411106 0.197272316 -0.0111063707 Jul 2 8.707787 1.214171158 0.0922131452 Aug 2 9.358290 0.692818917 -0.0582895730 Sep 2 9.337246 0.032550385 -0.0372457914 Oct 2 8.718675 -0.569688908 -0.0186751705 Nov 2 8.193788 -0.528250175 0.0062121004 Dec 2 8.275650 0.036044886 0.0243498718 Jan 3 8.486811 0.197960410 0.0131892538 Feb 3 8.606593 0.125609684 -0.0065932213 Mar 3 8.502791 -0.083634686 -0.0027906674 Apr 3 8.215310 -0.269995115 -0.0153098034 May 3 8.060847 -0.164543540 0.0391527393 Jun 3 8.006790 -0.063657951 -0.1067901578 Jul 3 8.464589 0.412516488 0.1354111138 Aug 3 8.744668 0.291582226 -0.0446676507 Sep 3 8.706978 -0.009084762 -0.0069781571 Oct 3 8.507064 -0.183337773 -0.0070640464 Nov 3 8.423496 -0.092233652 -0.0234955133 Dec 3 8.482545 0.045897687 0.0174547115 Jan 4 8.682091 0.186159733 0.0179094237 Feb 4 8.706676 0.038598812 -0.0066763381 Mar 4 8.591384 -0.100862299 0.0086162378 Apr 4 8.513852 -0.079671614 -0.0138520800 May 4 8.258906 -0.238637913 0.0410944608 Jun 4 8.149935 -0.121006421 -0.1499346134 Jul 4 8.065827 -0.087528554 0.1341728946 Aug 4 8.116175 0.037563817 -0.0161754099 Sep 4 8.090885 -0.019461246 0.0091149129 Oct 4 8.011865 -0.073497503 -0.0118646141 Nov 4 7.931211 -0.079990177 -0.0312108580 Dec 4 7.899634 -0.036071899 0.0003655977 Jan 5 7.966765 0.057544157 0.0332352742 Feb 5 7.996013 0.031877797 0.0039870377 Mar 5 7.907360 -0.076966587 -0.0073603247 Apr 5 7.983136 0.061209415 0.0168640391 May 5 7.689112 -0.259801448 0.0108880869 Jun 5 7.352600 -0.329126626 -0.1525998764 Jul 5 7.360126 -0.024839462 0.1398744944 Aug 5 7.314149 -0.043944149 -0.0141494167 Sep 5 7.002708 -0.285723115 -0.0027076340 Oct 5 6.991547 -0.037556793 0.0084528181 Nov 5 7.026933 0.028371934 -0.0269327098 Dec 5 7.202391 0.161294441 -0.0023913775 Jan 6 7.272887 0.079231867 0.0271131680 Feb 6 7.099853 -0.148686133 0.0001468148 Mar 6 6.842916 -0.246218775 -0.0429155977 Apr 6 6.360326 -0.459465532 0.0396740059 May 6 6.058975 -0.316922619 0.0410254733 Jun 6 6.635227 0.488227690 -0.1352270806 Jul 6 7.540674 0.864392093 0.1593259234 Aug 6 7.910566 0.418541278 -0.0105657417 Sep 6 7.573206 -0.262988182 -0.0732063272 Oct 6 6.916618 -0.617863481 -0.0166177472 Nov 6 6.623137 -0.325406663 -0.0231366116 Dec 6 6.867074 0.187822974 0.0329257574 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.30904054 -0.30437703 0.25949091 -0.52692225 0.77723668 2 0.35378064 -1.13787364 -1.20575177 -1.10993648 1.37659021 1.74730949 3 0.50401449 -0.22551529 -0.65424136 -0.57947818 0.32805402 0.31386634 4 0.43666040 -0.45928969 -0.43492708 0.06592860 -0.49445274 0.36597435 5 0.29143715 -0.07984661 -0.33906141 0.43000321 -0.99840201 -0.21568451 6 -0.25544483 -0.70892622 -0.30366264 -0.66367608 0.44333000 2.50494243 Jul Aug Sep Oct Nov Dec 1 2.12631693 -1.49835287 -0.57149361 0.02421644 0.01905685 0.56809830 2 3.16359847 -1.62193836 -2.05411003 -1.87357980 0.12891690 1.75553886 3 1.48138603 -0.37622942 -0.93538167 -0.54210501 0.28342794 0.42973244 4 0.10414984 0.38916545 -0.17740623 -0.16810809 -0.02019894 0.13663412 5 0.94664144 -0.05943506 -0.75217988 0.77205256 0.20510600 0.41354806 6 1.17026054 -1.38704799 -2.12025519 -1.10402931 0.90983740 1.59680928 > 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/100s21260547397.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/2crip1260547397.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/3vwg01260547397.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/4km8j1260547397.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/50esj1260547397.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/6dxdk1260547397.tab") > system("convert tmp/100s21260547397.ps tmp/100s21260547397.png") > system("convert tmp/2crip1260547397.ps tmp/2crip1260547397.png") > system("convert tmp/3vwg01260547397.ps tmp/3vwg01260547397.png") > system("convert tmp/4km8j1260547397.ps tmp/4km8j1260547397.png") > system("convert tmp/50esj1260547397.ps tmp/50esj1260547397.png") > > > proc.time() user system elapsed 1.458 0.816 2.717