Home

HTDV — Hypothesis Testing for Dependent Variables with Unbalanced Data

A unified R framework for inference on dependent and unbalanced data under strong-mixing conditions. Combines hierarchical Bayesian estimation via Hamiltonian Monte Carlo with frequentist (HAR, block bootstrap) and distribution-free (adaptive conformal) robustness anchors, with explicit finite-sample constants for the metric equivalence between three law-of-large-numbers regimes.

What HTDV does, in one paragraph
When to use HTDV (and when not)
Quick start
Mathematical foundations
The three convergence regimes: TAC, WSC, MPC
The metric equivalence theorem
The Bayesian core
Finite-sample frequentist anchors
Heavy-tail and Rényi refinements
Function reference
Validation evidence
Tutorial walkthroughs
Frequently asked questions
References

1. What HTDV does, in one paragraph

HTDV answers a single inferential question — do these dependent and possibly unequally-sized samples come from a common population? — under the worst combination of nuisance conditions an applied statistician encounters: temporal dependence with polynomial decay, heavy tails near the boundary of finite second moments, severe imbalance across groups, and finite samples that put asymptotic theory under stress. It does so by running three independent inferential layers in parallel — a hierarchical Bayesian fit via Hamiltonian Monte Carlo, a heteroskedasticity-and- autocorrelation-robust (HAR) Wald test with fixed-bandwidth critical values, and a stationary block bootstrap with automatic block length — and exposing the disagreement between them as a calibration signal. Where all three layers concur, the inferential conclusion is robust. Where they disagree, the framework points the user to the layer whose assumptions are most defensible for the data at hand. The choice of three layers is not a stylistic preference: a pre-registered factorial Monte Carlo (1,024 cells × 500 replications) shipped with the package shows that no single layer maintains nominal calibration across the full design, and an external benchmark suite reproduces three published references on public-source datasets to certify the framework against the literature.

2. When to use HTDV (and when not)

Use HTDV when

Your data are time-dependent, spatially dependent, or otherwise exhibit serial / cross-sectional correlation that breaks the independence assumption of textbook tests.
The samples you compare are of unequal size (panels with different follow-up windows; sectors with different aggregation levels; groups with different sampling intensity).
You suspect heavy tails but cannot rule them out a priori — i.e., you don't know whether the second moment is finite enough for Bernstein-style concentration.
You need a hypothesis test or interval estimate for a population mean, group-mean difference, or inflation/yield/return parameter that requires inference, not point estimation.
You want inferential answers that survive review — the framework ships its own validation evidence and external benchmarks.

Do not use HTDV (or use with caution) when

Your data are independent and identically distributed with finite second moments. Classical tools (Welch t-test, ordinary OLS) are simpler, equivalent, and computationally cheaper.
Your dependence is long-memory ($\alpha$-mixing decay slower than $t^{-1}$). The metric equivalence theorem of Section 6 requires $\gamma > 1$. Long-memory inference needs different tools (FELT, FARIMA).
Your data have unit roots or cointegration. HTDV's AR(1) Bayesian core can fit near-unit-root data but its conclusions are honest about the resulting uncertainty inflation; for explicit unit-root testing, use ADF / Phillips-Perron / KPSS.
Your data have structural breaks. The TAC and WSC routes are fragile under regime change; the framework's decision table flags this. If you know the break date, segment the sample first.
Your goal is prediction, not inference. HTDV is built for hypothesis testing and parameter estimation under uncertainty, not for forecasting accuracy.

3. Quick start

Installation

# Stable from GitHub
remotes::install_github("IsadoreNabi/HTDV")

# With the validation vignette
remotes::install_github("IsadoreNabi/HTDV", build_vignettes = TRUE)

rstan is required (compile-on-first-use Stan models). Optional backends: bridgesampling (Bayes factors), loo (WAIC / PSIS-LOO), posterior (draws utilities), bayesplot (visualization), readxl (empirical benchmarks vignette).

A minimal one-sample example

library(HTDV)

# Generate a dependent series
set.seed(1)
x <- as.numeric(stats::arima.sim(model = list(ar = 0.6), n = 200))

# Three inferential layers
fb <- htdv_fixedb(x, theta0 = 0, B = 0.1, sims = 5000, seed = 1)
fit <- htdv_fit(x, model = "tac", chains = 4, iter = 4000, seed = 1)
bb <- htdv_boot(x, mean, R = 1999, type = "stationary", seed = 1)

# Bayesian decision
draws <- as.numeric(rstan::extract(fit$stanfit, "theta")$theta)
htdv_rope(draws, rope = c(-0.1, 0.1))$decision

# Diagnostics
htdv_diagnostics(fit)$passed

A two-sample example via the simulation engine

res <- htdv_simstudy(
  n_grid     = c(80, 200),
  phi_grid   = c(0, 0.6),
  tail_grid  = c("normal", "t3"),
  imb_grid   = c(1, 3),
  delta_grid = c(0, 0.25),
  R          = 50,
  layers     = c("har", "boot", "bayes"),
  seed       = 20260427
)
htdv_simstudy_summary(res)
htdv_simstudy_warnings(htdv_simstudy_summary(res), threshold = 0.7)

Browse shipped validation

# 1024-cell factorial simulation summary
data(htdv_sim_summary)
str(htdv_sim_summary)

# Three-dataset external benchmark
data(htdv_empirical_benchmarks)
htdv_empirical_benchmarks

# Full narrative
vignette("HTDV-validation")

4. Mathematical foundations

4.1. Probability space and filtration

All random variables in HTDV live on a single complete filtered probability space

$$(\Omega,\mathcal{F},(\mathcal{F}_t)_{t\in\mathbb{N}},\mathbb{P}).$$

The filtration ${\mathcal F_t}$ is the natural device for encoding "information available up to time $t$" and is referenced explicitly whenever a martingale or mixingale construction is invoked. For each variable $X_i$ the framework treats its triangular-array realization ${X_{i,m}^{(k)}}$, its correlation-weighted realization ${X_i}$ with weights $w_i$, and its mixingale realization ${X_{i,t}}$ as adapted to ${\mathcal{F}_t}$; the three differ only in how the underlying observations are indexed and combined.

4.2. $\alpha$-mixing with polynomial decay

