R/BvarSV.R

In hoanguc3m/fatBVARS: Bayesian VAR with Stochastic volatility and fat tails

#' Bayesian inference of VAR model with RW-SV
#'
#' Bayesian inference of VAR model with RW-SV
#' \deqn{y_t = B x_t + SQRT(w_t) A^(-1) H_t eps_t}
#' @param y The input data as matrix T x K where T is the number of observations and K is the number of variables
#' @param K The number of variables in BVAR model.
#' @param p The number of lags in BVAR model.
#' @param dist The variable specifies the BVAR error distribution. It should be one of
#' c("Gaussian","Student","Skew.Student","Skew.Student", "MT","Skew.MT","MST").
#' @param y0 The number of observations
#' @param prior The prior specification of BVAR.
#' @param inits The initial values of BVAR.
#' @return A coda object of the posterior samples.
#' @export
#' @examples
#' \dontrun{
#' datagen <- sim.VAR.SV(dist="Gaussian")
#' y <- datagen$y
#' prior <- get_prior(y, p = 2, dist="Gaussian", SV = T)
#' inits <- get_init(prior)
#' Chain1 <- BVAR.SV(y, K = 5, p = 2, dist = "Gaussian", y0 = NULL, prior = prior, inits = inits)
#' plot(Chain1)
#' }
BVAR.SV <- function(y, K, p, dist, y0 = NULL, prior = NULL, inits = NULL){
  if (!(dist %in% c("Gaussian","Student","Skew.Student",
                    "MT","MST",
                    "OT","OST",
                    "dynSkew.Student", "dynMST", "dynOST") ))
    stop("dist is not implemented.")
  if (prior$SV == TRUE){
    Start = Sys.time()
    if (dist == "Gaussian") Chain <- BVAR.Gaussian.SV(y, K, p, y0, prior, inits)
    if (dist == "Student") Chain <- BVAR.Student.SV(y, K, p, y0, prior, inits)
    if (dist == "Skew.Student") Chain <- BVAR.Skew.Student.SV(y, K, p, y0, prior, inits)
    if (dist == "MT") Chain <- BVAR.MT.SV(y, K, p, y0, prior, inits)
    if (dist == "MST") Chain <- BVAR.MST.SV(y, K, p, y0, prior, inits)
    if (dist == "OT") Chain <- BVAR.OT.SV(y, K, p, y0, prior, inits)
    if (dist == "OST") Chain <- BVAR.OST.SV(y, K, p, y0, prior, inits)

    if (dist == "dynSkew.Student") Chain <- BVAR.dynSkew.Student.SV(y, K, p, y0, prior, inits)
    if (dist == "dynMST") Chain <- BVAR.dynMST.SV(y, K, p, y0, prior, inits)
    if (dist == "dynOST") Chain <- BVAR.dynOST.SV(y, K, p, y0, prior, inits)
    elapsedTime = Sys.time() - Start
    print(elapsedTime)
    out <- list(mcmc = Chain,
                y = y,
                y0 = y0,
                K = K,
                p = p,
                dist = dist,
                prior = prior,
                inits = inits,
                esttime = elapsedTime)
    class(out) <- c("fatBVARSV")
    return(out)
  } else {
    warning("prior$SV is TRUE")
  }
}
#' @export
BVAR.Gaussian.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "Gaussian", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior


  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)
  # new precompute
  theta.prior.precmean <- V_b_prior_inv %*% b_prior

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # Output
  mcmc <- matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K*t_max,
                     ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # Sample B, use reduce sum here
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i])) %*% A)
    #   b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% yt[,i])
    # }

    # V_b_post_inv <- V_b_prior_inv +
    #   Reduce( f = "+",
    #           x = lapply(1:t_max, function(i) kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A)),
    #           accumulate = FALSE)
    # b_post <- Reduce(f = "+",
    #                  x = lapply(1:t_max, function(i) kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A) %*% yt[,i])),
    #                  accumulate = FALSE)
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)

    # Sample B
    wt <- as.vector(exp(-h/2))
    y.tilde <- as.vector( A %*% yt ) * wt
    x.tilde <- kronecker( t(xt), A ) * wt
    theta.prec.chol <- chol( V_b_prior_inv + crossprod(x.tilde) )
    b_sample <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*m) )
    B <- matrix(b_sample,K,m)

    # Sample vol
    ytilde <- A%*% (yt - B %*% xt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h


    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*%xt)
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    if (K > 1) {
      for (i in c(2:K)){
        id_end <- i*(i-1)/2
        id_start <- id_end - i + 2
        a_sub <- a_prior[id_start:id_end]
        V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
        a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                     xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                     a_sub = a_sub,
                                                     V_a_sub = V_a_sub)
      }
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1
    if ((j > inits$burnin) & (j %% inits$thin == 0)){
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, diag(sigma_h), h0, as.numeric(h))
    }

    if (j %% 1000 == 0) { cat(" Iteration ", j, " \n")}
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K))
  )
  return(as.mcmc(t(mcmc)))
}

###########################################################################
#' @export
BVAR.Student.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "Student", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b

  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  nu <- inits$nu
  logsigma_nu <- 0
  acount_nu <- 0
  acount_w <- rep(0, t_max)
  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # new precompute
  theta.prior.precmean <- V_b_prior_inv %*% b_prior

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + 1 + K + K + K*t_max + t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i])
    #   b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i]) %*% yt[,i])
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)

    # Sample B
    wt <- as.vector(1/(exp(h/2)*w_sqrt))
    y.tilde <- as.vector( A %*% yt ) * wt
    x.tilde <- kronecker( t(xt), A ) * wt
    theta.prec.chol <- chol( V_b_prior_inv + crossprod(x.tilde) )
    b_sample <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*m) )
    B <- matrix(b_sample,K,m)

    # Sample vol
    ytilde <- A%*% ((yt - B %*% xt)/w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    #   h0[i,] <- as.numeric(latent0(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*%xt)/ w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    if (K > 1) {
      for (i in c(2:K)){
        id_end <- i*(i-1)/2
        id_start <- id_end - i + 2
        a_sub <- a_prior[id_start:id_end]
        V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
        a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                     xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                     a_sub = a_sub,
                                                     V_a_sub = V_a_sub)
      }
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # # Sample w
    # u <- (yt - B %*%xt)
    # for (i in c(1:t_max)){
    #   w_sample[i] <- rinvgamma(1, shape = nu*0.5 + K*0.5, rate = nu*0.5 + 0.5 * t(u[,i]) %*% (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A) %*% u[,i])
    # }
    # w <- reprow(w_sample, K)
    # w_sqrt <- sqrt(w)

    # Sample w
    u <-  A %*% (yt - B %*%xt)
    w_sample <- rinvgamma( n = t_max, shape = (nu+K)*0.5, rate = 0.5*( nu + colSums((u^2)/exp(h)) ) )
    w <- reprow(w_sample, K)
    w_sqrt <- sqrt(w)

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(1)
    if (nu_temp > 4 && nu_temp < 100){
      num_mh = dgamma(nu_temp, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu_temp*0.5, rate = nu_temp*0.5, log = T))
      denum_mh = dgamma(nu, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu*0.5, rate = nu*0.5, log = T))
      alpha = num_mh - denum_mh;
      temp = log(runif(1));
      if (alpha > temp){
        nu = nu_temp
        acount_nu = acount_nu + 1
      }

    }

    if(j %% batchlength == 0 ){
      if (acount_nu > batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu + adaptamount(j %/% batchlength);
      }
      if (acount_nu < batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu - adaptamount(j %/% batchlength);
      }
      acount_nu = 0
    }

    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w_sample))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", round(nu,2)," \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu"),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        paste("w",c(1:t_max), sep = ""))

  return(as.mcmc(t(mcmc)))
}

