I am trying to build a multivariate linear regression and the main goal is to understand how the various features impact the response by understanding the coefficients and their confidence intervals.
For this reason, I chose multivariate linear regression because the coefficients are intuitive to interpret and from the standard error and the degrees of freedom, I can get the 95% confidence interval of the coefficients. Therefore, I can tell what the impact of a unit increase of a predictor on the outcome is. I could use more complex models such as tree based models but, even if I can get the variable importance, it is not easy to quantify the coefficient of each variable.
My problem, though, is the data is time series data and hence it shows auto-correlation. I know using regression with ARIMA errors can handle the autocorrelation issue but I found it difficult to interpret the coefficients specially when d is non-zero in ARIMA(p, d, q). The attached image, for example, is model diagnosis for one of the models I have (I have to build 1000s of them).
How can I handle the autocorrelation issue and still get coefficients that are readily interpretable like in a multivariate linear regression? My residuals are not normally distributed and I am planning to use box-cox transform to see if it can solve that problem. But, I am not sure what I should use for the autocorrelation issue.