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(2.05,2.11,2.09,2.05,2.08,2.06,2.06,2.08,2.07,2.06,2.07,2.06,2.09,2.07,2.09,2.28,2.33,2.35,2.52,2.63,2.58,2.70,2.81,2.97,3.04,3.28,3.33,3.50,3.56,3.57,3.69,3.82,3.79,3.96,4.06,4.05,4.03,3.94,4.02,3.88,4.02,4.03,4.09,3.99,4.01,4.01,4.19,4.30,4.27,3.82,3.15,2.49,1.81,1.26,1.06,0.84,0.78,0.70,0.36,0.35) > 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.000000e+00 1.348242e-02 7.380607e-06 1.042469e-03 > m$fitted level slope sea Jan 1 2.0500000 0.0000000000 0.000000e+00 Feb 1 2.1055023 0.0550869014 2.544225e-04 Mar 1 2.0937942 -0.0041317927 1.814800e-04 Apr 1 2.0518948 -0.0378711409 3.122225e-04 May 1 2.0759384 0.0174651610 4.533523e-04 Jun 1 2.0615893 -0.0109702547 2.648784e-04 Jul 1 2.0591370 -0.0033571117 3.666298e-04 Aug 1 2.0783148 0.0167839726 3.718876e-04 Sep 1 2.0711049 -0.0046610384 2.934215e-04 Oct 1 2.0600335 -0.0103904328 3.400911e-04 Nov 1 2.0685271 0.0064875724 3.723539e-04 Dec 1 2.0605341 -0.0064547356 3.097912e-04 Jan 2 2.0885588 0.0241023938 -5.512563e-04 Feb 2 2.0721715 -0.0085319196 -3.312868e-04 Mar 2 2.0889268 0.0141235063 -3.588774e-04 Apr 2 2.2707186 0.1635207441 -3.983669e-04 May 2 2.3365978 0.0763974264 -9.326571e-04 Jun 2 2.3542667 0.0239862962 -8.586549e-04 Jul 2 2.5124840 0.1437730646 -2.731860e-04 Aug 2 2.6321912 0.1222972414 -7.947779e-04 Sep 2 2.5906911 -0.0238716260 -1.186446e-03 Oct 2 2.6929038 0.0886432216 -2.201075e-04 Nov 2 2.8090204 0.1131599948 -6.146527e-04 Dec 2 2.9680347 0.1540778584 -6.953873e-04 Jan 3 3.0423045 0.0832559090 2.302156e-03 Feb 3 3.2712365 0.2036929037 1.702875e-03 Mar 3 3.3380145 0.0813145616 -2.366062e-04 Apr 3 3.4934640 0.1472913675 2.269970e-03 May 3 3.5633893 0.0783613619 1.081085e-03 Jun 3 3.5737373 0.0177555847 1.931775e-04 Jul 3 3.6828776 0.0991844224 1.841457e-03 Aug 3 3.8159481 0.1293783744 2.093715e-03 Sep 3 3.7989291 -0.0010677976 -4.692075e-04 Oct 3 3.9487738 0.1334020026 2.505303e-03 Nov 3 4.0605453 0.1141280075 7.047198e-04 Dec 3 4.0557535 0.0081737651 1.118044e-03 Jan 4 4.0460569 -0.0076876192 -1.502753e-02 Feb 4 3.9426789 -0.0880079825 2.134429e-03 Mar 4 4.0119301 0.0522816868 -8.709764e-04 Apr 4 3.8894139 -0.1030881251 6.105411e-04 May 4 4.0047445 0.0912411772 2.683737e-03 Jun 4 4.0351353 0.0370948645 -1.632435e-03 Jul 4 4.0883995 0.0514823121 6.697221e-04 Aug 4 3.9963023 -0.0762730635 1.962699e-03 Sep 4 4.0064213 0.0005975477 -1.394339e-03 Oct 4 4.0084328 0.0018556897 1.485784e-03 Nov 4 4.1786813 0.1516923193 1.625236e-03 Dec 4 4.3016151 0.1261066227 4.009996e-05 Jan 5 4.2895995 0.0035419241 -1.166204e-02 Feb 5 3.8465356 -0.3750977024 -3.525607e-03 Mar 5 3.1686744 -0.6447344580 -1.465856e-03 Apr 5 2.4942010 -0.6711364925 -2.501130e-03 May 5 1.8090673 -0.6835737196 1.735477e-03 Jun 5 1.2554136 -0.5681182898 -2.865799e-03 Jul 5 1.0370136 -0.2573449521 2.926981e-03 Aug 5 0.8370900 -0.2063187788 -3.835399e-04 Sep 5 0.7730472 -0.0798884664 -1.207899e-03 Oct 5 0.7010303 -0.0728935772 -1.481757e-03 Nov 5 0.3764757 -0.2965315402 -2.040567e-03 Dec 5 0.3345733 -0.0703012559 8.244792e-04 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.49419190 -0.51950475 -0.29104269 0.47657006 -0.24489367 2 0.27550988 -0.27496134 0.19156559 1.28429829 -0.75033040 -0.45137497 3 -0.62957154 1.02218876 -1.03999433 0.56734125 -0.59364419 -0.52194915 4 -0.13988841 -0.68433888 1.19608689 -1.33624691 1.67361715 -0.46631898 5 -1.07542349 -3.23373607 -2.30416984 -0.22708737 -0.10711280 0.99432532 Jul Aug Sep Oct Nov Dec 1 0.06556622 0.17345983 -0.18468956 -0.04934291 0.14535741 -0.11146225 2 1.03163201 -0.18495492 -1.25884121 0.96900475 0.21114432 0.35239467 3 0.70128446 0.26003753 -1.12343360 1.15808605 -0.16599229 -0.91250760 4 0.12390811 -1.10025982 0.66202807 0.01083542 1.29042967 -0.22035334 5 2.67645371 0.43944959 1.08884806 0.06024166 -1.92602590 1.94841070 > 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/1zylj1259959161.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/24xot1259959161.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/3diob1259959161.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/4q5ij1259959161.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/5uuy51259959161.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/6nykj1259959161.tab") > > system("convert tmp/1zylj1259959161.ps tmp/1zylj1259959161.png") > system("convert tmp/24xot1259959161.ps tmp/24xot1259959161.png") > system("convert tmp/3diob1259959161.ps tmp/3diob1259959161.png") > system("convert tmp/4q5ij1259959161.ps tmp/4q5ij1259959161.png") > system("convert tmp/5uuy51259959161.ps tmp/5uuy51259959161.png") > > > proc.time() user system elapsed 1.421 0.796 1.628