#############################################################################################
#' @export
BVAR.Skew.Student.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "Skew.Student", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  nu <- inits$nu
  mu.xi <- nu/( nu - 2 )
  logsigma_nu <- 0
  acount_nu <- 0
  acount_w <- rep(0, t_max)
  gamma <- inits$gamma
  D <- diag(gamma, nrow = K)

  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # new precompute
  i <- K*m
  theta.prior.prec = matrix(0,i+K,i+K)
  theta.prior.prec[1:i,1:i] <- V_b_prior_inv
  theta.prior.prec[i+(1:K),i+(1:K)] <- solve( V_gamma_prior )
  theta.prior.precmean <- theta.prior.prec %*% c( b_prior, gamma_prior )

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + 1 + K + K + K*t_max + t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i])
    #   b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i]) %*% (yt[,i] - gamma * w_sample[i] ) )
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)
    #
    # # Sample gamma
    # gamma_post = rep(0, K)
    # V_gamma_post_inv = solve(V_gamma_prior)
    # for (i in c(1:t_max)){
    #   V_gamma_post_inv <- V_gamma_post_inv +  w_sample[i] * t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A # ww' sqrt(w_inv) Signma_inv sqrt(w_inv)
    #   gamma_post <- gamma_post + (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A) %*% (yt[,i] - B %*% xt[,i] )
    # }
    #
    # V_gamma_post <- solve(V_gamma_post_inv)
    # gamma_post <- V_gamma_post %*% ( solve(V_gamma_prior) %*% gamma_prior + gamma_post)
    # gamma <- as.numeric(gamma_post + t(chol(V_gamma_post)) %*% rnorm(K))
    # D <- diag(gamma, nrow = K)

    # Sample B and gamma
    wt <- as.vector(1/(exp(h/2)*w_sqrt))
    y.tilde <- as.vector( A %*% yt ) * wt
    x.tilde <- kronecker( cbind( t(xt), w_sample - mu.xi ), A ) * wt
    theta.prec.chol <- chol( theta.prior.prec + crossprod(x.tilde) )
    theta <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*(m+1)) )
    B <- matrix(theta[1:(m*K)],K,m)
    b_sample <- as.vector(B)
    gamma <- theta[(m*K+1):((m+1)*K)]
    D <- diag(gamma, K)

    # Sample vol
    ytilde <- A%*% ((yt - B %*% xt  - D%*% (w-mu.xi) )/w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # Sample A0
    u_std <- (yt - B %*% xt - D%*% (w-mu.xi))/ w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    if (K > 1) {
      for (i in c(2:K)){
        id_end <- i*(i-1)/2
        id_start <- id_end - i + 2
        a_sub <- a_prior[id_start:id_end]
        V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
        a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                     xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                     a_sub = a_sub,
                                                     V_a_sub = V_a_sub)
      }
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1


    # Sample w
    q2 <- colSums( ( ( A %*% ( yt - B%*%xt + mu.xi * gamma ) ) / sqrtvol )^2 )
    p2 <- colSums( ( c( A %*% gamma ) / sqrtvol )^2 )
    w_sample <- mapply( GIGrvg::rgig, n = 1, lambda = -(nu+K)*0.5, chi = nu + q2,
                         psi = p2 )
    w <- reprow(w_sample, K)
    w_sqrt <- sqrt(w)

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(1)
    if (nu_temp > 4 && nu_temp < 100){
      num_mh = dgamma(nu_temp, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu_temp*0.5, rate = nu_temp*0.5, log = T))
      denum_mh = dgamma(nu, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu*0.5, rate = nu*0.5, log = T))
      alpha = num_mh - denum_mh;
      temp = log(runif(1));
      if (alpha > temp){
        nu = nu_temp
        acount_nu = acount_nu + 1
      }

    }
    mu.xi <- nu/( nu - 2 )

    if(j %% batchlength == 0 ){
      if (acount_nu > batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu + adaptamount(j %/% batchlength);
      }
      if (acount_nu < batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu - adaptamount(j %/% batchlength);
      }
      acount_nu = 0
    }
    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, gamma, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w_sample))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2) , " ", gamma ,  " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("gamma",c(1:K), sep = ""),
                        paste("nu"),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        paste("w",c(1:t_max), sep = ""))

  return(as.mcmc(t(mcmc)))
}

###########################################################################
#' @export
BVAR.MT.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "MT", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b

  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)
