Fórmula de Breeden-Litzenberger para densidades de riesgo neutrales

Question

Basado en este tema: ¿Cómo derivar la distribución de probabilidad implícita a partir de las volatilidades de B-S?

Estoy tratando de implementar la fórmula de Breeden-Litzenberger para calcular las densidades de riesgo implícitas del mercado neutral para el S&P 500 para algunas fechas de comilla. Los pasos que sigo son los siguientes:

Paso 1: Extraer las call_strikes c_strikes para una madurez T dada y los precios de mercado correspondientes css.

Paso 2: Una vez que tengo las huelgas y los precios de mercado, calculo las volatilidades implícitas a través de la función ImplieVolatilities.m. Estoy 100% seguro de que esta función funciona.

Paso 3: A continuación, interpolar la curva de volatilidad implícita utilizando el comando Matlab interp1. Ahora, el vector xq es la cuadrícula de huelgas y el vector vq las volatilidades implícitas interp oladas correspondientes.

Paso 4: Dado que tengo la curva de volatilidad, puedo calcular la densidad implícita del mercado: \begin{align} f(K) &= e^{rT} \frac{\partial^2 C(K,T)}{\partial K^2} \\ &\approx e^{rT} \frac{C(K+\Delta_K,T)-2C(K,T)+C(K-\Delta_K, T)}{(\Delta_K)^2} \end{align} donde $\Delta_K = 0.2$ es el tamaño de la cuadrícula de huelga. Solo muestro la parte principal del programa de Matlab:

c_strikes = Call_r_strikes(find(imp_vols > 0));
css = Call_r_prices(find(imp_vols > 0)); 

impvolss = ImpliedVolatilities(S,c_strikes,r,q,Time,0,css); %volatilidades implícitas
 xq = (min(c_strikes):0.2:max(c_strikes)); %cuadrícula de huelgas
vq = interp1(c_strikes,impvolss,xq); %valores interpolados %volatilidades implícitas interpoladas

    f = zeros(1,length(xq)); %densidad neutral al riesgo
function [f] = secondDerivativeNonUniformMesh(x, y)
dx = diff(x); %tamaño de la cuadrícula (uniforme)
dxp = dx(2:end);
d2k = dxp.^2; %tamaño de cuadrícula al cuadrado
f = (1./d2k).*(y(1:end-2)-2*y(2:end-1)+y(3:end)); 
end

figure(2)
plot(xq,f);

Nota: blsprice es una función interna de Matlab por lo que definitivamente debería funcionar. Esta es la pdf neutral al riesgo para los puntos interiores:

Los datos La siguiente tabla muestra los datos de las opciones (columna izquierda: huelgas, medio: precios de mercado y columna derecha: volatilidades implícitas):

                   650                     387.5                  0.337024
                   700                    346.45                  0.325662
                   750                     306.8                  0.313846
                   800                    268.95                  0.302428
                   850                    232.85                  0.290759
                   900                    199.15                  0.279979
                   950                     168.1                  0.270041
                   975                     153.7                  0.265506
                  1000                     139.9                  0.260885
                  1025                    126.15                  0.255063
                  1050                     112.9                  0.248928
                  1075                     100.7                  0.243514
                  1100                     89.45                  0.238585
                  1125                     79.25                  0.234326
                  1150                      69.7                  0.229944
                  1175                      60.8                   0.22543
                  1200                      52.8                  0.221322
                  1225                      45.8                   0.21791
                  1250                      39.6                   0.21486
                  1275                      33.9                   0.21154
                  1300                     28.85                  0.208372
                  1325                      24.4                  0.205328
                  1350                      20.6                  0.202698
                  1375                      17.3                  0.200198
                  1400                     13.35                  0.193611
                  1450                      9.25                  0.190194
                  1500                      6.55                  0.188604
                  1550                       4.2                  0.184105
                  1600                      2.65                  0.180255
                  1650                     1.675                  0.177356
                  1700                     1.125                  0.176491
                  1800                     0.575                  0.177828

$T = 1.6329$ es el tiempo hasta el vencimiento, $r = 0.009779$ la tasa de interés libre de riesgo y $q = 0.02208$ el rendimiento por dividendo. El precio spot $S = 1036.2$.

