Replication

Alternative paradigms

Non-ADF-family approaches, principally the quantile-based global test and its recursive monitoring extension.

This page is the technical record behind the methods. To learn how to run them, start with Alternative paradigms in the guide.

This file covers methods that address the same problem, detecting explosive or bubble dynamics, from outside the ADF/SADF/GSADF/BSADF recursive-regression family around which exuber is built.

Method	Paper	Status
Quantile-based detection	Pavlidis (2025); Wu, Shi & Wu (2025)	the global test, `QPWY` and `QPSY` monitoring are done (`QPSY` carries a small-sample caveat away from the median); Pavlidis’s `Un`/`QKS` is not implemented
Noncausal / local explosive dynamics	Blasques, Koopman, Mingoli & Telg (2025)	evaluated, not implemented
Spectral fragility	Bhandari (arXiv)	out of scope
Stochastic tree asset pricing	Gourieroux & Jasiak (2025)	out of scope (pricing, not testing)

All papers are from the JTSA 46(5) special issue except Bhandari (arXiv). See references.md.

Quantile-based detection

Status: the global test is quantile_test(), and the recursive monitoring extensions QPWY and QPSY are monitor_quantile().

Sources

Pavlidis, E. G. (2025). Bubbles and crashes: A tale of quantiles. JTSA, 46(5), 884–907.
Wu, R., Shi, S. & Wu, J. (2025). Quantile analysis for financial bubble detection and surveillance. JTSA, 46(5), 908–931. Shi is a coauthor of PSY, so this is a quantile-based test from inside the same lineage.

Idea

These papers characterise explosive behaviour through quantile regression (QR) at chosen points of the conditional distribution, in place of a recursive mean-regression ADF statistic.

The global test of Wu, Shi & Wu (eq. 18) is the QR analogue of the DF $t$ -ratio:

