I am currently working on option pricing model and I'd like to include a method for maximizing the likelihood of returns under the P measure. I am using the Heston and Nandi (2000) model: \begin{align} ln S_{t+1} - ln S_t := r_{t+1} &= r_{ft+1} + \lambda h_{t+1} - \xi_{t+1} + \sqrt{h_{t+1}} z_{t+1}, \; z_{t+1} \sim N(0,1) \\ h_{t+1} &= \sigma^2 + \pi \left( h_t - \sigma^2 \right) + \alpha \left(z_t^2 - 1 - 2 \gamma \sqrt{h_t} z_t \right). \end{align}

Above, the frequency of the data will be daily. Moreover, $\xi_{t+1}$ is a convexity correction which ensures that expected gross returns $E_t(S_{t+1}/S_t) = E_t(\exp r_{t+1}) = \exp(r_{ft+1} + \lambda h_{t+1})$. Since $z_{t+1} \sim N(0,1)$, the logarithm of the conditional moment generating function is $\xi_{t+1} = h_{t+1}/2$.


The first thing I thought of doing for stabilizing the estimation is to make sure $h_{t+1}$ lies within certain bounds. I impose that at all times, $h_{t+1} > 0.01^2/N_{days}$ (i.e., I am excluding the possibility of seeing days with annualize volatility below 1%). Since I am working with the S&P500, I suppose it's not crazy. I also impose that it cannot be higher than 5 (i.e., annualized volatility in a day cannot exceed 500%). It's not crazy, especially since my sample stops in 2013. I enforce it directly in the filtration step of the optimization:

        for tt in range(0,T-1):
            z[tt]   = ( series[tt] - (lambda_-0.5)*h[tt] )/sqrt(h[tt])
            h[tt+1] = sigma2 + persistence*(h[tt] - sigma2) + alpha*(z[tt]**2 - 1 - \

        # To ensure smooth optimization, enforce bounds on h(t+1):
        h[tt+1] = max(self.h_min, min(h[tt+1], self.h_max))

And, obviously, I had a flag to tell me if I enforced the bounds.

The other thing I am doing is that I follow the literature in estimate $\sigma^2$ using the full sample and outside of the MLE: $\hat{\sigma}^2 := \frac{1}{T-1} \sum_{t=1}^{T} \left( r_{t+1} - \bar{r} \right)^2$. It's called "variance targeting" and it is typical in GARCH option pricing papers. The last thing I would do is enforce bounds on $\pi$, specifically $|\pi| < 1$. I could also put bounds on $(\alpha, \gamma)$ using previous results from the literature, but I am not sure it's going to entirely necessary. I think this should make sure nothing crazy happens, but if you have comments or other suggestions, I am very much open to them.


Now, for that part, I have no idea where to start. I suppose that if impose bounds on everything, I could pick a bunch of random points using uniform random variables and go with the solution which works best among them. I could also look at previous work and initialize at or close to their estimates.

Here, I would really appreciate some pointers for best practice.


Just so you know where this is going, the idea is to calibrate the HN2000 model for pricing European options on the S&P500 using a joint likelihood. What you see above is the P-measure part. The Q-measure part would use Q-dynamics to produce prices that would in turn be expressed as implied volatilities. The Q-likelihood would be a gaussian likelihood in the volatility surface, in other words.

So, you're looking at step 1 here and I need to make sure this works well before moving on to step 2. Thanks in advance.


I recommend you take a look at Christoffersen, Heston Jacobs (2013), since they conduct a joint and a sequential analysis and more importantly they include a non-linear pricing kernel which you can easily integrate in the HN2000 model.

Regarding $h_0$ and $h_t$:

For $h_t$ an initial value needs to be set. Heston and Nandi (2000) take the variance of the entire sample as a starting value, however, I found that only taking the variance of the previous let's say 15 trading days works very well too. This comes with a number of benefits. In a long return sample (which is what you want for a GARCH analysis), iterating all the way back to the starting value at every observation is computationally costly.

Therefore, at each point in time in the estimation, only the seven preceding values are considered, with the seventh not relying on the eight, but instead on the variance of the 15 previous days. This way 22 trading days of historical information are considered at every step of the optimization routine, which corresponds to roughly one month of time. Of course you can also go back further and you should analyse if that makes a significant difference for your data and your case.

In my case the log-likelihood did in fact increase with every previous period added, however, less so at each step. While this appears to be a restrictive assumption it turns out, results are very much in line with other empirical studies. A possible explanation is that the initial value implied by the 15-day variance is usually much closer to the level of conditional variance than the sample variance, therefore requiring less periods to converge. Furthermore, I found the starting value of the conditional variance, $h_0$, exerts a negligible influence on the results. As Heston and Nandi (2000) note, this is caused by the strong mean-reversion property of variance.

Regarding the return parameter ($\mu$):

Estimating the return parameter (in my case $\mu$) without any restriction leads to widely fluctuating values, even reaching negative territory at some points. Giving it a reasonable range of starting values avoids this issue without imposing any limitations. Moreover, as CHJ (2013) note, realistic values for $\mu$ can quickly be inferred from market data. An annual equity risk premium given by $\mu h_t$ = 10% and an annualized volatility of about 21% implies $\mu$ = ~2.26. These are approximate average return and volatility values for the years of the sample.

Two notes on this:

  • CHJ (2013) apparently impose $\mu$ = 0 after having obtained the original HN (2000) model parameters in their sequential estimation, without giving an explicit reason.
  • The p-value of $\mu$ is rather high and in the final results, a higher or lower parameter, even in the order of magnitude of 50%, is barely noticeable. So setting it to zero or giving it a reasonable starting value and restricting it should be a good approach.

Regarding $\omega$:

To a lesser extent, the same holds for $\omega$. Some authors, such as Byun (2011) and CHJ (2013), obtain slightly negative or very small values for $\omega$ and therefore set it equal to zero. While this lowers the computation time it lead to significantly lower log-likelihoods in my analysis. I propose to estimate $\omega$ freely and provide the alternative suggestion of setting it to a reasonably small number in the region of $1*10^-7$ to $1*10^-6$. This range is supported by results of a large number of reputable works (e.g., Heston and Nandi (2000), Christoffersen and Jacobs (2004) and Christoffersen et al. (2012)), as well as our own results and these parameters explains between $1.6$% and $5$% of annual conditional standard deviation.

Combining this insight with the suggestion to infer $\mu$ from market data, not just to check for plausibility, but to omit it from estimation, reduces the calibration problem from five to three parameters, allowing you to spend more time finding better values for the remaining parameters.

The remaining parameters $\alpha$, $\beta$ and $\gamma$:

These parameters, especially $\beta$ have a big impact on the results. It is therefore crucial to estimate them accurately and precisely, since any deviation has comparably large implications.

| improve this answer | |
  • 1
    $\begingroup$ I read CHJ2013 and it is indeed a good reference. I could use $\bar{r}_{ft}$ and $\bar{r}_t$, as well as $\hat{\sigma}^2$ as a way to infer a starting value for $\lambda$. For $\pi, \alpha$ and $\beta$, I could simply look at previous papers that use this model -- CDJW2010, CHJ2013, BCHJ2018, etc. But you have a good proposal with matching data moments. $\endgroup$ – Stéphane Apr 21 at 16:12
  • $\begingroup$ I am also wondering: the risk-free rate I would use for pricing options would in principle change across time, as well as across maturity. For the return likelihood under P, I am thinking I should be using excess returns as imputs $lnS_{t+1}-ln S_t - r_{ft+1}$, but this means picking a single interest rate series. Which one should I use? I mean, we're basically "cheating" here by not modeling the short rate process, but I need to take stock of changes and of its term structure. $\endgroup$ – Stéphane Apr 21 at 18:38
  • $\begingroup$ I agree, you could use an interest rate that is similar to your forecast horizon. $\endgroup$ – Andreas Apr 22 at 9:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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