Cálculo de la tasa de interés libre de riesgo y el rendimiento por dividendo a través de la paridad Put-Call $r$ es la tasa de interés libre de riesgo correspondiente al vencimiento $T$ y $q$ es el rendimiento por dividendo. En la implementación numérica, derivamos $r$ y $q$ nuevamente utilizando la paridad Put-Call. Con ese fin, asumimos la siguiente relación lineal: $$f(K) = \alpha K-\beta,$$ donde $\alpha = e^{-rT}$ y $-e^{-qT}S(0) = \beta$, y $f(K) = P(K,T)-C(K,T)$. Las constantes $\alpha$ y $\beta$ se calculan realizando una regresión lineal. En consecuencia, la tasa de interés libre de riesgo $r$ y el rendimiento por dividendos $q$ se obtienen entonces por \begin{align} r &= \frac{1}{T}\ln\left(\frac{1}{\alpha}\right), \nonumber \\ q &= \frac{1}{T} \ln \left(\frac{-S(0)}{\beta}\right). \nonumber \end{align>

Answer 1

Supongo que para aproximar la segunda derivada del precio de la llamada $C (K,T)$ en los límites del dominio de la huelga (ver los primeros 2 casos "if" del último bucle for de tu código) intentaste establecer condiciones límite.

En el límite derecho, tu aproximación $C(K+\Delta K, T) \approx 0$ podría tener sentido para $K \to \infty$ o al menos lo suficientemente grande.

En el límite izquierdo, sin embargo, $C(K-\Delta K, T)$ no se puede aproximar a cero, ya que el precio de la llamada es casi solo un valor intrínseco para huelgas bajas. En lugar de cero, podrías usar la diferencia descontada del precio a plazo menos la huelga, $C(K-\Delta K, T) \approx B(0,T)(F (0,T)-K)$, si $K \to 0$ o al menos lo suficientemente pequeño.

Esto explica por qué observas un comportamiento extraño en el LHS (con una pdf negativa), mientras supongo que si te acercas o eliminas ese valor atípico deberías estar bien.

Como primer paso, podrías simplemente calcular la pdf para los puntos interiores xq(2:end-1) (es decir, el último caso de tu bucle for) y ver qué da.

[Editar]

Mirando más de cerca, parece que tus precios de entrada permiten oportunidades de arbitraje estático, más específicamente arbitraje de mariposa. Como tal, no puedes esperar obtener una pdf razonable.

Aunque está lejos de ser perfecto, tal vez el script a continuación te ayude. En él, intento buscar las huelgas que conducen a arbitrajes y descartar las entradas correspondientes para obtener lo que se llama una pdf "limpia". Este es un enfoque no paramétrico. Como se discute en los comentarios, utilizar un enfoque paramétrico (por ejemplo, ajustar la sonrisa de volatilidad utilizando una representación sin arbitraje), producirá una pdf aún más suave.

function [ ] = someBreedenLitzenbergerCode( )

% Obtener datos e interpolar
[ K, C, ~ ] = getData();
dk = 0.5;
k = min(K):dk:max(K);
c = interp1(K, C, k, 'spline');

% Notificar oportunidades de arbitraje vertical
dC = diff(C);
if any(dC > 0)
    warning('Los precios de compra de la opción permiten arbitraje (spread vertical)'); 
end

% Notificar y tratar oportunidades de arbitraje de mariposas
f = secondDerivativeNonUniformMesh(k, c);
idxkArb = find(f < 0) + 1;                                     
if ~isempty(idxkArb)
    warning('Los precios de compra de la opción permiten arbitraje (mariposa)');
end

% Encontrar índices de huelgas/precios de entrada que llevan a arbitraje
idxKArb = [];
for i = 1:length(idxkArb)
    idxHi = find(K >= k(idxkArb(i)), 1, 'first');
    idxLo = find(K <= k(idxkArb(i)), 1, 'last');
    idxKArb = [idxKArb, idxHi, idxLo];
end
idxKArb = unique(idxKArb);
disp('Huelgas que llevan a arbitraje de mariposa');
disp(K(idxKArb));

% Limpiar entradas y calcular pdf
K_clean = K; K_clean(idxKArb) = [];
C_clean = C; C_clean(idxKArb) = [];
k_clean = min(K_clean):dk:max(K_clean);
c_clean = interp1(K_clean, C_clean, k_clean, 'spline');
f_clean = secondDerivativeNonUniformMesh(k_clean, c_clean);

% Graficar resultados
if ishandle(1), clf(1); end;
figure(1);
subplot(221);
plot(K, C, 'x'); hold on;
plot(k, c, 'r-');
plot(k(idxkArb), c(idxkArb), 'rd');
grid on;
title('Precios originales de compra, C(K)');
subplot(222);
plot(k(2:end-1), f); hold on;
plot(k(idxkArb), f(idxkArb - 1), 'rd');
grid on;
title('Pdf, \phi(K)');
subplot(223);
plot(k_clean, c_clean); hold on; 
lh.Box = 'off';
grid on;
title('Precios de compra limpios, C(K)');
subplot(224);
plot(k_clean(2:end-1), f_clean); hold on;
a = ksr(k_clean(2:end-1), f_clean, 60); % Disponible en Matlab Central
plot(a.x, a.f, 'k-');
lh = legend('FD','FD + Suavizado');
lh.Box = 'Off';
grid on;
title('Pdf, \phi(K)'); 

end

function d2y_dx2 = secondDerivativeNonUniformMesh(x, y)

dx = diff(x);
dxp = dx(2:end);
dxm = dx(1:end-1);
d2k = dxp .* dxm;
d2kp = dxp .* (dxm + dxp);
d2km = dxm .* (dxm + dxp);
d2y_dx2 =  2 ./ d2km .* y(1:end-2) - 2 ./ d2k .* y(2:end-1) + 2 ./ d2kp .* y(3:end);

end

function [strikes, prices, volatilities] = getData()

mat = [650        387.5      0.33702
      700       346.45      0.32566
      750        306.8      0.31385
      800       268.95      0.30243
      850       232.85      0.29076
      900       199.15      0.27998
      950        168.1      0.27004
      975        153.7      0.26551
     1000        139.9      0.26088
     1025       126.15      0.25506
     1050        112.9      0.24893
     1075        100.7      0.24351
     1100        89.45      0.23858
     1125        79.25      0.23433
     1150        69.7      0.22994
     1175        60.8      0.22543
     1200        52.8      0.22132
     1225        45.8      0.21791
     1250        39.6      0.21486
     1275        33.9      0.21154
     1300        28.85      0.20837
     1325        24.4      0.20533
     1350        20.6      0.2027
     1375        17.3      0.2002
     1400        13.35      0.19361
     1450         9.25      0.19019
     1500         6.55      0.1886
     1550         4.2      0.1841
     1600         2.65      0.18025
     1650         1.675      0.17736
     1700        1.125      0.17649
     1800        0.575      0.17783];
 strikes = mat(:,1);
 prices = mat(:,2);
 volatilities = mat(:,3);

end