--- title: "Estimation and inference with hdcce" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Estimation and inference with hdcce} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE, fig.width = 6.5, fig.height = 3.6, fig.align = "center") old <- options(width = 88) set.seed(20260810) ``` ```{r} library(hdcce) ``` # 1. Introduction **hdcce** fits panel data models in which the number of regressors may exceed the number of observations and the unobserved heterogeneity has an interactive structure. It handles two specifications, and offers estimation and inference for each. The **linear model**, for units $i = 1, \dots, n$ and periods $t = 1, \dots, T$, $$Y_{it} = X_{it}^\top \beta + \gamma_i^\top F_t + \varepsilon_{it},$$ with $X_{it} \in \mathbb{R}^p$, unobserved factors $F_t \in \mathbb{R}^K$ and unit-specific loadings $\gamma_i$. The vector $\beta$ is assumed sparse. The **additive model**, in which each covariate enters through an unknown function expanded in a known dictionary, $$Y_{it} = \sum_{j=1}^p m_j(X_{it,j}) + \gamma_i^\top F_t + \varepsilon_{it}, \qquad m_j(x) = \phi_j(x)^\top \beta_j ,$$ with $\phi_j = (\phi_{j1}, \dots, \phi_{jL_j})^\top$. Collecting the transformations of all covariates gives a design with $d = \sum_j L_j$ columns. In both cases the factors and loadings are never estimated. They are removed by projecting the data on the orthogonal complement of the space spanned by the cross-sectional averages of the regressors, after which a lasso is run on the projected data. Which specification is used is determined by a single argument: supplying `dictionaries` switches from the first model to the second, in both `hdcce_estimator` and `hdcce_inference`. ## Data layout **The panel must be sorted by unit.** Rows $(i-1)T + 1, \dots, iT$ of `x`, `y` and `Phi` belong to unit $i$. This is not checked, and a panel sorted by period returns numbers rather than an error. ## The shipped data Two simulated panels come with the package, both produced by `generate_data` with $n = T = 20$ and $p = 61$. ```{r} data("data_estimation") data("data_inference") obs_N <- 20 obs_T <- 20 p <- ncol(data_estimation$x) dim(data_estimation$x) ``` They are generated from a linear design whose coefficient vector has ten non-zero entries: the first regressor, and the first three regressors of each of three groups. ```{r} gsize <- (p - 1) / 3 beta_true <- c(1, rep(c(1, 1, 1, rep(0, gsize - 3)), 3)) which(beta_true != 0) ``` `data_inference$y` has three columns, differing only in the coefficient of the first regressor: $c^{**} = 0$, $0.1$ and $0.2$. The first column therefore has $\beta_1 = 0$, which is used below as a case where the null is true. Throughout, a dictionary with $\phi_j(x) = (x, x^2)$ serves as the example expansion. It is supplied as a matrix of stacked transformations together with a vector recording which covariate each column belongs to. ```{r} make_dict <- function(X) { Phi <- do.call(cbind, lapply(seq_len(ncol(X)), function(j) cbind(X[, j], X[, j]^2))) list(Phi = Phi, group = rep(seq_len(ncol(X)), each = 2)) } dict_est <- make_dict(data_estimation$x) dict_inf <- make_dict(data_inference$x) dim(dict_est$Phi) ``` Since the shipped panels come from a linear design, the true $m_j$ are linear and the quadratic columns of the dictionary have coefficient zero. # 2. Estimation ## 2.1 Linear specification `hdcce_estimator` estimates the factor space from the cross-sectional averages, projects the data on its orthogonal complement, and runs a lasso. ```{r} fit <- hdcce_estimator(data_estimation, obs_N = obs_N, obs_T = obs_T, NFOLDS = 5) fit$K_hat ``` The number of factors is the count of normalised eigenvalues of $\bar{\boldsymbol{X}}^\top \bar{\boldsymbol{X}} / T$ exceeding the truncation `TRUNC`. The returned `eigenvalues` show how clear-cut that decision was. ```{r scree} plot(fit$eigenvalues[1:15], type = "b", pch = 16, ylim = c(0, 1), xlab = "index", ylab = "normalised eigenvalue", main = sprintf("K_hat = %d", fit$K_hat)) abline(h = 0.01, col = "red", lty = 2) legend("topright", "TRUNC", lty = 2, col = "red", bty = "n") ``` A sharp drop after the third eigenvalue is what one hopes to see; a gradual decay means the factor number is not well identified and the result will be sensitive to `TRUNC`. Supplying `NFACTORS` fixes $\widehat{K}$ directly. The estimates recover the sparsity pattern, with the usual shrinkage: ```{r} est <- as.numeric(fit$coefs) round(head(est, 6), 3) c(selected = sum(est != 0), true_nonzero = sum(beta_true != 0), found = sum(est[beta_true != 0] != 0)) ``` All ten non-zero coefficients are selected, alongside a number of false positives. That is expected when the penalty is chosen by cross-validation, which targets prediction rather than selection, and it is the reason inference requires the debiasing step of Section 3 rather than reading significance off the lasso fit. ## 2.2 Dictionary specification Supplying `dictionaries` estimates the additive model instead. The returned coefficients are the $d$ dictionary coefficients rather than $p$ slopes, named by covariate and position within its block. ```{r} fit_d <- hdcce_estimator(data_estimation, obs_N = obs_N, obs_T = obs_T, dictionaries = dict_est, NFOLDS = 5) c(K_hat = fit_d$K_hat, n_coef = length(as.numeric(fit_d$coefs))) round(head(fit_d$coefs, 6), 3) ``` # 3. Inference ## 3.1 Linear specification: confidence intervals `hdcce_inference` returns a desparsified estimate for each coefficient in `COEF_INDEX_VEC`. The debiasing runs a nodewise lasso of the target regressor on the others and corrects the lasso estimate along the resulting residual direction, which restores asymptotic normality. ```{r} dat0 <- list(x = data_inference$x, y = data_inference$y[, 1]) # beta_1 = 0 inf0 <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 1, NFOLDS = 5) r <- inf0$results[["1"]] c(estimate = r$coef_despar, se = r$se, p_value = r$p_value) r$confidence_band ``` Rows of `confidence_band` follow the order of `alpha`, which defaults to `c(0.01, 0.05, 0.10)`. The interval covers zero, as it should. With a non-zero coefficient the interval moves away from zero: ```{r} dat2 <- list(x = data_inference$x, y = data_inference$y[, 3]) # beta_1 = 0.2 inf2 <- hdcce_inference(dat2, obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 1, NFOLDS = 5) inf2$results[["1"]]$confidence_band ``` The reported $p$-value tests $H_0 \colon \beta_j = 0$ and is exactly dual to the band: $p \le \alpha$ if and only if zero lies outside the $(1-\alpha)$ interval. ```{r} cb <- r$confidence_band cbind(alpha = inf0$alpha, p_le_alpha = r$p_value <= inf0$alpha, zero_excluded = !(cb[, 1] <= 0 & 0 <= cb[, 2])) ``` Several coefficients can be treated in one call; the result carries one entry per index, named by it. Coefficients 1 and 5 are truly zero here, coefficient 2 is truly one. ```{r} inf_multi <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = c(1, 2, 5), NFOLDS = 5) t(sapply(inf_multi$results, function(z) c(estimate = z$coef_despar, se = z$se, p = z$p_value))) ``` The `HAC` argument selects the variance estimator: `1` assumes homoscedastic errors without serial correlation, `2` (the default) allows heteroscedasticity, and `3` allows serial correlation as well. The choice matters when $T$ is small, since the unit-wise variance estimates underlying `HAC = 2` rest on few degrees of freedom. ```{r} sapply(1:3, function(k) hdcce_inference(dat2, obs_N, obs_T, COEF_INDEX_VEC = 1, NFOLDS = 5, HAC = k)$results[["1"]]$se) ``` ## 3.2 Dictionary specification: significance test With a dictionary the null of interest is $H_0 \colon m_j = 0$, a functional null. There is no scalar to invert, so the output is a test rather than a confidence interval: `se` and `confidence_band` are `NULL`, and a `statistic`, `critical_values`, `p_value` and `profile` take their place. The statistic is a maximum, over a grid of locations $w$, of self-normalised correlations between the residuals under the null and kernel weights centred at $w$. Its critical values come from a Gaussian coupling simulated `B` times. **A covariate with no effect.** In the first response column $\beta_1 = 0$, so $m_1 \equiv 0$ and the null is true. ```{r} tst1 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]), obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 1, dictionaries = dict_inf, NFOLDS = 5, B = 500) q1 <- tst1$results[["1"]] c(statistic = q1$statistic, p_value = q1$p_value) q1$critical_values ``` **A covariate with an effect.** The second covariate has coefficient one in the same response, so the null is false. ```{r} tst2 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]), obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 2, dictionaries = dict_inf, NFOLDS = 5, B = 500) q2 <- tst2$results[["2"]] c(statistic = q2$statistic, p_value = q2$p_value) q2$critical_values ``` Because the statistic is a maximum it is informative to see where the evidence sits. The `profile` records $\Psi_{w,h}$ across the grid. ```{r profiles, fig.height = 3.2} op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.4, 1)) for (z in list(list(q = q1, main = "covariate 1 (no effect)"), list(q = q2, main = "covariate 2 (effect)"))) { pr <- z$q$profile plot(pr$w, abs(pr$Psi_w), type = "b", pch = 16, main = z$main, xlab = "location w", ylab = expression(group("|", Psi[list(w, h)], "|")), ylim = range(0, abs(pr$Psi_w), z$q$critical_values)) abline(h = z$q$critical_values[2], col = "red", lty = 2) } par(op) ``` For the first covariate the curve stays well below the 5% critical value at every location; for the second it crosses it, and the crossing identifies the region of the nodewise residuals in which the deviation is detected. The `profile` also reports `n_eff`, the number of nodewise residuals falling in each bump $[w-h, w+h]$. This is the diagnostic to check before trusting the test: when the smallest count is low the Gaussian approximation is unreliable, and the bandwidth should be increased or the region narrowed. ```{r} q2$profile ``` Note that the projection is rebuilt from $\bar{\boldsymbol{X}}_{(-j)}$ for each tested index, so the grid, $\widehat{K}$ and the residuals differ slightly between the two tests above. # 4. Two things worth checking **The factor number.** Everything downstream depends on the projection removing the factors, and that in turn requires the cross-sectional averages to identify the factor space — the usual CCE rank condition, which asks the mean loading matrix to have rank $K$. When it fails, $\widehat{K}$ is driven below the truth, a factor component survives in the residuals, the estimated error variance inflates and the procedures lose power. The scree plot of Section 2.1 is the quick diagnostic. **The occupancy of the grid.** For the test, `min(profile$n_eff)` says how much data supports the sparsest direction. ```{r} rbind(covariate_1 = c(min = min(q1$profile$n_eff), max = max(q1$profile$n_eff)), covariate_2 = c(min = min(q2$profile$n_eff), max = max(q2$profile$n_eff))) c(nT = obs_N * obs_T) ``` # References Rücker, M., Vogt, M., Linton, O. and Walsh, C. (2025). Estimation and inference in high-dimensional panel data models with interactive fixed effects. *Quantitative Economics* **16**(4), 1457–1509. [doi:10.3982/QE2308](https://doi.org/10.3982/QE2308) Rücker, M., Vogt, M. and Linton, O. (2026). High-dimensional panel data models with interactive fixed effects: beyond the linear case. [arXiv:2608.02055](https://arxiv.org/abs/2608.02055) ```{r cleanup, include=FALSE} options(old) ```