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(12.610,10.862,52.929,56.902,81.776,87.876,82.103,72.846,60.632,33.521,15.342,7.758,8.668,13.082,38.157,58.263,81.153,88.476,72.329,75.845,61.108,37.665,12.755,2.793,12.935,19.533,33.404,52.074,70.735,69.702,61.656,82.993,53.990,32.283,15.686,2.713,12.842,19.244,48.488,54.464,84.192,84.458,85.793,75.163,68.212,49.233,24.302,5.402,15.058,33.559,70.358,85.934,94.452,129.305,113.882,107.256,94.274,57.842,26.611,14.521) > 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 -37.365440 48.52701 1.4484287 Feb 1 -30.983514 48.33852 -6.4930110 Mar 1 -2.023796 48.15004 6.8027579 Apr 1 10.705189 47.97356 -1.7767530 May 1 31.508174 47.79709 2.4707363 Jun 1 40.853703 47.63647 -0.6141730 Jul 1 31.886644 47.47585 2.7405062 Aug 1 31.190103 47.30396 -5.6480670 Sep 1 15.648150 47.13208 -2.1482277 Oct 1 -10.398977 46.89363 -2.9736497 Nov 1 -34.081308 46.65518 2.7681316 Dec 1 -46.938933 46.52030 8.1766298 Jan 2 -37.365440 46.38543 -0.3519902 Feb 2 -30.983514 46.30136 -2.2358454 Mar 2 -2.023796 46.21729 -6.0364918 Apr 2 10.705189 46.12923 1.4285765 May 2 31.508174 46.04118 3.6036450 Jun 2 40.853703 46.09161 1.5306901 Jul 2 31.886644 46.14203 -5.6996764 Aug 2 31.190103 46.21556 -1.5606585 Sep 2 15.648150 46.28908 -0.8292280 Oct 2 -10.398977 45.91105 2.1529275 Nov 2 -34.081308 45.53302 1.3032864 Dec 2 -46.938933 44.71685 5.0150840 Jan 3 -37.365440 43.90068 6.3997633 Feb 3 -30.983514 43.22831 7.2882046 Mar 3 -2.023796 42.55594 -7.1281455 Apr 3 10.705189 42.12019 -0.7513820 May 3 31.508174 41.68444 -2.4576182 Jun 3 40.853703 41.69026 -12.8419620 Jul 3 31.886644 41.69607 -11.9267173 Aug 3 31.190103 42.26581 9.5370875 Sep 3 15.648150 42.83554 -4.4936950 Oct 3 -10.398977 43.87279 -1.1908158 Nov 3 -34.081308 44.91004 4.8572668 Dec 3 -46.938933 45.96126 3.6906727 Jan 4 -37.365440 47.01248 3.1949603 Feb 4 -30.983514 47.77230 2.4552090 Mar 4 -2.023796 48.53213 1.9796664 Apr 4 10.705189 49.19099 -5.4321759 May 4 31.508174 49.84984 2.8339820 Jun 4 40.853703 50.47511 -6.8708121 Jul 4 31.886644 51.10037 2.8059822 Aug 4 31.190103 52.23647 -8.2635685 Sep 4 15.648150 53.37256 -0.8087066 Oct 4 -10.398977 55.26227 4.3697070 Nov 4 -34.081308 57.15198 1.2313240 Dec 4 -46.938933 59.65385 -7.3129119 Jan 5 -37.365440 62.15571 -9.7322661 Feb 5 -30.983514 64.58990 -0.0473842 Mar 5 -2.023796 67.02409 5.3577064 Apr 5 10.705189 67.96276 7.2660490 May 5 31.508174 68.90143 -5.9576082 Jun 5 40.853703 69.56716 18.8841402 Jul 5 31.886644 70.23288 11.7624770 Aug 5 31.190103 70.84381 5.2220901 Sep 5 15.648150 71.45473 7.1711158 Oct 5 -10.398977 71.88038 -3.6394064 Nov 5 -34.081308 72.30603 -11.6137253 Dec 5 -46.938933 72.52751 -11.0675723 > m$win s t l 601 19 13 > m$deg s t l 0 1 1 > m$jump s t l 61 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/rcomp/tmp/1c6u51259922693.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/2nbkd1259922693.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/3fl3c1259922693.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/41v0r1259922693.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/5c3v91259922693.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/63l2l1259922693.tab") > system("convert tmp/1c6u51259922693.ps tmp/1c6u51259922693.png") > system("convert tmp/2nbkd1259922693.ps tmp/2nbkd1259922693.png") > system("convert tmp/3fl3c1259922693.ps tmp/3fl3c1259922693.png") > system("convert tmp/41v0r1259922693.ps tmp/41v0r1259922693.png") > > > proc.time() user system elapsed 0.942 0.592 1.154