¿Cómo realizar simulaciones de Monte Carlo para fijar el precio de un contrato a plazo según el modelo de Schwartz de reversión a la media?

Question

Objetivo: (1) Implementar el Método Explícito de Euler para resolver la EDP de los precios de las opciones bajo el modelo de Schwartz de reversión a la media. (2) Comparar con una simulación Monte Carlo.

Me quedo con el punto 1 (ver esto pregunta ) así que mientras tanto estoy trabajando en el punto 2. Básicamente, tengo 200 precios spot mensuales de una materia prima (petróleo) siguiendo la SDE de Schwartz de reversión a la media $dS = \alpha(\mu-\log S)Sdt + \sigma S dW$ . Esta materia prima es el activo subyacente del contrato a plazo con precio a la entrega $F(S,T)$ donde $S$ es el precio al contado en el momento 0 y $T$ es el tiempo de expiración. Un contrato de este tipo equivale a tener una opción con $\text{payoff}(S_T)=S_T-K$ donde $K=F(S,T)$ es el precio de ejercicio.

Dejar $\tau=T-t$ sea el tiempo hasta el vencimiento, la fórmula para el precio del contrato a plazo es

$$\tag{1}F(S,\tau)=\exp\Big(e^{-\alpha\tau}\log S +(\mu-\frac{\sigma^2}{2\alpha}-\frac{\mu-r}{\alpha})(1-e^{-\alpha\tau})+\frac{\sigma^2}{4\alpha}(1-e^{-2\alpha\tau})\Big)$$

y el valor de la opción equivalente es $V(S,t)=e^{-r(T-t)}F(S,t)$ . ( $\alpha, \sigma, \mu$ se estiman a partir de los datos)

Conociendo esta información, ¿cómo ejecutar (estoy usando matlab pero otros lenguajes también valen) una simulación Monte Carlo para fijar el precio de la opción?

He leído que estos son los pasos

• generar un gran número de trayectorias de precios aleatorias para el activo subyacente
• para cada camino, calcule el beneficio de la opción
• calcula la media de todos los pagos
• descontar la media a hoy

pero como estoy trabajando con datos reales de precios al contado del petróleo, estoy un poco confuso: ¿cómo podemos obtener el valor de la opción teniendo una materia prima subyacente con precios al contado conocidos generando trayectorias aleatorias de precios al contado?

Además, ¿cómo realizar los siguientes pasos utilizando las fórmulas anteriores?

---

EDIT: Primer intento

He escrito el código que calcula lo que creo (¿puedes confirmarlo?) que es la solución exacta para la valoración de opciones ( V_exact ) y, a continuación, calcula la solución aproximada mediante simulación Monte Carlo ( V_Monte_Carlo ). La simulación MC utiliza el primer precio al contado (datos reales) y los parámetros estimados para calcular "aleatoriamente" los siguientes precios al contado. Para calcular los precios al contado, en lugar de utilizar la ecuación de $dS$ que escribí más arriba es más fácil usar este que se obtiene de $dS$ dejando $X=\log S$ .

$$\tag{2}X_t = X_0 e^{-\alpha t} + \mu^*(1-e^{-\alpha t}) + \underbrace{\sigma e^{-\alpha t} \int_0^t e^{\alpha s} dW_s}_Z$$

donde $\mu^* = \mu-\dfrac{\sigma^2}{2\alpha}$ y $Z \sim \sigma\sqrt{\dfrac{1-e^{-2\alpha t}}{2\alpha}} \cdot N(0,1)$

Estos son los pasos de la simulación Monte Carlo:

1. El primer precio F del contrato a plazo - y por tanto de la opción V también - se fija en 0, ya que no hay ningún coste por suscribir un contrato a plazo. Calcule el primer valor de $X$ tomando el logaritmo del primer precio al contado
2. calcula el siguiente valor de $X$ mediante la fórmula (2)
3. calcular el siguiente valor del precio al contado tomando el exponencial de $X$
4. calcular el siguiente precio F del contrato a plazo mediante la fórmula (1)
5. calcular el siguiente precio V de la opción mediante $V(S,t)=e^{-r(T-t)}F(S,t)$
6. calcular la media V_Monte_Carlo de los precios de las opciones
7. repita los pasos 2-6 hasta que se hayan calculado todos los valores

