¡Los resultados de Monte Carlo y PDE son diferentes para una Opción de Compra!

Question

Entonces, esta podría ser una pregunta bastante trivial, pero estoy teniendo problemas con el cálculo del valor de una opción de compra utilizando métodos de Monte Carlo y un método de Ecuaciones en Derivadas Parciales (PDE).

Cuando empecé, primero usé los parámetros:

Precio Spot = 0 a 20

Precio de Ejercicio = 10

Tasa de Interés = 0

Volatilidad = 0.25

Tiempo = 1

Y ambos métodos de Monte Carlo y PDE dieron resultados idénticos. Sin embargo, cuando cambié el rango de precio spot de 90 a 110 y el precio de ejercicio a 100, ¡ahora los métodos de Monte Carlo y PDE dan resultados diferentes! Por ejemplo, con un precio spot de 110, el método de Monte Carlo da (aproximadamente) 16.19 como precio de la opción. Sin embargo, el método de PDE en cambio da 10 como el precio. Y lo curioso es que esto sucede con cada código de Monte Carlo y PDE que pruebo (encontré algunos en internet y los he estado probando). ¿Alguien puede decirme por qué sucede esto?

Gracias de antemano.

Editar: Solo como un pequeño favor, ¿alguien puede decirme qué necesita ser modificado en este código para obtener el resultado correcto? Obtuve el código del enlace que mencioné a continuación

r=0.0; % Tasa de Interés
sigma=0.25; % Volatilidad del activo subyacente 
M=1600; % Número de puntos en el tiempo
N=160; % Número de puntos de precio de la acción
Smax=110; % Precio máximo de la acción considerado
Smin=90; % Precio mínimo de la acción considerado
T=1.; % Maduración (vencimiento) del contrato
E=100; % Precio de ejercicio del activo subyacente

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=(T/M); % Paso en el tiempo
ds=(Smax-Smin)/N; % Paso en el precio
% Inicializando la matriz del valor de la opción
v(1:N,1:M) = 0.0;
% Condiciones iniciales establecidas por la recompensa de la Opción de Compra Europea en el vencimiento: V(S,T)=máx(S-E,0);
v(1:N,1)=max((Smin+(0:N-1)*ds-E),zeros(size(1:N)))';
%v(1:N,1)=max((Smin+(0:N-1)*ds-E),(5))';

% Condiciones de frontera establecidas por la Opción de Compra Europea:
v(1,2:M)=zeros(M-1,1)'; % V(0,t)=0
v(N,2:M)=((N-1)*ds+Smin)-E*exp(-r*(1:M-1)*dt); % V(S,t)=S-Eexp[-r(T-t)] conforme S -> infinito.

% Determinando los coeficientes de la matriz del algoritmo explícito
aa=0.5*dt*(sigma*sigma*(1:N-2).*(1:N-2)-r*(1:N-2))';
bb=1-dt*(sigma*sigma*(1:N-2).*(1:N-2)+r)';
cc=0.5*dt*(sigma*sigma*(1:N-2).*(1:N-2)+r*(1:N-2))';

% Implementando el algoritmo explícito
for i=2:M,
v(2:N-1,i)=bb.*v(2:N-1,i-1)+cc.*v(3:N,i-1)+aa.*v(1:N-2,i-1);
end
% Reversión de los componentes de tiempo en la matriz ya que la solución de la ecuación de Black-Scholes fue realizada hacia atrás
v=fliplr(v);

Answer 1

Sin duda, recomendaría Volopta como una fuente confiable de códigos fuente de ingeniería financiera autocontenidos y comentados (útiles para prototipos/comprensión pero claramente no código de producción). Por ejemplo, he copiado y pegado el solucionador explícito de EDP que estás buscando (centrado en espacio, hacia atrás en el tiempo) a continuación (+ editado para mayor claridad + rendimiento mejorado, ver [Edits]).

Ejecutando el código en tu configuración: $S_0=110, K=100, T=1, r=q=0, \sigma=0.25$ da $16.1948$ (todavía puedes refinar las mallas de tiempo/punto, aunque ten en cuenta que los esquemas explícitos son solo condicionalmente estables), mientras que el precio de Black-Scholes es $16.1904$

% Opción de compra europea por Diferencias Finitas Explícitas
% Ejemplo tomado del libro de Clewlow y Strickland
% "Implementing Derivatives Models", Figura 3.8 (página 61)

