The R code

This files allows to estimate and predict using the Markov Walk model. The prediction part provides both point forecast and prediction intervals (as required to participate in the M5 competition). In the first chuck you find two functions (i.e. “Bernoulli”, “MW”). The “Bernoulli” function takes as input the time series (y) and the steps you wish to predict. This function returns, as output, a list containing:

  1. bern i.e. the “Bernoulli representation” of the time series, that is, taking 1 for the non-zero values and 0 for the zero values.
  2. p i.e. the proportion of the non zero values
  3. lambda i.e. the conditional probability of \(p(\chi_t=1|\chi_{t-1}=1)=\lambda\)
  4. MC i.e. \(E(\chi_{n+h}|\chi_{n+h-1}\cdots\chi_{n})=\sum_{i=1}^{h}\delta^{i-1}\xi+\delta^h\chi_{n}\)
  5. delta i.e. \(p(\chi_t=1|\chi_{t-1}=1)-p(\chi_t=1|\chi_{t-1}=0)\)
  6. xi i.e. \(p(\chi_t=1|\chi_{t-1}=0)\)

The MW is the function providing the Markov Walk prediction. It takes as input the data and the steps and it returns as output a list with the point forecast (the expected mean) and the prediction intervals for the percentages (50%,67%,95%,99%) as required in the M5 competition.

Bernoulli<-function(y,steps){alpha=-1;co=mean(y[y>0])
bern=y;bern[bern>0]=1;
su=bern[-1]+bern[-length(y)]
di=bern[-1]-bern[-length(y)]
n00=length(su[su==0]);n11=length(su[su==2]);n01=length(di[di==-1]);n10=length(di[di==1])
p00=n00/(n00+n10);xi=n10/(n10+n00);p=mean(bern);lambda=n11/(n11+n01);
if(n00==0&n10==0){lambda=1};if(n11==0&n01==0){lambda=0};if(lambda<=0){lambda=.0001}
delta=p00+lambda-1;if(lambda==1){p=1;xi=0;delta=1;}

MC<-c();MC[1]=xi+delta*tail(bern,1);for(s in 2:steps){MC[s]=xi+delta*MC[s-1]}

m<-v<-c();m[1]=0;k=(lambda^2-1+sqrt(1-lambda^2))/lambda
for (t in 1:length(y)){
  v[t]<-y[t]-bern[t]*co-m[t]
  m[t+1] <-lambda*m[t]+k*v[t]
}
fo<-c();for(s in 1:steps){fo[s]<-co*MC[s]+lambda^(s-1)*m[length(m)]}

list(bern,p,lambda,MC,delta,xi,fo,v,k)
}

MW=function(y,steps){prob1=.5;prob2=.67;prob3=.95;prob4=.99;
mysum<-function(x,steps){d<-0;for(f in 0:(steps-2)){d<-d+x^(f*2)};d}
# This applies the classical multiplicative seasonal adjustment.
# This is done if the number of non-zero observations are more than 5% of the sample
if((length(y[y==0])/length(y))<.95){
  s=7;h=steps;  cma<-matrix(NA,length(y),1);
  for(g in 1:(length(y)-s+1)){cma[g+((s+1)/2)-1]<-mean(y[g:(g+s-1)])};residuals<-y/cma
  sfactors<-c();for(seas in 1:s){
    sfactors[seas]<-mean(na.omit(residuals[seq(seas,length(y)-s+seas,by=s)]))}
    sfactout<-rep(sfactors,length(y)+h)[(length(y)+1):(length(y)+h)]
  y<-y/rep(sfactors,ceiling(length(y)/s))[1:length(y)]
  y[is.na(y)]=0;y[y==Inf]=0  
}else{sfactout=rep(1,steps)}  

h<-c();le=14
for(i in 0:1000){Y=tail(y,(length(y)-i*le));if(length(Y)<100){break};ins=head(Y,length(Y)-steps);
if((length(ins[ins==0])/length(ins))<.99&length(ins)>100){h[i+1]=mean((tail(Y,steps)-Bernoulli(ins,steps)[[7]])^2)}  
}
if(length(h)!=0){y=tail(y,(length(y)-which.min(h[!is.na(h)])*le)+steps)}


co=mean(y[y>0]);ma=Bernoulli(y,steps);bern=ma[[1]];lambda=ma[[3]];MC=ma[[4]]
p=ma[[2]];delta=ma[[5]];fo=ma[[7]];v=ma[[8]];k=ma[[9]]

fo=fo*sfactout
if(p<1){vari=((-1 + lambda)*(1 +lambda - 2*p)*p)/(-1 +p)}else{vari=0}

Interv<-c();Interv[1]=vari*co^2+var(v);for(j in 2:steps){Interv[j]=vari*mysum(delta,j)*co^2+(var(v)*(1+k^2*(mysum(lambda,j))))};
lower0=fo
lower50<-fo-qnorm((1+prob1)/2)*sqrt(Interv);
lower67<-fo-qnorm((1+prob2)/2)*sqrt(Interv);
lower95<-fo-qnorm((1+prob3)/2)*sqrt(Interv);
lower99<-fo-qnorm((1+prob4)/2)*sqrt(Interv);
upper0=fo
upper50<-fo+qnorm((1+prob1)/2)*sqrt(Interv);
upper67<-fo+qnorm((1+prob2)/2)*sqrt(Interv);
upper95<-fo+qnorm((1+prob3)/2)*sqrt(Interv);
upper99<-fo+qnorm((1+prob4)/2)*sqrt(Interv);
list(mean=fo,lower=cbind(lower0,lower50,lower67,lower95,lower99),upper=cbind(upper0,upper50,upper67,upper95,upper99))
}

A practical example

This example is taken from the Markov Walk paper. Suppose we observe the following series of observations over 30 consecutive days:

y=c(4,2,4,5,6,7,4,7,0,0,0,3,4,2,5,6,7,9,7,9,3,2,3,4,7,0,0,0,6,5)
y
##  [1] 4 2 4 5 6 7 4 7 0 0 0 3 4 2 5 6 7 9 7 9 3 2 3 4 7 0 0 0 6 5

This is the Bernoulli representation of the series:

Bernoulli(y,28)[[1]]
##  [1] 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 1

Note that, since \(n_{00}=4\), \(n_{11}=21\), \(n_{01}=2\) and \(n_{10}=2\), then \(\hat{\phi}=.8\) and \(\hat{\lambda}=0.913\). Moreover, \(\omega=E(y_t^{>0})=5.04\)

lambda=Bernoulli(y,28)[[3]];k=Bernoulli(y,28)[[9]]
bern=Bernoulli(y,28)[[1]];co=mean(y[y>0])
m<-v<-c();m[1]=0
k=(lambda^2-1+sqrt(1-lambda^2))/lambda
for (t in 1:length(y)){
  v[t]<-y[t]-bern[t]*co-m[t]
  m[t+1] <-lambda*m[t]+k*v[t]
}
plot(ts(y),ylab="items sold",xlab ="Days",col="grey")
lines(m[1:length(y)]+co*bern,col="black")
lines(rep(co,30),lty="dotted")

If you have any question please email to : \(\textbf{giacomo.sbrana@neoma-bs.fr}\)