I've searched thru dozens of papers and did not find in any of them satisfying and enough theoretical answers to my concerns. So I've combined everything what I found below. Please indicate if my understanding of the topic is proper and correct me if necessary. Beside theoretical side I've found also an obstacle in practical implementation.

Duan (1995) in his paper developed model for European options pricing with GARCH model. He introduced locally risk neutral measure (LRNVR) $\mathbb{Q} $ which is equivalent to physical market measure $\mathbb{P} $. Since price for either call or put option has no analytical solution in his framework it is necessary to run Monte Carlo simulations. My first concern for which I did not found explicit answer is the following:

(1) We have a time serie with some market data and we fit parameters of the model to this time serie under the physical measure $ \mathbb{P} $.

(2) Then we use estimated parameters to the transformed process under the LRNVR measure $\mathbb{Q} $ and we run Monte Carlo simulations to estimate price of an option.

Since MC simulation under measure $\mathbb{Q} $ is not a problem for me, let's now focus on theoretical derivation of MLE under physical measure $\mathbb{P} $.

Suppose we have a sample of $T$ log-returns of some financial asset. Let $X = (X_1, \dots , X_T) $ denote our sample and assume that $t$-th log-return follows GARCH(1,1)-M process proposed by Duan (1995): $$ X_t = \ln \frac{S_t}{S_{t-1}} = r + \lambda \sqrt{h_t} - 0.5h_t + \sqrt{h_t} z_t $$ under physical measure $\mathbb{P} $, where $z_t \overset{iid}{=} \mathcal{N} (0,1) $ and $h_t = \omega + \alpha h_{t-1} z_{t-1}^2 + \beta h_{t-1} $. Here we assume that $ \omega >0 $ and $ \alpha, \beta \geq 0 $. We would also require $ \alpha + \beta < 1$ to ensure stationarity. Parameter $r$ is a market risk-free interest rate (known parameter) and $\lambda$ is a risk-premium associated with a given financial asset (parameter to be estimated).

Since we have $ X_t \sim \mathcal{N} (r + \lambda \sqrt{h_t} - 0.5h_t, h_t )$, then likelihood function for $t$-th observation is: $$ l_t (X_t ; \theta ) = \frac{1}{\sqrt{2 \pi h_t}} \exp \left( - \frac{ ( X_t - r - \lambda \sqrt{h_t} + \frac{1}{2} h_t )^2 }{2h_t} \right) \text{,} $$

where $\theta = (\omega, \alpha, \beta, \lambda ) $ is a vector of parameters to be estimated. Likelihood function for a vector $X$ is: $$ \mathcal{L} (X; \theta ) = \prod_{t=1}^T l_t(X_t; \theta) \text{.} $$ Since it is easier to compute natural logs, then: $$ \ln \mathcal{L} (X; \theta ) = -\frac{T}{2} \ln \left( 2 \pi \right) - \frac{1}{2} \sum_{t=1}^T \ln \left( h_t \right) - \frac{1}{2} \sum_{t=1}^T \frac{ ( X_t - r - \lambda \sqrt{h_t} + \frac{1}{2} h_t )^2 }{2h_t} $$ to be maximized. We are seeking for $ \hat{\theta} $ which is: $$ \hat{\theta} = \arg \max_{\theta} \ln \mathcal{L} (X; \theta ) $$ with constrains for $\omega , \alpha, \beta $ as indicated above.

Now let's move to the practical implementation of the above:

Let our market data will be AAPL daily log-returns for the period 2016-2019 (or any other data, because the following problem does not disappear with the change of source data). Suppose risk-free interest rate is $r=0$. As initial variance $h_1$ we assume variance of our sample, i.e. $h_1 = Var(X) $. Function to be minimized is $ - \ln \mathcal{L} (X; \theta ) $ and is defined as follows:

loglike <- function(params, log_returns){

  omega <- params[1]
  alpha <- params[2]
  beta <- params[3]
  lambda <- params[4]

  bigT <- length(log_returns)  
  h <- c(var(log_returns))

  for (i in 2:bigT) {
    h[i] <- omega + alpha * (log_returns[i-1] - lambda * sqrt(h[i-1]) + 0.5 * h[i-1] )^2 + beta * h[i-1]

  likelihood <- 0.5 * bigT * log(2*pi) + 0.5 * sum(log(h) + ((log_returns - lambda * sqrt(h) + 0.5 * h)^2) / h )

Starting parameters and constrains for optimization problem are as follows:

params <- rep(0.01, 4)
lb <- c(0,0,0, -Inf)
A = cbind(0, 1, 1, 0)
b = 1

And if I try to use fmincon from pracma package for optimization I have the following error and warnings:

> fmincon(loglike, log_returns = log_returns, x0 = params, lb = lb, A = A, b = b)
Error in if (f < 0) { : missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In sqrt(h[i - 1]) : NaNs produced
2: In log(h) : NaNs produced
3: In sqrt(h) : NaNs produced

First of all I do not know why I get these warnings. According to the objective function defined above $ h \geq 0 $ so taking square roots and natural logs should not produce NaNs.

Secondly I do not understand the error returned by fmincon. What is wrong with my objective function?

  • $\begingroup$ A hint: In a standard GARCH(1,1)-model, the initial variance $h$ is commonly set to something like $h=\alpha/(1-\beta-\gamma)$ instead of $var(x)$. $\endgroup$ May 4, 2022 at 6:08

1 Answer 1


There are two things i recommend you to try, maybe that already fixes your problem.

First, drop the constant term $-\frac{T}{2}ln\left ( 2\pi \right )$ to eliminate one possible source for your error. (this part it is not needed for the optimization anyways)

Second, sometimes errors arise due to numerical issues when optimization routines do not strictly stick to the limits set by the user (even if they should). You can try setting your limit to an arbitrarily small positive number.

  • $\begingroup$ Thanks for reply. However the problem does not disappear. Even if I remove optimization boundaries, there is the same error. $\endgroup$ Apr 15, 2020 at 15:55
  • $\begingroup$ Thee might be a misunderstanding. I meant to suggest that you set the optimization boundaries not to zero but to a very small number such as 1e^-6, or whatever is reasonable for your parameters. $\endgroup$
    – Andreas
    Apr 15, 2020 at 18:03
  • $\begingroup$ I did it already. It does not help. $\endgroup$ Apr 15, 2020 at 18:49
  • 1
    $\begingroup$ Hi: I'm not saying that this will solve your problem but, when you have positivity constraints, it might be helpful take the exp of the variables so that they can't go negative. So try that and see if it makes a difference. Also, I'm not familiar with fmincon but it's usually useful to supply the gradient if the algorithm is a quasi-newton type procedure. I've had good experience Rvmmin in R. Finally, you might be better off sending that to the R-list or R stackoverflow. John Nash, the author of the of the majority of the optimization routines in R ( and generous too ), is on the R-list. $\endgroup$
    – mark leeds
    May 13, 2020 at 21:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.