Short-term gas demand forecasting I'm a 22 year old Statistics student with a big problem to solve. English isn't my first language, so I apologize in advance for any mistakes in my grammar.
I'm trying to make a short-term gas consumption forecast for a small city. In my country, gas consumption has its peak in the winter months (may, june, july, august and september), and in the rest of the year the consumption is much lower.
I'm working with daily data from 2015 up until 2022, but I can, if necessary, disregard the first few years.
Since I have daily data (one observation of the gas consumption per day), I have weekly and yearly seasonality. I also need to take into account some explanatory variables, the most important one being the daily mean temperature of the city.
Here I include a picture of the time series across the years.

To sum it up, I need a forecasting model that will handle daily data
with two types of seasonalities (weekly and yearly) and that also allows
me to incorporate covariates (the temperature and also some dummies for holidays and weekends).
In this context, with what type of model should I work with?
After applying a Box-Cox transformation to stabilize the variance of my variable (gas consumption), I tried using dynamic regression, following Rob Hyndman's book. I incorporated the following variables as covariates in the model: a dummy for saturdays, another for sundays and holidays, the mean daily temperature (in its linear and quadratic form) and the heating degree days (in its linear and quadratic form). The residuals plots look like this:

I also tried a dynamic harmonic regression with the following covariates: fourier terms, the mean daily temperature (in its linear and quadratic form) and the heating degree days (in its linear and quadratic form). The code I used was:
fit7 <- time_series %>%
  model(
    ARIMA( gasconsumption_transformed ~ PDQ(0, 0, 0) + pdq(d = 0) +
            mean_temp + I(mean_temp^2) + hdd + I(hdd^2)+
            fourier(period = "week", K = 3)+
            fourier(period = "year", K = 5))
  )
report(fit7)
fit7 %>% gg_tsresiduals()

And the residual plots were oddly similar to the previous fit with dynamic regression.

The residuals of the two previous models don't look particularly good. So I'm hoping to find another possible solution to my problem.
I've also read Rob J. Hyndman's paper titled "Short-term load forecasting based on a semi-parametric additive model", but I haven't found any package that may help me do this in the software R.
Another type of model that I've come across is state space models however in this case I haven't found any R package that implements this method either.
Any type of help with models that are easily implemented in R would be greatly appreciated! Also, if there's something hard to understand or information missing, I'll gladly re-check the cross-validated question.
One final comment: I am not allowed to share the data used.
 A: I've solved similar forecasting problems using Prophet. There are R and Python implementations, with many replicable examples on that webpage. It handles multiple seasonalities, holidays, and covariates.
A: First off, you have already implemented a lot of precisely those models one would recommend in this situation, congratulations!
One possibility for daily data would indeed be Prophet, as dimitry suggests. Another would be the models proposed by Hyndman et al. for multiple and complex seasonalities, like forecast::tbats() or forecast::mstl(), see the tag wiki. However, particularly the TBATS model is already very close to your own model with Fourier terms, so the results may not be very much better.
When I look at your residual plots, they indeed look rather similar. But they are not identical, and the similarity is probably a consequence of your models giving very similar fits, because I would suspect that both models extract about the same signal from the series. If that is so, then the other models suggested in the previous paragraph may not give very different fits (and residuals, and forecasts), either. One possible way forward would be to dig into what happened around the large residuals, e.g., the ones that seem to occur quite regularly around the end of the year. If you can understand these and model and forecast them, that might improve matters. However, since these seem to pertain to single days each, this may not be a systematic effect and just be a symptom of higher variance during "special times" - if so, there is probably little you can do in terms of an expectation forecast.
Indeed, there is probably a lot of residual variation in your data. At some point, forecasting algorithm will not be able to improve any further, simply because of this residual variation. See How to know that your machine learning problem is hopeless?