t_T(\tau) = \frac{\hat f(b_\tau)}{\sqrt{\tau(1-\tau)}}\, \bigl(Y_{-1}' P_Z Y_{-1}\bigr)^{1/2} \bigl(\hat\alpha(\tau) - 1\bigr).

Here $\hat\alpha(\tau)$ is the quantile-regression estimator of $y_t$ on an intercept and $y_{t-1}$ at quantile $\tau$ (eq. 13). The term $Y_{-1}' P_Z Y_{-1}$ is the demeaned sum of squares of the lagged level ( $P_Z$ is a demeaning projector and $Z$ a column of ones). The density $\hat f(b_\tau)$ is a kernel estimate of the density of the first-differenced series at its own $\tau$ -th sample quantile (eq. 19, $\hat f(b_\tau) = (Th)^{-1} \sum_t K\bigl((\hat b_\tau - \hat u_t)/h\bigr)$ with $\hat u_t = y_t - y_{t-1}$ ). It studentises the QR coefficient, as the residual variance does for the OLS DF ratio. The paper’s QPWY and QPSY statistics have the sup-scan structure of PWY and PSY, with this QR $t$ -ratio computed in every window.

Critical values (eq. 22–23). The limiting null distribution of $t_T(\tau)$ is

U(\tau) = \sqrt{1 - \delta(\tau)^2}\; z + \delta(\tau)\, Q, \qquad z \sim N(0,1),

where $\delta(\tau)$ is a correlation estimated from the data, $\operatorname{cor}\bigl(\hat u_t,\ \tau - \mathbf 1\{\hat u_t < \hat b_\tau\}\bigr)$ , between the innovation and its quantile-check score. $Q = (\int \bar W^2)^{-1/2} \int \bar W\, dW$ (eq. 23) is the standard demeaned Dickey–Fuller $t$ -distribution. Simulating $Q$ by the random-walk-plus-OLS- $t$ construction of radf_mc_cv() gives values identical, bit for bit, to the adf field of radf() on the same series. The quantiles of $U(\tau)$ therefore need one fresh normal draw and a critical value that the package already simulates.

Optimal quantile. The paper’s criterion (eq. 33) is $\tau^* = \arg\min_\tau \tau(1-\tau)/\hat f(b_\tau)^2$ , a grid search over the same $\hat f(b_\tau)$ .

Implementation: global test

quantile_test() implements the global test of Section 3.1, a single static QR fit at one quantile with no recursion. It adds quantreg (>= 5.9) as an estimation dependency of exuber. tau = "optimal" (the default) runs the grid search of eq. 33 over tau_grid (default seq(0.2, 0.8, by = 0.05), the practical range the paper recommends), and a fixed tau can be passed. quantreg::rq() fits the regression, and the density, the demeaned sum of squares and the critical-value simulation are plain R.

Validation. The empirical size under a pure random-walk null is 5.0% for tau = "optimal" and 3.0% for fixed tau = 0.5 (nominal 5%, 100 replications). Power under an explosive alternative is 100%, as for SADF on the same DGP. The selected $\tau$ lies inside the search grid across replications and does not collapse to a boundary. Replication script: replication/alternative-paradigms/radf_quantile_validation.R.

Implementation: QPWY

$\mathrm{QPWY}_r(\tau) := t_T^{0,r}(\tau)$ is a single recursion. The window start is fixed at 1 and only the end $r$ grows, which is the badf shape of radf(), so it needs $O(T)$ QR fits. QR has no closed form in cumulative sums, so each window needs a genuine fit (Corollaries 1–2, pages 10–11).

Corollary 1 decomposes the limit of the windowed statistic as

U'^{r_1, r_2}(\tau) = \sqrt{1 - \delta(\tau)^2}\; z + \delta(\tau)\, Q_{r_1, r_2},

the same decomposition used for quantile_test(). Corollary 2 identifies $Q_{0,r}$ with the badf sequence of radf() under a simulated null path.

monitor_quantile(data, tau = 0.5, minw, nrep, level, seed) is in exuber/R/monitor_quantile.R. qpwy_stat_path() is the $O(T)$ loop of quantreg::rq() fits and repeats the per-window $t$ -ratio of quantile_test() (eq. 18) over a growing window. The boundary is a single flat value, the quantile across replicates of each simulated path’s own maximum, as in the sadf_cv of radf_mc_cv(). A pointwise marginal quantile at each $r$ would not control the first crossing: it gave a false-alarm rate of 50% against a nominal 5%.

The limit in Theorem 1 is $\int \tilde W\, dB_\psi / \sqrt{\int \tilde W^2}$ , where $B_\psi$ is a Brownian motion with correlation $\delta$ to $W$ . Writing $B_\psi = \delta W + \sqrt{1 - \delta^2}\, V$ with $V$ independent of $W$ gives

\delta\, Q_{r_1, r_2} + \sqrt{1 - \delta^2}\, Z_{r_1, r_2}, \qquad Z = \frac{\int \tilde W\, dV}{\sqrt{\int \tilde W^2}} .

For a single window $Z$ is exactly $N(0,1)$ , which is all that quantile_test() needs. It varies across windows, however, and a monitoring boundary is a quantile of path suprema, so a single constant draw of $z$ understates it. At $n = 200$ (4000 replications, Gaussian innovations, $\tau = 0.5$ ) the single- $z$ 95% boundary and the correct one are:

$\delta$	single- $z$ boundary	correct boundary	true size of single- $z$
0.8	1.534	1.665	0.066
0.5	1.657	2.060	0.108
0.2	1.704	2.405	0.196

The distortion is largest where QPWY is meant to help, with heavy tails and non-central quantiles, because $\delta$ is small there. quantile_boundary_sim() simulates $Q$ and $Z$ jointly for every window from prefix sums of $(e_t, v_t)$ , at $O(1)$ per window, takes each path’s supremum for the data-estimated $\delta$ and matches a per-window brute force to $2 \times 10^{-15}$ .

Validation. The false-alarm rate under $H_0$ at a nominal 5% ( $n = 150$ , 200 replications) is 0.035 for Gaussian innovations and for $t_3$ at $\tau = 0.5$ . For $t_3$ at $\tau = 0.2$ , $0.8$ and $0.9$ it is 0.075, 0.085 and 0.125. The last is well above nominal, in line with the paper’s advice to avoid extreme quantiles in small samples (its Table II reports 0.07 for the global test at $\tau = 0.9$ under $t_3$ ). Power on a post-training explosive DGP (60 replications) is 50.0%, against 55.0% for SADF on the same DGP. QPWY gives up a little power for robustness to non-Gaussian innovations, which is the paper’s motivation. test-qpwy.R has 7 tests, including a check that the supremum-calibrated boundary stochastically dominates any single-column marginal quantile.

Implementation: QPSY

$\mathrm{QPSY}_r(\tau, r_0) = \sup_{r_1} t^{r_1, r}(\tau)$ (eq. 26) needs $O(T^2)$ QR fits, about $T^2/2$ per series. In R with rq.fit(method = "br") that takes 2.7 s at $n = 100$ and 12 s at $n = 200$ . The boundary needs none of them. The $\sup_{r_1} [\delta\, Q_{r_1, r} + \sqrt{1 - \delta^2}\, Z_{r_1, r}]$ of Corollary 2 comes from the same $Q$ / $Z$ simulation as QPWY, run over the full $(r_1, r)$ grid instead of the $r_1 = 0$ row.

monitor_quantile(..., type = "qpsy") uses windows $[r_1, r]$ with at least minw regression observations (the same first-window floor as QPWY) and the same flat, supremum-calibrated boundary.

Validation. qpsy_stat_path() matches a quantreg::rq() brute force over every window exactly. The grid simulation matches a per-window brute force to $2 \times 10^{-15}$ , and the QPSY suprema dominate those of QPWY on shared draws, as they must. At a nominal 5% ( $n = 100$ , 80 to 100 replications), Gaussian innovations give size 0.040 at $\tau = 0.5$ and 0.350 at $\tau = 0.9$ . With $t_3$ innovations the size is 0.037 at $\tau = 0.5$ , 0.212 at $\tau = 0.8$ and 0.440 at $\tau = 0.9$ . For power ( $n = 100$ , explosive from $t = 71$ , $\rho = 1.04$ , Gaussian innovations, 60 replications), QPWY gives 0.400, QPSY 0.433 and SADF 0.483. The OLS test is ahead with Gaussian errors, as in Table V of the paper.

Caveat. The asymptotic boundary is exact and well sized at the median, and does not hold in the small early windows (about 20 observations) away from the median. A double supremum over thousands of such windows amplifies the finite-sample error, which the single supremum of QPWY mostly avoids. The paper uses bootstrap critical values for monitoring (Algorithm 1, used for Table V) and advises against extreme quantiles. With tau away from 0.5, type = "qpsy" emits a caveat as a message and as attr(x, "caveat"), and ?monitor_quantile gives the numbers. pyexuber has monitor_quantile(type="qpsy") with a UserWarning for the caveat.

The bootstrap of Algorithm 1 is not implemented. It resamples the centred $\Delta y$ and recomputes the whole statistic path for every replicate, so each replicate repeats the full $O(T^2)$ QR sweep. That takes about 9 minutes per series at $n = 100$ with 199 replicates.

Tests are in test-qpwy.R. Replication script: replication/alternative-paradigms/radf_qpwy_validation.R. The sizes are Monte Carlo estimates from 80 to 200 replications, good to about 2 or 3 points.

Pavlidis’s quantile-autoregressive tests: not implemented

Pavlidis characterises bubbles through unit-root quantile-autoregressive models in which the largest autoregressive root may vary by quantile (below 1 at low quantiles and crashes, above 1 at high quantiles and expansions). Pages 6–7 (eq. 5–11) give the same ADF regression form as radf(), fitted by quantile regression at chosen $\tau$ . The statistics are $U_n(\tau) = n\,(\hat\alpha_1(\tau) - 1)$ (coefficient-based) and $\mathrm{QKS} = \sup_{\tau \in \mathcal T} U_n(\tau)$ (eq. 11). Critical values come from a residual or sieve bootstrap (page 7, steps 1–5) that is close to the Pedersen–Schütte bootstrap of radf_sb_(): fit an AR( $q$ ) to $\operatorname{diff}(y)$ under $H_0$ , resample the centred residuals and rebuild the series by cumulation. Table 2 of the paper gives empirical sizes (N(0,1), $t_3$ and $t_2$ errors, $n = 100, 200, 300, 400$ ).

The tests are not implemented because the bootstrap does not reproduce the published sizes. In a prototype, the size of $U_n(\tau = 0.5)$ at $n = 100$ with N(0,1) errors matched Table 2 (0.050 against 0.053), but $U_n$ at higher $\tau$ and $\mathrm{QKS}$ were oversized (0.075–0.100 against 0.052–0.063 published). Comparing the oracle null distribution of $U_n$ and $\mathrm{QKS}$ (1000 i.i.d. random walks, no bootstrap) with the critical values implied by the bootstrap on a single series, the bootstrap critical value stays substantially below the oracle at $\tau = 0.5, 0.8, 0.9$ even at 1999 bootstrap replications, so the gap is not a matter of replication count. Only $\tau = 0.95$ and $\mathrm{QKS}$ came close to the oracle at large nboot. The statistic $U_n(\tau) = n(\hat\alpha_1(\tau) - 1)$ is very sensitive to small differences in $\hat\alpha_1$ , because the coefficient is near 1 and $n \approx 99$ amplifies third-decimal differences by about 100. A next step would be to compare the bootstrap and oracle distributions of $\hat\alpha_1$ itself at several $\tau$ , or to implement $\mathrm{QKS}$ alone, the headline statistic of the paper, which calibrated well.

Noncausal / local explosive dynamics

Status: evaluated, not implemented.

Source

Blasques, F., Koopman, S. J., Mingoli, G. & Telg, S. (2025). A Novel Test for the Presence of Local Explosive Dynamics. JTSA, 46(5), 966–980, doi:10.1111/jtsa.70001.

Idea

The test is built for mixed causal-noncausal autoregressive processes, a model class in which part of the dynamics depends on future shocks (anticipative or noncausal terms) as well as past ones. The premise is that bubbles come from an extreme shock acting through the forward-looking component of the model, and not from a recursively estimated explosive AR root on the past alone. The distribution of the test statistic is derived analytically or approximated numerically, depending on the assumed error distribution. The application is a monthly oil price index, framed partly as a Value-at-Risk-style risk-assessment tool.

Fit with exuber

Mixed causal-noncausal AR models need their own estimation. There is no closed-form OLS or QR reduction. Noncausal components are typically estimated by approximate or simulated maximum likelihood under a specified non-Gaussian error distribution, because noncausal processes are identifiable only with non-Gaussian innovations. None of the recursive least squares of exubercore applies, nor does any transform-and-reuse approach of the kind used for STADF, the sign-based test or PDC. It would be a separate statistical framework with a new estimation dependency, and no mainstream R package for noncausal AR fitting is available.

Spectral fragility (out of scope)

Bhandari, A. Rational Bubbles at the Spectral Edge: An Operator-Spectral Theory of Fragility, Identification and Finite-Sample Certification. arXiv:2607.03933.

The paper detects factor and co-movement spectral fragility. It identifies market fragility through the strength of a dominant factor extracted from cross-sectional co-movement (fewer independent factors during crises than in calm periods). It is not a right-tailed unit-root test on a single series, and it detects fragility contemporaneously and not predictively. It is not built on ADF machinery, and it addresses market-wide fragility and not the explosiveness of a specific series, so it falls outside the scope of exuber.

Stochastic tree asset pricing (out of scope)

Gourieroux, C. & Jasiak, J. (2025). A Stochastic Tree for Bubble Asset Modelling and Pricing. JTSA, 46(5), 932–944.

The paper presents a stochastic-tree representation for modelling, forecasting and pricing bubbles, with closed-form option-pricing formulas. It is an asset-pricing model and not a test for the presence of a bubble, and exuber is concerned with testing.

Replication scripts

4 standalone scripts support the claims above, and none depends on another. The R scripts run from the exuber/ package root. The 2 Python scripts check that pyexuber reproduces the same numbers.

radf_qpwy_validation.R 130 lines

# Validation of monitor_quantile(), the QPWY and QPSY recursive quantile
# monitoring of Wu, Shi & Wu (2025). See docs/alternative-paradigms.md,
# "Quantile-based detection".
#
# The size checks show why the boundary is the quantile of each path's
# supremum and not a pointwise quantile, and why the independent Brownian
# component Z of the limit is simulated as a process over windows and not as
# one draw per replicate.
#
# Run from the exuber-project/ root. Sections 4-5 take about 30 minutes.

devtools::load_all("exuber", quiet = TRUE)

cat("=== 1. Point statistics vs quantreg::rq() brute force ===\n")
set.seed(3)
y <- cumsum(rnorm(60))
rq_stat <- function(yy, tau) {
  m <- length(yy)
  ylag <- yy[1:(m - 1)]
  yresp <- yy[2:m]
  a <- unname(coef(quantreg::rq(yresp ~ ylag, tau = tau))["ylag"])
  f <- exuber:::quantile_check_density(yresp - ylag, tau)$f_hat
  (f / sqrt(tau * (1 - tau))) * sqrt(sum((ylag - mean(ylag))^2)) * (a - 1)
}
d1 <- abs(exuber:::qpwy_stat_path(y, 0.5, 60) - rq_stat(y, 0.5))
minw <- 15
r_idx <- (minw + 1):60
qpsy <- exuber:::qpsy_stat_path(y, 0.7, r_idx, minw)
brute <- vapply(r_idx, function(r) max(vapply(1:(r - minw), function(r1) rq_stat(y[r1:r], 0.7), 0)), 0)
cat(sprintf("QPWY |diff| = %.2e   QPSY max|diff| = %.2e\n", d1, max(abs(qpsy - brute))))

cat("\n=== 2. quantile_boundary_sim vs per-window brute force ===\n")
n <- 30
minw <- 8
delta <- c(0.3, 0.9)
worst <- 0
for (type in c("qpwy", "qpsy")) {
  sim <- exuber:::quantile_boundary_sim(n, minw, 3, delta, type = type, seed = 11)
  set.seed(11)
  for (i in 1:3) {
    e <- rnorm(n - 1)
    v <- rnorm(n - 1)
    x <- c(0, cumsum(e))[1:(n - 1)]
    U <- NULL
    for (hi in minw:(n - 1)) {
      for (lo in if (type == "qpwy") 0 else 0:(hi - minw)) {
        k <- (lo + 1):hi
        xb <- x[k] - mean(x[k])
        U <- rbind(U, (delta * sum(xb * e[k]) + sqrt(1 - delta^2) * sum(xb * v[k])) / sqrt(sum(xb^2)))
      }
    }
    worst <- max(worst, abs(sim[i, ] - apply(U, 2, max)))
  }
}
cat(sprintf("max|diff| over qpwy/qpsy, 3 reps, 2 deltas = %.2e\n", worst))

cat("\n=== 3. One z per replicate vs Z as a process over windows ===\n")
# QPWY slice (lo = 0), n = 200, psy_minw, 4000 reps: 95% boundary under the
# old single-z construction, and that boundary's true size under the
# correct limiting process
n <- 200
minw <- psy_minw(n)
nrep <- 4000
set.seed(1)
hi <- minw:(n - 1)
Qm <- Zm <- matrix(NA_real_, nrep, length(hi))
for (i in seq_len(nrep)) {
  e <- rnorm(n - 1)
  v <- rnorm(n - 1)
  x <- c(0, cumsum(e))[1:(n - 1)]
  cs <- function(w) cumsum(w)[hi]
  sxx <- cs(x^2) - cs(x)^2 / hi
  Qm[i, ] <- (cs(x * e) - cs(x) * cs(e) / hi) / sqrt(sxx)
  Zm[i, ] <- (cs(x * v) - cs(x) * cs(v) / hi) / sqrt(sxx)
}
z1 <- rnorm(nrep)
for (d in c(0.8, 0.5, 0.2)) {
  right <- apply(d * Qm + sqrt(1 - d^2) * Zm, 1, max)
  single <- apply(d * Qm + sqrt(1 - d^2) * z1, 1, max)
  b_single <- quantile(single, 0.95)
  cat(sprintf(
    "delta = %.1f: single-z boundary %.3f, correct %.3f, true size of single-z boundary %.3f\n",
    d, b_single, quantile(right, 0.95), mean(right > b_single)
  ))
}

cat("\n=== 4. Empirical false-alarm rate under H0 (nominal 5%) ===\n")
fa_rate <- function(type, n, reps, innov, tau) {
  mean(vapply(seq_len(reps), function(i) {
    set.seed(5000 + i)
    yy <- cumsum(innov(n))
    !is.na(suppressMessages(monitor_quantile(yy, tau = tau, nrep = 300, seed = i, type = type))$alarm)
  }, logical(1)))
}
gauss <- function(n) rnorm(n)
t3 <- function(n) rt(n, df = 3)
for (cfg in list(
  list("qpwy", 150, 200, "gaussian", gauss, 0.5),
  list("qpwy", 150, 200, "t3", t3, 0.5),
  list("qpwy", 150, 200, "t3", t3, 0.9),
  list("qpwy", 150, 200, "t3", t3, 0.8),
  list("qpwy", 150, 200, "t3", t3, 0.2),
  list("qpsy", 100, 100, "gaussian", gauss, 0.5),
  list("qpsy", 100, 100, "t3", t3, 0.9),
  list("qpsy", 100, 80, "t3", t3, 0.5),
  list("qpsy", 100, 80, "t3", t3, 0.8),
  list("qpsy", 100, 80, "gaussian", gauss, 0.9)
)) {
  cat(sprintf(
    "%s n=%d reps=%d %s tau=%.1f: %.3f\n", cfg[[1]], cfg[[2]], cfg[[3]], cfg[[4]], cfg[[6]],
    fa_rate(cfg[[1]], cfg[[2]], cfg[[3]], cfg[[5]], cfg[[6]])
  ))
}

cat("\n=== 5. Detection power, explosive from t = 71 of n = 100, vs SADF ===\n")
reps <- 60
det <- c(qpwy = 0, qpsy = 0, sadf = 0)
cv <- radf_mc_cv(100, nrep = 2000, seed = 1)
for (i in seq_len(reps)) {
  set.seed(3000 + i)
  normal_part <- cumsum(rnorm(70))
  yy <- c(normal_part, normal_part[70] * 1.04^(1:30) + cumsum(rnorm(30)))
  det["qpwy"] <- det["qpwy"] + !is.na(monitor_quantile(yy, nrep = 300, seed = i)$alarm)
  det["qpsy"] <- det["qpsy"] + !is.na(monitor_quantile(yy, nrep = 300, seed = i, type = "qpsy")$alarm)
  det["sadf"] <- det["sadf"] + (radf(yy)$sadf > cv$sadf_cv[2])
}
print(round(det / reps, 3))

cat("\ndone\n")

radf_qpwy_validation.py 125 lines

"""Python cross-check of exuber's monitor_quantile() (the QPWY recursive
quantile monitoring of Wu, Shi & Wu 2025), mirroring
radf_qpwy_validation.R's formula-exact check of the per-window QR
t-ratio and the boundary-vs-marginal-quantile structural check (see
docs/alternative-paradigms.md, "Implementation: QPWY").

Run standalone: uv run --project pyexuber python
docs/replication/alternative-paradigms/radf_qpwy_validation.py

_qpwy_stat_path() needs no RNG and no radf()/C++ extension, so its check
runs fully offline and matches R bit-for-bit (within the IRLS-vs-simplex
QR-solver tolerance documented in monitor.py's module docstring). The
boundary simulation needs no radf() either (Q and Z come from prefix sums),
so sections 4-5 check it directly: the independent-BM term Z must be a
process over windows and not a single draw per path, and the QPSY grid
contains the QPWY path.

Reference values: a direct R run from the exuber-project/ root:
    Rscript -e 'Sys.setenv(NOT_CRAN="true");
    devtools::load_all("exuber", quiet=TRUE); set.seed(42);
    y <- cumsum(rnorm(80)); minw <- 15;
    r_idx <- (minw + 1L):length(y);
    stat <- exuber:::qpwy_stat_path(y, 0.5, r_idx);
    dput(round(tail(stat, 5), 10)); dput(length(stat))'
"""

import sys
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parents[3] / "pyexuber" / "src"))

