R version 2.12.1 (2010-12-16) Copyright (C) 2010 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(20995,17382,9367,31124,26551,30651,25859,25100,25778,20418,18688,20424,24776,19814,12738,31566,30111,30019,31934,25826,26835,20205,17789,20520,22518,15572,11509,25447,24090,27786,26195,20516,22759,19028,16971,20036,22485,18730,14538,27561,25985,34670,32066,27186,29586,21359,21553,19573,24256,22380,16167,27297,28287,33474,28229,28785,25597,18130,20198,22849,23118) > par8 = 'FALSE' > par7 = '1' > par6 = '' > par5 = '1' > par4 = '' > par3 = '0' > par2 = 'periodic' > par1 = '12' > main = 'Seasonal Decomposition by Loess' > #'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) #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 -295.1638 21633.47 -343.30570 Feb 1 -4491.9076 21834.85 39.06204 Mar 1 -10466.5723 22036.22 -2202.64937 Apr 1 5221.3996 22240.84 3661.76527 May 1 3579.9725 22445.45 525.57898 Jun 1 7856.9420 22658.85 135.20830 Jul 1 5355.3138 22872.25 -2368.56472 Aug 1 1946.9117 23081.50 71.58650 Sep 1 2540.9140 23290.75 -53.66661 Oct 1 -3756.4042 23476.19 698.21503 Nov 1 -4558.9237 23661.63 -414.70205 Dec 1 -2932.4822 23867.55 -511.06579 Jan 2 -295.1638 24073.47 997.69355 Feb 2 -4491.9076 24233.51 72.39372 Mar 2 -10466.5723 24393.56 -1188.98525 Apr 2 5221.3996 24419.03 1925.57066 May 2 3579.9725 24444.50 2086.52565 Jun 2 7856.9420 24348.59 -2186.52753 Jul 2 5355.3138 24252.67 2326.01696 Aug 2 1946.9117 24004.17 -125.08345 Sep 2 2540.9140 23755.67 538.41181 Oct 2 -3756.4042 23388.71 572.69222 Nov 2 -4558.9237 23021.75 -673.82608 Dec 2 -2932.4822 22611.27 841.21141 Jan 3 -295.1638 22200.79 612.37199 Feb 3 -4491.9076 21826.90 -1762.98780 Mar 3 -10466.5723 21453.00 522.57326 Apr 3 5221.3996 21234.52 -1008.91741 May 3 3579.9725 21016.04 -506.00902 Jun 3 7856.9420 21000.38 -1071.32028 Jul 3 5355.3138 20984.72 -145.03387 Aug 3 1946.9117 21156.36 -2587.27506 Sep 3 2540.9140 21328.01 -1109.92058 Oct 3 -3756.4042 21613.39 1171.01343 Nov 3 -4558.9237 21898.77 -368.85127 Dec 3 -2932.4822 22317.25 651.22731 Jan 4 -295.1638 22735.73 44.42898 Feb 4 -4491.9076 23210.15 11.76117 Mar 4 -10466.5723 23684.56 1320.01421 Apr 4 5221.3996 24052.54 -1712.94400 May 4 3579.9725 24420.53 -2015.50313 Jun 4 7856.9420 24631.78 2181.27812 Jul 4 5355.3138 24843.03 1867.65703 Aug 4 1946.9117 24996.79 242.30251 Sep 4 2540.9140 25150.54 1894.54367 Oct 4 -3756.4042 25205.50 -90.09764 Nov 4 -4558.9237 25260.46 851.46233 Dec 4 -2932.4822 25186.56 -2681.07638 Jan 5 -295.1638 25112.66 -561.49200 Feb 5 -4491.9076 24968.93 1902.97517 Mar 5 -10466.5723 24825.21 1808.36321 Apr 5 5221.3996 24710.76 -2635.16094 May 5 3579.9725 24596.31 110.71399 Jun 5 7856.9420 24517.10 1099.95319 Jul 5 5355.3138 24437.90 -1564.20995 Aug 5 1946.9117 24357.43 2480.65709 Sep 5 2540.9140 24276.97 -1220.88021 Oct 5 -3756.4042 24193.41 -2307.00259 Nov 5 -4558.9237 24109.85 647.07631 Dec 5 -2932.4822 24028.59 1752.88796 Jan 6 -295.1638 23947.34 -534.17730 > m$win s t l 611 19 13 > m$deg s t l 0 1 1 > m$jump s t l 62 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/rcomp/tmp/1f31u1322499506.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/www/rcomp/tmp/2kohb1322499506.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/www/rcomp/tmp/36npz1322499506.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/www/rcomp/tmp/44a221322499506.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/www/rcomp/tmp/5ja9w1322499506.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/www/rcomp/tmp/695ac1322499506.tab") > > try(system("convert tmp/1f31u1322499506.ps tmp/1f31u1322499506.png",intern=TRUE)) character(0) > try(system("convert tmp/2kohb1322499506.ps tmp/2kohb1322499506.png",intern=TRUE)) character(0) > try(system("convert tmp/36npz1322499506.ps tmp/36npz1322499506.png",intern=TRUE)) character(0) > try(system("convert tmp/44a221322499506.ps tmp/44a221322499506.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.680 0.272 2.038