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(10570,10297,10635,10872,10296,10383,10431,10574,10653,10805,10872,10625,10407,10463,10556,10646,10702,11353,11346,11451,11964,12574,13031,13812,14544,14931,14886,16005,17064,15168,16050,15839,15137,14954,15648,15305,15579,16348,15928,16171,15937,15713,15594,15683,16438,17032,17696,17745,19394,20148,20108,18584,18441,18391,19178,18079,18483,19644,19195,19650,20830,23595,22937,21814,21928,21777,21383,21467,22052,22680,24320,24977,25204) > 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 409.29150 10186.64 -2.592700e+01 Feb 1 819.53997 10253.57 -7.761053e+02 Mar 1 515.67567 10320.50 -2.011708e+02 Apr 1 185.76809 10373.90 3.123367e+02 May 1 61.52703 10427.30 -1.928223e+02 Jun 1 -361.28984 10469.34 2.749527e+02 Jul 1 -320.77272 10511.38 2.403938e+02 Aug 1 -643.35760 10549.52 6.678341e+02 Sep 1 -545.44246 10587.67 6.107744e+02 Oct 1 -251.20143 10598.06 4.581429e+02 Nov 1 61.53946 10608.45 2.020117e+02 Dec 1 68.72160 10629.42 -7.313870e+01 Jan 2 409.29150 10650.39 -6.526768e+02 Feb 2 819.53997 10730.92 -1.087461e+03 Mar 2 515.67567 10811.46 -7.711320e+02 Apr 2 185.76809 10982.12 -5.218892e+02 May 2 61.52703 11152.79 -5.123129e+02 Jun 2 -361.28984 11436.77 2.775201e+02 Jul 2 -320.77272 11720.75 -5.398079e+01 Aug 2 -643.35760 12085.79 8.570464e+00 Sep 2 -545.44246 12450.82 5.862170e+01 Oct 2 -251.20143 12861.02 -3.582042e+01 Nov 2 61.53946 13271.22 -3.017624e+02 Dec 2 68.72160 13676.82 6.645573e+01 Jan 3 409.29150 14082.42 5.228611e+01 Feb 3 819.53997 14433.38 -3.219165e+02 Mar 3 515.67567 14784.33 -4.140064e+02 Apr 3 185.76809 15036.36 7.828680e+02 May 3 61.52703 15288.40 1.714076e+03 Jun 3 -361.28984 15436.54 9.274596e+01 Jul 3 -320.77272 15584.69 7.860821e+02 Aug 3 -643.35760 15642.54 8.398182e+02 Sep 3 -545.44246 15700.39 -1.794573e+01 Oct 3 -251.20143 15693.57 -4.883730e+02 Nov 3 61.53946 15686.76 -1.003001e+02 Dec 3 68.72160 15668.60 -4.323261e+02 Jan 4 409.29150 15650.45 -4.807399e+02 Feb 4 819.53997 15695.15 -1.666926e+02 Mar 4 515.67567 15739.86 -3.275325e+02 Apr 4 185.76809 15879.17 1.060624e+02 May 4 61.52703 16018.48 -1.430093e+02 Jun 4 -361.28984 16247.96 -1.736734e+02 Jul 4 -320.77272 16477.44 -5.626715e+02 Aug 4 -643.35760 16770.56 -4.442046e+02 Sep 4 -545.44246 17063.68 -8.023780e+01 Oct 4 -251.20143 17349.02 -6.582309e+01 Nov 4 61.53946 17634.37 9.175024e-02 Dec 4 68.72160 17887.67 -2.113901e+02 Jan 5 409.29150 18140.97 8.437402e+02 Feb 5 819.53997 18353.24 9.752217e+02 Mar 5 515.67567 18565.51 1.026816e+03 Apr 5 185.76809 18724.52 -3.262908e+02 May 5 61.52703 18883.54 -5.040642e+02 Jun 5 -361.28984 19020.55 -2.682642e+02 Jul 5 -320.77272 19157.57 3.412018e+02 Aug 5 -643.35760 19360.99 -6.386315e+02 Sep 5 -545.44246 19564.41 -5.359649e+02 Oct 5 -251.20143 19844.14 5.106144e+01 Nov 5 61.53946 20123.87 -9.904121e+02 Dec 5 68.72160 20415.95 -8.346707e+02 Jan 6 409.29150 20708.03 -2.873171e+02 Feb 6 819.53997 20989.76 1.785696e+03 Mar 6 515.67567 21271.50 1.149822e+03 Apr 6 185.76809 21578.16 5.007106e+01 May 6 61.52703 21884.82 -1.834683e+01 Jun 6 -361.28984 22214.00 -7.570871e+01 Jul 6 -320.77272 22543.18 -8.394046e+02 Aug 6 -643.35760 22861.94 -7.515857e+02 Sep 6 -545.44246 23180.71 -5.832668e+02 Oct 6 -251.20143 23496.17 -5.649729e+02 Nov 6 61.53946 23811.64 4.468212e+02 Dec 6 68.72160 24134.51 7.737662e+02 Jan 7 409.29150 24457.39 3.373235e+02 > m$win s t l 731 19 13 > m$deg s t l 0 1 1 > m$jump s t l 74 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/freestat/rcomp/tmp/1utsr1292971629.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/252au1292971629.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/352au1292971629.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/4gurf1292971629.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/5ulo61292971629.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/6f45u1292971629.tab") > > try(system("convert tmp/1utsr1292971629.ps tmp/1utsr1292971629.png",intern=TRUE)) character(0) > try(system("convert tmp/252au1292971629.ps tmp/252au1292971629.png",intern=TRUE)) character(0) > try(system("convert tmp/352au1292971629.ps tmp/352au1292971629.png",intern=TRUE)) character(0) > try(system("convert tmp/4gurf1292971629.ps tmp/4gurf1292971629.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.478 0.876 1.604