from exuber.monitor_quantile import (  # noqa: E402
    _qpsy_stat_path,
    _qpwy_stat_path,
    _quantile_boundary_sim,
)
from exuber.quantile_test import quantile_test  # noqa: E402

# set.seed(42); y <- cumsum(rnorm(80))
Y42 = np.array(
    [
        1.3709584471, 0.8062602758, 1.1693886871, 1.802251292, 2.2065196152,
        2.1003950991, 3.6119170965, 3.5172580581, 5.535681772, 5.4729676729,
        6.7778373272, 9.0644827199, 7.6756220188, 7.3968332519, 7.2635119156,
        7.8994623136, 7.6152093922, 4.9587539713, 2.5182870427, 3.8384003885,
        3.5317617944, 1.7504533604, 1.5785360046, 2.7932107038, 4.6884041651,
        4.2579350335, 4.0006656507, 2.2375025655, 2.6975999203, 2.0576050444,
        2.5130551676, 3.2178925048, 4.2529960268, 3.6440696514, 4.1490247747,
        2.4320160956, 1.6475570873, 0.7966494931, -1.6175581569, -1.58143555,
        -1.3754369498, -1.7364942483, -0.9783310126, -1.7050358397, -3.0733168841,
        -2.6404988582, -3.4518920344, -2.0077907727, -2.4392369753, -1.7835890919,
        -1.4616638267, -2.2455027676, -0.6697752478, -0.0268759421, 0.0628847045,
        0.3394354518, 1.0187242679, 1.1085571544, -1.8845329287, -1.5996499752,
        -1.9668846179, -1.7816540531, -1.1998303257, 0.1999065016, -0.5273855579,
        0.7751570742, 1.1110051939, 2.1495112926, 3.0702398609, 3.7911180238,
        2.7479990852, 2.6578126986, 3.2813308606, 2.3278075028, 1.7849786883,
        2.3659751859, 3.1341539238, 3.5979215123, 2.7121452149, 1.6123643163,
    ]
)


