Understanding negative gamma value for the GJR-GARCH model:
$\gamma > 0$ is not a required condition to ensure a "valid" GJR-GARCH model. Let me explain why:
As you probably know, we need to impose some restrictions on the parameter space in order to obtain a proper volatility model. The two requirements we need to ensure, are positivity (positive estimates) and covariance stationarity. For simplicity, let us vaguely define the GJR-GARCH(1,1) model (I'm skipping the mean-model and thus imposing constant $\mu$):
\begin{align*}
r_t \vert \mathcal{F}_{t-1} &= \mu + \varepsilon_t\\
\varepsilon_t &= \sigma_t \cdot z_t\\
\sigma^2_t &= \omega + \alpha \varepsilon_{t-1}^2 + \beta \sigma_{t-1}^2 + \gamma I_{t-1} \varepsilon_{t-1}^2,
\end{align*}
where $z_t \overset{iid}{\sim} D(0,1)$ (which in your case is the skewed generalized distribution sged) and
$$I_{t-1} =\begin{cases}
1 & \text{if } \varepsilon_{t-1} < 0 \\
0 & \text{if } \varepsilon_{t-1} \geq 0
\end{cases}.$$
Here, positivity is still satisfied when we impose $\omega, \beta,\alpha > 0$ and $\alpha + \gamma > 0$. The latter condition is a broader statement than imposing $\alpha, \gamma >0$, since we can allow one of the parameters to become negative (in your case, $\gamma$). As seen from your parameter estimates, we have that $\alpha > \gamma$ and therefore positivity is still ensured. Again, it is trivial that positivity is still satisfied when $\gamma > 0$. Looking at the unconditional variance for the return process:
\begin{equation}
\mathbb{V}ar(r_t) := \sigma_t^2 = \frac{\omega}{1 - \alpha - \beta - \kappa \gamma},
\end{equation}
we can ensure covariance stationarity by restricting $0 < \alpha + \beta + \kappa \gamma < 1$ and $\omega>0$, where $\kappa = \mathbb{E}\left[I_{t-1} z_{t-1}^2\right] = \mathbb{P}(z_{t-1}<0)$ and is 0.5 for symmetric distributions.
Again, looking at your parameter estimates, we see that $0.2013 + 0.7828 - 0.0961 \cdot \kappa < 1$ is satisfied for $\kappa \in [0,1]$ (so, even though you are not working with a symmetric distribution, your model is still covariance stationary, in this particular scenario).
In general, when fitting the GJR-GARCH model on equities, you will often end up with a positive gamma parameter. When $\gamma > 0$ we observe asymmetrical effects in the volatility process, leading us to the conclusion that negative return-shocks causes larger variance. However, this does not imply that you'll get the same results for other asset classes. To provide some comfort, The V-lab at NYU have fitted a GJR-GARCH on the gold spot and likewise get a negative parameter for gamma.
Intuitive & technical reasons for a negative gamma parameter:
Here are my two cents on, why you are obtaining a negative gamma parameter:
Gold is a safe-haven asset and exhibit an opposite asymmetrical leverage effect, as opposed to equities. In times of crisis, many institutional investors reallocate large equity positions into gold and other "safe-haven" assets (in general, "safe-haven" assets are either uncorrelated (or slightly negatively correlated) with the equity market. In crisis they exhibit a stronger negative correlation, making them great for equity-portfolio hedges). This inevitably causes an opposite asymmetrical response (symmetrical response) in the gold prices, where future volatility is more affected by past positive returns (than negative). In essence, institutional investors collected moves from equity positions to gold positions, might result in an overall negative gamma parameter. This is further emphasised in the paper of Stavroyiannis (2018), where he also constructs a bivariate VAR model and finds short-run Granger causality from S&P 500 to the gold spot index (but not reverse), thus implying that past S&P 500 returns help us explain gold spot returns (empirical evidence of gold being a "safe-haven" asset for investors). This analysis was done from 2000 to 2016 and included several crises.
The model is trying to downscale persistence or volatility by letting $\gamma < 0$. A combination of high $\alpha$, $\beta$ together with a negative $\gamma$ might impose a better fit, than decreasing $\alpha$ and $\beta$ all together (and letting $\gamma \approx 0$).
The parameter estimations might be very dependent on the software you use. As described in Stavroyiannis (2018) you might get different parameter estimations depending on the software:
[...] the results of a GJR model depend highly on the software
used, affected by the inclusion or exclusion of certain constraint inequalities in the programming
approach of the optimization procedure. Eviews v.8.1 and OxMetrics v.7.1 allow for a negative 𝑎
parameter as far as 𝑎 + 𝛾 > 0, while both of them including Matlab v.2014a, R language v.3.3.2
using the rugarch package (Galanos, 2015), and Gretl v.2016d allow for a negative 𝛾 parameter
as far as 𝑎 + 𝛾 > 0, and 𝑎 > 0.
In conclusion, the term $\gamma I_{t-1} \varepsilon_{t-1}^2$, might not be adequate to explain the leverage-type effect for the gold return-process and changing the indicator function to allow for positive return variation, might improve the model.
Imposing specific bounds on $\gamma$ in the rugarch package:
Also, as @Kermittfrog wrote in the comments, if you want to impose a zero lower bound on $\gamma$ you can call setbounds() on the ugarchspec:
s_8 <- ugarchspec(variance.model = list(model="gjrGARCH",
garchOrder=c(1,1)),
mean.model = list(armaOrder=c(2,2)),
distribution.model = "sged")
setbounds(s_8)<-list(gamma1=c(0,1))
f_8<- ugarchfit(spec = s_8, data = ret)
f_8
But you will probably force $\gamma$ towards zero, by doing this. I hope this gives a little bit of insight and clarity on the GJR-GARCH model as-well as possible reasons for a negative $\gamma$ parameter.
