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(1775,2197,2920,4240,5415,6136,6719,6234,7152,3646,2165,2803,1615,2350,3350,3536,5834,6767,5993,7276,5641,3477,2247,2466,1567,2237,2598,3729,5715,5776,5852,6878,5488,3583,2054,2282,1552,2261,2446,3519,5161,5085,5711,6057,5224,3363,1899,2115,1491,2061,2419,3430,4778,4862,6176,5664,5529,3418,1941,2402,1579,2146,2462,3695,4831,5134,6250,5760,6249,2917,1741,2359,1511,2059,2635,2867,4403,5720,4502,5749,5627,2846,1762,2429,1169,2154,2249,2687,4359,5382,4459,6398,4596,3024,1887,2070,1351,2218,2461,3028,4784,4975,4607,6249,4809,3157,1910,2228,1594,2467,2222,3607,4685,4962,5770,5480,5000,3228,1993,2288,1588,2105,2191,3591,4668,4885,5822,5599,5340,3082,2010,2301) > 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 -2239.8282 4175.296 -160.467361 Feb 1 -1555.3589 4196.608 -444.249273 Mar 1 -1213.1624 4217.921 -84.758367 Apr 1 -301.5167 4232.529 308.988026 May 1 1221.7657 4247.137 -53.902263 Jun 1 1682.6615 4257.966 195.372631 Jul 1 1882.2844 4268.795 567.920437 Aug 1 2381.9805 4279.586 -427.566841 Sep 1 1775.1307 4290.377 1086.491837 Oct 1 -485.3838 4288.163 -156.779159 Nov 1 -1765.7166 4285.948 -355.231852 Dec 1 -1382.8563 4285.553 -99.696334 Jan 2 -2239.8282 4285.157 -430.328596 Feb 2 -1555.3589 4283.519 -378.160223 Mar 2 -1213.1624 4281.881 281.280968 Apr 2 -301.5167 4274.257 -436.739884 May 2 1221.7657 4266.632 345.602584 Jun 2 1682.6615 4254.979 829.359741 Jul 2 1882.2844 4243.326 -132.610190 Aug 2 2381.9805 4213.319 680.700896 Sep 2 1775.1307 4183.311 -317.442063 Oct 2 -485.3838 4143.735 -181.351100 Nov 2 -1765.7166 4104.158 -91.441834 Dec 2 -1382.8563 4070.234 -221.377470 Jan 3 -2239.8282 4036.309 -229.480885 Feb 3 -1555.3589 4023.678 -231.319593 Mar 3 -1213.1624 4011.048 -199.885484 Apr 3 -301.5167 4011.654 18.862497 May 3 1221.7657 4012.260 480.973797 Jun 3 1682.6615 4008.878 84.460033 Jul 3 1882.2844 4005.496 -35.780819 Aug 3 2381.9805 3990.068 505.951820 Sep 3 1775.1307 3974.639 -261.769586 Oct 3 -485.3838 3943.158 125.226096 Nov 3 -1765.7166 3911.677 -91.959918 Dec 3 -1382.8563 3871.524 -206.668128 Jan 4 -2239.8282 3831.372 -39.544117 Feb 4 -1555.3589 3797.632 18.727345 Mar 4 -1213.1624 3763.891 -104.728375 Apr 4 -301.5167 3742.899 77.617264 May 4 1221.7657 3721.908 217.326221 Jun 4 1682.6615 3708.837 -306.498906 Jul 4 1882.2844 3695.767 132.948880 Aug 4 2381.9805 3684.659 -9.639174 Sep 4 1775.1307 3673.551 -224.681272 Oct 4 -485.3838 3662.019 186.364494 Nov 4 -1765.7166 3650.488 14.228563 Dec 4 -1382.8563 3644.050 -146.194010 Jan 5 -2239.8282 3637.613 93.215638 Feb 5 -1555.3589 3637.957 -21.598597 Mar 5 -1213.1624 3638.302 -6.140014 Apr 5 -301.5167 3645.920 85.596493 May 5 1221.7657 3653.538 -97.303682 Jun 5 1682.6615 3665.092 -485.753489 Jul 5 1882.2844 3676.646 617.069615 Aug 5 2381.9805 3689.733 -407.713820 Sep 5 1775.1307 3702.821 51.048699 Oct 5 -485.3838 3715.574 187.810243 Nov 5 -1765.7166 3728.327 -21.609910 Dec 5 -1382.8563 3740.927 43.929033 Jan 6 -2239.8282 3753.528 65.300197 Feb 6 -1555.3589 3764.689 -63.329789 Mar 6 -1213.1624 3775.849 -100.686957 Apr 6 -301.5167 3777.450 219.066466 May 6 1221.7657 3779.051 -169.816793 Jun 6 1682.6615 3772.238 -320.899881 Jul 6 1882.2844 3765.426 602.289942 Aug 6 2381.9805 3755.492 -377.472777 Sep 6 1775.1307 3745.559 728.310460 Oct 6 -485.3838 3723.452 -321.067802 Nov 6 -1765.7166 3701.344 -194.627760 Dec 6 -1382.8563 3662.126 79.730204 Jan 7 -2239.8282 3622.908 127.920389 Feb 7 -1555.3589 3581.147 33.211966 Mar 7 -1213.1624 3539.386 308.776361 Apr 7 -301.5167 3515.903 -347.386096 May 7 1221.7657 3492.420 -311.185234 Jun 7 1682.6615 3484.563 552.775461 Jul 7 1882.2844 3476.707 -856.990933 Aug 7 2381.9805 3471.579 -104.559621 Sep 7 1775.1307 3466.452 385.417646 Oct 7 -485.3838 3460.963 -129.579658 Nov 7 -1765.7166 3455.475 72.241340 Dec 7 -1382.8563 3446.140 365.716017 Jan 8 -2239.8282 3436.805 -27.977087 Feb 8 -1555.3589 3422.545 286.813691 Mar 8 -1213.1624 3408.285 53.877287 Apr 8 -301.5167 3391.001 -402.484707 May 8 1221.7657 3373.718 -236.483382 Jun 8 1682.6615 3369.856 329.482495 Jul 8 1882.2844 3365.994 -789.278717 Aug 8 2381.9805 3382.358 633.661811 Sep 8 1775.1307 3398.721 -577.851706 Oct 8 -485.3838 3423.709 85.674836 Nov 8 -1765.7166 3448.697 204.019681 Dec 8 -1382.8563 3458.310 -5.453241 Jan 9 -2239.8282 3467.922 122.906058 Feb 9 -1555.3589 3465.603 307.756100 Mar 9 -1213.1624 3463.283 210.878961 Apr 9 -301.5167 3460.921 -131.404715 May 9 1221.7657 3458.559 103.674928 Jun 9 1682.6615 3463.559 -171.220154 Jul 9 1882.2844 3468.558 -743.842323 Aug 9 2381.9805 3485.057 381.962268 Sep 9 1775.1307 3501.556 -467.687186 Oct 9 -485.3838 3526.523 115.861006 Nov 9 -1765.7166 3551.489 124.227502 Dec 9 -1382.8563 3575.384 35.472273 Jan 10 -2239.8282 3599.279 234.549264 Feb 10 -1555.3589 3604.648 417.710579 Mar 10 -1213.1624 3610.018 -174.855288 Apr 10 -301.5167 3602.700 305.816878 May 10 1221.7657 3595.382 -132.147636 Jun 10 1682.6615 3587.807 -308.468341 Jul 10 1882.2844 3580.232 307.483866 Aug 10 2381.9805 3576.106 -478.086255 Sep 10 1775.1307 3571.980 -347.110422 Oct 10 -485.3838 3573.950 139.433700 Nov 10 -1765.7166 3575.921 182.796124 Dec 10 -1382.8563 3582.273 88.582855 Jan 11 -2239.8282 3588.626 239.201808 Feb 11 -1555.3589 3591.928 68.430568 Mar 11 -1213.1624 3595.230 -191.067854 Apr 11 -301.5167 3593.787 298.729355 May 11 1221.7657 3592.344 -146.110117 Jun 11 1682.6615 3592.483 -390.144629 Jul 11 1882.2844 3592.622 347.093771 Aug 11 2381.9805 3591.991 -374.971937 Sep 11 1775.1307 3591.361 -26.491689 Oct 11 -485.3838 3591.487 -24.103263 Nov 11 -1765.7166 3591.613 184.103465 Dec 11 -1382.8563 3593.339 90.516923 > m$win s t l 1321 19 13 > m$deg s t l 0 1 1 > m$jump s t l 133 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/freestat/rcomp/tmp/1rt1y1293536289.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/2rt1y1293536289.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/32k1j1293536289.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/42k1j1293536289.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/5q3xu1293536289.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/6jvxx1293536289.tab") > > try(system("convert tmp/1rt1y1293536289.ps tmp/1rt1y1293536289.png",intern=TRUE)) character(0) > try(system("convert tmp/2rt1y1293536289.ps tmp/2rt1y1293536289.png",intern=TRUE)) character(0) > try(system("convert tmp/32k1j1293536289.ps tmp/32k1j1293536289.png",intern=TRUE)) character(0) > try(system("convert tmp/42k1j1293536289.ps tmp/42k1j1293536289.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.818 0.907 2.038