Forecasting Seasonal Trends in pandas: Decomposing Time Series into Trend, Seasonality, and Residuals for Business Reporting

Every business has rhythms. Retail surges in November and December. HVAC service calls spike in July and January. E-commerce subscriptions dip in summer. The problem isn't knowing these patterns exist — every analyst who's spent five minutes with a sales chart knows something seasonal is happening. The problem is quantifying those patterns precisely enough to make useful forecasts, set realistic targets, and tell the difference between genuine business growth and the same seasonal bump you got last year.

Standard reporting tools handle this poorly. You can build a 12-month rolling average in Excel, but decomposing a series into its mathematical components — isolating the underlying trend from the seasonal cycle and from random noise — requires statistical machinery that Excel simply doesn't offer. pandas, combined with statsmodels, gives you that machinery, and once you understand how decomposition works at a conceptual level, you can apply it to almost any business reporting problem where time is a variable.

By the end of this lesson, you'll be able to load and prepare real time series data in pandas, apply classical and STL decomposition to separate trend, seasonality, and residuals, interpret each component for business audiences, use those components to build simple seasonal forecasts, and export the results in a format that goes straight into a management report. You'll also understand when decomposition breaks down and what to do about it.

What you'll learn:

• How additive and multiplicative decomposition models work and which to choose
• How to implement seasonal_decompose and STL decomposition using statsmodels
• How to extract and interpret the trend, seasonal, and residual components for business reporting
• How to build a naive seasonal forecast using decomposition outputs
• How to diagnose decomposition quality using residual analysis
• How to integrate decomposed data into formatted Excel reports for stakeholders

---

Prerequisites

You should be comfortable working with pandas DataFrames, including loading data, indexing by date, and performing basic aggregations. If you need a refresher on date handling specifically, the lesson on Working with Dates and Time Series in pandas: Parsing, Resampling, and Rolling Windows covers the foundational mechanics. You should also have statsmodels installed (pip install statsmodels) and a working Python environment — Setting Up Python for Data Analysis: Python, VS Code, Jupyter, and Virtual Environments walks through that setup if you're starting fresh.

---

What Is Time Series Decomposition, and Why Does It Matter?

Before touching any code, you need a solid mental model of what decomposition actually does — because the math is straightforward, but the business interpretation is where most analysts go wrong.

A time series is just a sequence of values measured at regular intervals. Monthly revenue, weekly website visits, daily order counts. Decomposition asserts that any such series can be broken into three additive (or multiplicative) parts:

Trend (T): The long-run direction of the series, stripped of seasonal fluctuation and noise. If your company is growing, the trend component rises over time regardless of which month it is.

Seasonal (S): The repeating, calendar-driven pattern. December is always higher than August for most retailers. This component captures that structure. In a well-decomposed series, the seasonal factors repeat identically (or near-identically) each year.

Residual (R): Everything left over after trend and seasonality are accounted for. Ideally this is small, random, and structurally uninteresting. If your residuals have obvious patterns, your model is missing something.

The two decomposition models — additive and multiplicative — differ in how they combine these components:

• Additive: Y = T + S + R — use this when the seasonal swings are roughly constant in absolute terms regardless of the trend level. If December is always about $500K above the trend, that's additive.
• Multiplicative: Y = T × S × R — use this when seasonal swings scale with the level of the series. If December is always about 40% above the trend, and the trend is growing, December's absolute premium grows too. This is far more common in real business data.

A practical way to check: take the log of your series and apply additive decomposition. log(T × S × R) = log(T) + log(S) + log(R), which means multiplicative decomposition on raw data is mathematically equivalent to additive decomposition on log-transformed data. This trick lets you use additive machinery everywhere.

---

Setting Up the Dataset

We'll use a realistic scenario throughout this lesson: monthly sales figures for a mid-sized consumer electronics retailer, spanning January 2018 through December 2023. Six years of monthly data gives us enough seasonal cycles to decompose reliably (you need at least two full cycles, and three or more is better).

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from statsmodels.tsa.seasonal import seasonal_decompose, STL
from statsmodels.graphics.tsaplots import plot_acf
import warnings
warnings.filterwarnings('ignore')

# Generate realistic monthly revenue data
# In practice, load this from CSV, Excel, or a database
np.random.seed(42)

date_range = pd.date_range(start='2018-01-01', end='2023-12-01', freq='MS')
n = len(date_range)

