Web Traffic Forecasting with Uncertainty

The pyINLA Model

Web traffic data typically consists of several overlapping patterns: a long-term trend (is interest growing or declining?), weekly cycles (do people browse more on weekends?), yearly cycles (are there seasonal topics?), and day-to-day momentum (if yesterday was busy, today probably will be too). Rather than trying to predict page views directly, we separate these patterns into individual components that we can model and interpret independently.

We model the log page views using this additive decomposition:

Here's what each symbol means:

• $y_t$ : the log page views on day $t$ (what we're trying to predict). We use the logarithm because page views are always positive and often have extreme spikes. Taking the log makes the data more symmetric and ensures that percentage changes are treated consistently: a doubling of traffic from 1,000 to 2,000 views is treated the same as a doubling from 10,000 to 20,000.
• $\mu$ : the overall average level of page views
• $f_{\text{trend}}(t)$ : a smooth function that captures long-term changes over time, where $t$ is the day number (1, 2, 3, ...). For example, overall interest in Peyton Manning grew during his career peak years and declined after retirement.
• $f_{\text{weekly}}(\text{dow}_t)$ : an adjustment based on which day of the week it is ($\text{dow}_t$ = day-of-week, where Monday=1 and Sunday=7)
• $f_{\text{yearly}}(\text{doy}_t)$ : an adjustment based on which day of the year it is ($\text{doy}_t$ = day-of-year, from 1 to 365). Unlike $t$ which increases forever, $\text{doy}_t$ resets to 1 every January 1st, so the same seasonal pattern repeats each year.
• $f_{\text{ar1}}(t)$ : captures day-to-day momentum (autocorrelation): if yesterday had unusually high traffic, today likely will too
• $\varepsilon_t$ : random noise that the model cannot explain

Here is how the indices $t$, $\text{dow}_t$, and $\text{doy}_t$ change over time:

Notice: $t$ always increases, $\text{dow}_t$ cycles 1-7, and $\text{doy}_t$ resets to 1 each January.

The "+" signs mean these effects are added together. So the predicted page views for any day is: average + trend adjustment + day-of-week adjustment + time-of-year adjustment + momentum + noise.

Model Components

Each component uses a random walk model (RW2), which assumes that values change smoothly from one point to the next. But there's an important distinction:

• Non-cyclic (for trend): Day 1 has no connection to day 2964. The trend can go up or down over time without any constraint to return to where it started.
• Cyclic (for weekly and yearly): The last value connects back to the first. For weekly effects, Sunday's value must smoothly connect to Monday's value because every week repeats. For yearly effects, December 31st connects to January 1st because every year repeats.

Calendar-Aligned Seasonality

A common mistake in time series modeling is to compute seasonal indices from a simple counter (1, 2, 3, ...) rather than from actual calendar dates. This can cause the model to misalign seasons: for example, if your data starts on a Wednesday, a naive approach might treat Wednesday as "day 1 of the week" instead of day 3.

We avoid this by extracting the day-of-week and day-of-year directly from the date. These values act as "lookup keys": for any given date, we use week_idx to find the corresponding weekly effect (one of 7 values) and year_idx to find the corresponding yearly effect (one of 365 values).

import numpy as np

# Calendar-aligned indices
dates = df['ds']  # datetime column from the dataframe
week_idx = dates.dt.dayofweek.values + 1   # Monday=1, Sunday=7
year_idx = dates.dt.dayofyear.values
year_idx = np.clip(year_idx, 1, 365)       # Map leap day 366 -> 365

What each line does:

• dates.dt.dayofweek extracts the day of week from each date (Python uses 0=Monday, 6=Sunday, so we add 1 to get 1=Monday, 7=Sunday)
• dates.dt.dayofyear extracts which day of the year it is (January 1 = 1, December 31 = 365)
• np.clip(..., 1, 365) handles leap years: February 29th would be day 366, but since it only occurs every 4 years, we map it to day 365 so our yearly pattern has a consistent length

Prior Specification

We use Penalized Complexity (PC) priors for all hyperparameters:

The AR1 prior is diffuse and centered at $\rho = 0$ (no autocorrelation), allowing the data to determine the degree of temporal dependence.

Interpretation of Yearly Pattern

The estimated yearly component shows a pattern with higher values during certain periods of the year. Attributing this pattern to specific domain events (e.g., NFL season) would require additional validation using event regressors or domain expertise; the pattern shown here is purely data-driven.

Component Decomposition

A key strength of pyINLA is the ability to examine learned components with full posterior uncertainty:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
random = result.summary_random

# (a) Trend
ax = axes[0, 0]
trend = random['t']
ax.plot(dates, trend['mean'].values, 'b-', linewidth=1.5)
ax.fill_between(dates, trend['0.025quant'].values, trend['0.975quant'].values,
                color='blue', alpha=0.2)
ax.axvline(dates.iloc[N_TRAIN], color='red', linestyle='--', alpha=0.5)
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_title('(a) Trend (RW2)')

# (b) Weekly seasonality
ax = axes[0, 1]
week = random['week']
days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
x = np.arange(7)
ax.bar(x, week['mean'].values[:7], color='green', alpha=0.7)
ax.errorbar(x, week['mean'].values[:7],
            yerr=[week['mean'].values[:7] - week['0.025quant'].values[:7],
                  week['0.975quant'].values[:7] - week['mean'].values[:7]],
            fmt='none', color='black', capsize=3)
ax.set_xticks(x)
ax.set_xticklabels(days)
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_title('(b) Weekly seasonality (cyclic RW2)')

# (c) Yearly seasonality
ax = axes[1, 0]
year = random['year_day']
ax.plot(np.arange(1, 366), year['mean'].values[:365], 'orange', linewidth=1)
ax.fill_between(np.arange(1, 366), year['0.025quant'].values[:365],
                year['0.975quant'].values[:365], color='orange', alpha=0.2)
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_title('(c) Yearly seasonality (cyclic RW2)')

# (d) AR1 residual component
ax = axes[1, 1]
ar = random['t_ar']
ax.plot(dates, ar['mean'].values, 'purple', linewidth=0.5, alpha=0.7)
ax.axvline(dates.iloc[N_TRAIN], color='red', linestyle='--', alpha=0.5)
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_title('(d) AR(1) residual component')

plt.tight_layout()

Interpreting the Components

• Trend: Shows the long-term evolution of page popularity, with uncertainty increasing in the forecast period
• Weekly: Day-of-week effects with error bars, revealing which days have higher or lower views
• Yearly: 365-day pattern capturing annual periodicity in the series
• AR1: The residual autocorrelation component with estimated $\hat{\rho} \approx 0.655$, meaning about 65% of yesterday's unexplained deviation carries over to today