knitr::opts_chunk$set(echo = TRUE)
This package contains the implementation of the algorithm in
Varando G, Hansen NR (2020) Graphical continuous Lyapunov models
gclm
contains methods to estimate a sparse parametrization
of covariance matrix as solution of a continuous time Lyapunov
equation (CLE):
[ B\Sigma + \Sigma B^t + C = 0 ]
Solving the following (\ell_1) penalized loss minimization problem:
[ \arg\min L(\Sigma(B,C)) + \lambda \rho_1(B) + \lambda_C ||C - C_0||^2_F ]
subject to (B) stable and (C) diagonal, where (\rho_1(B)) is the (\ell_1) norm of the off-diagonal elements of (B) and (||C - C_0||^2_F ) is the squared frobenius norm of the difference between (C) and a fixed diagonal matrix ( C_0 ) (usually the identity).
## version on CRAN install.packages("gclm") ## development version from github devtools::install_github("gherardovarando/gclm")
library(gclm) ### define coefficient matrices B <- matrix(nrow = 4, c(-4, 2, 0, 0, 0, -3, 1, 0, 0, 0, -2, 0.5, 0, 0, 0, -1), byrow = TRUE) C <- diag(c(1,4,1,4)) ### solve continuous Lyapunov equation ### to obtain real covariance matrix Sigma <- clyap(B, C) ### obtain observations sample <- MASS::mvrnorm(n = 1000, mu = rep(0,4), Sigma = Sigma) ### Solve minimization res <- gclm(cov(sample), lambda = 0.4, lambdac = 0.01) res$B res$C
The CLE can be freely multiplied by a scalar and thus the (B,C) parametrization can be rescaled. As an example we can impose ( C_{11} = 1) as in the true ( C ) matrix, obtaining the estimators:
C1 <- res$C / res$C[1] B1 <- res$B / res$C[1] B1 C1
path <- gclm.path(cov(sample), lambdac = 0.01, lambdas = 10^seq(0, -3, length = 10)) t(sapply(path, function(res) c(lambda = res$lambda, npar = sum(res$B!=0), fp = sum(res$B!=0 & B==0), tp = sum(res$B!=0 & B!=0) , fn = sum(res$B==0 & B!=0), tn = sum(res$B==0 & B==0), errs = sum(res$B!=0 & B==0) + sum(res$B==0 & B!=0))))
lyapunov
package (https://github.com/gragusa/lyapunov).Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.