How to add controls (regressors) to GARCH model in R?

How can I estimate a GARCH(1,1) model with control variables like this: $$Y_t=a_0+a_1X_t+e_t$$ where$$e_t\sim N(0,h_t)$$ $$h_t=b_0+b_1e_{t-1}^2+b_3h_{t-1}+b_3Z_t$$ I've checked some packages but can't fix it. Hope you guys can shed some lights!

• I'm voting to close this question as off-topic because it seems better suited for Cross Validated SE. – LocalVolatility Nov 23 '17 at 11:00
• @LocalVolatility, this is clearly off topic on Cross Validated because Cross Validated does not deal with software implementation of this kind. – Richard Hardy Nov 24 '17 at 19:19

library(rugarch)

The object fit contains the fitted model. If you uncomment the second line, you can try it out with randomly generated data.
(I suppose you meant $b_2$ -- not $b_3$ -- before $h_{t-1}$.)