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(11.1,10.9,10,9.2,9.2,9.5,9.6,9.5,9.1,8.9,9,10.1,10.3,10.2,9.6,9.2,9.3,9.4,9.4,9.2,9,9,9,9.8,10,9.8,9.3,9,9,9.1,9.1,9.1,9.2,8.8,8.3,8.4,8.1,7.7,7.9,7.9,8,7.9,7.6,7.1,6.8,6.5,6.9,8.2,8.7,8.3,7.9,7.5,7.8,8.3,8.4,8.2,7.7,7.2,7.3,8.1) > 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.000000000 0.146842615 0.001148952 0.000000000 > m$fitted level slope sea Jan 1 11.100000 0.000000000 0.0000000000 Feb 1 10.902314 -0.197877276 -0.0023139605 Mar 1 10.019231 -0.860783815 -0.0192310726 Apr 1 9.194827 -0.825853973 0.0051734651 May 1 9.178985 -0.047679531 0.0210152908 Jun 1 9.496613 0.303317255 0.0033865430 Jul 1 9.608752 0.119631753 -0.0087518179 Aug 1 9.505703 -0.094319525 -0.0057032555 Sep 1 9.107704 -0.386095917 -0.0077035236 Oct 1 8.894074 -0.220390061 0.0059261439 Nov 1 8.993383 0.086777247 0.0066171997 Dec 1 10.075875 1.043462157 0.0241245649 Jan 2 10.327244 0.282609815 -0.0272435729 Feb 2 10.193062 -0.118014592 0.0069375677 Mar 2 9.619067 -0.546087046 -0.0190671751 Apr 2 9.206940 -0.420475495 -0.0069396677 May 2 9.276231 0.038567885 0.0237690080 Jun 2 9.395143 0.113891419 0.0048571694 Jul 2 9.413780 0.024589848 -0.0137804576 Aug 2 9.202184 -0.196833787 -0.0021839783 Sep 2 9.006864 -0.195414658 -0.0068639367 Oct 2 8.973630 -0.043370782 0.0263700104 Nov 2 9.045622 0.064782357 -0.0456220277 Dec 2 9.728118 0.643878572 0.0718815678 Jan 3 10.041197 0.333836916 -0.0411966821 Feb 3 9.787153 -0.217578113 0.0128467301 Mar 3 9.317266 -0.451057299 -0.0172659525 Apr 3 9.022764 -0.305990234 -0.0227641075 May 3 8.972163 -0.069601753 0.0278366953 Jun 3 9.090532 0.104439149 0.0094676216 Jul 3 9.106463 0.022487413 -0.0064632585 Aug 3 9.100652 -0.003713984 -0.0006524485 Sep 3 9.206214 0.097463713 -0.0062139957 Oct 3 8.785127 -0.382662661 0.0148728247 Nov 3 8.391352 -0.392951971 -0.0913518939 Dec 3 8.310389 -0.104095019 0.0896106792 Jan 4 8.119105 -0.184804100 -0.0191048360 Feb 4 7.675970 -0.424043551 0.0240297316 Mar 4 7.882943 0.155900568 0.0170570248 Apr 4 7.941965 0.066705797 -0.0419647172 May 4 7.991337 0.050766260 0.0086627478 Jun 4 7.890051 -0.089097312 0.0099493333 Jul 4 7.611332 -0.263529391 -0.0113318137 Aug 4 7.137520 -0.456964842 -0.0375195839 Sep 4 6.759896 -0.383980961 0.0401037883 Oct 4 6.467508 -0.299726701 0.0324919456 Nov 4 6.978626 0.446157527 -0.0786258881 Dec 4 8.080070 1.048886109 0.1199296440 Jan 5 8.708368 0.662062153 -0.0083676482 Feb 5 8.383140 -0.246030203 -0.0831399056 Mar 5 7.883942 -0.477839829 0.0160582798 Apr 5 7.534192 -0.360374393 -0.0341917886 May 5 7.755099 0.172159058 0.0449014906 Jun 5 8.267971 0.484346706 0.0320294025 Jul 5 8.417286 0.177335317 -0.0172862775 Aug 5 8.256569 -0.132443127 -0.0565692057 Sep 5 7.653603 -0.563607654 0.0463972224 Oct 5 7.198769 -0.463931918 0.0012311188 Nov 5 7.414413 0.158801166 -0.1144127064 Dec 5 7.965757 0.518465225 0.1342430486 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.519268297 -1.757604254 0.091073383 2.030830003 0.915960879 2 -1.986254839 -1.049652127 -1.127976783 0.327427338 1.198063388 0.196563697 3 -0.809447418 -1.441598780 -0.612360285 0.378321768 0.616906128 0.454181548 4 -0.210734474 -0.624685550 1.517415880 -0.232720466 -0.041592466 -0.364999189 5 -1.010028217 -2.369917420 -0.605824340 0.306566532 1.389467349 0.814719275 Jul Aug Sep Oct Nov Dec 1 -0.479345722 -0.558327255 -0.761419674 0.432426001 0.801583809 2.496564945 2 -0.233041580 -0.577827100 0.003703359 0.396773704 0.282236524 1.511212181 3 -0.213861223 -0.068375181 0.264033327 -1.252937821 -0.026851040 0.753803795 4 -0.455196599 -0.504789479 0.190458740 0.219869964 1.946466962 1.572899777 5 -0.801173794 -0.808398085 -1.125166977 0.260114050 1.625088810 0.938609053 > 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/1fnv51259847034.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/207f21259847034.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/3wxfn1259847034.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/4dmpz1259847034.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/5q4pe1259847034.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/6rgdq1259847034.tab") > > system("convert tmp/1fnv51259847034.ps tmp/1fnv51259847034.png") > system("convert tmp/207f21259847034.ps tmp/207f21259847034.png") > system("convert tmp/3wxfn1259847034.ps tmp/3wxfn1259847034.png") > system("convert tmp/4dmpz1259847034.ps tmp/4dmpz1259847034.png") > system("convert tmp/5q4pe1259847034.ps tmp/5q4pe1259847034.png") > > > proc.time() user system elapsed 1.422 0.821 2.313