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.5,14.3,15.3,14.4,13.7,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.4517202 0.0000000 0.2605746 0.0000000 > m$fitted level slope sea Jan 1 14.50000 0.000000000 0.00000000 Feb 1 14.35322 -0.006502253 -0.05322182 Mar 1 14.90425 0.022151218 0.39575319 Apr 1 14.69348 0.015070672 -0.29348351 May 1 14.09893 0.005847680 -0.39893063 Jun 1 14.04481 0.005268923 0.15519059 Jul 1 13.72787 0.002387403 -0.22786845 Aug 1 12.64404 -0.007856503 -0.74403904 Sep 1 13.54436 0.001044583 1.05563778 Oct 1 14.86163 0.014055682 0.73836812 Nov 1 14.67178 0.012054124 -0.57178088 Dec 1 14.76127 0.012805650 0.13873098 Jan 2 14.55051 0.019392973 -0.35051497 Feb 2 14.73947 0.019564219 -0.13946722 Mar 2 15.85958 0.042461484 1.34041788 Apr 2 15.80206 0.040257908 -0.40206132 May 2 15.17575 0.029676502 -0.87574555 Jun 2 16.06779 0.038598358 1.43220762 Jul 2 15.26164 0.032152180 -0.76164283 Aug 2 15.33867 0.032465542 -0.93866591 Sep 2 15.74168 0.035128778 0.85832216 Oct 2 15.86474 0.035751657 0.83526468 Nov 2 16.57496 0.039314251 0.02504433 Dec 2 16.76615 0.039390137 0.13384916 Jan 3 16.58984 0.040541306 -0.88984059 Feb 3 16.83798 0.040887434 -0.43798003 Mar 3 16.94751 0.041658834 1.45248751 Apr 3 17.05232 0.042610858 -0.15231864 May 3 17.38969 0.046663012 -0.88968610 Jun 3 16.96071 0.041596933 1.33929488 Jul 3 16.36085 0.036363282 -1.26084999 Aug 3 16.52620 0.037241851 -0.82619590 Sep 3 16.99732 0.039873277 1.10268063 Oct 3 16.69837 0.038165450 0.10162504 Nov 3 17.83478 0.041568879 1.06521695 Dec 3 18.52044 0.041963963 0.47955764 Jan 4 18.95198 0.041857156 -0.85197564 Feb 4 18.70066 0.040974342 -0.90066362 Mar 4 19.45126 0.046639682 2.04874492 Apr 4 18.43232 0.034811944 -1.33232371 May 4 18.74684 0.038003186 -0.04684422 Jun 4 18.07572 0.030937321 0.92428344 Jul 4 17.81855 0.028594228 -1.41855411 Aug 4 17.82604 0.028454894 -0.92604461 Sep 4 17.66626 0.027450452 0.93373829 Oct 4 18.78064 0.031802484 0.51936472 Nov 4 18.89864 0.032017559 0.50136062 Dec 4 18.08161 0.030906428 -0.48160563 Jan 5 18.67548 0.031736148 -0.07548078 Feb 5 19.03825 0.032916586 -0.93824595 Mar 5 18.47228 0.028986092 1.92772464 Apr 5 18.92022 0.032686126 -0.82022447 May 5 19.16895 0.034749711 0.43105283 Jun 5 19.02946 0.033192000 0.87053690 Jul 5 19.84717 0.039204455 -0.64717348 Aug 5 19.42081 0.036294936 -1.62081040 Sep 5 19.03122 0.034214180 0.16878044 Oct 5 20.23765 0.038414487 1.76235176 Nov 5 20.44013 0.038816714 0.65986776 Dec 5 20.39496 0.038660949 -0.89496473 Jan 6 21.43690 0.040921462 0.76309749 Feb 6 21.69853 0.041747690 -0.79852767 Mar 6 21.12048 0.038180879 1.07951603 Apr 6 22.80500 0.050387815 0.69500164 May 6 22.13852 0.044544237 -0.63852227 Jun 6 22.88133 0.050075777 1.41867300 Jul 6 23.15466 0.051645505 -0.35466438 Aug 6 22.54186 0.047766990 -2.24186080 Sep 6 23.29007 0.050980136 0.40992593 Oct 6 22.53778 0.048216778 0.76221841 Nov 6 20.62854 0.043172767 -1.02853731 Dec 6 19.74265 0.041086156 -1.74265127 Jan 7 17.98953 0.036322510 -0.68953051 Feb 7 17.52403 0.034432938 -0.72402710 Mar 7 17.62154 0.034761232 0.57846057 Apr 7 16.56229 0.027711769 -0.06228989 May 7 16.66140 0.028217764 -0.66139976 Jun 7 16.88681 0.029602969 1.51318687 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.15182930 0.68030420 -0.33763062 -0.90856653 -0.08946153 2 -0.36079768 0.25173233 1.51102679 -0.14248623 -0.98244575 1.28323239 3 -0.32619434 0.30761791 0.09839646 0.09084005 0.43210402 -0.70532432 4 0.58118643 -0.43283724 1.03175369 -1.54675675 0.40996893 -1.04921918 5 0.83621511 0.48837487 -0.87616673 0.61178597 0.31713577 -0.25755292 6 1.48771109 0.32570836 -0.90993732 2.41351502 -1.05422109 1.03198624 7 -2.65871347 -0.74102542 0.09282572 -1.60820758 0.10516604 0.29149651 Jul Aug Sep Oct Nov Dec 1 -0.48020798 -1.61756646 1.35193603 1.95920790 -0.30352215 0.11526595 2 -1.25856165 0.06683560 0.55178282 0.13090789 1.00331810 0.22626898 3 -0.95489801 0.19217376 0.64635146 -0.50433493 1.63337617 0.95923317 4 -0.42836751 -0.03142939 -0.28032183 1.61730217 0.12818692 -1.26286987 5 1.16511728 -0.69277353 -0.63381815 1.74347338 0.24389846 -0.12478757 6 0.33127315 -0.98781571 1.04170245 -1.19418334 -2.90863296 -1.37943339 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/1g3jb1259939155.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/2vzdf1259939155.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/3n2p81259939155.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/4iqok1259939155.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/5quc91259939155.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/6vsnj1259939155.tab") > system("convert tmp/1g3jb1259939155.ps tmp/1g3jb1259939155.png") > system("convert tmp/2vzdf1259939155.ps tmp/2vzdf1259939155.png") > system("convert tmp/3n2p81259939155.ps tmp/3n2p81259939155.png") > system("convert tmp/4iqok1259939155.ps tmp/4iqok1259939155.png") > system("convert tmp/5quc91259939155.ps tmp/5quc91259939155.png") > > > proc.time() user system elapsed 1.472 0.822 1.674