# Build components manually so we know what "truth" looks like
trend_true = np.linspace(800_000, 1_400_000, n)  # Growing trend

# Seasonal pattern: strong Q4 (holiday), weak Q1 (post-holiday)
seasonal_pattern = np.array([
    -0.18, -0.12, -0.05, 0.02, 0.05, 0.08,
     0.04,  0.06,  0.07, 0.10, 0.22, 0.35
] * (n // 12 + 1))[:n]

# Multiplicative: seasonal impact scales with trend level
noise = np.random.normal(0, 0.03, n)
revenue = trend_true * (1 + seasonal_pattern + noise)

df = pd.DataFrame({
    'date': date_range,
    'revenue': revenue.round(2)
})
df = df.set_index('date')

print(df.head(12))
print(f"\nShape: {df.shape}")
print(f"\nRevenue stats:\n{df['revenue'].describe()}")

This produces a DataFrame with a proper DatetimeIndex, which is essential — statsmodels decomposition functions work directly with DatetimeIndex-based Series, not with date columns.

For loading real data, you'd typically do something like this:

# Loading from CSV (realistic version)
df = pd.read_csv('monthly_sales.csv', parse_dates=['date'])
df = df.set_index('date').sort_index()
df = df.asfreq('MS')  # Ensure monthly frequency is declared

# Check for gaps — decomposition can't handle missing periods
print(f"Expected periods: {len(pd.date_range(df.index.min(), df.index.max(), freq='MS'))}")
print(f"Actual periods: {len(df)}")
print(f"Missing months: {df['revenue'].isna().sum()}")

If you have missing months, you must fill them before decomposing. df['revenue'].interpolate(method='time') handles sparse gaps well; for larger gaps, you need domain judgment about whether interpolation makes sense at all. The lesson on Cleaning Messy Data with pandas: Missing Values, Duplicates, and Data Types covers the full toolkit for that step.

---

Classical Decomposition with `seasonal_decompose`

statsmodels.tsa.seasonal.seasonal_decompose implements classical decomposition, which works by:

1. Computing the trend using a centered moving average
2. Removing the trend from the original series
3. Averaging the de-trended values across the same month/quarter/week to estimate seasonal factors
4. Computing residuals as whatever is left over

It's fast, interpretable, and well-understood. It also has real limitations that we'll discuss. Let's run it first.

# Apply multiplicative decomposition (appropriate for growing revenue series)
result_mult = seasonal_decompose(
    df['revenue'],
    model='multiplicative',
    period=12,          # 12 months = 1 seasonal cycle
    extrapolate_trend='freq'  # Fill NaN at edges of trend component
)

# Inspect the components
components = pd.DataFrame({
    'observed': result_mult.observed,
    'trend': result_mult.trend,
    'seasonal': result_mult.seasonal,
    'residual': result_mult.resid
})

print(components.head(24).round(4))

The output gives you four columns. Let's understand what each means for a business audience:

# Trend: the underlying business trajectory
print("Trend range (first year):")
print(components['trend']['2018'].describe())

print("\nTrend range (last year):")
print(components['trend']['2023'].describe())

# Seasonal: multipliers that show calendar-driven variation
print("\nSeasonal factors by month (averaged across years):")
seasonal_by_month = components['seasonal'].groupby(components.index.month).mean()
seasonal_by_month.index = ['Jan','Feb','Mar','Apr','May','Jun',
                           'Jul','Aug','Sep','Oct','Nov','Dec']
print(seasonal_by_month.round(4))

The seasonal factors are multipliers in multiplicative decomposition. A December factor of 1.35 means December revenue is typically 35% above the trend line. A January factor of 0.82 means January runs about 18% below trend. These are the numbers your finance team actually wants when setting monthly targets.

# Convert seasonal factors to percentage deviation from trend
seasonal_pct = (seasonal_by_month - 1) * 100
print("\nSeasonal deviation from trend (%):")
print(seasonal_pct.round(1))

Now visualize all four components in a single chart:

fig, axes = plt.subplots(4, 1, figsize=(14, 12), sharex=True)
fig.suptitle('Monthly Revenue Decomposition (Multiplicative Model)', 
             fontsize=14, fontweight='bold', y=1.01)

# Format the x-axis for all subplots
for ax in axes:
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

components['observed'].plot(ax=axes[0], color='steelblue', linewidth=1.5)
axes[0].set_title('Observed Revenue', fontsize=11)
axes[0].yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'${x/1e6:.1f}M'))

