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(29.837,29.571,30.167,30.524,30.996,31.033,31.198,30.937,31.649,33.115,34.106,33.926,33.382,32.851,32.948,36.112,36.113,35.210,35.193,34.383,35.349,37.058,38.076,36.630,36.045,35.638,35.114,35.465,35.254,35.299,35.916,36.683,37.288,38.536,38.977,36.407,34.955,34.951,32.680,34.791,34.178,35.213,34.871,35.299,35.443,37.108,36.419,34.471,33.868,34.385,33.643,34.627,32.919,35.500,36.110,37.086,37.711,40.427,39.884,38.512,38.767) > 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 7.540313e-01 -2.700247e-17 3.490118e-02 0.000000e+00 > m$fitted level slope sea Jan 1 29.83700 0.0000000000 0.000000000 Feb 1 29.59041 -0.0133242641 -0.019410131 Mar 1 30.10473 -0.0059093844 0.062265129 Apr 1 30.50164 -0.0041790382 0.022358585 May 1 30.96004 -0.0021454240 0.035957500 Jun 1 31.04372 -0.0017382535 -0.010723917 Jul 1 31.18560 -0.0010532928 0.012401347 Aug 1 30.97037 -0.0020685124 -0.033372276 Sep 1 31.56833 0.0007606616 0.080674537 Oct 1 32.96697 0.0073207446 0.148025645 Nov 1 34.02695 0.0122377150 0.079048438 Dec 1 33.96629 0.0118987771 -0.040288600 Jan 2 33.64256 0.0253993901 -0.260563217 Feb 2 32.99635 0.0118802614 -0.145349415 Mar 2 32.88626 0.0104368688 0.061739088 Apr 2 35.60117 0.0204899103 0.510828307 May 2 36.12329 0.0216534240 -0.010286681 Jun 2 35.39599 0.0197223412 -0.185992964 Jul 2 35.12739 0.0189295563 0.065606961 Aug 2 34.51689 0.0171866650 -0.133894873 Sep 2 35.21295 0.0190679567 0.136045643 Oct 2 36.76444 0.0234988686 0.293562612 Nov 2 37.87891 0.0267545981 0.197085967 Dec 2 36.88816 0.0269270637 -0.258157832 Jan 3 36.32085 0.0360231237 -0.275849112 Feb 3 35.80284 0.0303057284 -0.164844123 Mar 3 35.49225 0.0269064084 -0.378245790 Apr 3 35.10784 0.0251432392 0.357158308 May 3 35.13095 0.0251389231 0.123051046 Jun 3 35.37861 0.0255792415 -0.079607167 Jul 3 35.66430 0.0261433222 0.251698146 Aug 3 36.68099 0.0283903192 0.002011696 Sep 3 37.30925 0.0298066248 -0.021251246 Oct 3 38.27833 0.0321700172 0.257672654 Nov 3 38.62875 0.0328204540 0.348254886 Dec 3 36.97027 0.0354773477 -0.563269616 Jan 4 35.46875 0.0472402792 -0.513754287 Feb 4 35.07879 0.0445433039 -0.127785084 Mar 4 33.31479 0.0296570254 -0.634792652 Apr 4 34.13840 0.0333005671 0.652601685 May 4 34.09061 0.0331142109 0.087392157 Jun 4 35.08981 0.0348343564 0.123189463 Jul 4 34.86651 0.0343493887 0.004494215 Aug 4 35.30579 0.0351771727 -0.006785450 Sep 4 35.61749 0.0357874675 -0.174493398 Oct 4 36.72534 0.0381806037 0.382659444 Nov 4 36.02351 0.0372140943 0.395488253 Dec 4 35.00223 0.0391053648 -0.531225222 Jan 5 34.36740 0.0419179791 -0.499397031 Feb 5 34.17557 0.0409393912 0.209434530 Mar 5 34.34824 0.0418350391 -0.705237660 Apr 5 34.02848 0.0401712742 0.598517659 May 5 33.17351 0.0379213263 -0.254510942 Jun 5 34.89107 0.0408777506 0.608929964 Jul 5 35.98349 0.0427021565 0.126511634 Aug 5 36.99888 0.0445614216 0.087124197 Sep 5 37.95197 0.0464443934 -0.240974096 Oct 5 39.61052 0.0495555496 0.816480863 Nov 5 39.43689 0.0493750270 0.447108463 Dec 5 39.05781 0.0500524704 -0.545813671 Jan 6 39.22681 0.0497749785 -0.459806089 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.172678862 0.605236884 0.465018942 0.533814452 0.099037846 2 -0.454947411 -0.682761335 -0.138488357 3.119824864 0.578398052 -0.863423460 3 -0.728019555 -0.604216926 -0.385678633 -0.474070437 -0.002347474 0.256492011 4 -1.822151527 -0.489271777 -2.047136705 0.913924958 -0.093497319 1.113574274 5 -0.787685748 -0.264536021 0.149409794 -0.415727421 -1.031962353 1.936102593 6 0.138046923 Jul Aug Sep Oct Nov Dec 1 0.165718174 -0.247137480 0.692367552 1.613013917 1.214642040 -0.084118058 2 -0.332379313 -0.725607017 0.782616550 1.766823560 1.258174382 -1.171986337 3 0.299795560 1.141628737 0.691426217 1.082898142 0.366842157 -1.949936720 4 -0.297500737 0.466664132 0.318735557 1.235932665 -0.852612568 -1.221328860 5 1.211897977 1.121000746 1.047229991 1.858254181 -0.257076126 -0.494451744 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/1kgi81260028851.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/2lr6v1260028851.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/3rgoq1260028851.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/432t61260028851.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/5y3yp1260028851.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/6axzs1260028851.tab") > system("convert tmp/1kgi81260028851.ps tmp/1kgi81260028851.png") > system("convert tmp/2lr6v1260028851.ps tmp/2lr6v1260028851.png") > system("convert tmp/3rgoq1260028851.ps tmp/3rgoq1260028851.png") > system("convert tmp/432t61260028851.ps tmp/432t61260028851.png") > system("convert tmp/5y3yp1260028851.ps tmp/5y3yp1260028851.png") > > > proc.time() user system elapsed 1.280 0.792 1.795