README.md

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

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 
## $N
## [1] 4
## 
## $Sigma
##            [,1]       [,2]       [,3]       [,4]
## [1,] 0.15601606 0.08321924 0.01558787 0.01211664
## [2,] 0.08321924 0.73264757 0.08683278 0.07235061
## [3,] 0.01558787 0.08683278 0.55492104 0.60110874
## [4,] 0.01211664 0.07235061 0.60110874 1.97583949
## 
## $B
##          [,1]       [,2]       [,3]       [,4]
## [1,] -3.58569  0.7140809  0.0000000  0.0000000
## [2,]  0.00000 -2.7795702  0.4197185  0.0000000
## [3,]  0.00000  0.1439541 -1.7630976  0.7750345
## [4,]  0.00000  0.0000000  0.3400480 -1.1156806
## 
## $C
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    4    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    4
## 
## $lambda
## [1] 0.05
## 
## $diff
## [1] 3.342472e-15
## 
## $objective
## [1] 1.461488
## 
## $iter
## [1] 95
## 
## $job
## [1] 1
### 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
##           [,1]        [,2]       [,3]          [,4]
## [1,] -3.190091  0.56017938  0.0000000  2.298023e-15
## [2,]  0.000000 -2.78017258  0.4427526  1.578973e-01
## [3,]  0.000000  0.07089838 -1.6362058  7.768397e-01
## [4,]  0.000000  0.00000000  0.3521275 -1.118297e+00

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))))
##       lambda npar fp tp fn tn errs
##  [1,]    0.1    8  1  7  0  8    1
##  [2,]    0.2    7  0  7  0  9    0
##  [3,]    0.3    5  0  5  2  9    2
##  [4,]    0.4    5  0  5  2  9    2
##  [5,]    0.5    5  0  5  2  9    2
##  [6,]    0.6    5  0  5  2  9    2
##  [7,]    0.7    5  0  5  2  9    2
##  [8,]    0.8    5  0  5  2  9    2
##  [9,]    0.9    5  0  5  2  9    2
## [10,]    1.0    5  0  5  2  9    2

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))
##      C      Cest
## [1,] 1 1.0290756
## [2,] 4 4.0162866
## [3,] 1 0.9962571
## [4,] 4 3.9106229

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.05, lambdac = 0.005, job = 0)
## B estimated
results$B
##           [,1]       [,2]        [,3]       [,4]
## [1,] -3.508704  0.5414103  0.00000000  0.0000000
## [2,]  0.000000 -0.6687548  0.07730917  0.0404080
## [3,]  0.000000  0.1426657 -1.20064232  0.4779548
## [4,]  0.000000  0.0000000  0.00000000 -0.2974951
## C estimated
cbind(C = diag(C), Cest = diag(results$C))
##      C      Cest
## [1,] 1 0.9690831
## [2,] 4 0.9396312
## [3,] 1 0.6999000
## [4,] 4 1.1927048

Related code



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