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(5.7,6.1,6,5.9,5.8,5.7,5.6,5.4,5.4,5.5,5.6,5.7,5.9,6.1,6,5.8,5.8,5.7,5.5,5.3,5.2,5.2,5,5.1,5.1,5.2,4.9,4.8,4.5,4.5,4.4,4.4,4.2,4.1,3.9,3.8,3.9,4.2,4.1,3.8,3.6,3.7,3.5,3.4,3.1,3.1,3.1,3.2,3.3,3.5,3.6,3.5,3.3,3.2,3.1,3.2,3,3,3.1,3.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 2.209261e-02 1.337643e-04 4.698744e-05 0.000000e+00 > m$fitted level slope sea Jan 1 5.700000 0.000000000 0.0000000000 Feb 1 6.078849 0.021506092 0.0211511918 Mar 1 5.980814 0.019621003 0.0191861075 Apr 1 5.880697 0.017123727 0.0193025559 May 1 5.781088 0.014052703 0.0189118483 Jun 1 5.681463 0.010486770 0.0185374290 Jul 1 5.581809 0.006511701 0.0181914739 Aug 1 5.383049 -0.001793076 0.0169508385 Sep 1 5.381589 -0.001778259 0.0184106500 Oct 1 5.481098 0.003093706 0.0189017443 Nov 1 5.581066 0.008071734 0.0189336381 Dec 1 5.680853 0.013054901 0.0191468731 Jan 2 5.977342 0.020912505 -0.0773419779 Feb 2 6.089175 0.027377371 0.0108248967 Mar 2 5.993510 0.019886785 0.0064896272 Apr 2 5.793285 0.006161774 0.0067151890 May 2 5.790403 0.005582290 0.0095966026 Jun 2 5.692607 -0.001193051 0.0073927501 Jul 2 5.492486 -0.014483134 0.0075138969 Aug 2 5.293722 -0.026998056 0.0062776466 Sep 2 5.191708 -0.032164431 0.0082920829 Oct 2 5.190634 -0.029997410 0.0093656237 Nov 2 4.993280 -0.041784374 0.0067197760 Dec 2 5.090564 -0.031915387 0.0094361509 Jan 3 5.169149 -0.024687178 -0.0691487381 Feb 3 5.190184 -0.021249754 0.0098155308 Mar 3 4.898281 -0.040875977 0.0017192240 Apr 3 4.795400 -0.045380868 0.0046004038 May 3 4.496139 -0.063908567 0.0038614191 Jun 3 4.490723 -0.059624654 0.0092770276 Jul 3 4.393499 -0.062386614 0.0065005677 Aug 3 4.393360 -0.057802254 0.0066400846 Sep 3 4.197242 -0.068011342 0.0027580923 Oct 3 4.089888 -0.070920721 0.0101120601 Nov 3 3.898152 -0.079870625 0.0018480172 Dec 3 3.793221 -0.081727882 0.0067787899 Jan 4 3.940154 -0.065155835 -0.0401539811 Feb 4 4.181872 -0.042124798 0.0181281523 Mar 4 4.101097 -0.045001756 -0.0010971958 Apr 4 3.795983 -0.064358036 0.0040173730 May 4 3.598207 -0.074293776 0.0017934132 Jun 4 3.685099 -0.062283240 0.0149014626 Jul 4 3.498796 -0.071528740 0.0012041736 Aug 4 3.391022 -0.074231893 0.0089784919 Sep 4 3.101357 -0.090304159 -0.0013569076 Oct 4 3.084576 -0.084817287 0.0154236743 Nov 4 3.094294 -0.077760073 0.0057062497 Dec 4 3.192986 -0.064593066 0.0070135416 Jan 5 3.352736 -0.047986244 -0.0527361886 Feb 5 3.480721 -0.034831549 0.0192793271 Mar 5 3.594402 -0.023738957 0.0055978996 Apr 5 3.495857 -0.029326050 0.0041430447 May 5 3.306722 -0.041263163 -0.0067215858 Jun 5 3.184197 -0.047333608 0.0158026625 Jul 5 3.100178 -0.050074373 -0.0001775521 Aug 5 3.182702 -0.040167516 0.0172979142 Sep 5 3.007278 -0.050273464 -0.0072780717 Oct 5 2.984312 -0.048233007 0.0156877178 Nov 5 3.093233 -0.036489362 0.0067671501 Dec 5 3.390049 -0.011593165 0.0099507856 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 1.529636038 -0.801081630 -0.799935279 -0.777645076 -0.755231688 2 2.215084192 0.507280713 -0.802957581 -1.434334208 -0.058875882 -0.672481894 3 0.781885054 0.271193289 -1.753526936 -0.401373510 -1.643316078 0.378560209 4 1.566517541 1.870044028 -0.250035773 -1.681359362 -0.862497210 1.042009655 5 1.512260977 1.088688007 0.960308896 -0.483341447 -1.032683735 -0.525128986 Jul Aug Sep Oct Nov Dec 1 -0.729844261 -1.356937535 0.002198096 0.666675182 0.636426541 0.601498344 2 -1.293076905 -1.197086405 -0.487029208 0.201751549 -1.085566915 0.901520025 3 -0.243303500 0.402756442 -0.894853207 -0.254511840 -0.781564316 -0.162036600 4 -0.801727117 -0.234306506 -1.392626468 0.475279264 0.611148466 1.140092338 5 -0.237075411 0.856887314 -0.874061654 0.176472542 1.015675249 2.152845446 > 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/1saqz1259917721.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/2jnhf1259917721.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/3ge7e1259917721.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/4h4ym1259917721.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/5nlbs1259917721.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/6m47u1259917721.tab") > system("convert tmp/1saqz1259917721.ps tmp/1saqz1259917721.png") > system("convert tmp/2jnhf1259917721.ps tmp/2jnhf1259917721.png") > system("convert tmp/3ge7e1259917721.ps tmp/3ge7e1259917721.png") > system("convert tmp/4h4ym1259917721.ps tmp/4h4ym1259917721.png") > system("convert tmp/5nlbs1259917721.ps tmp/5nlbs1259917721.png") > > > proc.time() user system elapsed 1.522 0.811 8.017