#  w_sqrt_inv <- 1/w_sqrt

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # new precompute
  theta.prior.precmean <- V_b_prior_inv %*% b_prior

  # Output
  mcmc <- matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K*t_max + K*t_max,
                ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = solve(V_b_prior)
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), diag(w_sqrt_inv[,i]) %*% t(A)%*% diag(1/exp(h[,i])) %*% A %*% diag(w_sqrt_inv[,i]))
    #   b_post <- b_post + kronecker(xt[,i], (diag(w_sqrt_inv[,i]) %*% t(A)%*% diag(1/exp(h[,i])) %*% A %*% diag(w_sqrt_inv[,i])) %*% yt[,i])
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)

    # Sample B
    wt <- as.vector(exp(-h/2))
    y.tilde <- as.vector( A %*% ( w^(-0.5) * yt ) ) * wt
    # make and stack X matrices
    x.tilde <- kronecker( t(xt), A )
    tmp <- kronecker( t(w^(-1/2)), matrix(1,K,1) ) * wt # W^(-1/2)
    ii <- 0
    for (i in 1:m) {
      x.tilde[,(ii+1):(ii+K)] <- x.tilde[,(ii+1):(ii+K)] * tmp
      ii <- ii+K
    }
    theta.prec.chol <- chol( V_b_prior_inv + crossprod(x.tilde) )
    b_sample <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*m) )
    B <- matrix(b_sample,K,m)

    # Sample vol
    ytilde <- A%*% ((yt - B %*%xt)/ w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*%xt) / w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # Sample w
    u <- (yt - B %*%xt)
    u_proposal <- A %*% (yt - B %*%xt)
    a_target <- (nu*0.5 + 1*0.5) * 0.75
    b_target <- (nu*0.5 + 0.5 * u_proposal^2 / sqrtvol^2)*0.75  # adjust by 0.75
    w_temp <- matrix( rinvgamma( n = K*t_max, shape = a_target, rate = b_target ), K, t_max ) # rinvgamma recycles arguments
    w_temp_sqrt <- sqrt(w_temp)
    # num_mh <- sapply(1:t_max, function (i) { Sigma2_inv_t = t(A)%*% diag(1/exp(h[,i])) %*% A; u.w <- w_temp_sqrt_inv[,i] * u[,i]
    #                                          sum(dinvgamma(w_temp[,i], shape = nu*0.5, rate = nu*0.5, log = T)) + #prior
    #                                           ( - 0.5 * sum(log_w_temp[,i]) - 0.5 * t(u.w) %*% Sigma2_inv_t %*% u.w ) - # posterior
    #                                           sum( dinvgamma(w_temp[,i], shape = a_target, rate = b_target[,i], log = T)) # proposal
    #                                        }
    #                   )
    # denum_mh <- sapply(1:t_max, function (i) { Sigma2_inv_t = t(A)%*% diag(1/exp(h[,i])) %*% A; u.w <- u[,i]/sqrt(w[,i])
    #                                            sum(dinvgamma(w[,i], shape = nu*0.5, rate = nu*0.5, log = T)) + #prior
    #                                            ( - 0.5 * sum(log(w[,i])) - 0.5 * t(u.w) %*% Sigma2_inv_t %*% u.w ) - # posterior
    #                                            sum( dinvgamma(w[,i], shape = a_target, rate = b_target[,i], log = T)) # proposal
    #                                          }
    #                   )
    A.w.u.s.prop <- ( A %*% ( u/w_temp_sqrt ) )/sqrtvol
    num_mh <- colSums(dinvgamma(w_temp, shape = nu*0.5, rate = nu*0.5, log = T)) - #prior
              0.5*( colSums(log(w_temp)) + colSums(A.w.u.s.prop^2) ) - # posterior
              colSums( dinvgamma(w_temp, shape = a_target, rate = b_target, log = T)) # proposal
    A.w.u.s.curr <- ( A %*% ( u/w_sqrt ) )/sqrtvol
    denum_mh <- colSums(dinvgamma(w, shape = nu*0.5, rate = nu*0.5, log = T)) - #prior
                0.5*( colSums(log(w)) + colSums(A.w.u.s.curr^2) ) - # posterior
                colSums( dinvgamma(w, shape = a_target, rate = b_target, log = T)) # proposal
    acc <- (num_mh-denum_mh) > log(runif(t_max))
    w[,acc] <- w_temp[,acc]
    w_sqrt[,acc] <- w_temp_sqrt[,acc]
    acount_w[acc] <- acount_w[acc] + 1

    # # Sample w
    # u <- (yt - B %*%xt)
    # u_proposal <- A %*% (yt - B %*%xt)
    # for (i in c(1:t_max)){
    #   a_target <- (nu*0.5 + 1*0.5) * 0.75 # adjust by 0.75
    #   b_target <- (nu*0.5 + 0.5 * u_proposal[,i]^2 / sqrtvol[,i]^2) * 0.75 # adjust by 0.75
    #   w_temp = w[,i]
    #   for (k in c(1:K)){
    #     w_temp[k] <- rinvgamma(1, shape = a_target[k], rate = b_target[k])
    #   }
    #   log_w_temp <- log(w_temp)
    #   w_temp_sqrt_inv = 1 / sqrt(w_temp)
    #
    #   #total_log_dens <- (-0.5) * sum(log(w[,i])) - 0.5 * t(u[,i]) %*% diag(w_sqrt_inv[,i]) %*% Sigma2_inv %*% diag(w_sqrt_inv[,i]) %*% u[,i]  + sum(dinvgamma(w[,i], shape = nu*0.5, rate = nu*0.5, log = T))
    #   Sigma2_inv_t = t(A)%*% diag(1/exp(h[,i])) %*% A
    #   num_mh =  sum(dinvgamma(w_temp, shape = nu*0.5, rate = nu*0.5, log = T)) +  # prior
    #     ( - 0.5 * sum(log_w_temp) - 0.5 * t(u[,i]) %*% diag(w_temp_sqrt_inv) %*% Sigma2_inv_t %*% diag(w_temp_sqrt_inv) %*% u[,i]) - # posterior
    #     sum( dinvgamma(w_temp, shape = a_target, rate = b_target, log = T)) # proposal
    #   denum_mh = sum(dinvgamma(w[,i], shape = nu*0.5, rate = nu*0.5, log = T)) + # prior
    #     (-0.5) * sum(log(w[,i])) - 0.5 * t(u[,i]) %*% diag(w_sqrt_inv[,i]) %*% Sigma2_inv_t %*% diag(w_sqrt_inv[,i]) %*% u[,i]  - # posterior
    #     sum( dinvgamma(w[,i], shape = a_target, rate = b_target, log = T))
    #   alpha = num_mh - denum_mh;
    #   temp = log(runif(1))
    #   if (alpha > temp){
    #     w[,i] <- w_temp
    #     w_sqrt[,i] <- sqrt(w_temp)
    #     w_sqrt_inv[,i] <- 1/sqrt(w_temp)
    #     acount_w[i] <- acount_w[i] + 1
    #   }
    # }



    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }
    }

    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }

    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w))

    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2), " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)))

  return(as.mcmc(t(mcmc)))
}

