R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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 17713.664 0.000 34615.617 807.299 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.738 -4.57503397 -334.7875196 Mar 1 9148.232 -22.22370942 -61.4690305 Apr 1 9343.348 -8.81409412 396.7931527 May 1 9095.716 -19.39005090 -504.9257038 Jun 1 9271.816 -13.67941076 455.8658809 Jul 1 9434.350 -10.03424633 125.6016347 Aug 1 9698.064 -5.32959008 295.1800125 Sep 1 9674.895 -5.62085046 -237.5847159 Oct 1 9807.963 -3.34319558 227.6316521 Nov 1 9889.114 -1.93755181 27.4222805 Dec 1 9667.914 -5.59677102 -412.1170850 Jan 2 9552.089 -3.60384830 186.8768243 Feb 2 9409.595 -2.67899180 -372.1339040 Mar 2 9311.859 -4.16496378 -177.3538882 Apr 2 9154.699 -8.33751561 334.5920685 May 2 9193.580 -6.97757668 -494.3039001 Jun 2 9253.457 -5.30704495 372.4759769 Jul 2 9171.221 -6.84965695 -222.9621120 Aug 2 9097.545 -7.92828444 186.5625588 Sep 2 9105.001 -7.72068337 -276.2576146 Oct 2 9330.412 -5.08131169 612.6975154 Nov 2 9424.017 -4.23441714 202.3342303 Dec 2 9518.239 -3.76445657 -201.8877552 Jan 3 9459.676 -3.82879669 146.2465359 Feb 3 9244.044 -4.56210827 -600.5076623 Mar 3 9214.619 -4.81080270 -0.2165604 Apr 3 9224.511 -4.57455968 342.2567606 May 3 9155.489 -5.80592525 -607.4743117 Jun 3 8982.938 -8.99015440 204.7233023 Jul 3 9195.954 -5.17658585 270.4452240 Aug 3 9183.677 -5.27971708 -60.5609705 Sep 3 9382.203 -2.87382164 -107.5684708 Oct 3 9524.104 -1.54094165 643.4976443 Nov 3 9459.836 -1.96760263 -24.7956415 Dec 3 9554.220 -1.49236275 99.1802315 Jan 4 9426.440 -2.04200723 4.6547065 Feb 4 9354.039 -2.44798556 -613.8803763 Mar 4 9400.487 -2.02101863 150.7185611 Apr 4 9364.464 -2.42525722 323.0822973 May 4 9439.237 -1.34166619 -421.4709486 Jun 4 9526.308 -0.04340074 144.2715824 Jul 4 9354.362 -2.44614244 -145.5749801 Aug 4 9291.784 -3.19031980 -221.8008299 Sep 4 9527.278 -0.71632542 256.7993941 Oct 4 9624.432 0.09970966 685.9547918 Nov 4 9836.135 1.50139419 265.3718236 Dec 4 9824.450 1.42816634 38.7673704 Jan 5 9747.286 1.00202237 -89.9909936 Feb 5 9801.812 1.34029273 -507.6898995 Mar 5 9816.199 1.44415573 129.5885070 Apr 5 9676.994 0.07080432 26.2818444 May 5 9569.029 -1.13267722 -518.2857008 Jun 5 9686.252 0.25506854 501.8347698 Jul 5 9795.081 1.49910586 -90.8449107 Aug 5 9964.849 3.26747802 -202.5970160 Sep 5 9928.379 2.90299521 -34.7277418 Oct 5 9736.872 1.39266123 260.3234979 Nov 5 9841.180 2.07165023 590.1264369 Dec 5 9924.933 2.55379055 146.7223251 Jan 6 10052.807 3.28444383 57.1328619 Feb 6 9964.623 2.70151742 -697.1240748 Mar 6 9805.676 1.51435467 16.9607326 Apr 6 9863.907 1.99248523 232.1709387 May 6 9841.333 1.76394392 -725.9339007 Jun 6 9901.656 2.33360760 508.3932940 Jul 6 9895.092 2.24798390 -216.9472566 Aug 6 10131.525 4.36862517 272.6490647 Sep 6 10180.009 4.73012028 -27.7317050 Oct 6 10200.940 4.84774720 166.7939296 Nov 6 10156.269 4.52735315 425.5441936 Dec 6 10266.987 5.16385994 328.2676204 Jan 7 10388.587 5.85328339 289.5017995 Feb 7 10404.403 5.91541211 -666.5658201 Mar 7 10104.982 3.83257015 -543.9927888 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94833334 -1.35932643 1.35800128 -1.70072006 1.44999185 2 -0.86796446 -1.07956785 -0.67263042 -1.05060412 0.33339663 0.48733594 3 -0.41302773 -1.58337128 -0.18133924 0.10542887 -0.46353080 -1.21609902 4 -0.94394921 -0.52260093 0.35935310 -0.24786787 0.56250889 0.64813494 5 -0.58584627 0.39728941 0.09629753 -1.03367140 -0.79342772 0.87169898 6 0.93326500 -0.67930953 -1.19659604 0.41881562 -0.18131315 0.43280512 7 0.86689732 0.07404471 -2.26492225 Jul Aug Sep Oct Nov Dec 1 1.31902541 2.05147093 -0.13357917 1.03724339 0.63141467 -1.63793928 2 -0.57080637 -0.49950793 0.11527125 1.74693469 0.73900553 0.73819043 3 1.64062703 -0.05289428 1.52369098 1.08335331 -0.46941634 0.72112740 4 -1.27009878 -0.44699684 1.78030332 0.73100211 1.58104226 -0.09847935 5 0.80323868 1.24993667 -0.29597802 -1.45002305 0.76792576 0.60922853 6 -0.06592761 1.73983213 0.32838972 0.12073877 -0.36918920 0.79146993 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/freestat/rcomp/tmp/1l0bo1291059515.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/freestat/rcomp/tmp/2l0bo1291059515.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/freestat/rcomp/tmp/3vrar1291059515.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/freestat/rcomp/tmp/4vrar1291059515.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/freestat/rcomp/tmp/5o0sc1291059515.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/6c2af1291059516.tab") > > try(system("convert tmp/1l0bo1291059515.ps tmp/1l0bo1291059515.png",intern=TRUE)) character(0) > try(system("convert tmp/2l0bo1291059515.ps tmp/2l0bo1291059515.png",intern=TRUE)) character(0) > try(system("convert tmp/3vrar1291059515.ps tmp/3vrar1291059515.png",intern=TRUE)) character(0) > try(system("convert tmp/4vrar1291059515.ps tmp/4vrar1291059515.png",intern=TRUE)) character(0) > try(system("convert tmp/5o0sc1291059515.ps tmp/5o0sc1291059515.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.156 1.162 2.394