BB_estimator <- function(X, sim, interval, method){
  N <- nrow(X)
  dirichlet_sample <- matrix(rexp(N * sim, 1) , ncol = N, byrow = TRUE)
  dirichlet_sample <- dirichlet_sample / rowSums(dirichlet_sample)
  off_diag <- list()
  diag_l <- list()

  if(method == "ridge"){

    for( i in 1:sim) {
      print(i)
      center <- colSums(X * dirichlet_sample[i,])
      x <- sqrt(dirichlet_sample[i,]) * sweep(X, 2, center, check.margin = FALSE)

      shrinkage.lambda <- lambda.TargetD(x = X)
      OPT <-  ridge.lambda <- shrinkage.lambda / (1.0 - shrinkage.lambda)
      mat_temp2 <- ridgeP(crossprod(x), lambda = OPT , type = "ArchII")

      off_diag[[i]] <- mat_temp2[upper.tri(mat_temp2)]
      diag_l[[i]] <- diag(mat_temp2)

    }
  }

  if(method == "lw"){
    for(i in 1:sim){
      n=nrow(X)
      p = ncol(X)
      print(i)
      # de-mean returns
      #Y = scale(Y, scale = FALSE)
      center <- colSums(X * dirichlet_sample[i,])
      x <- sqrt(dirichlet_sample[i,]) * sweep(X, 2, center, check.margin = FALSE)
      # compute S covariance matrix
      S = crossprod(X)/n
      S_w <- crossprod(x)
      # compute prior

      meanvar = mean(diag(S))

      prior = meanvar*diag(1,p)


      # what we call b
      T=X^2
      phiMat = (crossprod(T)/n) - S^2
      phi = sum(phiMat)

      # what we call c
      gamma = sum(abs(S-prior)^2)

      # compute shrinkage constant
      kappa = phi/gamma
      shrinkage = max(0,min(1,kappa/n))
      Sigma=shrinkage*prior+(1-shrinkage)*S_w
      Sigma <- solve(Sigma)
      off_diag[[i]] <- Sigma[upper.tri(Sigma)]
      diag_l[[i]] <- diag(Sigma)
    }
  }
if(method == "ml"){
  for(i in 1:sim){
  center <- colSums(X * dirichlet_sample[i,])
  x <- sqrt(dirichlet_sample[i,]) * sweep(X, 2, center, check.margin = FALSE)
  # compute S covariance matrix
  ml <- crossprod(x)
  ml <- solve(ml)
  off_diag[[i]] <- ml[upper.tri(ml)]
  diag_l[[i]] <- diag(ml)
  }

}
  mlt_dat <- reshape2::melt(diag_l)
  mlt_dat$position <- 1:ncol(X)
  t <- mlt_dat  %>% group_by(position) %>% summarise(mean(value))
  d <- t$`mean(value)`

  mlt <- reshape2::melt(off_diag)
  mlt$position <- 1:((ncol(X) * (ncol(X) - 1)) / 2)
  results  <- mlt %>%
    group_by(position) %>%
    summarise(mean(value),
              low = quantile(value, (1 - interval) / 2)[1],
              up = quantile(value, 1 -( (1 - interval) / 2))[1])

  mat_point <- matrix(0, ncol(X), ncol(X))
  mat_sig <- matrix(0, ncol(X), ncol(X))
  mat_point[upper.tri(mat_point)]   <-  results$`mean(value)`
  mat_point <- as.matrix(Matrix::forceSymmetric(mat_point))
  sig <- ifelse(results$low < 0 & results$up > 0, 0, 1)
  mat_sig[upper.tri(mat_sig)]   <-  sig
  mat_sig <- as.matrix(Matrix::forceSymmetric(mat_sig))
  mat_point_sparse <- mat_sig * mat_point
  diag(mat_point_sparse) <- d
  list(mat_point = mat_point, mat_sig = mat_sig, mat_point_sparse = mat_point_sparse)
}