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 = '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 216.659725 9517.959 -34.6185466 Feb 1 -524.326279 9520.554 84.7719805 Mar 1 -174.454657 9523.150 -264.6951176 Apr 1 148.134258 9522.729 72.1366729 May 1 -738.622861 9522.308 -196.6855021 Jun 1 217.654814 9519.207 -5.8619253 Jul 1 -166.900490 9516.106 213.7946308 Aug 1 4.899093 9513.709 479.3915984 Sep 1 -47.300641 9511.313 -27.0121160 Oct 1 519.336153 9505.144 13.5200682 Nov 1 389.139291 9498.975 29.8859085 Dec 1 155.781668 9473.355 -377.1368002 Jan 2 216.659725 9447.735 72.6048112 Feb 2 -524.326279 9410.349 148.9776698 Mar 2 -174.454657 9372.962 -65.5070969 Apr 2 148.134258 9346.342 -7.4758901 May 2 -738.622861 9319.722 118.9013513 Jun 2 217.654814 9303.515 105.8302088 Jul 2 -166.900490 9287.308 -173.4079543 Aug 2 4.899093 9275.856 2.2450035 Sep 2 -47.300641 9264.403 -388.1027207 Oct 2 519.336153 9259.546 168.1181644 Nov 2 389.139291 9254.688 -15.8272944 Dec 2 155.781668 9256.916 -94.6975841 Jan 3 216.659725 9259.144 129.1964463 Feb 3 -524.326279 9269.970 -105.6441395 Mar 3 -174.454657 9280.797 107.6576495 Apr 3 148.134258 9291.973 126.8923956 May 3 -738.622861 9303.150 -17.5268237 Jun 3 217.654814 9308.675 -341.3295245 Jul 3 -166.900490 9314.200 322.7007542 Aug 3 4.899093 9323.027 -204.9264246 Sep 3 -47.300641 9331.855 -6.5542853 Oct 3 519.336153 9352.527 298.1371071 Nov 3 389.139291 9373.199 -328.3378443 Dec 3 155.781668 9393.702 105.5167440 Jan 4 216.659725 9414.205 -201.8643476 Feb 4 -524.326279 9428.622 -165.2958062 Mar 4 -174.454657 9443.040 283.4151101 Apr 4 148.134258 9467.301 71.5651137 May 4 -738.622861 9491.562 266.0611519 Jun 4 217.654814 9520.373 -66.0278672 Jul 4 -166.900490 9549.184 -176.2839069 Aug 4 4.899093 9574.707 -510.6057282 Sep 4 -47.300641 9600.229 235.0717686 Oct 4 519.336153 9623.836 168.8282002 Nov 4 389.139291 9647.442 68.4182881 Dec 4 155.781668 9679.859 27.3591446 Jan 5 216.659725 9712.276 -272.9356787 Feb 5 -524.326279 9740.159 79.1673484 Mar 5 -174.454657 9768.042 352.4127502 Apr 5 148.134258 9781.412 -228.5463042 May 5 -738.622861 9794.782 -7.1593241 Jun 5 217.654814 9807.887 164.4579532 Jul 5 -166.900490 9820.992 51.9082100 Aug 5 4.899093 9832.330 -72.2292663 Sep 5 -47.300641 9843.668 96.6325755 Oct 5 519.336153 9853.119 -378.4548069 Nov 5 389.139291 9862.569 181.2914669 Dec 5 155.781668 9878.502 38.7158883 Jan 6 216.659725 9894.436 0.9046298 Feb 6 -524.326279 9919.817 -129.4906191 Mar 6 -174.454657 9945.198 49.2565068 Apr 6 148.134258 9974.150 -25.2840153 May 6 -738.622861 10003.101 -149.4785027 Jun 6 217.654814 10036.468 156.8772588 Jul 6 -166.900490 10069.834 -224.9340003 Aug 6 4.899093 10090.802 312.2986444 Sep 6 -47.300641 10111.770 88.5306071 Oct 6 519.336153 10129.561 -280.8967487 Nov 6 389.139291 10147.351 44.5095516 Dec 6 155.781668 10164.811 276.4076819 Jan 7 216.659725 10182.270 281.0701323 Feb 7 -524.326279 10198.831 63.4948913 Mar 7 -174.454657 10215.393 -484.9379749 > m$win s t l 751 19 13 > m$deg s t l 0 1 1 > m$jump s t l 76 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/freestat/rcomp/tmp/1v3221291927977.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/2v3221291927977.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/3v3221291927977.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/45c151291927977.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/51lze1291927977.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/654y21291927977.tab") > > try(system("convert tmp/1v3221291927977.ps tmp/1v3221291927977.png",intern=TRUE)) character(0) > try(system("convert tmp/2v3221291927977.ps tmp/2v3221291927977.png",intern=TRUE)) character(0) > try(system("convert tmp/3v3221291927977.ps tmp/3v3221291927977.png",intern=TRUE)) character(0) > try(system("convert tmp/45c151291927977.ps tmp/45c151291927977.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.477 0.887 1.598