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(283.042,276.687,277.915,277.128,277.103,275.037,270.150,267.140,264.993,287.259,291.186,292.300,288.186,281.477,282.656,280.190,280.408,276.836,275.216,274.352,271.311,289.802,290.726,292.300,278.506,269.826,265.861,269.034,264.176,255.198,253.353,246.057,235.372,258.556,260.993,254.663,250.643,243.422,247.105,248.541,245.039,237.080,237.085,225.554,226.839,247.934,248.333,246.969,245.098,246.263,255.765,264.319,268.347,273.046,273.963,267.430,271.993,292.710,295.881,293.299) > 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.6407728 7.4292011 4.8275843 0.0000000 > m$fitted level slope sea Jan 1 283.0420 0.0000000 0.00000000 Feb 1 279.6431 -2.9236168 -2.95609995 Mar 1 276.1406 -3.2617563 1.77444054 Apr 1 275.5598 -1.6553316 1.56824057 May 1 276.3389 -0.1800070 0.76411626 Jun 1 275.6012 -0.5114944 -0.56417322 Jul 1 271.6036 -2.5837714 -1.45358071 Aug 1 267.3033 -3.6080429 -0.16332093 Sep 1 264.3292 -3.2294219 0.66381602 Oct 1 278.8965 7.3984700 8.36249504 Nov 1 291.4611 10.4826558 -0.27510007 Dec 1 295.8636 6.8534805 -3.56362144 Jan 2 292.2557 0.6598808 -4.06974641 Feb 2 286.2519 -3.2940322 -4.77485520 Mar 2 281.8927 -3.9072077 0.76325478 Apr 2 279.1774 -3.2264196 1.01260871 May 2 278.1557 -1.9489068 2.25227522 Jun 2 275.4431 -2.3901821 1.39285281 Jul 2 273.8508 -1.9318786 1.36524037 Aug 2 274.1186 -0.6669680 0.23340611 Sep 2 277.0051 1.3809298 -5.69412491 Oct 2 282.6168 3.8203242 7.18521283 Nov 2 288.3960 4.9487235 2.33002919 Dec 2 292.3533 4.3783861 -0.05330597 Jan 3 285.3165 -2.1963309 -6.81053632 Feb 3 276.2565 -6.1477143 -6.43046758 Mar 3 266.6346 -8.1287944 -0.77364159 Apr 3 265.0870 -4.3924439 3.94702174 May 3 261.3843 -3.9986414 2.79168470 Jun 3 254.6891 -5.5408491 0.50888835 Jul 3 251.2377 -4.3497383 2.11533110 Aug 3 247.5679 -3.9625895 -1.51092363 Sep 3 243.7303 -3.8913906 -8.35826903 Oct 3 248.9557 1.3071531 9.60031263 Nov 3 256.2906 4.7432096 4.70242037 Dec 3 253.6671 0.5434434 0.99591452 Jan 4 253.6415 0.2187482 -2.99854519 Feb 4 250.6214 -1.6273064 -7.19936820 Mar 4 249.9836 -1.0659285 -2.87856062 Apr 4 246.0038 -2.7154338 2.53724195 May 4 241.4245 -3.7732408 3.61445043 Jun 4 236.9397 -4.1778680 0.14026226 Jul 4 233.9596 -3.4977602 3.12543024 Aug 4 228.6947 -4.4994757 -3.14066622 Sep 4 234.7026 1.4547501 -7.86356024 Oct 4 239.9354 3.5963144 7.99859013 Nov 4 241.8379 2.6358953 6.49506625 Dec 4 245.1634 3.0271715 1.80557409 Jan 5 247.5747 2.6775554 -2.47672781 Feb 5 252.6887 4.0588306 -6.42565313 Mar 5 257.3126 4.3783318 -1.54763391 Apr 5 261.0294 4.0047844 3.28960926 May 5 264.2460 3.5590730 4.10100068 Jun 5 271.0337 5.3876454 2.01231560 Jul 5 272.7624 3.3166038 1.20058417 Aug 5 275.1943 2.8166559 -7.76432278 Sep 5 280.1861 4.0446335 -8.19308340 Oct 5 284.2118 4.0339496 8.49820722 Nov 5 289.0686 4.4987710 6.81243857 Dec 5 292.1353 3.6891466 1.16366322 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.602608670 -0.125865217 0.625674701 0.536882430 -0.120727949 2 -2.311042496 -1.462259903 -0.222981967 0.253176943 0.472296550 -0.161077837 3 -2.423435738 -1.450033762 -0.724380149 1.376118853 0.145243434 -0.565358818 4 -0.119323540 -0.676709885 0.205616875 -0.606230501 -0.389366177 -0.148538252 5 -0.128339537 0.506263793 0.117099127 -0.137195331 -0.163895921 0.671488107 Jul Aug Sep Oct Nov Dec 1 -0.761464921 -0.376027614 0.138895033 3.899205711 1.131560095 -1.331485946 2 0.167747257 0.464367162 0.751506379 0.894942062 0.414038743 -0.209592665 3 0.435929106 0.141949478 0.026126051 1.907609365 1.261340796 -1.543293647 4 0.249143943 -0.367107800 2.184243871 0.786066647 -0.352698905 0.143736204 5 -0.759182632 -0.183216475 0.450382121 -0.003922001 0.170729131 -0.297349472 > 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/1riqv1260471343.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/27v431260471343.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/3vs1i1260471343.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/4lzrp1260471343.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/5vpav1260471343.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/6x2ot1260471343.tab") > system("convert tmp/1riqv1260471343.ps tmp/1riqv1260471343.png") > system("convert tmp/27v431260471343.ps tmp/27v431260471343.png") > system("convert tmp/3vs1i1260471343.ps tmp/3vs1i1260471343.png") > system("convert tmp/4lzrp1260471343.ps tmp/4lzrp1260471343.png") > system("convert tmp/5vpav1260471343.ps tmp/5vpav1260471343.png") > > > proc.time() user system elapsed 1.308 0.826 1.545