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(15579,16348,15928,16171,15937,15713,15594,15683,16438,17032,17696,17745,19394,20148,20108,18584,18441,18391,19178,18079,18483,19644,19195,19650,20830,23595,22937,21814,21928,21777,21383,21467,22052,22680,24320,24977,25204,25739,26434,27525,30695,32436,30160,30236,31293,31077,32226,33865,32810,32242,32700,32819,33947,34148,35261,39506,41591,39148,41216,40225,41126,42362,40740,40256,39804,41002,41702,42254,43605,43271) > 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 160.45195 15261.10 157.44352 Feb 1 701.33503 15464.81 181.85091 Mar 1 62.55186 15668.52 196.92455 Apr 1 -561.98896 15889.64 843.34855 May 1 -309.69542 16110.76 135.93820 Jun 1 -181.88617 16348.77 -453.88543 Jul 1 -538.07768 16586.79 -454.70828 Aug 1 -227.12593 16829.16 -919.02968 Sep 1 465.49237 17071.53 -1099.01762 Oct 1 -25.49122 17335.71 -278.22047 Nov 1 284.60683 17599.90 -188.50496 Dec 1 169.82852 17878.90 -303.73016 Jan 2 160.45195 18157.91 1075.64290 Feb 2 701.33503 18385.49 1061.17403 Mar 2 62.55186 18613.08 1432.37141 Apr 2 -561.98896 18764.21 381.77838 May 2 -309.69542 18915.34 -164.64902 Jun 2 -181.88617 19028.67 -455.78388 Jul 2 -538.07768 19142.00 574.08203 Aug 2 -227.12593 19330.04 -1023.91672 Sep 2 465.49237 19518.09 -1500.58201 Oct 2 -25.49122 19803.75 -134.25458 Nov 2 284.60683 20089.40 -1179.00879 Dec 2 169.82852 20407.18 -927.01078 Jan 3 160.45195 20724.96 -55.41450 Feb 3 701.33503 21022.02 1871.64848 Mar 3 62.55186 21319.07 1555.37770 Apr 3 -561.98896 21620.58 755.40696 May 3 -309.69542 21922.09 315.60185 Jun 3 -181.88617 22219.93 -261.04479 Jul 3 -538.07768 22517.77 -596.69067 Aug 3 -227.12593 22827.60 -1133.47628 Sep 3 465.49237 23137.44 -1550.92844 Oct 3 -25.49122 23640.25 -934.75848 Nov 3 284.60683 24143.06 -107.67015 Dec 3 169.82852 24896.20 -89.03052 Jan 4 160.45195 25649.34 -605.79262 Feb 4 701.33503 26463.35 -1425.68388 Mar 4 62.55186 27277.36 -905.90889 Apr 4 -561.98896 28029.22 57.77215 May 4 -309.69542 28781.08 2223.61883 Jun 4 -181.88617 29456.79 3161.09519 Jul 4 -538.07768 30132.51 565.57231 Aug 4 -227.12593 30668.33 -205.19941 Sep 4 465.49237 31204.15 -376.63767 Oct 4 -25.49122 31565.40 -462.90918 Nov 4 284.60683 31926.66 14.73768 Dec 4 169.82852 32273.96 1421.21045 Jan 5 160.45195 32621.27 28.28148 Feb 5 701.33503 33198.68 -1658.01120 Mar 5 62.55186 33776.09 -1138.63764 Apr 5 -561.98896 34499.38 -1118.39521 May 5 -309.69542 35222.68 -965.98715 Jun 5 -181.88617 35979.04 -1649.15767 Jul 5 -538.07768 36735.41 -936.32742 Aug 5 -227.12593 37491.94 2241.18972 Sep 5 465.49237 38248.47 2877.04030 Oct 5 -25.49122 38888.65 284.84444 Nov 5 284.60683 39528.83 1402.56694 Dec 5 169.82852 39978.53 76.64031 Jan 6 160.45195 40428.24 537.31194 Feb 6 701.33503 40707.82 952.84533 Mar 6 62.55186 40987.40 -309.95504 Apr 6 -561.98896 41266.34 -448.34805 May 6 -309.69542 41545.27 -1431.57541 Jun 6 -181.88617 41812.82 -628.92995 Jul 6 -538.07768 42080.36 159.71628 Aug 6 -227.12593 42348.30 132.82665 Sep 6 465.49237 42616.24 523.27047 Oct 6 -25.49122 42889.99 406.50512 > 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/1sppa1353682306.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/28lty1353682306.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/3blev1353682306.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/4xy3z1353682306.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/57jme1353682306.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/6ygnx1353682306.tab") > > try(system("convert tmp/1sppa1353682306.ps tmp/1sppa1353682306.png",intern=TRUE)) character(0) > try(system("convert tmp/28lty1353682306.ps tmp/28lty1353682306.png",intern=TRUE)) character(0) > try(system("convert tmp/3blev1353682306.ps tmp/3blev1353682306.png",intern=TRUE)) character(0) > try(system("convert tmp/4xy3z1353682306.ps tmp/4xy3z1353682306.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.734 0.439 3.154