# In gherardovarando/clggm: Continuous Lyapunov Gaussian Graphical Models

```knitr::opts_chunk\$set(echo = TRUE)
```

# clggm

`clggm` contatins methods to estimate the coefficients of the stochastic process (Ornstein–Uhlenbeck) \$dX_t = B(X_t - a)dt + DdW_t\$ from observations of the invariant distribution.

## Installation

```devtools::install_github("gherardovarando/clggm")
```

### Dependencies

• `MASS`
• `stats`
• `testthat` (only in dev)

## Usage

```library(clggm)

### define coefficient matrices
B <- matrix(nrow = 4, c(-4, 1,   0,  0,
0, -3,  1,  0,
0,  0, -2,  1,
0,  0,  0, -1), byrow = TRUE)
D <- diag(c(1,2,1,2))
C <- D %*% t(D)

### solve continuous Lyapunov equation
### to obtain covariance of invariant
### distribution
Sigma <- clyap(B, C)

### obtain observations from the invariant distribution
sample <- MASS::mvrnorm(n = 10000, mu = rep(0,4),Sigma =  Sigma)
```

### Estimate \$B\$ knowing \$C\$

```B0 <- - 0.5 * C %*% solve(cov(sample))
### penalized maximum-likelihood solved with proximal gradient
results <- proxgradllB(cov(sample), B= B0, C = C, eps = 1e-12,
maxIter = 1000, lambda = 0.05, job = 1)
results

### penalized least square 0.5*||S(B) - Sigma||_2^2 + lambda * ||B||_1,off
results2 <- proxgradlsB(cov(sample), B= B0, C = C, eps = 1e-12,
maxIter = 1000, lambda = 0.01, job = 0)
results2\$B
```

#### Solutions along a regularization path

```results <- llBpath(Sigma, eps = 1e-12, maxIter = 1000, job = 11)
t(sapply(results, 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))))
```

### Estimate \$C\$ knowing \$B\$

```results <- graddsllc(cov(sample), B, C = diag(4), C0 = diag(4),
eps = 1e-12, maxIter = 1000, lambda = 0.0005)
cbind(C = diag(C), Cest = diag(results\$C))
```

### Estimate \$B\$ and \$C\$

```B0 <- - 0.5 * diag(4) %*% solve(cov(sample))
results <- pnllbc(cov(sample), B0, C = diag(4), C0 = diag(4),
eps = 1e-14, maxIter = 1000, intitr = 1,
lambda = 0.2, lambdac = 0.005, job = 0)
## B estimated
results\$B

## C estimated
cbind(C = diag(C), Cest = diag(results\$C))
```

## Related code

• Some code is from the `lyapunov` package (https://github.com/gragusa/lyapunov).

• Using FORTRAN code of Algorithm 705 from the Collected Algorithms from ACM, TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL. 18, NO. 2, PP. 232-238.75.

gherardovarando/clggm documentation built on Feb. 14, 2020, 12:31 p.m.