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(9700,9081,9084,9743,8587,9731,9563,9998,9437,10038,9918,9252,9737,9035,9133,9487,8700,9627,8947,9283,8829,9947,9628,9318,9605,8640,9214,9567,8547,9185,9470,9123,9278,10170,9434,9655,9429,8739,9552,9687,9019,9672,9206,9069,9788,10312,10105,9863,9656,9295,9946,9701,9049,10190,9706,9765,9893,9994,10433,10073,10112,9266,9820,10097,9115,10411,9678,10408,10153,10368,10581,10597,10680,9738,9556)
> par8 = 'FALSE'
> par7 = '1'
> par6 = ''
> par5 = '1'
> par4 = ''
> par3 = '0'
> par2 = '6'
> 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  194.11396  9557.096  -51.210362
Feb 1 -515.86657  9549.139   47.727321
Mar 1 -156.08156  9541.182 -301.100524
Apr 1  198.58128  9533.225   11.193804
May 1 -730.20129  9522.906 -205.704727
Jun 1  161.03909  9512.587   57.373796
Jul 1  -96.19507  9502.268  156.926854
Aug 1   25.17304  9491.152  481.674483
Sep 1 -128.92519  9480.037   85.888444
Oct 1  660.85193  9468.921  -91.772932
Nov 1  337.05131  9457.377  123.572027
Dec 1   53.84131  9445.832 -247.673642
Jan 2  184.56515  9434.288  118.146859
Feb 2 -522.31619  9409.629  147.686996
Mar 2 -132.66236  9384.970 -119.308030
Apr 2  196.45726  9360.312  -69.768856
May 2 -715.16519  9341.720   73.445511
Jun 2  155.58277  9323.128  148.289471
Jul 2 -103.06854  9304.536 -254.467291
Aug 2  -17.02634  9294.827    5.199560
Sep 2 -109.83428  9285.118 -346.283440
Oct 2  658.73545  9275.409   12.855895
Nov 2  337.26585  9273.489   17.244771
Dec 2   75.08527  9271.570  -28.655361
Jan 3  173.60802  9269.651  161.741169
Feb 3 -530.13857  9276.465 -106.326210
Mar 3  -85.23625  9283.279   15.957513
Apr 3  188.36877  9290.093   88.538528
May 3 -662.94331  9301.313  -91.369577
Jun 3  149.25559  9312.533 -276.788651
Jul 3 -132.66478  9323.753  278.911542
Aug 3 -141.06881  9337.181  -73.112411
Sep 3  -81.16672  9350.609    8.557523
Oct 3  656.60145  9364.037  149.361374
Nov 3  330.47773  9380.206 -276.683578
Dec 3  137.63529  9396.375  120.990185
Jan 4  133.56915  9412.543 -117.112342
Feb 4 -572.23047  9433.165 -121.934681
Mar 4   27.67098  9453.787   70.541918
Apr 4  144.54344  9474.409   68.047499
May 4 -667.19778  9499.157  187.041266
Jun 4  184.69169  9523.904  -36.595647
Jul 4 -155.63406  9548.651 -187.017353
Aug 4 -178.96780  9575.617 -327.649678
Sep 4   12.27432  9602.584  173.142136
Oct 4  533.83602  9629.550  148.614370
Nov 4  369.42116  9654.124   81.455255
Dec 4  216.64286  9678.698  -32.340431
Jan 5  178.16700  9703.272 -225.438544
Feb 5 -568.06619  9724.217  138.849634
Mar 5  -81.23401  9745.162  282.072446
Apr 5  122.19804  9766.107 -187.304611
May 5 -706.39316  9782.350  -26.956551
Jun 5  233.23923  9798.593  158.167926
Jul 5 -199.42395  9814.836   90.587971
Aug 5  -96.69047  9829.622   32.068094
Sep 5   40.41327  9844.409    8.177964
Oct 5  451.52713  9859.195 -316.722291
Nov 5  403.08167  9877.293  152.625795
Dec 5  254.66536  9895.390  -77.055265
Jan 6  195.15825  9913.487    3.354468
Feb 6 -561.76235  9938.029 -110.266156
Mar 6 -129.59822  9962.570  -12.971515
Apr 6  107.72284  9987.111    2.166199
May 6 -719.53700 10009.689 -175.152208
Jun 6  253.18095 10032.267  125.551593
Jul 6 -214.67858 10054.846 -162.167126
Aug 6  -72.19758 10075.919  404.278703
Sep 6   50.95674 10096.992    5.051219
Oct 6  416.16835 10118.065 -166.233566
Nov 6  420.27056 10139.421   21.308147
Dec 6  264.26480 10160.777  171.957839
Jan 7  202.24222 10182.133  295.624340
Feb 7 -558.77230 10203.253   93.518855
Mar 7 -149.85299 10224.373 -518.520463
> m$win
 s  t  l 
 6 25 13 
> m$deg
s t l 
0 1 1 
> m$jump
s t l 
1 3 2 
> m$inner
[1] 2
> m$outer
[1] 0
> postscript(file="/var/www/html/freestat/rcomp/tmp/1mike1291064573.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/2xajh1291064573.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/3xajh1291064573.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/4qjik1291064573.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/54byb1291064573.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/6pbfz1291064573.tab") 
> postscript(file="/var/www/html/freestat/rcomp/tmp/7pbfz1291064573.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> myresid <- m$time.series[!is.na(m$time.series[,'remainder']),'remainder']
> hist(as.numeric(myresid), main='Residual Histogram', xlab='Residual Value')
> dev.off()
null device 
          1 
> try(system("convert tmp/1mike1291064573.ps tmp/1mike1291064573.png",intern=TRUE))
character(0)
> try(system("convert tmp/2xajh1291064573.ps tmp/2xajh1291064573.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xajh1291064573.ps tmp/3xajh1291064573.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qjik1291064573.ps tmp/4qjik1291064573.png",intern=TRUE))
character(0)
> try(system("convert tmp/7pbfz1291064573.ps tmp/7pbfz1291064573.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.683   1.081   1.885