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(220206,220115,218444,214912,210705,209673,237041,242081,241878,242621,238545,240337,244752,244576,241572,240541,236089,236997,264579,270349,269645,267037,258113,262813,267413,267366,264777,258863,254844,254868,277267,285351,286602,283042,276687,277915,277128,277103,275037,270150,267140,264993,287259,291186,292300,288186,281477,282656,280190,280408,276836,275216,274352,271311,289802,290726,292300,278506,269826,265861,269034,264176,255198,253353,246057,235372,258556,260993,254663,250643,243422,247105,248541,245039,237080,237085,225554,226839,247934,248333,246969,245098,246263,255765,264319,268347,273046,273963,267430,271993,292710,295881,293299,288576) > 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 2191.2067 216081.4 1933.36177 Feb 1 869.8936 218219.0 1026.11875 Mar 1 -3037.2958 220356.5 1124.75205 Apr 1 -5955.9043 222507.2 -1639.25833 May 1 -11875.1351 224657.8 -2077.64648 Jun 1 -13869.4003 226791.9 -3249.49696 Jul 1 8289.7110 228926.0 -174.72395 Aug 1 11142.3402 231038.5 -99.84565 Sep 1 9370.4705 233151.0 -643.46839 Oct 1 4100.9753 235317.5 3202.55310 Nov 1 -1510.3912 237483.9 2571.44583 Dec 1 283.5261 239650.4 403.07624 Jan 2 2191.2067 241816.8 743.94342 Feb 2 869.8936 243977.5 -271.37295 Mar 2 -3037.2958 246138.1 -1528.81299 Apr 2 -5955.9043 248196.0 -1699.13525 May 2 -11875.1351 250254.0 -2289.83527 Jun 2 -13869.4003 252208.9 -1342.51429 Jul 2 8289.7110 254163.9 2125.43021 Aug 2 11142.3402 256062.9 3143.80859 Sep 2 9370.4705 257961.8 2312.68594 Oct 2 4100.9753 259624.6 3311.42369 Nov 2 -1510.3912 261287.4 -1663.96733 Dec 2 283.5261 262606.4 -76.94697 Jan 3 2191.2067 263925.5 1296.31016 Feb 3 869.8936 265138.6 1357.50918 Mar 3 -3037.2958 266351.7 1462.58453 Apr 3 -5955.9043 267682.5 -2863.63748 May 3 -11875.1351 269013.4 -2294.23726 Jun 3 -13869.4003 270265.8 -1528.41786 Jul 3 8289.7110 271518.3 -2540.97496 Aug 3 11142.3402 272552.8 1655.89717 Sep 3 9370.4705 273587.3 3644.26827 Oct 3 4100.9753 274512.3 4428.72687 Nov 3 -1510.3912 275437.3 2760.05672 Dec 3 283.5261 276176.1 1455.41127 Jan 4 2191.2067 276914.8 -1977.99741 Feb 4 869.8936 277416.8 -1183.73207 Mar 4 -3037.2958 277918.9 155.40960 Apr 4 -5955.9043 278366.9 -2261.01814 May 4 -11875.1351 278815.0 200.17636 Jun 4 -13869.4003 279264.0 -401.61353 Jul 4 8289.7110 279713.1 -743.77992 Aug 4 11142.3402 280059.3 -15.61528 Sep 4 9370.4705 280405.5 2524.04833 Oct 4 4100.9753 280735.2 3349.83285 Nov 4 -1510.3912 281064.9 1922.48861 Dec 4 283.5261 281361.7 1010.74946 Jan 5 2191.2067 281658.5 -3659.75293 Feb 5 869.8936 281687.5 -2149.42403 Mar 5 -3037.2958 281716.5 -1843.21880 Apr 5 -5955.9043 281251.9 -79.97705 May 5 -11875.1351 280787.2 5439.88694 Jun 5 -13869.4003 279823.5 5356.92555 Jul 5 8289.7110 278859.7 2652.58767 Aug 5 11142.3402 277330.9 2252.76849 Sep 5 9370.4705 275802.1 7127.44828 Oct 5 4100.9753 273644.3 760.74440 Nov 5 -1510.3912 271486.5 -150.08824 Dec 5 283.5261 268848.9 -3271.45075 Jan 6 2191.2067 266211.4 631.42351 Feb 6 869.8936 263509.6 -203.52254 Mar 6 -3037.2958 260807.9 -2572.59227 Apr 6 -5955.9043 258386.1 922.78536 May 6 -11875.1351 255964.3 1967.78522 Jun 6 -13869.4003 254048.8 -4807.39107 Jul 6 8289.7110 252133.2 -1866.94386 Aug 6 11142.3402 250550.5 -699.79965 Sep 6 9370.4705 248967.7 -3675.15646 Oct 6 4100.9753 247637.4 -1095.34737 Nov 6 -1510.3912 246307.1 -1374.66702 Dec 6 283.5261 245248.0 1573.47267 Jan 7 2191.2067 244188.9 2160.84913 Feb 7 869.8936 243360.4 808.74637 Mar 7 -3037.2958 242531.8 -2414.48006 Apr 7 -5955.9043 242130.6 910.34924 May 7 -11875.1351 241729.3 -4300.19922 Jun 7 -13869.4003 242286.2 -1577.81636 Jul 7 8289.7110 242843.1 -3198.80998 Aug 7 11142.3402 244781.8 -7591.15198 Sep 7 9370.4705 246720.5 -9121.99500 Oct 7 4100.9753 249925.8 -8928.76282 Nov 7 -1510.3912 253131.1 -5357.65939 Dec 7 283.5261 257106.6 -1625.11378 Jan 8 2191.2067 261082.1 1045.66860 Feb 8 869.8936 264802.0 2675.06747 Mar 8 -3037.2958 268522.0 7561.34265 Apr 8 -5955.9043 272098.2 7820.66448 May 8 -11875.1351 275674.5 3630.60853 Jun 8 -13869.4003 279212.6 6649.79452 Jul 8 8289.7110 282750.7 1669.60400 Aug 8 11142.3402 286165.0 -1426.30336 Sep 8 9370.4705 289579.2 -5650.71174 Oct 8 4100.9753 292860.5 -8385.44386 > m$win s t l 941 19 13 > m$deg s t l 0 1 1 > m$jump s t l 95 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/rcomp/tmp/1gveu1260225950.ps",horizontal=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/281r51260225950.ps",horizontal=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/32l8u1260225950.ps",horizontal=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/431r91260225950.ps",horizontal=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/5i3cd1260225950.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/65hkd1260225951.tab") > system("convert tmp/1gveu1260225950.ps tmp/1gveu1260225950.png") > system("convert tmp/281r51260225950.ps tmp/281r51260225950.png") > system("convert tmp/32l8u1260225950.ps tmp/32l8u1260225950.png") > system("convert tmp/431r91260225950.ps tmp/431r91260225950.png") > > > proc.time() user system elapsed 1.076 0.628 1.725