components['trend'].plot(ax=axes[1], color='darkorange', linewidth=2)
axes[1].set_title('Trend Component', fontsize=11)
axes[1].yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'${x/1e6:.1f}M'))

components['seasonal'].plot(ax=axes[2], color='green', linewidth=1.5)
axes[2].set_title('Seasonal Component (multiplier)', fontsize=11)
axes[2].axhline(y=1, color='black', linestyle='--', linewidth=0.8)

components['residual'].plot(ax=axes[3], color='red', linewidth=1, alpha=0.7)
axes[3].set_title('Residual Component', fontsize=11)
axes[3].axhline(y=1, color='black', linestyle='--', linewidth=0.8)

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

For more on producing polished charts for reports, the lesson on Visualizing Data with matplotlib and seaborn: Charts That Explain Your Analysis has you covered on formatting and annotation.

---

Diagnosing Decomposition Quality: Residual Analysis

The residuals tell you how well your model fit. In a good decomposition, residuals should look like white noise — no obvious patterns, no autocorrelation, and a distribution centered near 1.0 (for multiplicative) or 0.0 (for additive).

residuals = components['residual'].dropna()

fig, axes = plt.subplots(1, 3, figsize=(16, 4))

# Distribution of residuals
axes[0].hist(residuals, bins=20, color='steelblue', edgecolor='white')
axes[0].axvline(x=1.0, color='red', linestyle='--', label='Expected mean (1.0)')
axes[0].set_title('Residual Distribution')
axes[0].set_xlabel('Residual Value')
axes[0].legend()

# Residuals over time
axes[1].plot(residuals.index, residuals.values, color='gray', linewidth=1)
axes[1].axhline(y=1, color='red', linestyle='--')
axes[1].set_title('Residuals Over Time')
axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

# Autocorrelation plot
plot_acf(residuals, ax=axes[2], lags=24, alpha=0.05)
axes[2].set_title('Residual Autocorrelation (ACF)')

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

# Statistical summary
print(f"Residual mean: {residuals.mean():.4f} (should be ~1.0 for multiplicative)")
print(f"Residual std:  {residuals.std():.4f} (smaller is better)")
print(f"Residual min:  {residuals.min():.4f}")
print(f"Residual max:  {residuals.max():.4f}")

# Check for outlier months
outlier_threshold = 3 * residuals.std()
distance_from_mean = (residuals - residuals.mean()).abs()
outliers = residuals[distance_from_mean > outlier_threshold]
print(f"\nOutlier months (>3σ from mean):")
print(outliers)

The autocorrelation plot is particularly important. If you see significant spikes at lag 12, 24, etc., the seasonal component wasn't fully captured. Spikes at lag 1 or 2 suggest the trend extraction is imperfect and there's remaining structure in the residuals that the model left on the table.

---

STL Decomposition: When Classical Isn't Enough

Classical seasonal_decompose has three well-known weaknesses:

1. The seasonal component is forced to repeat identically each year — it cannot evolve
2. It's sensitive to outliers because the moving average is not robust
3. It handles only integer-divisible periods cleanly

STL (Seasonal and Trend decomposition using Loess) fixes all three. It uses locally weighted regression (LOESS) to estimate both trend and seasonal components iteratively, with optional robustness weighting that downweights outliers. This makes STL significantly better for real business data, which often has structural shifts and exceptional events (like, say, a global pandemic in Q1 2020).

from statsmodels.tsa.seasonal import STL

# STL requires additive framing — use log transform for multiplicative data
df['log_revenue'] = np.log(df['revenue'])

stl = STL(
    df['log_revenue'],
    period=12,
    seasonal=13,         # Seasonal smoother window (odd number >= period + 2)
    trend=None,          # Let STL choose automatically (usually 2*period+1 rounded to odd)
    robust=True,         # Downweight outliers — strongly recommended for real data
    seasonal_deg=1,      # Degree of LOESS polynomial for seasonal (0=constant, 1=linear)
    trend_deg=1          # Degree of LOESS polynomial for trend
)

stl_result = stl.fit()

# Back-transform log components to original scale
stl_components = pd.DataFrame({
    'observed': np.exp(stl_result.observed),
    'trend': np.exp(stl_result.trend),
    'seasonal': np.exp(stl_result.seasonal),   # These are now multiplicative factors
    'residual': np.exp(stl_result.resid)        # Same
}, index=df.index)

