Skip to content

mfe.bootstrap

mfe.bootstrap

mfe.bootstrap — Dependent-data bootstrap methods.

wild_bootstrap_rv Wild bootstrap CI for realized volatility statistics wild_bootstrap_test Wild bootstrap p-value for hypothesis tests spa_test Hansen (2005) Superior Predictive Ability test step_m Romano-Wolf (2005) Stepdown Multiple Hypothesis Test (FWER control)

WildBootstrapResult dataclass 

WildBootstrapResult(statistic: float, ci_lower: float, ci_upper: float, ci_level: float, bootstrap_distribution: FloatArray, n_replications: int, multiplier: str)

Result from a wild bootstrap confidence interval computation.

StepMResult dataclass 

StepMResult(rejected: list[int], accepted: list[int], t_stats: FloatArray, p_values_raw: FloatArray, p_values_adjusted: FloatArray, n_models: int, n_obs: int, n_bootstrap: int, alpha: float)

StepM multiple hypothesis testing result.

wild_bootstrap_rv 

wild_bootstrap_rv(returns: FloatArray, statistic_fn: Callable[[FloatArray], float] | None = None, n_replications: int = 999, ci_level: float = 0.95, multiplier: str = 'rademacher', rng: Generator | None = None) -> WildBootstrapResult

Wild bootstrap confidence interval for a realized volatility statistic.

The bootstrap DGP is: r_t^* = w_t * r_t

where w_t is i.i.d. from the specified multiplier distribution. The statistic is re-evaluated on {r_t^*}.

Parameters:

Name Type Description Default
returns FloatArray
required
statistic_fn Callable[[FloatArray], float] | None
None
n_replications number of bootstrap replications
999
ci_level float
0.95
multiplier str
'rademacher'
rng Generator | None
None

Returns:

Type Description
WildBootstrapResult
Source code in src/mfe/bootstrap/wild.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def wild_bootstrap_rv(
    returns: FloatArray,
    statistic_fn: Callable[[FloatArray], float] | None = None,
    n_replications: int = 999,
    ci_level: float = 0.95,
    multiplier: str = "rademacher",
    rng: np.random.Generator | None = None,
) -> WildBootstrapResult:
    """
    Wild bootstrap confidence interval for a realized volatility statistic.

    The bootstrap DGP is:
        r_t^* = w_t * r_t

    where w_t is i.i.d. from the specified multiplier distribution.
    The statistic is re-evaluated on {r_t^*}.

    Parameters
    ----------
    returns        : (T,) log-return array
    statistic_fn   : function mapping returns -> float; default is sum(r^2) (RV)
    n_replications : number of bootstrap replications
    ci_level       : confidence level (e.g. 0.95 for 95% CI)
    multiplier     : "rademacher" | "mammen" | "normal"
    rng            : numpy random generator; if None, uses default_rng()

    Returns
    -------
    WildBootstrapResult
    """
    r = np.asarray(returns, dtype=np.float64)
    T = len(r)

    if rng is None:
        rng = np.random.default_rng()

    if statistic_fn is None:
        def statistic_fn(x: FloatArray) -> float:
            return float(np.sum(x ** 2))

    mult_func = _MULTIPLIER_FUNCS.get(multiplier)
    if mult_func is None:
        raise ValueError(f"multiplier must be one of {list(_MULTIPLIER_FUNCS)}, got '{multiplier}'")

    # Point estimate
    stat0 = statistic_fn(r)

    # Bootstrap distribution
    boot_stats = np.empty(n_replications, dtype=np.float64)
    for b in range(n_replications):
        w = mult_func(T, rng)
        r_star = w * r
        boot_stats[b] = statistic_fn(r_star)

    alpha = 1.0 - ci_level
    ci_lo = float(np.percentile(boot_stats, 100 * alpha / 2))
    ci_hi = float(np.percentile(boot_stats, 100 * (1 - alpha / 2)))

    return WildBootstrapResult(
        statistic=stat0,
        ci_lower=ci_lo,
        ci_upper=ci_hi,
        ci_level=ci_level,
        bootstrap_distribution=boot_stats,
        n_replications=n_replications,
        multiplier=multiplier,
    )

wild_bootstrap_test 

wild_bootstrap_test(returns: FloatArray, null_statistic: float, statistic_fn: Callable[[FloatArray], float] | None = None, n_replications: int = 999, multiplier: str = 'rademacher', rng: Generator | None = None) -> tuple[float, float]

Wild bootstrap p-value for a two-sided hypothesis test.

Parameters:

Name Type Description Default
null_statistic the value of the statistic under the null hypothesis
required

Returns:

Type Description
(observed_statistic, bootstrap_p_value)
Source code in src/mfe/bootstrap/wild.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def wild_bootstrap_test(
    returns: FloatArray,
    null_statistic: float,
    statistic_fn: Callable[[FloatArray], float] | None = None,
    n_replications: int = 999,
    multiplier: str = "rademacher",
    rng: np.random.Generator | None = None,
) -> tuple[float, float]:
    """
    Wild bootstrap p-value for a two-sided hypothesis test.

    Parameters
    ----------
    null_statistic : the value of the statistic under the null hypothesis

    Returns
    -------
    (observed_statistic, bootstrap_p_value)
    """
    result = wild_bootstrap_rv(
        returns,
        statistic_fn=statistic_fn,
        n_replications=n_replications,
        multiplier=multiplier,
        rng=rng,
    )
    p_val = float(np.mean(np.abs(result.bootstrap_distribution - null_statistic) >=
                          abs(result.statistic - null_statistic)))
    return result.statistic, p_val

spa_test 

spa_test(loss_benchmark: FloatArray, loss_models: FloatArray, n_bootstrap: int = 999, avg_block_len: float | None = None, bandwidth: int | None = None, rng: Generator | None = None) -> SPAResult

Hansen (2005) Superior Predictive Ability test.

Parameters:

Name Type Description Default
loss_benchmark (T,) loss series for the benchmark model (lower = better)
required
loss_models FloatArray
required
n_bootstrap int
999
avg_block_len float | None 
         if None, uses T^{1/3}
None
bandwidth int | None 
         if None, uses 1.2 * T^{1/3}
None
rng Generator | None
None

Returns:

Type Description
SPAResult

Report .p_value_consistent for the standard SPA p-value. Report .p_value_upper for White's Reality Check p-value.
Notes

Loss convention: LOWER is BETTER (e.g. MSE, MAE, negative log-lik). Loss differential d_{k,t} = L_benchmark_t - L_model_k_t. Positive d_bar_k means model k beats benchmark on average.

Source code in src/mfe/bootstrap/spa.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def spa_test(
    loss_benchmark: FloatArray,
    loss_models: FloatArray,
    n_bootstrap: int = 999,
    avg_block_len: float | None = None,
    bandwidth: int | None = None,
    rng: np.random.Generator | None = None,
) -> SPAResult:
    """
    Hansen (2005) Superior Predictive Ability test.

    Parameters
    ----------
    loss_benchmark : (T,) loss series for the benchmark model (lower = better)
    loss_models    : (T, M) loss series for M alternative models
    n_bootstrap    : stationary bootstrap replications (default 999)
    avg_block_len  : average block length for stationary bootstrap;
                     if None, uses T^{1/3}
    bandwidth      : Newey-West bandwidth for long-run variance;
                     if None, uses 1.2 * T^{1/3}
    rng            : numpy Generator; if None uses default_rng()

    Returns
    -------
    SPAResult
        Report .p_value_consistent for the standard SPA p-value.
        Report .p_value_upper for White's Reality Check p-value.

    Notes
    -----
    Loss convention: LOWER is BETTER (e.g. MSE, MAE, negative log-lik).
    Loss differential d_{k,t} = L_benchmark_t - L_model_k_t.
    Positive d_bar_k means model k beats benchmark on average.
    """
    lb = np.asarray(loss_benchmark, dtype=np.float64)
    lm = np.asarray(loss_models, dtype=np.float64)
    if lm.ndim == 1:
        lm = lm[:, None]

    T, M = lm.shape
    if rng is None:
        rng = np.random.default_rng()
    if avg_block_len is None:
        avg_block_len = max(2.0, float(T ** (1 / 3)))

    # Loss differentials: d_{k,t} = L_bench_t - L_model_k_t
    # d_bar_k > 0 means model k is better than benchmark
    d = lb[:, None] - lm   # (T, M)

    d_bar = d.mean(axis=0)  # (M,)

    # Long-run variance
    sigma2 = _long_run_variance(d, bandwidth=bandwidth)
    sigma = np.sqrt(sigma2)

    # Studentized statistics: t_k = sqrt(T) * d_bar_k / sigma_k
    t_stats = np.sqrt(T) * d_bar / sigma   # (M,)
    T_spa = float(np.max(t_stats))

    # Stationary bootstrap to get null distribution of max t_k
    idx = _stationary_bootstrap_indices(T, n_bootstrap, avg_block_len, rng)

    # Three null variants of Hansen (2005)
    # "consistent": zero out models with strongly negative d_bar (irrelevant)
    # "upper":      keep all models (White RC)
    # "lower":      only models with d_bar > 0 (conservative)

    # Threshold for "consistent": c_k = max(d_bar_k, -sqrt(sigma2_k * log(log(T)) / T))
    c_consistent = np.maximum(d_bar, -np.sqrt(sigma2 * np.log(np.log(T)) / T))
    c_upper = d_bar.copy()          # White Reality Check: mean-center on d_bar
    c_lower = np.maximum(d_bar, 0.0)

    boot_max_consistent = np.empty(n_bootstrap, dtype=np.float64)
    boot_max_upper = np.empty(n_bootstrap, dtype=np.float64)
    boot_max_lower = np.empty(n_bootstrap, dtype=np.float64)

    d_centered_consistent = d - c_consistent[None, :]
    d_centered_upper      = d - c_upper[None, :]
    d_centered_lower      = d - c_lower[None, :]

    for b in range(n_bootstrap):
        d_boot_c = d_centered_consistent[idx[b]].mean(axis=0)
        d_boot_u = d_centered_upper[idx[b]].mean(axis=0)
        d_boot_l = d_centered_lower[idx[b]].mean(axis=0)

        t_boot_c = np.sqrt(T) * d_boot_c / sigma
        t_boot_u = np.sqrt(T) * d_boot_u / sigma
        t_boot_l = np.sqrt(T) * d_boot_l / sigma

        boot_max_consistent[b] = float(np.max(t_boot_c))
        boot_max_upper[b]      = float(np.max(t_boot_u))
        boot_max_lower[b]      = float(np.max(t_boot_l))

    pval_consistent = float(np.mean(boot_max_consistent >= T_spa))
    pval_upper      = float(np.mean(boot_max_upper >= T_spa))
    pval_lower      = float(np.mean(boot_max_lower >= T_spa))

    return SPAResult(
        statistic=T_spa,
        p_value_consistent=pval_consistent,
        p_value_upper=pval_upper,
        p_value_lower=pval_lower,
        d_bar=d_bar,
        t_stats=t_stats,
        n_models=M,
        n_obs=T,
        n_bootstrap=n_bootstrap,
        bootstrap_distribution=boot_max_consistent,
    )

