Say you have a GARCH-M(1,1) model as follows:
$y_t = \beta y_{t-1} + \delta h_t + \epsilon_t, \quad \epsilon_t \sim N(0, h_t) $
$h_t = a_0 + a_1 \epsilon^2_{t-1} + b_1 h_{t-1}.$
How exactly does one estimate the parameter vector: $\Theta = (\delta, \beta, a_0, a_1, b_1)^\top$? This paper, shows the likelihood function for estimating the parameters in the model using the MLE method. For those who have applied this or a similar methodology, are there any nuances to take note of during the estimation process?