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(1.4,1.2,1,1.7,2.4,2,2.1,2,1.8,2.7,2.3,1.9,2,2.3,2.8,2.4,2.3,2.7,2.7,2.9,3,2.2,2.3,2.8,2.8,2.8,2.2,2.6,2.8,2.5,2.4,2.3,1.9,1.7,2,2.1,1.7,1.8,1.8,1.8,1.3,1.3,1.3,1.2,1.4,2.2,2.9,3.1,3.5,3.6,4.4,4.1,5.1,5.8,5.9,5.4,5.5,4.8,3.2,2.7) > 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.177451041 0.009938781 0.000000000 0.000000000 > m$fitted level slope sea Jan 1 1.400000 0.000000000 0.000000e+00 Feb 1 1.210215 -0.013779306 -1.021464e-02 Mar 1 1.010431 -0.032149922 -1.043086e-02 Apr 1 1.709695 0.065916366 -9.694610e-03 May 1 2.409159 0.167134887 -9.158758e-03 Jun 1 2.009553 0.066532635 -9.553016e-03 Jul 1 2.109534 0.072865664 -9.534153e-03 Aug 1 2.009612 0.038821321 -9.612396e-03 Sep 1 1.809699 -0.009392717 -9.698670e-03 Oct 1 2.709438 0.177050179 -9.437506e-03 Nov 1 2.309569 0.057603094 -9.568922e-03 Dec 1 1.909651 -0.037686912 -9.651435e-03 Jan 2 1.887415 -0.034599735 1.125853e-01 Feb 2 2.300153 0.056146654 -1.534731e-04 Mar 2 2.800492 0.149403094 -4.923038e-04 Apr 2 2.400161 0.033928269 -1.610485e-04 May 2 2.300097 0.005772572 -9.727436e-05 Jun 2 2.700246 0.088662522 -2.455275e-04 Jul 2 2.700219 0.070018776 -2.191962e-04 Aug 2 2.900250 0.097352414 -2.496808e-04 Sep 2 3.000250 0.097909191 -2.501712e-04 Oct 2 2.200119 -0.090921788 -1.188430e-04 Nov 2 2.300141 -0.050770260 -1.408947e-04 Dec 2 2.800191 0.065059632 -1.911307e-04 Jan 3 2.808968 0.053396871 -8.967593e-03 Feb 3 2.800085 0.040614627 -8.543282e-05 Mar 3 2.199732 -0.094267061 2.677832e-04 Apr 3 2.599947 0.009755996 5.261408e-05 May 3 2.800013 0.049783781 -1.277951e-05 Jun 3 2.499918 -0.023799145 8.216049e-05 Jul 3 2.399902 -0.039827597 9.849277e-05 Aug 3 2.299891 -0.052483718 1.086770e-04 Sep 3 1.899845 -0.125574176 1.551240e-04 Oct 3 1.699837 -0.141227195 1.629792e-04 Nov 3 1.999874 -0.048431074 1.262043e-04 Dec 3 2.099884 -0.017214296 1.164348e-04 Jan 4 1.744787 -0.087607103 -4.478651e-02 Feb 4 1.794464 -0.059252635 5.536058e-03 Mar 4 1.794488 -0.046780488 5.511720e-03 Apr 4 1.794503 -0.036936836 5.496547e-03 May 4 1.294385 -0.134355718 5.615141e-03 Jun 4 1.294412 -0.106093658 5.587970e-03 Jul 4 1.294429 -0.083778364 5.571027e-03 Aug 4 1.194427 -0.087190186 5.573072e-03 Sep 4 1.394456 -0.026788660 5.544472e-03 Oct 4 2.194521 0.147097281 5.479452e-03 Nov 4 2.894555 0.263379543 5.445115e-03 Dec 4 3.094552 0.250050143 5.448223e-03 Jan 5 3.541125 0.291092690 -4.112457e-02 Feb 5 3.598239 0.242584469 1.760938e-03 Mar 5 4.398421 0.359894499 1.579050e-03 Apr 5 4.098251 0.221053513 1.749074e-03 May 5 5.098409 0.384916322 1.590595e-03 Jun 5 5.798460 0.451192090 1.539975e-03 Jul 5 5.898415 0.377325814 1.584529e-03 Aug 5 5.398328 0.192805020 1.672423e-03 Sep 5 5.498320 0.173286575 1.679765e-03 Oct 5 4.798266 -0.010377367 1.734324e-03 Nov 5 3.198187 -0.344693252 1.812749e-03 Dec 5 2.698181 -0.377355849 1.818800e-03 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.267829279 -0.419696617 1.617461704 1.379884195 -1.221082256 2 0.038101678 0.795823002 0.936449825 -1.159084896 -0.282542240 0.831667719 3 -0.128662804 -0.118724164 -1.351200417 1.042580693 0.401305051 -0.737858471 4 -0.752598895 0.269399271 0.124981841 0.098678592 -0.976810310 0.283421970 5 0.431964445 -0.466574855 1.175787002 -1.391999703 1.643166745 0.664669528 Jul Aug Sep Oct Nov Dec 1 0.071537060 -0.367741955 -0.506485502 1.924821839 -1.219865618 -0.966600938 2 -0.187041520 0.274205023 0.005585254 -1.894192024 0.402759982 1.161878791 3 -0.160745685 -0.126934777 -0.733095901 -0.157004015 0.930786981 0.313121890 4 0.223805818 -0.034219958 0.605837552 1.744143520 1.166372112 -0.133702011 5 -0.740846503 -1.850745454 -0.195775649 -1.842234323 -3.353379268 -0.327626716 > 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/1n0nr1259949589.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/22ncu1259949589.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/3zb4a1259949589.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/4bgc21259949589.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/50vp51259949589.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/6slas1259949589.tab") > system("convert tmp/1n0nr1259949589.ps tmp/1n0nr1259949589.png") > system("convert tmp/22ncu1259949589.ps tmp/22ncu1259949589.png") > system("convert tmp/3zb4a1259949589.ps tmp/3zb4a1259949589.png") > system("convert tmp/4bgc21259949589.ps tmp/4bgc21259949589.png") > system("convert tmp/50vp51259949589.ps tmp/50vp51259949589.png") > > > proc.time() user system elapsed 1.307 0.806 1.530