--- title: "Inference for Sparse Spectral Precision Matrices with `cxreg`" author: - Younghoon Kim - Navonil Deb - Sumanta Basu date: "`r format(Sys.time(), '%B %d, %Y')`" output: pdf_document vignette: > %\VignetteIndexEntry{Inference for Sparse Spectral Precision Matrices with cxreg} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} link-citations: true bibliography: cxreg_refs.bib --- ```{r include=FALSE} hook_output <- knitr::knit_hooks$get("output") knitr::knit_hooks$set(output = function(x, options) { if (!is.null(n <- options$out.lines)) { x <- xfun::split_lines(x) if (length(x) > n) { x <- c(head(x, n), "....\n") } x <- paste(x, collapse = "\n") } hook_output(x, options) }) ``` ## Introduction This vignette describes the inference pipeline for sparse spectral precision matrices (SPMs) implemented in \texttt{cxreg}. The methodology is developed in [-@deb2026inference], referred to as DEB26 throughout. We focus on practical usage of the functions. The spectral precision matrix $\Theta(\omega)$ at frequency $\omega\in(-\pi,\pi]$ is the inverse of the spectral density matrix $f(\omega)$ at the same given frequency $\omega$, and plays the role of a precision matrix in the spectral domain: its zero and nonzero off-diagonal entries encode the absence or presence of frequency-specific conditional dependencies among time series components. The inference pipeline in \texttt{cxreg} consists of five steps: \begin{enumerate} \item \textbf{Data-driven bandwidth selection}: \texttt{select\_m()} selects the smoothing bandwidth $m$ via a generalised cross-validation (GCV) criterion. \item \textbf{Spectral density estimation}: \texttt{dft.all()} and \texttt{fhat\_at()} compute the discrete Fourier transforms (DFTs) and the smoothed periodogram $\hat{f}(\omega_j)$. \item \textbf{Debiased estimation}: \texttt{decglasso()} constructs the debiased complex graphical lasso (deCGLASSO) estimator $\tilde{\Theta}_j$ from the sparse complex graphical lasso (CGLASSO). \item \textbf{Variance estimation}: \texttt{var.cov()} estimates the asymptotic variance and pseudovariance of $\tilde{\Theta}_j$ using either a plug-in or heteroscedasticity and autocorrelation consistent (HAC) estimator. \item \textbf{Inference}: \texttt{spec.test()} constructs entry-wise confidence ellipsoids and test statistics; \texttt{spec.fdr()} applies multiple hypothesis testing with controlled false discovery rate (FDR). \end{enumerate} ## Background ### Spectral density and the DFT For a $p$-dimensional weakly stationary time series $\{X_t\}_{t=1}^n$, the discrete Fourier transform (DFT) at a Fourier frequency $\omega_j=2\pi j/n$ is $$ d_j := d(\omega_j) = \frac{1}{\sqrt{n}}\sum_{t=1}^n X_t \exp(-i t \omega_j). $$ Asymptotically, $d_j \sim \mathcal{N}_{\mathbb{C}}(0, f(\omega_j))$. The smoothed periodogram at $\omega_j$ with half-bandwidth $m$ is $$ \hat{f}(\omega_j) = \frac{1}{2\pi(2m+1)} \sum_{k \in W_m(j)} d_k d_k^*, $$ where $W_m(j)=\{k:|k-j|\leq m\}$ is the set of $2m+1$ neighboring frequencies around $\omega_j$. ### The CGLASSO and deCGLASSO estimators The CGLASSO estimator $\hat{\Theta}_j$ is the minimizer of the optimization problem $$ \hat{\Theta} = \arg\min_{\Theta \in \mathcal{H}_{++}^p} \left\{ -\log\det\Theta + \mathrm{trace}(\hat{f}_j\Theta) + \lambda \|\Theta^-\|_1 \right\}, $$ where $\|\Theta^-\|_1=\sum_{a \neq b} |\Theta_{a,b}|$ penalizes off-diagonal entries of $\Theta$. In high dimensions, $\hat{\Theta}_j$ has a regularization bias due to the penalty term. The deCGLASSO estimator removes this bias through the computation [-@deb2026inference]: $$ \tilde{\Theta}_j = 2\hat{\Theta}_j - \hat{\Theta}_j\hat{f}_j\hat{\Theta}_j. $$ Under suitable regularity conditions, [-@deb2026inference] showed that $\sqrt{2m+1}\,(\tilde{\Theta}_j-\Theta_j)_{ab}$ is asymptotically bivariate normal with a covariance matrix constructed by using $\Theta_j$. ### Step 1 — Generate a multivariate white noise process We illustrate the pipeline using a multivariate Gaussian white noise process $X_t\sim\mathrm{WN}(0,\Sigma)$, $t=1,\ldots,n$, with a sparse tridiagonal precision matrix $\Theta=\Sigma^{-1}$. Note that the spectral density of white noise process is constant across all frequencies: $$ f(\omega) = \frac{1}{2\pi}\Sigma,\quad \omega\in(-\pi,\pi], $$ so the true SPM is also constant over frequencies: $$ \Theta(\omega) = f(\omega)^{-1} = 2\pi \Sigma^{-1}. $$ ```{r} library(cxreg) library(mvtnorm) ``` ```{r} set.seed(1010) p <- 10 # number of variables n <- 500 # sample size # True precision matrix: tridiagonal C <- diag(0.7, p) C[row(C) == col(C) + 1] <- 0.3 C[row(C) == col(C) - 1] <- 0.3 Sigma <- solve(C) # True SPM: constant across all frequencies Theta_true <- 2*pi*C # Generate white noise observations X <- rmvnorm(n = n, mean = rep(0,p), sigma=Sigma) cat("Data dimensions:", dim(X), "\n") cat("True SPM (1:4, 1:4):\n") print(round(Re(Theta_true[1:4, 1:4]), 3)) ``` ### Step 2 — Bandwidth selection via GCV \texttt{select\_m()} searches over a grid of half-bandwidths and selects the one that minimizes a GCV criterion based on the deviance of the diagonal periodogram [-@ombao2001gcv]. The GCV criterion provides a principled data-driven choice of the half-bandwidth $m$. ```{r} bw_sel <- select_m(X, verbose = FALSE) m <- bw_sel$m_opt cat("Selected bandwidth: m =", m, "(full bandwidth:", 2*m+1,")\n") bw_sel$gcv_table ``` ### Step 3 — Compute the DFT and smoothed periodogram We target $j=\lfloor n/4 \rfloor$, corresponding to $\omega\approx\pi/2$. ```{r} j <- floor(n / 4) dft <- dft.all(X) # full DFT matrix: n x p fhat <- fhat_at(dft, j = j, m = m) # smoothed periodogram: p x p cat("Frequency index j =", j, "(omega =", round(2*pi*j/n, 3), ")\n") # For white noise: fhat should be close to Sigma / (2*pi) cat("Max deviation of Re(fhat) from Sigma/(2pi):", round(max(abs(Re(fhat) - Sigma / (2*pi))), 4), "\n") ``` ### Step 4 — Fit CGLASSO and compute the deCGLASSO estimator We fit the CGLASSO path using the EBIC stopping criterion, then apply one-step debiasing: ```{r} # Fit CGLASSO fit <- cglasso(S = fhat, m = m, type = "II") cat("EBIC-selected lambda index:", fit$min_index, "\n") cat("Lambda selected:", round(fit$lambda_grid[fit$min_index], 4), "\n") ``` ```{r} # One-step debiasing to produce deCGLASSO deb <- decglasso(object = fit, fhat = fhat) cat("Class of deCGLASSO output:", class(deb), "\n") Theta_tilde <- deb$Theta_tilde ``` ### Step 5 — Estimate the asymptotic variance ```{r} vc_plug <- var.cov(Theta = Theta_tilde, X = X, j = j, m = m, type = "plug-in") vc_hac <- var.cov(Theta = Theta_tilde, X = X, j = j, m = m, type = "HAC") cat("Fields:", names(vc_plug), "\n") ``` The output contains: - \texttt{Vre}: estimated variance of $\mathrm{Re}(\tilde{\Theta}_{ab})$, a $p \times p$ matrix. - \texttt{Vim}: estimated variance of $\mathrm{Im}(\tilde{\Theta}_{ab})$, a $p \times p$ matrix. - \texttt{Cov}: estimated covariance of $\mathrm{Re}$ and $\mathrm{Im}$ parts, a $p \times p$ matrix. - \texttt{V}: complex variance $\sigma^2$, a $p \times p$ complex matrix. - \texttt{PV}: pseudovariance $\delta^2$, a $p \times p$ complex matrix. These satisfy $\texttt{Vre} = \frac{1}{2}\mathrm{Re}(V + PV)$, $\texttt{Vim} = \frac{1}{2}\mathrm{Re}(V - PV)$, and $\texttt{Cov} = \frac{1}{2}\mathrm{Im}(PV)$. Compare the two estimators on the $(1,2)$ entry: ```{r} cat("Entry (1,2) variance (Re part):\n") cat("Plug-in Vre:", round(vc_plug$Vre[1,2], 6), "\n") cat("HAC Vre:", round(vc_hac$Vre[1,2], 6), "\n") ``` ### Step 6 — Entry-wise tests and confidence regions \texttt{spec.test()} computes for each off-diagonal entry $(a,b)$: - $Z$-statistics for $\mathrm{Re}(\tilde\Theta_{ab})$ and $\mathrm{Im}(\tilde\Theta_{ab})$. - The Mahalanobis chi-squared statistic $(2m+1)\hat{D}_{j,m}^{(a,b)}(0) \sim \chi^2_2$ under $H_0:\Theta_{ab}=0$. - Half-widths of the marginal confidence intervals (CI) for the real and imaginary parts, and areas of the joint 95% confidence region (ellipsoids) are produced. ```{r} alpha <- 0.2 st_plug <- spec.test(Est = Theta_tilde, varcov = vc_plug, m = m, alpha = alpha) st_hac <- spec.test(Est = Theta_tilde, varcov = vc_hac, m = m, alpha = alpha) cat("Chi-squared statistic at (1,2) [true edge]:\n") cat("Plug-in:", round(st_plug$Chi_sq[1,2], 3), "\n") cat("HAC:", round(st_hac$Chi_sq[1,2], 3), "\n") cat("Chi-sq(2, 0.95) critical value:", round(qchisq(0.95, 2), 3), "\n") cat("\n Chi-squared statistic at (1,3) [true null edge]:\n") cat("Plug-in:", round(st_plug$Chi_sq[1,3], 3), "\n") cat("HAC:", round(st_hac$Chi_sq[1,3], 3), "\n") ``` The CI half-widths for the real and imaginary parts at the $(1,2)$ entry: ```{r} cat("80% CI half-width at (1,2):\n") cat("Re part (plug-in):", round(st_plug$wing_re[1,2], 4), "\n") cat("Im part (plug-in):", round(st_plug$wing_im[1,2], 4), "\n") ``` For white noise process, the imaginary part of $\Theta(\omega)$ is zero, so we expect the imaginary CI to contain zero and the real CI to contain $2\pi (\Sigma^{-1})_{ab}$: ```{r} # Check if true value falls inside the 95% CI true_re_12 <- Re(Theta_true[1,2]) est_re_12 <- Re(Theta_tilde[1,2]) hw_re_12 <- st_plug$wing_re[1,2] cat("True Re(Theta_12):", round(true_re_12, 4), "\n") cat("80% CI: [",round(est_re_12 - hw_re_12, 4), ",", round(est_re_12 + hw_re_12, 4),"]\n") cat("True value inside CI:", abs(est_re_12 - true_re_12) <= hw_re_12, "\n") ``` ### Step 7 — FDR-controlled support recovery \texttt{spec.fdr()} applies the thresholding procedure of [-@liu2013gaussian] to control the FDR across all $p(p-1)/2$ off-diagonal entries simultaneously. ```{r} fdr_plug <- spec.fdr(Chi_sq = st_plug$Chi_sq, alpha = alpha, diag = FALSE) fdr_hac <- spec.fdr(Chi_sq = st_hac$Chi_sq, alpha = alpha, diag = FALSE) cat("Number of significant edges (plug-in):", sum(fdr_plug$Decision[upper.tri(fdr_plug$Decision)]), "\n") cat("Number of significant edges (HAC):", sum(fdr_hac$Decision[upper.tri(fdr_hac$Decision)]), "\n") cat("FDR threshold tau (plug-in):", round(fdr_plug$tau, 3), "\n") ``` For the white noise process, the true edge set of $\Theta(\omega) = 2\pi\Sigma^{-1}$ is the same at every frequency. Compare the discoveries to the truth: ```{r} # True support: nonzero off-diagonal entries of C (same as Theta_true up to scale) true_support <- (C!=0)*1L diag(true_support) <- 0 # Confusion matrix (plug-in) dec <- fdr_plug$Decision TP <- sum(dec == 1 & true_support == 1 & upper.tri(dec)) FP <- sum(dec == 1 & true_support == 0 & upper.tri(dec)) FN <- sum(dec == 0 & true_support == 1 & upper.tri(dec)) cat("True edges in upper triangle:", sum(true_support[upper.tri(true_support)]), "\n") cat("True Positives:", TP, "\n") cat("False Positives:", FP, "\n") cat("False Negatives:", FN, "\n") cat("Empirical FDR:", round(FP / max(TP + FP, 1), 3), "\n") cat("Power:", round(TP / max(TP + FN, 1), 3), "\n") ``` ## Summary of functions | Step | Function | Description | |------|----------|-------------| | 1 | `select_m()` | GCV bandwidth selection | | 2 | `dft.all()` | Full DFT matrix | | 2 | `fhat_at()` | Smoothed periodogram at frequency $j$ | | 3 | `cglasso()` | CGLASSO estimator path | | 4 | `decglasso()` | One-step debiasing | | 5 | `var.cov()` | Asymptotic variance estimation (plug-in or HAC) | | 6 | `spec.test()` | Chi-squared statistics and CI half-widths | | 7 | `spec.fdr()` | FDR-controlled support recovery | ## References