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(46.426,46.855,153.135,64,31,2.25,2.25,2.3,22.6,1.5,10.65,34,81.75,106.5,0.525,24.025,5.25,9,12.8,25.05,0.3,75.75,54.75,1.526) > par10 = '0.1' > par9 = '1' > par8 = 'dumresult' > par7 = 'dum' > par6 = '12' > par5 = 'ZZZ' > par4 = 'NA' > par3 = 'NA' > par2 = 'Croston' > par1 = 'Input box' > par10 <- '0.1' > par9 <- '1' > 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 25 30.99292 -21.37706 83.36290 -49.10005 111.0859 26 30.99292 -21.63826 83.62410 -49.49952 111.4854 27 30.99292 -21.89817 83.88401 -49.89702 111.8829 28 30.99292 -22.15681 84.14265 -50.29257 112.2784 29 30.99292 -22.41419 84.40003 -50.68621 112.6721 30 30.99292 -22.67034 84.65618 -51.07796 113.0638 31 30.99292 -22.92528 84.91112 -51.46784 113.4537 32 30.99292 -23.17901 85.16485 -51.85590 113.8417 33 30.99292 -23.43156 85.41740 -52.24214 114.2280 34 30.99292 -23.68294 85.66879 -52.62660 114.6124 $period Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 25 1 1 1 1 1 26 1 1 1 1 1 27 1 1 1 1 1 28 1 1 1 1 1 29 1 1 1 1 1 30 1 1 1 1 1 31 1 1 1 1 1 32 1 1 1 1 1 33 1 1 1 1 1 34 1 1 1 1 1 In-sample error measures: ME RMSE MAE MPE MAPE MASE -6.710034 40.525651 33.154766 -1329.170336 1356.553004 1.021489 Forecasts: Point Forecast 25 30.99292 26 30.99292 27 30.99292 28 30.99292 29 30.99292 30 30.99292 31 30.99292 32 30.99292 33 30.99292 34 30.99292 > postscript(file="/var/www/html/rcomp/tmp/1uggy1271420284.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/2n8yj1271420284.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/31hd91271420284.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/4u9dc1271420284.tab") > try(system("convert tmp/1uggy1271420284.ps tmp/1uggy1271420284.png",intern=TRUE)) character(0) > try(system("convert tmp/2n8yj1271420284.ps tmp/2n8yj1271420284.png",intern=TRUE)) character(0) > > #-SERVER-wessa.org > > > > proc.time() user system elapsed 1.647 0.335 1.805