I'm doing analysis with a model as below with autocorrelated residuals on time series data, and am trying to make predictions with pred intervals.
model<- gls(data=data,
outcome~variable1+variable2, correlation = corARMA(p = 1, q = 0),
method="ML")
I'm manually applying "arima" function on residuals of this model to make predictions with intervals, but am wondering if there is any R package to do this whole procedure with a single function.
Simply applying "predict" does not seem to reflect the AR(1) structure of residuals.