Using hac for arima models

Hi, sorry if this an obvious question, but I can't quite seem to figure it out.

I used 'estimate' in Matlab to estimate an ARIMA(2,0,2) model. From what I understand Matlab does this by MLE and then uses the inverse of the outer product of the score to get a variance-covariance matrix. I'm looking to get the full heteroskedastic-robust variance-covariance matrix though. In the MLE setting this would be the sandwich form of the Hessian and the outer product of the score (H^-1 ss' H^-1).

I noticed that Matlab has the function 'hac' that will give White or Newey-West standard errors using the full sandwich form of the variance-covariance matrix. This seems to only really work in a regression or GMM setting. Is there an analogue function for 'hac' that I could use for arima models?

Thanks!