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(3440,2678,2981,2260,2844,2546,2456,2295,2379,2479,2057,2280,2351,2276,2548,2311,2201,2725,2408,2139,1898,2539,2070,2063,2565,2442,2194,2798,2074,2628,2289,2154,2467,2137,1850,2075,1791,1755,2232,1952,1822,2522,2074,2366,2173,2094,1833,1858,2040,2133,2921,3252,3318,3554,2308,1621,1315,1501,1418,1657) > 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 90526.324 0.000 2152.175 14812.861 > m$fitted level slope sea Jan 1 3440.000 0.000000 0.0000000 Feb 1 2762.564 -38.429447 -46.9341401 Mar 1 2964.375 -34.499251 -15.2280897 Apr 1 2404.324 -38.129490 -72.4254750 May 1 2787.724 -35.639830 -1.4668169 Jun 1 2616.938 -36.431010 -52.4248638 Jul 1 2504.811 -36.872530 -38.4418938 Aug 1 2353.204 -37.538165 -42.4873016 Sep 1 2397.890 -37.063868 -30.1539361 Oct 1 2491.665 -36.313479 -30.5882479 Nov 1 2155.562 -38.023037 -57.4948203 Dec 1 2279.783 -37.103084 -22.0078043 Jan 2 2068.042 -28.511078 307.1086461 Feb 2 2284.244 -22.563764 -33.2335761 Mar 2 2475.324 -20.349011 45.1095490 Apr 2 2406.492 -20.540164 -89.0668256 May 2 2219.498 -21.035837 3.5891071 Jun 2 2647.626 -19.690708 17.7787827 Jul 2 2473.890 -20.157572 -45.4519325 Aug 2 2218.347 -20.871307 -48.1179165 Sep 2 1967.151 -21.569207 -38.5934376 Oct 2 2418.827 -20.119192 57.3875731 Nov 2 2184.379 -20.793463 -85.9432322 Dec 2 2075.520 -20.928724 -0.8312873 Jan 3 2251.961 -25.141107 286.4759373 Feb 3 2452.967 -21.913312 -36.8706840 Mar 3 2212.183 -23.799873 9.4941115 Apr 3 2724.158 -21.945959 4.2837065 May 3 2250.973 -22.949644 -118.2796460 Jun 3 2514.061 -22.326397 76.7358999 Jul 3 2360.766 -22.620270 -54.7334896 Aug 3 2208.933 -22.913192 -38.1297624 Sep 3 2470.624 -22.261697 -40.6361620 Oct 3 2151.322 -22.965928 24.3069111 Nov 3 1959.959 -23.356765 -88.0568967 Dec 3 2054.974 -23.369313 4.6245173 Jan 4 1674.146 -18.980097 164.0739301 Feb 4 1735.291 -18.228122 10.1131397 Mar 4 2156.698 -14.919468 20.2852888 Apr 4 1958.340 -15.548500 17.1849692 May 4 1956.504 -15.521272 -136.2667937 Jun 4 2320.047 -14.833819 153.2039444 Jul 4 2161.793 -15.103481 -69.3515063 Aug 4 2367.167 -14.677774 -29.5154432 Sep 4 2220.445 -14.940134 -30.4669249 Oct 4 2076.113 -15.208346 34.5248389 Nov 4 1951.513 -15.402977 -104.4445125 Dec 4 1836.780 -15.316432 33.9993629 Jan 5 1882.469 -15.790876 149.5759522 Feb 5 2115.767 -14.145975 -13.2382629 Mar 5 2713.746 -10.085442 130.9808431 Apr 5 3131.834 -8.600924 65.6555750 May 5 3426.501 -8.012549 -147.2469228 Jun 5 3399.449 -8.043553 156.9843132 Jul 5 2619.234 -9.333066 -212.5700924 Aug 5 1824.255 -10.705420 -102.8790547 Sep 5 1405.432 -11.453109 -38.2909543 Oct 5 1417.355 -11.409290 80.6581114 Nov 5 1477.440 -11.312818 -68.5766198 Dec 5 1592.242 -11.461681 48.6323267 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.39131405 0.78945166 -1.74972541 1.40372920 -0.45003948 2 -0.70171830 0.69841380 0.70214853 -0.16127212 -0.55378764 1.49426312 3 0.71548278 0.69360077 -0.71832699 1.78145326 -1.50093781 0.95134608 4 -1.24768538 0.25370845 1.44190298 -0.60965079 0.04560536 1.26064780 5 0.20890095 0.80218318 2.00792148 1.42225303 1.00862759 -0.06331902 Jul Aug Sep Oct Nov Dec 1 -0.25206377 -0.38205563 0.27379799 0.43567217 -0.99824293 0.54023511 2 -0.51246296 -0.78305053 -0.76621716 1.57439874 -0.71311953 -0.29276113 3 -0.43557530 -0.42973254 0.94655020 -0.98800646 -0.56016101 0.39346450 4 -0.47693577 0.73317685 -0.43912124 -0.43035973 -0.36380151 -0.33031293 5 -2.56771568 -2.61247733 -1.35719268 0.07774772 0.23771639 0.41958003 > 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/19e7c1322577468.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/2l4ao1322577468.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/3ckyj1322577468.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/4zq481322577468.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/5o0671322577468.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/640wi1322577468.tab") > > try(system("convert tmp/19e7c1322577468.ps tmp/19e7c1322577468.png",intern=TRUE)) character(0) > try(system("convert tmp/2l4ao1322577468.ps tmp/2l4ao1322577468.png",intern=TRUE)) character(0) > try(system("convert tmp/3ckyj1322577468.ps tmp/3ckyj1322577468.png",intern=TRUE)) character(0) > try(system("convert tmp/4zq481322577468.ps tmp/4zq481322577468.png",intern=TRUE)) character(0) > try(system("convert tmp/5o0671322577468.ps tmp/5o0671322577468.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.793 0.258 2.075