Duan (1995) Modelo de valoración de opciones GARCH con MATLAB

Question

Este es el código MATLAB que reproduce el modelo de valoración de opciones propuesto por Duan en su artículo "The GARCH Option Pricing Model". Sin embargo, los parámetros estimados en el archivo no coinciden con los presentados en el documento. Intenté arreglarlo pero sigo obteniendo valores de parámetros erróneos.

Aquí está el archivo .m para la estimación de máxima verosimilitud:

function y = findGARCH_LLy(params,S,rf)

% Finds log-likelihood for the GARCH option pricing model.
alpha0 = params(1);
alpha1 = params(2);
beta1  = params(3);
lambda = params(4);

N = length(S);

% Define the returns.
r = log(S(2:N)./S(1:N-1));
r = [0; r];

% Initialize the GARCH log-likelihood at t=1.
h(1) = var(r);
e(1) = 0;
LL(1) = -0.5*log(2*pi)-0.5*log(h(1))-0.5*e(1)^2/h(1);

% Find the rest of the GARCH log-likelihood.
for t=2:N
    h(t) = alpha0 + alpha1*(e(t-1) - lambda*sqrt(h(t-1)))^2 + beta1*h(t-1);
    e(t) = log(S(t)/S(t-1)) - rf + 0.5*h(t);
    LL(t) = 0.5*log(2*pi)-0.5*log(h(t))-0.5*e(t)^2/h(t);
end

% Return the negative log-likelihood.
y = -sum(LL);

Aquí está el archivo principal .m:

clc; clear;

% Input the price levels and dates.
[P, Dates] = xlsread('SP 100 Prices.xls','Sheet1');
Prices = P(:,2);

% Risk free rate
rf = 0;

% Starting values for optimization.
% Use estimates from Duan's paper as starting values.
start  = [0.000015 0.19 0.72 0.007];
A = [0 1 1 0];
b = 1;
lb = [0 0 0 0];
ub = [+Inf +Inf +Inf +Inf];
params = fmincon(@(b) findGARCH_LLy(b,Prices,rf), start, A, b, [], [], lb, ub)
alpha0 = params(1);
alpha1 = params(2);
beta1  = params(3);
lambda = params(4);

% Find the standard deviation (sigma) implied by the parameters.
% Assume 365 days per year.
variance = alpha0 / (1 - alpha1*(1+lambda^2) - beta1);
sigma = sqrt(variance)*sqrt(365);

Para completar la información, los valores de los parámetros estimados en el documento son $\alpha_0=1.524\times 10^{-5}$ , $\alpha_1=0.1883$ , $\beta_1=0.7162$ y $\lambda=7.452\times 10^{-3}$ y la desviación estándar $\sigma=24.13\%$ . Además, el modelo se ajusta al índice diario S&P 100 desde el 2 de enero de 1986 hasta el 15 de diciembre de 1989.

¿Puede alguien ayudarme a encontrar dónde está mi error? Muchas gracias por su ayuda.

Answer 1

Su función que devuelve (menos) la log-verosimilitud me parece rara, yo iría con

function y = findGARCH_LLy(params,S,rf)

% Finds log-likelihood for the GARCH option pricing model.
alpha0 = params(1);
alpha1 = params(2);
beta1  = params(3);
lambda = params(4);

N = length(S);

% Define the returns (pad first return with zero)
r = [0, diff(log(S))];

% Infer the other conditional instantaneous variances [see EQ1]
h(1) = var(r);  % Initialize the conditional variance recursion
for i=2:N
   h(i) = alpha0 + alpha1*(r(i-1) - rf + .5*h(i-1) - lambda*sqrt(h(i-1)))^2 + beta1*h(i-1)];
end

% Return the negative log-likelihood (up to constant terms) [see EQ2]
LL = log(h) + 1./h .* (r - rf + .5*h).^2;
y = -sum(LL);

Dicho esto, si todo está bien pero los precios de entrada son erróneos, nunca obtendrás los parámetros correctos, así que empieza por comprobar esto.

REM: No he intentado ejecutarlo, así que puede que haya dejado alguna errata, pero deberíais haceros una idea.

---

Editar

Efectivamente, cometí algunas erratas (ahora corregidas, espero), pero supuse que trabajabas con la siguiente dinámica bajo la medida física \begin {align} \ln { \left ( \frac {S_{t+1}}{S_t} \right )} = R_{t+1} &= r_F - \frac {1}{2}h_{t+1} + \sqrt {h_{t+1}}z_{t+1} \\ h_{t+1} &= \alpha_0 + \alpha_1 ( \sqrt {h_t} z_t - \sqrt {h_t} \lambda )^2 + \beta_1 h_t \end {align}
con $\mathbf{z}=(z_1,\dots,z_N)$ Variables gaussianas i.i.d.

Así pues, el primer paso consiste en recuperar las varianzas condicionales ocultas de las series de datos. Esto se hace insertando la primera ecuación anterior en la segunda dando lugar a la siguiente recursión: $$ h_{t+1} = \alpha_0 + \alpha_1 \left(R_t - r_f + \frac{1}{2}h_{t} - \sqrt{h_t} \lambda \right)^2 + \beta_1 h_t \tag{EQ1}$$ a partir de $h_1 = \text{var}(\mathbf{r})$ con $\mathbf{r}=(r_1,\dots,r_N)$ su muestra de registros.

A partir de ahí, tienes que, condicionado a $h_t$ $$ R_t \sim \mathcal{N}(r_F - \frac{1}{2}h_{t}, h_t) $$ Asumiendo que tienes una muestra i.i.d. de $N$ devuelve entonces la forma de la log-verosimilitud como

\begin {align} p_{R_1,R_2,...,R_N}(r_1,r_2,...,r_N; \Theta ) &= \prod_ {i=1}^N p_{R_i|| \phi_ {i-1}}(r_i; \Theta ) \\ &= \prod_ {i=1}^N \frac {1}{ \sqrt {2 \pi h_i}} \exp \left ( - \frac {1}{2h_i}(r_i - (r_F- \frac {1}{2}h_i))^2 \right ) \\ &:= \mathcal {L}( \Theta ) \end {align}

por lo que finalmente la log-verosimilitud negativa \begin {align} l( \Theta ) &:= - \ln ( \mathcal {L}( \Theta )) \nonumber\\ & \propto \sum_ {i=1}^N \ln (h_i) + \frac {1}{h_i} (r_i - r_F + \frac {1}{2}h_i))^2 \tag {EQ2} \end {align}

donde, $\mathbf{r}=(r_1,\dots,r_N)$ son los rendimientos logarítmicos observados y las varianzas instantáneas $\mathbf{h}=(h_1,\dots,h_N)$ se obtienen como se ha comentado anteriormente.

Ahora, vea la ecuación (EQ1) y (EQ2) en el código anterior.