step_m 

step_m(loss_benchmark: FloatArray, loss_models: FloatArray, alpha: float = 0.05, n_bootstrap: int = 999, avg_block_len: float | None = None, bandwidth: int | None = None, rng: Generator | None = None) -> StepMResult

Romano-Wolf StepM stepdown multiple hypothesis test.

Tests H_k: E[L_bench - L_model_k] <= 0 for k = 1..M. Rejects H_k (model k beats benchmark) for k in result.rejected. Controls FWER <= alpha across all M tests.

Parameters:

Name Type Description Default
loss_benchmark (T,) benchmark loss series (lower = better)
required
loss_models FloatArray
required
alpha float
0.05
n_bootstrap int
999
avg_block_len float | None
None
bandwidth int | None
None
rng Generator | None
None

Returns:

Type Description
StepMResult

.rejected — 0-based indices of models significantly beating benchmark .accepted — the rest
Source code in src/mfe/bootstrap/stepM.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def step_m(
    loss_benchmark: FloatArray,
    loss_models: FloatArray,
    alpha: float = 0.05,
    n_bootstrap: int = 999,
    avg_block_len: float | None = None,
    bandwidth: int | None = None,
    rng: np.random.Generator | None = None,
) -> StepMResult:
    """
    Romano-Wolf StepM stepdown multiple hypothesis test.

    Tests H_k: E[L_bench - L_model_k] <= 0  for k = 1..M.
    Rejects H_k (model k beats benchmark) for k in result.rejected.
    Controls FWER <= alpha across all M tests.

    Parameters
    ----------
    loss_benchmark : (T,) benchmark loss series (lower = better)
    loss_models    : (T, M) alternative model loss series
    alpha          : familywise error rate (default 0.05)
    n_bootstrap    : stationary bootstrap replications
    avg_block_len  : average block length; if None uses T^{1/3}
    bandwidth      : Newey-West bandwidth; if None uses 1.2 * T^{1/3}
    rng            : random generator

    Returns
    -------
    StepMResult
        .rejected  — 0-based indices of models significantly beating benchmark
        .accepted  — the rest
    """
    lb = np.asarray(loss_benchmark, dtype=np.float64)
    lm = np.asarray(loss_models, dtype=np.float64)
    if lm.ndim == 1:
        lm = lm[:, None]

    T, M = lm.shape
    if rng is None:
        rng = np.random.default_rng()
    if avg_block_len is None:
        avg_block_len = max(2.0, T ** (1 / 3))
    if bandwidth is None:
        bandwidth = max(1, int(1.2 * T ** (1 / 3)))

    # Loss differentials: d_{k,t} = L_bench_t - L_model_k_t
    d = lb[:, None] - lm                   # (T, M)
    d_bar = d.mean(axis=0)                 # (M,)
    sigma = _long_run_std(d, bandwidth)    # (M,)
    t_stats = np.sqrt(T) * d_bar / sigma  # (M,)

    # Unadjusted p-values (individual, no FWER control)
    # Use bootstrap max distribution over all M models
    boot = _stationary_bootstrap(d - d_bar[None, :], n_bootstrap, avg_block_len, rng)
    # boot: (n_boot, T, M) — resampled centered loss diffs

    p_raw = np.empty(M, dtype=np.float64)
    boot_max_all = (np.sqrt(T) * boot.mean(axis=1) / sigma[None, :]).max(axis=1)
    for k in range(M):
        boot_k = np.sqrt(T) * boot[:, :, k].mean(axis=1) / sigma[k]
        p_raw[k] = float(np.mean(boot_k >= t_stats[k]))

    # Stepdown procedure
    remaining = list(range(M))
    rejected = []
    p_adjusted = np.ones(M, dtype=np.float64)

    step = 0
    while remaining:
        # Bootstrap max over remaining models
        boot_max = np.empty(n_bootstrap, dtype=np.float64)
        for b in range(n_bootstrap):
            t_boot_remaining = np.sqrt(T) * boot[b, :, :][:, remaining].mean(axis=0) / sigma[remaining]
            boot_max[b] = float(np.max(t_boot_remaining))

        # Critical value at level alpha
        cv = float(np.quantile(boot_max, 1 - alpha))

        # Find the model with max t_stat among remaining
        t_remaining = t_stats[remaining]
        max_idx_in_remaining = int(np.argmax(t_remaining))
        max_k = remaining[max_idx_in_remaining]
        max_t = float(t_stats[max_k])

        if max_t > cv:
            # Reject this model
            p_adjusted[max_k] = float(np.mean(boot_max >= max_t))
            rejected.append(max_k)
            remaining.remove(max_k)
            step += 1
        else:
            # No more rejections possible
            break

    # Adjusted p-values for accepted: use the last step's distribution
    # (conservative: bound by the p-value from the last step)
    for k in remaining:
        p_adjusted[k] = float(np.mean(boot_max_all >= t_stats[k]))

    # Monotonise: stepdown p-values must be non-decreasing when sorted by t_stat descending
    order = np.argsort(-t_stats)
    p_mono = p_adjusted[order].copy()
    for i in range(1, M):
        p_mono[i] = max(p_mono[i], p_mono[i - 1])
    p_adjusted[order] = p_mono

    accepted = [k for k in range(M) if k not in rejected]

    return StepMResult(
        rejected=sorted(rejected),
        accepted=sorted(accepted),
        t_stats=t_stats,
        p_values_raw=p_raw,
        p_values_adjusted=p_adjusted,
        n_models=M,
        n_obs=T,
        n_bootstrap=n_bootstrap,
        alpha=alpha,
    )

