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

Question

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))

Answer 1

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.

It suggests that the choice of parameters for the Schwartz model are not consistent with the interest rate r, unless a non-zero convenience yield is expected.