Sin embargo, como se puede ver en la imagen inferior, la curva de los precios de las opciones obtenida con la simulación MC es suave, mientras que lo que yo esperaba era una curva similar a la de la solución exacta. Así que cometí algún error en alguna parte

% S = oil spot prices
S = [ 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 ];
r = .1;    % yearly instantaneous interest rate
T = 1/2;   % expiry time

alpha = 0.0692;
sigma = 0.0876; % values estimated from data
mu = 3.0582;

t = linspace(0,T,numel(S));
tau = T-t;
lambdahat = (mu-r)/alpha;
muhat = mu-sigma^2/2/alpha-lambdahat;

%% Is this V the exact solution?
% this F is the formula (1) on stackexchange
F = exp( exp(-alpha*tau).*log(S) + muhat*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) );
% K = F(end); % shouldn't we use also the strike price?
V_exact = exp( -r*tau ) .* F; % V contains the prices of the option
plot(t,S)
hold on
plot(t,V_exact)

%% Monte Carlo simulation
n = numel(S);
N = 10000; % number of samples
mustar = mu-sigma^2/2/alpha;
X = zeros(N,n);
X(:,1) = log(S(1));
S = zeros(N,n);
S(:,1) = S(1);
F = zeros(N,n);
V = zeros(N,n);
V_Monte_Carlo = zeros(1,n);
dt = 1; % which values for dt has to be pick?
for k = 2:n
    tau = T-t(k);
    % this X is the formula (2) on stackexchange
    X(:,k) = X(k-1)*exp(-alpha*dt) + mustar*(1-exp(-alpha*dt)) + sigma*sqrt(1-exp(-2*alpha*dt))/sqrt(2*alpha)*randn(N,1);
    S(:,k) = exp(X(:,k));
    F(:,k) = exp( exp(-alpha*tau).*log(S(:,k)) + muhat*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) );
    V(:,k) = exp(-r*tau)*F(:,k);
    V_Monte_Carlo(k) = mean(V(:,k));
end
plot(t,V_Monte_Carlo,'g')
legend('Spot prices (S)','Option prices (V) from exact solution','Option prices from Monte-Carlo simulation')

Answer 1

Veamos.

Se tiene la siguiente SDE para el precio de las acciones bajo la medida $P$ :

$$dS(t) = \alpha \cdot (\mu - \log S(t)) \cdot S(t) \cdot dt + \sigma \cdot S \cdot dW^P(t),$$

con condición inicial $S(0) = S_0$ . Además, definir $X(t) = \log S(t)$ suponiendo un precio de mercado constante del riesgo $\lambda = \mu - r$ y realizando un cambio de medida, se obtiene el proceso $X(t)$ bajo la medida de riesgo neutro $Q$ como:

$$dX(t) = \alpha \cdot \left( \mu - \frac{\sigma^2}{2\alpha} - \frac{(\mu - r) \cdot \sigma}{\alpha} - X(t) \right) \cdot dt + \sigma\cdot dW^Q(t).$$

Definamos las funciones de deriva y difusión:

$$\begin{aligned} f(X_t, t) &= \alpha \cdot \left( \mu - \frac{\sigma^2}{2\alpha} - \frac{(\mu - r) \cdot \sigma}{\alpha} - X_t \right),\\ g(X_t, t) &= \sigma. \end{aligned}$$

Ahora, basta con aplicar el esquema EulerMaruyama:

$$X(t_i) = X(t_{i-1}) + f(X(t_{i-1}), t_{i-1}) \cdot \Delta t_i + g(X(t_{i-1}), t_{i-1}) \cdot \left[ W(t_i) - W(t_{i-1}) \right],$$

donde $W(t_i) - W(t_{i-1}) = \mathcal{N}(0, \Delta t_i)$ .

No codifico en Python ni en Matlab, pero puedes encontrar esas implementaciones en wikipedia, ici .

Yo recomendaría utilizar el lenguaje de programación Julia para la simulación de ecuaciones diferenciales estocásticas.

Desea calcular el siguiente precio derivado:

$$\begin{aligned} V(0) &= E_t^Q \left[D(0, T) \cdot \left( S(T) - K \right) \right]\\ V(0) &= E_t^Q \left[D(0, T) \cdot \left( \exp \left( X(T) \right) - K \right) \right]\\ \end{aligned}$$

Dado que el factor de descuento es determinista y $K$ es constante:

$$\begin{aligned} V(0) &= D(0, T) \cdot \left( E_t^Q \left[ \exp \left( X(T) \right) \right] - K \right) \end{aligned}$$

Ahora, aproxima la expectativa mediante un método de Monte Carlo. Tenemos que simular el proceso $X(t)$ para $N$ evaluarlo en el momento $T$ y aplicar la función exponencial.

$$E_t^Q \left[ \exp \left( X(T) \right) \right] = \frac{1}{N} \cdot \sum_{n=1}^N \exp \left( X_n(T) \right)$$

El precio derivado es

$$V(0) = D(0, T) \cdot \left[ \frac{1}{N} \cdot \sum_{n=1}^N \exp \left( X_n(T) \right) - K \right]$$

para cualquier $K$ . En particular, si utiliza $K = F(S_0, T)$ obtendrá que el precio derivado $V(0) = 0$ que es lo que ocurre en la realidad, donde $K$ se supone que toma ese valor. Tenga en cuenta que su función $F(S_0, T)$ puede tener un error, hay un $\sigma$ falta multiplicar $\mu - r$ , ver mi código a continuación en la función F definición. De todos modos, la eliminación de ese $\sigma$ tanto en la SDE para $X(t)$ y $F(S_0, T)$ ceder a precios correctos. Pero quería destacar esta diferencia que vi en este en la sección 2.1.

Eso es todo. Aquí tienes el código en Julia:

using StochasticDiffEq
using Parameters
using Statistics
using Printf
using Plots

function drift(u, p, t)
  @unpack , , , r = p
  return  * (( - ^2 / (2) - ( - r) / ) - u)
end

function diffusion(u, p, t)
  @unpack  = p
  return 
end

function F(S0, , p)
  @unpack , , , r = p
  return exp(
    exp(- * ) * log(S0) + ( - ^2 / (2) - ( - r) / ) * (1 - exp(- * )) +
    ^2 / (4 * ) * (1 - exp(-2 *  * ))
  )
end

let
  params = (
     = 0.0692,
     = 0.0876,
     = 0.25,
    r = 0.10
  )

  S0 = 22.93
  X0 = log(S0)
  T  = 1/2
  trajectories = 100_000

  prob = SDEProblem(drift, diffusion, X0, (0., T), params)
  sde = EnsembleProblem(prob)
  sol = solve(sde, SRIW1(), trajectories = trajectories)

  evs = zeros(trajectories)
  for n in 1:trajectories
    evs[n] = exp(sol[n](T))
  end

  # fair value and standard deviation

  r = params.r
  K = sort(push!(collect(15:0.5:25), F(S0, T, params)))
  x = []
  y = []
  for k in K
     = mean(exp(-r * T) .* (evs .- k))
     = stdm(evs, ; corrected = true) / sqrt(trajectories)

    mcapprox = 
    analytic = exp(-r * T) * (F(S0, T, params) - k)

    push!(x, mcapprox)
    push!(y, analytic)

    # @printf("%e\t%e\n", mcapprox, analytic)
  end
  plot(K, x)
  plot!(K, y)
end

Obtengo un precio derivado igual a cero cuando $K = F(S_0, T)$ . Además, obtengo el mismo precio con el método de Monte Carlo que con la solución analítica cuando $K$ ¡es cualquier otra cosa!

Aquí puede consultar un gráfico con los resultados, donde los valores aproximados coinciden con los resultados analíticos, y $V(0)$ es igual a $0$ cuando se espera:

Por supuesto que los precios pueden ser positivos o negativos, depende del valor de la huelga $K$ . Si $K < F(S_0, T)$ el precio de la opción es positivo porque, para ser justos, el tenedor (el que va a recibir $S(T)$ ), debe pagar $V(0)$ por adelantado. Lo contrario ocurre si $K > F(S_0, T)$ .