The framework's primary dependence assumption is $\alpha$-mixing with polynomial decay rate $\gamma > 1$. Define the $\alpha$-mixing coefficient between two $\sigma$-fields $\mathcal{A}, \mathcal{B} \subset \mathcal{F}$ as

$$ \alpha(\mathcal{A}, \mathcal{B}) = \sup_{A\in\mathcal{A},,B\in\mathcal{B}} |\mathbb{P}(A\cap B) - \mathbb{P}(A)\mathbb{P}(B)|. $$

For a sequence ${X_t}$, define $\alpha(k) = \sup_{t\ge 1} \alpha(\sigma(X_s: s\le t),, \sigma(X_s: s\ge t+k))$. The series is strong-mixing if $\alpha(k)\to 0$ as $k\to\infty$. The framework requires the decay rate

$$ \alpha(k) \le C,k^{-\gamma}, \qquad \gamma > 1. $$

This rate is sharp enough to ensure that Rosenthal-type moment inequalities (Rio 1995; Shao 1995) and Hoffmann–Jørgensen-type maximal inequalities (Hoffmann–Jørgensen 1974; Peligrad 2002) deliver $O(n^{-1})$ variance control for the partial-sum / weighted-sum functionals on which the three regimes rely.

4.3. Modes of convergence used in the framework

Three modes of convergence are used explicitly. They are not interchangeable, and each theorem in the package's documentation commits to one:

Convergence in probability ($\xrightarrow{p}$): the default for uniform tightness statements (e.g., the proof of asymptotic normality of $\bar X_n$).
Convergence in $L^2$ ($\xrightarrow{L^2}$): the mode used for the WSC and MPC regimes, since they rely on second-moment control. Implies convergence in probability; on a subsequence implies almost-sure convergence.
Almost-sure convergence ($\xrightarrow{\text{a.s.}}$): claimed only when extractable from Kronecker's lemma + Birkhoff–Khinchin (TAC) or from a second Borel–Cantelli step after summability of tail probabilities.

4.4. Reference inequalities

Throughout the framework, the following inequalities are invoked. Each one is referenced to a specific source so that the dependence assumptions under which it applies are unambiguous.

Inequality	Form	Used for	Reference
Rosenthal (mixing)	$\mathbb{E}\|S_n\|^q \le C_q\big(\sum\mathbb{E}\|X_t\|^q + (\sum\mathrm{Var}(X_t))^{q/2}\big)$	TAC, WSC variance bounds	Rio (1995); Shao (1995)
Marcinkiewicz–Zygmund	$\mathbb{E}\|S_n\|^q \le C_q n^{q/2-1}\sum\mathbb{E}\|X_t\|^q$	Heavy-tail concentration	de la Peña & Giné (1999)
Hoffmann–Jørgensen (mixing)	$\Pr(\sup_n\|S_n\|>3M+\varepsilon) \le 4\Pr(\sup_n\|X_n\|>\varepsilon)+4(\Pr(\sup_n\|S_n\|>M))^2$	Mixingale tightness	Hoffmann-Jørgensen (1974); Peligrad (2002)
Concentration under $\tau$-dependence	bounded-difference under mixing	McDiarmid-type bounds	Rio (2000); Samson (2000); Dedecker–Prieur (2005)
Andrews bandwidth	optimal $h_n = c,n^{1/(2\beta+1)}$	HAC LRV estimation	Andrews (1991)
Kiefer–Vogelsang (fixed-$b$)	non-degenerate functional null distribution	Fixed-$b$ HAR	Kiefer & Vogelsang (2005)

5. The three convergence regimes: TAC, WSC, MPC

5.1. Triangular Arrays Convergence (TAC)

A triangular array is a doubly-indexed sequence ${X_{i,m}^{(k)} : 1\le m\le k,,k\in\mathbb{N}}$ where row $k$ has $k$ entries. The TAC regime claims

$$ \bar X_k = \frac{1}{k}\sum_{m=1}^{k} X_{i,m}^{(k)} \xrightarrow{p} \mu \quad\text{as } k\to\infty, $$

under $\alpha$-mixing with $\gamma>1$ and $\mathbb{E}|X_{i,m}^{(k)}|^q < \infty$ for $q > 2/(\gamma-1)+2$. The TAC regime is most natural for hierarchical panels where increasing aggregation introduces more observations per "row" — input-output disaggregation, sectoral subdivisions, hierarchical strata.

The HTDV likelihood for TAC data is the Stan model tac_hierarchical.stan, parameterized as $X_t \sim \mathcal{N}(\theta + \phi(X_{t-1}-\theta),,\sigma_\varepsilon^2)$ with $X_1 \sim \mathcal{N}(\theta, \sigma_\varepsilon^2/(1-\phi^2))$ at the stationary distribution.

5.2. Weighted Sums with Correlation (WSC)

The WSC regime considers a weighted sum

$$ S_n = \sum_{i=1}^n w_i X_i \xrightarrow{L^2} \mu, \qquad w_i = \frac{1}{\sqrt{n(1+\rho_{ij})}}, $$

where the weights ${w_i}$ encode the degree of contemporaneous correlation $\rho_{ij}$ between component variables. WSC is most natural for cross-sectional dependence with explicit covariance structure — input-output trade flows, networks of related markets, regional sectors that move together.

The convergence is claimed in $L^2$ rather than almost-sure, because the WSC machinery relies on Rosenthal-type variance bounds for the partial sum (Rio 1995, Theorem 2.1; Shao 1995). On a subsequence the convergence strengthens to a.s. by Borel–Cantelli applied to the tail probabilities.

5.3. Mixingale Processes Convergence (MPC)

A mixingale on a filtration ${\mathcal{F}t}$ is a sequence ${X_t}$ with $\mathbb{E}[X_t|\mathcal{F}{t-k}]\le c_t \psi_k$ and $|X_t - \mathbb{E}[X_t|\mathcal{F}_{t+k}]|_2 \le c_t \psi_k$ for some sequence $\psi_k\downarrow 0$ and constants $c_t$. The MPC regime considers the weighted sum

$$ \mathbb{E}\Bigg[\Bigg(\sum_{t=1}^T \psi_t (X_{i,t}-\mu)\Bigg)^2\Bigg] \le C_T $$

