Designing an EWMA Chart for Slow Drifts in Python

This is the step-by-step recipe for building an exponentially weighted moving average (EWMA) chart that reliably catches a slow drift in the process mean — a tool wearing a few microns per shift, a bath concentration creeping off target — the kind of shift a Shewhart chart takes dozens of samples to notice. It is the implementation companion to EWMA and CUSUM charts for small-shift detection: here we focus solely on the EWMA half, getting the time-varying limits exactly right, choosing the smoothing constant and limit width for a stated target, and verifying the chart against a fixture with a known injected shift.

Prerequisites

Before you compute a single EWMA point, confirm the upstream state:

  • Python 3.9+ with numpy and pandas installed (pip install numpy pandas).
  • A frozen Phase I baseline: an in-control mean $\mu_0$ and standard deviation $\sigma$ from at least 20–30 stable observations with no known assignable causes.
  • Observations are approximately independent — lag-1 autocorrelation below about 0.25. On correlated data the limits are too tight; model the series first, as you would before an I-MR chart.
  • A stated target shift $\delta$ (in σ units) you need to detect quickly — this drives the choice of $\lambda$.
  • Data arrives as a numeric pandas Series in observation order, with nulls already resolved upstream.

Step 1 — Fix the baseline and the design target

The EWMA is a Phase II monitor, so $\mu_0$ and $\sigma$ are inputs, never recomputed from the live stream. Capture them as immutable values alongside the design target — the shift $\delta$ you care about and the in-control run length you can tolerate. Writing these down first makes the parameter choice in Step 2 a derivation rather than a guess.

from dataclasses import dataclass


@dataclass(frozen=True)
class EwmaDesign:
    """Frozen EWMA design: baseline plus tuning target."""
    mu0: float          # in-control mean (Phase I)
    sigma: float        # in-control std dev (Phase I)
    lam: float = 0.2    # smoothing constant lambda in (0, 1]
    L: float = 3.0      # limit width in sigma units

    def __post_init__(self) -> None:
        if self.sigma <= 0:
            raise ValueError("sigma must be positive.")
        if not 0.0 < self.lam <= 1.0:
            raise ValueError("lambda must be in (0, 1].")
        if self.L <= 0:
            raise ValueError("L must be positive.")

Step 2 — Choose lambda and L for the target shift

The smoothing constant $\lambda$ sets the chart's memory: smaller $\lambda$ averages over more history and is more sensitive to small, slow drifts, while larger $\lambda$ reacts faster to big jumps and, at $\lambda = 1$, reduces to a Shewhart individuals chart. The limit width $L$ controls the false-alarm rate. The well-established ARL-optimal pairings from Montgomery's tables are the safe starting points: use a small $\lambda$ with $L$ slightly below 3 so the in-control run length stays near the Shewhart benchmark of about 370.

Target shift $\delta$ $\lambda$ $L$ Rationale
≈ 0.5σ (very slow drift) 0.05–0.10 2.6–2.8 Long memory; smaller L keeps in-control ARL from ballooning
≈ 1.0σ (moderate drift) 0.10–0.20 2.7–3.0 The workhorse default; in-control ARL ≈ 370
≈ 1.5–2σ (faster shift) 0.25–0.40 3.0 Shorter memory; approaches Shewhart behavior 
def recommend_lambda_L(delta_sigma: float) -> tuple[float, float]:
    """Return a (lambda, L) starting pair for a target shift in sigma units."""
    if delta_sigma <= 0:
        raise ValueError("delta_sigma must be positive.")
    if delta_sigma <= 0.75:
        return 0.10, 2.70
    if delta_sigma <= 1.25:
        return 0.20, 2.86
    return 0.40, 3.00

Treat these as informed defaults, not gospel; if you have an ARL calculator or a Markov-chain routine, confirm the in-control ARL for your exact pair. The key discipline is that $\lambda$ follows from the shift you want to detect, not from habit.

Step 3 — Compute the EWMA statistic

The statistic is the recursion $z_t = \lambda x_t + (1-\lambda)z_{t-1}$, seeded at $z_0 = \mu_0$. Each new value blends in with weight $\lambda$; older values decay geometrically.

import numpy as np
import pandas as pd


def ewma_statistic(values: pd.Series, design: EwmaDesign) -> np.ndarray:
    """Compute the EWMA statistic z_t seeded at mu0."""
    x = np.asarray(values, dtype=np.float64)
    if x.size == 0:
        raise ValueError("Input series is empty.")
    if np.isnan(x).any():
        raise ValueError("Resolve missing values upstream before charting.")

    z = np.empty_like(x)
    prev = design.mu0
    for t in range(x.size):
        prev = design.lam * x[t] + (1.0 - design.lam) * prev
        z[t] = prev
    return z

