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(235243,230354,227184,221678,217142,219452,256446,265845,248624,241114,229245,231805,219277,219313,212610,214771,211142,211457,240048,240636,230580,208795,197922,194596,194581,185686,178106,172608,167302,168053,202300,202388,182516,173476,166444,171297,169701,164182,161914,159612,151001,158114,186530,187069,174330,169362,166827,178037,186413,189226,191563,188906,186005,195309,223532,226899,214126,206903,204442,220375,214320,212588,205816,202196,195722,198563,229139,229527,211868,203555,195770) > 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 -1294.6913 235593.4 944.29251 Feb 1 -3809.2186 235420.9 -1257.65671 Mar 1 -7318.7443 235248.4 -745.60745 Apr 1 -9722.8344 234964.7 -3563.83160 May 1 -14466.0922 234681.0 -3072.88793 Jun 1 -10350.8923 234252.5 -4449.57818 Jul 1 21166.4717 233824.0 1455.56741 Aug 1 23856.0974 233320.4 8668.52622 Sep 1 9097.5600 232816.8 6709.64820 Oct 1 -293.0042 232123.4 9283.62112 Nov 1 -6969.5655 231430.0 4784.59115 Dec 1 104.9066 230177.4 1522.72763 Jan 2 -1294.6913 228924.8 -8353.06583 Feb 2 -3809.2186 227166.6 -4044.33247 Mar 2 -7318.7443 225408.3 -5479.60063 Apr 2 -9722.8344 223295.1 1198.76048 May 2 -14466.0922 221181.8 4426.28941 Jun 2 -10350.8923 218808.1 2999.77316 Jul 2 21166.4717 216434.4 2447.09277 Aug 2 23856.0974 213584.5 3195.42260 Sep 2 9097.5600 210734.5 10747.91561 Oct 2 -293.0042 207263.8 1824.23573 Nov 2 -6969.5655 203793.0 1098.55296 Dec 2 104.9066 200103.6 -5612.53482 Jan 3 -1294.6913 196414.2 -538.55255 Feb 3 -3809.2186 192950.3 -3455.08838 Mar 3 -7318.7443 189486.4 -4061.62572 Apr 3 -9722.8344 186568.1 -4237.21748 May 3 -14466.0922 183649.7 -1881.64140 Jun 3 -10350.8923 181469.0 -3065.09468 Jul 3 21166.4717 179288.2 1845.28788 Aug 3 23856.0974 177636.8 895.10820 Sep 3 9097.5600 175985.3 -2566.90832 Oct 3 -293.0042 174679.1 -910.11997 Nov 3 -6969.5655 173372.9 40.66549 Dec 3 104.9066 172230.7 -1038.59586 Jan 4 -1294.6913 171088.5 -92.78716 Feb 4 -3809.2186 170128.3 -2137.05581 Mar 4 -7318.7443 169168.1 64.67401 Apr 4 -9722.8344 168749.6 585.22426 May 4 -14466.0922 168331.1 -2864.05766 Jun 4 -10350.8923 168868.4 -403.52096 Jul 4 21166.4717 169405.7 -4042.14841 Aug 4 23856.0974 171159.0 -7946.12508 Sep 4 9097.5600 172912.4 -7679.93858 Oct 4 -293.0042 175645.5 -5990.44995 Nov 4 -6969.5655 178378.5 -4581.96422 Dec 4 104.9066 181702.4 -3770.26676 Jan 5 -1294.6913 185026.2 2681.50076 Feb 5 -3809.2186 188456.4 4578.81627 Mar 5 -7318.7443 191886.6 6995.13026 Apr 5 -9722.8344 195046.5 3582.36822 May 5 -14466.0922 198206.3 2264.77400 Jun 5 -10350.8923 200875.0 4784.90939 Jul 5 21166.4717 203543.6 -1178.11937 Aug 5 23856.0974 205484.4 -2441.50712 Sep 5 9097.5600 207425.2 -2396.73170 Oct 5 -293.0042 208651.4 -1455.39283 Nov 5 -6969.5655 209877.6 1533.94314 Dec 5 104.9066 210539.9 9730.22308 Jan 6 -1294.6913 211202.1 4412.57307 Feb 6 -3809.2186 210981.1 5416.13486 Mar 6 -7318.7443 210760.0 2374.69513 Apr 6 -9722.8344 210003.9 1914.92959 May 6 -14466.0922 209247.8 940.33188 Jun 6 -10350.8923 208400.8 513.11080 Jul 6 21166.4717 207553.8 418.72557 Aug 6 23856.0974 206598.5 -927.62991 Sep 6 9097.5600 205643.3 -2872.82222 Oct 6 -293.0042 204595.4 -747.37306 Nov 6 -6969.5655 203547.5 -807.92678 > m$win s t l 711 19 13 > m$deg s t l 0 1 1 > m$jump s t l 72 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/rcomp/tmp/18pgx1293569253.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/html/rcomp/tmp/28pgx1293569253.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/html/rcomp/tmp/3ojn71293569253.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/html/rcomp/tmp/4ojn71293569253.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/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/5f8dr1293569253.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/60qcf1293569253.tab") > > try(system("convert tmp/18pgx1293569253.ps tmp/18pgx1293569253.png",intern=TRUE)) character(0) > try(system("convert tmp/28pgx1293569253.ps tmp/28pgx1293569253.png",intern=TRUE)) character(0) > try(system("convert tmp/3ojn71293569253.ps tmp/3ojn71293569253.png",intern=TRUE)) character(0) > try(system("convert tmp/4ojn71293569253.ps tmp/4ojn71293569253.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.991 0.628 4.296