R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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(37,30,47,35,30,43,82,40,47,19,52,136,80,42,54,66,81,63,137,72,107,58,36,52,79,77,54,84,48,96,83,66,61,53,30,74,69,59,42,65,70,100,63,105,82,81,75,102,121,98,76,77,63,37,35,23,40,29,37,51,20,28,13,22,25,13,16,13,16,17,9,17,25,14,8,7,10,7,10,3) > par8 = 'TRUE' > 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 9.6800417 27.233988 0.08597015 Feb 1 -0.3438541 29.739676 0.60417828 Mar 1 -9.7456700 32.245364 24.50030639 Apr 1 2.9160997 34.768022 -2.68412159 May 1 1.9021097 37.290680 -9.19278993 Jun 1 -5.3175616 39.816646 8.50091585 Jul 1 -3.8556116 42.342611 43.51300047 Aug 1 -4.2221224 44.878224 -0.65610205 Sep 1 0.1629425 47.413838 -0.57678018 Oct 1 -6.0556588 50.050920 -24.99526137 Nov 1 -1.0385484 52.688003 0.35054592 Dec 1 15.9178320 55.363017 64.71915072 Jan 2 9.6800417 58.038032 12.28192627 Feb 2 -0.3438541 60.548931 -18.20507693 Mar 2 -9.7456700 63.059830 0.68583984 Apr 2 2.9160997 64.929413 -1.84551293 May 2 1.9021097 66.798996 12.29889394 Jun 2 -5.3175616 67.671264 0.64629780 Jul 2 -3.8556116 68.543531 72.31208050 Aug 2 -4.2221224 69.164833 7.05728952 Sep 2 0.1629425 69.786135 37.05092293 Oct 2 -6.0556588 70.213540 -6.15788144 Nov 2 -1.0385484 70.640946 -33.60239733 Dec 2 15.9178320 70.666411 -34.58424290 Jan 3 9.6800417 70.691876 -1.37191772 Feb 3 -0.3438541 70.467983 6.87587109 Mar 3 -9.7456700 70.244090 -6.49842011 Apr 3 2.9160997 69.673714 11.41018625 May 3 1.9021097 69.103338 -23.00544775 Jun 3 -5.3175616 68.010979 33.30658243 Jul 3 -3.8556116 66.918620 19.93699146 Aug 3 -4.2221224 65.772156 4.44996662 Sep 3 0.1629425 64.625691 -3.78863384 Oct 3 -6.0556588 63.546931 -4.49127269 Nov 3 -1.0385484 62.468171 -31.42962304 Dec 3 15.9178320 61.868351 -3.78618272 Jan 4 9.6800417 61.268530 -1.94857163 Feb 4 -0.3438541 61.986288 -2.64243425 Mar 4 -9.7456700 62.704047 -10.95837688 Apr 4 2.9160997 64.613473 -2.52957282 May 4 1.9021097 66.522899 1.57499086 Jun 4 -5.3175616 68.989441 36.32812028 Jul 4 -3.8556116 71.455983 -4.60037146 Aug 4 -4.2221224 73.719351 35.50277098 Sep 4 0.1629425 75.982720 5.85433780 Oct 4 -6.0556588 76.127672 10.92798694 Nov 4 -1.0385484 76.272624 -0.23407544 Dec 4 15.9178320 74.372753 11.70941490 Jan 5 9.6800417 72.472882 38.84707599 Feb 5 -0.3438541 68.990162 29.35369171 Mar 5 -9.7456700 65.507443 20.23822741 Apr 5 2.9160997 61.631653 12.45224736 May 5 1.9021097 57.755863 3.34202694 Jun 5 -5.3175616 53.446672 -11.12911091 Jul 5 -3.8556116 49.137482 -10.28186991 Aug 5 -4.2221224 44.805136 -17.58301331 Sep 5 0.1629425 40.472790 -0.63573233 Oct 5 -6.0556588 37.028872 -1.97321360 Nov 5 -1.0385484 33.584955 4.45359362 Dec 5 15.9178320 31.184379 3.89778861 Jan 6 9.6800417 28.783804 -18.46384565 Feb 6 -0.3438541 27.005811 1.33804287 Mar 6 -9.7456700 25.227819 -2.48214864 Apr 6 2.9160997 23.476928 -4.39302724 May 6 1.9021097 21.726037 1.37185378 Jun 6 -5.3175616 20.200032 -1.88247035 Jul 6 -3.8556116 18.674027 1.18158436 Aug 6 -4.2221224 17.632630 -0.41050733 Sep 6 0.1629425 16.591232 -0.75417464 Oct 6 -6.0556588 15.667595 7.38806362 Nov 6 -1.0385484 14.743958 -4.70540964 Dec 6 15.9178320 14.048869 -12.96670101 Jan 7 9.6800417 13.353780 1.96617838 Feb 7 -0.3438541 12.669451 1.67440293 Mar 7 -9.7456700 11.985123 5.76054746 Apr 7 2.9160997 11.323582 -7.23968188 May 7 1.9021097 10.662042 -2.56415159 Jun 7 -5.3175616 10.027909 2.28965206 Jul 7 -3.8556116 9.393777 4.46183456 Aug 7 -4.2221224 8.786751 -1.56462868 > m$win s t l 801 19 13 > m$deg s t l 0 1 1 > m$jump s t l 81 2 2 > m$inner [1] 1 > m$outer [1] 15 > postscript(file="/var/fisher/rcomp/tmp/17qwd1353683906.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/fisher/rcomp/tmp/2juhj1353683906.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/fisher/rcomp/tmp/3m6891353683906.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/fisher/rcomp/tmp/46foj1353683906.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/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,'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/fisher/rcomp/tmp/5wpd51353683906.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/fisher/rcomp/tmp/61eke1353683906.tab") > > try(system("convert tmp/17qwd1353683906.ps tmp/17qwd1353683906.png",intern=TRUE)) character(0) > try(system("convert tmp/2juhj1353683906.ps tmp/2juhj1353683906.png",intern=TRUE)) character(0) > try(system("convert tmp/3m6891353683906.ps tmp/3m6891353683906.png",intern=TRUE)) character(0) > try(system("convert tmp/46foj1353683906.ps tmp/46foj1353683906.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.502 0.582 3.071