# Numerically validating rates of convergence of approximations of expectations?

In applied mathematics it is standard practise to often validate theoretical approximations using numerical simulations. Since these simulations typically use numerical methods that convergence very fast , e.g. Gaussian quadrature, spectral methods, we can use them to check approximations featuring terms on the order of, say, $$O(n^{-1})$$ or $$O(n^{-2})$$, as $$n$$ gets large without any problems since the the discretization error of the numerical method itself decays much faster and thus is negligble.

Note we are talking about numerically validating theoretical approximations. Both the theoretical component and the numerical component feature error which can be viewed as approximation error. So to avoid confusion I have used the term 'discretization error' for the error of the numerical approximation, and and I will use 'approximation error' to refer to the error of the theoretical approximation.

Now, consider the fields of probability and statistics and suppose we want to validate approximations of expectations. One typically uses Monte Carlo simulations to compute approximations of expectations. But Monte Carlo integration converges slowly, that is, the discretization error of Monte Carlo integration converges at $$O(n^{-1/2})$$. So it seems there might be a problem when it comes to validating higher order terms in a theoretical approximation of an expectation.

Let's look at an example. Suppose we have a random variable $$X$$ that depends on the sample size $$n$$ such that its variance is $$\sigma^2 = C/n$$ for some $$C>0$$. Then we have the following Taylor approximation for the expectation of a function $$g$$ of $$X$$: \begin{align} E[g(X)] & \approx g(\mu) + \frac{1}{2}g''(\mu)\sigma^2 \\ &= g(\mu) + \frac{C}{2n}g''(\mu). \end{align} This is a simple case for the sake of example. In realty, we could be dealing with a very complicated $$X$$ which could be a function of many other random variables, and we could have derived an approximation of $$E[g(X)]$$ featuring several relatively complicated terms that each have different rates of convergence with respect to $$n$$.

In any case, I am wondering how we can numerically validate that we have a correct second order approximation of the expectation $$E[g(X)]$$. If we can validate the simple case above, we can validate more complicated real world cases.

Suppose we approximate $$E[g(X)]$$ using Monte Carlo integration to get $$E[g(X)] = E_\text{MC}(n) + \varepsilon(n),$$ where $$E_\text{MC}(n)$$ is our Monte Carlo approximation and $$\varepsilon(n) = O(n^{-1/2})$$ is the discretization error. We consider $$E_\text{MC}(n)$$ to be our (numerical) 'reference solution'.

To check our theoretical Taylor series approximation we compare it against our reference solution $$E_\text{MC}(n)$$. First, we should find that $$g(\mu)$$ is close to $$E_\text{MC}(n)$$, matching it to say one or two decimal places. Then, we would like to check that our second order term is correct, that is, $$g(\mu) + \frac{C}{2n}g''(\mu)$$ should be even closer to $$E_\text{MC}(n)$$, probably matching it perfectly to, say, four or five decimal places. However, it appears there is no point checking this term against $$E_\text{MC}(n)$$ as the Monte Carlo error that has been neglected $$\varepsilon(n) = O(n^{-1/2})$$ which dominates over the $$O(n^{-2})$$ term in our theoretical approximation.

So, in general, how can we numerically validate theoretical approximations of expectations that feature terms converging at rates such as $$O(n^{-1})$$ or $$O(n^{-2})$$ when the Monte Carlo error dominates over them? Is it simply not possible to numerically validate higher order terms in theoretical approximations in probability and statistics?