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(96.2,96.8,109.9,88,91.1,106.4,68.6,100.1,108,106,108.6,91.5,99.2,98,96.6,102.8,96.9,110,70.5,101.9,109.6,107.8,113,93.8,108,102.8,116.3,89.2,106.7,112.1,74.2,108.8,111.5,118.8,118.9,97.6,116.4,107.9,121.2,97.9,113.4,117.6,79.6,115.9,115.7,129.1,123.3,96.7,121.2,118.2,102.1,125.4,116.7,121.3,85.3,114.2,124.4,131,118.3,99.6) > 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 3.8137034 96.02648 -3.6401802 Feb 1 0.0225964 96.25981 0.5175887 Mar 1 4.1714973 96.49315 9.2353496 Apr 1 -4.6108956 96.72709 -4.1161923 May 1 -0.5332715 96.96102 -5.3277511 Jun 1 7.7903802 97.17829 1.4313272 Jul 1 -30.2459673 97.39556 1.4504048 Aug 1 2.0623729 97.58115 0.4564789 Sep 1 7.4907078 97.76673 2.7425583 Oct 1 11.9071980 97.96819 -3.8753893 Nov 1 9.5036885 98.16965 0.9266627 Dec 1 -11.3720131 98.48073 4.3912860 Jan 2 3.8137034 98.79181 -3.4055087 Feb 2 0.0225964 99.00641 -1.0290111 Mar 2 4.1714973 99.22102 -6.7925215 Apr 2 -4.6108956 99.45566 7.9552342 May 2 -0.5332715 99.69030 -2.2570271 Jun 2 7.7903802 100.08404 2.1255810 Jul 2 -30.2459673 100.47778 0.2681884 Aug 2 2.0623729 101.00774 -1.1701175 Sep 2 7.4907078 101.53771 0.5715818 Oct 2 11.9071980 101.91027 -6.0174721 Nov 2 9.5036885 102.28284 1.2134737 Dec 2 -11.3720131 102.61292 2.5590886 Jan 3 3.8137034 102.94301 1.2432854 Feb 3 0.0225964 103.35194 -0.5745393 Mar 3 4.1714973 103.76087 8.3676281 Apr 3 -4.6108956 104.21261 -10.4017127 May 3 -0.5332715 104.66434 2.5689296 Jun 3 7.7903802 105.13094 -0.8213246 Jul 3 -30.2459673 105.59755 -1.1515796 Aug 3 2.0623729 106.17201 0.5656213 Sep 3 7.4907078 106.74646 -2.7371726 Oct 3 11.9071980 107.35672 -0.4639159 Nov 3 9.5036885 107.96697 1.4293406 Dec 3 -11.3720131 108.51744 0.4545756 Jan 4 3.8137034 109.06790 3.5183926 Feb 4 0.0225964 109.55034 -1.6729386 Mar 4 4.1714973 110.03278 6.9957223 Apr 4 -4.6108956 110.44748 -7.9365867 May 4 -0.5332715 110.86218 3.0710872 Jun 4 7.7903802 111.17050 -1.3608818 Jul 4 -30.2459673 111.47882 -1.6328516 Aug 4 2.0623729 111.74603 2.0916014 Sep 4 7.4907078 112.01323 -3.8039404 Oct 4 11.9071980 112.48561 4.7071909 Nov 4 9.5036885 112.95799 0.8383218 Dec 4 -11.3720131 113.43102 -5.3590089 Jan 5 3.8137034 113.90405 3.4822423 Feb 5 0.0225964 114.23831 3.9390947 Mar 5 4.1714973 114.57256 -16.6440607 Apr 5 -4.6108956 114.57537 15.4355222 May 5 -0.5332715 114.57818 2.6550881 Jun 5 7.7903802 114.50825 -0.9986265 Jul 5 -30.2459673 114.43831 1.1076582 Aug 5 2.0623729 114.34075 -2.2031245 Sep 5 7.4907078 114.24319 2.6660979 Oct 5 11.9071980 114.10884 4.9839645 Nov 5 9.5036885 113.97448 -5.1781693 Dec 5 -11.3720131 113.79995 -2.8279406 > 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/19d2u1259957654.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/2wdbs1259957654.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/36p1l1259957654.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/4la1a1259957654.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/5l1b81259957654.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/60na01259957654.tab") > > system("convert tmp/19d2u1259957654.ps tmp/19d2u1259957654.png") > system("convert tmp/2wdbs1259957654.ps tmp/2wdbs1259957654.png") > system("convert tmp/36p1l1259957654.ps tmp/36p1l1259957654.png") > system("convert tmp/4la1a1259957654.ps tmp/4la1a1259957654.png") > > > proc.time() user system elapsed 0.972 0.612 1.136