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(1038.00,934.00,988.00,870.00,854.00,834.00,872.00,954.00,870.00,1238.00,1082.00,1053.00,934.00,787.00,1081.00,908.00,995.00,825.00,822.00,856.00,887.00,1094.00,990.00,936.00,1097.00,918.00,926.00,907.00,899.00,971.00,1087.00,1000.00,1071.00,1190.00,1116.00,1070.00,1314.00,1068.00,1185.00,1215.00,1145.00,1251.00,1363.00,1368.00,1535.00,1853.00,1866.00,2023.00,1373.00,1968.00,1424.00,1160.00,1243.00,1375.00,1539.00,1773.00,1906.00,2076.00,2004.00) > 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 18.582278 1007.5159 11.9018104 Feb 1 -9.750896 999.9102 -56.1593362 Mar 1 -36.083639 992.3046 31.7790857 Apr 1 -157.069829 985.2614 41.8084471 May 1 -154.056117 978.2182 29.8379069 Jun 1 -141.418453 972.1259 3.2925777 Jul 1 -67.380778 966.0335 -26.6527626 Aug 1 -24.565573 960.7366 17.8289342 Sep 1 28.249771 955.4397 -113.6895077 Oct 1 253.050264 954.2777 30.6720726 Nov 1 162.850831 953.1156 -33.9664199 Dec 1 127.592264 957.1433 -31.7355568 Jan 2 18.582278 961.1710 -45.7532752 Feb 2 -9.750896 960.7150 -163.9641476 Mar 2 -36.083639 960.2591 156.8245486 Apr 2 -157.069829 954.8304 110.2393816 May 2 -154.056117 949.4018 199.6543129 Jun 2 -141.418453 944.2945 22.1239551 Jul 2 -67.380778 939.1872 -49.8064139 Aug 2 -24.565573 935.0650 -54.4993965 Sep 2 28.249771 930.9427 -72.1925178 Oct 2 253.050264 929.6823 -88.7326101 Nov 2 162.850831 928.4219 -101.2727753 Dec 2 127.592264 939.0055 -130.5977204 Jan 3 18.582278 949.5890 128.8287530 Feb 3 -9.750896 967.8023 -40.0514094 Mar 3 -36.083639 986.0156 -23.9320032 Apr 3 -157.069829 999.7963 64.2734828 May 3 -154.056117 1013.5770 39.4790671 Jun 3 -141.418453 1024.1124 88.3060397 Jul 3 -67.380778 1034.6478 119.7330012 Aug 3 -24.565573 1047.1315 -22.5658789 Sep 3 28.249771 1059.6151 -16.8648977 Oct 3 253.050264 1078.4355 -141.4857779 Nov 3 162.850831 1097.2559 -144.1067310 Dec 3 127.592264 1123.0467 -180.6389479 Jan 4 18.582278 1148.8375 146.5802536 Feb 4 -9.750896 1186.2071 -108.4561832 Mar 4 -36.083639 1223.5767 -2.4930514 Apr 4 -157.069829 1274.5727 97.4971476 May 4 -154.056117 1325.5687 -26.5125550 Jun 4 -141.418453 1378.2788 14.1396389 Jul 4 -67.380778 1430.9890 -0.6081783 Aug 4 -24.565573 1470.8889 -78.3233731 Sep 4 28.249771 1510.7889 -4.0387066 Oct 4 253.050264 1530.2368 69.7129235 Nov 4 162.850831 1549.6847 153.4644808 Dec 4 127.592264 1559.8262 335.5815384 Jan 5 18.582278 1569.9677 -215.5499856 Feb 5 -9.750896 1586.3986 391.3522486 Mar 5 -36.083639 1602.8296 -142.7459487 Apr 5 -157.069829 1623.8121 -306.7422352 May 5 -154.056117 1644.7945 -247.7384232 Jun 5 -141.418453 1663.8969 -147.4784710 Jul 5 -67.380778 1682.9993 -76.6185299 Aug 5 -24.565573 1703.8633 93.7023190 Sep 5 28.249771 1724.7272 153.0230292 Oct 5 253.050264 1749.1335 73.8161992 Nov 5 162.850831 1773.5399 67.6092964 > m$win s t l 591 19 13 > m$deg s t l 0 1 1 > m$jump s t l 60 2 2 > m$inner [1] 2 > m$outer [1] 0 > postscript(file="/var/www/html/freestat/rcomp/tmp/1q1nx1292883368.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/2q1nx1292883368.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/3q1nx1292883368.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/41bm01292883368.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/5fkkr1292883368.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/6030x1292883368.tab") > > try(system("convert tmp/1q1nx1292883368.ps tmp/1q1nx1292883368.png",intern=TRUE)) character(0) > try(system("convert tmp/2q1nx1292883368.ps tmp/2q1nx1292883368.png",intern=TRUE)) character(0) > try(system("convert tmp/3q1nx1292883368.ps tmp/3q1nx1292883368.png",intern=TRUE)) character(0) > try(system("convert tmp/41bm01292883368.ps tmp/41bm01292883368.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.414 0.878 1.537