How to model a GARCH(1,1) with covariate? The purpose of my study is to understand if changes in environment policy or changes in people concerns about climate change affects volatility or if they can help in the prediction of volatility. In this case I check the volatility of the Dow Jones from January 2003 to June 2018. In order to do it, I’m using a GARCH(1,1) model with the addition of a covariate. The covariates I chose are two indeces, one called EnvP, which measures the salience of US environmental policy, and the other one called Media climate change concerns, which captures unexpected increases in climate change concerns. I use the software gretl for my analysis. 
Anyway, I got a lot of doubts about how to estimate this model.
For a start, I estimate a first GARCH(1,1) to give an initial interpretation of the parameters and to use it as a sort of benchmark:

The first doubt I have is: do I have to insert a lag of the returns in the conditional mean equation of the GARCH? 
In that case I’d have:

Is there any significant change? Is there a motive I should or shouldn’t insert this lag?
Now, I insert the covariate (EnvP) in the model:

Same model but with the returns lagged(1):

I need help understanding these outputs, to understand if it’s methodologically correct to insert those lags and, in general, if I can say something regarding the aim of my study. 
Any help will be very appreciated.
 A: 
[D]o I have to insert a lag of the returns in the conditional mean equation of the GARCH? <...> Is there any significant change? Is there a motive I should or shouldn’t insert this lag?

If Dow Jones were an AR(1)-GARCH(1,1) process, it would make sense to include the first lag of Dow Jones in the conditional mean equation. Otherwise the residuals would be autocorrelated and would violate the assumptions of the model. Of course, Dow Jones is not exactly an AR(1)-GARCH(1,1) process, but perhaps the conditional mean is decently approximated by that. To find out, examine the standardized residuals of the model and their ACF and PACF in particular. The standardized residuals should not be autocorrelated at any lag. (Note that an uncorrelated series will produce significant estimated autocorrelations in 5% of all lags under examination given a 95% confidence level. Thus, significant autocorrelation at a small proportion of lags does not necessarily mean trouble.)
Also, some of the information criteria (AIC, HQ) suggest that including the first lag improves the fit to a sufficient degree to warrant the inclusion, while BIC is fairly indifferent regarding the two models.

I need help understanding these outputs, to understand if it’s methodologically correct to insert those lags and, in general, if I can say something regarding the aim of my study.

In both const-GARCH(1,1) and AR(1)-GARCH(1,1), the estimated effect of the EnvPIndex is virtually the same. (That makes your life a little easier, since whichever model you choose, the variable of interest contributes to it in the same way.) Since it is statistically significant, it seems that when EnvPIndex is larger, the conditional variance is larger, and this is unlikely to be due to chance. Whether the "effect" is economically large or small depends on the scale of the index relative to the scale of $\sigma^2_t$. E.g. if the variance of EnvP is much larger than the variance of $\sigma_t^2$, the economic effect could be large, while if not it may be small.
