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(539,548,563,581,572,519,521,531,540,548,556,551,549,564,586,604,601,545,537,552,563,575,580,575,558,564,581,597,587,536,524,537,536,533,528,516,502,506,518,534,528,478,469,490,493,508,517,514,510,527,542,565,555,499,511,526,532,549,561,557,566,588,620,626,620,573,573,574,580,590) > par8 = 'FALSE' > par7 = '1' > par6 = '' > par5 = '1' > par4 = '' > par3 = '0' > par2 = 'periodic' > par1 = '12' > main = 'Seasonal Decomposition by Loess' > par1 <- as.numeric(par1) #seasonal period > if (par2 != 'periodic') par2 <- as.numeric(par2) #s.window > par3 <- as.numeric(par3) #s.degree > if (par4 == '') par4 <- NULL else par4 <- as.numeric(par4)#t.window > par5 <- as.numeric(par5)#t.degree > if (par6 != '') par6 <- as.numeric(par6)#l.window > par7 <- as.numeric(par7)#l.degree > if (par8 == 'FALSE') par8 <- FALSE else par9 <- TRUE #robust > nx <- length(x) > x <- ts(x,frequency=par1) > if (par6 != '') { + m <- stl(x,s.window=par2, s.degree=par3, t.window=par4, t.degre=par5, l.window=par6, l.degree=par7, robust=par8) + } else { + m <- stl(x,s.window=par2, s.degree=par3, t.window=par4, t.degre=par5, l.degree=par7, robust=par8) + } > m$time.series seasonal trend remainder Jan 1 -8.532077 541.6063 5.9257390 Feb 1 3.047411 542.6707 2.2818747 Mar 1 21.293565 543.7351 -2.0286563 Apr 1 36.809297 544.8746 -0.6838855 May 1 28.825027 546.0141 -2.8391117 Jun 1 -24.195797 547.2063 -4.0104788 Jul 1 -27.549959 548.3985 0.1514917 Aug 1 -16.093350 549.6569 -2.5635956 Sep 1 -11.470084 550.9154 0.5546596 Oct 1 -2.796746 552.6891 -1.8923604 Nov 1 3.688487 554.4628 -2.1512750 Dec 1 -3.025778 556.4991 -2.4732725 Jan 2 -8.532077 558.5353 -1.0032369 Feb 2 3.047411 560.4676 0.4849858 Mar 2 21.293565 562.3999 2.3065417 Apr 2 36.809297 564.3504 2.8402798 May 2 28.825027 566.3010 5.8740209 Jun 2 -24.195797 567.8105 1.3852976 Jul 2 -27.549959 569.3200 -4.7700881 Aug 2 -16.093350 569.7090 -1.6156981 Sep 2 -11.470084 570.0980 4.3720345 Oct 2 -2.796746 569.5405 8.2562015 Nov 2 3.688487 568.9830 7.3284740 Dec 2 -3.025778 567.8653 10.1604522 Jan 3 -8.532077 566.7476 -0.2155366 Feb 3 3.047411 564.8762 -3.9236136 Mar 3 21.293565 563.0048 -3.2983573 Apr 3 36.809297 559.8510 0.3397088 May 3 28.825027 556.6972 1.4777778 Jun 3 -24.195797 552.4660 7.7297988 Jul 3 -27.549959 548.2348 3.3151573 Aug 3 -16.093350 543.2352 9.8581472 Sep 3 -11.470084 538.2356 9.2344796 Oct 3 -2.796746 532.8035 2.9932670 Nov 3 3.688487 527.3714 -3.0598401 Dec 3 -3.025778 522.1984 -3.1726451 Jan 4 -8.532077 517.0255 -6.4934172 Feb 4 3.047411 513.0325 -10.0799499 Mar 4 21.293565 509.0396 -12.3331493 Apr 4 36.809297 506.9759 -9.7852066 May 4 28.825027 504.9122 -5.7372610 Jun 4 -24.195797 504.9586 -2.7628403 Jul 4 -27.549959 505.0050 -8.4550820 Aug 4 -16.093350 506.5890 -0.4956901 Sep 4 -11.470084 508.1730 -3.7029557 Oct 4 -2.796746 510.4848 0.3119420 Nov 4 3.688487 512.7966 0.5149450 Dec 4 -3.025778 515.3719 1.6538698 Jan 5 -8.532077 517.9472 0.5848275 Feb 5 3.047411 520.8784 3.0742356 Mar 5 21.293565 523.8095 -3.1030229 Apr 5 36.809297 527.1267 1.0640056 May 5 28.825027 530.4439 -4.2689629 Jun 5 -24.195797 534.4197 -11.2238752 Jul 5 -27.549959 538.3954 0.1545501 Aug 5 -16.093350 543.4960 -1.4026258 Sep 5 -11.470084 548.5965 -5.1264591 Oct 5 -2.796746 554.4041 -2.6073568 Nov 5 3.688487 560.2117 -2.9001491 Dec 5 -3.025778 565.8994 -5.8736584 Jan 6 -8.532077 571.5872 2.9448653 Feb 6 3.047411 575.4872 9.4653946 Mar 6 21.293565 579.3872 19.3192571 Apr 6 36.809297 583.0625 6.1281902 May 6 28.825027 586.7378 4.4371262 Jun 6 -24.195797 590.2836 6.9121520 Jul 6 -27.549959 593.8294 6.7205154 Aug 6 -16.093350 597.1331 -7.0397514 Sep 6 -11.470084 600.4368 -8.9666755 Oct 6 -2.796746 603.5136 -10.7168173 > m$win s t l 701 19 13 > m$deg s t l 0 1 1 > m$jump s t l 71 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/wessaorg/rcomp/tmp/1m6wl1324644955.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(m,main=main) > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/wessaorg/rcomp/tmp/2mdfm1324644955.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(as.numeric(m$time.series[,'trend']),na.action=na.pass,lag.max = mylagmax,main='Trend') > acf(as.numeric(m$time.series[,'seasonal']),na.action=na.pass,lag.max = mylagmax,main='Seasonal') > acf(as.numeric(m$time.series[,'remainder']),na.action=na.pass,lag.max = mylagmax,main='Remainder') > par(op) > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3h4ws1324644955.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(as.numeric(m$time.series[!is.na(m$time.series[,'trend']),'trend']),main='Trend') > spectrum(as.numeric(m$time.series[!is.na(m$time.series[,'seasonal']),'seasonal']),main='Seasonal') > spectrum(as.numeric(m$time.series[!is.na(m$time.series[,'remainder']),'remainder']),main='Remainder') > par(op) > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/4ldhr1324644955.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(as.numeric(m$time.series[!is.na(m$time.series[,'trend']),'trend']),main='Trend') > cpgram(as.numeric(m$time.series[!is.na(m$time.series[,'seasonal']),'seasonal']),main='Seasonal') > cpgram(as.numeric(m$time.series[!is.na(m$time.series[,'remainder']),'remainder']),main='Remainder') > 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,'Seasonal Decomposition by Loess - Parameters',4,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Component',header=TRUE) > a<-table.element(a,'Window',header=TRUE) > a<-table.element(a,'Degree',header=TRUE) > a<-table.element(a,'Jump',header=TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Seasonal',header=TRUE) > a<-table.element(a,m$win['s']) > a<-table.element(a,m$deg['s']) > a<-table.element(a,m$jump['s']) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Trend',header=TRUE) > a<-table.element(a,m$win['t']) > a<-table.element(a,m$deg['t']) > a<-table.element(a,m$jump['t']) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Low-pass',header=TRUE) > a<-table.element(a,m$win['l']) > a<-table.element(a,m$deg['l']) > a<-table.element(a,m$jump['l']) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/wessaorg/rcomp/tmp/5kxrm1324644955.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Seasonal Decomposition by Loess - Time Series Components',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,'Fitted',header=TRUE) > a<-table.element(a,'Seasonal',header=TRUE) > a<-table.element(a,'Trend',header=TRUE) > a<-table.element(a,'Remainder',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,x[i]+m$time.series[i,'remainder']) + a<-table.element(a,m$time.series[i,'seasonal']) + a<-table.element(a,m$time.series[i,'trend']) + a<-table.element(a,m$time.series[i,'remainder']) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/wessaorg/rcomp/tmp/6xs0j1324644955.tab") > > try(system("convert tmp/1m6wl1324644955.ps tmp/1m6wl1324644955.png",intern=TRUE)) character(0) > try(system("convert tmp/2mdfm1324644955.ps tmp/2mdfm1324644955.png",intern=TRUE)) character(0) > try(system("convert tmp/3h4ws1324644955.ps tmp/3h4ws1324644955.png",intern=TRUE)) character(0) > try(system("convert tmp/4ldhr1324644955.ps tmp/4ldhr1324644955.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.358 0.266 1.629