  We consider the problem of Bayesian estimation of static parameters
associated to a partially and discretely observed diffusion process. We assume
that the exact transition dynamics of the diffusion process are unavailable,
even up-to an unbiased estimator and that one must time-discretize the
diffusion process. In such scenarios it has been shown how one can introduce
the multilevel Monte Carlo method to reduce the cost to compute posterior
expected values of the parameters for a pre-specified mean square error (MSE).
These afore-mentioned methods rely on upon the Euler-Maruyama discretization
scheme which is well-known in numerical analysis to have slow convergence
properties. We adapt stochastic Runge-Kutta (SRK) methods for Bayesian
parameter estimation of static parameters for diffusions. This can be
implemented in high-dimensions of the diffusion and seemingly under-appreciated
in the uncertainty quantification and statistics fields. For a class of
diffusions and SRK methods, we consider the estimation of the posterior
expectation of the parameters. We prove that to achieve a MSE of
$\mathcal{O}(\epsilon^2)$, for $\epsilon>0$ given, the associated work is
$\mathcal{O}(\epsilon^{-2})$. Whilst the latter is achievable for the Milstein
scheme, this method is often not applicable for diffusions in dimension larger
than two. We also illustrate our methodology in several numerical examples.