print(stl_components.describe().round(2))

The seasonal parameter controls the width of the LOESS smoother applied to the seasonal component. Larger values make the seasonal component smoother and more stable across years; smaller values allow it to evolve more rapidly. seasonal=13 is a good default for monthly data — it's the smallest odd number that's at least period + 1. If you know your seasonal pattern is genuinely evolving (e.g., the December bump has grown year-over-year as your business matures), you might experiment with values as high as seasonal=25 or seasonal=35.

Let's compare classical and STL decomposition on the same data:

# Side-by-side trend comparison
fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

# Trend comparison
axes[0].plot(components['trend'], color='darkorange', linewidth=2, 
             label='Classical (MA-based)', linestyle='--')
axes[0].plot(stl_components['trend'], color='steelblue', linewidth=2, 
             label='STL (LOESS-based)')
axes[0].set_title('Trend Component: Classical vs. STL', fontsize=12)
axes[0].yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'${x/1e6:.1f}M'))
axes[0].legend()
axes[0].xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

# Residual comparison
axes[1].plot(components['residual'], color='darkorange', linewidth=1, 
             alpha=0.7, label='Classical residual')
axes[1].plot(stl_components['residual'], color='steelblue', linewidth=1, 
             alpha=0.7, label='STL residual')
axes[1].axhline(y=1, color='black', linestyle='--', linewidth=0.8)
axes[1].set_title('Residuals: Classical vs. STL', fontsize=12)
axes[1].legend()

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

---

Extracting Seasonal Factors for Business Reporting

The seasonal component is the most immediately useful output for business reporting. It tells you exactly how much of each month's performance is simply calendar-driven versus actual business dynamics.

# Extract clean seasonal factors from STL result
# Use the last full year of seasonal factors (most current estimates)
last_year_seasonal = stl_components['seasonal']['2023']

# Seasonal index: how does each month compare to the annual average?
seasonal_index = pd.DataFrame({
    'month': ['Jan','Feb','Mar','Apr','May','Jun',
              'Jul','Aug','Sep','Oct','Nov','Dec'],
    'seasonal_factor': last_year_seasonal.values,
    'pct_above_below_trend': ((last_year_seasonal.values - 1) * 100).round(1)
})

# Add ranking
seasonal_index['rank'] = seasonal_index['seasonal_factor'].rank(ascending=False).astype(int)
seasonal_index = seasonal_index.sort_values('rank')

print("Seasonal Index by Month:")
print(seasonal_index.to_string(index=False))

For business audiences, expressing this as "December is typically 32% above trend while January runs 19% below trend" is immediately actionable. Finance teams can use these factors to set realistic monthly targets on a growing trend line.

# Build monthly target model for 2024
# Step 1: Project the trend forward using linear extrapolation from STL trend

from scipy.stats import linregress

trend_series = stl_components['trend'].dropna()
# Use last 24 months for the projection basis
recent_trend = trend_series[-24:]
x = np.arange(len(recent_trend))
slope, intercept, r_value, p_value, std_err = linregress(x, recent_trend)

print(f"Trend growth rate: ${slope:,.0f} per month")
print(f"Trend R²: {r_value**2:.4f}")
print(f"Trend p-value: {p_value:.6f}")

# Step 2: Project 12 months forward
future_x = np.arange(len(recent_trend), len(recent_trend) + 12)
projected_trend = intercept + slope * future_x

# Step 3: Apply seasonal factors to projected trend
forecast_dates = pd.date_range('2024-01-01', periods=12, freq='MS')
forecast_df = pd.DataFrame({
    'date': forecast_dates,
    'projected_trend': projected_trend,
    'seasonal_factor': seasonal_index.sort_values('month')['seasonal_factor'].values
    # Note: this assumes seasonal_index is in Jan-Dec order; adjust if not
})

# Reorder seasonal_index properly
monthly_seasonal = stl_components['seasonal'].groupby(
    stl_components.index.month
).mean()

forecast_df['seasonal_factor'] = [monthly_seasonal[m] for m in forecast_df['date'].dt.month]
forecast_df['forecast_revenue'] = forecast_df['projected_trend'] * forecast_df['seasonal_factor']
forecast_df['date'] = forecast_df['date'].dt.strftime('%b %Y')

