Modeling Airline Passenger Seasonality

Step 3: Interpret the Results

Year Effects (Trend)

The year coefficients $\alpha_{\text{year}}$ are estimated on the log scale. To convert back to the original passenger scale, we exponentiate:

For example, for 1949 we have $\alpha_{1949} = 4.855$, so:

# Convert log-scale year means to original scale
year_means = result.summary_fixed['mean']
print("Average monthly passengers by year (thousands):")
for year_label, log_mean in year_means.items():
    print(f"  {year_label.replace('year', '')}: {np.exp(log_mean):.0f}")

Average monthly passengers by year (thousands):
  1949: 128
  1950: 141
  1951: 172
  1952: 200
  1953: 228
  1954: 241
  1955: 287
  1956: 331
  1957: 371
  1958: 383
  1959: 431
  1960: 480

Passengers nearly quadrupled from 1949 to 1960 (3.8x growth), reflecting the rapid growth of commercial aviation.

Seasonal Effects

The seasonal random effects capture the monthly pattern within each year:

# Extract seasonal effects (first 12 values = monthly pattern)
seasonal = result.summary_random['month']
first_year = seasonal.iloc[:12][['mean', '0.025quant', '0.975quant']].copy()
first_year.columns = ['Effect', 'Lower', 'Upper']

month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
               'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
first_year['Month'] = month_names

print("Seasonal Effects (deviations from year mean):")
print(first_year[['Month', 'Effect', 'Lower', 'Upper']].to_string(index=False))

Seasonal Effects (deviations from year mean):
Month    Effect     Lower     Upper
  Jan -0.159375 -0.227468 -0.091198
  Feb -0.171319 -0.239346 -0.103270
  Mar -0.031091 -0.098980  0.036921
  Apr -0.052377 -0.120124  0.015620
  May -0.044802 -0.112390  0.023216
  Jun  0.087265  0.019847  0.155333
  Jul  0.201182  0.133886  0.269269
  Aug  0.201991  0.134712  0.270008
  Sep  0.067588  0.000224  0.135450
  Oct -0.060295 -0.127786  0.007387
  Nov -0.193776 -0.261374 -0.126237
  Dec -0.069798 -0.137420 -0.002301

Trend and Seasonality Together

Now we can visualize both components side by side:

# Plot year effects (trend) and seasonal pattern side by side
year_effects = result.summary_fixed
years_list = [int(y.replace('year', '')) for y in year_effects.index]  # strips 'year' prefix from e.g. 'year1949'
year_means = np.exp(year_effects['mean'].values)
year_lower = np.exp(year_effects['0.025quant'].values)
year_upper = np.exp(year_effects['0.975quant'].values)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Left: Year effects on original scale with 95% CI
axes[0].fill_between(years_list, year_lower, year_upper, alpha=0.3, color='#3498db')
axes[0].plot(years_list, year_means, 'o-', color='#2c3e50', lw=2, markersize=8)
axes[0].set_xlabel('Year')
axes[0].set_ylabel('Average Passengers (thousands)')
axes[0].set_title('Estimated Yearly Mean (with 95% CI)')
axes[0].grid(True, alpha=0.3, linestyle='--')
axes[0].set_xticks(years_list)
axes[0].set_xticklabels(years_list, rotation=45)
axes[0].annotate(f'{year_means[0]:.0f}k', xy=(years_list[0], year_means[0]),
                 xytext=(years_list[0]+0.5, year_means[0]+30), fontsize=9)
axes[0].annotate(f'{year_means[-1]:.0f}k', xy=(years_list[-1], year_means[-1]),
                 xytext=(years_list[-1]-1.5, year_means[-1]+30), fontsize=9, ha='right')