where $C_T = O(T^\rho)$ for $\rho < 2$. The MPC regime is the natural fit for temporally-dependent series with decreasing memory — forecast errors, residuals from time-series models, prediction intervals that gradually lose information about the past.

The corresponding HTDV likelihood is mpc_hierarchical.stan, which fits the mixingale dependence weights $\psi_t$ jointly with the long-run mean $\theta$ and innovation scale.

5.4. Why three regimes, not one?

Each regime captures a structurally different way that a real-world dataset can fail i.i.d. assumptions:

TAC: hierarchical aggregation. Information accumulates as you slice the index more finely.
WSC: cross-sectional correlation. Information shares a common covariance structure.
MPC: temporal memory. Information decays smoothly as you reach back in time.

A single-regime framework would be either too narrow (missing the other two failure modes) or too general (so weak that the finite-sample constants are vacuous). The three-regime architecture lets the framework dispatch the correct likelihood to the correct data-generating hypothesis while preserving a unified inferential pipeline downstream.

6. The metric equivalence theorem

The decisive theoretical claim of HTDV is that the three regimes above induce three strictly equivalent metrics on the space of hypothesis-testing problems, under explicit finite-sample constants. This is what makes the unified framework defensible: a result obtained under one regime translates, with computable error, to either of the other two.

6.1. The unified metric

For $j\in{\text{TAC}, \text{WSC}, \text{MPC}}$ and a sample of size $T$, define

$$ D_j^{(T)}(X_i, \mu) = \begin{cases} \bigg|\dfrac{1}{T}\sum_{m=1}^{T} X_{i,m}^{(T)}-\mu\bigg| & j=\text{TAC},\[1.0em] \bigg|\sum_{i=1}^T w_i X_i - \mu\bigg| & j=\text{WSC},\[1.0em] \sqrt{\mathbb{E}\bigg[\big(\sum_{t=1}^T \psi_t(X_{i,t}-\mu)\big)^2\bigg]} & j=\text{MPC}. \end{cases} $$

Each $D_j^{(T)}$ is a non-negative scalar measuring the "distance" from the sample functional to the hypothesized population mean.

6.2. The strict equivalence theorem (informal statement)

Under $\alpha$-mixing with $\gamma > 1$, finite moments $\mathbb{E}|X_i|^q$ for $q\in{4, 6, \infty}$, and finite sample size $T$, there exist explicit constants $K_{\text{TAC}}, K_{\text{WSC}}, K_{\text{MPC}}, c_L, c_U$ — computable as functions of $\gamma, q, T$ — such that for all admissible $X_i$ and $\mu$,

$$ c_L \cdot D_{\text{WSC}}^{(T)}(X_i,\mu) \le D_{\text{TAC}}^{(T)}(X_i,\mu) \le c_U \cdot D_{\text{WSC}}^{(T)}(X_i,\mu), $$

$$ K_{\text{TAC}} \cdot D_{\text{MPC}}^{(T)}(X_i,\mu) \le D_{\text{TAC}}^{(T)}(X_i,\mu), $$

$$ D_{\text{TAC}}^{(T)}(X_i,\mu) \le K_{\text{MPC}} \cdot D_{\text{MPC}}^{(T)}(X_i,\mu). $$

The constants $K_\cdot, c_L, c_U$ are exposed by the function htdv_equivalence_constants(), and they satisfy:

$K_{\text{TAC}}, K_{\text{MPC}}$ decreasing in $\gamma$ (faster mixing decay tightens the bound).
$K_{\text{TAC}}, K_{\text{MPC}}$ increasing in $q$ (more moments loosen the bound).
$c_L \to 1$ and $c_U \to 1$ as $T \to \infty$ (the equivalence becomes asymptotically isometric).

6.3. Why this matters for the user

The metric equivalence theorem is the reason HTDV's three-layer inferential architecture is mathematically well-defined. Without it, running HAR-Wald on a TAC dataset and a Bayesian fit on a WSC dataset would compare apples and oranges. The theorem certifies that the metrics are coercible to one another with computable error, so all three layers can be evaluated against the same null hypothesis.

The function htdv_equivalence_constants() returns the constants for any admissible $(\gamma, q, T)$ configuration:

htdv_equivalence_constants(gamma = 2, q = 6, n = 500)
# $K_TAC, $K_MPC, $c_L, $c_U

For $\gamma=2, q=6, n=500$, $K_{\text{TAC}}\approx 1.45$ and $c_U\approx 1.18$, so a TAC-derived inferential conclusion translates to a WSC-equivalent conclusion within an 18% slack, a margin that is usually irrelevant for hypothesis testing decisions.

7. The Bayesian core

7.1. Likelihood layer: Whittle and composite

HTDV supports two parametric likelihood families for the underlying dependence structure:

Whittle likelihood (frequency domain). For a stationary process with spectral density $f_\theta(\lambda)$ and observed periodogram $I_n(\lambda) = (2\pi n)^{-1}|\sum_{t=1}^n X_t e^{-i\lambda t}|^2$, the Whittle log-likelihood is

$$ \ell_W(\theta) = -\frac{1}{2}\sum_{j=1}^{n/2} \bigg(\log f_\theta(\lambda_j) + \frac{I_n(\lambda_j)}{f_\theta(\lambda_j)}\bigg). $$

Whittle is consistent and asymptotically Gaussian under regularity, with the advantage that it is not vulnerable to the boundary effects that plague time-domain Gaussian likelihoods on short series. It is the canonical likelihood for stationary dependence with a known spectral form. Available in HTDV via model = "whittle" in htdv_fit().

Composite likelihood (time domain). For a process with conditional density $f_\theta(X_t | X_{t-1}, \dots, X_{t-k})$ over a $k$-block, the composite log-likelihood is

$$ \ell_C(\theta) = \sum_{t=k+1}^n \log f_\theta(X_t | X_{t-1}, \dots, X_{t-k}). $$

Composite is a robust choice when the spectral form is uncertain or when the data have non-Gaussian innovations; the trade-off is mild inefficiency relative to the full likelihood (Varin, Reid & Firth 2011). Available in HTDV via model = "composite".

The two likelihoods can be computed individually via htdv_whittle() and htdv_composite(), or stacked in the Berger-robust envelope of Section 7.5.