Step 4 — Build the exact time-varying limits

The variance of $z_t$ grows from zero toward a steady state, so the limits widen over the first several samples. Use the exact form with the $\left(1-(1-\lambda)^{2t}\right)$ factor rather than the flat asymptote — otherwise the earliest points are held to a band that is too wide and a fast early drift slips through:

$$\text{UCL}_t,\ \text{LCL}_t = \mu_0 \pm L\,\sigma\sqrt{\frac{\lambda}{2-\lambda}\left(1-(1-\lambda)^{2t}\right)}$$

def ewma_limits(n: int, design: EwmaDesign) -> tuple[np.ndarray, np.ndarray]:
    """Exact time-varying UCL/LCL arrays for n points (t = 1..n)."""
    if n < 1:
        raise ValueError("n must be at least 1.")
    t = np.arange(1, n + 1)
    steady = np.sqrt(design.lam / (2.0 - design.lam))
    width = steady * np.sqrt(1.0 - (1.0 - design.lam) ** (2 * t))
    half = design.L * design.sigma * width
    return design.mu0 + half, design.mu0 - half

The limits reach 99% of their steady-state width after roughly $t \approx \ln(0.005)/\big(2\ln(1-\lambda)\big)$ samples — about 12 points at $\lambda = 0.2$ — after which the flat asymptotic band is a fine approximation. Automating the exact form costs nothing and removes the early-sample edge case entirely.

Step 5 — Flag signals

A signal is any point where the EWMA statistic falls outside its own time-matched limit. Because $z_t$ is smoothed, a signal reflects sustained drift, not a lone spike.

def ewma_chart(values: pd.Series, design: EwmaDesign) -> pd.DataFrame:
    """Assemble the EWMA statistic, limits, and per-point signal flags."""
    z = ewma_statistic(values, design)
    ucl, lcl = ewma_limits(len(z), design)
    out = pd.DataFrame(
        {"x": np.asarray(values, dtype=np.float64), "z": z, "ucl": ucl, "lcl": lcl},
        index=getattr(values, "index", pd.RangeIndex(len(z))),
    )
    out["signal"] = (out["z"] > out["ucl"]) | (out["z"] < out["lcl"])
    return out

Verification

Prove the chart on a fixture with a known injected shift: a constant in-control run must stay silent, and a sustained drift must signal within a few samples of crossing. Assert both directions.

import numpy as np
import pandas as pd


def test_ewma_flags_known_shift():
    design = EwmaDesign(mu0=100.0, sigma=2.0, lam=0.2, L=2.86)

    # 20 in-control points, then a sustained +1.5 sigma (=+3.0 unit) shift.
    x = np.concatenate([np.full(20, 100.0), np.full(20, 103.0)])
    result = ewma_chart(pd.Series(x), design)

    # No signal while in control.
    assert not result["signal"].iloc[:20].any()
    # Signals within ~10 samples of the shift onset.
    first = int(result["signal"].values.argmax())
    assert 20 <= first <= 30

    # Limits widen monotonically toward the steady-state band, then flatten.
    ucl = result["ucl"].to_numpy()
    assert ucl[0] < ucl[5] < ucl[15]
    assert np.isclose(ucl[-1], ucl[-2], atol=1e-6)


def test_deadon_series_is_silent():
    design = EwmaDesign(mu0=50.0, sigma=1.5)
    result = ewma_chart(pd.Series(np.full(60, 50.0)), design)
    assert result["signal"].sum() == 0
    assert np.allclose(result["z"], 50.0)     # z stays pinned to mu0

The monotonic-limit assertion is the one people forget: if ucl[0] is already at the steady-state value, you have silently used the asymptotic band and the first-sample sensitivity is wrong. The dead-on test confirms the seed $z_0 = \mu_0$ keeps the statistic pinned when nothing drifts.

Root-Cause Table

Symptom Cause Fix
Slow drift never signals $\lambda$ too large — chart behaves like a Shewhart individuals chart Lower $\lambda$ toward 0.1–0.2 to lengthen the memory
Frequent false alarms in control $L$ too small, or observations autocorrelated Raise $L$ toward 3.0; check lag-1 ACF and model the series if needed
Earliest points over-sensitive Flat asymptotic limits used from $t=1$ Use the exact $\left(1-(1-\lambda)^{2t}\right)$ time-varying form
Statistic drifts and never returns Baseline $\mu_0$ recomputed on shifted live data Freeze $\mu_0,\sigma$ in Phase I; recalibrate only after a verified change
z diverges or is NaN Missing values propagated into the recursion Resolve nulls upstream; the function raises on NaN

For chart selection across every variable and attribute chart, see SPC Fundamentals & Control Chart Taxonomy.