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(37,30,47,35,30,43,82,40,47,19,52,136,80,42,54,66,81,63,137,72,107,58,36,52,79,77,54,84,48,96,83,66,61,53,30,74,69,59,42,65,70,100,63,105,82,81,75,102,121,98,76,77,63,37,35,23,40,29,37,51,20,28,13,22,25,13,16,13,16,17,9,17,25,14,8,7,10,7,10,3) > 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 66.90099 0.00000 0.00000 308.92244 > m$fitted level slope sea Jan 1 37.000000 0.0000000000 0.0000000000 Feb 1 34.067122 -0.3156475648 -0.3156475642 Mar 1 39.275416 0.1130184938 0.1130184938 Apr 1 37.598824 0.0060913332 0.0060913332 May 1 34.632000 -0.1387019757 -0.1387019757 Jun 1 37.829326 0.0002859076 0.0002859076 Jul 1 54.734506 0.6263501415 0.6263501415 Aug 1 49.275525 0.4201618229 0.4201618229 Sep 1 48.514007 0.3827671799 0.3827671799 Oct 1 37.441873 0.0393381287 0.0393381287 Nov 1 42.949323 0.1961723047 0.1961723047 Dec 1 78.093308 1.1614128174 1.1614128174 Jan 2 82.405038 0.8036406695 -8.8400473598 Feb 2 65.051623 0.0542762837 0.0542762801 Mar 2 60.629162 -0.0882405584 -0.0882405584 Apr 2 62.690879 -0.0332794524 -0.0332794524 May 2 69.659962 0.1173098829 0.1173098829 Jun 2 67.173084 0.0681303204 0.0681303204 Jul 2 93.436365 0.5179984697 0.5179984697 Aug 2 85.531341 0.3829514108 0.3829514108 Sep 2 93.660794 0.5010822848 0.5010822848 Oct 2 80.458151 0.2997739527 0.2997739527 Nov 2 63.928272 0.0595876340 0.0595876340 Dec 2 59.490134 -0.0031229268 -0.0031229268 Jan 3 65.169043 -0.2969538192 3.2664920132 Feb 3 69.932722 -0.1633123201 -0.1633123225 Mar 3 63.696644 -0.2869490594 -0.2869490594 Apr 3 71.361924 -0.1564251609 -0.1564251609 May 3 62.518572 -0.2770143166 -0.2770143166 Jun 3 74.997892 -0.1206427913 -0.1206427913 Jul 3 77.957327 -0.0861183040 -0.0861183040 Aug 3 73.475654 -0.1323596319 -0.1323596319 Sep 3 68.793795 -0.1781142508 -0.1781142508 Oct 3 62.867856 -0.2340744045 -0.2340744045 Nov 3 50.575671 -0.3486724708 -0.3486724708 Dec 3 59.201727 -0.2649565248 -0.2649565248 Jan 4 61.189337 -0.3403683766 3.7440521439 Feb 4 60.250713 -0.3519846967 -0.3519846982 Mar 4 53.165493 -0.4527514787 -0.4527514787 Apr 4 57.526352 -0.3945665441 -0.3945665441 May 4 62.104399 -0.3436029335 -0.3436029335 Jun 4 76.167961 -0.2128651651 -0.2128651651 Jul 4 71.207738 -0.2523612601 -0.2523612601 Aug 4 83.716725 -0.1524613105 -0.1524613105 Sep 4 83.039365 -0.1563993677 -0.1563993677 Oct 4 82.241351 -0.1610706984 -0.1610706984 Nov 4 79.509800 -0.1793825794 -0.1793825794 Dec 4 87.817070 -0.1199051041 -0.1199051041 Jan 5 98.177637 -0.3793057676 4.1723634461 Feb 5 98.025247 -0.3758220022 -0.3758220051 Mar 5 89.551339 -0.4716322246 -0.4716322246 Apr 5 84.714285 -0.5133931937 -0.5133931937 May 5 76.468043 -0.5762076058 -0.5762076058 Jun 5 61.615472 -0.6790478345 -0.6790478345 Jul 5 51.542235 -0.7411861056 -0.7411861056 Aug 5 40.746767 -0.8038586130 -0.8038586130 Sep 5 40.262604 -0.8019459242 -0.8019459242 Oct 5 35.875456 -0.8227912675 -0.8227912675 Nov 5 36.080657 -0.8169345447 -0.8169345447 Dec 5 41.405656 -0.7824620547 -0.7824620547 Jan 6 30.815421 -0.5899077067 6.4889847728 Feb 6 29.592923 -0.5979386197 -0.5979386189 Mar 6 23.164359 -0.6549761831 -0.6549761831 Apr 6 22.563968 -0.6545441183 -0.6545441183 May 6 23.306545 -0.6451448059 -0.6451448059 Jun 6 19.306742 -0.6651794757 -0.6651794757 Jul 6 17.907289 -0.6692102689 -0.6692102689 Aug 6 15.913281 -0.6760705909 -0.6760705909 Sep 6 15.771059 -0.6734146373 -0.6734146373 Oct 6 16.053014 -0.6687910465 -0.6687910465 Nov 6 13.264747 -0.6788514581 -0.6788514581 Dec 6 14.474297 -0.6700127533 -0.6700127533 Jan 7 15.108447 -0.6912298405 7.6035282435 Feb 7 14.522470 -0.6900907828 -0.6900907812 Mar 7 11.892182 -0.7062690061 -0.7062690061 Apr 7 9.884129 -0.7150564370 -0.7150564370 May 7 9.744995 -0.7117498664 -0.7117498664 Jun 7 8.542055 -0.7142551131 -0.7142551131 Jul 7 8.899339 -0.7092279478 -0.7092279478 Aug 7 6.527869 -0.7165896890 -0.7165896890 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.29155821 0.57810794 -0.19482668 -0.33272002 0.38039034 2 0.53070282 -1.77265834 -0.48610517 0.24498872 0.81704541 -0.30766236 3 0.81616565 0.53927214 -0.68946181 0.92987301 -1.03081505 1.52551439 4 0.30676469 -0.06624858 -0.78005252 0.56963710 0.59475104 1.73294072 5 1.38956946 0.02567701 -0.94906246 -0.52017383 -0.92915174 -1.72303018 6 -1.27958777 -0.07258097 -0.68844640 0.00653360 0.16837829 -0.40577759 7 0.16832709 0.01219341 -0.23028752 -0.15630744 0.06955623 -0.05950736 Jul Aug Sep Oct Nov Dec 1 1.95163955 -0.70835479 -0.13831563 -1.34590728 0.64421889 4.12555502 2 3.11621286 -1.00607463 0.92755597 -1.64344570 -2.02043377 -0.54035233 3 0.36994385 -0.52928211 -0.54866103 -0.69384602 -1.45650351 1.08452121 4 -0.57283348 1.54266792 -0.06352345 -0.07770251 -0.31143571 1.02853627 5 -1.13661231 -1.21823775 0.03877004 -0.43501897 0.12477702 0.74567274 6 -0.08899897 -0.16076592 0.06482961 0.11606987 -0.25757200 0.22953254 7 0.13004202 -0.20193383 > 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 > postscript(file="/var/www/html/rcomp/tmp/1bzj51291132007.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/html/rcomp/tmp/2m80p1291132007.ps",horizontal=F,onefile=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/3m80p1291132007.ps",horizontal=F,onefile=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/4xh0s1291132007.ps",horizontal=F,onefile=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/57qzw1291132007.ps",horizontal=F,onefile=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/6b9y11291132007.tab") > > try(system("convert tmp/1bzj51291132007.ps tmp/1bzj51291132007.png",intern=TRUE)) character(0) > try(system("convert tmp/2m80p1291132007.ps tmp/2m80p1291132007.png",intern=TRUE)) character(0) > try(system("convert tmp/3m80p1291132007.ps tmp/3m80p1291132007.png",intern=TRUE)) character(0) > try(system("convert tmp/4xh0s1291132007.ps tmp/4xh0s1291132007.png",intern=TRUE)) character(0) > try(system("convert tmp/57qzw1291132007.ps tmp/57qzw1291132007.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.580 0.835 3.585