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(6.3,6.2,6.1,6.3,6.5,6.6,6.5,6.2,6.2,5.9,6.1,6.1,6.1,6.1,6.1,6.4,6.7,6.9,7,7,6.8,6.4,5.9,5.5,5.5,5.6,5.8,5.9,6.1,6.1,6,6,5.9,5.5,5.6,5.4,5.2,5.2,5.2,5.5,5.8,5.8,5.5,5.3,5.1,5.2,5.8,5.8,5.5,5,4.9,5.3,6.1,6.5,6.8,6.6,6.4,6.4,6.6) > 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.0000000000 0.0583044545 0.0001326111 0.0000000000 > m$fitted level slope sea Jan 1 6.300000 0.000000e+00 0.000000e+00 Feb 1 6.200339 -9.968924e-02 -3.387563e-04 Mar 1 6.099767 -1.005625e-01 2.325517e-04 Apr 1 6.297453 1.937710e-01 2.546944e-03 May 1 6.500628 2.030540e-01 -6.282820e-04 Jun 1 6.600972 1.016653e-01 -9.717137e-04 Jul 1 6.501583 -9.680078e-02 -1.582665e-03 Aug 1 6.201387 -2.975773e-01 -1.386981e-03 Sep 1 6.197052 -8.109658e-03 2.947695e-03 Oct 1 5.903191 -2.901828e-01 -3.190970e-03 Nov 1 6.095224 1.858253e-01 4.776051e-03 Dec 1 6.102670 9.742496e-03 -2.670127e-03 Jan 2 6.099611 -2.893737e-03 3.894338e-04 Feb 2 6.097851 -1.774059e-03 2.149064e-03 Mar 2 6.103277 5.270743e-03 -3.277349e-03 Apr 2 6.393558 2.835247e-01 6.442207e-03 May 2 6.701417 3.072866e-01 -1.417049e-03 Jun 2 6.902393 2.034741e-01 -2.393080e-03 Jul 2 6.997518 9.767195e-02 2.481573e-03 Aug 2 7.009323 1.382265e-02 -9.322807e-03 Sep 2 6.792473 -2.114284e-01 7.527321e-03 Oct 2 6.412704 -3.758119e-01 -1.270406e-02 Nov 2 5.894675 -5.146865e-01 5.324975e-03 Dec 2 5.500161 -3.973390e-01 -1.607017e-04 Jan 3 5.494000 -1.538731e-02 6.000400e-03 Feb 3 5.592981 9.631440e-02 7.018590e-03 Mar 3 5.807498 2.109315e-01 -7.498166e-03 Apr 3 5.898435 9.470249e-02 1.564568e-03 May 3 6.097479 1.957742e-01 2.520525e-03 Jun 3 6.104744 1.316191e-02 -4.744150e-03 Jul 3 6.000045 -1.010117e-01 -4.468809e-05 Aug 3 6.003218 -8.700507e-05 -3.217785e-03 Sep 3 5.893785 -1.060116e-01 6.215387e-03 Oct 3 5.513964 -3.712532e-01 -1.396360e-02 Nov 3 5.581733 5.403306e-02 1.826693e-02 Dec 3 5.414214 -1.605858e-01 -1.421379e-02 Jan 4 5.195623 -2.167707e-01 4.376981e-03 Feb 4 5.193606 -8.698717e-03 6.394428e-03 Mar 4 5.201907 7.684707e-03 -1.907476e-03 Apr 4 5.498504 2.860154e-01 1.495823e-03 May 4 5.794699 2.958213e-01 5.300630e-03 Jun 4 5.805822 2.155808e-02 -5.822432e-03 Jul 4 5.511914 -2.823456e-01 -1.191357e-02 Aug 4 5.302005 -2.125639e-01 -2.005160e-03 Sep 4 5.076915 -2.246308e-01 2.308487e-02 Oct 4 5.223093 1.325849e-01 -2.309338e-02 Nov 4 5.767855 5.296534e-01 3.214489e-02 Dec 4 5.822502 7.206416e-02 -2.250194e-02 Jan 5 5.510274 -2.981069e-01 -1.027369e-02 Feb 5 4.991254 -5.109495e-01 8.746056e-03 Mar 5 4.898855 -1.092616e-01 1.145324e-03 Apr 5 5.292385 3.733239e-01 7.615124e-03 May 5 6.082276 7.730637e-01 1.772431e-02 Jun 5 6.504813 4.366675e-01 -4.812872e-03 Jul 5 6.818672 3.188096e-01 -1.867230e-02 Aug 5 6.606041 -1.912064e-01 -6.041108e-03 Sep 5 6.379662 -2.249615e-01 2.033842e-02 Oct 5 6.438705 4.759396e-02 -3.870484e-02 Nov 5 6.557813 1.162252e-01 4.218725e-02 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.413526898 -0.003634748 1.218836123 0.038444902 -0.419893279 2 -0.052338210 0.004643509 0.029299828 1.152087498 0.098409103 -0.429931210 3 1.582025870 0.463081881 0.476226005 -0.481200502 0.418591076 -0.756273520 4 -0.232720493 0.862275194 0.068010707 1.152327772 0.040611883 -1.135837870 5 -1.533307882 -0.881785298 1.666340877 1.998084490 1.655528407 -1.393157229 Jul Aug Sep Oct Nov Dec 1 -0.821931127 -0.831499670 1.198806743 -1.168182892 1.971348722 -0.729232585 2 -0.438170815 -0.347255019 -0.932858742 -0.680781089 -0.575137685 0.485985066 3 -0.472840849 0.417971521 -0.438678121 -1.098476609 1.761288726 -0.888827211 4 -1.258592626 0.288995202 -0.049974202 1.479380004 1.644427276 -1.895070266 5 -0.488098870 -2.112190030 -0.139793929 1.128766312 0.284230858 > 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/1p9mh1259947365.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/2ypll1259947365.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/383ic1259947365.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/4lik11259947365.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/5ohd91259947365.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/69eyx1259947365.tab") > system("convert tmp/1p9mh1259947365.ps tmp/1p9mh1259947365.png") > system("convert tmp/2ypll1259947365.ps tmp/2ypll1259947365.png") > system("convert tmp/383ic1259947365.ps tmp/383ic1259947365.png") > system("convert tmp/4lik11259947365.ps tmp/4lik11259947365.png") > system("convert tmp/5ohd91259947365.ps tmp/5ohd91259947365.png") > > > proc.time() user system elapsed 1.489 0.783 2.317