This files allows to reproduce the Monte Carlo results as in Table 3 of Sbrana and Silvestrini (IJF, 2024).
The code generate two sets of processes (DGP1 and DGP2). You can choose which one to select.
If you have any question please email to : giacomo.sbrana@neoma-bs.fr
library("forecast")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary("Mcomp")
library("M4comp2018")
THETAstructural<-function(y,steps){
likelihood<-function(parameters){q=abs(parameters[1]);co=abs(parameters[2])
state<-p<-rep(0,length(y)+1)
state[1]=y[1]
p[1]=10000
k<-inn<-rep(0,length(y))
k[1]=p[1]/(p[1]+1);
sigmae=0
p[2]=p[1]-p[1]*k[1]+q;
state[2]=co+state[1]+k[1]*inn[1]
for(t in (2):length(y)){
k[t]=p[t]/(p[t]+1);
inn[t]=y[t]-state[t]
sigmae=sigmae+inn[t]^2/(p[t]+1)
p[t+1]=p[t]-p[t]*k[t]+q;
state[t+1]=co+state[t]+k[t]*inn[t]
}
sum(log(p+1))+length(y)*log(sigmae/length(y))}
results<-optim(c(.01,.001),likelihood)
q=abs(results[[1]][1]);co=abs(results[[1]][2])
state<-p<-rep(0,length(y)+1)
state[1]=y[1]
p[1]=10000
k<-inn<-rep(0,length(y))
k[1]=p[1]/(p[1]+1);
sigmae=0
p[2]=p[1]-p[1]*k[1]+q;
state[2]=co+state[1]+k[1]*inn[1]
for(t in (2):length(y)){
k[t]=p[t]/(p[t]+1);
inn[t]=y[t]-state[t]
sigmae=sigmae+inn[t]^2/(p[t]+1)
p[t+1]=p[t]-p[t]*k[t]+q;
state[t+1]=co+state[t]+k[t]*inn[t]
}
PointF<-rep(0,steps);PointF[1]=state[length(y)+1];
for( t in 2:steps){PointF[t]=co*(t-1)+state[length(y)+1]}
H<-sigmae/(length(y)-1);P<-tail(p,1)*H;Eta<-q*H;
Inter<-c();Inter[1]=P;
for(j in 2:steps){Inter[j]=Inter[j-1]+Eta};
Interv<-c();
for(j in 1:steps){Interv[j]=Inter[j]+H};
prob1=.85;prob2=.95;
lower1<-PointF-qnorm((1+prob1)/2)*sqrt(Interv);
#for(d in 1:steps){if(lower1[d]<0){lower1[d]=0}};
upper1<-PointF+qnorm((1+prob1)/2)*sqrt(Interv);
lower2<-PointF-qnorm((1+prob2)/2)*sqrt(Interv);
#for(d in 1:steps){if(lower2[d]<0){lower2[d]=0}};
upper2<-PointF+qnorm((1+prob2)/2)*sqrt(Interv);
list(mean=PointF,lower=cbind(lower1,lower2),upper=cbind(upper1,upper2))
}
AUTO<-function(y,steps){forecast::forecast(auto.arima(y),h=steps,level = c(.85,.95))}
THETA<-function(y,steps){forecast::forecast(thetaf(y,h=steps,level = c(.85,.95)))}
ETS<-function(y,steps){forecast::forecast(ets(y),h=steps)}
competitors<-list(THETAstructural,ETS,AUTO)
####################### data Monte Carlo ##################################################
steps<-6
freq<-frequ<-1
times=10000
set.seed(1)
DGP1=list()
for(u in 1:times){
n=60+steps
eta=runif(1,.1,.5)*rnorm(n)
mu=y=trend=c()
alpha=runif(1,.1,5);beta=runif(1,.01,.1)
ar=runif(1,.7,.99)
delta=runif(1,-.7,-.1)
trend[1]=alpha+beta
mu[1]=eta[1];y[1]=trend[1]+mu[1]
for(t in 2:n){mu[t]=ar*mu[t-1]+eta[t]+delta*eta[t-1]
trend[t]=alpha+beta*t
y[t]=trend[t]+mu[t]
}
DGP1[[u]]=y
}
DGP2=list()
for(u in 1:times){
n=60+steps
eta=runif(1,.1,.5)*rnorm(n)
mu=y=c()
ar=runif(1,.1,.7)
delta=runif(1,-.7,0)
co=runif(1,.01,.1);
trend=somma=c();trend[1]=co
mu[1]=eta[1];y[1]=co;somma[1]=mu[1]
for(t in 2:n){mu[t]=ar*mu[t-1]+eta[t]+delta*eta[t-1]
trend[t]=trend[t-1]+co
somma[t]=mu[t]+somma[t-1]
y[t]=somma[t]+trend[t]
}
DGP2[[u]]=y
}
######################## Select the dataset to use! ##############################
DATA=DGP1
goods=length(DATA)
#############################################################################################
rmsse<-function(y,yf,steps){OutSample=tail(y,length(yf));sqrt(mean((OutSample[1:steps]-yf[1:steps])^2)/mean(diff(head(y,(length(y)-length(yf))),frequ)^2))}
mase<-function(y,yf,steps){OutSample=tail(y,length(yf));(mean(abs(OutSample[1:steps]-yf[1:steps])))/mean(abs(diff(head(y,(length(y)-length(yf))),frequ)))}
mis<-function(yout,lower,upper,prob,yin){
a<-b<-d<-rep(0,length(yout));
for(t in 1:length(yout)){
if(yout[t]<lower[t]){a[t]=(2/(1-prob))*(lower[t]-yout[t])}
if(yout[t]>upper[t]){b[t]=(2/(1-prob))*(yout[t]-upper[t])}
d[t]=upper[t]-lower[t]}
mean(a+b+d)/mean(abs(diff(yin,freq)))
}
par(mfrow=c(1,3))
resMASE<-resRMSSE<-resSMAPE<-matrix(NA,steps,length(competitors))
resCOVER85<-resCOVER95<-resMIS85<-resMIS95<-matrix(NA,1,length(competitors))
results1MASE<-resultsMASE<-results1RMSSE<-resultsRMSSE<-results1SMAPE<-resultsSMAPE<-forec<-list()
resultsCOVER85<-resultsMIS85<-resultsCOVER95<-resultsMIS95<-matrix(NA,goods,length(competitors))
results1COVER85<-results1MIS85<-results1COVER95<-results1MIS95<-list()
for(h in 1:length(DATA)){
y=ts(DATA[[h]],frequency = freq);
for(j in 1:length(competitors)){
forec[[j]]<-competitors[[j]](head(y,(length(y)-steps)),steps)
for(s in 1:steps){
resRMSSE[s,j]<-rmsse(y,forec[[j]]$mean,s)
resMASE[s,j]<-mase(y,forec[[j]]$mean,s)
}
resMIS85[1,j]<-mis(tail(y,steps),forec[[j]]$lower[,1],forec[[j]]$upper[,1],.85,head(y,length(y)-steps));
resMIS95[1,j]<-mis(tail(y,steps),forec[[j]]$lower[,2],forec[[j]]$upper[,2],.95,head(y,length(y)-steps));
}
results1RMSSE[[1]]=resRMSSE
results1MASE[[1]]=resMASE
resultsRMSSE[[h]]<-Reduce("+", results1RMSSE)
resultsMASE[[h]]<-Reduce("+", results1MASE)
resultsMIS85[h,]<-resMIS85
resultsMIS95[h,]<-resMIS95
data_frame <- data.frame(col1 = 1:steps)
}
print(round(Reduce("+", resultsRMSSE)/h,3))## [,1] [,2] [,3]
## [1,] 0.716 0.728 0.729
## [2,] 0.840 0.856 0.857
## [3,] 0.914 0.935 0.933
## [4,] 0.972 0.999 0.994
## [5,] 1.020 1.054 1.046
## [6,] 1.062 1.103 1.090 print(rbind(round(colMeans(resultsMIS85,na.rm = TRUE),3),round(colMeans(resultsMIS95,na.rm = TRUE),3)))## [,1] [,2] [,3]
## [1,] 5.554 5.957 5.729
## [2,] 7.115 7.569 7.294