R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-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(6827,6178,7084,8162,8462,9644,10466,10748,9963,8194,6848,7027,7269,6775,7819,8371,9069,10248,11030,10882,10333,9109,7685,7602,8350,7829,8829,9948,10638,11253,11424,11391,10665,9396,7775,7933,8186,7444,8484,9864,10252,12282,11637,11577,12417,9637,8094,9280,8334,7899,9994,10078,10801,12950,12222,12246,13281,10366,8730,9614,8639,8772,10894,10455,11179,10588,10794,12770,13812,10857,9290,10925,9491,8919,11607,8852,12537,14759,13667,13731,15110,12185,10645,12161,10840,10436,13589,13402,13103,14933,14147,14057,16234,12389,11595,12772) > 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 -1623.7651 8150.612 300.15270 Feb 1 -2138.5373 8183.557 132.98021 Mar 1 -436.8098 8216.502 -695.69199 Apr 1 -380.1758 8248.449 293.72653 May 1 436.0833 8280.397 -254.47996 Jun 1 1720.9295 8311.334 -388.26345 Jul 1 1520.0256 8342.271 603.70314 Aug 1 1726.7507 8377.255 643.99435 Sep 1 2233.2257 8412.239 -682.46428 Oct 1 -283.7440 8454.123 23.62128 Nov 1 -1774.3390 8496.007 126.33219 Dec 1 -999.6436 8533.710 -507.06608 Jan 2 -1623.7651 8571.412 321.35265 Feb 2 -2138.5373 8611.870 301.66694 Mar 2 -436.8098 8652.328 -396.51847 Apr 2 -380.1758 8707.966 43.20956 May 2 436.0833 8763.604 -130.68741 Jun 2 1720.9295 8829.101 -302.03028 Jul 2 1520.0256 8894.597 615.37693 Aug 2 1726.7507 8981.915 173.33450 Sep 2 2233.2257 9069.232 -969.45778 Oct 2 -283.7440 9175.760 216.98416 Nov 2 -1774.3390 9282.288 177.05145 Dec 2 -999.6436 9376.156 -774.51267 Jan 3 -1623.7651 9470.025 503.74023 Feb 3 -2138.5373 9525.585 441.95189 Mar 3 -436.8098 9581.146 -315.33616 Apr 3 -380.1758 9603.042 725.13403 May 3 436.0833 9624.937 576.97921 Jun 3 1720.9295 9614.492 -82.42150 Jul 3 1520.0256 9604.047 299.92787 Aug 3 1726.7507 9573.056 91.19342 Sep 3 2233.2257 9542.065 -1110.29087 Oct 3 -283.7440 9526.009 153.73472 Nov 3 -1774.3390 9509.953 39.38565 Dec 3 -999.6436 9539.269 -606.62511 Jan 4 -1623.7651 9568.584 241.18114 Feb 4 -2138.5373 9633.706 -51.16860 Mar 4 -436.8098 9698.828 -778.01803 Apr 4 -380.1758 9766.455 477.72075 May 4 436.0833 9834.082 -18.16547 Jun 4 1720.9295 9890.236 670.83481 Jul 4 1520.0256 9946.389 170.58517 Aug 4 1726.7507 9995.536 -145.28669 Sep 4 2233.2257 10044.683 139.09161 Oct 4 -283.7440 10088.205 -167.46093 Nov 4 -1774.3390 10131.727 -263.38812 Dec 4 -999.6436 10181.099 98.54462 Jan 5 -1623.7651 10230.471 -272.70563 Feb 5 -2138.5373 10297.034 -259.49684 Mar 5 -436.8098 10363.598 67.21225 Apr 5 -380.1758 10429.098 29.07739 May 5 436.0833 10494.599 -129.68247 Jun 5 1720.9295 10543.248 685.82297 Jul 5 1520.0256 10591.896 110.07849 Aug 5 1726.7507 10633.990 -114.74053 Sep 5 2233.2257 10676.084 371.69061 Oct 5 -283.7440 10688.615 -38.87118 Nov 5 -1774.3390 10701.147 -196.80762 Dec 5 -999.6436 10651.705 -38.06189 Jan 6 -1623.7651 10602.264 -339.49915 Feb 6 -2138.5373 10569.092 341.44522 Mar 6 -436.8098 10535.920 794.88989 Apr 6 -380.1758 10563.579 271.59628 May 6 436.0833 10591.239 151.67766 Jun 6 1720.9295 10651.896 -1784.82539 Jul 6 1520.0256 10712.553 -1438.57836 Aug 6 1726.7507 10765.989 277.26074 Sep 6 2233.2257 10819.424 759.35000 Oct 6 -283.7440 10898.237 242.50697 Nov 6 -1774.3390 10977.050 87.28929 Dec 6 -999.6436 11125.072 799.57166 Jan 7 -1623.7651 11273.094 -158.32897 Feb 7 -2138.5373 11415.235 -357.69774 Mar 7 -436.8098 11557.376 486.43379 Apr 7 -380.1758 11673.469 -2441.29307 May 7 436.0833 11789.562 311.35508 Jun 7 1720.9295 11919.243 1118.82766 Jul 7 1520.0256 12048.924 98.05031 Aug 7 1726.7507 12210.751 -206.50176 Sep 7 2233.2257 12372.578 504.19633 Oct 7 -283.7440 12530.193 -61.44901 Nov 7 -1774.3390 12687.808 -268.46901 Dec 7 -999.6436 12774.019 386.62482 Jan 8 -1623.7651 12860.229 -396.46434 Feb 8 -2138.5373 12909.111 -334.57406 Mar 8 -436.8098 12957.993 1067.81653 Apr 8 -380.1758 12999.190 782.98556 May 8 436.0833 13040.387 -373.47040 Jun 8 1720.9295 13089.116 122.95481 Jul 8 1520.0256 13137.844 -510.86990 Aug 8 1726.7507 13181.107 -850.85745 Sep 8 2233.2257 13224.369 776.40518 Oct 8 -283.7440 13264.776 -592.03207 Nov 8 -1774.3390 13305.183 64.15603 Dec 8 -999.6436 13345.787 425.85657 > m$win s t l 961 19 13 > m$deg s t l 0 1 1 > m$jump s t l 97 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/wessaorg/rcomp/tmp/1ju251324215486.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/2tg0w1324215486.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/3mrnj1324215486.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/4o50f1324215486.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/5aguh1324215486.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/6m43z1324215486.tab") > > try(system("convert tmp/1ju251324215486.ps tmp/1ju251324215486.png",intern=TRUE)) character(0) > try(system("convert tmp/2tg0w1324215486.ps tmp/2tg0w1324215486.png",intern=TRUE)) character(0) > try(system("convert tmp/3mrnj1324215486.ps tmp/3mrnj1324215486.png",intern=TRUE)) character(0) > try(system("convert tmp/4o50f1324215486.ps tmp/4o50f1324215486.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.392 0.256 1.650