Modeling ARCH effects in Vector Autoregression

Hi all,

I have constructed a VAR(p) model with 3 time series in Stata. When I predict the residuals and run an AR(p) model (with p=1, 2, 3, and 4) on the residuals, I get significant lag coefficients (i.e. ARCH effects). I learned that I could now model the residuals using those AR(p) coefficients but I don't know how to implement this in Stata. The ultimate goal is go get model (VAR+ARCH) residuals that are free of ARCH effects.

My code so far:
var IV GDP UNEMP, lags(1/4)
predict resIV, res equation(IV)
gen redIVsqr=resIV^2
arima resIVsqr, arima(1,0,0)
arima resIVsqr, arima(2,0,0)
arima resIVsqr, arima(3,0,0)
arima resIVsqr, arima(4,0,0)