Hello, PyMC community!

Inspired by of Prophet-way seasonality modelling with airport example (here) and an interesting case by @juanitorduz here, I started modelling seasonality in PyMC with Fourier-series as well.

To create weekly, monthly and yearly seasonality, I created datasets with Fourier series and applied them in a model like this.

```
# === Creating Fourier features here. n_order is as recommended by Prophet for each type of seasonality
# Fourier features annual
n_order_annual = 10
periods_annual = df["date_day"].dt.dayofyear / 365.25
fourier_features_annual = pd.DataFrame(
{
f"{func}_order_{order}": getattr(np, func)(2 * np.pi * periods_annual * order)
for order in range(1, n_order_annual + 1)
for func in ("sin", "cos")
}
)
# Fourier features monthly
n_order_monthly = 5
periods_monthly = df["date_day"].dt.dayofyear / 30.5
fourier_features_monthly = pd.DataFrame(
{
f"{func}_order_{order}": getattr(np, func)(2 * np.pi * periods_monthly * order)
for order in range(1, n_order_monthly + 1)
for func in ("sin", "cos")
}
)
# Fourier features weekly
n_order_weekly = 3
periods_weekly = df["date_day"].dt.dayofyear / 7
fourier_features_weekly = pd.DataFrame(
{
f"{func}_order_{order}": getattr(np, func)(2 * np.pi * periods_weekly * order)
for order in range(1, n_order_weekly + 1)
for func in ("sin", "cos")
}
)
coords = {
"date":date,
"fourier_features_annual": np.arange(2 * n_order_annual),
"fourier_features_monthly": np.arange(2 * n_order_monthly),
"fourier_features_weekly": np.arange(2 * n_order_weekly)
}
with pm.Model(coords=coords) as model:
sales = pm.MutableData(name="sales", value=df['sales'])
intercept = pm.Normal(name="intercept", mu=0, sigma=200)
beta_trend = pm.Normal(name="beta_trend", mu=150, sigma=200)
trend = pm.MutableData(name="trend", value=df['trend'])
beta_fourier_annual = pm.Normal("beta_fourier_annual", mu=0, sigma=10, dims="fourier_features_annual")
seasonality_annual = pm.Deterministic("seasonality_annual", pm.math.dot(beta_fourier_annual, fourier_features_annual.to_numpy().T))
beta_fourier_monthly = pm.Normal("beta_fourier_monthly", mu=0, sigma=10, dims="fourier_features_monthly")
seasonality_monthly = pm.Deterministic("seasonality_monthly", pm.math.dot(beta_fourier_monthly, fourier_features_monthly.to_numpy().T))
beta_fourier_weekly = pm.Normal("beta_fourier_weekly", mu=0, sigma=10, dims="fourier_features_weekly")
seasonality_weekly = pm.Deterministic("seasonality_weekly", pm.math.dot(beta_fourier_weekly, fourier_features_weekly.to_numpy().T))
mu = pm.Deterministic(name="mu", var=intercept + trend*beta_trend +
seasonality_weekly + seasonality_monthly + seasonality_annual)
sigma = pm.HalfNormal(name="sigma", sigma=200)
pm.Normal(name="likelihood", mu=mu, sigma=sigma, observed=sales)
```

And this seems to work!

However, in my dataset I quickly found out that there are **interaction effects between types of seasonality**.

Specifically, weekends have an even larger uplift when they coincide with monthly uplift. (So weekly seasonality interacts with monthly).

Do you think thereâ€™s optimal way to model it?

One approach to modelling interaction effects is to multiply the two variables with a separate beta and include them to the model, like:

```
beta_interaction*variable1*variable2
```

But I am not sure if itâ€™s possible with Prophet-style modelling given both variables are matrices with different orders.

Would appreciate any recommendation!