spa

Superior Predictive Ability (SPA) test — Hansen (2005).

Hansen, P.R. (2005): "A Test for Superior Predictive Ability", Journal of Business & Economic Statistics, 23(4), 365-380.

Also implements White (2000) Reality Check as a special case.

Setup

Given M models and a benchmark, let d_{k,t} = L(y_t, f_{0,t}) - L(y_t, f_{k,t}) be the loss differential at time t for model k vs. benchmark (model 0). d_{k,t} > 0 means model k is better than benchmark at time t.

H0: max_k E[d_{k,t}] <= 0 (no model beats the benchmark on average) H1: max_k E[d_{k,t}] > 0 (at least one model is strictly better)

Test statistic

T_SPA = max_k ( sqrt(T) * d_bar_k / sigma_k ) where d_bar_k = mean(d_{k,t}) and sigma_k^2 is the long-run variance of d_{k,t}.

Under H0, T_SPA has a distribution that depends on the correlation structure of {d_{k,t}} across k. P-values are computed by the stationary bootstrap.

Hansen's SPA uses a "studentized" version with three variants of the null: - "consistent" (default): removes irrelevant models from the null (d_bar_k << 0) - "upper": retains all models (equivalent to White's Reality Check) - "lower": most conservative, all models treated as tied with benchmark

References

White, H. (2000): "A Reality Check for Data Snooping", Econometrica.

spa_test 

spa_test(loss_benchmark: FloatArray, loss_models: FloatArray, n_bootstrap: int = 999, avg_block_len: float | None = None, bandwidth: int | None = None, rng: Generator | None = None) -> SPAResult

