R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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(1038.00,934.00,988.00,870.00,854.00,834.00,872.00,954.00,870.00,1238.00,1082.00,1053.00,934.00,787.00,1081.00,908.00,995.00,825.00,822.00,856.00,887.00,1094.00,990.00,936.00,1097.00,918.00,926.00,907.00,899.00,971.00,1087.00,1000.00,1071.00,1190.00,1116.00,1070.00,1314.00,1068.00,1185.00,1215.00,1145.00,1251.00,1363.00,1368.00,1535.00,1853.00,1866.00,2023.00,1373.00,1968.00,1424.00,1160.00,1243.00,1375.00,1539.00,1773.00,1906.00,2076.00,2004.00) > par8 = 'FALSE' > par7 = '1' > par6 = '' > par5 = '1' > par4 = '' > par3 = '0' > par2 = '6' > par1 = '12' > main = 'Seasonal Decomposition by Loess' > 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 64.707689 964.0392 9.253089 Feb 1 -72.157452 963.4273 42.730175 Mar 1 22.070221 962.8153 3.114446 Apr 1 -77.968060 962.2034 -14.235329 May 1 -76.394046 961.5594 -31.165393 Jun 1 -103.917952 960.9155 -22.997538 Jul 1 -53.748781 960.2715 -34.522760 Aug 1 -42.356313 959.3669 36.989418 Sep 1 -23.019961 958.4622 -65.442288 Oct 1 214.368709 957.5576 66.073689 Nov 1 109.748194 956.8245 15.427349 Dec 1 46.331793 956.0913 50.576896 Jan 2 65.200296 955.3582 -86.558462 Feb 2 -66.603416 952.2573 -98.653858 Mar 2 11.936916 949.1564 119.906703 Apr 2 -84.793011 946.0555 46.737523 May 2 -87.747128 944.5314 138.215739 Jun 2 -105.833448 943.0073 -12.173844 Jul 2 -51.125712 941.4832 -68.357482 Aug 2 -42.404507 944.0535 -45.648982 Sep 2 -12.897379 946.6238 -46.726406 Oct 2 218.368997 949.1941 -73.563078 Nov 2 119.605243 954.7929 -84.398097 Dec 2 55.621721 960.3916 -80.013349 Jan 3 60.490148 965.9904 70.519450 Feb 3 -55.139151 974.0344 -0.895215 Mar 3 -9.095788 982.0783 -46.982542 Apr 3 -101.979514 990.1223 18.857220 May 3 -110.162724 1004.3862 4.776547 Jun 3 -111.515758 1018.6501 63.865699 Jul 3 -46.875910 1032.9139 100.961970 Aug 3 -46.366415 1054.2459 -7.879509 Sep 3 13.559079 1075.5779 -18.136987 Oct 3 220.029115 1096.9099 -126.939007 Nov 3 137.642032 1124.0636 -145.705662 Dec 3 118.939685 1151.2174 -200.157054 Jan 4 25.160906 1178.3711 110.467988 Feb 4 17.864939 1214.6597 -164.524647 Mar 4 -47.535464 1250.9483 -18.412845 Apr 4 -163.905867 1287.2369 91.668956 May 4 -168.503392 1327.4055 -13.902069 Jun 4 -139.781054 1367.5740 23.207043 Jul 4 -57.602157 1407.7426 12.859596 Aug 4 -27.317507 1436.9599 -41.642392 Sep 4 51.307481 1466.1772 17.515283 Oct 4 245.993392 1495.3946 111.612034 Nov 4 169.193784 1516.8313 179.974894 Dec 4 143.739210 1538.2681 340.992721 Jan 5 13.646218 1559.7048 -200.351034 Feb 5 33.055981 1582.2874 352.656643 Mar 5 -65.502493 1604.8699 -115.367442 Apr 5 -185.279339 1627.4525 -282.173156 May 5 -192.091227 1648.0343 -212.943081 Jun 5 -151.149403 1668.6161 -142.466718 Jul 5 -62.729939 1689.1979 -87.467996 Aug 5 -24.927704 1709.0907 88.836956 Sep 5 61.770855 1728.9836 115.245584 Oct 5 254.854635 1748.8764 72.268990 Nov 5 182.097495 1769.1857 52.716852 > m$win s t l 6 25 13 > m$deg s t l 0 1 1 > m$jump s t l 1 3 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/freestat/rcomp/tmp/1zfym1292884111.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/freestat/rcomp/tmp/2soyp1292884111.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/freestat/rcomp/tmp/3soyp1292884111.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/freestat/rcomp/tmp/4lyfs1292884111.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/5z7vj1292884111.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/freestat/rcomp/tmp/6k8t71292884111.tab") > postscript(file="/var/www/html/freestat/rcomp/tmp/7dhta1292884111.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > myresid <- m$time.series[!is.na(m$time.series[,'remainder']),'remainder'] > hist(as.numeric(myresid), main='Residual Histogram', xlab='Residual Value') > dev.off() null device 1 > > try(system("convert tmp/1zfym1292884111.ps tmp/1zfym1292884111.png",intern=TRUE)) character(0) > try(system("convert tmp/2soyp1292884111.ps tmp/2soyp1292884111.png",intern=TRUE)) character(0) > try(system("convert tmp/3soyp1292884111.ps tmp/3soyp1292884111.png",intern=TRUE)) character(0) > try(system("convert tmp/4lyfs1292884111.ps tmp/4lyfs1292884111.png",intern=TRUE)) character(0) > try(system("convert tmp/7dhta1292884111.ps tmp/7dhta1292884111.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.620 1.046 1.761