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(464,460,467,460,448,443,436,431,484,510,513,503,471,471,476,475,470,461,455,456,517,525,523,519,509,512,519,517,510,509,501,507,569,580,578,565,547,555,562,561,555,544,537,543,594,611,613,611,594,595,591,589,584,573,567,569,621,629,628,612,595,597,593,590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516,528,533,536,537,524,536,587,597,581,564,558,575,580,575,563,552,537,545,601,604,586,564,549) > 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 -4.885406 460.4942 8.39118156 Feb 1 -1.823171 461.6339 0.18930370 Mar 1 1.305063 462.7735 2.92142675 Apr 1 -3.197795 464.0188 -0.82099030 May 1 -12.600652 465.2641 -4.66340746 Jun 1 -19.361145 466.5411 -4.17991883 Jul 1 -29.821637 467.8181 -1.99643229 Aug 1 -27.783047 469.0508 -10.26772643 Sep 1 25.455542 470.2835 -11.73902022 Oct 1 35.119009 471.7224 3.15855158 Nov 1 26.382479 473.1614 13.45612076 Dec 1 11.210755 474.8461 16.94311653 Jan 2 -4.885406 476.5309 -0.64545000 Feb 2 -1.823171 478.1689 -5.34570858 Mar 2 1.305063 479.8069 -5.11196626 Apr 2 -3.197795 481.2200 -3.02216073 May 2 -12.600652 482.6330 -0.03235531 Jun 2 -19.361145 484.6037 -4.24259674 Jul 2 -29.821637 486.5745 -1.75284024 Aug 2 -27.783047 489.7272 -5.94415590 Sep 2 25.455542 492.8799 -1.33547121 Oct 2 35.119009 496.5635 -6.68248502 Nov 2 26.382479 500.2470 -3.62950146 Dec 2 11.210755 504.1701 3.61916570 Jan 3 -4.885406 508.0931 5.79227055 Feb 3 -1.823171 512.2779 1.54532029 Mar 3 1.305063 516.4626 1.23237094 Apr 3 -3.197795 520.6793 -0.48146356 May 3 -12.600652 524.8960 -2.29529817 Jun 3 -19.361145 528.8141 -0.45297398 Jul 3 -29.821637 532.7323 -1.91065188 Aug 3 -27.783047 536.4043 -1.62123075 Sep 3 25.455542 540.0763 3.46819074 Oct 3 35.119009 543.6671 1.21391732 Nov 3 26.382479 547.2579 4.35964128 Dec 3 11.210755 550.4866 3.30260390 Jan 4 -4.885406 553.7154 -1.82999578 Feb 4 -1.823171 556.4164 0.40672854 Mar 4 1.305063 559.1175 1.57745377 Apr 4 -3.197795 561.7492 2.44856959 May 4 -12.600652 564.3810 3.21968529 Jun 4 -19.361145 567.5598 -4.19868585 Jul 4 -29.821637 570.7387 -3.91705907 Aug 4 -27.783047 574.0228 -3.23970610 Sep 4 25.455542 577.3068 -8.76235277 Oct 4 35.119009 580.1859 -4.30489152 Nov 4 26.382479 583.0650 3.55256712 Dec 4 11.210755 585.6329 14.15634374 Jan 5 -4.885406 588.2008 10.68455806 Feb 5 -1.823171 590.2921 6.53106077 Mar 5 1.305063 592.3834 -2.68843560 Apr 5 -3.197795 593.6110 -1.41321023 May 5 -12.600652 594.8386 1.76201503 Jun 5 -19.361145 595.3464 -2.98524089 Jul 5 -29.821637 595.8541 0.96750111 Aug 5 -27.783047 596.1342 0.64884976 Sep 5 25.455542 596.4143 -0.86980124 Oct 5 35.119009 596.5235 -2.64250287 Nov 5 26.382479 596.6327 4.98479290 Dec 5 11.210755 596.6606 4.12861258 Jan 6 -4.885406 596.6885 3.19686995 Feb 6 -1.823171 596.7145 2.10862149 Mar 6 1.305063 596.7406 -5.04562606 Apr 6 -3.197795 596.3392 -3.14137570 May 6 -12.600652 595.9378 -3.33712544 Jun 6 -19.361145 594.5178 -1.15668168 Jul 6 -29.821637 593.0979 9.72376000 Aug 6 -27.783047 590.5701 10.21290009 Sep 6 25.455542 588.0424 6.50204054 Oct 6 35.119009 584.4054 6.47562268 Nov 6 26.382479 580.7683 12.84920221 Dec 6 11.210755 576.0006 0.78862003 Jan 7 -4.885406 571.2329 -0.34752445 Feb 7 -1.823171 565.6599 -6.83672073 Mar 7 1.305063 560.0869 -0.39191611 Apr 7 -3.197795 554.4862 -2.28842213 May 7 -12.600652 548.8856 -4.28492826 Jun 7 -19.361145 543.8294 1.53170636 Jul 7 -29.821637 538.7733 2.04833890 Aug 7 -27.783047 534.6393 -7.85621437 Sep 7 25.455542 530.5052 -0.96076729 Oct 7 35.119009 527.0326 2.84840105 Nov 7 26.382479 523.5600 -7.94243322 Dec 7 11.210755 520.5765 -4.78721667 Jan 8 -4.885406 517.5930 -2.70756241 Feb 8 -1.823171 515.1324 0.69080723 Mar 8 1.305063 512.6718 3.02317778 Apr 8 -3.197795 510.5669 0.63087654 May 8 -12.600652 508.4621 -2.86142480 Jun 8 -19.361145 506.9225 2.43864214 Jul 8 -29.821637 505.3829 -6.56129300 Aug 8 -27.783047 505.1782 0.60482077 Sep 8 25.455542 504.9735 -2.42906511 Oct 8 35.119009 506.8818 -8.00081102 Nov 8 26.382479 508.7901 -17.17255954 Dec 8 11.210755 512.7306 -17.94140230 Jan 9 -4.885406 516.6712 -9.78580736 Feb 9 -1.823171 521.9589 -4.13570795 Mar 9 1.305063 527.2465 -0.55160764 Apr 9 -3.197795 532.8856 3.31220148 May 9 -12.600652 538.5246 10.07601049 Jun 9 -19.361145 543.5687 12.79247958 Jul 9 -29.821637 548.6127 5.20894660 Aug 9 -27.783047 552.6856 11.09746007 Sep 9 25.455542 556.7585 4.78597390 Oct 9 35.119009 559.7569 2.12410450 Nov 9 26.382479 562.7553 -8.13776752 Dec 9 11.210755 564.5743 -11.78506611 Jan 10 -4.885406 566.3933 -3.50792700 Feb 10 -1.823171 567.6258 9.19738910 Mar 10 1.305063 568.8582 9.83670611 Apr 10 -3.197795 569.1377 9.06005877 May 10 -12.600652 569.4172 6.18341131 Jun 10 -19.361145 568.6045 2.75663904 Jul 10 -29.821637 567.7918 -0.97013531 Aug 10 -27.783047 566.8219 5.96117395 Sep 10 25.455542 565.8520 9.69248356 Oct 10 35.119009 564.7040 4.17698327 Nov 10 26.382479 563.5560 -3.93851964 Dec 10 11.210755 562.1823 -9.39306587 Jan 11 -4.885406 560.8086 -6.92317439 > m$win s t l 1211 19 13 > m$deg s t l 0 1 1 > m$jump s t l 122 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/rcomp/tmp/1ofo21293480365.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/rcomp/tmp/2ofo21293480365.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/rcomp/tmp/3hp551293480365.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/rcomp/tmp/4hp551293480365.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/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/5o82y1293480365.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/6hh211293480365.tab") > > try(system("convert tmp/1ofo21293480365.ps tmp/1ofo21293480365.png",intern=TRUE)) character(0) > try(system("convert tmp/2ofo21293480365.ps tmp/2ofo21293480365.png",intern=TRUE)) character(0) > try(system("convert tmp/3hp551293480365.ps tmp/3hp551293480365.png",intern=TRUE)) character(0) > try(system("convert tmp/4hp551293480365.ps tmp/4hp551293480365.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.290 0.647 2.798