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(17,14,15,16,16,15,13,12,13,13,12,10,14,14,15,16,16,15,15,13,15,15,15,13,16,16,14,16,15,14,15,15,14,13,12,13,12,9,10,8,11,8,8,8,4,6,8,10,5,6,5,9,8,6,9,11,11,8,11,11,13) > 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.965577 0.000000 0.000000 1.367741 > m$fitted level slope sea Jan 1 17.000000 0.00000000 0.00000000 Feb 1 15.033059 -0.14454806 -0.14454806 Mar 1 15.036681 -0.13830763 -0.13830763 Apr 1 15.596131 -0.11941742 -0.11941742 May 1 15.835756 -0.11211844 -0.11211844 Jun 1 15.382783 -0.11800519 -0.11800519 Jul 1 14.069124 -0.13691827 -0.13691827 Aug 1 12.932596 -0.15197898 -0.15197898 Sep 1 12.987371 -0.14895565 -0.14895565 Oct 1 13.011304 -0.14648052 -0.14648052 Nov 1 12.464889 -0.15210857 -0.15210857 Dec 1 11.109800 -0.16878063 -0.16878063 Jan 2 11.375642 -0.20525124 2.25776363 Feb 2 13.038397 -0.14218902 -0.14218902 Mar 2 14.177455 -0.11744830 -0.11744830 Apr 2 15.214141 -0.10298988 -0.10298988 May 2 15.665204 -0.09772867 -0.09772867 Jun 2 15.305389 -0.09986134 -0.09986134 Jul 2 15.146598 -0.10030414 -0.10030414 Aug 2 13.961178 -0.10813117 -0.10813117 Sep 2 14.552863 -0.10319463 -0.10319463 Oct 2 14.814058 -0.10065938 -0.10065938 Nov 2 14.929325 -0.09917191 -0.09917191 Dec 2 13.865000 -0.10576642 -0.10576642 Jan 3 14.236578 -0.12489941 1.37389351 Feb 3 15.308478 -0.09847276 -0.09847276 Mar 3 14.580701 -0.10635954 -0.10635954 Apr 3 15.388350 -0.09891449 -0.09891449 May 3 15.183024 -0.09957345 -0.09957345 Jun 3 14.534366 -0.10249541 -0.10249541 Jul 3 14.806033 -0.10065202 -0.10065202 Aug 3 14.925892 -0.09960650 -0.09960650 Sep 3 14.420906 -0.10149073 -0.10149073 Oct 3 13.640003 -0.10461289 -0.10461289 Nov 3 12.737243 -0.10825297 -0.10825297 Dec 3 12.896343 -0.10704081 -0.10704081 Jan 4 11.851433 -0.08234735 0.90582089 Feb 4 10.198007 -0.10812237 -0.10812237 Mar 4 10.100173 -0.10802694 -0.10802694 Apr 4 8.937770 -0.11438969 -0.11438969 May 4 10.102878 -0.10851130 -0.10851130 Jun 4 8.941765 -0.11267274 -0.11267274 Jul 4 8.429299 -0.11413793 -0.11413793 Aug 4 8.202990 -0.11453403 -0.11453403 Sep 4 5.871073 -0.12221930 -0.12221930 Oct 4 5.957183 -0.12150459 -0.12150459 Nov 4 7.111112 -0.11715706 -0.11715706 Dec 4 8.736623 -0.11124492 -0.11124492 Jan 5 6.223469 -0.06415994 0.70575937 Feb 5 6.104203 -0.06487917 -0.06487917 Mar 5 5.490790 -0.06891836 -0.06891836 Apr 5 7.461868 -0.05914141 -0.05914141 May 5 7.769310 -0.05780272 -0.05780272 Jun 5 6.788410 -0.06070584 -0.06070584 Jul 5 8.029775 -0.05690694 -0.05690694 Aug 5 9.693994 -0.05206496 -0.05206496 Sep 5 10.428872 -0.04989051 -0.04989051 Oct 5 9.079218 -0.05344790 -0.05344790 Nov 5 10.157324 -0.05036861 -0.05036861 Dec 5 10.633431 -0.04894170 -0.04894170 Jan 6 11.536167 -0.06380147 0.70181615 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.39601525 0.13792829 0.68875869 0.36036526 -0.34406780 2 0.56769426 1.56498830 1.24240555 1.15631407 0.56004645 -0.26568395 3 0.55391490 1.07834053 -0.62087590 0.92079318 -0.10781667 -0.55737250 4 -1.04350229 -1.46319107 0.01023350 -1.06500714 1.29787006 -1.06918732 5 -2.61528643 -0.05230325 -0.54823485 2.06375594 0.37209581 -0.93800338 6 1.02265694 Jul Aug Sep Oct Nov Dec 1 -1.20961452 -1.01224113 0.20945813 0.17518912 -0.40530699 -1.21924732 2 -0.05980345 -1.10167800 0.71063547 0.37005724 0.21929512 -0.98023884 3 0.38007397 0.22405929 -0.41195292 -0.69046820 -0.81115894 0.27171585 4 -0.40630148 -0.11402124 -2.25415898 0.21179309 1.29666645 1.77169793 5 1.32363668 1.74992448 0.80016776 -1.32164912 1.15062568 0.53535269 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/1t07p1259928514.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/2lj9r1259928514.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/33q6n1259928514.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/435031259928514.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/5ukjq1259928514.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/64qwd1259928515.tab") > system("convert tmp/1t07p1259928514.ps tmp/1t07p1259928514.png") > system("convert tmp/2lj9r1259928514.ps tmp/2lj9r1259928514.png") > system("convert tmp/33q6n1259928514.ps tmp/33q6n1259928514.png") > system("convert tmp/435031259928514.ps tmp/435031259928514.png") > system("convert tmp/5ukjq1259928514.ps tmp/5ukjq1259928514.png") > > > proc.time() user system elapsed 1.371 0.791 2.654