# How to reduce variance in Monte Carlo using Control Variates when spot prices are decreasing?

I'm trying to use the Control Variates technique to reduce the variance of the estimate obtained from a Monte Carlo simulation for option pricing. As suggested in the book by Glasserman I'm using this control variate estimator

$$\text{"option price at time 0"} \approx \hat Y = \frac 1n\sum_{i=1}^n Z_i$$

where $$Z_i$$ are the components of the vector $$Z = Y-\theta(X-\mathbb E[X])$$, with $$V=e^{-rT}(S(T)-K)$$ vector of discounted payoffs (outputs of the Monte Carlo simulation), $$X=e^{-rT}S(T)$$ and $$S(T)$$ is the vector of spot prices at expiry time $$T$$ generated in the simulation, $$\theta$$ is a constant chosen to be the minimzer of $$Z$$ that is $$\theta=\dfrac{\text{cov}(Y,X)}{\text{var}(X)}$$. Finally, under the risk-neutral measure $$X$$ is a martingale and $$\mathbb E[X]=S(0)$$.

The last identity comes from the previous book "the absence of arbitrage is essentially equivalent to the requirement that appropriately discounted asset prices be martingales. Any martingale with a known initial value provides a potential control variate precisely because its expectation at any future time is its initial value".

What I don't get is the basic assumption $$\mathbb E[S(T)]=e^{rT}S(0)$$ which implies that the spot prices will keep growing in the future ($$e^{rT}$$ is strictly bigger than $$1$$).

In the example I'm working on - option under the Schwartz model $$dS = \alpha(\mu-\log S)Sdt + \sigma S dW$$ - the initial spot price is $$S(0)=22.93$$ but almost all (98.5%) the spot prices $$S(T)$$ computed with the Monte Carlo simulation are smaller than $$S(0)$$, hence $$\mathbb E[S(T)] and $$\hat Y$$ is a bad estimator of the option price (exact solution is 2.08 while the control variate estimator is 5.88).

So I guess that a different $$X$$ has to be chosen, any idea on possible candidates?

This is the output of the Matlab code used to compute the price V of the option at time 0 using Monte Carlo simulations with the suggestion by jherek

V_MC_standard = 0.070141, std = 0.000144
V_MC_controlv = 0.070216, std = 0.000074


and this is the code

S0 = 1; % spot price at time 0
K = 1; % strike prices
T = 1/2; % expiry time
r = .1; % risk-free interest rate
alpha = .2;
sigma = 0.4;
mu = 0.3;

%% Standard Monte Carlo
N = 1e6;
X = log(S0)*exp(-alpha*T) + (mu-sigma^2/2/alpha-(mu-r)/alpha)*(1-exp(-alpha*T)) + sigma*sqrt(1-exp(-2*alpha*T))/sqrt(2*alpha)*randn(N,1);
S = exp(X);
V = exp( -r*T ) * max(0,S-K);
V0 = mean(V);
fprintf('V_MC_standard = %f, std = %f\n' , V0 , std(V)/sqrt(N) );

%% Control Variates
VC = exp(-r*T)*S; % mean(VC) == S0
C = cov(V,VC); % the covariance matrix
theta = C(1,2)/C(2,2); % the optimal theta
F = exp( exp(-alpha*T)*log(S0) + (mu-sigma^2/2/alpha-(mu-r)/alpha)*(1-exp(-alpha*T)) + sigma^2/4/alpha*(1-exp(-2*alpha*T)) );
V = V-theta*(VC-exp(-r*T)*F);
V0 = mean(V); % Controlled Monte Carlo estimate of the option value
fprintf('V_MC_controlv = %f, std = %f\n' , V0 , std(V)/sqrt(N))

• Just a minor comment: I believe you are using this equation to compute X but I suspect the normal should be $N(0, T)$. What do you think? Then when you compute the mean, that term vanishes since its expected value is zero. – rvignolo Oct 2 '20 at 16:18
• @rvignolo The explicit formula for the $\log$ of the spot price is $$X_t = X_0 e^{-\alpha t} + \Big(\mu-\frac{\sigma^2}{2\alpha}-\lambda\Big)(1-e^{-\alpha t}) + \sigma e^{-\alpha t} \underbrace{\int_0^t e^{\alpha s} dW_s}_{I}$$ $I$ is the "sum" of independent Normals so it is also Normal with mean the sum of the means, ie $0$, and variance the sum of the variances, ie $\displaystyle\int_0^t e^{2\alpha s}ds=\frac{e^{2\alpha t}-1}{2\alpha}$ hence – sound wave Oct 2 '20 at 18:00
• hence $I \sim N\Bigg(0,\dfrac{e^{2\alpha t}-1}{2\alpha}\Bigg) \sim \sqrt{\dfrac{e^{2\alpha t}-1}{2\alpha}} N(0,1)$ and $$\sigma e^{-\alpha t} I \sim \sigma \sqrt{e^{-2\alpha t}\frac{e^{2\alpha t}-1}{2\alpha}} N(0,1) \sim \sigma\sqrt{\dfrac{1-e^{-2\alpha t}}{2\alpha}} N(0,1)$$ – sound wave Oct 2 '20 at 18:00
• Obvious first place to look is the solution of the PDE...did you derive it yourself or get one from a Numerical Recipes sort of text? – Chris Oct 4 '20 at 5:49
• @Chris The solution is derived by Schwartz himself at page 5 here sci-hub.st/10.1111/j.1540-6261.1997.tb02721.x – sound wave Oct 4 '20 at 6:00

Some of the assumptions here are wrong. The issue here is that $$S_0 \neq e^{-rT} E[S],$$ but $$F = E[S].$$
And thus Z should be Z=V-theta*(VC-exp(-rT)*F). If you output mean(VC) it's very clear.
• Thank you for help, I tried your formula for Z but then mean(Z) is exactly equal (up to 10 decimal digits) to V_exact, and by lowering the number of simulations N to 1e1 the precision increases (equal up to 14 decimal digits). It seems quite strange, is this how it is supposed to work? – sound wave Oct 5 '20 at 11:50
• The formula $S_0 = e^{-rT} E[S]$ is written in the book by Glasserman look here i.imgur.com/zlBFzTB.png – sound wave Oct 5 '20 at 11:51
• The formula $F=S_0 e^{rT}$ implies that spot prices increase in time, but usng the parameters estmated from market data, they decrease if $T$ is low, increase if $T$ is higher, I think this is due to the fact that the market data that I use to estimate the parameters don't always increase, they are the first 200 monthly oil spot prices downloadable from here eia.gov/dnav/pet/pet_pri_spt_s1_m.htm (sheet Data 1) – sound wave Oct 5 '20 at 11:59