Time series forecasting - Residuals not white noise This is my first message on CrossValidated to get some insights on an issue I am facing while trying to model properly a time series. I am relatively new to this science so please brace with me.
My time series is related to energy consumption with a daily seasonality (96 values per day). Here a description:


For the model definition, I used an auto.arima as a start but the model returned (SARIMA (1,0,0)(0,1,0)[96]) has residuals that are not white noise, this means that the stochastic part of my series is not entirely modeled. Below a summary of residuals checking:

Ljung-Box test

data:  Residuals from ARIMA(5,1,1)(0,1,1)[96]
Q* = 248.44, df = 185, p-value = 0.001282

Model df: 7.   Total lags used: 192 

I tried to create my own model using differentiation but still end up with residuals that are not white noise which is frustrating. I differentiated twice to remove seasonal and trend patterns, this is what I obtained :
After first differentiation:

After second differentiation:

I chose a SARIMA(0,1,4)(0,1,1)[96] because:

*

*we can see a very significant spike on lag 96 on acf with an exponential decay on pacf seasonal lags which suggests a seasonal MA1.


*And we observe another significant spike on lag 4 on acf for which I chose a non-seasonal MA4.
At the end, I obtained these residuals:

Do you think it is reasonable to choose such models in regard to the significance of residuals correlations ? Do you have any advice on how should I carry on from here ?
Thank you in advance for your help ?
 A: Well, 96 values per day means you have data every 15 minutes. I think the best first step would be to aggregate to hourly if you can, then you would have a pretty standard hourly seasonal shape.  Then you would still have multiple seasonal patterns to worry about but it would be much simpler IMO.  You would have 24 for the daily pattern then potentially 24 * 7 to account for a weekly pattern.
With that said, if you have patterns in your residuals that just means you are not accounting for those patterns in your model.  Now should we change our order to pick up this signal in our residuals?  No, I would fit to produce the best forecasts on validation sets, which is our ultimate goal anyway. If the residuals contain real signal then that would mean a different order could produce better forecasts, and if it is just noise then the new order would produce worse forecasts.  There are tons of statistical tests for white noise and all that but the best test for us is to just test the forecast.
If you are new to time series then I always recommend more plug and play models like fbprophet.  It also handles multiple seasonality quite well which seems like your biggest driver.
A: While the SARIMA model includes a seasonal component, it is a seasonal component where the period is stochastic rather than deterministic.  In the present case, where you have a fixed number of data points per day (occurring in 15 minute intervals) you would expect to have daily signals with a fixed period of 96 time units.  Your diagnostic plots confirm this, since they show a clear "spike" in autocorrelation at this lag.
In view of this, my first recommendation would be to create a Fourier intensity plot for your residuals to confirm that there is a fixed signal in the frequency domain, at the period $\Delta = 96$.  To incorporate this signal into your model, there are two choices: (1) you could add 95 model terms to act as effects for each of the time-periods in the day (if you a large enough volume of data to deal with the loss of degrees-of-freedom) or (2) you could add a relatively small number of sinusoidal terms at this frequency and its harmonics.  Adding a periodic term with a fixed period out to deal with the clear daily effect in your data, which will allow you to model the remaining aspects of your data more appropriately.  As a flow-on effect, this is likely to reduce the excess kurtosis in your residual plot.