#############################################################################################
#' @export
BVAR.MST.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "MST", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  mu.xi <- nu/( nu - 2 )
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  gamma <- inits$gamma
  D <- diag(gamma)

  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # new precompute
  i <- K*m
  theta.prior.prec = matrix(0,i+K,i+K)
  theta.prior.prec[1:i,1:i] <- V_b_prior_inv
  theta.prior.prec[i+(1:K),i+(1:K)] <- solve( V_gamma_prior )
  theta.prior.precmean <- theta.prior.prec %*% c( b_prior, gamma_prior )

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K + K*t_max + K*t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), diag(w_sqrt_inv[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt_inv[,i]) )
    #   b_post <- b_post + kronecker(xt[,i], (diag(w_sqrt_inv[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt_inv[,i])) %*% (yt[,i] - gamma * w[,i]))
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)
    #
    # gamma_post = rep(0, K)
    # V_gamma_post_inv = solve(V_gamma_prior)
    # for (i in c(1:t_max)){
    #   V_gamma_post_inv <- V_gamma_post_inv +  diag(w_sqrt[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt[,i]) # ww' sqrt(w_inv) Signma_inv sqrt(w_inv)
    #   gamma_post <- gamma_post + (diag(w_sqrt[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt_inv[,i])) %*% (yt[,i] - B %*% xt[,i] )
    # }
    #
    # V_gamma_post <- solve(V_gamma_post_inv)
    # gamma_post <- V_gamma_post %*% ( solve(V_gamma_prior) %*% gamma_prior + gamma_post)
    # gamma <- as.numeric(gamma_post + t(chol(V_gamma_post)) %*% rnorm(K))
    # D <- diag(gamma)


    # Sample B and gamma
    wt <- as.vector(exp(-h/2))
    y.tilde <- as.vector( A %*% ( w^(-0.5) * yt ) ) * wt
    # make and stack X matrices
    x1 <- kronecker( t(xt), A )
    tmp <- kronecker( t(w^(-1/2)), matrix(1,K,1) ) * wt # W^(-1/2)
    ii <- 0
    for (i in 1:m) {
      x1[,(ii+1):(ii+K)] <- x1[,(ii+1):(ii+K)] * tmp
      ii <- ii+K
    }
    W.mat <- matrix(0,K*t_max,K)
    idx <- (0:(t_max-1))*K
    for ( i in 1:K ) {
      W.mat[idx + (i-1)*t_max*K + i] <- (w[i,] - mu.xi[i])/sqrt(w[i,]) # W^(-1/2) * (W - bar(W))
    }
    x2 <- matrix(0,K*t_max,K)
    for ( i in 1:K ) {
      x2[,i] <- as.vector( A%*%matrix(W.mat[,i],K,t_max) ) * wt
    }
    x.tilde <- cbind( x1, x2 )
    theta.prec.chol <- chol( theta.prior.prec + crossprod(x.tilde) )
    theta <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*(m+1)) )
    B <- matrix(theta[1:(m*K)],K,m)
    b_sample <- as.vector(B)
    gamma <- theta[(m*K+1):((m+1)*K)]
    D <- diag(gamma)

    # Sample vol
    ytilde <- A%*% ((yt - B %*%xt - D %*% ( w - mu.xi ))/ w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # Sample A0
    u_std <- (yt - B %*% xt - D%*% ( w - mu.xi )) / w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # Sample w
    u <- (yt - B %*%xt + mu.xi * gamma  )
    u_proposal <- A %*% (yt - B %*%xt) + mu.xi * gamma
    a_target <- (nu*0.5 + 1*0.5) * 0.75
    b_target <- (nu*0.5 + 0.5 * u_proposal^2 / sqrtvol^2)*0.75  # adjust by 0.75
    w_temp <- matrix( rinvgamma( n = K*t_max, shape = a_target, rate = b_target ), K, t_max ) # rinvgamma recycles arguments
    log_w_temp <- log(w_temp)
    w_temp_sqrt <- sqrt(w_temp)

    A.w.u.s.prop <- ( A %*% ( u/w_temp_sqrt - w_temp_sqrt*gamma ) )/sqrtvol
    num_mh <- colSums(dinvgamma(w_temp, shape = nu*0.5, rate = nu*0.5, log = T)) - #prior
      0.5*( colSums(log(w_temp)) + colSums(A.w.u.s.prop^2) ) - # posterior
      colSums( dinvgamma(w_temp, shape = a_target, rate = b_target, log = T)) # proposal
    A.w.u.s.curr <- ( A %*% ( u/w_sqrt - w_sqrt*gamma ) )/sqrtvol
    denum_mh <- colSums(dinvgamma(w, shape = nu*0.5, rate = nu*0.5, log = T)) - #prior
      0.5*( colSums(log(w)) + colSums(A.w.u.s.curr^2) ) - # posterior
      colSums( dinvgamma(w, shape = a_target, rate = b_target, log = T)) # proposal
    acc <- (num_mh-denum_mh) > log(runif(t_max))
    w[,acc] <- w_temp[,acc]
    w_sqrt[,acc] <- w_temp_sqrt[,acc]
    acount_w[acc] <- acount_w[acc] + 1

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }

    }
    mu.xi <- nu / ( nu - 2 )

    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }
    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, gamma, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2), " gamma ", round(gamma,2) ," \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("gamma",c(1:K), sep = ""),
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K))
  )

  return(as.mcmc(t(mcmc)))
}

###########################################################################
#' @export
BVAR.OT.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "OT", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b

  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # new precompute
  theta.prior.precmean <- V_b_prior_inv %*% b_prior

  # Output
  mcmc <- matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K*t_max + K*t_max,
                 ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]) / w[,i] ) %*% A )
    #   b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i])/ w[,i] ) %*% A) %*% yt[,i])
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)

    # Sample B
    wt <- as.vector(1/(exp(h/2)*w_sqrt))
    y.tilde <- as.vector( A %*% yt ) * wt
    x.tilde <- kronecker( t(xt), A ) * wt
    theta.prec.chol <- chol( V_b_prior_inv + crossprod(x.tilde) )
    b_sample <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                            + rnorm(K*m) )
    B <- matrix(b_sample,K,m)

    # Sample vol
    ytilde <- (A %*% (yt - B %*% xt)) / w_sqrt
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*%xt) # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,] / w_sqrt[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / reprow(sqrtvol[i,],i-1) / reprow(w_sqrt[i,], i-1), nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # # Sample w
    # u <-  A %*% (yt - B %*%xt)
    # for (i in c(1:t_max)){
    #   w[,i] <- mapply(rinvgamma, n = 1, shape = nu*0.5 + 1*0.5, rate = nu*0.5 + 0.5 * (u[,i]^2) / exp(h[,i]))
    # }
    # w_sqrt <- sqrt(w)

    # Sample w
    u <-  A %*% (yt - B %*%xt)
    w <- matrix( rinvgamma( n = K*t_max, shape = nu*0.5 + 1*0.5, rate = nu*0.5 + 0.5 * (u^2) / exp(h) ),
                 K, t_max )
    w_sqrt <- sqrt(w)


    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }
    }


    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }

    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w))

    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", round(nu,2)," \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)))

  return(as.mcmc(t(mcmc)))
}

