  In this article we consider computing expectations w.r.t.~probability laws
associated to a certain class of stochastic systems. In order to achieve such a
task, one must not only resort to numerical approximation of the expectation,
but also to a biased discretization of the associated probability. We are
concerned with the situation for which the discretization is required in
multiple dimensions, for instance in space and time. In such contexts, it is
known that the multi-index Monte Carlo (MIMC) method can improve upon
i.i.d.~sampling from the most accurate approximation of the probability law.
Indeed by a non-trivial modification of the multilevel Monte Carlo (MLMC)
method and it can reduce the work to obtain a given level of error, relative to
the afore mentioned i.i.d.~sampling and relative even to MLMC. In this article
we consider the case when such probability laws are too complex to sampled
independently. We develop a modification of the MIMC method which allows one to
use standard Markov chain Monte Carlo (MCMC) algorithms to replace independent
and coupled sampling, in certain contexts. We prove a variance theorem which
shows that using our MIMCMC method is preferable, in the sense above, to
i.i.d.~sampling from the most accurate approximation, under assumptions. The
method is numerically illustrated on a problem associated to a stochastic
partial differential equation (SPDE).
