I have been asked a few time to provide an example how to use MAPAx, from the MAPA package for R, so I prepared this blog post. I admit the documentation could be better, so I put together this example from a retailing case – the original setting MAPAx was developed for (see paper here). This provides a step by step use of MAPAx.
First we load the package and some data.
y <- cbind(c(1030.829,893.551,1084.09,1278.436,936.708,915.322,885.713,2364.399,774.88,977.506,831.616,813.656,1569.956,967.925,806.146,1063.117,1123.787,906.686,996.498,1088.464,977.414,1128.328,896.594,1007.172,1046.379,1514.648,1626.115,2959.558,838.506,949.377,1433.307,805.048,1218.907,872.43,1730.103,865.734,1845.713,919.291,1003.363,1102.969,847.38,1965.26,809.673,953.193,1066.089,991.352,1115.694,1003.333,1090.48,930.749,1006.184,1239.068,873.707,728.583,881.316,1302.468,997.442,3481.118,841.042,997.601,1830.194,909.693,2358.406,2573.673,777.08,773.781,945.424,968.646,1074.589,1046.22,1155.559,990.627,931.943,786.285,2297.025,628.166,889.238,937.631,1113.925,870.384,1018.375,799.458,1542.328,1879.587,750.307,1087.948,1247.803,1052.352,883.899,793.126,913.736,1082.142,968.823,2099.176,841.224,964.227,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)) x <- cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0), c(0,0,0,1,1,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,1,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,1,1,1,0,1,0,0,0,1,1,1,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,1,1,0,0)) y <- y[!is.na(y)] y <- ts(y,frequency=12)
This toy set includes sales of a heavily promoted item, with three different promotion. Let’s visualise what we have:
plot(y)
cmp <- brewer.pal(3,"Set1")
for (i in 1:3){
points(time(y)[which(x[,i]==1)],y[which(x[,i]==1)],col=cmp[i],pch=i,cex=1.5)
}
legend("topright",c("Promo1","Promo2","Promo3"),col=cmp,pch=1:3)
Now let’s fit a MAPA and MAPAx and produce forecasts. I am modelling the sales in logs to capture the multiplicative promotional effects. For other types of applications this is not needed.
library(MAPA) mapafit.x <- mapaest(log(y),type="es",display=1,outplot=1,xreg=x) mapafit <- mapaest(log(y),type="es",outplot=1) frc.x <- mapafor(log(y),mapafit.x,ifh=13,fh=13,xreg=x,conf.lvl=c(0.8,0.9,0.95),comb="w.mean") frc <- mapafor(log(y),mapafit,ifh=13,fh=13,conf.lvl=c(0.8,0.9,0.95),comb="w.mean")
Let’s plot the results
par(mfrow=c(1,2))
plot(1:96, y, type="l",xlim=c(1,109),xaxs="i",main="MAPA",xlab="Period",ylab="Sales")
for (i in 1:96){
lines(i:(12+i),exp(frc$infor[,i]),col="red")
}
cmp <- brewer.pal(9,"Reds")[4:2]
for (i in 1:3){
polygon(c(97:109,109:97),exp(c(frc$PI[i,],rev(frc$PI[7-i,]))),col=cmp[i],border=NA)
}
lines(97:109,exp(frc$outfor),col="red")
plot(1:96, y, type="l",xlim=c(1,109),xaxs="i",main="MAPAx",xlab="Period",ylab="Sales")
for (i in 1:96){
lines(i:(12+i),exp(frc.x$infor[,i]),col="red")
}
cmp <- brewer.pal(9,"Reds")[4:2]
for (i in 1:3){
polygon(c(97:109,109:97),exp(c(frc.x$PI[i,],rev(frc.x$PI[7-i,]))),col=cmp[i],border=NA)
}
lines(97:109,exp(frc.x$outfor),col="red")
As you can see, MAPA provides a flat forecast. This is “correct” in the sense that apart from promotions this time series contains like information that could be captured by exponential smoothing. On the other hand, MAPAx is provided with the promotional information and makes use of it.
Really appreciate the how-to, but can you fix the formatting? Also, I think you may have forgot the library(MAPA) line (should be obvious, but you never know if a new R user is trying the code). Thanks again for the example!
Thanks! It should be fixed now – added library as well!
Hi Nikos,
I have this issue when I select “MAPA” package. i know this is due to “smooth” package but raised it because I use only MAPA and smooth is required for “MAPA”. I have raised it as an issue in github for forecast package. I have cross posted it below. Appreciate all your help to the forecasting community.
I get following error when I use forecast function. The problem arises due to conflict with “smooth” package. I need to select smooth package because “MAPA” package requires it.
Error in UseMethod(“forecast”) :
no applicable method for ‘forecast’ applied to an object of class “c(‘ARIMA’, ‘Arima’)”
Please see: https://github.com/trnnick/mapa/issues/2
Hi Nikos,
When I run
mapafit_x = mapaest(y,type=”es”,display=1,outplot=1,xreg=as.matrix(x))
I get this error
Aggregation level: 80/365 (21.92%)Error in if (results[[besti]]$ICs[ic] <= results[[i]]$ICs[ic]) { : missing value where TRUE/FALSE needed
Do you know what may be the reason?
I suspect the time series is very short, but I cannot be sure.
Hello Nikos,
Thanks a lot for making MAPA available!
Just like in your example, I’m using 3 regressors and the function mapafit.x is working fine.
However, when I run frc.x I get the following:
“Error in mapacalc(inobs, mapafit, fh = ifh.c, comb, outplot = 0, hybrid, :
Number of xreg input variables does not match mapafit specification.”
Do you have an hypothesis for why this is happening? Many thanks in advance
Not sure. The easiest way to get the xreg part working is to imput the complete matrix in both the estimation and forecast functions. It should be able to find what sample to use on its own from the target variable. It does not use dates to do that, but it counts observations, so the X’s do not have to be time series.
Hi Nikos
MAPAx sounds like an interesting alternative that would be worth trying out.
My experience with using ARIMAX and ETS so far in predicting 7-day forecasts for SKUs is that both models are very sensitive to trends just before the forecasting horizon: e.g. if the last 1-2 days see an increase in sales, both models assume an increase in every single day of the horizon leading to bad forecasting accuracy especially in the latest 3-7 days of the horizon where they extrapolate that initial trend too much .
If you have such experience yourself, do you know how MAPAx deals with ‘last-minute’ trends that do not last?
Thanks,
Ioannis
By construction MAPA and MAPAx would be resistant to such changes, as short term movements would be filtered out in the temporally aggregate views of the time series. Therefore, even if on the original sampling frequency the ETS model might have poor performance, at an aggregate level it would be more robust and eventually the final MAPA forecast that would combine information from both would be resistant to these effects.