def main() -> None:
    print("=== 1. Formula-exact check: qpwy_stat_path -- cross-check vs R ===")
    minw = 15
    r_idx = np.arange(minw + 1, len(Y42) + 1)
    stat = _qpwy_stat_path(Y42, 0.5, r_idx)
    assert len(stat) == 65
    expected_tail = np.array(
        [-1.128241465, -1.2210267401, -1.2384527008, -1.2255226955, -1.1087521367]
    )
    np.testing.assert_allclose(stat[-5:], expected_tail, atol=1e-4)
    print(f"len={len(stat)}, tail={stat[-5:]}")

    print("\n=== 2. Structural check: Q_{0,r} at r=n equals quantile_test()'s tstat ===")
    # Corollary 2's own claim: QPWY_n(tau) (the full-window recursion
    # endpoint) is exactly quantile_test()'s own single-shot tstat.
    qt = quantile_test(Y42, tau=0.5, nrep=10, sig_lvl=95, seed=1)
    assert abs(stat[-1] - qt.tstat[0]) < 1e-6
    print(f"stat[-1]={stat[-1]}, quantile_test tstat={qt.tstat[0]}: match")

    print(
        "\n=== 3. Boundary-vs-marginal-quantile sanity ==="
    )
    # A per-r MARGINAL quantile of simulated paths badly inflates the
    # false-alarm rate relative to a SUPREMUM-calibrated one (see
    # docs/alternative-paradigms.md). Demonstrated here directly on
    # simulated standard-normal paths.
    rng = np.random.default_rng(0)
    nrep, n_mon = 500, 40
    paths = rng.normal(size=(nrep, n_mon))  # stand-in for the Q_{0,r} paths
    marginal_boundary = np.quantile(paths, 0.95, axis=0)  # WRONG (per-r)
    sup_boundary = np.quantile(paths.max(axis=1), 0.95)  # RIGHT (supremum)

    fpr_marginal = np.mean(np.any(paths > marginal_boundary, axis=1))
    fpr_sup = np.mean(paths.max(axis=1) > sup_boundary)
    print(f"per-r marginal boundary false-alarm rate: {fpr_marginal:.3f} (badly inflated)")
    print(f"supremum-calibrated boundary false-alarm rate: {fpr_sup:.3f} (~nominal 0.05)")
    assert fpr_marginal > 3 * fpr_sup
    assert abs(fpr_sup - 0.05) < 0.03

    print("\n=== 4. Z is a process over windows, not one z per path ===")
    # delta = 0 leaves only Z: one shared z per path would make sup_r Z
    # exactly N(0,1), 95% quantile 1.645
    sim = _quantile_boundary_sim(150, 20, 400, np.array([0.0]), False, np.random.default_rng(3))
    q95 = np.quantile(sim[:, 0], 0.95)
    print(f"95% quantile of sup_r Z_r: {q95:.3f} (single-z construction: 1.645)")
    assert q95 > 2

    print("\n=== 5. QPSY: its grid contains QPWY's path ===")
    sy = _qpsy_stat_path(Y42, 0.5, r_idx, minw)
    assert np.all(sy >= stat - 1e-12) and abs(sy[0] - stat[0]) < 1e-12
    d = np.array([0.2, 0.8])
    wy = _quantile_boundary_sim(60, 12, 20, d, False, np.random.default_rng(5))
    sb = _quantile_boundary_sim(60, 12, 20, d, True, np.random.default_rng(5))
    assert np.all(sb >= wy - 1e-12)
    print(f"QPSY path >= QPWY path everywhere; max QPSY = {sy.max():.4f}")

    print("\nAll checks passed.")


