Covariance of forecasts from Python/statsmodels SARIMAXResults object

Question

I have a SARIMAX model fitted at daily frequency using statsmodels.tsa.statepsace.sarimax. It is a "full" SARIMAX model in the sense that it has AR, MA, seasonal and exogenous components.

If I have correctly understood the statsmodels implementation of (linear Gaussian) state space models, then the forecasts of the next 31 days after the historical data form a multivariate Gaussian distribution. I would like know the parameters of that distribution.

I can find the mean vector easily enough with

fittedModel.forecast(31, exog=fwd_exog)

I can estimate the covariance matrix with

fittedModel.simulate(31, anchor='end', repetitions=10000, exog=fwd_exog).transpose().cov()

but that is a little slow and is only a numerical estimate.

It seems like maybe I could derive this from fittedModel.states.filtered_cov but I'm not clear on how to do that.

I have two questions:

1. Am I correct that the documented structure of state space models in statsmodels guarantees that the forecasts will be jointly Gaussian?
2. Does statsmodels.tsa.statespace.sarimax.SARIMAXResults have a method that will provide the covariances of the forecasts? If not, is there a straightforward way to compute that matrix from the methods that it does provide?