R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(27951,29781,32914,33488,35652,36488,35387,35676,34844,32447,31068,29010,29812,30951,32974,32936,34012,32946,31948,30599,27691,25073,23406,22248,22896,25317,26558,26471,27543,26198,24725,25005,23462,20780,19815,19761,21454,23899,24939,23580,24562,24696,23785,23812,21917,19713,19282,18788,21453,24482,27474,27264,27349,30632,29429,30084,26290,24379,23335,21346,21106,24514,28353,30805,31348,34556,33855,34787,32529,29998,29257,28155,30466,35704,39327,39351,42234,43630,43722,43121,37985,37135,34646,33026,35087,38846,42013,43908,42868,44423,44167,43636,44382,42142,43452,36912,42413,45344,44873,47510,49554,47369,45998,48140,48441,44928,40454,38661,37246,36843,36424,37594,38144,38737,34560,36080,33508,35462,33374,32110,35533,35532,37903,36763,40399,44164,44496,43110,43880,43930,44327) > 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 -3271.3448 31678.69 -456.34678 Feb 1 -1022.0822 31888.14 -1085.05769 Mar 1 931.0883 32097.59 -114.67646 Apr 1 1379.5985 32281.05 -172.64557 May 1 2562.3802 32464.51 625.11391 Jun 1 3400.3664 32624.43 463.19957 Jul 1 2243.7162 32784.36 358.92157 Aug 1 2316.1628 32927.18 432.65239 Sep 1 470.5170 33070.01 1303.47557 Oct 1 -1384.9630 33067.90 764.06492 Nov 1 -2752.1696 33065.79 754.38095 Dec 1 -4873.2694 32813.57 1069.70279 Jan 2 -3271.3448 32561.34 522.00013 Feb 2 -1022.0822 32115.69 -142.61216 Mar 2 931.0883 31670.04 372.86770 Apr 2 1379.5985 31076.93 479.47024 May 2 2562.3802 30483.82 965.80136 Jun 2 3400.3664 29875.06 -329.42879 Jul 2 2243.7162 29266.31 437.97741 Aug 2 2316.1628 28706.54 -423.70739 Sep 2 470.5170 28146.78 -926.29984 Oct 2 -1384.9630 27620.84 -1162.87461 Nov 2 -2752.1696 27094.89 -936.72270 Dec 2 -4873.2694 26586.09 535.18140 Jan 3 -3271.3448 26077.28 90.06097 Feb 3 -1022.0822 25623.71 715.37305 Mar 3 931.0883 25170.13 456.77728 Apr 3 1379.5985 24792.45 298.95044 May 3 2562.3802 24414.77 565.85218 Jun 3 3400.3664 24149.33 -1351.69576 Jul 3 2243.7162 23883.89 -1402.60736 Aug 3 2316.1628 23730.19 -1041.35428 Sep 3 470.5170 23576.49 -585.00884 Oct 3 -1384.9630 23465.54 -1300.57542 Nov 3 -2752.1696 23354.58 -787.41533 Dec 3 -4873.2694 23250.24 1384.02984 Jan 4 -3271.3448 23145.89 1579.45049 Feb 4 -1022.0822 23029.03 1892.05026 Mar 4 931.0883 22912.17 1095.74218 Apr 4 1379.5985 22761.76 -561.35958 May 4 2562.3802 22611.35 -611.73277 Jun 4 3400.3664 22511.00 -1215.36342 Jul 4 2243.7162 22410.64 -869.35774 Aug 4 2316.1628 22498.31 -1002.47414 Sep 4 470.5170 22585.98 -1139.49818 Oct 4 -1384.9630 22890.46 -1792.49915 Nov 4 -2752.1696 23194.94 -1160.77343 Dec 4 -4873.2694 23646.12 15.14884 Jan 5 -3271.3448 24097.30 627.04660 Feb 5 -1022.0822 24576.99 927.09595 Mar 5 931.0883 25056.67 1486.23745 Apr 5 1379.5985 25428.78 455.61804 May 5 2562.3802 25800.89 -1014.27278 Jun 5 3400.3664 25964.67 1266.96528 Jul 5 2243.7162 26128.44 1056.83969 Aug 5 2316.1628 26186.58 1581.25900 Sep 5 470.5170 26244.71 -425.22933 Oct 5 -1384.9630 26390.08 -626.11381 Nov 5 -2752.1696 26535.44 -448.27160 Dec 5 -4873.2694 26843.65 -624.38455 Jan 6 -3271.3448 27151.87 -2774.52201 Feb 6 -1022.0822 27606.54 -2070.45433 Mar 6 931.0883 28061.21 -639.29448 Apr 6 1379.5985 28633.36 792.03767 May 6 2562.3802 29205.52 -419.90158 Jun 6 3400.3664 29862.55 1293.07912 Jul 6 2243.7162 30519.59 1091.69616 Aug 6 2316.1628 31277.75 1193.08562 Sep 6 470.5170 32035.92 22.56743 Oct 6 -1384.9630 32824.24 -1441.27698 Nov 6 -2752.1696 33612.56 -1603.39471 Dec 6 -4873.2694 34422.82 -1394.54571 Jan 7 -3271.3448 35233.07 -1495.72122 Feb 7 -1022.0822 35987.64 738.43894 Mar 7 931.0883 36742.22 1653.69126 Apr 7 1379.5985 37342.34 629.06573 May 7 2562.3802 37942.45 1729.16878 Jun 7 3400.3664 38318.91 1910.72306 Jul 7 2243.7162 38695.37 2782.91370 Aug 7 2316.1628 38911.30 1893.53372 Sep 7 470.5170 39127.24 -1612.75390 Oct 7 -1384.9630 39266.13 -746.16590 Nov 7 -2752.1696 39405.02 -2006.85123 Dec 7 -4873.2694 39534.25 -1634.97874 Jan 8 -3271.3448 39663.48 -1305.13078 Feb 8 -1022.0822 39925.34 -57.25962 Mar 8 931.0883 40187.21 894.70370 Apr 8 1379.5985 40656.71 1871.69108 May 8 2562.3802 41126.21 -820.59295 Jun 8 3400.3664 41633.28 -610.64909 Jul 8 2243.7162 42140.35 -217.06887 Aug 8 2316.1628 42582.27 -1262.43305 Sep 8 470.5170 43024.19 887.29514 Oct 8 -1384.9630 43422.66 104.29869 Nov 8 -2752.1696 43821.14 2383.02893 Dec 8 -4873.2694 44154.43 -2369.15604 Jan 9 -3271.3448 44487.71 1196.63448 Feb 9 -1022.0822 44743.88 1622.19972 Mar 9 931.0883 45000.05 -1058.14288 Apr 9 1379.5985 45130.65 999.74775 May 9 2562.3802 45261.25 1730.36696 Jun 9 3400.3664 45141.64 -1173.00872 Jul 9 2243.7162 45022.03 -1267.74806 Aug 9 2316.1628 44512.14 1311.69345 Sep 9 470.5170 44002.26 3968.22732 Oct 9 -1384.9630 43204.55 3108.41148 Nov 9 -2752.1696 42406.85 799.32231 Dec 9 -4873.2694 41422.36 2111.90858 Jan 10 -3271.3448 40437.87 79.47033 Feb 10 -1022.0822 39364.58 -1499.49401 Mar 10 931.0883 38291.28 -2798.36619 Apr 10 1379.5985 37400.93 -1186.53050 May 10 2562.3802 36510.59 -928.96622 Jun 10 3400.3664 36082.67 -746.03603 Jul 10 2243.7162 35654.75 -3338.46950 Aug 10 2316.1628 35599.87 -1836.02924 Sep 10 470.5170 35544.98 -2507.49662 Oct 10 -1384.9630 35719.30 1127.66020 Nov 10 -2752.1696 35893.63 232.54369 Dec 10 -4873.2694 36325.07 658.20072 Jan 11 -3271.3448 36756.51 2047.83323 Feb 11 -1022.0822 37454.14 -900.05677 Mar 11 931.0883 38151.77 -1179.85461 Apr 11 1379.5985 38988.93 -3605.52663 May 11 2562.3802 39826.09 -1989.47006 Jun 11 3400.3664 40671.77 91.86652 Jul 11 2243.7162 41517.44 734.83944 Aug 11 2316.1628 42390.59 -1596.75119 Sep 11 470.5170 43263.73 145.75053 Oct 11 -1384.9630 44175.75 1139.21563 Nov 11 -2752.1696 45087.76 1991.40740 > m$win s t l 1311 19 13 > m$deg s t l 0 1 1 > m$jump s t l 132 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/rcomp/tmp/1feke1292703215.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/rcomp/tmp/2feke1292703215.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/rcomp/tmp/3q51g1292703215.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/rcomp/tmp/40wi11292703215.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/5f6ga1292703215.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/rcomp/tmp/6pgyv1292703215.tab") > > try(system("convert tmp/1feke1292703215.ps tmp/1feke1292703215.png",intern=TRUE)) character(0) > try(system("convert tmp/2feke1292703215.ps tmp/2feke1292703215.png",intern=TRUE)) character(0) > try(system("convert tmp/3q51g1292703215.ps tmp/3q51g1292703215.png",intern=TRUE)) character(0) > try(system("convert tmp/40wi11292703215.ps tmp/40wi11292703215.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.319 0.675 3.149