R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing 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(56,55,54,52,72,71,56,46,47,47,48,50,44,38,33,33,52,54,39,22,31,31,38,42,41,31,36,34,51,47,31,19,30,33,36,40,32,25,28,29,55,55,40,38,44,41,49,59,61,47,43,39,66,68,63,68,67,59,68,78,82,70,62,68,94,102,100,104,103,93,110,114,120,102,95,103,122,139,135,135,137,130,148,148,145,128,131,133,146,163,151,157,152,149,172,167,160,150,160,165,171,179,171,176,170,169,194,196,188,174,186,191,197,206,197,204,201,190,213,213) > par8 = 'FALSE' > par7 = '1' > par6 = '' > par5 = '1' > par4 = '' > par3 = '0' > par2 = 'periodic' > par1 = '12' > main = 'Seasonal Decomposition by Loess' > par8 <- 'FALSE' > par7 <- '1' > par6 <- '' > par5 <- '1' > par4 <- '' > par3 <- '0' > par2 <- 'periodic' > par1 <- '12' > #'GNU S' R Code compiled by R2WASP v. 1.2.327 () > #Author: root > #To cite this work: Wessa P., (2013), Decomposition by Loess (v1.0.2) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_decomposeloess.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > # > 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 2.1004552 64.07227 -10.17272215 Feb 1 -9.9402932 62.41446 2.52583121 Mar 1 -10.2810431 60.75666 3.52438619 Apr 1 -9.3688316 59.10796 2.26087403 May 1 7.5433860 57.45926 6.99735583 Jun 1 12.4106282 55.86571 2.72365687 Jul 1 1.3778719 54.27217 0.34995640 Aug 1 -1.0425457 52.66456 -5.62201561 Sep 1 -0.7629647 51.05695 -3.29398622 Oct 1 -6.0124822 49.33423 3.67825530 Nov 1 6.1380015 47.61150 -5.74950443 Dec 1 7.8378200 46.19623 -4.03404631 Jan 2 2.1004552 44.78095 -2.88140497 Feb 2 -9.9402932 43.39959 4.54070739 Mar 2 -10.2810431 42.01822 1.26282137 Apr 2 -9.3688316 40.76879 1.60004332 May 2 7.5433860 39.51935 4.93725924 Jun 2 12.4106282 38.63485 2.95452224 Jul 2 1.3778719 37.75034 -0.12821628 Aug 2 -1.0425457 37.31560 -14.27305722 Sep 2 -0.7629647 36.88086 -5.11789677 Oct 2 -6.0124822 36.88680 0.12567725 Nov 2 6.1380015 36.89275 -5.03074998 Dec 2 7.8378200 36.95415 -2.79197297 Jan 3 2.1004552 37.01556 1.88398727 Feb 3 -9.9402932 36.90522 4.03507416 Mar 3 -10.2810431 36.79488 9.48616266 Apr 3 -9.3688316 36.49103 6.87780412 May 3 7.5433860 36.18717 7.26943956 Jun 3 12.4106282 35.60967 -1.02029500 Jul 3 1.3778719 35.03216 -5.41003106 Aug 3 -1.0425457 34.34156 -14.29901858 Sep 3 -0.7629647 33.65097 -2.88800472 Oct 3 -6.0124822 33.54436 5.46812477 Nov 3 6.1380015 33.43775 -3.57574699 Dec 3 7.8378200 34.18359 -2.02140930 Jan 4 2.1004552 34.92943 -5.02988838 Feb 4 -9.9402932 36.06511 -1.12481712 Mar 4 -10.2810431 37.20079 1.08025576 Apr 4 -9.3688316 38.34838 0.02045007 May 4 7.5433860 39.49598 7.96063836 Jun 4 12.4106282 40.98244 1.60693253 Jul 4 1.3778719 42.46890 -3.84677480 Aug 4 -1.0425457 44.01715 -4.97459970 Sep 4 -0.7629647 45.56539 -0.80242319 Oct 4 -6.0124822 46.86442 0.14805879 Nov 4 6.1380015 48.16346 -5.30146048 Dec 4 7.8378200 49.59624 1.56593759 Jan 5 2.1004552 51.02903 7.87051890 Feb 5 -9.9402932 52.86228 4.07801606 Mar 5 -10.2810431 54.69553 -1.41448515 Apr 5 -9.3688316 56.39500 -8.02616380 May 5 7.5433860 58.09446 0.36215152 Jun 5 12.4106282 59.75611 -4.16674092 Jul 5 1.3778719 61.41776 0.20436514 Aug 5 -1.0425457 63.28490 5.75764795 Sep 5 -0.7629647 65.15203 2.61093216 Oct 5 -6.0124822 67.36385 -2.35137086 Nov 5 6.1380015 69.57567 -7.71367514 Dec 5 7.8378200 72.11696 -1.95478162 Jan 6 2.1004552 74.65825 5.24129512 Feb 6 -9.9402932 77.66383 2.27646109 Mar 6 -10.2810431 80.66941 -8.38837133 Apr 6 -9.3688316 83.87548 -6.50664549 May 6 7.5433860 87.08154 -0.62492568 Jun 6 12.4106282 90.29635 -0.70697579 Jul 6 1.3778719 93.51116 5.11097261 Aug 6 -1.0425457 96.49285 8.54969820 Sep 6 -0.7629647 99.47454 4.28842519 Oct 6 -6.0124822 102.15092 -3.13843818 Nov 6 6.1380015 104.82730 -0.96530280 Dec 6 7.8378200 107.33150 -1.16931537 Jan 7 2.1004552 109.83569 8.06385529 Feb 7 -9.9402932 112.60221 -0.66191366 Mar 7 -10.2810431 115.36872 -10.08768099 Apr 7 -9.3688316 118.38817 -6.01933374 May 7 7.5433860 121.40761 -6.95099252 Jun 7 12.4106282 124.34116 2.24820853 Jul 7 1.3778719 127.27472 6.34740806 Aug 7 -1.0425457 129.96407 6.07847526 Sep 7 -0.7629647 132.65342 5.10954385 Oct 7 -6.0124822 134.91240 1.10007862 Nov 7 6.1380015 137.17139 4.69061213 Dec 7 7.8378200 138.82518 1.33699554 Jan 8 2.1004552 140.47898 2.42056218 Feb 8 -9.9402932 141.92151 -3.98121212 Mar 8 -10.2810431 143.36403 -2.08298480 Apr 8 -9.3688316 144.99595 -2.62712210 May 8 7.5433860 146.62788 -8.17126542 Jun 8 12.4106282 148.39909 2.19027823 Jul 8 1.3778719 150.17031 -0.54817962 Aug 8 -1.0425457 152.19865 5.84389560 Sep 8 -0.7629647 154.22699 -1.46402778 Oct 8 -6.0124822 156.33531 -1.32282521 Nov 8 6.1380015 158.44362 7.41837611 Dec 8 7.8378200 160.24681 -1.08462627 Jan 9 2.1004552 162.04999 -4.15044541 Feb 9 -9.9402932 163.57987 -3.63957315 Mar 9 -10.2810431 165.10974 5.17130073 Apr 9 -9.3688316 166.83101 7.53781830 May 9 7.5433860 168.55228 -5.09567016 Jun 9 12.4106282 170.56146 -3.97208661 Jul 9 1.3778719 172.57063 -2.94850458 Aug 9 -1.0425457 174.78208 2.26046924 Sep 9 -0.7629647 176.99352 -6.23055555 Oct 9 -6.0124822 179.26386 -4.25137861 Nov 9 6.1380015 181.53420 6.32779707 Dec 9 7.8378200 183.88523 4.27694768 Jan 10 2.1004552 186.23626 -0.33671849 Feb 10 -9.9402932 188.43459 -4.49429442 Mar 10 -10.2810431 190.63291 5.64813126 Apr 10 -9.3688316 192.37542 7.99341623 May 10 7.5433860 194.11792 -4.66130484 Jun 10 12.4106282 195.82528 -2.23590380 Jul 10 1.3778719 197.53263 -1.91050427 Aug 10 -1.0425457 199.20413 5.83841253 Sep 10 -0.7629647 200.87563 0.88733072 Oct 10 -6.0124822 202.51579 -6.50330281 Nov 10 6.1380015 204.15594 2.70606242 Dec 10 7.8378200 205.77580 -0.61361723 > m$win s t l 1201 19 13 > m$deg s t l 0 1 1 > m$jump s t l 121 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/wessaorg/rcomp/tmp/1nm3f1387554871.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/2sjz51387554871.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/3yv5c1387554871.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/4udd51387554871.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/5i2ks1387554871.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/60wyb1387554871.tab") > > try(system("convert tmp/1nm3f1387554871.ps tmp/1nm3f1387554871.png",intern=TRUE)) character(0) > try(system("convert tmp/2sjz51387554871.ps tmp/2sjz51387554871.png",intern=TRUE)) character(0) > try(system("convert tmp/3yv5c1387554871.ps tmp/3yv5c1387554871.png",intern=TRUE)) character(0) > try(system("convert tmp/4udd51387554871.ps tmp/4udd51387554871.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 5.243 0.899 6.121