Discretización de Euler de la SDE de Heston en Mathematica

Question

A continuación se muestra una implementación de la solución numérica de la EDE de Heston utilizando la discretización de Euler. Tarda menos de un segundo en ejecutarse en Mathematica.

Los parámetros de calibración dan un buen ajuste a la superficie de volatilidad utilizando la técnica de la función característica/transformación de Fourier.

Estoy tratando de utilizar el siguiente código para fijar el precio de un derivado exótico mediante simulación MC, pero soy incapaz de igualar la superficie de volatilidad como primer paso. Sospecho que el código está simplemente tomando demasiado tiempo para converger.

¿Hay alguna solución rápida que pueda acelerar este código de forma significativa?

\[Rho] = -0.4042;
v0 = 0.2577^2;
\[Kappa] = 0.2656;
vbar = .1851;
\[Sigma]v = 0.2992;
n = 100;
NPaths = 100;
Tmax = 574/365;
dt = Tmax/n;
dw = RandomVariate[
BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]], {NPaths, n}];

HestonPaths[G0_] :=
 Module[{XPaths, X, v},
  XPaths = {};
  Do[
   X = {Log[G0]};
   v = {v0};
   Do[
    v = Append[v, 
      Abs[Last[v] + \[Kappa] (vbar - Last[v]) dt + 
        Sqrt[Last[v]] \[Sigma]v dw[[idx, i]][[1]]]];
    X = Append[X, 
      Last[X] - 1/2 Last[v] dt + Sqrt[Last[v]] dw[[idx, i]][[2]]];
    , {i, 1, n}];
   XPaths = Append[XPaths, X];
   , {idx, 1, NPaths}];
  Exp[XPaths]
  ]

ListLinePlot[HestonPaths[500]]

Answer 1

Cambiar la discretización y utilizar el enfoque QE-M: Andersen (2006) el sesgo es mucho menor que el del Euler simple. además, se puede intentar utilizar variantes de control/números anatómicos para reducir la varianza de la muestra.

Answer 2

Algunas mejoras sencillas:

1) Sustituir la aproximación de discretización de Euler de la volatilidad por una aproximación de discretización de Milstein. Véase, por ejemplo estas notas por Rouah.

2) 100 caminos es un número muy bajo de caminos, y conduce a un gran error estándar en su estimación. Por lo tanto, debería multiplicarse por un factor de ~100.

3) Debe utilizar alguna forma de reducción de la varianza. Las variables antitéticas son fáciles de implementar y proporcionan una gran mejora en su error estándar.

Como mencionó Phun, hay una serie de enfoques más complicados que podrías utilizar para construir tus rutas. Estos reducen su sesgo, ya que estos enfoques pueden normalmente evitar que la varianza sea negativa.

Pero yo sugeriría probar primero las simplificaciones anteriores, porque el sesgo no debería ser extremadamente grande. Con este esquema deberías poder acercarte bastante a las volatilidades implícitas determinadas por el método de precios de Fourier.

Y, por último, algunas comprobaciones de cordura: ¿calibraste el modelo Heston utilizando también una tasa cero? ¿Exponen los HestonPaths finales para obtener los Heston Paths "reales"? (ya que está aproximando $d\log[S]$ ).

Answer 3

Gracias por las respuestas. Todavía estoy desconcertada con esto.

Aquí hay una implementación que utiliza tanto la discretización de Milstein como las variables antitéticas.

El código construye el sesgo de la volatilidad para una opción de compra de T = 574 días con un precio inicial a plazo G0 = 570,856 y un tipo de interés r = 0,05327. El sesgo de la volatilidad es claramente disparatado: es cóncavo con un sesgo positivo mientras que el modelo de Heston debería dar una sonrisa convexa y un sesgo negativo.

He puesto el número de incrementos a ser n = 100 y NPaths = 100 sólo para demostrar el código. Aumentar a n = 1000 incrementos con NPaths = 10000 lleva una hora de ejecución y no parece mejorar el resultado.

\[Rho] = -0.4042;
v0 = 0.2577^2;
\[Kappa] = 0.2656;
vbar = .1851;
\[Sigma]v = 0.2992;
n = 100;
NPaths = 100;

HestonPaths[G0_, T_] :=
 Module[{XPaths, X1, v1, X2, v2, dt, dw},
  dt = T/n;
  dw = RandomVariate[
    BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]], {NPaths, n}];
  XPaths = {};
  Do[
   X1 = {Log[G0]};
   v1 = {v0};
   X2 = {Log[G0]};
   v2 = {v0};
   Do[
    v1 = Append[v1, 
      Abs[\[Kappa] (vbar - Last[v1]) dt - 
        1/4 \[Sigma]v^2 dt + (Sqrt[Last[v1]] + 
          1/2 \[Sigma]v dw[[idx, i]][[1]])^2]];
    X1 = Append[X1, 
      Last[X1] - 1/2 Last[v1] dt + Sqrt[Last[v1]] dw[[idx, i]][[2]]];
    v2 = Append[v2, 
      Abs[\[Kappa] (vbar - Last[v2]) dt - 
        1/4 \[Sigma]v^2 dt + (Sqrt[Last[v1]] - 
          1/2 \[Sigma]v dw[[idx, i]][[1]])^2]];
    X2 = Append[X2, 
      Last[X2] - 1/2 Last[v2] dt - Sqrt[Last[v2]] dw[[idx, i]][[2]]];
    , {i, 1, n}];
   XPaths = Append[XPaths, X1];
   XPaths = Append[XPaths, X2];
   , {idx, 1, NPaths}];
  Exp[XPaths]
  ]

T = 574/365;
G0 = 570.856;
r = 0.05327;

Paths = HestonPaths[G0, T];

CallPrice[K_] :=
 Block[{CF, \[CurlyPhi]VanillaCall, G},
  G = Table[Part[Paths, i, n], {i, 1, NPaths}];
  \[CurlyPhi]VanillaCall[S_] := Max[0, S - K];
  CF = Map[\[CurlyPhi]VanillaCall, G];
  Exp[-r T] Mean[CF]
  ]

BCallIV[G_, K_, r_, T_, value_] :=  
  FinancialDerivative[{"European", "Call"}, {"StrikePrice" -> K, 
    "Expiration" -> T, "Value" -> value}, {"InterestRate" -> r, 
    "CurrentPrice" -> G, "Dividend" -> r}, "ImpliedVolatility"
   ];

ListPlot[Table[{k, BCallIV[G0, k, r, T, CallPrice[k]]}, {k, 300, 700, 
   10}]]

Answer 4

Parece que obtengo mejores resultados si reemplazo

BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]]

con

BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho] dt]

Esto es muy extraño porque la documentación de Mathematica indica claramente que BinormalDistribution toma la correlación $\rho$ como entrada, mientras que la covarianza viene dada por $\rho \, dt$ .

Intentaré ejecutar este código modificado con $10^4$ rutas de muestreo y $10^4$ incrementos de tiempo.

¿Es esto excesivo para obtener estadísticas decentes?