Estoy intentando utilizar la técnica de las Variables de Control para reducir la varianza de la estimación obtenida de una simulación de Monte Carlo para la valoración de opciones. Como se sugiere en el libro de Glasserman estoy utilizando este estimador de variantes de control
$$ \text{"option price at time 0"} \approx \hat Y = \frac 1n\sum_{i=1}^n Z_i $$
donde $Z_i$ son los componentes del vector $Z = Y-\theta(X-\mathbb E[X])$ con $V=e^{-rT}(S(T)-K)$ vector de resultados descontados (resultados de la simulación de Montecarlo), $X=e^{-rT}S(T)$ y $S(T)$ es el vector de precios al contado en el momento del vencimiento $T$ generado en la simulación, $\theta$ es una constante elegida para ser el minimizador de $Z$ es decir $\theta=\dfrac{\text{cov}(Y,X)}{\text{var}(X)}$ . Finalmente, bajo la medida de riesgo neutral $X$ es una martingala y $\mathbb E[X]=S(0)$ .
La última identidad proviene del libro anterior "la ausencia de arbitraje es esencialmente equivalente al requisito de que los precios de los activos debidamente descontados sean martingalas. Cualquier martingala con un valor inicial conocido proporciona una variante de control potencial precisamente porque su expectativa en cualquier momento futuro es su valor inicial".
Lo que no entiendo es la suposición básica $\mathbb E[S(T)]=e^{rT}S(0)$ lo que implica que los precios al contado seguirán creciendo en el futuro ( $e^{rT}$ es estrictamente mayor que $1$ ).
En el ejemplo en el que estoy trabajando - opción bajo el modelo de Schwartz $dS = \alpha(\mu-\log S)Sdt + \sigma S dW$ - el precio inicial al contado es $S(0)=22.93$ pero casi todos (98,5%) los precios al contado $S(T)$ calculados con la simulación de Monte Carlo son menores que $S(0)$ Por lo tanto $\mathbb E[S(T)]<e^{rT}S(0)$ y $\hat Y$ es un mal estimador del precio de la opción (la solución exacta es 2,08 mientras que el estimador de la variable de control es 5,88).
Así que supongo que una $X$ tiene que ser elegido, ¿alguna idea sobre posibles candidatos?
Este es el resultado del código Matlab utilizado para calcular la solución exacta y la simulación Monte Carlo ( V denotan el valor de la opción)
V_exact = 2.032381
V_MC = 2.034583, error = 0.002202, std = 0.001169 <-- Monte Carlo
V_MC_CV = 5.807870, error = 3.775489, std = 0.000000 <-- Control Variatey este es el código
S0 = 22.93; % spot price at time 0
r = .1; % risk-free interest rate
T = 1/2; % expiry time
K = 18; % strike price
alpha = 0.069217; %
sigma = 0.087598; % estimated from data
mu = 3.058244; %
%% Analytical Solution of the PDE for option pricing
F = exp( exp(-alpha*T)*log(S0) + (mu-sigma^2/2/alpha-sigma*(mu-r)/alpha)*(1-exp(-alpha*T)) + sigma^2/4/alpha*(1-exp(-2*alpha*T)) );
V_exact = exp( -r*T ) * (F-K); % true value of the option at time 0
%% Standard Monte Carlo
N = 1e6; % number of simulations
X = log(S0)*exp(-alpha*T) + (mu-sigma^2/2/alpha-sigma*(mu-r)/alpha)*(1-exp(-alpha*T)) + sigma*sqrt(1-exp(-2*alpha*T))/sqrt(2*alpha)*randn(N,1); % analytical solution of the Schwartz SDE
S = exp(X);
V = exp( -r*T ) * (S-K);
V_MC = mean(V); % estimate of the option value at time 0
%% Control Variates
VC = exp(-r*T)*S; % mean(VC) should be almost equal to S0
C = cov(V,VC); % the covariance matrix
theta = C(1,2)/C(2,2); % the optimal theta
Z = V-theta*(VC-S0);
V_MC_CV = mean(Z); % controlled estimate of the option value at time 0
fprintf( 'V_exact = %f\nV_MC = %f, error = %f, std = %f\n' , V_exact , V_MC , abs(V_exact-V_MC) , std(V)/sqrt(N) )
fprintf( 'V_MC_CV = %f, error = %f, std = %f\n' , V_MC_CV , abs(V_exact-V_MC_CV) , std(Z)/sqrt(N) )Este es el código utilizado para estimar los parámetros a partir de los datos del mercado
oil_monthly_spot_prices = [ 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11 18.65 17.75 18.30 18.68 19.44 20.07 21.34 20.31 19.53 19.86 18.85 17.27 17.13 16.80 16.20 17.86 17.42 16.53 15.50 15.52 14.54 13.77 14.14 16.38 18.02 17.94 19.48 21.07 20.12 20.05 19.78 18.58 19.59 20.10 19.86 21.10 22.86 22.11 20.39 18.43 18.20 16.70 18.45 27.31 33.51 36.04 32.33 27.28 25.23 20.48 19.90 20.83 21.23 20.19 21.40 21.69 21.89 23.23 22.46 19.50 18.79 19.01 18.92 20.23 20.98 22.38 21.78 21.34 21.88 21.69 20.34 19.41 19.03 20.09 20.32 20.25 19.95 19.09 17.89 18.01 17.50 18.15 16.61 14.51 15.03 14.78 14.68 16.42 17.89 19.06 19.65 18.38 17.45 17.72 18.07 17.16 18.04 18.57 18.54 19.90 19.74 18.45 17.33 18.02 18.23 17.43 17.99 19.03 18.85 19.09 21.33 23.50 21.17 20.42 21.30 21.90 23.97 24.88 23.71 25.23 25.13 22.18 20.97 19.70 20.82 19.26 19.66 19.95 19.80 21.33 20.19 18.33 16.72 16.06 15.12 15.35 14.91 13.72 14.17 13.47 15.03 14.46 13.00 11.35 12.51 12.01 14.68 17.31 17.72 17.92 20.10 21.28 23.80 22.69 25.00 26.10 27.26 29.37 29.84 25.72 28.79 31.82 29.70 31.26 33.88 33.11 34.42 28.44 29.59 29.61 27.24 27.49 28.63 27.60 26.42 27.37 26.20 22.17 19.64 19.39 19.71 20.72 24.53 26.18 27.04 25.52 26.97 28.39 ];
% now calibrate this data to the Schwartz model
S = oil_monthly_spot_prices(1:200);
X = log(S);
n = length(X);
M = [X(1:end-1)',ones(n-1,1)];
b = X(2:end)';
c = M\b;
r = M*c-b;
h = 1; % set h to 1 for raw data, and to dt for simulated data
EstAlpha = -log(c(1))/h;
EstSigma = std(r)*sqrt(2*EstAlpha)/sqrt(1-c(1)^2);
EstMu = c(2)/(1-c(1))+EstSigma^2/2/EstAlpha;
fprintf('alpha = %f\nsigma = %f\nmu = %f\n',EstAlpha,EstSigma,EstMu)
