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(102.86,102.55,102.28,102.26,102.57,103.08,102.76,102.51,102.87,103.14,103.12,103.16,102.48,102.57,102.88,102.63,102.38,101.69,101.96,102.19,101.87,101.6,101.63,101.22,101.21,101.49,101.64,101.66,101.77,101.82,101.78,101.28,101.29,101.37,101.12,101.51,102.24,102.94,103.09,103.46,103.64,104.39,104.15,105.21,105.8,105.91,105.39,105.46,104.72,103.14,102.63,102.32,101.93,100.62,100.6,99.63,98.9,98.32,99.22,98.81) > 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.17299278 0.01490638 0.00000000 0.00000000 > m$fitted level slope sea Jan 1 102.86000 0.000000000 0.0000000000 Feb 1 102.56568 -0.024101798 -0.0156821794 Mar 1 102.29565 -0.059019789 -0.0156542556 Apr 1 102.27566 -0.051770688 -0.0156578635 May 1 102.58568 0.025577406 -0.0156841616 Jun 1 103.09571 0.137359171 -0.0157112499 Jul 1 102.77569 0.027294510 -0.0156918296 Aug 1 102.52568 -0.041008229 -0.0156829555 Sep 1 102.88569 0.059052589 -0.0156925866 Oct 1 103.15570 0.112068677 -0.0156963797 Nov 1 103.13569 0.078743572 -0.0156946041 Dec 1 103.17569 0.068945531 -0.0156942150 Jan 2 102.57586 -0.094256794 -0.0958649243 Feb 2 102.55954 -0.075416433 0.0104587084 Mar 2 102.86997 0.022529692 0.0100253447 Apr 2 102.61975 -0.046664898 0.0102539759 May 2 102.36962 -0.098264758 0.0103812700 Jun 2 101.67934 -0.248385199 0.0106577344 Jul 2 101.94952 -0.116894197 0.0104769737 Aug 2 102.17961 -0.028910530 0.0103866917 Sep 2 101.85956 -0.102736662 0.0104432362 Oct 2 101.58953 -0.145156838 0.0104674872 Nov 2 101.61955 -0.100735447 0.0104485323 Dec 2 101.20953 -0.179167124 0.0104735123 Jan 3 101.27836 -0.117221366 -0.0683597098 Feb 3 101.47921 -0.039032483 0.0107883265 Mar 3 101.62935 0.008991094 0.0106465196 Apr 3 101.64936 0.011785729 0.0106403574 May 3 101.75940 0.036706853 0.0105993321 Jun 3 101.80940 0.040079075 0.0105951880 Jul 3 101.76939 0.019767243 0.0106138202 Aug 3 101.26930 -0.112060396 0.0107040833 Sep 3 101.27931 -0.081103725 0.0106882622 Oct 3 101.35933 -0.040245926 0.0106726762 Nov 3 101.10931 -0.093441338 0.0106878224 Dec 3 101.49934 0.029162601 0.0106617666 Jan 4 102.26620 0.214411009 -0.0262032251 Feb 4 102.93282 0.326494682 0.0071768503 Mar 4 103.08273 0.281676616 0.0072749420 Apr 4 103.45276 0.304092045 0.0072383121 May 4 103.63272 0.272609057 0.0072767193 Jun 4 104.38283 0.393704707 0.0071664437 Jul 4 104.14272 0.232974300 0.0072756991 Aug 4 105.20283 0.442726256 0.0071692765 Sep 4 105.79284 0.480076971 0.0071551315 Oct 4 105.90282 0.386221807 0.0071816614 Nov 4 105.38277 0.156397239 0.0072301507 Dec 4 105.45277 0.134486376 0.0072336012 Jan 5 104.86573 -0.047162470 -0.1457256429 Feb 5 103.14084 -0.465004434 -0.0008411156 Mar 5 102.63082 -0.476427294 -0.0008212632 Apr 5 102.32088 -0.434196322 -0.0008760578 May 5 101.93089 -0.422984350 -0.0008869175 Jun 5 100.62072 -0.647975546 -0.0007242471 Jul 5 100.60081 -0.488702372 -0.0008102028 Aug 5 99.63076 -0.610768310 -0.0007610321 Sep 5 98.90075 -0.641006919 -0.0007519403 Oct 5 98.32076 -0.625535033 -0.0007554125 Nov 5 99.22082 -0.238649001 -0.0008202182 Dec 5 98.81081 -0.282104592 -0.0008147852 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.41844744 -0.54762662 0.08465289 0.77123147 1.02151189 2 -1.61460401 0.13558246 0.80029043 -0.56597742 -0.42231394 -1.22905563 3 0.55664424 0.59461745 0.39253420 0.02286350 0.20398818 0.02761061 4 1.61463477 0.87131756 -0.36652715 0.18343947 -0.25774172 0.99158087 5 -1.55911142 -3.28718271 -0.09344672 0.34566281 0.09179796 -1.84242012 Jul Aug Sep Oct Nov Dec 1 -0.95821655 -0.57879029 0.83523360 0.43884105 -0.27456194 -0.08051496 2 1.07673397 0.72054214 -0.60463389 -0.34743116 0.36382810 -0.64239169 3 -0.16633251 -1.07962398 0.25353707 0.33463702 -0.43569218 1.00418487 4 -1.31627990 1.71784674 0.30590954 -0.76870733 -1.88236647 -0.17946094 5 1.30438547 -0.99972417 -0.24766235 0.12672094 3.16877964 -0.35592362 > 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/1wsqi1260465045.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/2w0fx1260465045.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/39ov91260465045.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/4rlce1260465045.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/5lqsx1260465045.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/6y4sj1260465045.tab") > system("convert tmp/1wsqi1260465045.ps tmp/1wsqi1260465045.png") > system("convert tmp/2w0fx1260465045.ps tmp/2w0fx1260465045.png") > system("convert tmp/39ov91260465045.ps tmp/39ov91260465045.png") > system("convert tmp/4rlce1260465045.ps tmp/4rlce1260465045.png") > system("convert tmp/5lqsx1260465045.ps tmp/5lqsx1260465045.png") > > > proc.time() user system elapsed 1.271 0.796 1.526