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(62,30,31,50,33,12,20,30,21.5,23,13.5,0.5,12,10,70.5,30,20.5,12,20,45,11.505,0,10,5.5,27.5,0.5,7,0,2.5,0,0,6.025,1,0,0,0,0,2,0,6,20,0,0,0,7,35,0,0,0,1) > par10 = '0.1' > par9 = '3' > par8 = 'dumresult' > par7 = 'dum' > par6 = '12' > par5 = 'ZZZ' > par4 = 'NA' > par3 = 'NA' > par2 = 'Croston' > par1 = 'Input box' > par10 <- '0.1' > par9 <- '3' > par8 <- 'dumresult' > par7 <- 'dum' > par6 <- '12' > par5 <- 'ZZZ' > par4 <- 'NA' > par3 <- 'NA' > par2 <- 'Croston' > par1 <- 'Input box' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Dr. Ian E. Holliday > #To cite this work: Ian E. Holliday, 2009, 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: > #Technical description: > if(par3!='NA') par3 <- as.numeric(par3) else par3 <- NA > if(par4!='NA') par4 <- as.numeric(par4) else par4 <- NA > par6 <- as.numeric(par6) #Seasonal Period > par9 <- as.numeric(par9) #Forecast Horizon > par10 <- as.numeric(par10) #Alpha > library(forecast) Loading required package: tseries Loading required package: quadprog Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from package:base : as.Date.numeric This is forecast 2.03 > if (par1 == 'CSV') { + xarr <- read.csv(file=paste('tmp/',par7,'.csv',sep=''),header=T) + numseries <- length(xarr[1,])-1 + n <- length(xarr[,1]) + nmh <- n - par9 + nmhp1 <- nmh + 1 + rarr <- array(NA,dim=c(n,numseries)) + farr <- array(NA,dim=c(n,numseries)) + parr <- array(NA,dim=c(numseries,8)) + colnames(parr) = list('ME','RMSE','MAE','MPE','MAPE','MASE','ACF1','TheilU') + for(i in 1:numseries) { + sindex <- i+1 + x <- xarr[,sindex] + if(par2=='Croston') { + if (i==1) m <- croston(x,alpha=par10) + if (i==1) mydemand <- m$model$demand[] + fit <- croston(x[1:nmh],h=par9,alpha=par10) + } + if(par2=='ARIMA') { + m <- auto.arima(ts(x,freq=par6),d=par3,D=par4) + mydemand <- forecast(m) + fit <- auto.arima(ts(x[1:nmh],freq=par6),d=par3,D=par4) + } + if(par2=='ETS') { + m <- ets(ts(x,freq=par6),model=par5) + mydemand <- forecast(m) + fit <- ets(ts(x[1:nmh],freq=par6),model=par5) + } + try(rarr[,i] <- mydemand$resid,silent=T) + try(farr[,i] <- mydemand$mean,silent=T) + if (par2!='Croston') parr[i,] <- accuracy(forecast(fit,par9),x[nmhp1:n]) + if (par2=='Croston') parr[i,] <- accuracy(fit,x[nmhp1:n]) + } + write.csv(farr,file=paste('tmp/',par8,'_f.csv',sep='')) + write.csv(rarr,file=paste('tmp/',par8,'_r.csv',sep='')) + write.csv(parr,file=paste('tmp/',par8,'_p.csv',sep='')) + } > if (par1 == 'Input box') { + numseries <- 1 + n <- length(x) + if(par2=='Croston') { + m <- croston(x) + mydemand <- m$model$demand[] + } + if(par2=='ARIMA') { + m <- auto.arima(ts(x,freq=par6),d=par3,D=par4) + mydemand <- forecast(m) + } + if(par2=='ETS') { + m <- ets(ts(x,freq=par6),model=par5) + mydemand <- forecast(m) + } + summary(m) + } Forecast method: Croston's method Model Information: $demand Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 36 14.36090 -6.657339 35.37915 -17.78372 46.50553 37 14.36090 -6.762168 35.48398 -17.94405 46.66586 38 14.36090 -6.866481 35.58829 -18.10358 46.82539 39 14.36090 -6.970283 35.69209 -18.26233 46.98414 40 14.36090 -7.073582 35.79539 -18.42031 47.14212 41 14.36090 -7.176386 35.89819 -18.57754 47.29935 42 14.36090 -7.278701 36.00051 -18.73402 47.45582 43 14.36090 -7.380535 36.10234 -18.88976 47.61157 44 14.36090 -7.481895 36.20370 -19.04477 47.76658 45 14.36090 -7.582786 36.30459 -19.19907 47.92088 $period Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 36 2.008931 0.8309275 3.186934 0.2073303 3.810531 37 2.008931 0.8250522 3.192809 0.1983447 3.819517 38 2.008931 0.8192058 3.198656 0.1894035 3.828458 39 2.008931 0.8133881 3.204473 0.1805060 3.837356 40 2.008931 0.8075985 3.210263 0.1716516 3.846210 41 2.008931 0.8018366 3.216025 0.1628396 3.855022 42 2.008931 0.7961022 3.221759 0.1540695 3.863792 43 2.008931 0.7903947 3.227467 0.1453407 3.872521 44 2.008931 0.7847139 3.233148 0.1366526 3.881209 45 2.008931 0.7790593 3.238802 0.1280046 3.889857 In-sample error measures: ME RMSE MAE MPE MAPE MASE -12.612959 19.192634 16.750877 -Inf Inf 1.465571 Forecasts: Point Forecast 51 7.148531 52 7.148531 53 7.148531 54 7.148531 55 7.148531 56 7.148531 57 7.148531 58 7.148531 59 7.148531 60 7.148531 > postscript(file="/var/www/html/rcomp/tmp/193hl1275297248.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,1)) > if (par2=='Croston') plot(m) > if ((par2=='ARIMA') | par2=='ETS') plot(forecast(m)) > plot(mydemand$resid,type='l',main='Residuals', ylab='residual value', xlab='time') > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2kuzo1275297248.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,2)) > acf(mydemand$resid, lag.max=n/3, main='Residual ACF', ylab='autocorrelation', xlab='time lag') > pacf(mydemand$resid,lag.max=n/3, main='Residual PACF', ylab='partial autocorrelation', xlab='time lag') > cpgram(mydemand$resid, main='Cumulative Periodogram of Residuals') > qqnorm(mydemand$resid); qqline(mydemand$resid, col=2) > 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,'Demand Forecast',6,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Point',header=TRUE) > a<-table.element(a,'Forecast',header=TRUE) > a<-table.element(a,'95% LB',header=TRUE) > a<-table.element(a,'80% LB',header=TRUE) > a<-table.element(a,'80% UB',header=TRUE) > a<-table.element(a,'95% UB',header=TRUE) > a<-table.row.end(a) > for (i in 1:length(mydemand$mean)) { + a<-table.row.start(a) + a<-table.element(a,i+n,header=TRUE) + a<-table.element(a,as.numeric(mydemand$mean[i])) + a<-table.element(a,as.numeric(mydemand$lower[i,2])) + a<-table.element(a,as.numeric(mydemand$lower[i,1])) + a<-table.element(a,as.numeric(mydemand$upper[i,1])) + a<-table.element(a,as.numeric(mydemand$upper[i,2])) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/3g4ew1275297248.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Actuals and Interpolation',3,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Time',header=TRUE) > a<-table.element(a,'Actual',header=TRUE) > a<-table.element(a,'Forecast',header=TRUE) > a<-table.row.end(a) > for (i in 1:n) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,x[i]) + a<-table.element(a,x[i] - as.numeric(m$resid[i])) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/48vvh1275297248.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'What is next?',1,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink(paste('http://www.wessa.net/Patrick.Wessa/rwasp_demand_forecasting_simulate.wasp',sep=''),'Simulate Time Series','',target='')) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink(paste('http://www.wessa.net/Patrick.Wessa/rwasp_demand_forecasting_croston.wasp',sep=''),'Generate Forecasts','',target='')) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink(paste('http://www.wessa.net/Patrick.Wessa/rwasp_demand_forecasting_analysis.wasp',sep=''),'Forecast Analysis','',target='')) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/51mdk1275297248.tab") > try(system("convert tmp/193hl1275297248.ps tmp/193hl1275297248.png",intern=TRUE)) character(0) > try(system("convert tmp/2kuzo1275297248.ps tmp/2kuzo1275297248.png",intern=TRUE)) character(0) > > #-SERVER-wessa.org > > > > proc.time() user system elapsed 2.563 0.338 2.771