Pawel's Helper Functions

The original file was written in R, thus syntax follows the same convention:

Author: Pawel Szerszen (pawel.j.szerszen@frb.gov)

Summary Statistics

mean_exp <- function(x, log_weights, scaled = FALSE) {
    max_x <- max(x)
    w <- exp(log_weights-max(log_weights))
    mean_exp_scaled <- sum(w * exp(x - max_x)) / sum(w)
    if (scaled == FALSE) {
        res_out <- mean_exp_scaled * exp(max_x)
    } else {
        res_out <- mean_exp_scaled
    }
    return(res_out)
}

Sigma_exp <- function(x, log_weights) {
    max_x      <- max(x)
    w <- exp(log_weights-max(log_weights))
    Sigma_x <- cov.wt(x = x, wt = w)
    return(Sigma_x)
}

Effective Sample Size

ESS_calc <- function(logx) {
    x_scaled <- exp(logx - max(logx))
    ESS <- sum(x_scaled)^2 / sum(x_scaled^2)
    if(is.nan(ESS)) ESS <- 0
    return(ESS)
}

ESS_calc_mat <- function(logx) {
    x_scaled <- exp(logx - matrix(rep(apply(logx,2,max),nrow(logx)),ncol=ncol(logx), byrow = T))
    ESS <- colSums(x_scaled)^2 / colSums(x_scaled^2)
    if(is.nan(ESS)) ESS <- 0
    return(ESS)
}

Truncated Distributions

# draw from proposal denisty of transition to x_prop
r_trunc_mult_norm <- function(mean_par, cov_par, lowerbnd, upperbnd, fixpar = NULL) {
    Sigma_list = lapply(
        1:length(mean_par),
        function(j_par) {
            Sig_mix = solve(cov_par[-j_par, -j_par], cov_par[-j_par, j_par])
            cov_par_cond = cov_par[j_par,j_par] - cov_par[j_par, -j_par] %*% Sig_mix
            return(list(Sig_mix = Sig_mix, sigma_par_cond = sqrt(cov_par_cond)))
        }
    )

    #initialize the algorithm
    x <- mean_par
    for (j_par in (1:length(mean_par))) {
        if (!(j_par %in% fixpar)) {
            mean_par_cond  = mean_par[j_par] + t(Sigma_list[[j_par]]$Sig_mix) %*% (x[-j_par] - mean_par[-j_par])
            lowerbnd_scaled = (lowerbnd[[j_par]] - mean_par_cond) / Sigma_list[[j_par]]$sigma_par_cond
            upperbnd_scaled = (upperbnd[[j_par]] - mean_par_cond) / Sigma_list[[j_par]]$sigma_par_cond
            x[j_par] = qnorm(runif(1, min = pnorm(lowerbnd_scaled), max = pnorm(upperbnd_scaled))) * Sigma_list[[j_par]]$sigma_par_cond + mean_par_cond 
        }
    }
    return(x)
}

# evaluate proposal denisty of transition to x_prop
D_trunc_mult_norm <- function(x_prop, mean_par, cov_par, lowerbnd, upperbnd, fixpar = NULL) {
    Sigma_list = lapply(
        1:length(mean_par),
        function(j_par) {
            Sig_mix = solve(cov_par[-j_par, -j_par], cov_par[-j_par, j_par])
            cov_par_cond = cov_par[j_par,j_par] - cov_par[j_par, -j_par] %*% Sig_mix
            return(list(Sig_mix = Sig_mix, sigma_par_cond = sqrt(cov_par_cond)))
        }
    )

    #initialize the algorithm
    loglik <- 0
    x_current <- mean_par
    for (j_par in (1:length(mean_par))) {
        if (!(j_par %in% fixpar)) {
            mean_par_cond  = mean_par[j_par] + t(Sigma_list[[j_par]]$Sig_mix) %*% (x_current[-j_par] - mean_par[-j_par])
            lowerbnd_in = lowerbnd[[j_par]]
            upperbnd_in = upperbnd[[j_par]]
            loglik <- loglik + dnorm(x_prop[j_par], mean = mean_par_cond, sd = Sigma_list[[j_par]]$sigma_par_cond, log = T) - log(pnorm(upperbnd_in) - pnorm(lowerbnd_in))
            x_current[j_par]    = x_prop[j_par]
        }
    }
    return(loglik)
}