Skip to content

Commit

Permalink
Better handling of unsuccessful Gibbs samplers
Browse files Browse the repository at this point in the history
  • Loading branch information
@franzmohr
franzmohr committed Jul 7, 2023
1 parent 2f78384 commit eb23a4e
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 101 deletions.
15 changes: 6 additions & 9 deletions R/select_submodels.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@
#' specified information criterion per country is selected. If \code{"rank"}, the
#' function selects the model, after which the selected information criterion increases
#' for the first time.
#' @param teststats optional. A named list of test statistics for each model in `object`. Usually,
#' the result of a call to \code{\link{submodel_test_statistics}[["teststats"]]}.
#' @param teststats optional. An object of class "ctryvartest", usually,
#' the result of a call to \code{\link{submodel_test_statistics}}.
#'
#' @examples
#' # Load data
Expand Down Expand Up @@ -90,14 +90,11 @@ select_submodels <- function(object, ic = "BIC", select = "order", teststats = N

# Obtain test statistics
if (is.null(teststats)) {
criteria <- submodel_test_statistics(object)[["teststats"]]
criteria <- submodel_test_statistics(object)
} else {
# Basic checks
if (!"data.frame" %in% class(teststats)) {
stop("If provided, argument 'teststats' must be of class data.frame.")
}
if (length(teststats) != length(unique(names(object)))) {
stop("Number of countries in argument 'object' differs from number of countries in argument 'teststats'.")
if (!"ctryvartest" %in% class(teststats)) {
stop("If provided, argument 'teststats' must be of class ctryvartest.")
}
criteria <- teststats
}
Expand All @@ -107,7 +104,7 @@ select_submodels <- function(object, ic = "BIC", select = "order", teststats = N
crit <- list()
ctry <- unique(names(object))
for (i in 1:length(ctry)) {
crit[[i]] <- criteria[criteria[, "ctry"] == ctry[i],]
crit[[i]] <- criteria[[i]][["teststats"]]
}
names(crit) <- ctry
criteria <- crit
Expand Down
181 changes: 91 additions & 90 deletions R/submodel_test_statistics.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,20 +107,11 @@ submodel_test_statistics <- function(object, ...){
current_ctry <- names(object)[i]
n_models_curr_ctry <- sum(names(object) == current_ctry)
n_completed_curr_ctry <- sum(names(object)[1:i] == current_ctry)
update_result_element <- n_models_curr_ctry == n_completed_curr_ctry
pos <- cumsum(as.numeric(names(object) == current_ctry))
pos[as.numeric(names(object) == current_ctry) == 0] <- 0

teststats[i, "ctry"] <- names(object)[i]

# Skip tests if posterior simulation was not successful
if (!is.null(object[[i]][["error"]])) {
if (object[[i]][["error"]]) {
loglik[[pos]] <- NA
next
}
}

# Get model specifications
structural <- object[[i]][["model"]][["structural"]]
tvp <- object[[i]][["model"]][["tvp"]]
Expand All @@ -144,113 +135,123 @@ submodel_test_statistics <- function(object, ...){
}
if (!is.null(object[[i]][["model"]][["rank"]])) {
teststats[i, "r"] <- object[[i]][["model"]][["rank"]]
}
}

if (tvp) {
temp_pars <- list()
length(temp_pars) <- tt
# Skip tests if posterior simulation was not successful
if (!is.null(object[[i]][["error"]])) {
if (object[[i]][["error"]]) {
loglik[[pos[i]]] <- NA
next
}
} else {
temp_pars <- NULL
}
x <- NULL

if ("ctryvarest" %in% class(object[[i]])) {

x <- t(object[[i]][["data"]][["Z"]])
tot_pars <- NCOL(object[[i]][["data"]][["Z"]])

vars <- c("domestic", "foreign", "global", "deterministic")
for (j in vars) {
if (!is.null(object[[i]][["posteriors"]][[j]])) {
if (is.list(object[[i]][["posteriors"]][[j]])) {
for (period in 1:tt) {
temp_pars[[period]] <- cbind(temp_pars[[period]], object[[i]][["posteriors"]][[j]][[period]])
if (tvp) {
temp_pars <- list()
length(temp_pars) <- tt
} else {
temp_pars <- NULL
}
x <- NULL

if ("ctryvarest" %in% class(object[[i]])) {

x <- t(object[[i]][["data"]][["Z"]])
tot_pars <- NCOL(object[[i]][["data"]][["Z"]])

vars <- c("domestic", "foreign", "global", "deterministic")
for (j in vars) {
if (!is.null(object[[i]][["posteriors"]][[j]])) {
if (is.list(object[[i]][["posteriors"]][[j]])) {
for (period in 1:tt) {
temp_pars[[period]] <- cbind(temp_pars[[period]], object[[i]][["posteriors"]][[j]][[period]])
}
} else {
temp_pars <- cbind(temp_pars, object[[i]][["posteriors"]][[j]])
}
} else {
temp_pars <- cbind(temp_pars, object[[i]][["posteriors"]][[j]])
}
}
}
}

if ("ctryvecest" %in% class(object[[i]])) {

object[[i]] <- .create_pi_matrices(object[[i]])

x <- t(object[[i]][["data"]][["X"]])
tot_pars <- object[[i]][["model"]][["rank"]] + nrow(x)
if (object[[i]][["model"]][["rank"]] > 0) {
x <- rbind(t(object[[i]][["data"]][["W"]]), x)
}

vars <- c("pi_domestic", "pi_foreign", "pi_global", "pi_deterministic",
"gamma_domestic", "gamma_foreign", "gamma_global", "gamma_deterministic")
for (j in vars) {
if (!is.null(object[[i]][["posteriors"]][[j]])) {
if (is.list(object[[i]][["posteriors"]][[j]])) {
for (period in 1:tt) {
temp_pars[[period]] <- cbind(temp_pars[[period]], object[[i]][["posteriors"]][[j]][[period]])
if ("ctryvecest" %in% class(object[[i]])) {

object[[i]] <- .create_pi_matrices(object[[i]])

x <- t(object[[i]][["data"]][["X"]])
tot_pars <- object[[i]][["model"]][["rank"]] + nrow(x)
if (object[[i]][["model"]][["rank"]] > 0) {
x <- rbind(t(object[[i]][["data"]][["W"]]), x)
}

vars <- c("pi_domestic", "pi_foreign", "pi_global", "pi_deterministic",
"gamma_domestic", "gamma_foreign", "gamma_global", "gamma_deterministic")
for (j in vars) {
if (!is.null(object[[i]][["posteriors"]][[j]])) {
if (is.list(object[[i]][["posteriors"]][[j]])) {
for (period in 1:tt) {
temp_pars[[period]] <- cbind(temp_pars[[period]], object[[i]][["posteriors"]][[j]][[period]])
}
} else {
temp_pars <- cbind(temp_pars, object[[i]][["posteriors"]][[j]])
}
} else {
temp_pars <- cbind(temp_pars, object[[i]][["posteriors"]][[j]])
}
}

}

}

if (tvp) {
draws <- nrow(temp_pars[[1]])
} else {
draws <- nrow(temp_pars)
}
loglik[[pos[i]]] <- matrix(NA, tt, draws)
y <- t(object[[i]][["data"]][["Y"]])
u <- y * 0
if (sv) {
sigma <- matrix(NA_real_, k_domestic * tt, k_domestic)
} else {
sigma <- matrix(NA_real_, k_domestic, k_domestic)
}

for (j in 1:draws) {
# Residuals
if (tvp) {
for (period in 1:tt) {
draws <- nrow(temp_pars[[1]])
} else {
draws <- nrow(temp_pars)
}
loglik[[pos[i]]] <- matrix(NA, tt, draws)
y <- t(object[[i]][["data"]][["Y"]])
u <- y * 0
if (sv) {
sigma <- matrix(NA_real_, k_domestic * tt, k_domestic)
} else {
sigma <- matrix(NA_real_, k_domestic, k_domestic)
}

for (j in 1:draws) {
# Residuals
if (tvp) {
for (period in 1:tt) {
if (structural) {
A0 <- matrix(object[[i]][["posteriors"]][["a0"]][[period]][j, ], k_domestic)
} else {
A0 <- diag(k_domestic)
}
u[, period] <- y[, period] - matrix(temp_pars[[period]][j, ], k_domestic) %*% x[, period]
}
} else {
if (structural) {
A0 <- matrix(object[[i]][["posteriors"]][["a0"]][[period]][j, ], k_domestic)
A0 <- matrix(object[[i]][["posteriors"]][["a0"]][j, ], k_domestic)
} else {
A0 <- diag(k_domestic)
}
u[, period] <- y[, period] - matrix(temp_pars[[period]][j, ], k_domestic) %*% x[, period]
u <- A0 %*% y - matrix(temp_pars[j, ], k_domestic) %*% x
}
} else {
if (structural) {
A0 <- matrix(object[[i]][["posteriors"]][["a0"]][j, ], k_domestic)

if (sv) {
for (period in 1:tt) {
sigma[(period - 1) * k_domestic + 1:k_domestic,] <- matrix(object[[i]][["posteriors"]][["sigma"]][[period]][j,], k_domestic)
}
} else {
A0 <- diag(k_domestic)
sigma <- matrix(object[[i]][["posteriors"]][["sigma"]][j,], k_domestic)
}
u <- A0 %*% y - matrix(temp_pars[j, ], k_domestic) %*% x

# Calculate log-likelihood
loglik[[pos[i]]][, j] <- bvartools::loglik_normal(u, sigma)
}

if (sv) {
for (period in 1:tt) {
sigma[(period - 1) * k_domestic + 1:k_domestic,] <- matrix(object[[i]][["posteriors"]][["sigma"]][[period]][j,], k_domestic)
}
} else {
sigma <- matrix(object[[i]][["posteriors"]][["sigma"]][j,], k_domestic)
}

# Calculate log-likelihood
loglik[[pos[i]]][, j] <- bvartools::loglik_normal(u, sigma)
ll <- sum(rowMeans(loglik[[pos[i]]]))
teststats[i, "LL"] <- ll
teststats[i, "AIC"] <- 2 * tot_pars - 2 * ll
teststats[i, "BIC"] <- log(tt) * tot_pars - 2 * ll
teststats[i, "HQ"] <- 2 * log(log(tt)) * tot_pars - 2 * ll
}

ll <- sum(rowMeans(loglik[[pos[i]]]))
teststats[i, "LL"] <- ll
teststats[i, "AIC"] <- 2 * tot_pars - 2 * ll
teststats[i, "BIC"] <- log(tt) * tot_pars - 2 * ll
teststats[i, "HQ"] <- 2 * log(log(tt)) * tot_pars - 2 * ll

# If all model of a country are finished, updated the result object
if (n_models_curr_ctry == n_completed_curr_ctry) {
model_pos <- which(teststats[, "ctry"] == current_ctry)
Expand Down
4 changes: 2 additions & 2 deletions man/select_submodels.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit eb23a4e

Please sign in to comment.