Hansen (2005) Superior Predictive Ability test.

Parameters:

Name Type Description Default
loss_benchmark (T,) loss series for the benchmark model (lower = better)
required
loss_models FloatArray
required
n_bootstrap int
999
avg_block_len float | None 
         if None, uses T^{1/3}
None
bandwidth int | None 
         if None, uses 1.2 * T^{1/3}
None
rng Generator | None
None

Returns:

Type Description
SPAResult

Report .p_value_consistent for the standard SPA p-value. Report .p_value_upper for White's Reality Check p-value.
Notes

Loss convention: LOWER is BETTER (e.g. MSE, MAE, negative log-lik). Loss differential d_{k,t} = L_benchmark_t - L_model_k_t. Positive d_bar_k means model k beats benchmark on average.

Source code in src/mfe/bootstrap/spa.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def spa_test(
    loss_benchmark: FloatArray,
    loss_models: FloatArray,
    n_bootstrap: int = 999,
    avg_block_len: float | None = None,
    bandwidth: int | None = None,
    rng: np.random.Generator | None = None,
) -> SPAResult:
    """
    Hansen (2005) Superior Predictive Ability test.

    Parameters
    ----------
    loss_benchmark : (T,) loss series for the benchmark model (lower = better)
    loss_models    : (T, M) loss series for M alternative models
    n_bootstrap    : stationary bootstrap replications (default 999)
    avg_block_len  : average block length for stationary bootstrap;
                     if None, uses T^{1/3}
    bandwidth      : Newey-West bandwidth for long-run variance;
                     if None, uses 1.2 * T^{1/3}
    rng            : numpy Generator; if None uses default_rng()

    Returns
    -------
    SPAResult
        Report .p_value_consistent for the standard SPA p-value.
        Report .p_value_upper for White's Reality Check p-value.

    Notes
    -----
    Loss convention: LOWER is BETTER (e.g. MSE, MAE, negative log-lik).
    Loss differential d_{k,t} = L_benchmark_t - L_model_k_t.
    Positive d_bar_k means model k beats benchmark on average.
    """
    lb = np.asarray(loss_benchmark, dtype=np.float64)
    lm = np.asarray(loss_models, dtype=np.float64)
    if lm.ndim == 1:
        lm = lm[:, None]

    T, M = lm.shape
    if rng is None:
        rng = np.random.default_rng()
    if avg_block_len is None:
        avg_block_len = max(2.0, float(T ** (1 / 3)))

    # Loss differentials: d_{k,t} = L_bench_t - L_model_k_t
    # d_bar_k > 0 means model k is better than benchmark
    d = lb[:, None] - lm   # (T, M)

    d_bar = d.mean(axis=0)  # (M,)

    # Long-run variance
    sigma2 = _long_run_variance(d, bandwidth=bandwidth)
    sigma = np.sqrt(sigma2)

    # Studentized statistics: t_k = sqrt(T) * d_bar_k / sigma_k
    t_stats = np.sqrt(T) * d_bar / sigma   # (M,)
    T_spa = float(np.max(t_stats))

    # Stationary bootstrap to get null distribution of max t_k
    idx = _stationary_bootstrap_indices(T, n_bootstrap, avg_block_len, rng)

    # Three null variants of Hansen (2005)
    # "consistent": zero out models with strongly negative d_bar (irrelevant)
    # "upper":      keep all models (White RC)
    # "lower":      only models with d_bar > 0 (conservative)

    # Threshold for "consistent": c_k = max(d_bar_k, -sqrt(sigma2_k * log(log(T)) / T))
    c_consistent = np.maximum(d_bar, -np.sqrt(sigma2 * np.log(np.log(T)) / T))
    c_upper = d_bar.copy()          # White Reality Check: mean-center on d_bar
    c_lower = np.maximum(d_bar, 0.0)

    boot_max_consistent = np.empty(n_bootstrap, dtype=np.float64)
    boot_max_upper = np.empty(n_bootstrap, dtype=np.float64)
    boot_max_lower = np.empty(n_bootstrap, dtype=np.float64)

    d_centered_consistent = d - c_consistent[None, :]
    d_centered_upper      = d - c_upper[None, :]
    d_centered_lower      = d - c_lower[None, :]

    for b in range(n_bootstrap):
        d_boot_c = d_centered_consistent[idx[b]].mean(axis=0)
        d_boot_u = d_centered_upper[idx[b]].mean(axis=0)
        d_boot_l = d_centered_lower[idx[b]].mean(axis=0)

        t_boot_c = np.sqrt(T) * d_boot_c / sigma
        t_boot_u = np.sqrt(T) * d_boot_u / sigma
        t_boot_l = np.sqrt(T) * d_boot_l / sigma

        boot_max_consistent[b] = float(np.max(t_boot_c))
        boot_max_upper[b]      = float(np.max(t_boot_u))
        boot_max_lower[b]      = float(np.max(t_boot_l))

    pval_consistent = float(np.mean(boot_max_consistent >= T_spa))
    pval_upper      = float(np.mean(boot_max_upper >= T_spa))
    pval_lower      = float(np.mean(boot_max_lower >= T_spa))

    return SPAResult(
        statistic=T_spa,
        p_value_consistent=pval_consistent,
        p_value_upper=pval_upper,
        p_value_lower=pval_lower,
        d_bar=d_bar,
        t_stats=t_stats,
        n_models=M,
        n_obs=T,
        n_bootstrap=n_bootstrap,
        bootstrap_distribution=boot_max_consistent,
    )

