I have time series data recorded at multiple locations, stored in a matrix $Y$. I have fit a Vector Autoregressive Model to it which forecasts the data pretty well on a test set. However, if I plot the residuals there is definitely some time series structure left; ideally the residuals would look like noise.
I have been learning about boosted trees (package 'gbm' in R) that fits a response vector $Y$ to some input features $X$, pulls the residuals $E = Y - f(X)$, and then fits $E$ to the same input features $X$ again (repeat until the process starts to overfit). It's my belief that a tree won't be able to model time series structure well, so I am wondering if the same process can be done using an AR model as the weak learner.
After fitting my VAR model, I take the residuals and fit an AR process to it, and forecast it using the 'forecast' package in R. I then take my VAR forecasts and add them to my AR forecasts to get a final prediction. The error at some locations decrease, but for most locations the error increases.
Below are the plots of the residuals having only fit the VAR model - there is definitely some structure that can still be explained using some model. Any suggestions?
Update: The time series represent energy consumption at various locations. During the night we expect energy consumption to be very low and flat, during the day there may be high volatility so the variance is definitely not constant. My first step was to remove the daily seasonal component using STL (Seasonal Trend Decomposition Using Loess). After that, I fit my VAR model, and subsequently another AR model to the residuals. The time series after removing the Seasonal Component look like this:
The error terms from the VAR look like