Modelling Prophet-way seasonality, but with interaction effects

1

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!

2

I think an interaction effect will work fine, basically as you wrote it. The prophet model is an additive linear model, so you can think of the seasonal components as “features”. All the usual tricks apply.

2 Likes