#############################################################################################
#' @export
BVAR.OST.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "OST", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  mu.xi <- nu/( nu - 2 )
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  gamma <- inits$gamma
  D <- diag(gamma)

  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # new precompute
  i <- K*m
  theta.prior.prec = matrix(0,i+K,i+K)
  theta.prior.prec[1:i,1:i] <- V_b_prior_inv
  theta.prior.prec[i+(1:K),i+(1:K)] <- solve( V_gamma_prior )
  theta.prior.precmean <- theta.prior.prec %*% c( b_prior, gamma_prior )

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K + K*t_max + K*t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # # Sample B
    # b_post = rep(0, m*K)
    # V_b_post_inv = V_b_prior_inv
    # inv_A <- solve(A)
    # for (i in c(1:t_max)){
    #   V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]) / w[,i] ) %*% A )
    #   b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i])/ w[,i] ) %*% A) %*% (yt[,i] - inv_A %*% (gamma * w[,i])) )
    # }
    #
    # V_b_post <- solve(V_b_post_inv)
    # b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    # b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    # B <- Vec_to_Mat(b_sample, K,p)
    #
    # # Sample gamma
    # gamma_post = rep(0, K)
    # V_gamma_post_inv = solve(V_gamma_prior)
    # u <- A %*% (yt - B %*%xt)
    #
    # V_gamma_post_inv <- solve(V_gamma_prior) + diag(apply(w / exp(h), MARGIN = 1, FUN = sum))
    # gamma_post <- apply(u / exp(h), MARGIN = 1, FUN = sum)
    #
    # V_gamma_post <- solve(V_gamma_post_inv)
    # gamma_post <- V_gamma_post %*% ( solve(V_gamma_prior) %*% gamma_prior + gamma_post)
    # gamma <- as.numeric(gamma_post + t(chol(V_gamma_post)) %*% rnorm(K))
    # D <- diag(gamma)

    # Sample B and gamma
    wt <- as.vector( 1/sqrt(w*exp(h)) )
    y.tilde <- as.vector( A %*% yt ) * wt
    # make and stack diagonal W matrices
    W.mat <- matrix(0,K*t_max,K)
    idx <- (0:(t_max-1))*K
    for ( i in 1:K ) {
      W.mat[idx + (i-1)*t_max*K + i] <- w[i,] -  mu.xi[i]
    }
    x.tilde <- cbind( kronecker( t(xt), A ), W.mat ) * wt
    theta.prec.chol <- chol( theta.prior.prec + crossprod(x.tilde) )
    theta <- backsolve( theta.prec.chol,
                        backsolve( theta.prec.chol, theta.prior.precmean + crossprod( x.tilde, y.tilde ),
                                   upper.tri = T, transpose = T )
                        + rnorm(K*(m+1)) )
    B <- matrix(theta[1:(m*K)],K,m)
    b_sample <- as.vector(B)
    gamma <- theta[(m*K+1):((m+1)*K)]
    D <- diag(gamma)

    # Sample vol
    ytilde <- (A %*% (yt - B %*% xt) - D %*% ( w - mu.xi ) ) / w_sqrt
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # Sample A0
    u_std <- (yt - B %*% xt) # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = (u_std[i,] - (w[i,]  - mu.xi[i] ) * gamma[i]) / sqrtvol[i,] / w_sqrt[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / reprow(sqrtvol[i,],i-1) / reprow(w_sqrt[i,], i-1), nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # Sample w
    u <-  A %*% (yt - B %*%xt) + gamma * mu.xi
    w <- matrix( mapply( GIGrvg::rgig, n = 1, lambda = -(nu+1)*0.5, chi = nu + (u^2)/exp(h),
                         psi = gamma^2/exp(h) ), K, t_max )
    w_sqrt <- sqrt(w)

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }
    }
    mu.xi <- nu / ( nu - 2 )

    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }

    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, gamma, nu, diag(sigma_h), h0, as.numeric(h), as.numeric(w))
    if (j %% 1000 == 0) {
      cat(" Iteration- ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2) , " ", round(gamma,2) ,  " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("gamma",c(1:K), sep = ""),
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)))

  return(as.mcmc(t(mcmc)))
}

#############################################################################################
#' @export
BVAR.dynSkew.Student.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "dynSkew.Student", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  nu <- inits$nu
  logsigma_nu <- 0
  acount_nu <- 0
  acount_w <- rep(0, t_max)
  Gamma_T <- inits$gamma
  if (length(Gamma_T) <= K){
    Gamma_T = repcol(Gamma_T, t_max)
  }
  gamma0 = Gamma_T[,1]
  sigma_G = diag(1e-4, nrow = K)
  #D <- diag(gamma, nrow = K)

  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + 1 + K + K + K + K + K * t_max + K*t_max + t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # Sample B
    b_post = rep(0, m*K)
    V_b_post_inv = V_b_prior_inv
    for (i in c(1:t_max)){
      V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i])
      b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A /w_sample[i]) %*% (yt[,i] - Gamma_T[,i] * w_sample[i] ) )
    }

    V_b_post <- solve(V_b_post_inv)
    b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    B <- Vec_to_Mat(b_sample, K,p)

    # Sample gamma : carterkohn(y,Z,Ht)
    # y: n*T
    # Z: (nT)*(n_max*K)
    # Ht: (nT)*n
    # ysim <-  Z %*% B_t + t(chol(Hsim)) %*% rnorm(M)
    Sigma_T <- matrix(0, nrow = K * t_max, ncol = K)
    A_inv <- solve(A)
    for (i in c(1:t_max)){
      Sigma_T[ ((i-1)*K + 1):(i*K) , ] <- A_inv %*% diag(exp(h[,i]), nrow = K) %*% t(A_inv) * w_sample[i]
    }
    Gamma_filter <- fatBVARS::carterkohn(yt - B %*% xt, w_sample %x% diag(K),
                         Sigma_T, sigma_G, K, K, t_max,
                         gamma_prior,V_gamma_prior)

    Gamma_T <- Gamma_filter$bdraws
    gamma0 <- Gamma_filter$b0draws

    # Draw sigma_G, the covariance of B(t)
    if (K>1) {
      sse_2 <- apply( (Gamma_T[,1:t_max] - cbind(gamma0,Gamma_T[,1:(t_max-1)]) )^2, MARGIN = 1, FUN = sum)
    } else {
      sse_2 <- sum( (Gamma_T[,1:t_max] - c(gamma0,Gamma_T[,1:(t_max-1)]) )^2)
    }
    sigma_G_a <- rep(t_max,K) # prior of sigma_G
    sigma_G_b <- sse_2 # prior of sigma_G

    for (i in c(1:K)){
      sigma_new <- rinvgamma(1, shape = sigma_G_a[i] * 0.5, rate = sigma_G_b[i] * 0.5)
      alpha = (sigma_G[i,i] - sigma_new) / 2 + 0.5 * (log(sigma_new) - log(sigma_G[i,i])) # B_sigma = 1
      temp = log(runif(1))
      if (alpha > temp){
        sigma_G[i,i] <- sigma_new
      }
    }


    # Sample vol
    ytilde <- A%*% ((yt - B %*% xt  - Gamma_T * w)/w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*% xt - Gamma_T * w)/ w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    if (K > 1) {
      for (i in c(2:K)){
        id_end <- i*(i-1)/2
        id_start <- id_end - i + 2
        a_sub <- a_prior[id_start:id_end]
        V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
        a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                     xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                     a_sub = a_sub,
                                                     V_a_sub = V_a_sub)
      }
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1


    # Sample w
    u <- (yt - B %*%xt)
    for (i in c(1:t_max)){
      ran_stu = rt(1, df = 4)
      Sigma2_inv_t = (t(A)%*% diag(1/exp(h[,i]), nrow = K) %*% A)
      a_eq <- -0.5* (t(Gamma_T[,i]) %*% Sigma2_inv_t %*% Gamma_T[,i])
      b_eq <- -((0.5*K + nu/2+1))
      c_eq =  0.5 * t(u[,i]) %*% Sigma2_inv_t %*% u[,i] + nu/2
      w_mode <- (- b_eq - sqrt(b_eq^2-4*a_eq*c_eq))/(2*a_eq)
      w_mode <- ifelse(abs(w_mode)<1e-5, - c_eq/b_eq, w_mode)

      mode_target <- log( w_mode )
      cov_target <- 1/( c_eq / w_mode - a_eq * w_mode) * 1.2 # scale by 1.2

      log_w_temp = mode_target + sqrt(cov_target)*ran_stu
      w_temp = exp(log_w_temp)

      num_mh =  dinvgamma(w_temp, shape = nu*0.5, rate = nu*0.5, log = T) +  # prior
        ( - 0.5 * K * log_w_temp - 0.5 * t((u[,i] - Gamma_T[,i]*w_temp)) %*% Sigma2_inv_t %*% (u[,i] - Gamma_T[,i]*w_temp) / w_temp) - # posterior
        dt( (log_w_temp - mode_target)/sqrt(cov_target), df = 4, log = T) + log_w_temp # proposal
      denum_mh = dinvgamma(w_sample[i], shape = nu*0.5, rate = nu*0.5, log = T) +  # prior
        ( - 0.5 * K * log(w_sample[i]) - 0.5 * t((u[,i] - Gamma_T[,i]*w_sample[i])) %*% Sigma2_inv_t %*% (u[,i] - Gamma_T[,i]*w_sample[i]) / w_sample[i]) - # posterior
        dt( (log(w_sample[i]) - mode_target)/sqrt(cov_target), df = 4, log = T) + log(w_sample[i]) # proposal
      alpha = num_mh - denum_mh;
      temp = log(runif(1))
      if (alpha > temp){
        w_sample[i] <- w_temp
        w[,i] <- w_temp
        w_sqrt[,i] <- sqrt(w_temp)
        acount_w[i] <- acount_w[i] + 1
      }
    }

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(1)
    if (nu_temp > 4 && nu_temp < 100){
      num_mh = dgamma(nu_temp, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu_temp*0.5, rate = nu_temp*0.5, log = T))
      denum_mh = dgamma(nu, shape = nu_gam_a, rate = nu_gam_b, log = T) +
        sum(dinvgamma(w_sample, shape = nu*0.5, rate = nu*0.5, log = T))
      alpha = num_mh - denum_mh;
      temp = log(runif(1));
      if (alpha > temp){
        nu = nu_temp
        acount_nu = acount_nu + 1
      }

    }

    if(j %% batchlength == 0 ){
      if (acount_nu > batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu + adaptamount(j %/% batchlength);
      }
      if (acount_nu < batchlength * TARGACCEPT){
        logsigma_nu = logsigma_nu - adaptamount(j %/% batchlength);
      }
      acount_nu = 0
    }
    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, diag(sigma_G), gamma0, as.numeric(Gamma_T), as.numeric(h), as.numeric(w_sample))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2) , " ",
          round(apply(Gamma_T,MARGIN = 1,FUN = mean),2) ,  " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu"),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        paste("sigma_G",c(1:K), sep = ""),
                        paste("G0",c(1:K), sep = ""),
                        sprintf("gamma_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        paste("w",c(1:t_max), sep = ""))

  return(as.mcmc(t(mcmc)))
}