# Right: Seasonal pattern (% deviation from year mean)
seasonal_vals = first_year['Effect'].values
multipliers = np.exp(seasonal_vals)
pct_dev = (multipliers - 1) * 100
colors = ['#3498db' if m > 1 else '#e74c3c' for m in multipliers]
bars = axes[1].bar(month_names, pct_dev, color=colors, edgecolor='white', linewidth=1)
axes[1].axhline(0, color='k', lw=0.8)
axes[1].set_xlabel('Month')
axes[1].set_ylabel('% Deviation from Year Mean')
axes[1].set_title('Seasonal Pattern (% above/below yearly average)')
axes[1].grid(True, alpha=0.3, linestyle='--', axis='y')
axes[1].tick_params(axis='x', rotation=45)

for bar, m in zip(bars, multipliers):
    if abs(m - 1) > 0.15:
        height = (m - 1) * 100
        axes[1].annotate(f'{height:+.0f}%',
                         xy=(bar.get_x() + bar.get_width() / 2, height),
                         xytext=(0, 3 if height >= 0 else -12),
                         textcoords="offset points", ha='center',
                         va='bottom' if height >= 0 else 'top', fontsize=9, fontweight='bold')

plt.tight_layout()
plt.savefig('trend_and_seasonality.png', dpi=150, bbox_inches='tight')

# Bar chart of seasonal effects with value labels
seasonal_vals = first_year['Effect'].values
colors = ['#3498db' if e > 0 else '#e74c3c' for e in seasonal_vals]

fig, ax = plt.subplots(figsize=(10, 5))
bars = ax.bar(month_names, seasonal_vals, color=colors, edgecolor='white', linewidth=1.5)

for bar, val in zip(bars, seasonal_vals):
    height = bar.get_height()
    ax.annotate(f'{val:+.2f}',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3 if height >= 0 else -12),
                textcoords="offset points", ha='center',
                va='bottom' if height >= 0 else 'top', fontsize=9, fontweight='bold')

ax.axhline(0, color='k', lw=0.8)
ax.set_ylabel('Effect (log scale)')
ax.set_xlabel('Month')
ax.set_title('Seasonal Effects (deviations from year mean)')
ax.grid(True, alpha=0.3, linestyle='--', axis='y')
ax.set_ylim(-0.3, 0.35)

plt.tight_layout()
plt.savefig('seasonal_effects.png', dpi=150, bbox_inches='tight')

Interpretation:

• July-August peak: Highest travel (~22% above yearly average)
• November trough: Lowest travel (~18% below yearly average)
• Summer dominance: June-September all show positive effects
• Winter lull: January-February and November show the largest negative effects

Multiplicative Interpretation

Since the model is on the log scale, a month effect $u_m$ corresponds to a multiplicative factor $\exp(u_m)$ on the original scale:

The term $\exp(u_m)$ is a multiplier relative to the yearly baseline. To express as a percentage change:

Why subtract 1? When there is no seasonal effect ($u_m = 0$), the multiplier is $\exp(0) = 1$, meaning "100% of baseline" or "no change." Subtracting 1 converts the multiplier to a deviation: $1.22 - 1 = 0.22 = +22\%$ above baseline.

For example, if $u_{\text{Jul}} \approx 0.20$, then July is associated with an expected multiplier of $\exp(0.20) \approx 1.22$, i.e., about a 22% increase relative to the year baseline.

Comparing Two Months

Using the Seasonal model results, a direct comparison between two months $m_1$ and $m_2$ on the original scale is given by:

# Compare July (index 6) vs November (index 10) using seasonal model
u_jul = first_year.loc[6, 'Effect']  # 0.201
u_nov = first_year.loc[10, 'Effect']  # -0.194

ratio = np.exp(u_jul - u_nov)
print(f"July effect: {u_jul:+.3f}")
print(f"November effect: {u_nov:+.3f}")
print(f"Difference: {u_jul - u_nov:.3f}")
print(f"Ratio exp({u_jul - u_nov:.3f}) = {ratio:.2f}")
print(f"July demand is roughly {(ratio - 1) * 100:.0f}% higher than November")

July effect: +0.201
November effect: -0.194
Difference: 0.395
Ratio exp(0.395) = 1.48
July demand is roughly 48% higher than November

