I would fit an AR(1)-GARCH(1,1) with arch in mean effect and with the square of the serie at lag 1 as explanatory variable in the mean process :
The AR(1)-GARCH(1,1) component for :
$$
y_{t}=\alpha_{0}+\alpha_{1}y_{t-1}+\varepsilon_{t}
$$
$$
h_{t}=\beta_{0}+\beta_{1}\varepsilon_{t-1}^{2}+\beta_{2}h_{t-1}
$$
+ the arch in mean effect ($\lambda h_{t}$ ):
$$
y_{t}=\alpha_{0}+\alpha_{1}y_{t-1} +\lambda h_{t}+\varepsilon_{t}
$$
$$
h_{t}=\beta_{0}+\beta_{1}\varepsilon_{t-1}^{2}+\beta_{2}h_{t-1}
$$
+ the explanatory variable in the mean process ($ y_{t-1}^{2} $):
$$
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}
$$
All econometrical packages allowing to add 1)an arch in mean effect and 2) an explanatory variable can be used (the G@rch package (Ox) for example).