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(7291,6820,8031,7862,7357,7213,7079,7012,7319,8148,7599,6908,7878,7407,7911,7323,7179,6758,6934,6696,7688,8296,7697,7907,7592,7710,9011,8225,7733,8062,7859,8221,8330,8868,9053,8811,8120,7953,8878,8601,8361,9116,9310,9891,10147,10317,10682,10276,10614,9413,11068,9772,10350,10541,10049,10714,10759,11684,11462,10485,11056,10184,11082,10554,11315,10847,11104,11026,11073,12073,12328,11172) > 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 35884.73559 60.98047 33842.96309 40587.55547 > m$fitted level slope sea Jan 1 7291.000 0.000000 0.00000 Feb 1 7048.684 -13.261165 -134.84061 Mar 1 7493.268 18.112521 390.70038 Apr 1 7713.329 28.806651 68.17160 May 1 7594.138 23.079846 -171.85157 Jun 1 7412.844 16.705993 -106.71810 Jul 1 7240.416 11.334001 -74.49667 Aug 1 7113.355 7.468740 -37.69585 Sep 1 7177.860 9.086174 114.94337 Oct 1 7600.439 21.100522 358.03195 Nov 1 7671.126 22.577434 -94.80645 Dec 1 7357.737 12.338774 -296.37782 Jan 2 7511.471 8.513654 298.78742 Feb 2 7659.670 11.847862 -311.03203 Mar 2 7671.234 11.835801 239.87142 Apr 2 7480.772 2.989743 -80.53806 May 2 7361.826 -1.857577 -133.96279 Jun 2 7105.353 -10.998255 -242.66877 Jul 2 7000.519 -14.155146 -27.57278 Aug 2 6900.204 -16.979717 -168.35807 Sep 2 7226.205 -5.761113 319.13191 Oct 2 7567.813 5.616955 583.83396 Nov 2 7675.534 8.886221 -20.92764 Dec 2 7951.361 16.589070 -155.28898 Jan 3 7772.575 12.670787 -97.83754 Feb 3 7835.491 14.187525 -146.08139 Mar 3 8181.058 26.900153 701.20063 Apr 3 8249.253 28.590285 -40.20751 May 3 8054.824 19.596558 -234.01274 Jun 3 8106.607 20.845241 -57.48651 Jul 3 8035.276 17.391586 -139.09793 Aug 3 8256.762 24.867681 -118.30560 Sep 3 8277.293 24.711006 54.46186 Oct 3 8326.535 25.584860 531.54670 Nov 3 8685.637 37.171792 232.67447 Dec 3 8842.150 41.145223 -79.37670 Jan 4 8678.049 34.578519 -474.56165 Feb 4 8507.431 27.421205 -471.79887 Mar 4 8350.259 20.390220 600.39917 Apr 4 8403.540 21.695479 184.61768 May 4 8495.062 24.479544 -161.51029 Jun 4 8780.554 34.753685 231.89056 Jul 4 9143.921 47.469968 34.95615 Aug 4 9572.856 62.007369 165.61551 Sep 4 9905.727 72.196286 132.97738 Oct 4 10016.546 73.630471 285.02841 Nov 4 10249.387 79.455713 369.09112 Dec 4 10306.331 78.643766 -21.34102 Jan 5 10610.389 86.755263 -86.66974 Feb 5 10373.687 74.779080 -831.62006 Mar 5 10428.399 74.010231 647.52610 Apr 5 10169.496 60.984670 -266.77891 May 5 10377.954 66.789656 -85.98597 Jun 5 10505.521 69.173204 11.43735 Jul 5 10455.983 64.554887 -359.84438 Aug 5 10576.959 66.730172 114.60682 Sep 5 10665.792 67.574870 84.42375 Oct 5 11037.520 79.101021 525.75122 Nov 5 11154.642 80.530856 292.27666 Dec 5 10971.811 70.675054 -382.24862 Jan 6 10973.862 68.102015 109.43140 Feb 6 10996.625 66.382446 -794.61989 Mar 6 10748.352 54.278014 457.96544 Apr 6 10796.133 54.025552 -239.57587 May 6 11072.330 62.693912 155.15242 Jun 6 11027.843 58.515667 -138.49816 Jul 6 11242.885 64.595027 -200.87952 Aug 6 11197.467 60.343122 -127.85992 Sep 6 11193.467 57.869122 -94.97906 Oct 6 11347.073 61.532450 688.04402 Nov 6 11610.551 69.230075 637.57413 Dec 6 11632.872 67.445216 -442.30300 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.043554792 1.835929028 0.950525853 -0.745883742 -1.050983666 2 0.832768071 0.699422629 -0.001309988 -0.963915942 -0.603078426 -1.281671879 3 -1.017989284 0.252633678 1.599720239 0.199606300 -1.096392312 0.160177160 4 -1.035592559 -1.024139899 -0.904397560 0.160550047 0.343326489 1.293210684 5 1.124749769 -1.607197490 -0.098833695 -1.634108677 0.726126703 0.300654684 6 -0.340815453 -0.224752029 -1.552745575 -0.031985041 1.095287425 -0.529881670 Jul Aug Sep Oct Nov Dec 1 -0.977067470 -0.714894154 0.294204356 2.129232932 0.254901741 -1.724262206 2 -0.475030151 -0.436603799 1.737029424 1.757298146 0.515692754 1.348386909 3 -0.461110356 1.022305834 -0.021714924 0.122682031 1.665571721 0.596870614 4 1.634663900 1.899317628 1.347628975 0.191897580 0.790404977 -0.111945867 5 -0.588897115 0.280032877 0.109597201 1.506237596 0.188225193 -1.305651695 6 0.775309708 -0.545048847 -0.318490297 0.473401534 0.998473944 -0.232186314 > 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/12xjt1259694170.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/2dpqd1259694170.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/3uivy1259694170.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/4psd31259694170.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/57s9s1259694170.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/6jqeb1259694170.tab") > > system("convert tmp/12xjt1259694170.ps tmp/12xjt1259694170.png") > system("convert tmp/2dpqd1259694170.ps tmp/2dpqd1259694170.png") > system("convert tmp/3uivy1259694170.ps tmp/3uivy1259694170.png") > system("convert tmp/4psd31259694170.ps tmp/4psd31259694170.png") > system("convert tmp/57s9s1259694170.ps tmp/57s9s1259694170.png") > > > proc.time() user system elapsed 1.649 0.828 2.000