covgpenal: covgpenal

In VecDep: Measuring Copula-Based Dependence Between Random Vectors

covgpenal

Description

This function computes the empirical penalized Gaussian copula covariance matrix with the Gaussian log-likelihood plus a lasso-type penalty as objective function. Model selection is done by choosing omega such that BIC is maximal.

Usage

covgpenal(
  S,
  n,
  omegas,
  derpenal = function(t, omega) {
     derSCAD(t, omega, 3.7)
 },
  nsteps = 1
)

Arguments

S	The sample matrix of normal scores covariances.
n	The sample size.
omegas	The candidate values for the tuning parameter in [0,\infty).
derpenal	The derivative of the penalty function to be used (default = scad with parameter a = 3.7).
nsteps	The number of weighted covariance graphical lasso iterations (default = 1).

Details

The aim is to solve/compute

\widehat{\boldsymbol{\Sigma}}_{\text{LT},n} \in \text{arg min}_{\boldsymbol{\Sigma} > 0} \left \{ \ln \left | \boldsymbol{\Sigma} \right | + \text{tr} \left (\boldsymbol{\Sigma}^{-1} \widehat{\boldsymbol{\Sigma}}_{n} \right ) + P_{\text{LT}}\left (\boldsymbol{\Sigma},\omega_{n} \right ) \right \},

where the penalty function P_{\text{LT}} is of lasso-type:

P_{\text{LT}} \left (\boldsymbol{\Sigma},\omega_{n} \right ) = \sum_{ij} p_{\omega_{n}} \left (\Delta_{ij} \left |\sigma_{ij} \right | \right ),

for a certain penalty function p_{\omega_{n}} with penalty parameter \omega_{n}, and \sigma_{ij} the (i,j)'th entry of \boldsymbol{\Sigma} with \Delta_{ij} = 1 if i \neq j and \Delta_{ij} = 0 if i = j (in order to not shrink the variances). The matrix \widehat{\boldsymbol{\Sigma}}_{n} is the matrix of sample normal scores covariances.

In case p_{\omega_{n}}(t) = \omega_{n} t is the lasso penalty, the implementation for the (weighted) covariance graphical lasso is available in the R package ‘covglasso’ (see the manual for further explanations). For general penalty functions, we perform a local linear approximation to the penalty function and iteratively do (nsteps, default = 1) weighted covariance graphical lasso optimizations.

The default for the penalty function is the scad (derpenal = derivative of scad penalty), i.e.,

p_{\omega_{n},\text{scad}}^{\prime}(t) = \omega_{n} \left [1 \left (t \leq \omega_{n} \right ) + \frac{\max \left (a \omega_{n} - t,0 \right )}{\omega_{n} (a-1)} 1 \left (t > \omega_{n} \right ) \right ],

with a = 3.7 by default.

For tuning \omega_{n}, we maximize (over a grid of candidate values) the BIC criterion

\text{BIC} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right ) = -n \left [\ln \left |\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right | + \text{tr} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}}^{-1} \widehat{\boldsymbol{\Sigma}}_{n} \right ) \right ] - \ln(n) \text{df} \left (\widehat{\boldsymbol{\Sigma}}_{\omega_{n}} \right ),

where \widehat{\boldsymbol{\Sigma}}_{\omega_{n}} is the estimated candidate covariance matrix using \omega_{n} and df (degrees of freedom) equals the number of non-zero entries in \widehat{\boldsymbol{\Sigma}}_{\omega_{n}}, not taking the elements under the diagonal into account.

Value

A list with elements "est" containing the (lasso-type) penalized matrix of sample normal scores rank correlations (output as provided by the function “covglasso.R”), and "omega" containing the optimal tuning parameter.

References

De Keyser, S. & Gijbels, I. (2024). High-dimensional copula-based Wasserstein dependence. doi: https://doi.org/10.48550/arXiv.2404.07141.

Fop, M. (2021). covglasso: sparse covariance matrix estimation, R package version 1.0.3. url: https://CRAN.R-project.org/package=covglasso.

Wang, H. (2014). Coordinate descent algorithm for covariance graphical lasso. Statistics and Computing 24:521-529. doi: https://doi.org/10.1007/s11222-013-9385-5.

See Also

grouplasso for group lasso estimation of the normal scores rank correlation matrix.

Examples

q = 10
dim = c(5,5)
n = 100

# AR(1) correlation matrix with correlation 0.5
R = 0.5^(abs(matrix(1:q-1,nrow = q, ncol = q, byrow = TRUE) - (1:q-1)))

# Sparsity on off-diagonal blocks
R0 = createR0(R,dim)

# Sample from multivariate normal distribution
sample = mvtnorm::rmvnorm(n,rep(0,q),R0,method = "chol")

# Normal scores
scores = matrix(0,n,q)
for(j in 1:q){scores[,j] = qnorm((n/(n+1)) * ecdf(sample[,j])(sample[,j]))}

# Sample matrix of normal scores covariances
Sigma_est = cov(scores) * ((n-1)/n)

# Candidate tuning parameters
omega = seq(0.01, 0.6, length = 50)

Sigma_est_penal = covgpenal(Sigma_est,n,omega)

VecDep documentation built on April 4, 2025, 5:14 a.m.