print("\n2024 Revenue Forecast:")
print(forecast_df[['date','projected_trend','seasonal_factor','forecast_revenue']].
      to_string(index=False, float_format='${:,.0f}'.format))

---

Handling Structural Breaks: What Happens When the Pattern Shifts

Real business data has structural breaks — points where the underlying data-generating process changed permanently. COVID-19 in Q1-Q2 2020 is the obvious example, but structural breaks also happen due to major product launches, pricing changes, market entries, or acquisitions.

Classical decomposition handles these badly: the unusual observations get averaged into the seasonal factors, distorting them. STL with robust=True handles them better because outlier observations get down-weighted.

But the best approach is explicit: identify the break, decompose the pre-break and post-break periods separately, and only use post-break factors for future projections.

# Detect anomalous periods using residual magnitude
residuals_stl = stl_components['residual']
residual_zscore = (residuals_stl - residuals_stl.mean()) / residuals_stl.std()

anomalies = residual_zscore[residual_zscore.abs() > 2.0]
print("Anomalous months (|z-score| > 2.0):")
print(pd.DataFrame({
    'month': anomalies.index.strftime('%b %Y'),
    'residual': anomalies.values.round(4),
    'z_score': residual_zscore[anomalies.index].round(2)
}))

# If you know there's a structural break (e.g., April 2020), 
# decompose post-break separately

post_break = df['log_revenue']['2020-07-01':]  # 6-month grace period after the shock

if len(post_break) >= 24:  # Need at least 2 full cycles
    stl_post = STL(post_break, period=12, seasonal=13, robust=True).fit()
    
    # Post-break seasonal factors (more relevant for future forecasting)
    post_break_seasonal = pd.Series(
        stl_post.seasonal, 
        index=post_break.index
    )
    
    # Use only these for future projections
    current_seasonal = np.exp(
        post_break_seasonal.groupby(post_break_seasonal.index.month).mean()
    )
    print("Post-break seasonal factors:")
    print(current_seasonal.round(4))
else:
    print(f"Only {len(post_break)} periods post-break — not enough for reliable decomposition")
    print("Using full-series factors with robust=True as fallback")

---

Seasonally Adjusting a Series for Apples-to-Apples Comparisons

One of the most powerful applications of decomposition in business reporting is seasonal adjustment — removing the seasonal component to see the "pure" underlying performance. This lets you compare months that would otherwise be incomparable.

# Seasonally adjusted revenue: removes the seasonal multiplier
# Observed = Trend × Seasonal × Residual
# Seasonally adjusted = Trend × Residual = Observed / Seasonal

df['seasonally_adjusted'] = df['revenue'] / stl_components['seasonal']
df['trend_only'] = stl_components['trend']

# Now compare Q4 vs Q1 without seasonal distortion
q4_2023 = df.loc['2023-10':'2023-12', 'seasonally_adjusted'].mean()
q1_2024_actual = df.loc['2024-01':'2024-03', 'revenue'].mean() if '2024' in str(df.index.max()) else None

print(f"Q4 2023 seasonally adjusted avg monthly revenue: ${q4_2023:,.0f}")

# MoM growth on seasonally adjusted series
df['sa_mom_growth'] = df['seasonally_adjusted'].pct_change() * 100

print("\nSeasonally Adjusted Month-over-Month Growth:")
print(df['sa_mom_growth'].dropna().tail(12).round(2).to_string())

The pct_change() on seasonally adjusted figures gives you the "real" momentum of the business, independent of calendar effects. A December that shows +2% seasonally adjusted growth is genuinely stronger than a December showing +25% raw growth if the seasonal factor is normally +30%.

For more on period-over-period comparisons in pandas, the lesson on Calculating Month-over-Month and Year-over-Year Changes in pandas: pct_change, shift, and Period Comparisons for Business Reporting covers the full toolkit with shift() and pct_change().

---

Scaling to Multiple Product Lines or Regions

In practice, you're rarely decomposing a single aggregate series. You need to run decomposition across multiple dimensions — product categories, regions, sales channels — and then aggregate the results coherently.

# Simulate multi-category data
categories = ['Consumer Electronics', 'Accessories', 'Services']
category_data = []

