R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i686-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) > par20 = '' > par19 = '' > par18 = '' > par17 = '' > par16 = '' > par15 = '' > par14 = '' > par13 = '' > par12 = '' > par11 = '' > par10 = '' > par9 = '' > par8 = '' > par7 = '' > par6 = '' > par5 = '' > par4 = '' > par3 = '' > par2 = '' > par1 = '12' > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 17712.7499 0.0000 34617.8219 806.4681 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.755 -4.57309005 -334.8096561 Mar 1 9148.237 -22.22283242 -61.4764080 Apr 1 9343.340 -8.81343040 396.8037191 May 1 9095.718 -19.38977054 -504.9320859 Jun 1 9271.809 -13.67890532 455.8762416 Jul 1 9434.342 -10.03351191 125.6135385 Aug 1 9698.053 -5.32867043 295.1963833 Sep 1 9674.894 -5.61978427 -237.5842080 Oct 1 9807.960 -3.34211104 227.6369264 Nov 1 9889.112 -1.93643007 27.4257593 Dec 1 9667.922 -5.59561438 -412.1294553 Jan 2 9552.097 -3.60291303 186.8673313 Feb 2 9409.605 -2.67799039 -372.1463240 Mar 2 9311.867 -4.16388138 -177.3639328 Apr 2 9154.705 -8.33649044 334.5839639 May 2 9193.583 -6.97662343 -494.3062243 Jun 2 9253.457 -5.30611167 372.4775320 Jul 2 9171.222 -6.84882861 -222.9641323 Aug 2 9097.544 -7.92757698 186.5626922 Sep 2 9104.999 -7.71999510 -276.2545180 Oct 2 9330.401 -5.08062386 612.7133194 Nov 2 9424.007 -4.23368923 202.3460114 Dec 2 9518.233 -3.76366411 -201.8808635 Jan 3 9459.679 -3.82804163 146.2424875 Feb 3 9244.056 -4.56140679 -600.5230688 Mar 3 9214.628 -4.81013763 -0.2251594 Apr 3 9224.517 -4.57394589 342.2510525 May 3 9155.494 -5.80533755 -607.4807652 Jun 3 8982.944 -8.98964453 204.7144422 Jul 3 9195.950 -5.17610190 270.4534013 Aug 3 9183.678 -5.27917748 -60.5614924 Sep 3 9382.197 -2.87322832 -107.5589979 Oct 3 9524.094 -1.54031565 643.5103384 Nov 3 9459.830 -1.96698428 -24.7905174 Dec 3 9554.211 -1.49171292 99.1913631 Jan 4 9426.440 -2.04138556 4.6521080 Feb 4 9354.045 -2.44735856 -613.8878195 Mar 4 9400.491 -2.02040768 150.7148452 Apr 4 9364.470 -2.42462996 323.0754178 May 4 9439.241 -1.34107661 -421.4740886 Jun 4 9526.313 -0.04280025 144.2681458 Jul 4 9354.371 -2.44554986 -145.5868988 Aug 4 9291.789 -3.18980848 -221.8070516 Sep 4 9527.269 -0.71585058 256.8129433 Oct 4 9624.420 0.10021147 685.9681053 Nov 4 9836.120 1.50198135 265.3900692 Dec 4 9824.442 1.42878011 38.7751318 Jan 5 9747.287 1.00264192 -89.9925766 Feb 5 9801.813 1.34093348 -507.6904521 Mar 5 9816.203 1.44481250 129.5848499 Apr 5 9677.004 0.07151484 26.2692339 May 5 9569.039 -1.13197880 -518.2971269 Jun 5 9686.254 0.25567406 501.8348458 Jul 5 9795.084 1.49973166 -90.8459645 Aug 5 9964.850 3.26812001 -202.5947024 Sep 5 9928.380 2.90361451 -34.7286819 Oct 5 9736.873 1.39319798 260.3191051 Nov 5 9841.169 2.07215315 590.1389971 Dec 5 9924.922 2.55432680 146.7344593 Jan 6 10052.799 3.28504224 57.1429537 Feb 6 9964.624 2.70214316 -697.1267158 Mar 6 9805.684 1.51499737 16.9504548 Apr 6 9863.913 1.99311293 232.1659131 May 6 9841.342 1.76460063 -725.9441148 Jun 6 9901.662 2.33422210 508.3884058 Jul 6 9895.096 2.24857312 -216.9513563 Aug 6 10131.518 4.36915097 272.6603758 Sep 6 10180.004 4.73067643 -27.7264851 Oct 6 10200.938 4.84832578 166.7958202 Nov 6 10156.267 4.52790590 425.5454959 Dec 6 10266.978 5.16440960 328.2783905 Jan 7 10388.576 5.85385733 289.5146811 Feb 7 10404.399 5.91603110 -666.5622412 Mar 7 10104.995 3.83324770 -544.0110753 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94834233 -1.35941203 1.35792559 -1.70068213 1.44995450 2 -0.86799729 -1.07958715 -0.67266837 -1.05065570 0.33337676 0.48731361 3 -0.41298127 -1.58335184 -0.18137442 0.10541058 -0.46355025 -1.21612600 4 -0.94390451 -0.52257533 0.35934589 -0.24786600 0.56250478 0.64814967 5 -0.58580422 0.39730441 0.09631312 -1.03365617 -0.79345717 0.87165892 6 0.93330476 -0.67926414 -1.19658315 0.41881039 -0.18129460 0.43278283 7 0.86690006 0.07409722 -2.26486245 Jul Aug Sep Oct Nov Dec 1 1.31903979 2.05149724 -0.13351080 1.03724835 0.63142850 -1.63791153 2 -0.57081806 -0.49954096 0.11525693 1.74690713 0.73903111 0.73823999 3 1.64058374 -0.05286422 1.52367746 1.08334431 -0.46940320 0.72111678 4 -1.27010286 -0.44703988 1.78023393 0.73099950 1.58105541 -0.09843705 5 0.80326015 1.24994536 -0.29599167 -1.45005712 0.76785277 0.60923914 6 -0.06594807 1.73978727 0.32841426 0.12075899 -0.36920605 0.79143511 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/fisher/rcomp/tmp/1naqq1385606018.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/fisher/rcomp/tmp/2z4hn1385606018.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/fisher/rcomp/tmp/3vime1385606018.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/fisher/rcomp/tmp/4p38s1385606018.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/fisher/rcomp/tmp/5litz1385606018.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/6nmj41385606018.tab") > > try(system("convert tmp/1naqq1385606018.ps tmp/1naqq1385606018.png",intern=TRUE)) character(0) > try(system("convert tmp/2z4hn1385606018.ps tmp/2z4hn1385606018.png",intern=TRUE)) character(0) > try(system("convert tmp/3vime1385606018.ps tmp/3vime1385606018.png",intern=TRUE)) character(0) > try(system("convert tmp/4p38s1385606018.ps tmp/4p38s1385606018.png",intern=TRUE)) character(0) > try(system("convert tmp/5litz1385606018.ps tmp/5litz1385606018.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 4.114 0.809 4.897