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(548,563,581,572,519,521,531,540,548,556,551,549,564,586,604,601,545,537,552,563,575,580,575,558,564,581,597,587,536,524,537,536,533,528,516,502,506,518,534,528,478,469,490,493,508,517,514,510,527,542,565,555,499,511,526,532,549,561,557,566,588,620,626,620,573,573,574,580,590,593,597,595,612,628,629,621,569,567,573,584,589,591,595,594,611,613,611,594,543,537,544,555,561,562,555,547,565,578,580,569,507,501,509,510,517,519,512,509,519) > 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 8.472976 539.7053 -0.17831570 Feb 1 24.448758 541.2921 -2.74084670 Mar 1 35.464288 542.8788 2.65687512 Apr 1 26.758264 544.4787 0.76307111 May 1 -26.169979 546.0785 -0.90851248 Jun 1 -29.226942 547.7132 2.51371247 Jul 1 -18.395017 549.3480 0.04704865 Aug 1 -11.813731 551.0276 0.78610871 Sep 1 -3.010221 552.7073 -1.69705600 Oct 1 1.510182 554.5266 -0.03679920 Nov 1 -1.969413 556.3460 -3.37654397 Dec 1 -6.069163 558.3033 -3.23410320 Jan 2 8.472976 560.2606 -4.73355174 Feb 2 24.448758 562.2708 -0.71956387 Mar 2 35.464288 564.2810 4.25467684 Apr 2 26.758264 566.2387 8.00302208 May 2 -26.169979 568.1964 2.97358774 Jun 2 -29.226942 569.2283 -3.00132556 Jul 2 -18.395017 570.2601 0.13487240 Aug 2 -11.813731 570.0164 4.79728639 Sep 2 -3.010221 569.7727 8.23747562 Oct 2 1.510182 568.7884 9.70146089 Nov 2 -1.969413 567.8040 9.16544459 Dec 2 -6.069163 566.3725 -2.30334604 Jan 3 8.472976 564.9411 -9.41402597 Feb 3 24.448758 562.5174 -5.96616352 Mar 3 35.464288 560.0938 1.44195177 Apr 3 26.758264 556.4709 3.77079061 May 3 -26.169979 552.8481 9.32184987 Jun 3 -29.226942 548.2481 4.97880721 Jul 3 -18.395017 543.6481 11.74687579 Aug 3 -11.813731 538.2773 9.53641563 Sep 3 -3.010221 532.9065 3.10373070 Oct 3 1.510182 527.3669 -0.87706613 Nov 3 -1.969413 521.8273 -3.85786453 Dec 3 -6.069163 517.1227 -9.05353637 Jan 4 8.472976 512.4181 -14.89109752 Feb 4 24.448758 509.3822 -15.83097167 Mar 4 35.464288 506.3463 -7.81059301 Apr 4 26.758264 505.4943 -4.25256743 May 4 -26.169979 504.6423 -0.47232144 Jun 4 -29.226942 505.5945 -7.36754407 Jul 4 -18.395017 506.5467 1.84834454 Aug 4 -11.813731 508.5111 -3.69735308 Sep 4 -3.010221 510.4755 0.53472454 Oct 4 1.510182 512.7897 2.70014092 Nov 4 -1.969413 515.1039 0.86555574 Dec 4 -6.069163 517.8393 -1.77013137 Jan 5 8.472976 520.5747 -2.04770778 Feb 5 24.448758 523.7546 -6.20338950 Mar 5 35.464288 526.9345 2.60118160 Apr 5 26.758264 530.6397 -2.39791461 May 5 -26.169979 534.3448 -9.17479041 Jun 5 -29.226942 538.9106 1.31630486 Jul 5 -18.395017 543.4765 0.91851138 Aug 5 -11.813731 548.9321 -5.11832193 Sep 5 -3.010221 554.3876 -2.37738000 Oct 5 1.510182 560.1245 -0.63469790 Nov 5 -1.969413 565.8614 -6.89201737 Dec 5 -6.069163 571.1302 0.93892909 Jan 6 8.472976 576.3990 3.12798625 Feb 6 24.448758 580.5634 14.98787047 Mar 6 35.464288 584.7277 5.80800750 Apr 6 26.758264 587.7360 5.50570782 May 6 -26.169979 590.7444 8.42562854 Jun 6 -29.226942 592.7123 9.51465806 Jul 6 -18.395017 594.6802 -2.28520119 Aug 6 -11.813731 595.5859 -3.77220327 Sep 6 -3.010221 596.4917 -3.48143011 Oct 6 1.510182 596.6793 -5.18947714 Nov 6 -1.969413 596.8669 2.10247426 Dec 6 -6.069163 596.8103 4.25884770 Jan 7 8.472976 596.7537 6.77333182 Feb 7 24.448758 596.6894 6.86182453 Mar 7 35.464288 596.6251 -3.08942995 Apr 7 26.758264 596.3437 -2.10195863 May 7 -26.169979 596.0622 -0.89226689 Jun 7 -29.226942 595.6925 0.53444228 Jul 7 -18.395017 595.3228 -3.92773731 Aug 7 -11.813731 594.5620 1.25170699 Sep 7 -3.010221 593.8013 -1.79107348 Oct 7 1.510182 592.1953 -2.70543399 Nov 7 -1.969413 590.5892 6.38020393 Dec 7 -6.069163 588.2006 11.86853448 Jan 8 8.472976 585.8120 16.71497571 Feb 8 24.448758 583.0720 5.47925202 Mar 8 35.464288 580.3319 -4.79621885 Apr 8 26.758264 577.1671 -9.92536563 May 8 -26.169979 574.0023 -4.83229200 Jun 8 -29.226942 570.6495 -4.42253099 Jul 8 -18.395017 567.2967 -4.90165873 Aug 8 -11.813731 564.4576 2.35611923 Sep 8 -3.010221 561.6185 2.39167242 Oct 8 1.510182 559.1277 1.36216465 Nov 8 -1.969413 556.6368 0.33265531 Dec 8 -6.069163 553.6820 -0.61279741 Jan 9 8.472976 550.7272 5.79986056 Feb 9 24.448758 547.2434 6.30786908 Mar 9 35.464288 543.7596 0.77613042 Apr 9 26.758264 540.0587 2.18299610 May 9 -26.169979 536.3579 -3.18791780 Jun 9 -29.226942 532.7638 -2.53684398 Jul 9 -18.395017 529.1697 -1.77465891 Aug 9 -11.813731 525.5563 -3.74261354 Sep 9 -3.010221 521.9430 -1.93279293 Oct 9 1.510182 518.3424 -0.85260287 Nov 9 -1.969413 514.7418 -0.77241439 Dec 9 -6.069163 511.1947 3.87444467 Jan 10 8.472976 507.6476 2.87941441 > m$win s t l 1091 19 13 > m$deg s t l 0 1 1 > m$jump s t l 110 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/wessaorg/rcomp/tmp/1s84i1322768048.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/2tk9j1322768048.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/3gyss1322768048.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/4dzqb1322768048.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/53fbm1322768048.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/60oki1322768048.tab") > > try(system("convert tmp/1s84i1322768048.ps tmp/1s84i1322768048.png",intern=TRUE)) character(0) > try(system("convert tmp/2tk9j1322768048.ps tmp/2tk9j1322768048.png",intern=TRUE)) character(0) > try(system("convert tmp/3gyss1322768048.ps tmp/3gyss1322768048.png",intern=TRUE)) character(0) > try(system("convert tmp/4dzqb1322768048.ps tmp/4dzqb1322768048.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.526 0.206 1.728