---
title: "The HCbeta algorithm and its consistency"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{The HCbeta algorithm and its consistency}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(hcinfer.use_emoji = FALSE)
```

This vignette presents the HCbeta adjustment factor as an explicit algorithm and
records the argument that it yields a consistent covariance estimator. It
complements `vignette("hcinfer-methodology", package = "hcinfer")`, which places
HCbeta within the broader family of HC estimators, and
`vignette("hcinfer-hcbeta", package = "hcinfer")`, which is usage oriented.

## Motivation

The classical leverage corrections write the adjustment factor as a power of the
leverage complement $u_t = 1 - h_t$. Because $u_t \in [0, 1]$, it can be read as
the value at $1 - h_t$ of the cumulative distribution function of a
Uniform$(0, 1)$ random variable. HCbeta keeps this reading but replaces the
uniform cdf with a Beta cdf evaluated at a truncated complement $w_t$. The two
shape parameters give the correction curve an extra degree of flexibility, so the
adjustment can adapt to the observed leverage configuration instead of following a
fixed exponent.

This flexibility aims to temper the overshooting that HC3, HC4, and HC4m can
display under strong leverage, where small values of $1 - h_t$ raised to a fixed
negative exponent produce very large factors. HCbeta calibrates the Beta shape
parameters to the first two sample moments of $w_t$, so the resulting factor
reflects the actual spread of the leverages rather than a worst case exponent.

HCbeta does not assume that the leverage complements follow a Beta distribution.
The Beta cdf is a calibration device whose shape parameters are matched to the
sample mean and variance of the truncated complements, and the construction would
be unchanged if the leverages arose from another mechanism. The Beta family is
used because its support is $[0, 1]$ and its two parameters are enough to bend
the uniform baseline in either direction.

## The algorithm

Let $X$ be the $n \times p$ design matrix of an OLS fit with residuals
$\hat e_t$ and leverages $h_t = x_t'(X'X)^{-1}x_t$. HCbeta computes the
adjustment factor $g_t$ in four steps.

1. **Truncated leverage complements.** Form $u_t = 1 - h_t$ and truncate it to
   $w_t = \max\{\text{lower}, \min(u_t, \text{upper})\}$, with defaults
   $\text{lower} = 0.01$ and $\text{upper} = 0.99$. The lower bound prevents the
   adjustment from exploding when $h_t$ approaches one, and the upper bound keeps
   the values away from the boundary when leverages are small.

2. **Method of moments.** Estimate the mean and variance of the truncated
   complements, $\hat\mu = n^{-1}\sum_{t=1}^{n} w_t$ and
   $s_w^2 = (n - 1)^{-1}\sum_{t=1}^{n}(w_t - \hat\mu)^2$, and map them to Beta
   shape parameters through the mean precision parameterization,
   $\hat\phi = \hat\mu(1 - \hat\mu)/s_w^2 - 1$, $\hat a = \hat\mu\hat\phi$, and
   $\hat b = (1 - \hat\mu)\hat\phi$.

3. **Shrinkage with a fixed floor and caps.** Shrink the estimates toward the
   uniform case $a = b = 1$ using $\zeta = n/(n + 50)$, then set
   $\tilde a = \min\{\max\{(1 - \zeta) + \zeta\hat a,\, \varepsilon\},\,
   A_{\max}\}$ and
   $\tilde b = \min\{\max\{(1 - \zeta) + \zeta\hat b,\, \varepsilon\},\,
   B_{\max}\}$. The floor is fixed at $\varepsilon = 0.01$, while the defaults
   $A_{\max} = B_{\max} = 10000$ are set through `a_max` and `b_max`, each of
   which accepts finite values in the inclusive range `[50, 25000]`. The floor
   is applied after shrinkage and before the caps. It is distinct from `lower`,
   which truncates $w_t$, and it is not a method argument.

4. **Adjustment factor and sandwich.** Combine the Beta cdf with a decaying
   exponent,
   $g_t = \frac{n}{n - p}\left\{1/F_B(w_t; \tilde a, \tilde b)\right\}^{c_1/n^{c_2}}$,
   with defaults $c_1 = 7$ and $c_2 = 0.75$. The covariance estimator is the
   sandwich
   $\widehat\Psi_{HC\beta} = (X'X)^{-1} X'\widehat\Omega_\beta X (X'X)^{-1}$ with
   $\widehat\Omega_\beta = \operatorname{diag}(\hat e_t^2 g_t)$.

The degrees of freedom factor $n/(n - p)$ enters $g_t$ in the same way as in
HC1. When $\tilde a = \tilde b = 1$ the Beta cdf returns $w_t$, and with
$c_1 = 0$ the factor reduces to $n/(n - p)$, so HCbeta coincides with HC1. This
identity is a convenient implementation check and is exercised in the package
tests.

## Consistency

The consistency of HCbeta follows from the standard argument for HC0 given by
White (1980) for fixed regressors, extended to stochastic regressors in the
usual way. The role of the HCbeta construction is only to show that its
adjustment factor does not disturb that argument asymptotically. The reasoning
below is presented at an accessible level and does not introduce a new theorem.

Two features keep the adjustment under control. First, the truncation to
$[0.01, 0.99]$ and the bounds on the shape parameters place $\tilde a$ and
$\tilde b$ in a compact interval $[\varepsilon, \max(A_{\max}, B_{\max})]$ with
the fixed value $\varepsilon = 0.01$. On that compact set, the Beta cdf
evaluated at $w_t$ is bounded away from zero, so that
$\inf_t F_B(w_t; \tilde a, \tilde b) \geq \delta$
for some $\delta > 0$, and therefore
$\sup_t |\log F_B(w_t; \tilde a, \tilde b)| \leq -\log\delta < \infty$. The caps
make the parameter space compact, keeping the log adjustments uniformly bounded.

The caps are essential when the leverage variance vanishes. Under the standard
fixed-design condition $\max_t h_t \to 0$, the complements eventually reach the
common upper truncation, so $s_w^2 = 0$ and the uncapped method-of-moments shape
estimates diverge. HCbeta then uses $A_{\max}$ and $B_{\max}$ as the limiting
shape parameters instead of aborting. This keeps the Beta cdf and $g_t$
well-defined throughout the sequence; with the default caps, the common cdf is
numerically one and $g_t$ equals the HC1 scale $n/(n - p)$.

Second, a first-order Taylor expansion of $g_t$ around one gives
$\sup_t |g_t - 1| = O(n^{-c_2}) = o(1)$, because the exponent $c_1/n^{c_2}$ tends
to zero while it multiplies the bounded log cdf. It follows that
$\widehat\Omega_\beta$ and the HC0 residual matrix
$\operatorname{diag}(\hat e_t^2)$ are asymptotically equivalent, so
$\widehat\Psi_{HC\beta}$ and $\widehat\Psi_{HC0}$ share the same limit. HCbeta
is asymptotically equivalent to HC0, so it inherits the same consistency, while
the finite sample factor $g_t$ aims to temper the overshooting that fixed
exponent methods can produce under strong leverage.

## Numerical implementation

Four safeguards keep the computation of $g_t$ numerically stable. First, the
post-shrinkage shape parameters are clamped to the fixed floor 0.01 and their
respective caps. Second, the moment branch is used only when
`is.finite(s2_w) && s2_w > .Machine$double.eps`; otherwise, the shapes are set
directly to their caps. Third, the Beta cdf is evaluated with
`pbeta(w, a_tilde, b_tilde, log.p = TRUE)`, which returns $\log F_B$ directly
and avoids the $\log(0) = -\infty$ produced when a cdf underflows. Finally, the
exponent $-(c_1/n^{c_2})\log F_B(w_t; \tilde a, \tilde b)$ is capped at 700
before exponentiation so that `exp()` cannot overflow to infinity.

In the public-schools model shown here, the default caps
$A_{\max} = B_{\max} = 10000$ are far from binding.

```{r}
library(hcinfer)

schools <- transform(
  PublicSchools,
  income_scaled = income / 10000,
  income_scaled_sq = (income / 10000)^2
)
fit <- lm(expenditure ~ income_scaled + income_scaled_sq, data = schools)
cov <- vcov_hc(fit, type = "hcbeta")

cov$method_params[c("a_tilde", "b_tilde", "a_max", "b_max")]
max(cov$weights)
```

The adjusted shape parameters are well below the caps, and the largest
adjustment factor matches the value reported for this model. By construction each
factor is at least the HC1 scale $n/(n - p)$, because the capped exponent is
nonnegative and the leading factor is exactly $n/(n - p)$. The code below
verifies this lower-bound property.

```{r}
n <- cov$n
p <- cov$p

all(cov$weights >= n / (n - p))
```