if __name__ == "__main__":
    main()

radf_quantile_validation.R 67 lines

# Validation script for quantile_test() (Wu, Shi & Wu 2025, "Quantile
# analysis for financial bubble detection and surveillance", the
# "global test" of their Section 3.1). See docs/
# alternative-paradigms.md, "Quantile-based detection", for the full
# write-up. Run from the exuber-project/ root (or adjust the
# devtools::load_all() path below).

Sys.setenv(NOT_CRAN = "true")
options(exuber.parallel = FALSE, exuber.show_progress = FALSE)
devtools::load_all("exuber", quiet = TRUE)

cat("=== 1. Structural check: the critical value's Q functional is\n")
cat("exactly radf()'s own single-shot adf t-statistic, not a new\n")
cat("computation (bit-for-bit) ===\n\n")
set.seed(1)
y <- cumsum(rnorm(100))
full <- radf(y, minw = 90)
q_manual <- exuber:::quantile_adf_tstat(y)
cat("radf()$adf:", full$adf, " manual Q:", q_manual, " match:",
    isTRUE(all.equal(unname(full$adf), q_manual, tolerance = 1e-8)), "\n\n")

cat("=== 2. Empirical size under H0 (pure random walk), tau='optimal' ===\n")
run_null <- function(seed) {
  set.seed(seed)
  y <- cumsum(rnorm(100))
  quantile_test(y, nrep = 200, seed = 1)$detected[["series1"]]
}
rate_h0 <- mean(sapply(1:100, run_null))
cat(sprintf("False-detection rate (100 reps, nominal 5%%): %.3f\n\n", rate_h0))