7.2. Hierarchical priors on the nuisance triple

The dependence-aware likelihood is parameterized by a triple $(\sigma_\infty, \gamma, b)$ where $\sigma_\infty$ is the long-run standard deviation, $\gamma$ is the mixing decay rate, and $b$ is the bandwidth fraction (for HAR-related calculations). The hierarchical priors used by HTDV are:

$$ \sigma_\infty \sim \text{half-Cauchy}(0, 5), \quad \gamma \sim \text{Gamma}(2, 0.5),\text{shifted to support},(1, \infty), $$

$$ b \sim \text{Beta}(2, 8) ,\text{(concentrated near small bandwidth fractions)}. $$

These are weakly informative priors (Gelman 2006) that respect the admissible parameter ranges without overwhelming the data. The practitioner can override them for problem-specific knowledge.

7.3. HMC inference and the diagnostic gate

Posterior sampling uses the No-U-Turn Sampler (Hoffman & Gelman 2014) implemented in Stan via rstan. After sampling, every fit must pass five diagnostic gates before its posterior is admissible:

$\hat R < 1.01$ (split-chain $\hat R$ from Vehtari et al. 2021).
Bulk ESS $> 400$.
Tail ESS $> 400$.
Zero post-warmup divergences.
E-BFMI $> 0.3$ (energy Bayesian fraction of missing information, Betancourt 2016).

The function htdv_diagnostics() returns a boolean passed flag and a character vector of failed checks.

7.4. Decision layer

Given a clean posterior, the framework provides four decision tools:

ROPE (Region of Practical Equivalence; Kruschke 2018). Given a parameter posterior $p(\theta|x)$ and a tolerance interval $[a, b]$, the decision rule

$$ d(\text{ROPE}) = \begin{cases} \text{accept} & \text{if HDI}_{1-\alpha} \subset [a,b] \\ \text{reject} & \text{if HDI}_{1-\alpha} \cap [a,b] = \emptyset \\ \text{undecided} & \text{otherwise} \end{cases} $$

is implemented as htdv_rope().

Bridge-sampling Bayes factor (Meng & Wong 1996; Gronau et al. 2017). For two model fits $M_1, M_0$ over the same data, the marginal likelihood ratio $\text{BF}_{10} = p(x|M_1)/p(x|M_0)$ is approximated by bridge sampling and exposed via htdv_bf(). Requires the bridgesampling package.

WAIC and PSIS-LOO-LFO. Watanabe-Akaike Information Criterion (Watanabe 2010) and leave-future-out cross-validation (Bürkner, Gabry & Vehtari 2020) for predictive comparison are delivered by htdv_waic_lfo(). Requires the loo package.

Predictive stacking (Yao, Vehtari, Simpson & Gelman 2018). Given a list of fits, optimal convex stacking weights for predictive distributions are computed by htdv_stack().

7.5. The Berger-robust envelope

The function htdv_envelope() takes a list of fitted models (typically Whittle and composite, or TAC/WSC/MPC) and combines their posteriors into a robust envelope:

$$ \theta_{\text{env}} \mid x \sim \text{worst-case}, \bigl[p_1(\theta \mid x),, p_2(\theta \mid x),, \dots\bigr] $$

where "worst-case" is operationalized as the union of pre-specified credible regions across the component posteriors (Berger 1985, 1994). The envelope is wider than any single posterior but tighter than the naive max-min, and it gives the practitioner a hedge against the worst-modelled-likelihood scenario.

8. Finite-sample frequentist anchors

These anchors do not depend on a likelihood. They provide a parallel inferential channel that the user can compare against the Bayesian core to detect specification problems.

8.1. HAC long-run variance: `htdv_lrv()`

For a stationary ${X_t}$ with autocovariances $\gamma(k)=\mathrm{Cov}(X_t, X_{t+k})$, the long-run variance is

$$ \sigma_\infty^2 = \sum_{k=-\infty}^{\infty} \gamma(k). $$

The Newey–West (1987) and Andrews (1991) estimators replace the infinite sum with a kernel-weighted finite sum:

$$ \hat\sigma_\infty^2 = \hat\gamma(0) + 2\sum_{k=1}^{n-1} K!\left(\tfrac{k}{h}\right)\hat\gamma(k), $$

where $K$ is a kernel (Bartlett, Parzen, Quadratic-Spectral) and $h$ is a bandwidth. htdv_lrv() implements this with three kernels and the data-driven Andrews bandwidth (which is asymptotically optimal for the mean-squared error of $\hat\sigma_\infty^2$).

out <- htdv_lrv(x, kernel = "qs", bandwidth = "andrews")
out$lrv      # estimate
out$bandwidth  # bandwidth used
out$kernel     # kernel name

8.2. Fixed-bandwidth HAR Wald: `htdv_fixedb()`

The classical HAR-Wald statistic is

$$W = n(\bar X - \theta_0)^2 / \hat\sigma_\infty^2,$$

with asymptotic $\chi^2_1$ critical values. Kiefer & Vogelsang (2005) showed that under fixed-$b$ asymptotics — i.e., letting the bandwidth grow with the sample as $h_n = b,n$ for a fixed fraction $b\in(0,1]$ — the statistic converges to a non-degenerate functional of Brownian motion, $\Lambda_B$, whose critical values can be tabulated by simulation.

The advantage of fixed-$b$ is that the asymptotic distribution is non-pivotal in a way that automatically corrects for the bandwidth- dependent finite-sample distortion of standard HAR. htdv_fixedb() simulates $\Lambda_B$ critical values from $5,000$ Monte Carlo draws by default and returns the test statistic, critical value, $p$-value, and rejection flag.

fb <- htdv_fixedb(x, theta0 = 0, B = 0.1, kernel = "bartlett",
                  sims = 5000, alpha = 0.05, seed = 1)
fb$statistic
fb$p_value
fb$reject

8.3. Stationary block bootstrap: `htdv_boot()`

The stationary bootstrap of Politis & Romano (1994) draws blocks of random geometric length to preserve the dependence structure of the original series. The optimal block length under MSE-optimality (Patton, Politis & White 2009) is computed automatically when block_length = "auto".

The function returns three confidence intervals:

Percentile:

