How to fit a GAM with double seasonality to a daily time series? (mgcv package)

Question

I'm trying to adjust a Generalized Additive Model to a daily time series. My goal is to do a short-term forecast for the gas demand of my city. I have data since 2015, including information about the weather (minimum and maximum temperatures). The data with which I'm working is daily data, and since I have information from many years back I have double seasonality: yearly and weekly.

Here is an image of the historical data throught the years. We can see that every year the demand rises during winter months (may, june, july, aug, sept):

And here is an image of the data weekly. We can see that during the end of the week the demand decreases:

The real question is the following: what is the right sintax for taking into account this two seasonalities when trying to adjust a gam() model from the mgcv package in Rstudio? I know that it should be something like this:

library(mgcv)
gam_1 <- gam(gas_demand ~ s(x1, bs = "cr", k = 7) +
               s(x2, bs = "ps", k = 365),
             data = df,
             family = gaussian)
#This is just an example, x1 and x2 have not been defined, they
#are supposed to be the covariates of the model
#The daily data should be stored in the data frame 'df'
#The response variable that I'm trying to forecast is 'gas_demand'

Should I create a column in my data frame that goes from 1 to 7, depending on the day of the week of the observation to take into account the weekly seasonality?

And for the yearly seasonality, should I create a column with values from 1 to 365 (depending on the day of the year)? Or a colum with values from 1 to 12, depending on the month of the year? I'm not sure what would be the right way to do it.

And my final question: which type of basis function is recommended for each type of seasonality?

I'm really desperate for a response since I'm struggling to find examples that work with this exact type of data. Thanks in advance! :)

Answer 1

As the comment from @user11852 suggests, you can incorporate multiple forms of seasonality (and how those seasonalities may change over time) using the cyclic cubic regression basis in the {mgcv} package (bs = 'cc'). You can see more information on the types of basis classes that are available in the help file for smooth.terms, and you can find particular information on the cyclic smooths in the help file on cubic regression splines.

Be aware however that forecasting from GAMs is a tricky business, as the splines will generally extrapolate linearly forever. Using the default 2nd derivative penalty from {mgcv}, these extrapolations are extremely sensitive to whatever wiggles the function happens to do at the boundary. There are better ways to handle nonlinear splines and generate good forecasts, which I have made available in the {mvgam} R package. This uses Stan for full Bayesian inference and incorporates latent dynamic processes in univariate or multivariate time series models. See this post for links to a brief webinar and summary of the issues of fitting GAMs to time series