cat("=== 3. Empirical size under H0, fixed tau=0.5 ===\n")
run_null_fixed <- function(seed) {
  set.seed(seed)
  y <- cumsum(rnorm(100))
  quantile_test(y, tau = 0.5, nrep = 200, seed = 1)$detected[["series1"]]
}
rate_h0_fixed <- mean(sapply(1:100, run_null_fixed))
cat(sprintf("False-detection rate at tau=0.5 (nominal 5%%): %.3f\n\n", rate_h0_fixed))

cat("=== 4. Power under a genuine explosive alternative, vs standard SADF ===\n")
run_power <- function(seed) {
  set.seed(seed)
  n1 <- 60
  y <- 100 * 1.03^(1:n1) + cumsum(rnorm(n1, sd = 1))
  quantile_test(y, nrep = 200, seed = 1)$detected[["series1"]]
}
run_power_sadf <- function(seed) {
  set.seed(seed)
  n1 <- 60
  y <- 100 * 1.03^(1:n1) + cumsum(rnorm(n1, sd = 1))
  full <- radf(y, minw = 20)
  cv <- radf_mc_cv(n = n1, minw = 20, nrep = 200, seed = 1)
  full$sadf > cv$sadf_cv["95%"]
}
cat(sprintf("Detection rate, quantile_test: %.3f\n", mean(sapply(1:60, run_power))))
cat(sprintf("Detection rate, standard SADF (same DGP, rough cross-check): %.3f\n\n",
            mean(sapply(1:60, run_power_sadf))))