#############################################################################################
#' @export
BVAR.dynMST.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "dynMST", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  Gamma_T <- inits$gamma
  if (length(Gamma_T) <= K){
    Gamma_T = repcol(Gamma_T, t_max)
  }
  gamma0 = Gamma_T[,1]
  sigma_G = diag(1e-4, nrow = K)
  #D <- diag(gamma, nrow = K)


  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)
  w_sqrt_inv <- 1/w_sqrt

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K + K + K * t_max + K*t_max + K*t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # Sample B
    b_post = rep(0, m*K)
    V_b_post_inv = V_b_prior_inv
    for (i in c(1:t_max)){
      V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), diag(w_sqrt_inv[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt_inv[,i]) )
      b_post <- b_post + kronecker(xt[,i], (diag(w_sqrt_inv[,i]) %*% (t(A)%*% diag(1/exp(h[,i])) %*% A) %*% diag(w_sqrt_inv[,i])) %*% (yt[,i] - Gamma_T[,i] * w[,i]))
    }

    V_b_post <- solve(V_b_post_inv)
    b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    B <- Vec_to_Mat(b_sample, K,p)

    # Sample gamma : carterkohn(y,Z,Ht)
    # y: n*T
    # Z: (nT)*(n_max*K)
    # Ht: (nT)*n
    # ysim <-  Z %*% B_t + t(chol(Hsim)) %*% rnorm(M)
    Sigma_T <- matrix(0, nrow = K * t_max, ncol = K)
    Z_T <- matrix(0, nrow = K * t_max, ncol = K)
    A_inv <- solve(A)
    for (i in c(1:t_max)){
      Sigma_T[ ((i-1)*K + 1):(i*K) , ] <- diag(w_sqrt[,i], nrow = K) %*% (A_inv %*% diag(exp(h[,i]), nrow = K) %*% t(A_inv)) %*% diag(w_sqrt[,i], nrow = K)
      Z_T[ ((i-1)*K + 1):(i*K) , ] <- diag(w[,i], nrow = K)
    }
    Gamma_filter <- fatBVARS::carterkohn(yt - B %*% xt, Z_T,
                                         Sigma_T, sigma_G, K, K, t_max,
                                         gamma_prior,V_gamma_prior)

    Gamma_T <- Gamma_filter$bdraws
    gamma0 <- Gamma_filter$b0draws

    # Draw sigma_G, the covariance of B(t)
    if (K>1) {
      sse_2 <- apply( (Gamma_T[,1:t_max] - cbind(gamma0,Gamma_T[,1:(t_max-1)]) )^2, MARGIN = 1, FUN = sum)
    } else {
      sse_2 <- sum( (Gamma_T[,1:t_max] - c(gamma0,Gamma_T[,1:(t_max-1)]) )^2)
    }
    sigma_G_a <- rep(t_max,K) # prior of sigma_G
    sigma_G_b <- sse_2 # prior of sigma_G

    for (i in c(1:K)){
      sigma_new <- rinvgamma(1, shape = sigma_G_a[i] * 0.5, rate = sigma_G_b[i] * 0.5)
      alpha = (sigma_G[i,i] - sigma_new) / 2 + 0.5 * (log(sigma_new) - log(sigma_G[i,i])) # B_sigma = 1
      temp = log(runif(1))
      if (alpha > temp){
        sigma_G[i,i] <- sigma_new
      }
    }


    # Sample vol
    ytilde <- A%*% ((yt - B %*%xt - Gamma_T * w)/ w_sqrt)
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # for (i in c(1:K)){
    #   svdraw[[i]] <- svsample2(ytilde[i,], startpara = para(svdraw[[i]]),
    #                            startlatent = latent(svdraw[[i]]), priormu = priormu,
    #                            priorphi = priorphi, priorsigma = priorsigma)
    #   paravol[i,] <- para(svdraw[[i]])
    #   h[i,] <- as.numeric(latent(svdraw[[i]]))
    # }
    # sqrtvol <- exp(h/2)

    # Sample A0
    u_std <- (yt - B %*% xt - Gamma_T * w) / w_sqrt # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = u_std[i,] / sqrtvol[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / sqrtvol[i,], nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # Sample w
    u <- (yt - B %*%xt)
    u_proposal <- A %*% (yt - B %*%xt)
    for (i in c(1:t_max)){
      a_target <- (nu*0.5 + 1*0.5) * 0.75
      b_target <- (nu*0.5 + 0.5 * u_proposal[,i]^2 / sqrtvol[,i]^2)*0.75  # adjust by 0.75
      w_temp = w[,i]
      for (k in c(1:K)){
        w_temp[k] <- rinvgamma(1, shape = a_target[k], rate = b_target[k])
      }
      log_w_temp <- log(w_temp)
      w_temp_sqrt_inv = 1 / sqrt(w_temp)

      #total_log_dens <- (-0.5) * sum(log(w[,i])) - 0.5 * t(u[,i]) %*% diag(w_sqrt_inv[,i]) %*% Sigma2_inv %*% diag(w_sqrt_inv[,i]) %*% u[,i]  + sum(dinvgamma(w[,i], shape = nu*0.5, rate = nu*0.5, log = T))
      Sigma2_inv_t = t(A)%*% diag(1/exp(h[,i])) %*% A
      num_mh =  sum(dinvgamma(w_temp, shape = nu*0.5, rate = nu*0.5, log = T)) +  # prior
        ( - 0.5 * sum(log_w_temp) - 0.5 * t(u[,i] - Gamma_T[,i]*w_temp) %*% diag(w_temp_sqrt_inv) %*% Sigma2_inv_t %*% diag(w_temp_sqrt_inv) %*% (u[,i] - Gamma_T[,i]*w_temp)) - # posterior
        sum( dinvgamma(w_temp, shape = a_target, rate = b_target, log = T)) # proposal
      denum_mh = sum(dinvgamma(w[,i], shape = nu*0.5, rate = nu*0.5, log = T)) + # prior
        (-0.5) * sum(log(w[,i])) - 0.5 * t(u[,i] - Gamma_T[,i]*w[,i]) %*% diag(w_sqrt_inv[,i]) %*% Sigma2_inv_t %*% diag(w_sqrt_inv[,i]) %*% (u[,i] - Gamma_T[,i]*w[,i])  - # posterior
        sum( dinvgamma(w[,i], shape = a_target, rate = b_target, log = T))
      alpha = num_mh - denum_mh;
      temp = log(runif(1))
      if (alpha > temp){
        w[,i] <- w_temp
        w_sqrt[,i] <- sqrt(w_temp)
        w_sqrt_inv[,i] <- 1/sqrt(w_temp)
        acount_w[i] <- acount_w[i] + 1
      }
    }

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }

    }

    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }
    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, diag(sigma_G), gamma0, as.numeric(Gamma_T), as.numeric(h), as.numeric(w))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2),
          round(apply(Gamma_T,MARGIN = 1,FUN = mean),2) ,  " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        paste("sigma_G",c(1:K), sep = ""),
                        paste("G0",c(1:K), sep = ""),
                        sprintf("gamma_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K))
  )

  return(as.mcmc(t(mcmc)))
}

