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(7.3,7.6,7.5,7.6,7.9,7.9,8.1,8.2,8,7.5,6.8,6.5,6.6,7.6,8,8.1,7.7,7.5,7.6,7.8,7.8,7.8,7.5,7.5,7.1,7.5,7.5,7.6,7.7,7.7,7.9,8.1,8.2,8.2,8.2,7.9,7.3,6.9,6.6,6.7,6.9,7,7.1,7.2,7.1,6.9,7,6.8,6.4,6.7,6.6,6.4,6.3,6.2,6.5,6.8,6.8,6.4,6.1,5.8,6.1,7.2,7.3,6.9,6.1,5.8,6.2,7.1,7.7,7.9,7.7,7.4,7.5) > 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.008854744 0.090947153 0.000000000 0.000000000 > m$fitted level slope sea Jan 1 7.300000 0.0000000000 0.000000000 Feb 1 7.596395 0.2342854223 0.003604598 Mar 1 7.493578 -0.0750679612 0.006421723 Apr 1 7.593699 0.0856459617 0.006300757 May 1 7.893711 0.2824141660 0.006288605 Jun 1 7.893710 0.0231696542 0.006289919 Jul 1 8.093710 0.1854925901 0.006289851 Aug 1 8.193710 0.1070139322 0.006289854 Sep 1 7.993710 -0.1748121458 0.006289855 Oct 1 7.493710 -0.4733211317 0.006289855 Nov 1 6.793710 -0.6814029472 0.006289855 Dec 1 6.493710 -0.3312908335 0.006289855 Jan 2 6.622494 0.0871420214 -0.022494025 Feb 2 7.587536 0.7972277514 0.012463893 Mar 2 7.985591 0.4309751285 0.014409441 Apr 2 8.085458 0.1271446158 0.014542391 May 2 7.685440 -0.3567524123 0.014559763 Jun 2 7.485441 -0.2128601879 0.014559339 Jul 2 7.585441 0.0743325107 0.014559270 Aug 2 7.785441 0.1896900627 0.014559268 Sep 2 7.785441 0.0155624392 0.014559268 Oct 2 7.785441 0.0012767644 0.014559268 Nov 2 7.485441 -0.2752828311 0.014559268 Dec 2 7.485441 -0.0225845901 0.014559268 Jan 3 7.272263 -0.1965301673 -0.172263359 Feb 3 7.480906 0.1438345236 0.019093848 Mar 3 7.480394 0.0113757243 0.019605746 Apr 3 7.580420 0.0927309007 0.019579876 May 3 7.680420 0.0994036339 0.019579702 Jun 3 7.680420 0.0081552137 0.019579897 Jul 3 7.880420 0.1842607843 0.019579866 Aug 3 8.080420 0.1987087326 0.019579866 Sep 3 8.180420 0.1080982031 0.019579866 Oct 3 8.180420 0.0088685284 0.019579866 Nov 3 8.180420 0.0007275865 0.019579866 Dec 3 7.880420 -0.2753278864 0.019579866 Jan 4 7.519345 -0.3537097158 -0.219344589 Feb 4 6.881929 -0.5969513647 0.018070855 Mar 4 6.582759 -0.3236737728 0.017240880 Apr 4 6.682856 0.0652478194 0.017143749 May 4 6.882859 0.1889447560 0.017141214 Jun 4 6.982859 0.1072971527 0.017141351 Jul 4 7.082859 0.1005986687 0.017141352 Aug 4 7.182859 0.1000491156 0.017141352 Sep 4 7.082859 -0.0835876896 0.017141352 Oct 4 6.882859 -0.1904493705 0.017141352 Nov 4 6.982859 0.0761711257 0.017141352 Dec 4 6.782859 -0.1773425329 0.017141352 Jan 5 6.589150 -0.1923163587 -0.189149927 Feb 5 6.681271 0.0546877203 0.018728858 Mar 5 6.580915 -0.0876044754 0.019084788 Apr 5 6.380894 -0.1907803572 0.019106003 May 5 6.280895 -0.1074477402 0.019104597 Jun 5 6.180895 -0.1006110231 0.019104588 Jul 5 6.480895 0.2671333090 0.019104546 Aug 5 6.780895 0.2973035705 0.019104546 Sep 5 6.780895 0.0243912025 0.019104546 Oct 5 6.380895 -0.3651823496 0.019104546 Nov 5 6.080895 -0.3053476515 0.019104546 Dec 5 5.780895 -0.3004387288 0.019104546 Jan 6 6.286023 0.4370238505 -0.186022512 Feb 6 7.181000 0.8382020034 0.018999978 Mar 6 7.279557 0.1593656406 0.020443462 Apr 6 6.879467 -0.3541149639 0.020533186 May 6 6.079461 -0.7634189981 0.020539054 Jun 6 5.779461 -0.3380195453 0.020538554 Jul 6 6.179462 0.3394518399 0.020538488 Aug 6 7.079462 0.8540118416 0.020538484 Sep 6 7.679462 0.6208394883 0.020538485 Oct 6 7.879462 0.2345262627 0.020538485 Nov 6 7.679462 -0.1643508551 0.020538485 Dec 6 7.379462 -0.2888711536 0.020538485 Jan 7 7.710147 0.2785397814 -0.210147142 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.872311949 -0.970173911 0.532724147 0.652468147 -0.859636822 2 1.462968388 2.310335514 -1.178795029 -1.007283046 -1.604566891 0.477136641 3 -0.596561262 1.117024918 -0.429865863 0.269729623 0.022126299 -0.302573436 4 -0.266412201 -0.801098917 0.891022042 1.289492485 0.410170160 -0.270737793 5 -0.050630485 0.814982755 -0.465341983 -0.342092500 0.276325012 0.022670080 6 2.485110900 1.325130668 -2.224677896 -1.702532562 -1.357223097 1.410595102 7 1.907508766 Jul Aug Sep Oct Nov Dec 1 0.538251607 -0.260229791 -0.934515745 -0.989835112 -0.689984881 1.160947507 2 0.952311087 0.382517648 -0.577395133 -0.047370308 -0.917052454 0.837929856 3 0.583953869 0.047908396 -0.300458237 -0.329038725 -0.026994799 -0.915380818 4 -0.022211709 -0.001822280 -0.608926921 -0.354345929 0.884095089 -0.840633724 5 1.219414721 0.100042497 -0.904958500 -1.291798903 0.198407713 0.016277647 6 2.246448167 1.706245308 -0.773183365 -1.280987886 -1.322648882 -0.412900681 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/1ktuf1259954782.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/2eepu1259954782.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/3ls7k1259954782.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/4cx3q1259954782.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/5xc6t1259954782.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/6s0ar1259954782.tab") > system("convert tmp/1ktuf1259954782.ps tmp/1ktuf1259954782.png") > system("convert tmp/2eepu1259954782.ps tmp/2eepu1259954782.png") > system("convert tmp/3ls7k1259954782.ps tmp/3ls7k1259954782.png") > system("convert tmp/4cx3q1259954782.ps tmp/4cx3q1259954782.png") > system("convert tmp/5xc6t1259954782.ps tmp/5xc6t1259954782.png") > > > proc.time() user system elapsed 1.409 0.793 1.619