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.4816,1.4562,1.4268,1.4088,1.4016,1.3650,1.3190,1.3050,1.2785,1.3239,1.3449,1.2732,1.3322,1.4369,1.4975,1.5770,1.5553,1.5557,1.5750,1.5527,1.4748,1.4718,1.4570,1.4684,1.4227,1.3896,1.3622,1.3716,1.3419,1.3511,1.3516,1.3242,1.3074,1.2999,1.3213,1.2881,1.2611,1.2727,1.2811,1.2684,1.2650,1.2770,1.2271,1.2020,1.1938,1.2103,1.1856,1.1786,1.2015,1.2256,1.2292,1.2037,1.2165,1.2694,1.2938,1.3201,1.3014,1.3119,1.3408,1.2991) > 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.001149493 0.000000000 0.000000000 0.000000000 > m$fitted level slope sea Jan 1 1.481600 0.0000000000 0.0000000000 Feb 1 1.457521 -0.0013207999 -0.0013208000 Mar 1 1.428233 -0.0014326693 -0.0014326693 Apr 1 1.410298 -0.0014984126 -0.0014984126 May 1 1.403121 -0.0015209486 -0.0015209486 Jun 1 1.366659 -0.0016590551 -0.0016590551 Jul 1 1.320833 -0.0018329411 -0.0018329411 Aug 1 1.306880 -0.0018804687 -0.0018804687 Sep 1 1.280476 -0.0019762645 -0.0019762645 Oct 1 1.325693 -0.0017926356 -0.0017926356 Nov 1 1.346605 -0.0017046331 -0.0017046331 Dec 1 1.275174 -0.0019738461 -0.0019738461 Jan 2 1.298646 -0.0030503936 0.0335543298 Feb 2 1.436538 0.0003616363 0.0003616365 Mar 2 1.497029 0.0004709619 0.0004709619 Apr 2 1.576386 0.0006141304 0.0006141304 May 2 1.554726 0.0005737794 0.0005737794 Jun 2 1.555127 0.0005734657 0.0005734657 Jul 2 1.574393 0.0006072072 0.0006072072 Aug 2 1.552134 0.0005660072 0.0005660072 Sep 2 1.474375 0.0004251346 0.0004251346 Oct 2 1.471381 0.0004189964 0.0004189964 Nov 2 1.456608 0.0003917710 0.0003917710 Dec 2 1.467989 0.0004114286 0.0004114286 Jan 3 1.434542 0.0010765051 -0.0118415563 Feb 3 1.389244 0.0003562353 0.0003562353 Mar 3 1.361876 0.0003236193 0.0003236193 Apr 3 1.371266 0.0003342723 0.0003342723 May 3 1.341601 0.0002990621 0.0002990621 Jun 3 1.350791 0.0003094848 0.0003094848 Jul 3 1.351290 0.0003097076 0.0003097076 Aug 3 1.323923 0.0002773364 0.0002773364 Sep 3 1.307143 0.0002574096 0.0002574096 Oct 3 1.299652 0.0002483683 0.0002483683 Nov 3 1.321027 0.0002729918 0.0002729918 Dec 3 1.287866 0.0002340698 0.0002340698 Jan 4 1.266687 0.0005079511 -0.0055874617 Feb 4 1.272136 0.0005644348 0.0005644349 Mar 4 1.280529 0.0005712424 0.0005712424 Apr 4 1.267840 0.0005597222 0.0005597222 May 4 1.264444 0.0005562879 0.0005562879 Jun 4 1.276434 0.0005662045 0.0005662045 Jul 4 1.226577 0.0005225108 0.0005225108 Aug 4 1.201500 0.0005003460 0.0005003460 Sep 4 1.193307 0.0004928263 0.0004928263 Oct 4 1.209793 0.0005066494 0.0005066494 Nov 4 1.185115 0.0004849008 0.0004849008 Dec 4 1.178122 0.0004784483 0.0004784483 Jan 5 1.204103 0.0002366120 -0.0026027322 Feb 5 1.225175 0.0004251034 0.0004251035 Mar 5 1.228773 0.0004272915 0.0004272915 Apr 5 1.203291 0.0004094352 0.0004094352 May 5 1.216082 0.0004179628 0.0004179628 Jun 5 1.268946 0.0004540578 0.0004540578 Jul 5 1.293329 0.0004705154 0.0004705154 Aug 5 1.319612 0.0004882555 0.0004882555 Sep 5 1.300925 0.0004750858 0.0004750858 Oct 5 1.311418 0.0004819616 0.0004819616 Nov 5 1.340299 0.0005014393 0.0005014393 Dec 5 1.298627 0.0004725342 0.0004725342 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.426435187 -0.826541507 -0.487681101 -0.167835020 -1.032614066 2 0.908524197 3.543725884 1.775110959 2.328840524 -0.657557496 -0.005120975 3 -1.101268776 -1.243316484 -0.818185983 0.267549708 -0.885338170 0.262378439 4 -0.674528375 0.136116288 0.231008727 -0.391263915 -0.116740916 0.337384762 5 0.790458499 0.582845714 0.093610957 -0.764459253 0.365332573 1.547419557 Jul Aug Sep Oct Nov Dec 1 -1.305264436 -0.358164478 -0.724736268 1.394647063 0.670967232 -2.060531676 2 0.551839260 -0.675037707 -2.312267872 -0.100933408 -0.448480952 0.324396484 3 0.005615939 -0.816817269 -0.503400217 -0.228670573 0.623501828 -0.986708072 4 -1.487850682 -0.755406214 -0.256504835 0.471925942 -0.743146566 -0.220671152 5 0.706040201 0.761576162 -0.565761248 0.295582333 0.837899644 -1.244300496 > 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/1xayz1259929534.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/2uikb1259929534.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/3451r1259929534.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/4cjym1259929534.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/53fxe1259929534.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/6qcr91259929534.tab") > system("convert tmp/1xayz1259929534.ps tmp/1xayz1259929534.png") > system("convert tmp/2uikb1259929534.ps tmp/2uikb1259929534.png") > system("convert tmp/3451r1259929534.ps tmp/3451r1259929534.png") > system("convert tmp/4cjym1259929534.ps tmp/4cjym1259929534.png") > system("convert tmp/53fxe1259929534.ps tmp/53fxe1259929534.png") > > > proc.time() user system elapsed 1.503 0.825 3.423