for cat in categories:
    # Different trend slopes and seasonal patterns per category
    scale = {'Consumer Electronics': 1.0, 'Accessories': 0.3, 'Services': 0.15}[cat]
    seasonal_shift = {'Consumer Electronics': 0, 'Accessories': 1, 'Services': -2}[cat]
    
    cat_revenue = df['revenue'].values * scale
    # Add category-specific noise
    cat_revenue *= (1 + np.random.normal(0, 0.04, len(df)))
    
    category_data.append(pd.DataFrame({
        'date': df.index,
        'category': cat,
        'revenue': cat_revenue
    }))

multi_df = pd.concat(category_data).set_index(['date', 'category'])

# Decompose each category separately
def decompose_category(series, period=12):
    """Run STL decomposition on a single time series, return components DataFrame."""
    if series.isna().any():
        series = series.interpolate(method='time')
    
    log_series = np.log(series)
    stl_result = STL(log_series, period=period, seasonal=13, robust=True).fit()
    
    return pd.DataFrame({
        'trend': np.exp(stl_result.trend),
        'seasonal': np.exp(stl_result.seasonal),
        'residual': np.exp(stl_result.resid),
        'seasonally_adjusted': series / np.exp(stl_result.seasonal)
    }, index=series.index)

# Apply across all categories
all_components = {}
for cat in categories:
    cat_series = multi_df.xs(cat, level='category')['revenue']
    all_components[cat] = decompose_category(cat_series)
    print(f"Decomposed {cat}: {len(cat_series)} periods")

# Combine into a single multi-level DataFrame
combined = pd.concat(all_components, axis=1)
combined.columns = pd.MultiIndex.from_tuples(
    [(cat, comp) for cat in categories for comp in ['trend','seasonal','residual','seasonally_adjusted']],
    names=['category','component']
)

print("\nMulti-category decomposition complete")
print(f"Shape: {combined.shape}")
print(f"Categories: {categories}")

# Robust multi-series decomposition with error logging
def batch_decompose(grouped_series_dict, period=12, min_periods=24):
    """
    Decompose multiple time series, logging failures.
    Returns dict of successful decompositions and list of failures.
    """
    results = {}
    failures = []
    
    for name, series in grouped_series_dict.items():
        try:
            series = series.dropna()
            if len(series) < min_periods:
                raise ValueError(f"Insufficient data: {len(series)} periods (need {min_periods})")
            if (series <= 0).any():
                raise ValueError("Series contains non-positive values — log transform will fail")
            
            results[name] = decompose_category(series, period=period)
            
        except Exception as e:
            failures.append({'series': name, 'error': str(e)})
            print(f"  FAILED: {name} — {e}")
    
    if failures:
        failures_df = pd.DataFrame(failures)
        print(f"\n{len(failures)} decomposition failures:")
        print(failures_df)
    
    return results, failures

category_series_dict = {
    cat: multi_df.xs(cat, level='category')['revenue'] 
    for cat in categories
}

results, failures = batch_decompose(category_series_dict)
print(f"\nSuccessful decompositions: {len(results)}")

---

Building the Report Output

Decomposition results are only useful if they reach decision-makers. Let's build a formatted Excel output that a finance team can actually use, combining the raw decomposition data with the business-ready summary stats.

import openpyxl
from openpyxl.styles import PatternFill, Font, Alignment, Border, Side
from openpyxl.utils import get_column_letter