stepM

StepM: Stepdown Multiple Hypothesis Testing with FWER control.

Romano, J.P. & Wolf, M. (2005): "Stepwise Multiple Testing as Formalized Data Snooping", Econometrica, 73(4), 1237-1282.

Setup

M null hypotheses H_k: mu_k <= 0 for k = 1..M, where mu_k = E[d_{k,t}] is the mean performance differential of model k vs. the benchmark.

The algorithm controls the familywise error rate (FWER): FWER = P(reject at least one true H_k) <= alpha

This is more powerful than Bonferroni and more interpretable than SPA: it returns which models are significantly better, not just whether any is.

Algorithm (Algorithm 4.1 of Romano & Wolf 2005)
  1. Start with all M models.
  2. Compute test statistics t_k = sqrt(T) * d_bar_k / sigma_k.
  3. Use the stationary bootstrap to get the joint null distribution of max_k t_k (over the remaining models).
  4. Reject the model with the largest t_k if it exceeds the bootstrap critical value at level alpha.
  5. Remove rejected models from the set and repeat.
  6. Stop when no more rejections occur.

The result is a set of models significantly better than the benchmark.

This matches the MFE MATLAB implementation under bootstrap/stepm.m.

StepMResult dataclass 

StepMResult(rejected: list[int], accepted: list[int], t_stats: FloatArray, p_values_raw: FloatArray, p_values_adjusted: FloatArray, n_models: int, n_obs: int, n_bootstrap: int, alpha: float)

StepM multiple hypothesis testing result.

step_m 

step_m(loss_benchmark: FloatArray, loss_models: FloatArray, alpha: float = 0.05, n_bootstrap: int = 999, avg_block_len: float | None = None, bandwidth: int | None = None, rng: Generator | None = None) -> StepMResult

Romano-Wolf StepM stepdown multiple hypothesis test.

Tests H_k: E[L_bench - L_model_k] <= 0 for k = 1..M. Rejects H_k (model k beats benchmark) for k in result.rejected. Controls FWER <= alpha across all M tests.

Parameters:

Name Type Description Default
loss_benchmark (T,) benchmark loss series (lower = better)
required
loss_models FloatArray
required
alpha float
0.05
n_bootstrap int
999
avg_block_len float | None
None
bandwidth int | None
None
rng Generator | None
None

Returns:

Type Description
StepMResult