clc; clear;

% Configuración
T = 1;                                  % Maduración
Spot = 110;                             % Precio del activo subyacente
K = 100;                                % Precio de ejercicio
v = 0.25;                               % Volatilidad
r = 0.0;                                % Tasa libre de riesgo
q = 0.0;                                % Rendimiento de dividendos
N = 250;                                % Número de pasos en el tiempo
M = 250;                                % Número de pasos en el precio de la acción

% Mallas de tiempo y punto
J = 0:1:N;                              % Índices para paso en el tiempo
dt = T/(length(J)-1);                   % Incremento de tiempo
I = M:-1:-M;                            % Índices para paso en el precio de la acción
dx = v*sqrt(3*dt);                      % Incremento para el precio de la acción
S = Spot*exp(I.*dx);                    % Inicializar malla de puntos (no uniforme)

% Discretización de la EDP
mu = r - q - v^2/2;                     % Deriva para el proceso de la acción
pu = dt*(v^2/2/dx^2 + mu/2/dx);         % Probabilidad hacia arriba
pm = 1 - dt*v^2/dx^2 - r*dt;            % Probabilidad media
pd = dt*(v^2/2/dx^2 - mu/2/dx);         % Probabilidad hacia abajo
V = zeros(length(I), length(J));        % Inicializar precio de la opción

% Condición terminal (ausencia de arbitraje)
V(:,end) = max(S - K, 0);

// Trabajar hacia atrás a través de la rejilla
for j=N:-1:1        

    % Nodos espaciales internos en el punto t_j = (j-1)*dt
    I = 2:2*M;
    V(I,j) = pu*V(I-1,j+1) + pm*V(I,j+1) + pd*V(I+1,j+1);

    % Condición de límite inferior en S_{min} = S(2*M+1)
    V(2*M+1,j) = V(2*M,j);

    % Condición de límite superior en S_{max} = S(1)
    V(1,j) = V(2,j) + (S(1) - S(2));        

end

% Precio de Diferencias Finitas
fdPrice = V(M+1,1);

% Precio analítico de BS
DF0T = exp(-r*T);
F0T = Spot*exp((r-q)*T);
bsPrice = bsPrice(v, DF0T, F0T, K, T, 'C');

fprintf('Diferencias Finitas: %0.4f\n', fdPrice)
fprintf('Analítico: %0.4f\n', bsPrice)

[Edits]

Disculpas por no haber leído el código que proporcioné en mi respuesta: solo me aseguré de que produjera el resultado correcto. De hecho, este código era bastante desordenado (mezclaba los índices $i$ y $j$ utilizados una vez para espacio/tiempo respectivamente y luego para tiempo/espacio), además tenía un conjunto adicional de líneas dentro del bucle principal hacia atrás

for i=1:2*M
    V(i,j) = max(K - S(j), V(i,j));
end

eso no tenía sentido. He agregado mi propia versión corregida arriba. Debería ejecutarse más rápido también ya que he eliminado bucles "for" inútiles.

En cuanto a tu propia implementación, hay tantos problemas que te aconsejaría que comiences desde el que te di y hagas preguntas. Si necesitas una revisión de código adecuada, publícalo en Code Review Stack Exchange. Ejemplos de problemas con tu implementación incluyen:

• Ya sea que insertes tu código de MATLAB dentro de una función adecuada, o te asegures de borrar todas las variables previamente utilizadas antes de ejecutarlo. Si no lo haces, tus resultados pueden ser incorrectos (especialmente si cambias el número de nodos de tiempo/espacio en tu caso, de grandes M y N, a pequeños).

• Seleccionas $S_{\text{min}}=90$ y $S_{\text{max}}=110$ para fijar un precio de opción en $E=100$... $S_{\text{min}}$ y $S_{\text{max}}$ deberían corresponder a los límites del dominio espacial, es decir, idealmente, $S_{\text{min}}\approx 0$ y $S_{\text{max}}\approx +\infty$. Al restringir excesivamente el dominio espacial, no puedes esperar aproximarte razonablemente a la solución real (especialmente si tus condiciones de límite son incorrectas).

• Tienes un problema con tus condiciones de límite, pero también con tus coeficientes de matriz en mi opinión.