def build_decomposition_report(stl_components, forecast_df, seasonal_index_df, filename):
    """
    Export decomposition results to a formatted multi-sheet Excel workbook.
    """
    
    with pd.ExcelWriter(filename, engine='openpyxl') as writer:
        
        # Sheet 1: Full decomposition data
        full_data = pd.DataFrame({
            'Month': stl_components.index.strftime('%b %Y'),
            'Observed Revenue': stl_components['observed'].round(0),
            'Trend': stl_components['trend'].round(0),
            'Seasonal Factor': stl_components['seasonal'].round(4),
            'Seasonal % Effect': ((stl_components['seasonal'] - 1) * 100).round(1),
            'Residual': stl_components['residual'].round(4),
            'Seasonally Adjusted': (
                stl_components['observed'] / stl_components['seasonal']
            ).round(0)
        })
        full_data.to_excel(writer, sheet_name='Decomposition Detail', index=False)
        
        # Sheet 2: Seasonal index summary
        seasonal_summary = pd.DataFrame({
            'Month': ['Jan','Feb','Mar','Apr','May','Jun',
                     'Jul','Aug','Sep','Oct','Nov','Dec'],
            'Seasonal Factor': [
                stl_components['seasonal'].groupby(
                    stl_components.index.month
                ).mean()[m] for m in range(1, 13)
            ],
        })
        seasonal_summary['Seasonal Factor'] = seasonal_summary['Seasonal Factor'].round(4)
        seasonal_summary['% Above/Below Trend'] = (
            (seasonal_summary['Seasonal Factor'] - 1) * 100
        ).round(1)
        seasonal_summary['Interpretation'] = seasonal_summary['% Above/Below Trend'].apply(
            lambda x: f"{'+'}{x:.1f}% above trend" if x > 0 else f"{x:.1f}% below trend"
        )
        seasonal_summary.to_excel(writer, sheet_name='Seasonal Index', index=False)
        
        # Sheet 3: 2024 Forecast
        forecast_export = forecast_df.copy()
        forecast_export['projected_trend'] = forecast_export['projected_trend'].round(0)
        forecast_export['forecast_revenue'] = forecast_export['forecast_revenue'].round(0)
        forecast_export.to_excel(writer, sheet_name='2024 Forecast', index=False)
    
    # Apply formatting using openpyxl directly
    wb = openpyxl.load_workbook(filename)
    
    header_fill = PatternFill(start_color='1F4E79', end_color='1F4E79', fill_type='solid')
    header_font = Font(color='FFFFFF', bold=True, size=11)
    
    for sheet_name in wb.sheetnames:
        ws = wb[sheet_name]
        
        # Format header row
        for cell in ws[1]:
            cell.fill = header_fill
            cell.font = header_font
            cell.alignment = Alignment(horizontal='center')
        
        # Auto-size columns
        for col in ws.columns:
            max_length = max(len(str(cell.value or '')) for cell in col)
            ws.column_dimensions[get_column_letter(col[0].column)].width = min(max_length + 4, 30)
        
        # Freeze top row
        ws.freeze_panes = 'A2'
    
    wb.save(filename)
    print(f"Report saved: {filename}")

# Build the forecast DataFrame for export
forecast_dates = pd.date_range('2024-01-01', periods=12, freq='MS')
projected_trend_values = intercept + slope * np.arange(
    len(recent_trend), len(recent_trend) + 12
)
monthly_seasonal_factors = [monthly_seasonal[m] for m in range(1, 13)]

forecast_export_df = pd.DataFrame({
    'date': forecast_dates.strftime('%b %Y'),
    'projected_trend': projected_trend_values,
    'seasonal_factor': monthly_seasonal_factors,
    'forecast_revenue': projected_trend_values * monthly_seasonal_factors
})

build_decomposition_report(
    stl_components=stl_components,
    forecast_df=forecast_export_df,
    seasonal_index_df=seasonal_summary if 'seasonal_summary' in dir() else pd.DataFrame(),
    filename='seasonal_decomposition_report.xlsx'
)

For more advanced Excel formatting options, including conditional formatting and named ranges, the lesson on Automating Excel Reports with pandas and openpyxl: Formatted Workbooks Without Manual Work goes deep on the openpyxl API.

---

Hands-On Exercise

Work through this exercise using a publicly available dataset. The U.S. Census Bureau publishes retail sales data monthly (FRED: RSXFS — Advance Retail Sales). Download it as a CSV and complete the following:

Part 1: Data Preparation (15 min)

1. Load the retail sales CSV and set the date column as the DatetimeIndex with monthly frequency
2. Validate that no months are missing using the gap-detection code from earlier in this lesson
3. Create a chart of the raw series and annotate it with your observations about trend direction and seasonal behavior

Part 2: Decomposition (20 min) 4. Run seasonal_decompose with model='multiplicative' and period=12 5. Extract the seasonal factors and express them as percentage deviation from trend 6. Identify the three strongest and three weakest months based on seasonal factors 7. Run STL decomposition with robust=True and compare the residuals between the two methods

Part 3: Residual Analysis (15 min) 8. Plot the ACF of residuals from both methods 9. Identify any months where the absolute residual z-score exceeds 2.5 — these are your "surprise" months where actual performance significantly diverged from model expectations 10. Look up what happened in those months in the retail industry context