$$[\hat F^{-1}(\alpha/2),\hat F^{-1}(1-\alpha/2)]$$

Basic:

$$[2\hat\theta - \hat F^{-1}(1-\alpha/2),,2\hat\theta - \hat F^{-1}(\alpha/2)]$$

Studentized:

$$\hat\theta \pm z_{1-\alpha/2}\hat\sigma^*$$

bb <- htdv_boot(x, statistic = mean, R = 1999L,
                type = "stationary", block_length = "auto",
                level = 0.95, seed = 1)
bb$ci_percentile
bb$block_length

8.4. Adaptive conformal inference: `htdv_conformal()`

For an online predictor $\hat f$ and a target coverage $1-\alpha$, the adaptive conformal inference of Gibbs & Candès (2021) updates the significance level $\alpha_t$ at each step by

$$ \alpha_{t+1} = \alpha_t + \lambda(\alpha - \mathbf{1}{Y_t \notin C_t}), $$

where $C_t$ is the prediction set at time $t$ and $\lambda$ is a learning-rate hyperparameter. The procedure delivers $(1-\alpha-O(n^{-1/2}))$-coverage prediction sets even under distribution shift, with no parametric assumption.

pred <- function(history) if (length(history) >= 1L) history[length(history)] else 0
conf <- htdv_conformal(x, pred, alpha_target = 0.1, lambda = 0.05, burn_in = 20L)
conf$coverage

9. Heavy-tail and Rényi refinements

9.1. Marcinkiewicz–Zygmund concentration

When the second moment is at the boundary of existence (e.g., Student-$t$ with $\nu$ slightly above $2$), the standard Bernstein-type concentration inequality

$$ \Pr(|S_n| > n\varepsilon) \le \exp!\left(-\frac{n\varepsilon^2}{2(\sigma^2+M\varepsilon/3)}\right) $$

becomes vacuous because $\sigma^2$ is large or undefined. The Marcinkiewicz–Zygmund inequality (de la Peña & Giné 1999) replaces this with a polynomial bound

$$ \Pr(|S_n| > n\varepsilon) \le \frac{C_q}{n^{q/2-1}\varepsilon^q}\sum_{t=1}^{n}\mathbb{E}|X_t|^q $$

valid for $q\in(1, 2]$ even when the variance is infinite. HTDV applies this refinement automatically when the user's data fail the finite-variance check inside the framework's preprocessing.

9.2. Rényi $- \alpha$ divergence tightenings

The KL divergence between two distributions is the limiting case $\alpha\to 1$ of the Rényi $- \alpha$ family:

$$ D_\alpha(P|Q) = \frac{1}{\alpha-1}\log\int p^\alpha q^{1-\alpha}. $$

Bounds expressed in $D_\alpha$ for $\alpha\ne 1$ can be tighter than KL for some inferential decisions, and they relate cleanly to the total- variation and Hellinger distances. The framework uses Rényi $- \alpha$ tightenings (Van Erven & Harremoës 2014) in the bridge-sampling Bayes factor when KL would be infinite or undefined.

10. Function reference

Core fitting

`htdv_fit()`

Hierarchical Bayesian fit via HMC. Compiles a Stan model on first call and caches it for the session.

htdv_fit(x,
         model = c("tac", "wsc", "mpc", "whittle", "composite"),
         chains = 4, iter = 4000, warmup = floor(iter / 2),
         refresh = 0, seed = NULL,
         priors = NULL, control = NULL, ...)

Arguments:

x: numeric vector or matrix of observations.
model: which dependence likelihood to fit.
chains, iter, warmup: HMC hyperparameters.
priors: named list of prior overrides (defaults are weakly informative per Section 7.2).
control: passed to rstan::sampling() (e.g., adapt_delta, max_treedepth).

Returns: an htdv_fit S3 object with $stanfit, $model, $n, $priors. Diagnostic-pass status is computable via htdv_diagnostics().

fit <- htdv_fit(x, model = "tac", chains = 4, iter = 4000, seed = 1)
htdv_diagnostics(fit)$passed

`htdv_envelope()`

Berger-robust envelope across fitted models.

htdv_envelope(fits, target = "theta", level = 0.95, ...)

fits is a list of htdv_fit objects; the function returns the envelope posterior on the parameter target at credibility level.

`htdv_summary()`

Summary of an htdv_fit posterior.

htdv_summary(fit, target = "theta", rope = NULL, level = 0.95)

Returns posterior mean, HDI, and (optionally) ROPE-decision flag.

Likelihoods (standalone)

`htdv_whittle()`

Whittle log-likelihood evaluation.

htdv_whittle(x, spec_fun, theta, taper = FALSE)

spec_fun is a user-supplied function returning the spectral density $f_\theta(\lambda)$ at the Fourier frequencies. Returns $loglik and $periodogram.

`htdv_composite()`

Composite log-likelihood evaluation.

htdv_composite(x, density_fun, theta, k = 1L, log = FALSE)

density_fun(x_t, x_{t-1}, ..., x_{t-k}, theta) is the user-supplied $k$-block conditional density. Returns $loglik.

Long-run variance and HAR

`htdv_lrv()`

HAC long-run variance with Andrews (1991) optimal bandwidth.

htdv_lrv(x, kernel = c("bartlett", "parzen", "qs"),
         bandwidth = "andrews", prewhiten = FALSE)

Returns $lrv, $bandwidth, $kernel, and (if prewhiten = TRUE) $ar_coef for the AR(1) prewhitening step.

`htdv_fixedb()`

Fixed-bandwidth HAR Wald test (Kiefer & Vogelsang 2005) with simulated null distribution.

htdv_fixedb(x, theta0 = 0, B = 0.1,
            kernel = c("bartlett", "parzen", "qs"),
            sims = 5000L, alpha = 0.05, seed = NULL)

B is the bandwidth fraction; sims is the number of Monte Carlo draws from $\Lambda_B$. Returns $statistic, $critical_value, $p_value, $reject.

Bootstrap

`htdv_boot()`

Stationary or circular block bootstrap.

htdv_boot(x, statistic, R = 1999L,
          type = c("circular", "stationary"),
          block_length = "auto",
          level = 0.95, seed = NULL)

block_length = "auto" invokes the Patton-Politis-White rule. Returns $t0, $t_star, three CI types, and $block_length.

