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(12008,9169,8788,8417,8247,8197,8236,8253,7733,8366,8626,8863,10102,8463,9114,8563,8872,8301,8301,8278,7736,7973,8268,9476,11100,8962,9173,8738,8459,8078,8411,8291,7810,8616,8312,9692,9911,8915,9452,9112,8472,8230,8384,8625,8221,8649,8625,10443,10357,8586,8892,8329,8101,7922,8120,7838,7735,8406,8209,9451,10041,9411,10405,8467,8464,8102,7627,7513,7510,8291,8064,9383) > 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 189898.6321 272.8315 10729.0724 0.0000 > m$fitted level slope sea Jan 1 12008.000 0.000000 0.000000 Feb 1 9388.892 -142.267204 -219.891874 Mar 1 8928.993 -148.056580 -140.993391 Apr 1 8577.330 -149.799459 -160.329701 May 1 8389.369 -150.171333 -142.369069 Jun 1 8329.194 -149.137135 -132.194399 Jul 1 8360.962 -146.813587 -124.961988 Aug 1 8383.297 -144.432430 -130.296519 Sep 1 7930.669 -149.140299 -197.668956 Oct 1 8407.379 -138.851646 -41.379120 Nov 1 8733.471 -130.685260 -107.471065 Dec 1 8974.048 -123.761886 -111.048035 Jan 2 8665.556 -120.111824 1436.444107 Feb 2 8627.572 -117.349492 -164.572147 Mar 2 9088.719 -100.673708 25.280996 Apr 2 8771.400 -105.519165 -208.400387 May 2 8931.367 -99.760154 -59.367410 Jun 2 8514.209 -106.981777 -213.208704 Jul 2 8407.990 -106.963617 -106.990359 Aug 2 8324.270 -106.390963 -46.270335 Sep 2 8018.044 -111.471056 -282.043788 Oct 2 8011.008 -108.728664 -38.007948 Nov 2 8322.920 -97.385585 -54.920372 Dec 2 9346.911 -69.120662 129.089417 Jan 3 9661.960 -62.854668 1438.040050 Feb 3 9318.009 -72.219673 -356.008535 Mar 3 9137.042 -75.945558 35.957713 Apr 3 9012.539 -77.456162 -274.539354 May 3 8558.465 -88.746526 -99.464906 Jun 3 8317.649 -93.346908 -239.648957 Jul 3 8463.139 -85.988483 -52.139262 Aug 3 8322.964 -87.684906 -31.964065 Sep 3 8119.738 -91.355182 -309.737921 Oct 3 8581.599 -73.541096 34.401171 Nov 3 8583.754 -71.102874 -271.753656 Dec 3 9426.619 -43.036138 265.381486 Jan 4 8653.418 -63.920586 1257.581790 Feb 4 9118.270 -45.455696 -203.269913 Mar 4 9358.736 -35.071407 93.264037 Apr 4 9320.671 -35.176218 -208.671137 May 4 8684.023 -55.707971 -212.022784 Jun 4 8505.369 -59.893582 -275.369393 Jul 4 8401.976 -61.384149 -17.976363 Aug 4 8557.483 -53.894729 67.516632 Sep 4 8600.159 -50.536599 -379.158758 Oct 4 8616.980 -48.184377 32.020094 Nov 4 8967.667 -34.343588 -342.667254 Dec 4 9844.356 -3.501413 598.644128 Jan 5 9384.051 -18.790692 972.948537 Feb 5 8948.493 -33.796988 -362.493167 Mar 5 8805.332 -37.837753 86.667648 Apr 5 8479.938 -48.322056 -150.938463 May 5 8301.402 -52.999666 -200.401503 Jun 5 8191.261 -55.042583 -269.261486 Jul 5 8162.713 -54.093813 -42.713025 Aug 5 7842.729 -63.648866 -4.729298 Sep 5 8046.132 -54.022144 -311.131609 Oct 5 8385.113 -39.848749 20.887131 Nov 5 8678.079 -27.916192 -469.079221 Dec 5 8757.616 -24.108324 693.383513 Jan 6 8938.586 -16.837846 1102.414353 Feb 6 9597.383 7.852230 -186.382694 Mar 6 10160.207 28.430079 244.793311 Apr 6 8935.016 -17.845154 -468.015659 May 6 8675.486 -26.695265 -211.486440 Jun 6 8408.071 -35.475794 -306.070597 Jul 6 7744.849 -58.370189 -117.848663 Aug 6 7563.586 -62.858407 -50.585523 Sep 6 7805.095 -51.729570 -295.094701 Oct 6 8218.663 -34.728579 72.336647 Nov 6 8492.003 -23.517407 -428.003166 Dec 6 8702.994 -15.028432 680.006275 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -3.663850783 -0.722911171 -0.467289400 -0.087466501 0.206089529 2 -0.490615788 0.165819693 1.289982463 -0.492000875 0.602130350 -0.719305696 3 0.917071342 -0.605346919 -0.240653723 -0.109458120 -0.848507729 -0.342359597 4 -1.684176764 1.163173173 0.632355164 -0.006718301 -1.350480001 -0.275819163 5 -1.038749331 -0.924423480 -0.242260994 -0.643842691 -0.291901610 -0.127977612 6 0.463214323 1.504043815 1.231665133 -2.803287523 -0.541349304 -0.538760070 Jul Aug Sep Oct Nov Dec 1 0.413929007 0.386718887 -0.704062662 1.428633472 1.060531318 0.846224807 2 0.001729066 0.052613135 -0.452096222 0.236179881 0.950954699 2.530410394 3 0.537526130 -0.121919195 -0.259917963 1.244420281 0.170138674 2.054590101 4 -0.097569314 0.486473644 0.216610906 0.151071278 0.893628843 2.044929586 5 0.059326270 -0.595439390 0.598105137 0.879948596 0.744474208 0.241013864 6 -1.404523863 -0.274983568 0.681087404 1.040761324 0.688663991 0.525732422 > 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/1znc71259944980.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/2zxz81259944980.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/35gaf1259944980.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/4rtpg1259944980.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/54u9p1259944980.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/6ve4z1259944980.tab") > > system("convert tmp/1znc71259944980.ps tmp/1znc71259944980.png") > system("convert tmp/2zxz81259944980.ps tmp/2zxz81259944980.png") > system("convert tmp/35gaf1259944980.ps tmp/35gaf1259944980.png") > system("convert tmp/4rtpg1259944980.ps tmp/4rtpg1259944980.png") > system("convert tmp/54u9p1259944980.ps tmp/54u9p1259944980.png") > > > proc.time() user system elapsed 1.583 0.836 2.038