R version 2.12.1 (2010-12-16) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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(89924,31795,27922,59954,52150,39964,34604,51106,52593,68794,47124,32315,42248,36088,52744,72586,92334,80761,71078,63713,57122,55243,62143,62708,62474,64250,71866,69886,58724,55298,52594,54854,54694,49298,44659,43657,47002,47042,48959,49750,54048,60067,68929,74617,75940,72762,75621,73008,74196,78878,83812,91624,89388,110410,113857,112060,117236,132810,137699,146409) > 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 137856037.5 156065.3 0.0 0.0 > m$fitted level slope sea Jan 1 89924.00 0.00000 0.0000 Feb 1 34816.60 -3042.91316 -3021.5994 Mar 1 30946.85 -3047.98628 -3024.8535 Apr 1 62842.33 -2795.70701 -2888.3251 May 1 55057.66 -2837.05008 -2907.6561 Jun 1 42907.41 -2923.99265 -2943.4055 Jul 1 37556.62 -2949.14224 -2952.6243 Aug 1 53985.85 -2728.83012 -2879.8469 Sep 1 55457.27 -2676.95340 -2864.2672 Oct 1 71589.43 -2426.64794 -2795.4286 Nov 1 49988.60 -2699.65456 -2864.6041 Dec 1 35222.48 -2882.36055 -2907.4778 Jan 2 17251.98 -2492.49143 24996.0237 Feb 2 37806.15 -1680.62812 -1718.1531 Mar 2 54433.22 -1386.71480 -1689.2218 Apr 2 74242.29 -1028.55841 -1656.2923 May 2 93958.63 -661.12463 -1624.6343 Jun 2 82401.95 -862.66219 -1640.9540 Jul 2 72731.89 -1032.23153 -1653.8921 Aug 2 65376.00 -1158.57511 -1662.9960 Sep 2 58792.64 -1270.74135 -1670.6443 Oct 2 56914.48 -1283.70663 -1671.4824 Nov 2 63803.45 -1104.03616 -1660.4537 Dec 2 64366.26 -1066.37307 -1658.2552 Jan 3 47635.72 -1186.59176 14838.2837 Feb 3 65311.07 -552.17647 -1061.0661 Mar 3 72920.59 -356.71121 -1054.5920 Apr 3 70941.85 -396.39970 -1055.8471 May 3 59787.96 -664.93068 -1063.9637 Jun 3 56363.99 -735.09736 -1065.9924 Jul 3 53661.40 -786.01031 -1067.4017 Aug 3 55919.28 -705.95726 -1065.2788 Sep 3 55758.91 -691.39006 -1064.9084 Oct 3 50366.01 -818.70308 -1068.0135 Nov 3 45729.47 -923.45695 -1070.4658 Dec 3 44727.51 -925.63729 -1070.5148 Jan 4 36304.81 -1080.67899 10697.1912 Feb 4 47873.87 -652.99373 -831.8717 Mar 4 49789.80 -579.52226 -830.8022 Apr 4 50580.25 -539.96559 -830.2483 May 4 54876.35 -399.07875 -828.3501 Jun 4 60892.91 -210.61314 -825.9058 Jul 4 69751.55 57.88500 -822.5528 Aug 4 75437.53 225.72140 -820.5341 Sep 4 76760.15 258.65496 -820.1524 Oct 4 73583.31 154.85019 -821.3116 Nov 4 76441.43 237.01600 -820.4272 Dec 4 73829.33 149.93607 -821.3308 Jan 5 66023.78 -63.73305 8172.2204 Feb 5 79502.83 394.96457 -624.8306 Mar 5 84435.83 535.57322 -623.8257 Apr 5 92246.26 761.94853 -622.2649 May 5 90010.89 668.30539 -622.8878 Jun 5 111028.79 1306.45476 -618.7909 Jul 5 114475.37 1373.80275 -618.3736 Aug 5 112678.97 1273.71257 -618.9723 Sep 5 117854.26 1397.27054 -618.2588 Oct 5 133425.75 1847.43395 -615.7493 Nov 5 138314.23 1944.27356 -615.2281 Dec 5 147023.11 2160.22502 -614.1056 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -2.817557992 -0.070482201 2.977002868 -0.424792518 -0.792539802 2 -1.541667227 1.662196392 1.549165388 1.792734558 1.753824360 -0.920742919 3 -1.451578952 1.446424609 0.687310022 -0.136555057 -0.905399526 -0.232153241 4 -0.671284463 0.995225184 0.215735123 0.115030648 0.406005964 0.538542726 5 -0.700133806 1.080665572 0.380555123 0.610017551 -0.251316111 1.706133960 Jul Aug Sep Oct Nov Dec 1 -0.206401190 1.647246810 0.356846619 1.597142145 -1.627290492 -1.023543922 2 -0.743972801 -0.533955005 -0.457878596 -0.051250343 0.689313290 0.140538923 3 -0.165509284 0.255997210 0.045874690 -0.395233236 -0.320883076 -0.006596015 4 0.761195398 0.472316004 0.092042155 -0.288245786 0.226786656 -0.238999048 5 0.179419492 -0.265761032 0.327054125 1.188108444 0.254893176 0.566966647 > 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/rcomp/tmp/1nv0u1322218760.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/rcomp/tmp/2gg951322218760.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/rcomp/tmp/3qeti1322218760.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/rcomp/tmp/4cnnk1322218760.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',type='b') > grid() > dev.off() null device 1 > postscript(file="/var/www/rcomp/tmp/524091322218760.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/6i7ft1322218760.tab") > > try(system("convert tmp/1nv0u1322218760.ps tmp/1nv0u1322218760.png",intern=TRUE)) character(0) > try(system("convert tmp/2gg951322218760.ps tmp/2gg951322218760.png",intern=TRUE)) character(0) > try(system("convert tmp/3qeti1322218760.ps tmp/3qeti1322218760.png",intern=TRUE)) character(0) > try(system("convert tmp/4cnnk1322218760.ps tmp/4cnnk1322218760.png",intern=TRUE)) character(0) > try(system("convert tmp/524091322218760.ps tmp/524091322218760.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.188 0.332 2.546