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(14.9,18.6,19.1,18.8,18.2,18,19,20.7,21.2,20.7,19.6,18.6,18.7,23.8,24.9,24.8,23.8,22.3,21.7,20.7,19.7,18.4,17.4,17,18,23.8,25.5,25.6,23.7,22,21.3,20.7,20.4,20.3,20.4,19.8,19.5,23.1,23.5,23.5,22.9,21.9,21.5,20.5,20.2,19.4,19.2,18.8,18.8,22.6,23.3,23,21.4,19.9,18.8,18.6,18.4,18.6,19.9,19.2,18.4) > 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.306016620 0.963726631 0.002974124 0.000000000 > m$fitted level slope sea Jan 1 14.90000 0.00000000 0.00000000 Feb 1 18.49637 1.95952201 0.10363366 Mar 1 19.00180 0.78918907 0.09820006 Apr 1 18.68933 -0.08197203 0.11066739 May 1 18.08526 -0.49540241 0.11473653 Jun 1 17.88049 -0.26518039 0.11951262 Jul 1 18.87504 0.73280572 0.12496043 Aug 1 20.57900 1.50219178 0.12099595 Sep 1 21.09268 0.71905399 0.10731713 Oct 1 20.59060 -0.24837507 0.10939607 Nov 1 19.48750 -0.92552856 0.11250491 Dec 1 18.48257 -0.98843139 0.11743227 Jan 2 19.58025 0.63599789 -0.88024503 Feb 2 23.63367 2.96176953 0.16632953 Mar 2 24.78996 1.53109007 0.11004467 Apr 2 24.66723 0.23069543 0.13276572 May 2 23.65972 -0.74386970 0.14027803 Jun 2 22.15965 -1.33916482 0.14034848 Jul 2 21.54331 -0.77013665 0.15669352 Aug 2 20.54784 -0.94752408 0.15215852 Sep 2 19.55262 -0.98507401 0.14738076 Oct 2 18.25560 -1.23065006 0.14439947 Nov 2 17.23483 -1.06541244 0.16516766 Dec 2 16.85706 -0.52423050 0.14293972 Jan 3 19.35464 1.83813004 -1.35464241 Feb 3 23.59061 3.56245780 0.20938714 Mar 3 25.40347 2.18612133 0.09653337 Apr 3 25.47615 0.53269589 0.12384511 May 3 23.56598 -1.38012051 0.13401946 Jun 3 21.86128 -1.63432121 0.13872280 Jul 3 21.11221 -0.94102677 0.18779003 Aug 3 20.53934 -0.65270379 0.16066346 Sep 3 20.23411 -0.38057134 0.16588728 Oct 3 20.14368 -0.15334034 0.15631579 Nov 3 20.20303 0.01324441 0.19696822 Dec 3 19.70470 -0.38717631 0.09530322 Jan 4 21.12973 1.02794310 -1.62972660 Feb 4 22.86846 1.55189471 0.23153806 Mar 4 23.40663 0.75880252 0.09336524 Apr 4 23.33316 0.11021818 0.16684326 May 4 22.75831 -0.42402720 0.14168573 Jun 4 21.78908 -0.84930710 0.11092390 Jul 4 21.30698 -0.56287854 0.19302294 Aug 4 20.35552 -0.86597277 0.14447787 Sep 4 20.02757 -0.44630779 0.17242687 Oct 4 19.26709 -0.69137183 0.13290550 Nov 4 18.95592 -0.39478512 0.24408280 Dec 4 18.74315 -0.25292495 0.05685135 Jan 5 20.43098 1.25958563 -1.63097653 Feb 5 22.34426 1.74804911 0.25573618 Mar 5 23.22001 1.06846114 0.07998713 Apr 5 22.85269 -0.04719809 0.14731388 May 5 21.26792 -1.24258392 0.13208440 Jun 5 19.80295 -1.41552605 0.09705172 Jul 5 18.58016 -1.26564078 0.21984006 Aug 5 18.44992 -0.38269370 0.15007930 Sep 5 18.21659 -0.26653736 0.18340523 Oct 5 18.46512 0.13402324 0.13487714 Nov 5 19.60742 0.91817783 0.29257577 Dec 5 19.23388 -0.08534274 -0.03388363 Jan 6 20.08258 0.64114214 -1.68257648 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 2.67574541 -1.09926466 -0.88351553 -0.42108708 0.23451341 2 1.81717331 2.25596284 -1.41917815 -1.32112553 -0.99269215 -0.60639399 3 2.52028625 1.71548337 -1.38271892 -1.67984582 -1.94841303 -0.25894039 4 1.48271640 0.52597407 -0.80197511 -0.65889056 -0.54418400 -0.43320939 5 1.56953415 0.49258384 -0.68978078 -1.13333102 -1.21760145 -0.17616655 6 0.74942427 Jul Aug Sep Oct Nov Dec 1 1.01659380 0.78373164 -0.79773979 -0.98546725 -0.68977936 -0.06407568 2 0.57963791 -0.18069490 -0.03825007 -0.25015493 0.16831894 0.55128848 3 0.70622121 0.29369890 0.27720650 0.23146785 0.16969228 -0.40794516 4 0.29176917 -0.30874557 0.42748987 -0.24963347 0.30212102 0.14455400 5 0.15267996 0.89941005 0.11832213 0.40802994 0.79878638 -1.02281489 6 > 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/1m21n1259953504.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/2dbij1259953504.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/3aigr1259953504.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/47qxz1259953504.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/5ptz61259953505.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/66mto1259953505.tab") > system("convert tmp/1m21n1259953504.ps tmp/1m21n1259953504.png") > system("convert tmp/2dbij1259953504.ps tmp/2dbij1259953504.png") > system("convert tmp/3aigr1259953504.ps tmp/3aigr1259953504.png") > system("convert tmp/47qxz1259953504.ps tmp/47qxz1259953504.png") > system("convert tmp/5ptz61259953505.ps tmp/5ptz61259953505.png") > > > proc.time() user system elapsed 1.461 0.831 5.583