--- title: "The HCbeta algorithm and its consistency" 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)) ```