Control Variates, Monte Carlo integration

Question

Exercise: Calculate $P(N>2.5)$ where $N$~$N(0,1)$ through simple monte carlo integration, and then use control variables to reduce the variance of my estimator.

I did

#Simulate Normal Random Variables
 randnormal<- function(n){
        Z<-NULL
        for(i in 1:n){
            U1=runif(1)
            U2=runif(1)
            Z[i]=sqrt(-2*log(U1))*cos(2*U2*pi)
        }
    return(Z)}

Then I realized the simple monte carlo integration through

#Monte Carlo Integration
montecarlo<-function(n){
    Z=randnormal(n)
    sum(Z>2.5)/n}

My question is how can I choose my control variable in this case, there is a criteria for selecting a control variable?

Answer 1

The idea behind "control variates" is to use the same Monte-Carlo simulation to estimate the known expectation of a function $g$ that is similar to the one, $f$, whose expectation you are trying to estimate. When $g$ and $f$ are strongly correlated, you can use the errors in estimating $g$ to correct the estimate of $f$.

Therefore there are two principal criteria governing your choice of control variate:

1. It must be a function whose expectation you know (or can somehow estimate to a high accuracy).

2. It must be correlated with the target function $f$.

In the present case, $f(z)$ is zero for all $z \lt 2.5$ and jumps to unity for all $z \gt 2.5$. Any function that will be correlated with it should either tend to increase with increasing $z$ or it should to tend to decrease: that will take care of (2).

Finding a $g$ whose expectation you can easily compute is more of a challenge. But in writing its formula as

$$\mathbb{E}(g) = \frac{1}{\sqrt{2\pi}}\int_\mathbb{R} g(z) \exp(-z^2/2) dz,$$

ask yourself whether there might be a relatively simple form of $g$ that enables this integral to be calculated. Indeed there is: one of the simplest choices would be $g(z) = z$, for then the integral could easily be computed using the substitution $z^2/2 = y, z dz = dy$.

Although $g(z)=z$ does work, following criterion (1) we might elect to modify it so that it most closely emulates $f$. I would suggest doing this by cutting $g$ off so that it's zero when $f$ is zero, whence

$$g(z) = z I(z \ge 2.5)$$

ought to do the trick. Its expectation is easily found in closed form as

$$\mathbb{E}(g) = \mathbb{E}(z I(z \ge 2.5)) = \frac{1}{\sqrt{2\pi}} \exp(-2.5^2/2).$$

You will find that this reduces the estimation variance by a factor of almost $200$--which means that with this technique, you need only $1/200$ as many iterations. That's a pretty good improvement!