I was trying to reproduce the partial non-parametric garch model from paper on page 1764 https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1540-6261.1993.tb05127.x
The problem is how to derive the likelihood function? For example, I divide the range of $\varepsilon_t$ into 6 intervals using unconditional variance $\sigma$ calculated as sample variance. So far I derived the likelihood function is the following.
$L(\omega, \beta, \theta_0,\theta_1,\theta_2, \delta_0,\delta_1,\delta_2| \varepsilon) = \prod_{t=1}^{n}\frac{1}{\sqrt{h_t}}f_z(\frac{\varepsilon_t-\mu}{\sqrt{h_t}})$
$$ h_t = \left\{ \begin{array}{ll} A+\delta_0 \varepsilon_{t-1} + \delta_1(\varepsilon_{t-1}+\sigma) + \delta_2(\varepsilon_{t-1}+2\sigma)& \quad \varepsilon_{t-1} \leq -2\sigma \\ A+\delta_0 \varepsilon_{t-1} + \delta_1(\varepsilon_{t-1}+\sigma) & \quad -2\sigma < \varepsilon_{t-1} \leq -\sigma \\ A+\delta_0 \varepsilon_{t-1}& \quad-\sigma< \varepsilon_{t-1} \leq 0 \\ A+\theta_0 \varepsilon_{t-1}& \quad0< \varepsilon_{t-1} \leq \sigma \\ A+\theta_0 \varepsilon_{t-1} + \theta_1(\varepsilon_{t-1}-\sigma) & \quad \sigma < \varepsilon_{t-1} \leq 2\sigma \\ A+\theta_0 \varepsilon_{t-1} + \theta_1(\varepsilon_{t-1}-\sigma) + \theta_2(\varepsilon_{t-1}-2\sigma)& \quad \varepsilon_{t-1} >2\sigma \\ \end{array} \right. $$ where $A = \omega + \beta\sigma^2$
The likelihood function is a piecewise function, how can I fit this model using R? Is there any way to change this likelihood function into a simplified version? The author used QML by BHHH algorithms, how can I implement that? I think I can coup with that if it's not piecewise function but how to change it?