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(58608,46865,51378,46235,47206,45382,41227,33795,31295,42625,33625,21538,56421,53152,53536,52408,41454,38271,35306,26414,31917,38030,27534,18387,50556,43901,48572,43899,37532,40357,35489,29027,34485,42598,30306,26451,47460,50104,61465,53726,39477,43895,31481,29896,33842,39120,33702,25094,51442,45594,52518,48564,41745,49585,32747,33379,35645,37034,35681,20972) > 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 25218066 0 5646430 0 > m$fitted level slope sea Jan 1 58608.00 0.0000 0.00000 Feb 1 48584.94 -509.6645 -1719.93799 Mar 1 50205.04 -434.3353 1172.95778 Apr 1 47804.35 -460.6843 -1569.34627 May 1 47353.72 -460.6185 -147.72400 Jun 1 46153.71 -465.1682 -771.71005 Jul 1 42789.19 -484.6541 -1562.18558 Aug 1 36325.67 -526.2441 -2530.67439 Sep 1 32328.96 -550.3762 -1033.95521 Oct 1 39422.96 -497.6860 3202.04159 Nov 1 36393.59 -514.9908 -2768.59473 Dec 1 25852.13 -583.0328 -4314.13028 Jan 2 38303.58 -1043.0550 18117.42192 Feb 2 50964.03 -902.7612 2187.97331 Mar 2 53164.64 -839.2545 371.36310 Apr 2 54346.71 -811.4051 -1938.70663 May 2 46231.12 -863.8639 -4777.12037 Jun 2 40168.55 -887.8654 -1897.54757 Jul 2 36039.85 -901.8264 -733.85077 Aug 2 30342.01 -924.0153 -3928.00577 Sep 2 32953.98 -906.3240 -1036.98006 Oct 2 33523.34 -898.6499 4506.66484 Nov 2 29375.99 -912.0751 -1841.98536 Dec 2 27329.27 -911.0815 -8942.26917 Jan 3 33416.86 -978.4488 17139.14458 Feb 3 40130.21 -948.6101 3770.78720 Mar 3 46427.82 -856.8207 2144.18109 Apr 3 44902.94 -864.7226 -1003.94166 May 3 41824.13 -882.0424 -4292.12794 Jun 3 40763.25 -882.9444 -406.25072 Jul 3 36525.76 -896.5568 -1036.75690 Aug 3 34374.63 -901.6715 -5347.62524 Sep 3 34710.77 -896.3249 -225.76565 Oct 3 35638.57 -888.9721 6959.43273 Nov 3 32973.41 -892.9971 -2667.40918 Dec 3 35902.76 -897.6798 -9451.75718 Jan 4 34480.45 -896.0157 12979.55293 Feb 4 43690.63 -866.6060 6413.37070 Mar 4 54438.14 -765.1858 7026.86443 Apr 4 54934.69 -752.9589 -1208.68545 May 4 47494.80 -804.0397 -8017.80469 Jun 4 43872.76 -819.3394 22.24271 Jul 4 35689.97 -850.4913 -4208.97378 Aug 4 34642.12 -851.2551 -4746.11714 Sep 4 33952.33 -850.6562 -110.32837 Oct 4 31894.90 -854.3343 7225.09763 Nov 4 34985.87 -848.6772 -1283.86784 Dec 4 35033.00 -849.1898 -9938.99681 Jan 5 40022.89 -853.4197 11419.11443 Feb 5 42548.60 -843.9604 3045.39814 Mar 5 44616.98 -824.5329 7901.02097 Apr 5 46420.91 -803.4584 2143.08839 May 5 47626.04 -789.1769 -5881.04458 Jun 5 47096.02 -787.7409 2488.97879 Jul 5 40183.90 -814.4969 -7436.89676 Aug 5 37971.10 -819.7105 -4592.09691 Sep 5 35860.13 -823.8842 -215.12837 Oct 5 32297.90 -830.4893 4736.10455 Nov 5 35139.63 -826.3126 541.37484 Dec 5 33796.45 -826.3332 -12824.45080 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.272042063 0.393769644 -0.392172175 0.002009666 -0.147558518 2 2.913306395 2.584201728 0.580228633 0.398484153 -1.456859116 -1.037122154 3 1.436712367 1.507629597 1.387513580 -0.131044411 -0.440454980 -0.035665379 4 -0.105505542 1.992744813 2.254813934 0.247508970 -1.327150212 -0.561605487 5 1.165925496 0.667326304 0.569605284 0.516459588 0.398109050 0.051605672 Jul Aug Sep Oct Nov Dec 1 -0.578387495 -1.192592715 -0.692251320 1.524821271 -0.504989342 -1.999940899 2 -0.646205167 -0.956001214 0.704852946 0.294250759 -0.647534710 -0.226266504 3 -0.668955563 -0.250135534 0.246795404 0.363642226 -0.353767297 0.762914087 4 -1.468496047 -0.039360568 0.032195824 -0.240480897 0.785652713 0.178674558 5 -1.221257103 -0.278881537 -0.257444136 -0.545538991 0.731221879 -0.102996599 > 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/16hwg1259943886.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/2z3231259943886.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/3j9lo1259943886.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/4vsxp1259943886.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/5n2nr1259943886.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/6cjv21259943886.tab") > system("convert tmp/16hwg1259943886.ps tmp/16hwg1259943886.png") > system("convert tmp/2z3231259943886.ps tmp/2z3231259943886.png") > system("convert tmp/3j9lo1259943886.ps tmp/3j9lo1259943886.png") > system("convert tmp/4vsxp1259943886.ps tmp/4vsxp1259943886.png") > system("convert tmp/5n2nr1259943886.ps tmp/5n2nr1259943886.png") > > > proc.time() user system elapsed 1.373 0.802 1.574