R version 2.12.0 (2010-10-15) 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(9700,9081,9084,9743,8587,9731,9563,9998,9437,10038,9918,9252,9737,9035,9133,9487,8700,9627,8947,9283,8829,9947,9628,9318,9605,8640,9214,9567,8547,9185,9470,9123,9278,10170,9434,9655,9429,8739,9552,9687,9019,9672,9206,9069,9788,10312,10105,9863,9656,9295,9946,9701,9049,10190,9706,9765,9893,9994,10433,10073,10112,9266,9820,10097,9115,10411,9678,10408,10153,10368,10581,10597,10680,9738,9556) > 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 17712.8169 0.0000 34617.7004 806.6183 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.755 -4.57323719 -334.8078991 Mar 1 9148.237 -22.22291463 -61.4760003 Apr 1 9343.341 -8.81350400 396.8025128 May 1 9095.718 -19.38982861 -504.9313567 Jun 1 9271.810 -13.67897096 455.8751656 Jul 1 9434.342 -10.03358572 125.6125544 Aug 1 9698.053 -5.32874492 295.1950700 Sep 1 9674.893 -5.61986802 -237.5839450 Oct 1 9807.960 -3.34218992 227.6365197 Nov 1 9889.112 -1.93650787 27.4255474 Dec 1 9667.922 -5.59569782 -412.1282138 Jan 2 9552.096 -3.60297503 186.8679564 Feb 2 9409.604 -2.67806138 -372.1454296 Mar 2 9311.867 -4.16396482 -177.3632522 Apr 2 9154.705 -8.33657129 334.5844148 May 2 9193.583 -6.97670755 -494.3061702 Jun 2 9253.457 -5.30619486 372.4773654 Jul 2 9171.222 -6.84890395 -222.9638873 Aug 2 9097.544 -7.92764327 186.5626435 Sep 2 9104.999 -7.72006153 -276.2547673 Oct 2 9330.401 -5.08068647 612.7119651 Nov 2 9424.008 -4.23375213 202.3451698 Dec 2 9518.233 -3.76373048 -201.8812644 Jan 3 9459.679 -3.82810473 146.2430037 Feb 3 9244.055 -4.56146645 -600.5217144 Mar 3 9214.627 -4.81019516 -0.2246411 Apr 3 9224.517 -4.57400113 342.2514175 May 3 9155.494 -5.80539090 -607.4802818 Jun 3 8982.944 -8.98969150 204.7150848 Jul 3 9195.950 -5.17614872 270.4524616 Aug 3 9183.678 -5.27922955 -60.5613090 Sep 3 9382.197 -2.87328272 -107.5597204 Oct 3 9524.095 -1.54036934 643.5094490 Nov 3 9459.830 -1.96703686 -24.7907015 Dec 3 9554.211 -1.49176665 99.1904466 Jan 4 9426.440 -2.04143744 4.6524625 Feb 4 9354.045 -2.44741156 -613.8871818 Mar 4 9400.491 -2.02045921 150.7150492 Apr 4 9364.469 -2.42468219 323.0759130 May 4 9439.241 -1.34112591 -421.4739698 Jun 4 9526.312 -0.04285135 144.2684353 Jul 4 9354.370 -2.44559839 -145.5858431 Aug 4 9291.789 -3.18985193 -221.8065965 Sep 4 9527.269 -0.71589069 256.8116732 Oct 4 9624.421 0.10017039 685.9671427 Nov 4 9836.121 1.50193541 265.3887170 Dec 4 9824.443 1.42873329 38.7747351 Jan 5 9747.286 1.00259474 -89.9922998 Feb 5 9801.813 1.34088482 -507.6904216 Mar 5 9816.203 1.44476309 129.5851444 Apr 5 9677.003 0.07146133 26.2702973 May 5 9569.038 -1.13203120 -518.2963496 Jun 5 9686.254 0.25562814 501.8345364 Jul 5 9795.084 1.49968240 -90.8459013 Aug 5 9964.850 3.26806929 -202.5947919 Sep 5 9928.379 2.90356738 -34.7284203 Oct 5 9736.874 1.39315731 260.3195914 Nov 5 9841.170 2.07211550 590.1377810 Dec 5 9924.923 2.55428631 146.7335004 Jan 6 10052.800 3.28499710 57.1422559 Feb 6 9964.623 2.70209669 -697.1263161 Mar 6 9805.684 1.51495021 16.9512712 Apr 6 9863.913 1.99306605 232.1661217 May 6 9841.342 1.76455060 -725.9432978 Jun 6 9901.661 2.33417543 508.3887166 Jul 6 9895.096 2.24852865 -216.9510276 Aug 6 10131.519 4.36911179 272.6593023 Sep 6 10180.005 4.73063532 -27.7267648 Oct 6 10200.938 4.84828293 166.7958682 Nov 6 10156.268 4.52786575 425.5455170 Dec 6 10266.979 5.16437023 328.2774800 Jan 7 10388.577 5.85381682 289.5136305 Feb 7 10404.399 5.91598737 -666.5624233 Mar 7 10104.994 3.83319950 -544.0093998 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94834014 -1.35940558 1.35792881 -1.70068401 1.44995581 2 -0.86799351 -1.07958565 -0.67266591 -1.05065058 0.33337576 0.48731453 3 -0.41298340 -1.58335109 -0.18137177 0.10541107 -0.46354821 -1.21612243 4 -0.94390669 -0.52257735 0.35934606 -0.24786545 0.56250444 0.64814736 5 -0.58580630 0.39730284 0.09631224 -1.03365607 -0.79345392 0.87166071 6 0.93330042 -0.67926581 -1.19658207 0.41880957 -0.18129657 0.43278422 7 0.86689958 0.07409339 -2.26486453 Jul Aug Sep Oct Nov Dec 1 1.31903751 2.05149394 -0.13351473 1.03724817 0.63142775 -1.63791143 2 -0.57081628 -0.49953750 0.11525687 1.74690753 0.73902948 0.73823580 3 1.64058433 -0.05286695 1.52367693 1.08334526 -0.46940242 0.72111764 4 -1.27010016 -0.44703674 1.78023708 0.73099940 1.58105312 -0.09843854 5 0.80325655 1.24994319 -0.29598875 -1.45005228 0.76785759 0.60923696 6 -0.06594624 1.73978933 0.32841238 0.12075726 -0.36920343 0.79143735 7 > 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/11vgc1322317754.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/2pxuq1322317754.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/3f46n1322317754.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/4yaxt1322317754.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/5fcpj1322317754.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/6zczw1322317754.tab") > > try(system("convert tmp/11vgc1322317754.ps tmp/11vgc1322317754.png",intern=TRUE)) character(0) > try(system("convert tmp/2pxuq1322317754.ps tmp/2pxuq1322317754.png",intern=TRUE)) character(0) > try(system("convert tmp/3f46n1322317754.ps tmp/3f46n1322317754.png",intern=TRUE)) character(0) > try(system("convert tmp/4yaxt1322317754.ps tmp/4yaxt1322317754.png",intern=TRUE)) character(0) > try(system("convert tmp/5fcpj1322317754.ps tmp/5fcpj1322317754.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.060 0.120 2.188