invGlasso <- function(x){
if (all(eigen(x)$values > sqrt(.Machine$double.eps))){
Xinv <- solve(x)
} else {
Xglas <- glasso(x,0.05,penalize.diagonal=FALSE)
Xinv <- Xglas$wi
}
Xinv
}
generate_lambdas <- function(
X,
Y,
nLambda_kappa = 10,
nLambda_beta = 10,
lambda_min_kappa = 0.05,
lambda_min_beta = 0.05,
penalize.diagonal=TRUE,
version0.1.4 = FALSE
){
N <- nrow(Y)
P <- ncol(Y)
#### Lambda sequence for Kappa:
corY <- cov2cor(t(Y)%*%Y/nrow(Y))
if (version0.1.4){
lam_K_max = max(abs(corY))
} else {
lam_K_max = max(abs(corY[upper.tri(corY)]))
}
lam_K_min = lambda_min_kappa*lam_K_max
lam_K = exp(seq(log(lam_K_max), log(lam_K_min), length = nLambda_kappa))
#### Lambda sequence for Beta
# Initial estimate for Kappa:
# Yinv <- invGlasso(t(Y) %*% Y / N)
# Xinv <- invGlasso(t(X) %*% X / N)
# beta <- t(Y) %*% X %*% Xinv
# S <- 1/(nrow(Y)) * (
# t(Y) %*% Y -
# t(Y) %*% X %*% t(beta) -
# beta %*% t(X) %*% Y +
# beta %*% t(X) %*% X %*% t(beta)
# )
# S <- (S + t(S)) / 2
# if (any(eigen(S)$value < -sqrt(.Machine$double.eps))) stop("Residual covariances not postive definite")
# kappa <- invGlasso(S)
# kappa <- (kappa + t(kappa)) / 2
# lam_B_max = max(abs((1/N)*t(X)%*%Y%*%Yinv))
Yinv <- invGlasso(t(Y) %*% Y)
lam_B_max = max(abs(t(X)%*%Y%*%Yinv))
lam_B_min = lambda_min_beta*lam_B_max
lam_B = exp(seq(log(lam_B_max), log(lam_B_min), length = nLambda_beta))
return(list(lambda_kappa = lam_K, lambda_beta = lam_B))
}
### BACKWARD COMPETABILITY ###
SparseTSCGM_lambdas <-
function(X, Y, nlambda = 100, lambda.min.ratio = 0.01){
# From SparseTSCGM package:
lambda.seq <- function(SS, SA,nlambda)
{
if (length(nlambda)==1) nlambda <- rep(nlambda,2)
# lambda.min.ratio=0.1
d= dim(SS)[2]
lambda.max1 = max(max(SS-diag(d)),-min(SS-diag(d)))
lambda.min1 = lambda.min.ratio*lambda.max1
lambda1 = exp(seq(log(lambda.max1), log(lambda.min1), length = nlambda[1]))
lambda.min.ratio2=0.15
lambda.max2 = max(max(SA),-min(SA))
lambda.min2 = lambda.min.ratio2*lambda.max2
lambda2 = exp(seq(log(lambda.max2), log(lambda.min2), length = nlambda[2]))
return(list(lambda1=lambda1, lambda2=lambda2))
}
T <- dim(Y)[1]
p <- dim(X)[2]
n <- 1
q <- dim(Y)[2]
xtyi <- array(NA, c(p,q,n))
xtxi <- array(NA, c(p,p,n))
ytyi <- array(NA, c(q,q,n))
XX <- X
YY <- Y
XX2 <- X^2
YY2 <- Y^2
xtyi <- crossprod(XX,YY)
xtxi <- crossprod(XX)
ytyi <- crossprod(YY)
xty=apply(xtyi, c(1,2), sum)
xtx=apply(xtxi, c(1,2), sum)
yty=apply(ytyi, c(1,2), sum)
xtxt=apply(xtxi, c(1,2), sum)/(n*T)
xtx2=(n*T)*colMeans(apply(XX2, c(1,2), sum))
yty2=(n*T)*colMeans(apply(YY2, c(1,2), sum))
SX <- xtx/(n*T)
mSX <- glasso(SX,0.05,penalize.diagonal=FALSE)
SX <- xtx/(n*T)
mSX <- glasso(SX,0.05,penalize.diagonal=FALSE)
SXi <- mSX$wi
SS =(yty)/(n*T)
SS = cov2cor(SS)
SAs = xty/(n*T)
SA = t(SAs) %*% SXi
lambda <- lambda.seq(SS=SS,SA=SA, nlambda=nlambda)
lam1 <- round(lambda$lambda1,3)
lam2 <- round(lambda$lambda2,3)
lam2 <- round(lam2/max(lam2),3)
return(list(lambda_kappa = lam1, lambda_beta = lam2))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.