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(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 17719.925 0.000 34600.134 815.119 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.625 -4.58904263 -334.6262719 Mar 1 9148.208 -22.23036189 -61.4190689 Apr 1 9343.404 -8.81937153 396.7083849 May 1 9095.702 -19.39285047 -504.8744767 Jun 1 9271.864 -13.68363223 455.7850361 Jul 1 9434.408 -10.03992396 125.5140933 Aug 1 9698.138 -5.33632031 295.0604656 Sep 1 9674.896 -5.62859802 -237.5819788 Oct 1 9807.978 -3.35094275 227.5934810 Nov 1 9889.124 -1.94548511 27.3982919 Dec 1 9667.855 -5.60501880 -412.0214721 Jan 2 9552.041 -3.61040227 186.9432523 Feb 2 9409.530 -2.68610766 -372.0453594 Mar 2 9311.802 -4.17279617 -177.2831957 Apr 2 9154.666 -8.34497094 334.6470081 May 2 9193.558 -6.98469901 -494.2897316 Jun 2 9253.459 -5.31403602 372.4637230 Jul 2 9171.217 -6.85589376 -222.9456110 Aug 2 9097.557 -7.93365259 186.5607721 Sep 2 9105.021 -7.72594706 -276.2802348 Oct 2 9330.491 -5.08649186 612.5802750 Nov 2 9424.085 -4.23982837 202.2503908 Dec 2 9518.270 -3.77030289 -201.9348515 Jan 3 9459.653 -3.83436376 146.2801629 Feb 3 9243.964 -4.56730217 -600.3926398 Mar 3 9214.564 -4.81574726 -0.1574085 Apr 3 9224.469 -4.57916826 342.2964655 May 3 9155.452 -5.81034756 -607.4279227 Jun 3 8982.900 -8.99400505 204.7865499 Jul 3 9195.984 -5.18028867 270.3794909 Aug 3 9183.672 -5.28384289 -60.5541585 Sep 3 9382.239 -2.87829800 -107.6368611 Oct 3 9524.171 -1.54558637 643.4076574 Nov 3 9459.879 -1.97218082 -24.8282626 Dec 3 9554.286 -1.49714182 99.0984421 Jan 4 9426.439 -2.04658906 4.6767932 Feb 4 9353.996 -2.45262077 -613.8251408 Mar 4 9400.454 -2.02553276 150.7437093 Apr 4 9364.420 -2.42987756 323.1313324 May 4 9439.204 -1.34601322 -421.4508420 Jun 4 9526.269 -0.04784727 144.2969671 Jul 4 9354.300 -2.45049009 -145.4858465 Aug 4 9291.748 -3.19410451 -221.7563347 Sep 4 9527.342 -0.71983533 256.6965501 Oct 4 9624.511 0.09602796 685.8597871 Nov 4 9836.232 1.49713127 265.2409343 Dec 4 9824.505 1.42373578 38.7154690 Jan 5 9747.285 0.99755164 -89.9762518 Feb 5 9801.800 1.33567199 -507.6861450 Mar 5 9816.170 1.43942943 129.6152476 Apr 5 9676.923 0.06568997 26.3750423 May 5 9568.965 -1.13769999 -518.2052104 Jun 5 9686.240 0.25070343 501.8277898 Jul 5 9795.057 1.49455615 -90.8376573 Aug 5 9964.838 3.26280651 -202.6118803 Sep 5 9928.375 2.89852693 -34.7169593 Oct 5 9736.868 1.38879256 260.3583828 Nov 5 9841.258 2.06803444 590.0302674 Dec 5 9925.009 2.54992804 146.6340137 Jan 6 10052.859 3.28013451 57.0615418 Feb 6 9964.615 2.69702650 -697.1007935 Mar 6 9805.627 1.50975423 17.0356307 Apr 6 9863.866 1.98797486 232.2035348 May 6 9841.262 1.75920451 -725.8593462 Jun 6 9901.612 2.32917574 508.4272680 Jul 6 9895.063 2.24374033 -216.9173238 Aug 6 10131.574 4.36484885 272.5629338 Sep 6 10180.038 4.72613312 -27.7668886 Oct 6 10200.947 4.84359730 166.7843357 Nov 6 10156.284 4.52340528 425.5373128 Dec 6 10267.050 5.15994571 328.1879560 Jan 7 10388.663 5.84920868 289.4073554 Feb 7 10404.426 5.91101680 -666.5897164 Mar 7 10104.894 3.82775191 -543.8548044 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94825201 -1.35872117 1.35850437 -1.70098600 1.45024067 2 -0.86770754 -1.07943661 -0.67237091 -1.05021515 0.33348931 0.48748447 3 -0.41333638 -1.58347704 -0.18108752 0.10554256 -0.46338244 -1.21588208 4 -0.94425360 -0.52279134 0.35940018 -0.24786871 0.56252922 0.64800878 5 -0.58613071 0.39717523 0.09619240 -1.03376293 -0.79320034 0.87196827 6 0.93295767 -0.67960485 -1.19665516 0.41883112 -0.18146023 0.43296299 7 0.86687876 0.07366984 -2.26531949 Jul Aug Sep Oct Nov Dec 1 1.31890702 2.05127008 -0.13404660 1.03721999 0.63132643 -1.63810581 2 -0.57070790 -0.49925374 0.11535120 1.74711149 0.73883354 0.73783000 3 1.64089527 -0.05312115 1.52376780 1.08343245 -0.46947738 0.72121094 4 -1.27002893 -0.44669269 1.78077302 0.73102041 1.58093275 -0.09874863 5 0.80304840 1.24985211 -0.29584215 -1.44974138 0.76844300 0.60912792 6 -0.06577464 1.74014105 0.32821505 0.12058985 -0.36904207 0.79171879 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 > postscript(file="/var/www/html/rcomp/tmp/1hjgl1291896156.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') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/html/rcomp/tmp/2hjgl1291896156.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/html/rcomp/tmp/3saf61291896156.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/html/rcomp/tmp/421w91291896156.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/html/rcomp/tmp/5vsdu1291896156.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/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/6gtci1291896156.tab") > > try(system("convert tmp/1hjgl1291896156.ps tmp/1hjgl1291896156.png",intern=TRUE)) character(0) > try(system("convert tmp/2hjgl1291896156.ps tmp/2hjgl1291896156.png",intern=TRUE)) character(0) > try(system("convert tmp/3saf61291896156.ps tmp/3saf61291896156.png",intern=TRUE)) character(0) > try(system("convert tmp/421w91291896156.ps tmp/421w91291896156.png",intern=TRUE)) character(0) > try(system("convert tmp/5vsdu1291896156.ps tmp/5vsdu1291896156.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.607 0.869 4.260