R version 2.13.0 (2011-04-13) Copyright (C) 2011 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 17713.0690 0.0000 34617.1005 806.4073 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.749 -4.57366027 -334.8034490 Mar 1 9148.234 -22.22303413 -61.4737147 Apr 1 9343.342 -8.81354309 396.8020472 May 1 9095.717 -19.38972344 -504.9310781 Jun 1 9271.812 -13.67895746 455.8742269 Jul 1 9434.345 -10.03366337 125.6103384 Aug 1 9698.058 -5.32892305 295.1918418 Sep 1 9674.895 -5.62008630 -237.5854155 Oct 1 9807.962 -3.34244047 227.6354054 Nov 1 9889.113 -1.93678404 27.4245585 Dec 1 9667.919 -5.59596809 -412.1268983 Jan 2 9552.093 -3.60321825 186.8704441 Feb 2 9409.601 -2.67830098 -372.1425205 Mar 2 9311.864 -4.16419357 -177.3607064 Apr 2 9154.702 -8.33677975 334.5869121 May 2 9193.582 -6.97686102 -494.3051148 Jun 2 9253.457 -5.30634136 372.4772478 Jul 2 9171.222 -6.84902720 -222.9638628 Aug 2 9097.543 -7.92773976 186.5627888 Sep 2 9104.999 -7.72014666 -276.2553742 Oct 2 9330.405 -5.08078862 612.7092385 Nov 2 9424.011 -4.23387533 202.3423797 Dec 2 9518.236 -3.76387412 -201.8833090 Jan 3 9459.679 -3.82824221 146.2429385 Feb 3 9244.052 -4.56158952 -600.5192107 Mar 3 9214.625 -4.81030747 -0.2221712 Apr 3 9224.515 -4.57409511 342.2529625 May 3 9155.492 -5.80547880 -607.4788543 Jun 3 8982.941 -8.98976467 204.7171396 Jul 3 9195.951 -5.17620808 270.4521289 Aug 3 9183.678 -5.27929635 -60.5618446 Sep 3 9382.200 -2.87336903 -107.5617572 Oct 3 9524.098 -1.54047711 643.5063644 Nov 3 9459.832 -1.96714539 -24.7927357 Dec 3 9554.214 -1.49188759 99.1883573 Jan 4 9426.440 -2.04155098 4.6523159 Feb 4 9354.043 -2.44752304 -613.8858972 Mar 4 9400.490 -2.02056833 150.7162064 Apr 4 9364.468 -2.42479728 323.0775255 May 4 9439.240 -1.34123315 -421.4727502 Jun 4 9526.312 -0.04295657 144.2690507 Jul 4 9354.368 -2.44571076 -145.5839419 Aug 4 9291.787 -3.18994179 -221.8051712 Sep 4 9527.272 -0.71597523 256.8098279 Oct 4 9624.424 0.10007515 685.9640416 Nov 4 9836.126 1.50181437 265.3846159 Dec 4 9824.445 1.42860135 38.7721854 Jan 5 9747.287 1.00246122 -89.9926634 Feb 5 9801.813 1.34074611 -507.6902501 Mar 5 9816.202 1.44461879 129.5858603 Apr 5 9677.000 0.07130546 26.2725502 May 5 9569.035 -1.13218532 -518.2934679 Jun 5 9686.253 0.25549674 501.8358897 Jul 5 9795.084 1.49955475 -90.8455975 Aug 5 9964.851 3.26793917 -202.5956815 Sep 5 9928.380 2.90343374 -34.7290747 Oct 5 9736.872 1.39304067 260.3198540 Nov 5 9841.172 2.07200420 590.1362449 Dec 5 9924.926 2.55416915 146.7310422 Jan 6 10052.803 3.28486616 57.1397617 Feb 6 9964.623 2.70195652 -697.1266429 Mar 6 9805.681 1.51480364 16.9533344 Apr 6 9863.911 1.99292660 232.1679903 May 6 9841.339 1.76440911 -725.9412728 Jun 6 9901.660 2.33404241 508.3900473 Jul 6 9895.095 2.24839983 -216.9502183 Aug 6 10131.521 4.36899457 272.6578212 Sep 6 10180.006 4.73050982 -27.7284212 Oct 6 10200.940 4.84815277 166.7945951 Nov 6 10156.268 4.52773785 425.5446938 Dec 6 10266.981 5.16424053 328.2755663 Jan 7 10388.580 5.85367872 289.5111680 Feb 7 10404.400 5.91583858 -666.5636030 Mar 7 10104.991 3.83303815 -544.0067345 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94834264 -1.35938515 1.35795496 -1.70069468 1.44996901 2 -0.86799113 -1.07958031 -0.67265498 -1.05064355 0.33339100 0.48732227 3 -0.41299927 -1.58336341 -0.18136382 0.10541896 -0.46354597 -1.21612197 4 -0.94392068 -0.52258191 0.35934888 -0.24786867 0.56250760 0.64814874 5 -0.58581971 0.39730107 0.09630737 -1.03366377 -0.79345113 0.87167398 6 0.93329663 -0.67928281 -1.19659267 0.41881565 -0.18129783 0.43278988 7 0.86689913 0.07408123 -2.26488584 Jul Aug Sep Oct Nov Dec 1 1.31903816 2.05149161 -0.13353506 1.03724496 0.63142262 -1.63792541 2 -0.57081705 -0.49953401 0.11526499 1.74691902 0.73902164 0.73822607 3 1.64060386 -0.05287137 1.52368500 1.08334469 -0.46941270 0.72111880 4 -1.27010851 -0.44702631 1.78025952 0.73100041 1.58105416 -0.09845537 5 0.80325983 1.24994659 -0.29599391 -1.45005328 0.76787572 0.60924009 6 -0.06594299 1.73980300 0.32840673 0.12075362 -0.36920539 0.79144574 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/wessaorg/rcomp/tmp/1y1ud1322214735.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/wessaorg/rcomp/tmp/24o1p1322214735.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/wessaorg/rcomp/tmp/330hy1322214735.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/wessaorg/rcomp/tmp/41bgq1322214736.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/wessaorg/rcomp/tmp/5hf3y1322214736.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/6stxm1322214736.tab") > > try(system("convert tmp/1y1ud1322214735.ps tmp/1y1ud1322214735.png",intern=TRUE)) character(0) > try(system("convert tmp/24o1p1322214735.ps tmp/24o1p1322214735.png",intern=TRUE)) character(0) > try(system("convert tmp/330hy1322214735.ps tmp/330hy1322214735.png",intern=TRUE)) character(0) > try(system("convert tmp/41bgq1322214736.ps tmp/41bgq1322214736.png",intern=TRUE)) character(0) > try(system("convert tmp/5hf3y1322214736.ps tmp/5hf3y1322214736.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.802 0.237 2.047