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(96.96,93.11,95.62,98.30,96.38,100.82,99.06,94.03,102.07,99.31,98.64,101.82,99.14,97.63,100.06,101.32,101.49,105.43,105.09,99.48,108.53,104.34,106.10,107.35,103.00,104.50,105.17,104.84,106.18,108.86,107.77,102.74,112.63,106.26,108.86,111.38,106.85,107.86,107.94,111.38,111.29,113.72,111.88,109.87,113.72,111.71,114.81,112.05,111.54,110.87,110.87,115.48,111.63,116.24,113.56,106.01,110.45,107.77,108.61,108.19) > 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.37017880 0.01290006 1.16041900 0.00000000 > m$fitted level slope sea Jan 1 96.96000 0.000000000 0.000000000 Feb 1 95.69137 -0.000980021 -2.581366947 Mar 1 94.87369 -0.071519014 0.746314327 Apr 1 96.01534 0.052764979 2.284664081 May 1 96.44760 0.094045031 -0.067603171 Jun 1 98.16786 0.279267966 2.652140244 Jul 1 99.11862 0.359601884 -0.058621968 Aug 1 97.59607 0.123310476 -3.566066097 Sep 1 98.76500 0.260186108 3.305003207 Oct 1 99.39728 0.310521604 -0.087284825 Nov 1 99.36433 0.262934355 -0.724332243 Dec 1 100.32543 0.361272711 1.494568404 Jan 2 100.09352 0.281118250 -0.953522871 Feb 2 100.21969 0.259560722 -2.589687155 Mar 2 100.44422 0.254555729 -0.384224941 Apr 2 100.33000 0.201763459 0.989998438 May 2 101.21857 0.299678804 0.271428025 Jun 2 102.12512 0.386183344 3.304883714 Jul 2 103.31372 0.500789107 1.776278721 Aug 2 103.99209 0.526174495 -4.512085306 Sep 2 104.89593 0.580125421 3.634074882 Oct 2 105.13944 0.532121431 -0.799437528 Nov 2 106.10892 0.594354402 -0.008922449 Dec 2 106.31430 0.539123196 1.035695258 Jan 3 105.60964 0.362575165 -2.609643496 Feb 3 106.07637 0.377400403 -1.576373914 Mar 3 106.17651 0.337866658 -1.006510699 Apr 3 105.70708 0.222975900 -0.867077939 May 3 105.98498 0.230770091 0.195015682 Jun 3 106.21875 0.231194269 2.641253145 Jul 3 106.36914 0.219728657 1.400858631 Aug 3 106.89861 0.263739396 -4.158607067 Sep 3 107.88462 0.366387079 4.745378235 Oct 3 107.97990 0.327900004 -1.719904910 Nov 3 108.21150 0.314249277 0.648500391 Dec 3 108.88428 0.365019864 2.495716690 Jan 4 109.36292 0.381113108 -2.512924487 Feb 4 109.48080 0.343789774 -1.620796057 Mar 4 109.23067 0.259540915 -1.290671965 Apr 4 110.41000 0.389919243 0.969999592 May 4 111.14879 0.439307416 0.141209385 Jun 4 111.51767 0.429344262 2.202327670 Jul 4 111.51625 0.368385500 0.363745419 Aug 4 112.73347 0.488587857 -2.863466043 Sep 4 111.75376 0.280610168 1.966235729 Oct 4 112.22877 0.308132938 -0.518765533 Nov 4 113.29048 0.414747754 1.519520161 Dec 4 112.28772 0.214298741 -0.237721551 Jan 5 112.67075 0.238160947 -1.130748902 Feb 5 112.79709 0.222339536 -1.927086023 Mar 5 112.93762 0.210762398 -2.067618710 Apr 5 113.71911 0.291510083 1.760885466 May 5 113.14394 0.168957211 -1.513943927 Jun 5 113.31652 0.169469429 2.923475267 Jul 5 113.42974 0.161517491 0.130256475 Aug 5 111.48250 -0.136677085 -5.472499539 Sep 5 109.84610 -0.348782323 0.603903112 Oct 5 108.73696 -0.456301469 -0.966963979 Nov 5 107.36886 -0.585193834 1.241137762 Dec 5 107.27065 -0.516369493 0.919346461 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -2.354398700 -0.843449423 1.435496038 0.514396722 2.268052813 2 -0.779133951 -0.207945452 -0.044283221 -0.449635778 0.849420627 0.767563496 3 -1.573881484 0.132346390 -0.348798959 -1.003391964 0.068209721 0.003743282 4 0.142362656 -0.330251127 -0.742603906 1.144866038 0.433707370 -0.087767594 5 0.210503334 -0.139596124 -0.101994701 0.710247001 -1.077761346 0.004509969 Jul Aug Sep Oct Nov Dec 1 0.921625369 -2.530729566 1.384990018 0.487976288 -0.447416128 0.905238792 2 1.021987404 0.225129445 0.475465675 -0.421640736 0.546819406 -0.487828636 3 -0.101655098 0.389490496 0.903999093 -0.337764564 -0.119835733 0.447353214 4 -0.538398849 1.061399709 -1.832172611 0.241957838 0.937191751 -1.765090808 5 -0.070101489 -2.629064211 -1.868261768 -0.946090585 -1.133975809 0.605936835 > 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/1xoj71259861147.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/2kgkz1259861147.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/3gp6u1259861147.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/4o6kq1259861147.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/5km221259861147.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/64lh41259861147.tab") > system("convert tmp/1xoj71259861147.ps tmp/1xoj71259861147.png") > system("convert tmp/2kgkz1259861147.ps tmp/2kgkz1259861147.png") > system("convert tmp/3gp6u1259861147.ps tmp/3gp6u1259861147.png") > system("convert tmp/4o6kq1259861147.ps tmp/4o6kq1259861147.png") > system("convert tmp/5km221259861147.ps tmp/5km221259861147.png") > > > proc.time() user system elapsed 1.400 0.846 1.631