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$.
CONSTRAINTS
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 - \
2*gamma*sqrt(h[tt])*z[tt])
# 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.
INITIALIZATION
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.
WHERE IT ALL FITS
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.