R/bdfa.R

In vsartor/bnsa: Bayesian Non-Structured Analysis

#' Bayesian Dynamic Factor Analysis
#'
#' Simulates from the posterior distribution of the parameters of a
#' simple Bayesian Dynamic Factor Analysis using the Gibbs algorithm,
#' considering vague prior distributions for all parameters and a random
#' walk evolution for the factors.
#'
#' @param Y A T by m matrix containing each m-dimensional observation
#'          in each row.
#' @param k The number of factors to be considered.
#' @param N The number of iterations for the MCMC algorithm to run.
#'
#' @return A list with the chains `factor`, `beta`, `vars` and `lambda`.
#'
#' @importFrom stats rgamma
#' @importFrom MASS mvrnorm
#'
#' @export
#'
bdfa = function(Y, k, N = 5000) {
    # Dimensions of the problem
    n = nrow(Y)
    m = ncol(Y)

    # List to be returned
    ch         = list()
    ch$factors = array(0, dim = c(n, k, N))
    ch$beta    = array(0, dim = c(m, k, N))
    ch$vars    = array(0, dim = c(m, N))
    ch$lambda  = array(0, dim = c(k, N))

    # Hyperparameters
    nu = 1
    s  = 1
    m0 = 0
    C0 = 1e8

    # Set up non-zero initial values
    if (k == 1) {
        ch$beta[ ,1,1] = 1
    } else {
        diag(ch$beta[ , ,1]) = 1
    }
    ch$vars[ ,1]   = 1
    ch$lambda[ ,1] = 1

    # Initialize loop variables
    factors = matrix(ch$factors[ , ,1], nrow = n)
    beta    = matrix(ch$beta[ , ,1], nrow = m)
    vars    = ch$vars[ ,1]
    lambda  = ch$lambda[ ,1]

    # Initialize constants
    I = vector("list", k)
    for (i in 1:k) I[[i]] = diag(i)

    for (it in 2:N) {
        # Update user on progress
        if (it %% 250 == 0)
            cat(sprintf("[%d/%d] Computing...\n", it, N))

        # Draw for the factors
        kfs = kalmanfs(Y, beta, I[[k]], diag(vars), diag(lambda))
        for (j in 1:n) {
            # Positivize covariance matrix
            tmp = eigen(kfs$S[ , ,j])
            V = tmp$vectors
            L = ifelse(tmp$values <= 0, 0.0001, tmp$values)
            Cov = V %*% diag(L) %*% solve(V)
            # Actually simulate
            factors[j, ] = mvrnorm(1, kfs$s[ ,j], Cov)
        }
        ch$factors[ , ,it] = factors

        # Draw for the variances
        for (i in 1:m) {
            di      = sum((Y[ ,i] - factors %*% beta[i, ])^2)
            vars[i] = 1 / rgamma(1, (nu + n)/2, (nu*s + di)/2)
        }
        ch$vars[ ,it] = vars

        # Draw for the factor variances
        for (i in 1:k) {
            lambda.shape = n - 1 + 2e-5
            lambda.rate  = sum((factors[2:n,i] - factors[1:(n-1),i])^2) + 2e-5
            lambda[i]    = rgamma(1, lambda.shape / 2, lambda.rate / 2)
        }
        ch$lambda[ ,it] = lambda

        # Draw for loadings
        for (i in 1:k) {
            Fac         = factors[ ,1:i]
            Ci          = solve(I[[i]] / C0 + t(Fac) %*% Fac / vars[i])
            mi          = Ci %*% (m0 * rep(1, i) / C0 + t(Fac) %*% Y[ ,i] / vars[i])
            beta[i,1:i] = mvrnorm(1, mi, Ci)
            beta[i,i]   = 1
        }
        for (i in (k+1):m) {
            Fac       = factors
            Ci        = solve(I[[k]] / C0 + t(Fac) %*% Fac / vars[i])
            mi        = Ci %*% (m0 * rep(1, k) / C0 + t(Fac) %*% Y[ ,i] / vars[i])
            beta[i, ] = mvrnorm(1, mi, Ci)
        }
        ch$beta[ , ,it] = beta
    }

    # Return the chains
    ch
}