Adaptive conformal

`htdv_conformal()`

Online adaptive conformal inference for dependent data (Gibbs & Candès 2021).

htdv_conformal(x, pred, alpha_target = 0.1, lambda = 0.05,
               burn_in = 50L)

pred(history) is a user-supplied scalar predictor; the function returns the per-step prediction set, the empirical coverage, and the trajectory of $\alpha_t$.

Decision tools

`htdv_rope()`

Region of Practical Equivalence decision.

htdv_rope(draws, rope, level = 0.95)

draws is a numeric vector of posterior draws of the parameter of interest; rope is a length-$2$ tolerance interval. Returns $decision ("accept" | "reject" | "undecided"), $hdi, $rope, $prob_in_rope.

`htdv_bf()`

Bridge-sampling Bayes factor between two htdv_fit objects.

htdv_bf(fit1, fit0, ...)

Requires the bridgesampling package. Returns $bf10, $log_bf10, $logml1, $logml0, $error_percentage.

`htdv_waic_lfo()`

WAIC and leave-future-out CV for an htdv_fit object.

htdv_waic_lfo(fit, L = NULL, k_threshold = 0.7)

Requires the loo package. L is the leave-future-out window length (default $\lfloor n/5 \rfloor$). Returns $waic, $elpd_lfo, $elpd_se, $pareto_k, $refresh_times.

`htdv_stack()`

Predictive stacking weights across fits.

htdv_stack(fits, method = c("log_score"))

Requires the loo package. Returns $weights, $model_names, $elpd_per_model.

Diagnostics and posterior-predictive checks

`htdv_diagnostics()`

HMC diagnostic gate.

htdv_diagnostics(fit, rhat_threshold = 1.01,
                 ess_threshold = 400, bfmi_threshold = 0.3)

Returns $rhat, $ess_bulk, $ess_tail, $bfmi, $divergences, $max_treedepth, $passed (logical), $failures (character vector of failed checks).

`htdv_ppc()`

Posterior-predictive checks on dependence statistics (autocorrelation, Ljung-Box, spectral discrepancy, Hurst exponent).

htdv_ppc(fit, x_obs, lag_max = 10L)

Returns $p_acf, $p_ljung_box, $p_spectral, $p_hurst. The posterior-predictive $p$-values should not be too close to $0$ or $1$ (values in $[0.05, 0.95]$ indicate good model fit on the dependence structure).

Equivalence constants

`htdv_equivalence_constants()`

The explicit constants $K_{\text{TAC}}, K_{\text{WSC}}, K_{\text{MPC}}, c_L, c_U$ for the metric equivalence theorem.

htdv_equivalence_constants(gamma, q, n)

gamma is the polynomial mixing decay rate; q is the moment index; n is the sample size. Returns the named list of constants.

htdv_equivalence_constants(gamma = 2, q = 6, n = 500)
# $K_TAC, $K_WSC, $K_MPC, $c_L, $c_U

Simulation engine

`htdv_simstudy()`

Factorial Monte Carlo runner. Five factors crossed at four levels each:

htdv_simstudy(n_grid     = c(40L, 80L, 200L, 500L),
              phi_grid   = c(0, 0.3, 0.6, 0.85),
              tail_grid  = c("normal", "t5", "t3", "t2_1"),
              imb_grid   = c(1, 1.5, 3, 6),
              delta_grid = c(0, 0.1, 0.25, 0.5),
              R          = 500L,
              seed       = 20260422L,
              n_cores    = 1L,
              layers     = c("har", "boot", "bayes"),
              bayes_chains = 2L, bayes_iter = 600L, bayes_warmup = 300L,
              boot_R     = 999L, rope = c(-0.1, 0.1), alpha = 0.05,
              out_dir    = NULL, progress = TRUE)

out_dir = NULL returns results in memory; passing a directory path writes one RDS per cell atomically (cell-level resume support; survives worker crashes). Per-cell parallelism via parallel::mclapply when n_cores > 1 on Unix.

Returns a tidy data frame with one row per (cell × layer × replication).

`htdv_simstudy_summary()`

Aggregate per-cell results.

htdv_simstudy_summary(res, alpha = 0.05)

Returns one row per (cell × layer) with columns reject_rate, coverage, mean_ci_length, mean_runtime, diag_pass_rate.

`htdv_simstudy_warnings()`

Flag cells in the limit-of-identification zone.

htdv_simstudy_warnings(summ, threshold = 0.7)

Returns the subset of cells where Bayesian diag_pass_rate < threshold. Default $0.7$ catches the high-persistence × small-sample corner ($\phi=0.85, n\le 80$) plus a fringe of $\phi=0.6, n=40$. Practitioners should treat flagged cells with caution: increase HMC iterations, tighten $\phi$ priors, or fall back on the conformal anchor.

11. Validation evidence

11.1. The factorial simulation: `htdv_sim_summary`

A pre-registered factorial Monte Carlo of 1,024 cells × 500 replications × 3 inferential layers, executed in 31.1 hours on a 16-core Linux workstation. The five factors are:

$n \in {40, 80, 200, 500}$
$\phi \in {0, 0.3, 0.6, 0.85}$
innovation tail $\in {\mathcal{N}(0,1), t_5, t_3, t_{2.1}}$
imbalance ratio $n_1/n_2 \in {1, 1.5, 3, 6}$
location shift $\Delta \in {0, 0.1, 0.25, 0.5}\sigma_\infty$

Headline finding (sign-corrected):

Layer	Mean size	Mean coverage	Min coverage
bayes	0.056	0.944	0.884
HAR	0.180	0.820	0.286
bootstrap	0.175	0.824	0.386

The Bayesian envelope is the only layer with empirical size and coverage close to nominal across the entire grid. HAR and bootstrap inflate dramatically under high persistence (size up to 0.71; coverage down to 0.29). The full report is in vignette("HTDV-validation").

library(HTDV)
data(htdv_sim_summary)
str(htdv_sim_summary)   # 3069 rows x 12 columns

11.2. The empirical benchmarks: `htdv_empirical_benchmarks`

Three external validation runs against published references on public- source data:

Dataset	$\widehat\phi$	Reference	Verdict
FRED-MD CPI inflation post-1984	$\approx 0.45$	Stock & Watson (2007)	agreement
Shiller log-CAPE 1881-latest	$\approx 0.97$	Campbell & Shiller (1998)	agreement
US-Canada 10y yield differential	$\approx 0.99$	Naive iid Welch baseline	agreement

All three references reproduce within the framework's 95% credible interval. The Bayesian envelope half-width grows monotonically with the underlying persistence (0.46 → 0.70 → 7.35), consistent with the finite-sample inflation predicted by the metric equivalence theorem.

data(htdv_empirical_benchmarks)
htdv_empirical_benchmarks[, c("dataset", "agreement_bayes")]

12. Tutorial walkthroughs

12.1. For the novice (early-career student)

You have a time series of monthly inflation rates for the United States from 1984 to today, and you want to know whether the long-run mean is zero or not. The classical $t$-test is wrong because the observations are autocorrelated. Here is what HTDV does for you:

library(HTDV)
# (Assume `infl` is a numeric vector of monthly inflation rates,
#  e.g., from FRED-MD's CPIAUCSL transformed to 1200 * Δlog(CPI))

# Step 1: compute the long-run variance correctly
lrv <- htdv_lrv(infl, kernel = "qs", bandwidth = "andrews")$lrv
# This is what replaces the classical sample variance when data is dependent.

# Step 2: HAR-Wald test — is the mean equal to zero?
fb <- htdv_fixedb(infl, theta0 = 0, B = 0.1, sims = 5000, seed = 1)
fb$p_value           # < 0.05 means reject H0: mean = 0
fb$reject

# Step 3: Bayesian fit — what's the posterior on the mean?
fit <- htdv_fit(infl, model = "tac", chains = 4, iter = 4000, seed = 1)
htdv_summary(fit, target = "theta", level = 0.95)
# Posterior mean, 95% credible interval, ROPE decision

# Step 4: Sanity check — did the chains converge?
htdv_diagnostics(fit)$passed   # TRUE means OK

If steps 2 and 3 agree (e.g., both reject H0 with similar effect-size estimates), you have robust evidence. If they disagree, run step 5:

# Step 5: A bootstrap second opinion
bb <- htdv_boot(infl, mean, R = 1999, type = "stationary",
                block_length = "auto", seed = 1)
bb$ci_percentile

If two of the three layers agree, that is your conclusion. If the Bayesian gives a much wider interval than HAR or bootstrap, the data likely have strong persistence: trust the Bayesian.

12.2. For the applied statistician

You have a panel of sectoral profitability rates with $T=200$ observations per sector and $N=20$ sectors of unequal size (some have $T_i=200$, others $T_i = 50$ due to data availability). You want to test whether the sectors converge to a common mean.

library(HTDV)
# `panel` is a list of length 20, each a numeric vector of length T_i

# Run the simulation engine on a small calibration grid first
calib <- htdv_simstudy(
  n_grid     = c(50, 200),
  phi_grid   = c(0, 0.6),     # informed by panel autocorrelation
  tail_grid  = c("normal", "t5"),
  imb_grid   = c(1, 4),       # span your imbalance
  delta_grid = c(0, 0.25),
  R          = 100,
  layers     = c("har", "boot", "bayes"),
  seed       = 1
)
htdv_simstudy_warnings(htdv_simstudy_summary(calib), threshold = 0.7)
# This tells you whether your (n, phi) configuration is in the
# limit-of-identification zone.

# Now test each pair: is sector i's mean equal to sector j's?
# Use the Berger-robust envelope:
fits <- lapply(panel, function(x) htdv_fit(x, model = "tac",
                                            chains = 4, iter = 4000, seed = 1))
env <- htdv_envelope(fits, target = "theta", level = 0.95)
env$intervals

12.3. For the mathematician

You want to verify the metric equivalence constants for your specific $(\gamma, q, n)$ triple before deploying the framework on a critical dataset.

library(HTDV)

# Compute the constants for several (gamma, q, n)
configs <- expand.grid(gamma = c(1.5, 2, 3, 4),
                       q = c(4, 6),
                       n = c(100, 500, 5000))
constants <- t(apply(configs, 1, function(r) {
  out <- htdv_equivalence_constants(gamma = r["gamma"], q = r["q"],
                                     n = r["n"])
  c(out$K_TAC, out$K_MPC, out$c_L, out$c_U)
}))
colnames(constants) <- c("K_TAC", "K_MPC", "c_L", "c_U")
cbind(configs, constants)
# Verify: K_TAC decreasing in gamma, increasing in q;
#         c_L -> 1, c_U -> 1 as n grows.

For the proof apparatus, see Sections 4–6 of the companion paper vignette("HTDV-validation"). The Rosenthal-type inequalities underlying the WSC bound use Rio (1995) and Shao (1995); the Hoffmann–Jørgensen variant for mixingales uses Peligrad (2002) with Rio (2000) for the polynomial-decay sharpening; the heavy-tail refinements rely on de la Peña & Giné (1999); the bounded-difference McDiarmid extensions to the dependent setting use Rio (2000), Samson (2000), and Dedecker & Prieur (2005).

13. Frequently asked questions

Q. Why three regimes (TAC, WSC, MPC) instead of one general theorem? A. Because each regime captures a structurally different way that real-world data fail i.i.d.: hierarchical aggregation (TAC), cross-sectional correlation (WSC), or temporal memory (MPC). A single-regime framework would either miss two of these failure modes or be too general for finite-sample constants to be useful. The metric equivalence theorem of Section 6 unifies them at the asymptotic level while preserving the finite-sample specificity.

Q. Is the Bayesian envelope just a heavier hammer? A. No — see Section 11.1. The simulation shows that under high persistence ($\phi \ge 0.6$) and small samples ($n \le 200$), the asymptotic critical values used by HAR and bootstrap are unreliable, and only the Bayesian envelope maintains nominal coverage. The Bayesian intervals are wider, but the frequentist intervals achieve their narrower width by miscalibrating, not by being more informative.

Q. What if my data have unit roots? A. The framework's AR(1) Stan model will fit near-unit-root data with $\hat\phi$ close to one. The posterior on the long-run mean becomes correspondingly wide, which is the honest answer. For an explicit unit-root test, use the standard ADF or Phillips-Perron statistics; the framework is for inference under stationary mixing, not for testing against the unit-root null.