Note: This is a model-based point estimate from the Seasonal model. For rigorous inference, report the posterior credible interval for $\exp(u_{\text{Jul}} - u_{\text{Nov}})$.

# Convert seasonal effects to multiplicative factors
print("Monthly multipliers (relative to year average):")
for month, effect in zip(month_names, first_year['Effect']):
    multiplier = np.exp(effect)
    print(f"  {month}: {multiplier:.3f} ({(multiplier-1)*100:+.1f}%)")

Monthly multipliers (relative to year average):
  Jan: 0.853 (-14.7%)
  Feb: 0.843 (-15.7%)
  Mar: 0.969 (-3.1%)
  Apr: 0.949 (-5.1%)
  May: 0.956 (-4.4%)
  Jun: 1.091 (+9.1%)
  Jul: 1.223 (+22.3%)
  Aug: 1.224 (+22.4%)
  Sep: 1.070 (+7.0%)
  Oct: 0.941 (-5.9%)
  Nov: 0.824 (-17.6%)
  Dec: 0.933 (-6.7%)

Step 4: Model Fit Assessment

How well does the seasonal model fit the observed data? We compare fitted values $\hat{y}_t = \exp(\hat{\eta}_t)$ to observed passengers using descriptive metrics:

Coefficient of determination (R²) measures the proportion of variance explained by the model:

$R^2 = 1$ means perfect fit; $R^2 = 0$ means the model does no better than predicting the mean. Note: This is an in-sample measure and should not be interpreted as out-of-sample forecast accuracy.

Root mean squared error (RMSE) measures average prediction error in original units:

# Get fitted values on original scale (use .values to avoid index mismatch)
fitted_log = result.summary_fitted_values['mean'].values
fitted = np.exp(fitted_log)  # Back-transform from log scale
observed = df['passengers'].values

# Compute fit statistics
residuals = observed - fitted
ss_res = np.sum(residuals**2)
ss_tot = np.sum((observed - observed.mean())**2)
r_squared = 1 - ss_res / ss_tot
rmse = np.sqrt(np.mean(residuals**2))

print(f"R² = {r_squared:.4f}")
print(f"RMSE = {rmse:.1f} thousand passengers")

R² = 0.9925
RMSE = 10.4 thousand passengers

Visualizing the Fit

import matplotlib.pyplot as plt

# Create time index for x-axis
dates = pd.date_range('1949-01', periods=144, freq='M')

fig, axes = plt.subplots(2, 1, figsize=(10, 6), gridspec_kw={'height_ratios': [2, 1]})

# Top: Observed vs Fitted
axes[0].scatter(dates, observed, s=25, alpha=0.6, c='#2c3e50', label='Observed', zorder=3)
axes[0].plot(dates, fitted, color='#e74c3c', lw=2, label='Fitted', zorder=4)
axes[0].fill_between(dates, fitted, alpha=0.1, color='#e74c3c')
axes[0].set_ylabel('Passengers (thousands)')
axes[0].legend(loc='upper left', framealpha=0.95)
axes[0].set_title('Seasonal Model: Observed vs Fitted')
axes[0].grid(True, alpha=0.3, linestyle='--')

# Bottom: Residuals
colors = ['#27ae60' if r >= 0 else '#e74c3c' for r in residuals]
axes[1].bar(dates, residuals, width=25, color=colors, alpha=0.7, edgecolor='none')
axes[1].axhline(0, color='k', lw=0.8)
axes[1].set_ylabel('Residual')
axes[1].set_xlabel('Date')
axes[1].set_title('Residuals over Time')
axes[1].grid(True, alpha=0.3, linestyle='--')

plt.tight_layout()
plt.savefig('model_fit.png', dpi=150)

The seasonal model achieves R² ≈ 0.99 for this dataset, and fitted values closely track observations (residuals on the order of ±10 thousand passengers). These summaries describe in-sample agreement only; predictive performance should be assessed with out-of-sample validation or posterior predictive checks.