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.59,1.26,1.13,1.92,2.61,2.26,2.41,2.26,2.03,2.86,2.55,2.27,2.26,2.57,3.07,2.76,2.51,2.87,3.14,3.11,3.16,2.47,2.57,2.89,2.63,2.38,1.69,1.96,2.19,1.87,1.6,1.63,1.22,1.21,1.49,1.64,1.66,1.77,1.82,1.78,1.28,1.29,1.37,1.12,1.51,2.24,2.94,3.09,3.46,3.64,4.39,4.15,5.21,5.8,5.91,5.39,5.46,4.72,3.14,2.63) > 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.15559467 0.01496316 0.00000000 0.00000000 > m$fitted level slope sea Jan 1 1.590000 0.00000000 0.000000e+00 Feb 1 1.276641 -0.02661302 -1.664148e-02 Mar 1 1.146589 -0.04267860 -1.658931e-02 Apr 1 1.936925 0.12468757 -1.692502e-02 May 1 2.627101 0.25419385 -1.710073e-02 Jun 1 2.276959 0.10590744 -1.695903e-02 Jul 1 2.426967 0.11713425 -1.696673e-02 Aug 1 2.276932 0.04776197 -1.693216e-02 Sep 1 2.046906 -0.02513937 -1.690565e-02 Oct 1 2.876966 0.20058577 -1.696573e-02 Nov 1 2.566939 0.06539468 -1.693935e-02 Dec 1 2.286926 -0.02620949 -1.692625e-02 Jan 2 2.126972 -0.06046113 1.330276e-01 Feb 2 2.570832 0.06686880 -8.316638e-04 Mar 2 3.071362 0.18226226 -1.362313e-03 Apr 2 2.760920 0.05127101 -9.197024e-04 May 2 2.510721 -0.02884602 -7.208550e-04 Jun 2 2.870909 0.07452413 -9.092772e-04 Jul 2 3.140979 0.12647941 -9.788226e-04 Aug 2 3.110938 0.08489312 -9.379465e-04 Sep 2 3.160931 0.07562036 -9.312538e-04 Oct 2 2.470823 -0.12783552 -8.234285e-04 Nov 2 2.570847 -0.06729149 -8.469888e-04 Dec 2 2.890876 0.03562482 -8.763959e-04 Jan 3 2.668534 -0.03188621 -3.853366e-02 Feb 3 2.380178 -0.09781989 -1.780648e-04 Mar 3 1.689698 -0.25548614 3.015178e-04 Apr 3 1.960011 -0.11570205 -1.084793e-05 May 3 2.190162 -0.02378651 -1.617067e-04 Jun 3 1.870067 -0.10252352 -6.680350e-05 Jul 3 1.600027 -0.14703448 -2.740683e-05 Aug 3 1.630058 -0.09998680 -5.798454e-05 Sep 3 1.220019 -0.18237079 -1.866802e-05 Oct 3 1.210035 -0.13656538 -3.471931e-05 Nov 3 1.490063 -0.02586958 -6.320209e-05 Dec 3 1.640072 0.02086474 -7.203175e-05 Jan 4 1.659396 0.02045926 6.043374e-04 Feb 4 1.769110 0.04360774 8.895334e-04 Mar 4 1.819114 0.04530881 8.857012e-04 Apr 4 1.779077 0.02262204 9.232550e-04 May 4 1.278908 -0.11631283 1.092158e-03 Jun 4 1.288938 -0.08274005 1.062186e-03 Jul 4 1.368966 -0.03948975 1.033833e-03 Aug 4 1.118939 -0.09543248 1.060762e-03 Sep 4 1.508985 0.03356646 1.015166e-03 Oct 4 2.239033 0.21863384 9.671336e-04 Nov 4 2.939057 0.34654898 9.427565e-04 Dec 4 3.089050 0.29431957 9.500650e-04 Jan 5 3.463217 0.31538200 -3.217428e-03 Feb 5 3.640862 0.27947681 -8.617459e-04 Mar 5 4.391086 0.40465219 -1.085718e-03 Apr 5 4.150860 0.23324272 -8.604507e-04 May 5 5.211073 0.45301056 -1.072557e-03 Jun 5 5.801098 0.48941939 -1.098361e-03 Jul 5 5.911046 0.38858623 -1.045885e-03 Aug 5 5.390954 0.14713430 -9.536158e-04 Sep 5 5.460948 0.12663676 -9.478642e-04 Oct 5 4.720900 -0.10365882 -9.004148e-04 Nov 5 3.140841 -0.49597139 -8.410631e-04 Dec 5 2.630841 -0.49969924 -8.406489e-04 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.46894988 -0.24087732 1.88692291 1.25832719 -1.33054365 2 -0.33720710 0.91581390 0.94047396 -1.06909787 -0.65437801 0.84464880 3 -0.60507642 -0.50086801 -1.28597156 1.14132478 0.75090999 -0.64344467 4 -0.00352566 0.17972277 0.01388262 -0.18529477 -1.13523333 0.27438480 5 0.18036708 -0.28207680 1.02193085 -1.40025567 1.79590028 0.29757970 Jul Aug Sep Oct Nov Dec 1 0.09650662 -0.58269168 -0.60473942 1.85989903 -1.10989284 -0.75058201 2 0.42462574 -0.33992121 -0.07579929 -1.66318697 0.49493740 0.84133209 3 -0.36380715 0.38457513 -0.67345212 0.37444837 0.90492455 0.38205020 4 0.35352097 -0.45729704 1.05452368 1.51289389 1.04569500 -0.42697324 5 -0.82421825 -1.97375094 -0.16756203 -1.88263533 -3.20712862 -0.03047507 > 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/1e44x1259865983.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/2xqvi1259865983.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/3kdee1259865983.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/471ak1259865983.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/59toq1259865983.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/6yjst1259865983.tab") > system("convert tmp/1e44x1259865983.ps tmp/1e44x1259865983.png") > system("convert tmp/2xqvi1259865983.ps tmp/2xqvi1259865983.png") > system("convert tmp/3kdee1259865983.ps tmp/3kdee1259865983.png") > system("convert tmp/471ak1259865983.ps tmp/471ak1259865983.png") > system("convert tmp/59toq1259865983.ps tmp/59toq1259865983.png") > > > proc.time() user system elapsed 1.247 0.837 1.781