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)]<e^{rT}S(0)$ 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))
X
but I suspect the normal should be $N(0, T)$. What do you think? Then when you compute themean
, that term vanishes since its expected value is zero. $\endgroup$