.rejected — 0-based indices of models significantly beating benchmark .accepted — the rest
Source code in src/mfe/bootstrap/stepM.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def step_m(
    loss_benchmark: FloatArray,
    loss_models: FloatArray,
    alpha: float = 0.05,
    n_bootstrap: int = 999,
    avg_block_len: float | None = None,
    bandwidth: int | None = None,
    rng: np.random.Generator | None = None,
) -> StepMResult:
    """
    Romano-Wolf StepM stepdown multiple hypothesis test.

    Tests H_k: E[L_bench - L_model_k] <= 0  for k = 1..M.
    Rejects H_k (model k beats benchmark) for k in result.rejected.
    Controls FWER <= alpha across all M tests.

    Parameters
    ----------
    loss_benchmark : (T,) benchmark loss series (lower = better)
    loss_models    : (T, M) alternative model loss series
    alpha          : familywise error rate (default 0.05)
    n_bootstrap    : stationary bootstrap replications
    avg_block_len  : average block length; if None uses T^{1/3}
    bandwidth      : Newey-West bandwidth; if None uses 1.2 * T^{1/3}
    rng            : random generator

    Returns
    -------
    StepMResult
        .rejected  — 0-based indices of models significantly beating benchmark
        .accepted  — the rest
    """
    lb = np.asarray(loss_benchmark, dtype=np.float64)
    lm = np.asarray(loss_models, dtype=np.float64)
    if lm.ndim == 1:
        lm = lm[:, None]

    T, M = lm.shape
    if rng is None:
        rng = np.random.default_rng()
    if avg_block_len is None:
        avg_block_len = max(2.0, T ** (1 / 3))
    if bandwidth is None:
        bandwidth = max(1, int(1.2 * T ** (1 / 3)))

    # Loss differentials: d_{k,t} = L_bench_t - L_model_k_t
    d = lb[:, None] - lm                   # (T, M)
    d_bar = d.mean(axis=0)                 # (M,)
    sigma = _long_run_std(d, bandwidth)    # (M,)
    t_stats = np.sqrt(T) * d_bar / sigma  # (M,)

    # Unadjusted p-values (individual, no FWER control)
    # Use bootstrap max distribution over all M models
    boot = _stationary_bootstrap(d - d_bar[None, :], n_bootstrap, avg_block_len, rng)
    # boot: (n_boot, T, M) — resampled centered loss diffs

    p_raw = np.empty(M, dtype=np.float64)
    boot_max_all = (np.sqrt(T) * boot.mean(axis=1) / sigma[None, :]).max(axis=1)
    for k in range(M):
        boot_k = np.sqrt(T) * boot[:, :, k].mean(axis=1) / sigma[k]
        p_raw[k] = float(np.mean(boot_k >= t_stats[k]))

    # Stepdown procedure
    remaining = list(range(M))
    rejected = []
    p_adjusted = np.ones(M, dtype=np.float64)

    step = 0
    while remaining:
        # Bootstrap max over remaining models
        boot_max = np.empty(n_bootstrap, dtype=np.float64)
        for b in range(n_bootstrap):
            t_boot_remaining = np.sqrt(T) * boot[b, :, :][:, remaining].mean(axis=0) / sigma[remaining]
            boot_max[b] = float(np.max(t_boot_remaining))

        # Critical value at level alpha
        cv = float(np.quantile(boot_max, 1 - alpha))

        # Find the model with max t_stat among remaining
        t_remaining = t_stats[remaining]
        max_idx_in_remaining = int(np.argmax(t_remaining))
        max_k = remaining[max_idx_in_remaining]
        max_t = float(t_stats[max_k])

        if max_t > cv:
            # Reject this model
            p_adjusted[max_k] = float(np.mean(boot_max >= max_t))
            rejected.append(max_k)
            remaining.remove(max_k)
            step += 1
        else:
            # No more rejections possible
            break

    # Adjusted p-values for accepted: use the last step's distribution
    # (conservative: bound by the p-value from the last step)
    for k in remaining:
        p_adjusted[k] = float(np.mean(boot_max_all >= t_stats[k]))

    # Monotonise: stepdown p-values must be non-decreasing when sorted by t_stat descending
    order = np.argsort(-t_stats)
    p_mono = p_adjusted[order].copy()
    for i in range(1, M):
        p_mono[i] = max(p_mono[i], p_mono[i - 1])
    p_adjusted[order] = p_mono

    accepted = [k for k in range(M) if k not in rejected]

    return StepMResult(
        rejected=sorted(rejected),
        accepted=sorted(accepted),
        t_stats=t_stats,
        p_values_raw=p_raw,
        p_values_adjusted=p_adjusted,
        n_models=M,
        n_obs=T,
        n_bootstrap=n_bootstrap,
        alpha=alpha,
    )

wild

Wild bootstrap for realized volatility and related statistics.

Gonçalves, S. & Meddahi, N. (2009): "Bootstrapping Realized Volatility", Econometrica, 77(1), 283-306.

The wild bootstrap resamples by multiplying each squared return by an i.i.d. multiplier w_t drawn from a two-point distribution that matches the first two moments of the standard normal.

This is appropriate for realized volatility statistics because: 1. The squared-return sequence has heterogeneous conditional variance. 2. Block resampling destroys the i.i.d.-ness of squared returns under the null. 3. The wild bootstrap is consistent for RV-based test statistics even in the presence of microstructure noise (with appropriate pre-averaging).

Two-point Rademacher multiplier: w_t = +1 or -1 with prob 1/2. Mammen (1993) multiplier: w_t = -(sqrt(5)-1)/2 or (sqrt(5)+1)/2.

WildBootstrapResult dataclass 

WildBootstrapResult(statistic: float, ci_lower: float, ci_upper: float, ci_level: float, bootstrap_distribution: FloatArray, n_replications: int, multiplier: str)

Result from a wild bootstrap confidence interval computation.

