I want to fit the following model to a time series:
$$ y_{t}=\alpha_{0}+\alpha_{1}y_{t-1}+\alpha_{2}y_{t-1}^{2}+\lambda h_{t}+\varepsilon_{t} $$
$$ h_{t}=\beta_{0}+\beta_{1}\varepsilon_{t-1}^{2}+\beta_{2}h_{t-1} $$
How can I do this with R or with any other statistical software?
Thanks