Week Over Week: A Comparison That Moves Five Percent on Its Own

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np

RNG = np.random.default_rng(83)
DAYS = 560
LEVEL = 1000.0
WEEKDAY = np.array([1.05, 1.08, 1.06, 1.04, 1.10, 0.82, 0.85])   # Mon to Sun
RHO, SD = 0.55, 0.035                                            # day-to-day noise


def series(days=DAYS, step_at=None, step=0.0, rng=RNG):
    """A daily metric with weekday seasonality, autocorrelated noise and an
    optional permanent step change."""
    eps = rng.normal(0, SD * np.sqrt(1 - RHO ** 2), days)
    z = np.zeros(days)
    for i in range(1, days):
        z[i] = RHO * z[i - 1] + eps[i]
    level = np.full(days, LEVEL)
    if step_at is not None:
        level[step_at:] *= (1 + step)
    return level * WEEKDAY[np.arange(days) % 7] * (1 + z)


def wow_day(y, t):
    """Today against the same weekday last week."""
    return y[t] / y[t - 7] - 1


def wow_mean(y, t, w=7):
    """The last seven days against the seven before them."""
    return y[t - w + 1:t + 1].mean() / y[t - 2 * w + 1:t - w + 1].mean() - 1


def yoy_day(y, t):
    """Today against the nearest same weekday a year earlier."""
    return y[t] / y[t - 364] - 1


r = np.random.default_rng(3)
d1, d7, dyy = [], [], []
for _ in range(400):
    y = series(rng=r)
    for t in range(400, DAYS):
        d1.append(wow_day(y, t))
        d7.append(wow_mean(y, t))
        dyy.append(yoy_day(y, t))
for name, v in (("same weekday last week", d1), ("seven-day averages", d7),
                ("same weekday last year", dyy)):
    v = np.array(v)
    print(f"{name:24s} spread {v.std():6.2%}, "
          f"exceeds 5% in {np.mean(np.abs(v) > 0.05):5.1%} of days, "
          f"exceeds 10% in {np.mean(np.abs(v) > 0.10):5.1%}")

1
2
3
4
5
6
for name, v in (("same weekday last week", d1), ("seven-day averages", d7),
                ("same weekday last year", dyy)):
    v = np.array(v)
    print(f"{name:24s} a 5% false alarm rate needs a threshold of "
          f"{np.quantile(np.abs(v), 0.95):5.2%}, a 1% rate needs "
          f"{np.quantile(np.abs(v), 0.99):5.2%}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
STEP_AT = 450
WINDOW = 28


def rules(y, thresholds):
    """Each rule as a function of the day index, at its own threshold."""
    pattern = np.array([y[:STEP_AT - 20][np.arange(STEP_AT - 20) % 7 == k].mean()
                        for k in range(7)])
    adj = y / pattern[np.arange(len(y)) % 7]
    base = adj[200:STEP_AT - 20]
    return {
        "same weekday": lambda t: abs(wow_day(y, t)) > thresholds["same weekday"],
        "seven-day means": lambda t: abs(wow_mean(y, t)) > thresholds["seven-day means"],
        "seasonal residual": lambda t: (abs(adj[t] - base.mean())
                                        > thresholds["seasonal residual"] * base.std()),
    }


def false_alarm_rate(thresholds, runs=200, seed=7):
    r = np.random.default_rng(seed)
    hits = {k: 0 for k in thresholds}
    for _ in range(runs):
        y = series(rng=r)
        fns = rules(y, thresholds)
        for k, fn in fns.items():
            hits[k] += any(fn(t) for t in range(STEP_AT, STEP_AT + WINDOW))
    return {k: v / runs for k, v in hits.items()}


# Tune each threshold until a quiet four-week window trips it one time in ten.
thresholds = {"same weekday": 0.10, "seven-day means": 0.06, "seasonal residual": 3.0}
for key, lo, hi in (("same weekday", 0.05, 0.40), ("seven-day means", 0.02, 0.30),
                    ("seasonal residual", 1.5, 6.0)):
    for _ in range(12):
        mid = (lo + hi) / 2
        trial = dict(thresholds, **{key: mid})
        if false_alarm_rate(trial, runs=120)[key] > 0.10:
            lo = mid
        else:
            hi = mid
    thresholds[key] = (lo + hi) / 2
rates = false_alarm_rate(thresholds, runs=300)
for k, v in thresholds.items():
    unit = "standard deviations" if k == "seasonal residual" else "percent"
    shown = f"{v:.2f} {unit}" if k == "seasonal residual" else f"{v:.2%}"
    print(f"{k:18s} threshold {shown:24s} false alarms {rates[k]:5.1%} per quiet month")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
for step in (0.02, 0.03, 0.05, 0.10):
    r = np.random.default_rng(5)
    found = {k: [] for k in thresholds}
    for _ in range(300):
        y = series(step_at=STEP_AT, step=step, rng=r)
        fns = rules(y, thresholds)
        for k, fn in fns.items():
            hit = next((t - STEP_AT for t in range(STEP_AT, DAYS) if fn(t)), None)
            found[k].append(hit if hit is not None else DAYS - STEP_AT)
    parts = []
    for k, v in found.items():
        v = np.array(v)
        seen = v < 28
        median = np.median(v[seen]) if seen.any() else float("nan")
        parts.append(f"{k} {seen.mean():5.1%} within four weeks"
                     + (f", median {median:4.1f} days" if seen.any() else ""))
    print(f"step of {step:4.0%}: " + "; ".join(parts))

1
2
3
4
5
6
7
8
9
10
r = np.random.default_rng(11)
singles, means = [], []
for _ in range(2000):
    y = series(days=120, rng=r)
    singles.append(wow_day(y, 100))
    means.append(wow_mean(y, 100))
print(f"single day comparison spread {np.std(singles):.3%}")
print(f"seven-day average comparison spread {np.std(means):.3%}")
print(f"ratio {np.std(singles) / np.std(means):.2f}, against {np.sqrt(7):.2f} "
      f"if the days were independent")