  In this work, we extend the hybrid Chernoff tau-leap method to the multilevel
Monte Carlo (MLMC) setting. Inspired by the work of Anderson and Higham on the
tau-leap MLMC method with uniform time steps, we develop a novel algorithm that
is able to couple two hybrid Chernoff tau-leap paths at different levels. Using
dual-weighted residual expansion techniques, we also develop a new way to
estimate the variance of the difference of two consecutive levels and the bias.
This is crucial because the computational work required to stabilize the
coefficient of variation of the sample estimators of both quantities is often
unaffordable for the deepest levels of the MLMC hierarchy. Our method bounds
the global computational error to be below a prescribed tolerance, $TOL$,
within a given confidence level. This is achieved with nearly optimal
computational work. Indeed, the computational complexity of our method is of
order $\mathcal{O}\left(TOL^{-2}\right)$, the same as with an exact method, but
with a smaller constant. Our numerical examples show substantial gains with
respect to the previous single-level approach and the Stochastic Simulation
Algorithm.