cat("=== 5. Optimal-tau selection lands sensibly within the search grid ===\n")
taus <- sapply(1:20, function(s) {
  set.seed(s)
  y <- cumsum(rnorm(120))
  quantile_test(y, nrep = 100, seed = 1)$tau[["series1"]]
})
cat("selected taus:", paste(taus, collapse = ", "), "\n")
cat("all within [0.2, 0.8]:", all(taus >= 0.2 & taus <= 0.8), "\n")

radf_quantile_validation.py 114 lines

"""Python cross-check of exuber's quantile_test() (Wu, Shi & Wu 2025's
quantile-based global test), mirroring radf_quantile_validation.R's
structural check that Q matches radf()'s own adf statistic, plus a
direct bit-for-bit cross-check of the deterministic point statistics
(tstat/tau/delta) against R.

Run standalone: uv run --project pyexuber python
docs/replication/alternative-paradigms/radf_quantile_validation.py

Reference values: a direct R run from the exuber-project/ root,
    Rscript -e 'Sys.setenv(NOT_CRAN="true");
    devtools::load_all("exuber", quiet=TRUE); set.seed(42);
    y <- cumsum(rnorm(80));
    dput(exuber:::quantile_adf_tstat(y));
    qcd <- exuber:::quantile_check_density(diff(y), 0.5);
    dput(qcd$b_tau); dput(qcd$f_hat);
    qt <- quantile_test(y, tau=0.5, nrep=50, sig_lvl=95, seed=7);
    dput(unname(qt$tstat)); dput(round(unname(qt$delta), 10));
    qt_opt <- quantile_test(y, tau="optimal", nrep=50, sig_lvl=95, seed=7);
    dput(unname(qt_opt$tau)); dput(unname(qt_opt$tstat))'

Note on tolerance: pyexuber's QR fit uses an IRLS solver (see monitor.py's
module docstring) rather than R's quantreg::rq() simplex/interior-point
solver, so tstat matches to ~1e-6, not to machine precision the way the
fully closed-form pieces (badf/cusum/lbi) do elsewhere in this family.
"""