Q. How much CPU time does htdv_simstudy() need? A. The full default grid (1024 cells × 500 replications × 3 layers) runs in 31.1 hours on a 16-core Linux workstation. Smaller grids scale roughly linearly. Bayesian fits at $n=500$ are the bottleneck (~25s each); HAR and bootstrap are sub-second.

Q. Why are my HMC diagnostics failing for $n=40, \phi=0.85$? A. The AR(1) likelihood at $n=40$ with $\phi$ near 1 is at its limit of identification. The chains can struggle to explore the joint posterior on $(\theta, \phi)$. Increase iter to 8000+, tighten the prior on $\phi$ toward zero, or fall back on the conformal anchor whose validity does not require identifiability. See htdv_simstudy_warnings() for an automated diagnostic-pass-rate filter.

Q. Can I use HTDV for spatial data? A. The current implementation focuses on temporal $\alpha$-mixing. Spatial mixing introduces additional complexities (anisotropy, irregular sampling) that are not yet fully wired through the package. The MPC likelihood with custom $\psi_t$ weights can be adapted to a 1-D spatial index; for full 2-D spatial inference, use complementary tools.

Q. The HTDV validation card says HAR has bad coverage. Why ship HAR at all? A. HAR remains the right tool when persistence is low to moderate ($\phi \le 0.5$) and sample sizes are large enough for the asymptotics to bite ($n \ge 500$). It is also computationally an order of magnitude faster than the Bayesian fit. The framework's three-layer design lets the user exploit HAR's speed and simplicity in the regime where it is calibrated, while flagging the cells where it is not.

Q. Where is cran-comments.md referenced from the package? A. It isn't — cran-comments.md is a CRAN-internal communication file not intended for end users. The CRAN release notes are in NEWS.md.

14. References

Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation consistent covariance matrix estimation. Econometrica, 59(3), 817–858.

Berger, J. O. (1994). An overview of robust Bayesian analysis. Test, 3(1), 5–124.

Betancourt, M. (2016). Diagnosing suboptimal cotangent disintegrations in Hamiltonian Monte Carlo. arXiv:1604.00695.

Brown, B. M. (1971). Martingale central limit theorems. Annals of Mathematical Statistics, 42(1), 59–66.

Bürkner, P.-C., Gabry, J., & Vehtari, A. (2020). Approximate leave-future-out cross-validation for Bayesian time series models. Journal of Statistical Computation and Simulation, 90(14), 2499–2523.

Campbell, J. Y., & Shiller, R. J. (1998). Valuation ratios and the long-run stock market outlook. Journal of Portfolio Management, 24(2), 11–26.

Davidson, J. (1994). Stochastic Limit Theory. Oxford University Press.

Dedecker, J., & Prieur, C. (2005). New dependence coefficients. Probability Theory and Related Fields, 132(2), 203–236.

de la Peña, V. H., & Giné, E. (1999). Decoupling: From Dependence to Independence. Springer.

Durbin, J., & Koopman, S. J. (2001). Time Series Analysis by State Space Methods. Oxford University Press.

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515–534.

Gibbs, I., & Candès, E. (2021). Adaptive conformal inference under distribution shift. NeurIPS, 34, 1660–1672.

Gronau, Q. F. et al. (2017). A tutorial on bridge sampling. Journal of Mathematical Psychology, 81, 80–97.

Hall, P., & Heyde, C. C. (1980). Martingale Limit Theory and Its Application. Academic Press.

Hoffmann-Jørgensen, J. (1974). Sums of independent Banach space valued random variables. Studia Mathematica, 52(2), 159–186.

Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn sampler. JMLR, 15(1), 1593–1623.

Kiefer, N. M., & Vogelsang, T. J. (2005). A new asymptotic theory for HAR tests. Econometric Theory, 21(6), 1130–1164.

Kruschke, J. K. (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270–280.

McCracken, M. W., & Ng, S. (2016). FRED-MD: A monthly database for macroeconomic research. JBES, 34(4), 574–589.

Newey, W. K., & West, K. D. (1987). A simple, positive semi-definite HAC covariance matrix. Econometrica, 55(3), 703–708.

Patton, A., Politis, D. N., & White, H. (2009). Correction to "Automatic block-length selection for the dependent bootstrap". Econometric Reviews, 28(4), 372–375.

Peligrad, M. (2002). Some remarks on coupling of dependent random variables. Statistics & Probability Letters, 60(2), 201–209.

Politis, D. N., & Romano, J. P. (1994). The stationary bootstrap. JASA, 89(428), 1303–1313.

Rao, C. R. (1973). Linear Statistical Inference and Its Applications (2nd ed.). Wiley.

Rio, E. (1995). A maximal inequality and dependent Marcinkiewicz–Zygmund strong laws. Annals of Probability, 23(2), 918–937.

Rio, E. (2000). Théorie asymptotique des processus aléatoires faiblement dépendants. Springer.

Samson, P.-M. (2000). Concentration of measure inequalities for Markov chains and $\Phi$-mixing processes. Annals of Probability, 28(1), 416–461.

Shao, Q.-M. (1995). Maximal inequalities for partial sums of $\rho$-mixing sequences. Annals of Probability, 23(2), 948–965.

Stock, J. H., & Watson, M. W. (2007). Why has US inflation become harder to forecast? Journal of Money, Credit and Banking, 39, 3–33.

Van Erven, T., & Harremoës, P. (2014). Rényi divergence and Kullback–Leibler divergence. IEEE Transactions on Information Theory, 60(7), 3797–3820.

Varin, C., Reid, N., & Firth, D. (2011). An overview of composite likelihood methods. Statistica Sinica, 21(1), 5–42.

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432.

Vehtari, A. et al. (2021). Rank-normalization, folding, and localization: An improved $\hat R$ for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and WAIC. JMLR, 11, 3571–3594.

Whittle, P. (1953). Estimation and information in stationary time series. Arkiv för Matematik, 2(5), 423–434.

Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using stacking to average Bayesian predictive distributions. Bayesian Analysis, 13(3), 917–1007.

License

MIT.

Citation

citation("HTDV")

Issues and contributions

github.com/IsadoreNabi/HTDV/issues