R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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(41086,39690,43129,37863,35953,29133,24693,22205,21725,27192,21790,13253,37702,30364,32609,30212,29965,28352,25814,22414,20506,28806,22228,13971,36845,35338,35022,34777,26887,23970,22780,17351,21382,24561,17409,11514,31514,27071,29462,26105,22397,23843,21705,18089,20764,25316,17704,15548,28029,29383,36438,32034,22679,24319,18004,17537,20366,22782,19169,13807,29743,25591,29096,26482,22405,27044,17970,18730,19684,19785,18479,10698) > par8 = 'TRUE' > 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 8106.6521 34441.92 -1462.56763 Feb 1 5114.1533 33619.33 956.51839 Mar 1 8230.8769 32796.74 2101.38193 Apr 1 5525.9447 32009.75 327.30719 May 1 1087.1584 31222.75 3643.08670 Jun 1 638.9884 30473.50 -1979.48889 Jul 1 -3282.5715 29724.25 -1748.67448 Aug 1 -5559.2879 28974.69 -1210.40117 Sep 1 -3873.7462 28225.13 -2626.38598 Oct 1 216.8232 27598.45 -623.27720 Nov 1 -4965.8237 26971.78 -215.95203 Dec 1 -11239.1658 26802.93 -2310.76819 Jan 2 8106.6521 26634.09 2961.25566 Feb 2 5114.1533 26675.37 -1425.52409 Mar 2 8230.8769 26716.65 -2338.52632 Apr 2 5525.9447 26786.82 -2100.76757 May 2 1087.1584 26857.00 2020.84543 Jun 2 638.9884 26956.51 756.49832 Jul 2 -3282.5715 27056.03 2040.54123 Aug 2 -5559.2879 27214.25 759.03673 Sep 2 -3873.7462 27372.47 -2992.72589 Oct 2 216.8232 27433.93 1155.24861 Nov 2 -4965.8237 27495.38 -301.56050 Dec 2 -11239.1658 27326.99 -2116.82332 Jan 3 8106.6521 27158.59 1579.75388 Feb 3 5114.1533 26909.52 3314.32528 Mar 3 8230.8769 26660.45 130.67421 Apr 3 5525.9447 26362.74 2888.31316 May 3 1087.1584 26065.04 -265.19364 Jun 3 638.9884 25657.11 -2326.09418 Jul 3 -3282.5715 25249.18 813.39529 Aug 3 -5559.2879 24741.83 -1831.54166 Sep 3 -3873.7462 24234.48 1021.26328 Oct 3 216.8232 23787.18 556.99693 Nov 3 -4965.8237 23339.88 -965.05302 Dec 3 -11239.1658 23103.64 -350.47704 Jan 4 8106.6521 22867.41 539.93895 Feb 4 5114.1533 22823.66 -866.80949 Mar 4 8230.8769 22779.90 -1548.78042 Apr 4 5525.9447 22855.84 -2276.78728 May 4 1087.1584 22931.78 -1621.93990 Jun 4 638.9884 23087.57 116.44506 Jul 4 -3282.5715 23243.35 1744.22002 Aug 4 -5559.2879 23493.80 154.48757 Sep 4 -3873.7462 23744.25 893.49701 Oct 4 216.8232 23967.34 1131.83563 Nov 4 -4965.8237 24190.43 -1520.60935 Dec 4 -11239.1658 24196.18 2590.98393 Jan 5 8106.6521 24201.93 -4279.58279 Feb 5 5114.1533 24089.83 179.01595 Mar 5 8230.8769 23977.73 4229.39221 Apr 5 5525.9447 23893.50 2614.55366 May 5 1087.1584 23809.27 -2217.43065 Jun 5 638.9884 23707.32 -27.31149 Jul 5 -3282.5715 23605.37 -2318.80232 Aug 5 -5559.2879 23377.61 -281.31799 Sep 5 -3873.7462 23149.84 1089.90822 Oct 5 216.8232 22904.34 -339.16302 Nov 5 -4965.8237 22658.84 1475.98215 Dec 5 -11239.1658 22572.63 2473.53300 Jan 6 8106.6521 22486.42 -850.07613 Feb 6 5114.1533 22450.96 -1974.11445 Mar 6 8230.8769 22415.50 -1550.37525 Apr 6 5525.9447 22393.04 -1436.98940 May 6 1087.1584 22370.59 -1052.74930 Jun 6 638.9884 22350.15 4054.86619 Jul 6 -3282.5715 22329.70 -1077.12832 Aug 6 -5559.2879 22332.94 1956.34971 Sep 6 -3873.7462 22336.18 1221.56962 Oct 6 216.8232 22353.58 -2785.40512 Nov 6 -4965.8237 22370.99 1073.83654 Dec 6 -11239.1658 22392.11 -454.94447 > m$win s t l 721 19 13 > m$deg s t l 0 1 1 > m$jump s t l 73 2 2 > m$inner [1] 1 > m$outer [1] 15 > postscript(file="/var/www/html/rcomp/tmp/1avqa1243876019.ps",horizontal=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/html/rcomp/tmp/2xqwk1243876019.ps",horizontal=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/html/rcomp/tmp/38vhg1243876019.ps",horizontal=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/html/rcomp/tmp/47jub1243876019.ps",horizontal=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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/html/rcomp/tmp/5ys201243876019.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/html/rcomp/tmp/6k1wn1243876019.tab") > > system("convert tmp/1avqa1243876019.ps tmp/1avqa1243876019.png") > system("convert tmp/2xqwk1243876019.ps tmp/2xqwk1243876019.png") > system("convert tmp/38vhg1243876019.ps tmp/38vhg1243876019.png") > system("convert tmp/47jub1243876019.ps tmp/47jub1243876019.png") > > > proc.time() user system elapsed 1.019 0.624 1.541