import sys
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parents[3] / "pyexuber" / "src"))

from exuber.quantile_test import (  # noqa: E402
    _quantile_adf_tstat,
    _quantile_check_density,
    quantile_test,
)

# set.seed(42); y <- cumsum(rnorm(80)) -- same series as
# docs/replication/monitoring/radf_monitor_kurozumi_boundary_validation.py.
Y42 = np.array(
    [
        1.3709584471, 0.8062602758, 1.1693886871, 1.802251292, 2.2065196152,
        2.1003950991, 3.6119170965, 3.5172580581, 5.535681772, 5.4729676729,
        6.7778373272, 9.0644827199, 7.6756220188, 7.3968332519, 7.2635119156,
        7.8994623136, 7.6152093922, 4.9587539713, 2.5182870427, 3.8384003885,
        3.5317617944, 1.7504533604, 1.5785360046, 2.7932107038, 4.6884041651,
        4.2579350335, 4.0006656507, 2.2375025655, 2.6975999203, 2.0576050444,
        2.5130551676, 3.2178925048, 4.2529960268, 3.6440696514, 4.1490247747,
        2.4320160956, 1.6475570873, 0.7966494931, -1.6175581569, -1.58143555,
        -1.3754369498, -1.7364942483, -0.9783310126, -1.7050358397, -3.0733168841,
        -2.6404988582, -3.4518920344, -2.0077907727, -2.4392369753, -1.7835890919,
        -1.4616638267, -2.2455027676, -0.6697752478, -0.0268759421, 0.0628847045,
        0.3394354518, 1.0187242679, 1.1085571544, -1.8845329287, -1.5996499752,
        -1.9668846179, -1.7816540531, -1.1998303257, 0.1999065016, -0.5273855579,
        0.7751570742, 1.1110051939, 2.1495112926, 3.0702398609, 3.7911180238,
        2.7479990852, 2.6578126986, 3.2813308606, 2.3278075028, 1.7849786883,
        2.3659751859, 3.1341539238, 3.5979215123, 2.7121452149, 1.6123643163,
    ]
)


def main() -> None:
    print("=== 1. quantile_adf_tstat / quantile_check_density -- cross-check vs R ===")
    tstat = _quantile_adf_tstat(Y42)
    np.testing.assert_allclose(tstat, -1.66093811600216, atol=1e-6)
    print(f"quantile_adf_tstat(y) = {tstat} (expect -1.66093811600216)")

    b_tau, f_hat = _quantile_check_density(np.diff(Y42), 0.5)
    np.testing.assert_allclose(b_tau, 0.0898328865790818, atol=1e-8)
    np.testing.assert_allclose(f_hat, 0.369113058891683, atol=1e-6)
    print(f"b_tau={b_tau}, f_hat={f_hat}")

    print("\n=== 2. quantile_test(tau=0.5) -- cross-check vs R ===")
    qt = quantile_test(Y42, tau=0.5, nrep=50, sig_lvl=95, seed=7)
    np.testing.assert_allclose(qt.tstat[0], -1.1087521367351, atol=1e-4)
    np.testing.assert_allclose(qt.delta[0], 0.7832341878, atol=1e-6)
    print(f"tstat={qt.tstat[0]} (expect -1.1087521367351), delta={qt.delta[0]}")

    print("\n=== 3. quantile_test(tau='optimal') -- cross-check vs R ===")
    qt_opt = quantile_test(Y42, tau="optimal", nrep=50, sig_lvl=95, seed=7)
    assert qt_opt.tau[0] == 0.8
    np.testing.assert_allclose(qt_opt.tstat[0], -0.440485749709493, atol=1e-4)
    print(f"tau={qt_opt.tau[0]} (expect 0.8), tstat={qt_opt.tstat[0]}")

    print("\n=== 4. Empirical size under H0 (100 reps, tau=0.5) ===")
    rng = np.random.default_rng(11)
    rejections = []
    for _ in range(100):
        y_null = np.cumsum(rng.normal(size=100))
        res = quantile_test(y_null, tau=0.5, nrep=200, sig_lvl=95, seed=1)
        rejections.append(bool(res.detected[0]))
    fpr = np.mean(rejections)
    print(f"empirical size: {fpr:.3f} (nominal 0.05)")

    print("\n=== 5. Detection power under a genuine explosive alternative ===")
    rng = np.random.default_rng(12)
    detect = []
    for _ in range(30):
        n_bubble = 100
        e = rng.normal(size=n_bubble - 1)
        y_bubble = np.concatenate(([0.0], np.cumsum(1.03 ** np.arange(n_bubble - 1) + e)))
        res = quantile_test(y_bubble, tau=0.5, nrep=200, sig_lvl=95, seed=1)
        detect.append(bool(res.detected[0]))
    print(f"detection rate: {np.mean(detect):.3f}")

    print("\nAll checks passed.")


if __name__ == "__main__":
    main()