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(3440,2678,2981,2260,2844,2546,2456,2295,2379,2479,2057,2280,2351,2276,2548,2311,2201,2725,2408,2139,1898,2539,2070,2063,2565,2442,2194,2798,2074,2628,2289,2154,2467,2137,1850,2075,1791,1755,2232,1952,1822,2522,2074,2366,2173,2094,1833,1858,2040,2133,2921,3252,3318,3554,2308,1621,1315,1501,1418,1657) > 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 50.76456 2800.459 588.776669 Feb 1 -112.86959 2751.922 38.947667 Mar 1 222.49589 2703.385 55.119045 Apr 1 185.05107 2658.058 -583.108578 May 1 145.40567 2612.730 85.864383 Jun 1 505.60712 2567.152 -526.759467 Jul 1 34.60842 2521.575 -100.183157 Aug 1 -143.26982 2477.364 -39.094227 Sep 1 -197.74816 2433.153 143.594810 Oct 1 -84.97671 2412.768 151.208677 Nov 1 -380.20521 2392.383 44.822495 Dec 1 -224.86335 2375.940 128.923540 Jan 2 50.76456 2359.497 -59.261470 Feb 2 -112.86959 2340.136 48.733650 Mar 2 222.49589 2320.775 4.729151 Apr 2 185.05107 2308.352 -182.402620 May 2 145.40567 2295.928 -240.333809 Jun 2 505.60712 2297.496 -78.102848 Jul 2 34.60842 2299.063 74.328272 Aug 2 -143.26982 2309.209 -26.939253 Sep 2 -197.74816 2319.355 -223.606672 Oct 2 -84.97671 2326.641 297.335341 Nov 2 -380.20521 2333.928 116.277305 Dec 2 -224.86335 2331.590 -43.726377 Jan 3 50.76456 2329.252 184.983885 Feb 3 -112.86959 2325.160 229.709727 Mar 3 222.49589 2321.068 -349.564051 Apr 3 185.05107 2313.931 299.017775 May 3 145.40567 2306.794 -378.199816 Jun 3 505.60712 2285.664 -163.271005 Jul 3 34.60842 2264.534 -10.142033 Aug 3 -143.26982 2232.893 64.376656 Sep 3 -197.74816 2201.253 463.495452 Oct 3 -84.97671 2167.006 54.970493 Nov 3 -380.20521 2132.760 97.445486 Dec 3 -224.86335 2103.387 196.476784 Jan 4 50.76456 2074.013 -333.777974 Feb 4 -112.86959 2057.577 -189.707055 Mar 4 222.49589 2041.140 -31.635756 Apr 4 185.05107 2040.662 -273.712836 May 4 145.40567 2040.184 -363.589334 Jun 4 505.60712 2054.306 -37.913507 Jul 4 34.60842 2068.429 -29.037519 Aug 4 -143.26982 2108.742 400.528150 Sep 4 -197.74816 2149.054 221.693925 Oct 4 -84.97671 2223.220 -44.243780 Nov 4 -380.20521 2297.387 -84.181535 Dec 4 -224.86335 2361.068 -278.204729 Jan 5 50.76456 2424.749 -435.513978 Feb 5 -112.86959 2425.841 -179.971326 Mar 5 222.49589 2426.932 271.571707 Apr 5 185.05107 2359.467 707.481616 May 5 145.40567 2292.002 880.592109 Jun 5 505.60712 2234.846 813.546720 Jul 5 34.60842 2177.690 95.701491 Aug 5 -143.26982 2112.921 -348.650943 Sep 5 -197.74816 2048.151 -535.403271 Oct 5 -84.97671 1970.149 -384.171975 Nov 5 -380.20521 1892.146 -93.940728 Dec 5 -224.86335 1805.719 76.143909 > 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/wessaorg/rcomp/tmp/1bfxz1322577404.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/2pq261322577404.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/3qnit1322577404.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/4bbe51322577404.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/5xl0e1322577404.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/60xol1322577404.tab") > > try(system("convert tmp/1bfxz1322577404.ps tmp/1bfxz1322577404.png",intern=TRUE)) character(0) > try(system("convert tmp/2pq261322577404.ps tmp/2pq261322577404.png",intern=TRUE)) character(0) > try(system("convert tmp/3qnit1322577404.ps tmp/3qnit1322577404.png",intern=TRUE)) character(0) > try(system("convert tmp/4bbe51322577404.ps tmp/4bbe51322577404.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.272 0.205 1.495