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(14.2,13.5,11.9,14.6,15.6,14.1,14.9,14.2,14.6,17.2,15.4,14.3,17.5,14.5,14.4,16.6,16.7,16.6,16.9,15.7,16.4,18.4,16.9,16.5,18.3,15.1,15.7,18.1,16.8,18.9,19,18.1,17.8,21.5,17.1,18.7,19,16.4,16.9,18.6,19.3,19.4,17.6,18.6,18.1,20.4,18.1,19.6,19.9,19.2,17.8,19.2,22,21.1,19.5,22.2,20.9,22.2,23.5,21.5,24.3,22.8,20.3,23.7,23.3,19.6,18,17.3,16.8,18.2,16.5,16,18.4) > 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.4805987 0.0000000 0.3006777 0.0000000 > m$fitted level slope sea Jan 1 14.20000 0.0000000000 0.00000000 Feb 1 13.69629 -0.0218460906 -0.19629362 Mar 1 12.57026 -0.0799509748 -0.67026247 Apr 1 13.53704 -0.0461511237 1.06295613 May 1 14.84642 -0.0239254681 0.75357512 Jun 1 14.69334 -0.0252516531 -0.59334251 Jul 1 14.78468 -0.0241722368 0.11532413 Aug 1 14.48843 -0.0268099521 -0.28842854 Sep 1 14.50094 -0.0264128347 0.09906114 Oct 1 15.99069 -0.0109231755 1.20930660 Nov 1 15.95316 -0.0111935145 -0.55315688 Dec 1 15.04284 -0.0202241713 -0.74283517 Jan 2 15.84458 -0.0439215904 1.65542250 Feb 2 15.13869 -0.0440762184 -0.63868989 Mar 2 15.31251 -0.0395886555 -0.91250854 Apr 2 15.75959 -0.0285755811 0.84041127 May 2 15.88332 -0.0260199958 0.81667934 Jun 2 16.56633 -0.0181442951 0.03366949 Jul 2 16.71382 -0.0167963714 0.18618360 Aug 2 16.45806 -0.0185438448 -0.75805785 Sep 2 16.68822 -0.0166973841 -0.28822002 Oct 2 16.93194 -0.0148072291 1.46806324 Nov 2 17.07679 -0.0139431355 -0.17679358 Dec 2 17.34915 -0.0137442623 -0.84914570 Jan 3 16.91939 -0.0116906701 1.38061400 Feb 3 16.35344 -0.0125862935 -1.25343770 Mar 3 16.54187 -0.0103600919 -0.84187425 Apr 3 17.02122 -0.0029033436 1.07878008 May 3 16.71424 -0.0072301156 0.08575549 Jun 3 17.78385 0.0048765877 1.11615253 Jul 3 18.42086 0.0103573765 0.57913889 Aug 3 18.82179 0.0131587939 -0.72179043 Sep 3 18.59667 0.0116609696 -0.79666721 Oct 3 19.37172 0.0156081322 2.12827923 Nov 3 18.44814 0.0125830048 -1.34814321 Dec 3 18.74566 0.0128163869 -0.04565775 Jan 4 18.10602 0.0128228094 0.89397619 Feb 4 17.86389 0.0120311021 -1.46389160 Mar 4 17.87196 0.0119996193 -0.97196296 Apr 4 17.69897 0.0099353279 0.90103052 May 4 18.77535 0.0223781798 0.52464560 Jun 4 18.86005 0.0230238900 0.53994686 Jul 4 18.02975 0.0157359946 -0.42975079 Aug 4 18.58096 0.0194460942 0.01903989 Sep 4 18.95412 0.0214121214 -0.85411578 Oct 4 18.44427 0.0192026570 1.95572969 Nov 4 18.90655 0.0203767241 -0.80655249 Dec 4 19.17616 0.0207527011 0.42383995 Jan 5 19.07053 0.0205402530 0.82947195 Feb 5 19.89807 0.0235132400 -0.69807076 Mar 5 19.49015 0.0206687710 -1.69015068 Apr 5 19.08184 0.0168603130 0.11816335 May 5 20.22482 0.0277871044 1.77518306 Jun 5 20.39669 0.0291123026 0.70331144 Jul 5 20.33149 0.0283612771 -0.83148614 Aug 5 21.32434 0.0346525445 0.87566271 Sep 5 21.59882 0.0358768769 -0.69881553 Oct 5 21.08109 0.0337927449 1.11890673 Nov 5 22.74526 0.0380830958 0.75474387 Dec 5 22.16712 0.0368199501 -0.66711613 Jan 6 22.93006 0.0385877169 1.36993880 Feb 6 23.23088 0.0396008318 -0.43088318 Mar 6 22.63757 0.0359232148 -2.33756747 Apr 6 23.35172 0.0409833005 0.34828115 May 6 22.56563 0.0341608227 0.73437078 Jun 6 20.63886 0.0182954705 -1.03885845 Jul 6 19.71386 0.0114594679 -1.71385626 Aug 6 17.96302 0.0007780645 -0.66302264 Sep 6 17.48063 -0.0015389524 -0.68062597 Oct 6 17.58043 -0.0011706745 0.61956714 Nov 6 16.55663 -0.0040021564 -0.05662945 Dec 6 16.65275 -0.0037581770 -0.65275068 Jan 7 16.87693 -0.0031147890 1.52307297 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.51096130 -1.28905028 1.46329924 1.95715040 -0.18679128 2 1.28248842 -0.95650384 0.29004144 0.66998995 0.21717761 1.02223378 3 -0.60930543 -0.79657962 0.27951867 0.68202789 -0.43139290 1.54656751 4 -0.94318038 -0.36490987 -0.00558289 -0.26018779 1.51362961 0.08930831 5 -0.18193801 1.15412229 -0.61195918 -0.60712352 1.60147036 0.20631930 6 1.04361331 0.37514617 -0.90085460 0.96376954 -1.17859558 -2.80788066 7 0.32738422 Jul Aug Sep Oct Nov Dec 1 0.16845068 -0.39277340 0.05673822 2.18754213 -0.03839854 -1.29732710 2 0.23920688 -0.34506955 0.35904019 0.37583875 0.23023539 0.41346328 3 0.91199350 0.56410683 -0.34411642 1.10155995 -1.35410045 0.41132085 4 -1.22942482 0.77294592 0.51059397 -0.76628887 0.63871706 0.35931565 5 -0.13571590 1.39088119 0.34594720 -0.79813563 2.34932923 -0.88747680 6 -1.35616079 -2.53892633 -0.69644836 0.14603766 -1.47289973 0.14409105 7 > 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/1j98d1259922178.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/2djn81259922178.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/3uxum1259922178.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/4dija1259922178.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/5zuv81259922178.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/6vmf41259922178.tab") > system("convert tmp/1j98d1259922178.ps tmp/1j98d1259922178.png") > system("convert tmp/2djn81259922178.ps tmp/2djn81259922178.png") > system("convert tmp/3uxum1259922178.ps tmp/3uxum1259922178.png") > system("convert tmp/4dija1259922178.ps tmp/4dija1259922178.png") > system("convert tmp/5zuv81259922178.ps tmp/5zuv81259922178.png") > > > proc.time() user system elapsed 1.525 0.798 1.764