Antithetic and Control Variables

Antithetic Variables

Suppose we want to estimate a parameter q=E[X] and suppose we have generated X₁ and X₂, two rv's with mean θ. Then

Obviously we want to do this so that we minimize our simulation error, that is the variance of our estimator. If we generate X₁ and X₂ indpendently we have Cov(X₁,X₂)=0, but we could do even better if we generated them so that Cov(X₁,X₂)<0

Example (Estimating an Integral) Let's say we want to use simulation to estimate the integral
.
(not that we need simulation to do this, but it illustrates the idea). Now

where U~U[0,1]

Now the straightforward approach would be to generate U₁,U₂ iid U[0,1] and to estimate the integral by (exp(U₁)+exp(U₂))/2. This would have variance

How can we generate X₁ and X₂ so that E[X₁]=E[X₂]=θ and Cov(X₁,X₂)<0? One idea is to use X₁=exp(U) and X₂=exp(1-U). 1-U~U[0,1], so E[X₁]=E[X₂]=θ. Now

and so

This is illustrated in vr(1)

Control Variables

Again we want to estimate a parameter q=E[X], where X is the result of a simulation. Now suppose that for some other output variable, say Y, we have E[Y]=μ, and μ is known. Then for any constant c

E[X+c·(y-μ)

so X+c·(y-μ)] is an unbiased estimate of θ. Of course the optimal c is the one that minimizes the variance, and so

and with this value we have

The quantity Y is called a control variable for the simulation variable X. Let's see how this works. Let's assume X and Y are positively correlated, that is large values of X go with large values of Y. So if we see a large value of Y (relative to μ) then it is probably true that X is also large (relative to θ) and we should lower our estimate a bit. But this is exactly what happens because now c^* is negative.

Example let's do again the integral example above. We have X=exp(U). An obvious choice for Y is U itself. We know E[Y]=½, V[Y]=1/12 and

This is illustrated in vr(2)

Example Let's consider the following problem: as part of an "online" computer program we need to find the following integral, for different values of (t₁,..,t_k):

It is necessary to find this integral with a precision of ±0.01, that is a 95% CI for g_k(t) should have a length of no more than 0.002. Each time this integral needs to be found the computer program has to wait and so solving it as fast as possible is important.

Let's first do this using a simple simulation. Generate U₁,..,U_k~U[0,1] and calculate X₁=sin(∏[U_i+t_i]) Repeat this n times and then use as an estimate of g_k(t). Also calculate the sample standard deviation s. By the CLT should be approximately normal with mean g_k(t) and standard deviation σ, so a 95% CI for g_k(t) is given by ±1.96σ/√n, so the error is 2·1.96σ/√n=3.92σ/√n=0.002, so σ/√n=0.0005. In other words we need n large enough to ensure σ/√n≤0.0005. But we don't know σ! We can (approximately) do this as follows:

1) Do a "small" trial run, say 1000 simulations. Based on these you can calculate an estimate of σ=s
2) If 1000 was already enough (s≤0.0005) you are done, otherwise
σ/√n=0.0005 means n=(σ/0.0005)² and now run the simulation again. Note: you only need n-1000 runs, you already have 1000.

This is implemented in vr6.

As a numerical example, say t=(0.3,0.3,0.3). Then it turns out that n=110,000

Now, can we speed this up? let's try the control variable approach. the obvious choice here is Y=∏U_i, but in order to use it we need to know μ=E[Y]:

so our new estimator is as -cov(X,Y)/Var[Y]*(Y-2^-k)
This is implemented in vr6(which=2), and for the example above now n=10000 is enough!