Part 4: Seasonal Forecast (20 min) 11. Fit a linear trend to the STL trend component using scipy.stats.linregress 12. Project the trend 12 months forward and apply monthly seasonal factors 13. Compare your naive seasonal forecast to the actual data for a holdout period (if your dataset extends far enough, use the last 12 months as holdout: fit on all data before, forecast those 12 months, compute MAPE)

Challenge: Run the full analysis for two different retail sub-categories (e.g., food and beverage vs. electronics). Compare their seasonal indices side by side. Write a 3-paragraph narrative explaining the key differences in their seasonal patterns and what this would imply for inventory planning.

---

Common Mistakes & Troubleshooting

"My trend component is all NaN"
You're almost certainly missing the extrapolate_trend='freq' argument, or your data has genuine missing periods that weren't filled before decomposition. Check df.isna().sum() and df.index.freq.

"Seasonal factors sum to zero / multiply to one but my data clearly has a strong seasonal pattern"
The seasonal period might be wrong. If you have weekly data but daily observations, period should be 7, not 30. If you have monthly data with quarterly seasonality (e.g., B2B contracts that renew quarterly), set period=3, not period=12. The period argument defines the length of one seasonal cycle in your observation units.

"STL raises 'Input series is not monotonic' or ValueError about length"
Your DatetimeIndex has either duplicates or non-uniform spacing. Use df.index.is_monotonic_increasing to check monotonicity and df.index.duplicated().sum() to check for duplicates. Fix with df = df[~df.index.duplicated(keep='last')].sort_index().

"My residuals show a strong pattern at lag 12 in the ACF"
The seasonal component wasn't fully extracted. This usually means your data has a more complex seasonal structure — possibly a 52-week pattern masquerading as a 12-month pattern, or seasonal factors that genuinely evolve over time. Try increasing the seasonal parameter in STL to allow more flexibility, or consider a sub-period seasonal window.

"Decomposition looks great for older data but terrible for the most recent 6 months"
Remember that classical decomposition loses half a period at each end of the trend. The most recent 6 months of your trend estimate may be extrapolated. STL doesn't have this exact problem but still has boundary effects. When forecasting, always use trend estimates from the interior of your series for the slope calculation, not from the most recent extrapolated endpoints.

"My multiplicative decomposition produces negative seasonal factors"
This is a sign that your series has zero or near-zero values, or actual sign changes. Multiplicative models assume all values are positive — they're undefined for zero. Check for zeros with (df['revenue'] == 0).sum() and decide whether to replace them with a small positive constant, impute from neighbors, or use additive decomposition instead.

"Running decomposition on 500 product series is too slow"
Consider using joblib.Parallel to parallelize across series. For very large numbers of series, also consider whether you need full STL (which involves iterative LOESS fitting) or whether a faster approximation suffices. The x13arima-seats approach via statsmodels.tsa.x13 is an alternative that may handle batch processing differently. If memory is the constraint rather than compute, the lesson on Handling Large Datasets in Python: Chunked Reading, Efficient dtypes, and When to Use Polars covers strategies for chunked and memory-efficient processing.

---

Summary & Next Steps

You've covered the full arc of time series decomposition for business reporting: understanding the mathematical framework, implementing classical and STL decomposition in pandas, diagnosing result quality via residual analysis, scaling to multi-series scenarios, and producing formatted report outputs.

The core workflow you should commit to muscle memory:

1. Prepare the series (DatetimeIndex, declared frequency, filled gaps)
2. Choose the model (additive vs. multiplicative; classical vs. STL)
3. Run decomposition and extract components
4. Validate with residual analysis (distribution, time plot, ACF)
5. Extract seasonal factors and express them in business language
6. Apply to trend projection for naive seasonal forecasts
7. Export cleanly for stakeholder consumption

Decomposition is a foundational skill that underlies nearly all business time series work. Once you understand it, you start seeing it everywhere — in how the Federal Reserve seasonally adjusts economic data, in how retailers set monthly targets, in how operations teams plan staffing.

Where to go next:

The natural extension of this work is incorporating it into automated recurring report pipelines. The lesson on Building and Automating Recurring Reports with pandas: Scheduling Scripts to Run Without You shows you how to schedule the full workflow — data loading, decomposition, and report export — to run on a cron job or task scheduler without manual intervention. When your seasonal forecast needs to compare against actuals as each new month comes in, the patterns in Comparing Periods and Calculating Month-over-Month, Year-over-Year, and Rolling Changes in pandas with shift and pct_change will complete the picture.