Duan (1995) GARCH Option Pricing Model with MATLAB

Question

This is the MATLAB code that replicates the option pricing model proposed by Duan in his paper "The GARCH Option Pricing Model". However, the parameters estimated in the file do not match with the ones presented in the paper. I tried to fix it but I still get wrong parameter values.

Here's the .m file for the Maximum Likelihood Estimation:

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);

Here's the main .m file:

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);

For completeness of information, the estimated parameter values in the paper are $\alpha_0=1.524\times 10^{-5}$, $\alpha_1=0.1883$, $\beta_1=0.7162$ and $\lambda=7.452\times 10^{-3}$, and the standard deviation $\sigma=24.13\%$. Moreover, the model is fitted to the S&P 100 daily index from Jan 2, 1986 to Dec 15, 1989.

Can someone help me to find where my error lies? Thank you very much for your help.

Answer 1

Your function returning (minus) the log-likelihood seems weird to me, I would go with

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);

That being said, if everything is OK but your input prices are wrong you'll never get the right parameters, so start by checking this.

REM: I didn't try to run it, so I might have left some typos but you should get the idea.

---

Edit

I made some typos indeed (corrected now hopefully) but I assumed you are working with the following dynamics under the physical measure \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}
with $\mathbf{z}=(z_1,\dots,z_N)$ i.i.d. Gaussian variates.

So the first step is to recover the hidden conditional variances from the data series. This is done by inserting the first equation above in the second yielding the following recursion: $$ 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}$$ starting from $h_1 = \text{var}(\mathbf{r})$ with $\mathbf{r}=(r_1,\dots,r_N)$ your sample log-returns.

From there onwards, you have that, conditionally on $h_t$ $$ R_t \sim \mathcal{N}(r_F - \frac{1}{2}h_{t}, h_t) $$ Assuming you have an i.i.d. sample of $N$ returns you then form the log-likelihood as

\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}

hence finally the negative log-likelihood \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}

where, $\mathbf{r}=(r_1,\dots,r_N)$ are the observed log-returns and the instantaneous variances $\mathbf{h}=(h_1,\dots,h_N)$ are obtained as discussed above.

Now, see equation (EQ1) and (EQ2) in the above code.