#############################################################################################
#' @export
BVAR.dynOST.SV <- function(y, K, p, y0 = NULL, prior = NULL, inits = NULL){
  # Init regressors in the right hand side
  t_max <- nrow(y)
  yt = t(y)
  xt <- makeRegressor(y, y0, t_max, K, p)

  # Init prior and initial values
  m = K * p + 1
  if (is.null(prior)){
    prior <- get_prior(y, p, priorStyle = "Minnesota", dist = "MST", SV = TRUE)
  }
  # prior B
  b_prior = prior$b_prior
  V_b_prior = prior$V_b_prior
  # prior sigma
  sigma0_T0 <- prior$sigma_T0
  sigma0_S0 <- prior$sigma_S0
  # prior h0_mean
  h0_mean <- 2 * log(prior$sigma)
  # prior A
  a_prior = prior$a_prior
  V_a_prior = prior$V_a_prior
  # prior nu
  nu_gam_a = prior$nu_gam_a
  nu_gam_b = prior$nu_gam_b
  # prior gamma
  gamma_prior = prior$gamma_prior
  V_gamma_prior = prior$V_gamma_prior
  # Initial values
  if (is.null(inits)){
    inits <- get_init(prior)
  }
  samples <- inits$samples
  A <- inits$A0
  B <- inits$B0
  h <- inits$h
  sigma_h <- inits$sigma_h

  V_b_prior_inv <- solve(V_b_prior)

  # Multi degrees of freedom
  nu <- inits$nu
  logsigma_nu <- rep(0,K)
  acount_nu <- rep(0,K)
  acount_w <- rep(0, t_max)
  Gamma_T <- inits$gamma
  if (length(Gamma_T) <= K){
    Gamma_T = repcol(Gamma_T, t_max)
  }
  gamma0 = Gamma_T[,1]
  sigma_G = diag(1e-4, nrow = K)
  #D <- diag(gamma, nrow = K)

  # Init w as Gaussian
  w_sample <- rep(1, t_max)
  w <- reprow(w_sample, K)
  w_sqrt <- sqrt(w)

  # svdraw <- list()
  # paravol <- matrix(0, ncol = 3, nrow = K)
  # for (i in c(1:K)){
  #   svdraw[[i]] <- list(para = c(mu = 0, phi = 0.95, sigma = 0.2),
  #                       latent = h[i,])
  # }

  # Output
  mcmc <-  matrix(NA, nrow = m*K + 0.5*K*(K-1) + K + K + K + K + K + K * t_max + K*t_max + K*t_max,
                  ncol = (samples - inits$burnin)%/% inits$thin)
  for (j in c(1:samples)){
    # Sample B
    b_post = rep(0, m*K)
    V_b_post_inv = V_b_prior_inv
    inv_A <- solve(A)
    for (i in c(1:t_max)){
      V_b_post_inv <- V_b_post_inv + kronecker(xt[,i] %*% t(xt[,i]), t(A)%*% diag(1/exp(h[,i]) / w[,i] ) %*% A )
      b_post <- b_post + kronecker(xt[,i], (t(A)%*% diag(1/exp(h[,i])/ w[,i] ) %*% A)  ) %*% (yt[,i] - inv_A %*% (Gamma_T[,i] * w[,i]))
    }

    V_b_post <- solve(V_b_post_inv)
    b_post <- V_b_post %*% ( solve(V_b_prior) %*% b_prior + b_post)
    b_sample <- b_post + t(chol(V_b_post)) %*% rnorm(m*K)
    B <- Vec_to_Mat(b_sample, K,p)

    # Sample gamma : carterkohn(y,Z,Ht)
    # y: n*T
    # Z: (nT)*(n_max*K)
    # Ht: (nT)*n
    # ysim <-  Z %*% B_t + t(chol(Hsim)) %*% rnorm(M)
    u <- A %*% (yt - B %*%xt)
    Sigma_T <- matrix(0, nrow = K * t_max, ncol = K)
    Z_T <- matrix(0, nrow = K * t_max, ncol = K)
    A_inv <- solve(A)
    for (i in c(1:t_max)){
      Sigma_T[ ((i-1)*K + 1):(i*K) , ] <- diag(w[,i] * exp(h[,i]), nrow = K)
      Z_T[ ((i-1)*K + 1):(i*K) , ] <- diag(w[,i], nrow = K)
    }
    Gamma_filter <- fatBVARS::carterkohn(u, Z_T,
                                         Sigma_T, sigma_G, K, K, t_max,
                                         gamma_prior,V_gamma_prior)

    Gamma_T <- Gamma_filter$bdraws
    gamma0 <- Gamma_filter$b0draws

    # Draw sigma_G, the covariance of B(t)
    if (K>1) {
      sse_2 <- apply( (Gamma_T[,1:t_max] - cbind(gamma0,Gamma_T[,1:(t_max-1)]) )^2, MARGIN = 1, FUN = sum)
    } else {
      sse_2 <- sum( (Gamma_T[,1:t_max] - c(gamma0,Gamma_T[,1:(t_max-1)]) )^2)
    }
    sigma_G_a <- rep(t_max,K) # prior of sigma_G
    sigma_G_b <- sse_2 # prior of sigma_G

    for (i in c(1:K)){
      sigma_new <- rinvgamma(1, shape = sigma_G_a[i] * 0.5, rate = sigma_G_b[i] * 0.5)
      alpha = (sigma_G[i,i] - sigma_new) / 2 + 0.5 * (log(sigma_new) - log(sigma_G[i,i])) # B_sigma = 1
      temp = log(runif(1))
      if (alpha > temp){
        sigma_G[i,i] <- sigma_new
      }
    }

    # Sample vol
    ytilde <- (A %*% (yt - B %*% xt) - Gamma_T * w) / w_sqrt
    aux <- sample_h_ele(ytilde = ytilde, sigma_h = sigma_h, h0_mean = h0_mean, h = h, K = K, t_max = t_max, prior = prior)
    h <- aux$Sigtdraw
    h0 <- as.numeric(aux$h0)
    sqrtvol <- aux$sigt
    sigma_h <- aux$sigma_h

    # Sample A0
    u_std <- (yt - B %*% xt) # change from Gaussian
    u_neg <- - u_std
    a_sample <- rep(0, K * (K - 1) /2)
    for (i in c(2:K)){
      id_end <- i*(i-1)/2
      id_start <- id_end - i + 2
      a_sub <- a_prior[id_start:id_end]
      V_a_sub <- V_a_prior[id_start:id_end, id_start:id_end]
      a_sample[c(id_start:id_end)] <- sample_A_ele(ysub = (u_std[i,] - w[i,] * Gamma_T[i,]) / sqrtvol[i,] / w_sqrt[i,],
                                                   xsub = matrix(u_neg[1:(i-1),] / reprow(sqrtvol[i,],i-1) / reprow(w_sqrt[i,], i-1), nrow = i-1),
                                                   a_sub = a_sub,
                                                   V_a_sub = V_a_sub)
    }
    A_post <- matrix(0, nrow = K, ncol = K)
    A_post[upper.tri(A)] <- a_sample
    A <- t(A_post)
    diag(A) <- 1

    # Sample w
    u <- A %*% (yt - B %*%xt)
    for (i in c(1:t_max)){
      ran_stu = rt(K, df = 4)
      Sigma2_inv_t = 1/exp(h[,i])
      a_eq <- -0.5* (Gamma_T[,i]^2 * Sigma2_inv_t)
      b_eq <- -((0.5*1 + nu/2+1))
      c_eq =  0.5 * (u[,i]^2 * Sigma2_inv_t) + nu/2
      w_mode <- (- b_eq - sqrt(b_eq^2-4*a_eq*c_eq))/(2*a_eq)
      w_mode <- ifelse(abs(w_mode)<1e-5, - c_eq/b_eq, w_mode)
      mode_target <- log( w_mode )
      cov_target <- 1/( c_eq / w_mode - a_eq * w_mode) * 1.2 # scale by 1.2

      log_w_temp = mode_target + sqrt(cov_target)*ran_stu
      w_temp = exp(log_w_temp)

      num_mh =  mapply( dinvgamma, x = w_temp, shape = nu*0.5, rate = nu*0.5, log = T) +  # prior
        ( - 0.5 * 1 * log_w_temp - 0.5 * (u[,i] - Gamma_T[,i]*w_temp)^2 * Sigma2_inv_t / w_temp) - # posterior
        dt( (log_w_temp - mode_target)/sqrt(cov_target), df = 4, log = T) + log_w_temp # proposal
      denum_mh = mapply( dinvgamma, x = w[,i], shape = nu*0.5, rate = nu*0.5, log = T) +  # prior
        ( - 0.5 * 1 * log(w[,i]) - 0.5 * (u[,i] - Gamma_T[,i]*w[,i])^2 * Sigma2_inv_t / w[,i]) - # posterior
        dt( (log(w[,i]) - mode_target)/sqrt(cov_target), df = 4, log = T) + log(w[,i]) # proposal
      alpha = num_mh - denum_mh
      temp = log(runif(K))
      acre = alpha > temp
      w[,i] <- ifelse(acre, w_temp, w[,i])
      w_sqrt[,i] <- sqrt(w[,i])
      acount_w[i] <- acount_w[i] + mean(acre)
    }

    # Sample nu
    nu_temp = nu + exp(logsigma_nu)*rnorm(K)
    for (k in c(1:K)){
      if (nu_temp[k] > 4 && nu_temp[k] < 100){
        num_mh = dgamma(nu_temp[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu_temp[k]*0.5, rate = nu_temp[k]*0.5, log = T))
        denum_mh = dgamma(nu[k], shape = nu_gam_a, rate = nu_gam_b, log = T) +
          sum(dinvgamma(w[k,], shape = nu[k]*0.5, rate = nu[k]*0.5, log = T))
        alpha = num_mh - denum_mh;
        temp = log(runif(1));
        if (alpha > temp){
          nu[k] = nu_temp[k]
          acount_nu[k] = acount_nu[k] + 1
        }

      }
    }

    if(j %% batchlength == 0 ){
      for (jj in c(1:K)) {
        if (acount_nu[jj] > batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] + adaptamount(j %/% batchlength);
        }
        if (acount_nu[jj] < batchlength * TARGACCEPT){
          logsigma_nu[jj] = logsigma_nu[jj] - adaptamount(j %/% batchlength);
        }
        acount_nu[jj] = 0
      }
    }

    if ((j > inits$burnin) & (j %% inits$thin == 0))
      mcmc[, (j - inits$burnin) %/% inits$thin] <- c(b_sample, a_sample, nu, diag(sigma_h), h0, diag(sigma_G), gamma0, as.numeric(Gamma_T), as.numeric(h), as.numeric(w))
    if (j %% 1000 == 0) {
      cat(" Iteration ", j, " ", logsigma_nu," ", min(acount_w)," ", max(acount_w)," ", mean(acount_w), " ", round(nu,2) , " ",
          round(apply(Gamma_T,MARGIN = 1,FUN = mean),2) ,  " \n")
      acount_w <- rep(0,t_max)
    }
  }
  nameA <- matrix(paste("a", reprow(c(1:K),K), repcol(c(1:K),K), sep = "_"), ncol = K)
  nameA <- nameA[upper.tri(nameA, diag = F)]
  row.names(mcmc) <- c( paste("B0",c(1:K), sep = ""),
                        sprintf("B%d_%d_%d",reprow(c(1:p),K*K), rep(repcol(c(1:K),K), p), rep(reprow(c(1:K),K), p)),
                        nameA,
                        paste("nu",c(1:K), sep = ""),
                        paste("sigma_h",c(1:K), sep = ""),
                        paste("lh0",c(1:K), sep = ""),
                        paste("sigma_G",c(1:K), sep = ""),
                        paste("G0",c(1:K), sep = ""),
                        sprintf("gamma_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("h_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)),
                        sprintf("w_%d_%d", repcol(c(1:K),t_max), reprow(c(1:t_max),K)))

  return(as.mcmc(t(mcmc)))
}