R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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.0622,1.0773,1.0807,1.0848,1.1582,1.1663,1.1372,1.1139,1.1222,1.1692,1.1702,1.2286,1.2613,1.2646,1.2262,1.1985,1.2007,1.2138,1.2266,1.2176,1.2218,1.249,1.2991,1.3408,1.3119,1.3014,1.3201,1.2938,1.2694,1.2165,1.2037,1.2292,1.2256,1.2015,1.1786,1.1856,1.2103,1.1938,1.202,1.2271,1.277,1.265,1.2684,1.2811,1.2727,1.2611,1.2881,1.3213,1.2999,1.3074,1.3242,1.3516,1.3511,1.3419,1.3716,1.3622,1.3896,1.4227,1.4684,1.457) > 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 6.853887e-04 8.690119e-07 0.000000e+00 0.000000e+00 > m$fitted level slope sea Jan 1 1.062200 0.0000000000 0.000000e+00 Feb 1 1.076515 0.0007910781 7.848775e-04 Mar 1 1.079905 0.0008076869 7.950848e-04 Apr 1 1.083992 0.0008326301 8.078683e-04 May 1 1.157113 0.0014687981 1.087164e-03 Jun 1 1.165188 0.0015346778 1.112432e-03 Jul 1 1.136203 0.0011952883 9.969909e-04 Aug 1 1.112994 0.0008965452 9.058105e-04 Sep 1 1.121267 0.0009949000 9.330027e-04 Oct 1 1.168100 0.0016548039 1.099552e-03 Nov 1 1.169103 0.0016447382 1.097218e-03 Dec 1 1.227304 0.0025736869 1.296215e-03 Jan 2 1.261107 0.0018242590 1.932781e-04 Feb 2 1.264578 0.0018843074 2.200655e-05 Mar 2 1.226241 0.0011753565 -4.051394e-05 Apr 2 1.198584 0.0006406646 -8.449812e-05 May 2 1.200782 0.0006709178 -8.216896e-05 Jun 2 1.213864 0.0009226146 -6.397972e-05 Jul 2 1.226647 0.0011728148 -4.696401e-05 Aug 2 1.217661 0.0009505867 -6.121936e-05 Sep 2 1.221857 0.0010239943 -5.676877e-05 Oct 2 1.249022 0.0016339733 -2.175198e-05 Nov 2 1.299058 0.0027962681 4.152828e-05 Dec 2 1.340709 0.0037543780 9.107245e-05 Jan 3 1.318536 0.0034874635 -6.635729e-03 Feb 3 1.301115 0.0027472415 2.846595e-04 Mar 3 1.319803 0.0031635038 2.967625e-04 Apr 3 1.293525 0.0023788127 2.750046e-04 May 3 1.269144 0.0016519488 2.557661e-04 Jun 3 1.216282 0.0001448823 2.176574e-04 Jul 3 1.203491 -0.0002186456 2.088684e-04 Aug 3 1.228974 0.0005146897 2.258324e-04 Sep 3 1.225377 0.0003956914 2.231968e-04 Oct 3 1.201292 -0.0003221220 2.079665e-04 Nov 3 1.178406 -0.0009918874 1.943451e-04 Dec 3 1.185401 -0.0007520951 1.990220e-04 Jan 4 1.208984 -0.0001834807 1.315753e-03 Feb 4 1.194083 -0.0007090616 -2.828815e-04 Mar 4 1.202279 -0.0004343554 -2.794486e-04 Apr 4 1.227370 0.0003598619 -2.699157e-04 May 4 1.277252 0.0019132777 -2.520004e-04 Jun 4 1.265257 0.0014737032 -2.568729e-04 Jul 4 1.268656 0.0015349908 -2.562198e-04 Aug 4 1.281353 0.0018925464 -2.525554e-04 Sep 4 1.272956 0.0015609217 -2.558246e-04 Oct 4 1.261360 0.0011344723 -2.598695e-04 Nov 4 1.288352 0.0019770116 -2.521789e-04 Dec 4 1.321543 0.0029990640 -2.431994e-04 Jan 5 1.299900 0.0022757534 -3.770936e-08 Feb 5 1.307355 0.0024609642 4.508282e-05 Mar 5 1.324152 0.0029360672 4.791491e-05 Apr 5 1.351547 0.0037496729 5.258610e-05 May 5 1.351048 0.0036078483 5.180174e-05 Jun 5 1.341850 0.0031790269 4.951695e-05 Jul 5 1.371546 0.0040696571 5.408914e-05 Aug 5 1.362148 0.0036160476 5.184518e-05 Sep 5 1.389544 0.0044190992 5.567365e-05 Oct 5 1.422640 0.0053898507 6.013410e-05 Nov 5 1.468334 0.0067573050 6.619048e-05 Dec 5 1.456936 0.0061400490 6.355518e-05 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.32824024 0.09933580 0.12527986 2.75969375 0.25203200 2 1.42961647 0.05323408 -1.52514772 -1.09269821 0.05898162 0.46992451 3 -1.07613398 -0.71837094 0.60134815 -1.11033574 -1.00889811 -2.05474731 4 0.97593954 -0.51879636 0.33501363 0.96005265 1.86239025 -0.52298563 5 -0.97141140 0.18516266 0.53856080 0.91878364 -0.15959397 -0.48096435 Jul Aug Sep Oct Nov Dec 1 -1.16365903 -0.92992837 0.28090653 1.74461729 -0.02481869 2.15007733 2 0.44887897 -0.38430537 0.12270850 0.98813174 1.82893365 1.46760104 3 -0.48746596 0.96827331 -0.15488025 -0.92185568 -0.84952386 0.30065315 4 0.07239918 0.41958772 -0.38676104 -0.49449872 0.97176622 1.17294858 5 0.99587189 -0.50576570 0.89301155 1.07683055 1.51339308 -0.68166757 > 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/wessaorg/rcomp/tmp/197bz1322430404.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2walv1322430404.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3d61n1322430404.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4z2w01322430404.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/50v9c1322430404.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/62xvc1322430404.tab") > > try(system("convert tmp/197bz1322430404.ps tmp/197bz1322430404.png",intern=TRUE)) character(0) > try(system("convert tmp/2walv1322430404.ps tmp/2walv1322430404.png",intern=TRUE)) character(0) > try(system("convert tmp/3d61n1322430404.ps tmp/3d61n1322430404.png",intern=TRUE)) character(0) > try(system("convert tmp/4z2w01322430404.ps tmp/4z2w01322430404.png",intern=TRUE)) character(0) > try(system("convert tmp/50v9c1322430404.ps tmp/50v9c1322430404.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.770 0.311 2.096