wild_bootstrap_rv 

wild_bootstrap_rv(returns: FloatArray, statistic_fn: Callable[[FloatArray], float] | None = None, n_replications: int = 999, ci_level: float = 0.95, multiplier: str = 'rademacher', rng: Generator | None = None) -> WildBootstrapResult

Wild bootstrap confidence interval for a realized volatility statistic.

The bootstrap DGP is: r_t^* = w_t * r_t

where w_t is i.i.d. from the specified multiplier distribution. The statistic is re-evaluated on {r_t^*}.

Parameters:

Name Type Description Default
returns FloatArray
required
statistic_fn Callable[[FloatArray], float] | None
None
n_replications number of bootstrap replications
999
ci_level float
0.95
multiplier str
'rademacher'
rng Generator | None
None

Returns:

Type Description
WildBootstrapResult
Source code in src/mfe/bootstrap/wild.py 
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def wild_bootstrap_rv(
    returns: FloatArray,
    statistic_fn: Callable[[FloatArray], float] | None = None,
    n_replications: int = 999,
    ci_level: float = 0.95,
    multiplier: str = "rademacher",
    rng: np.random.Generator | None = None,
) -> WildBootstrapResult:
    """
    Wild bootstrap confidence interval for a realized volatility statistic.

    The bootstrap DGP is:
        r_t^* = w_t * r_t

    where w_t is i.i.d. from the specified multiplier distribution.
    The statistic is re-evaluated on {r_t^*}.

    Parameters
    ----------
    returns        : (T,) log-return array
    statistic_fn   : function mapping returns -> float; default is sum(r^2) (RV)
    n_replications : number of bootstrap replications
    ci_level       : confidence level (e.g. 0.95 for 95% CI)
    multiplier     : "rademacher" | "mammen" | "normal"
    rng            : numpy random generator; if None, uses default_rng()

    Returns
    -------
    WildBootstrapResult
    """
    r = np.asarray(returns, dtype=np.float64)
    T = len(r)

    if rng is None:
        rng = np.random.default_rng()

    if statistic_fn is None:
        def statistic_fn(x: FloatArray) -> float:
            return float(np.sum(x ** 2))

    mult_func = _MULTIPLIER_FUNCS.get(multiplier)
    if mult_func is None:
        raise ValueError(f"multiplier must be one of {list(_MULTIPLIER_FUNCS)}, got '{multiplier}'")

    # Point estimate
    stat0 = statistic_fn(r)

    # Bootstrap distribution
    boot_stats = np.empty(n_replications, dtype=np.float64)
    for b in range(n_replications):
        w = mult_func(T, rng)
        r_star = w * r
        boot_stats[b] = statistic_fn(r_star)

    alpha = 1.0 - ci_level
    ci_lo = float(np.percentile(boot_stats, 100 * alpha / 2))
    ci_hi = float(np.percentile(boot_stats, 100 * (1 - alpha / 2)))

    return WildBootstrapResult(
        statistic=stat0,
        ci_lower=ci_lo,
        ci_upper=ci_hi,
        ci_level=ci_level,
        bootstrap_distribution=boot_stats,
        n_replications=n_replications,
        multiplier=multiplier,
    )

wild_bootstrap_test 

wild_bootstrap_test(returns: FloatArray, null_statistic: float, statistic_fn: Callable[[FloatArray], float] | None = None, n_replications: int = 999, multiplier: str = 'rademacher', rng: Generator | None = None) -> tuple[float, float]

Wild bootstrap p-value for a two-sided hypothesis test.

Parameters:

Name Type Description Default
null_statistic the value of the statistic under the null hypothesis
required

Returns:

Type Description
(observed_statistic, bootstrap_p_value)
Source code in src/mfe/bootstrap/wild.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def wild_bootstrap_test(
    returns: FloatArray,
    null_statistic: float,
    statistic_fn: Callable[[FloatArray], float] | None = None,
    n_replications: int = 999,
    multiplier: str = "rademacher",
    rng: np.random.Generator | None = None,
) -> tuple[float, float]:
    """
    Wild bootstrap p-value for a two-sided hypothesis test.

    Parameters
    ----------
    null_statistic : the value of the statistic under the null hypothesis

    Returns
    -------
    (observed_statistic, bootstrap_p_value)
    """
    result = wild_bootstrap_rv(
        returns,
        statistic_fn=statistic_fn,
        n_replications=n_replications,
        multiplier=multiplier,
        rng=rng,
    )
    p_val = float(np.mean(np.abs(result.bootstrap_distribution - null_statistic) >=
                          abs(result.statistic - null_statistic)))
    return result.statistic, p_val