#' Transform a VEC Model to a VAR in Levels
#'
#' An object of class \code{"bvec"} is transformed to a VAR in level representation.
#'
#' @param object an object of class \code{"bvec"}.
#'
#' @return An object of class \code{"bvar"}.
#'
#' @examples
#'
#' # Load data
#' data("e6")
#'
#' # Generate model
#' data <- gen_vec(e6, p = 4, r = 1, const = "unrestricted", season = "unrestricted")
#'
#' # Obtain data matrices
#' y <- t(data$data$Y)
#' w <- t(data$data$W)
#' x <- t(data$data$X)
#'
#' # Reset random number generator for reproducibility
#' set.seed(1234567)
#'
#' iterations <- 100 # Number of iterations of the Gibbs sampler
#' # Chosen number of iterations should be much higher, e.g. 30000.
#'
#' burnin <- 100 # Number of burn-in draws
#' draws <- iterations + burnin
#'
#' r <- 1 # Set rank
#'
#' tt <- ncol(y) # Number of observations
#' k <- nrow(y) # Number of endogenous variables
#' k_w <- nrow(w) # Number of regressors in error correction term
#' k_x <- nrow(x) # Number of differenced regressors and unrestrictec deterministic terms
#'
#' k_alpha <- k * r # Number of elements in alpha
#' k_beta <- k_w * r # Number of elements in beta
#' k_gamma <- k * k_x
#'
#' # Set uninformative priors
#' a_mu_prior <- matrix(0, k_x * k) # Vector of prior parameter means
#' a_v_i_prior <- diag(0, k_x * k) # Inverse of the prior covariance matrix
#'
#' v_i <- 0
#' p_tau_i <- diag(1, k_w)
#'
#' u_sigma_df_prior <- r # Prior degrees of freedom
#' u_sigma_scale_prior <- diag(0, k) # Prior covariance matrix
#' u_sigma_df_post <- tt + u_sigma_df_prior # Posterior degrees of freedom
#'
#' # Initial values
#' beta <- matrix(c(1, -4), k_w, r)
#' u_sigma_i <- diag(1 / .0001, k)
#' g_i <- u_sigma_i
#'
#' # Data containers
#' draws_alpha <- matrix(NA, k_alpha, iterations)
#' draws_beta <- matrix(NA, k_beta, iterations)
#' draws_pi <- matrix(NA, k * k_w, iterations)
#' draws_gamma <- matrix(NA, k_gamma, iterations)
#' draws_sigma <- matrix(NA, k^2, iterations)
#'
#' # Start Gibbs sampler
#' for (draw in 1:draws) {
#' # Draw conditional mean parameters
#' temp <- post_coint_kls(y = y, beta = beta, w = w, x = x, sigma_i = u_sigma_i,
#' v_i = v_i, p_tau_i = p_tau_i, g_i = g_i,
#' gamma_mu_prior = a_mu_prior,
#' gamma_v_i_prior = a_v_i_prior)
#' alpha <- temp$alpha
#' beta <- temp$beta
#' Pi <- temp$Pi
#' gamma <- temp$Gamma
#'
#' # Draw variance-covariance matrix
#' u <- y - Pi %*% w - matrix(gamma, k) %*% x
#' u_sigma_scale_post <- solve(tcrossprod(u) +
#' v_i * alpha %*% tcrossprod(crossprod(beta, p_tau_i) %*% beta, alpha))
#' u_sigma_i <- matrix(rWishart(1, u_sigma_df_post, u_sigma_scale_post)[,, 1], k)
#' u_sigma <- solve(u_sigma_i)
#'
#' # Update g_i
#' g_i <- u_sigma_i
#'
#' # Store draws
#' if (draw > burnin) {
#' draws_alpha[, draw - burnin] <- alpha
#' draws_beta[, draw - burnin] <- beta
#' draws_pi[, draw - burnin] <- Pi
#' draws_gamma[, draw - burnin] <- gamma
#' draws_sigma[, draw - burnin] <- u_sigma
#' }
#' }
#'
#' # Number of non-deterministic coefficients
#' k_nondet <- (k_x - 4) * k
#'
#' # Generate bvec object
#' bvec_est <- bvec(y = data$data$Y, w = data$data$W,
#' x = data$data$X[, 1:6],
#' x_d = data$data$X[, 7:10],
#' Pi = draws_pi,
#' Gamma = draws_gamma[1:k_nondet,],
#' C = draws_gamma[(k_nondet + 1):nrow(draws_gamma),],
#' Sigma = draws_sigma)
#'
#' # Thin posterior draws
#' bvec_est <- thin(bvec_est, thin = 5)
#'
#' # Transfrom VEC output to VAR output
#' bvar_form <- bvec_to_bvar(bvec_est)
#'
#' @references
#'
#' Lütkepohl, H. (2006). \emph{New introduction to multiple time series analysis} (2nd ed.). Berlin: Springer.
#'
#' @export
bvec_to_bvar <- function(object) {
if (!any(class(object) %in% "bvec")) {
stop("Argument 'object' must be of class 'bvec'.")
}
draws <- NULL
specs <- NULL
vars <- c("Pi", "Pi_x", "Pi_d", "Gamma", "Upsilon", "C", "A0")
for (i in vars) {
if (is.null(draws)) {
if (!is.null(object[[i]])) {
if (is.list(object[[i]])) {
draws <- nrow(object[[i]][[1]])
} else {
draws <- nrow(object[[i]])
}
}
}
if (is.null(specs)) {
if (is.list(object[[i]])) {
specs <- attr(object[[i]][[1]], "mcpar")
} else {
specs <- attr(object[[i]], "mcpar")
}
}
}
k <- NCOL(object[["y"]])
tt <- NROW(object[["y"]])
tvp <- object[["specifications"]][["tvp"]]
p <- object[["specifications"]][["lags"]][["p"]]
r <- object[["specifications"]][["rank"]]
A <- NULL
if (!is.null(object[["Gamma"]])) {
W <- diag(-1, k * p)
W[1:k, 1:k] <- diag(1, k)
W[-(1:k), -(k * (p - 1) + 1:k)] <- W[-(1:k),-(k * (p - 1) + 1:k)] + diag(k * (p - 1))
J <- matrix(0, k, k * p)
J[1:k, 1:k] <- diag(1, k)
n_gamma <- k * k * p
if (!tvp[["Pi"]] & !tvp[["Gamma"]]) {
A <- matrix(NA, n_gamma, draws)
for (draw in 1:draws) {
if (is.null(object[["Pi"]])) {
A[, draw] <- cbind(matrix(0, k, k), matrix(object[["Gamma"]][draw, ], k)) %*% W + J
} else {
A[, draw] <- cbind(matrix(object[["Pi"]][draw, ], k), matrix(object[["Gamma"]][draw, ], k)) %*% W + J
}
}
} else {
A <- matrix(NA, n_gamma * tt, draws)
pi_temp <- matrix(0, k, k)
for (draw in 1:draws) {
for (i in 1:tt) {
if (!is.null(object[["Pi"]])) {
if (tvp[["Pi"]]) {
pi_temp <- matrix(object[["Pi"]][[i]][draw, ], k)
} else {
pi_temp <- matrix(object[["Pi"]][draw, ], k)
}
}
if (tvp[["Gamma"]]) {
gamma_temp <- matrix(object[["Gamma"]][[i]][draw, ], k)
} else {
gamma_temp <- matrix(object[["Gamma"]][draw, ], k)
}
A[(i - 1) * n_gamma + 1:n_gamma, draw] <- cbind(pi_temp, gamma_temp) %*% W + J
}
}
}
} else {
if (!is.null(object[["Pi"]])) {
n_a <- k * k
if (tvp[["Pi"]]) {
A <- matrix(NA, n_a * tt, draws)
for (draw in 1:draws) {
for (i in 1:tt) {
A[(i - 1) * n_a + 1:n_a, draw] <- matrix(object[["Pi"]][[i]][draw, ], k) + matrix(diag(1, k), k)
}
}
} else {
A <- matrix(NA, n_a, draws)
for (draw in 1:draws) {
A[, draw] <- matrix(object[["Pi"]][draw, ], k) + matrix(diag(1, k), k)
}
}
} else {
A <- matrix(0, k * k, draws)
}
}
B <- NULL
m <- 0
s <- 0
if (!is.null(object[["Upsilon"]])) {
m <- object[["specifications"]][["dims"]][["M"]]
s <- object[["specifications"]][["lags"]][["s"]]
W <- diag(-1, m * (s + 1))
W[1:m, 1:m] <- 0
W[1:m, m + 1:m] <- diag(1, m)
W[-(1:m), 1:(m * s)] <- W[-(1:m), 1:(m * s)] + diag(1, m * s)
n_b <- k * m * (s + 1)
if (!tvp[["Pi_x"]] & !tvp[["Upsilon"]]) {
B <- matrix(NA, n_b, draws)
for (draw in 1:draws){
if (!is.null(object[["Pi_x"]])) {
B[, draw] <- cbind(matrix(object[["Pi_x"]][draw, ], k), matrix(object[["Upsilon"]][draw, ], k)) %*% W
} else {
B[, draw] <- cbind(matrix(0, k, m), matrix(object[["Upsilon"]][draw, ], k)) %*% W
}
}
} else {
B <- matrix(NA, n_b * tt, draws)
pix_temp <- matrix(0, k, m)
for (draw in 1:draws){
for (i in 1:tt) {
if (!is.null(object[["Pi_x"]])) {
if (tvp[["Pi_x"]]) {
pix_temp <- matrix(object[["Pi_x"]][[i]][draw, ], k)
} else {
pix_temp <- matrix(object[["Pi_x"]][draw, ], k)
}
}
if (tvp[["Upsilon"]]) {
ups_temp <- matrix(object[["Upsilon"]][[i]][draw, ], k)
} else {
ups_temp <- matrix(object[["Upsilon"]][draw, ], k)
}
B[(i - 1) * n_b + 1:n_b, draw] <- cbind(pix_temp, ups_temp) %*% W
}
}
}
}
y <- object[["y"]]
y_names <- gsub("d.", "", dimnames(y)[[2]])
dimnames(y) <- list(NULL, y_names)
# Reconstruct y if w is available
if (!is.null(object[["w"]])) {
y <- y + as.matrix(object[["w"]][, 1:k])
dimnames(y) <- list(NULL, y_names)
}
tsp_temp <- stats::tsp(y)
if (p > 1) {
y_diff <- matrix(object[["x"]][1,], k) # Get first observations
y_diff <- matrix(y_diff[,(p - 1):1], k) # Reverse order
y_diff <- cbind(y_diff, matrix(object[["y"]][1, ], k))
y_level <- y_diff * NA
y_level[, p] <- y[1, ] - y_diff[, p]
for (i in (p - 1):1){
y_level[, i] <- y_level[, i + 1] - y_diff[, i]
}
y <- rbind(t(y_level), y)
} else {
if (!is.null(object[["w"]])) {
y <- rbind(object[["w"]][1, ], y)
}
}
if (m > 0) {
x_names <- dimnames(object[["x_x"]])[[2]][1:m]
x_names <- gsub("d.", "", x_names)
x_names <- substring(x_names, 1, unlist(gregexpr(".l", x_names))[1] - 1)
if (s > 0) {
x0 <- matrix(object[["x_x"]][, 1:m], NROW(object[["x_x"]]))
x0 <- x0 + object[["w_x"]]
x_diff <- matrix(object[["x_x"]][1, ], m)
x_diff <- matrix(x_diff[, s:1], m)
x_level <- x_diff * NA
x_level[, s] <- x0[1, ] - x_diff[, s]
if (s > 1) {
for (i in (s - 1):1){
x_level[, i] <- x_level[, i + 1] - x_diff[, i]
}
}
x <- rbind(t(x_level), x0)
}
} else {
x <- NULL
}
if (!is.null(x)) {
if (NROW(y) < NROW(x)) {
y <- rbind(matrix(NA, NROW(x) - NROW(y), k), y)
}
if (NROW(y) > NROW(x)) {
x <- rbind(matrix(NA, NROW(y) - NROW(x), m), x)
}
}
y <- stats::ts(y)
if (!is.null(x)) {
x <- stats::ts(x)
}
# Generate bvar data matrices
temp <- y
temp_names <- y_names
if (!is.null(p)) {
for (i in 1:p) {
temp <- cbind(temp, stats::lag(y, -i))
if (nchar(p) > 2) {
i_temp <- paste0(c(rep(0, nchar(p) - nchar(i)), i), collapse = "")
} else {
i_temp <- paste0(c(rep(0, 2 - nchar(i)), i), collapse = "")
}
temp_names <- c(temp_names, paste0(y_names, ".l", i_temp))
}
}
if (m > 0) {
temp <- cbind(temp, x)
if (nchar(s) > 2) {
i_temp <- paste0(rep(0, nchar(s)), collapse = "")
} else {
i_temp <- paste0(rep(0, 2), collapse = "")
}
temp_names <- c(temp_names, paste0("exo.", x_names, ".l", i_temp))
for (i in 1:s) {
temp <- cbind(temp, stats::lag(x, -i))
if (nchar(s) > 2) {
i_temp <- paste0(c(rep(0, nchar(s) - nchar(i)), i), collapse = "")
} else {
i_temp <- paste0(c(rep(0, 2 - nchar(i)), i), collapse = "")
}
temp_names <- c(temp_names, paste0("exo.", x_names, ".l", i_temp))
}
}
temp <- stats::na.omit(temp)
dimnames(temp)[[2]] <- temp_names
y <- stats::ts(as.matrix(temp[, 1:k]), class = c("mts", "ts", "matrix"),
end = tsp_temp[2], frequency = tsp_temp[3])
dimnames(y) <- list(NULL, y_names)
x <- NULL
if (p > 0 | m > 0) {
x <- stats::ts(as.matrix(temp[, -(1:k)]), class = c("mts", "ts", "matrix"),
end = tsp_temp[2], frequency = tsp_temp[3])
dimnames(x) <- list(NULL, temp_names[-(1:k)])
}
if (!is.null(object[["A0"]])) {
if (tvp[["A0"]]) {
n_a0 <- k * k
A0 <- matrix(NA, n_a0 * tt, draws)
for (i in 1:tt) {
for (draw in 1:draws) {
A0[(i - 1) * n_a0 + 1:n_a0, draw] <- object[["A0"]][[i]][draw, ]
}
}
} else {
A0 <- t(object[["A0"]])
}
} else {
A0 <- NULL
}
n <- 0
n_Pi_d <- 0
C <- NULL
x_det <- NULL
x_det_names <- NULL
if (!is.null(object[["Pi_d"]])) {
if (tvp[["Pi_d"]]) {
n <- NCOL(object[["Pi_d"]][[1]]) / k
} else {
n <- NCOL(object[["Pi_d"]]) / k
}
n_Pi_d <- n
if (tvp[["Pi_d"]]) {
n_pi <- k * n_Pi_d
C <- matrix(NA, n_pi * tt, draws)
for (i in 1:tt) {
for (draw in 1:draws) {
C[(i - 1) * n_pi + 1:n_pi, draw] <- object[["Pi_d"]][[i]][draw, ]
}
}
} else {
C <- t(object[["Pi_d"]])
}
if (!is.null(object[["w_d"]])) {
x_det_names <- dimnames(object[["w_d"]])[[2]]
x_det <- object[["w_d"]]
}
}
n_c <- 0
C_ur <- NULL
if (!is.null(object[["C"]])) {
if (tvp[["C"]]) {
n_c <- NCOL(object[["C"]][[1]]) / k
} else {
n_c <- NCOL(object[["C"]]) / k
}
n <- n + n_c
if (tvp[["C"]]) {
n_c_ur <- k * n_c
C_ur <- matrix(NA, n_c_ur * tt, draws)
for (i in 1:tt) {
for (draw in 1:draws) {
C_ur[(i - 1) * n_c_ur + 1:n_c_ur, draw] <- object[["C"]][[i]][draw, ]
}
}
} else {
C_ur <- t(object[["C"]])
}
if (!is.null(C)) {
C <- rbind(C, C_ur)
} else {
C <- C_ur
}
}
if (!is.null(object[["x_d"]]) & n_c > 0) {
x_det_names <- c(x_det_names, dimnames(object[["x_d"]])[[2]])
x_temp <- object[["x_d"]]
if (is.null(x_det)) {
x_det <- x_temp
} else {
x_det <- cbind(x_det, x_temp)
}
}
if (!is.null(x_det)) {
x_names <- c(dimnames(x)[[2]], x_det_names)
x <- cbind(x, x_det)
dimnames(x)[[2]] <- x_names
}
if (!is.null(object[["Sigma"]])) {
if (tvp[["Sigma"]]) {
n_sigma <- k * k
Sigma <- matrix(NA, n_sigma * tt, draws)
for (i in 1:tt) {
for (draw in 1:draws) {
Sigma[(i - 1) * n_sigma + 1:n_sigma, draw] <- object[["Sigma"]][[i]][draw, ]
}
}
} else {
Sigma <- t(object[["Sigma"]])
}
} else {
Sigma <- NULL
}
if (!is.null(object[["data"]])) {
data <- object[["data"]]
} else {
data <- NULL
}
if (!is.null(object[["exogen"]])) {
exogen <- object[["exogen"]]
} else {
exogen <- NULL
}
object <- bvar(data = data, exogen = exogen, y = y, x = x, A0 = A0, A = A, B = B, C = C, Sigma = Sigma)
vars <- c("A", "B", "C", "Sigma", "A0")
for (i in vars) {
if(!is.null(object[[i]])) {
if (object[["specifications"]][["tvp"]][[i]]) {
for (j in 1:tt) {
attr(object[[i]][[j]], "mcpar") <- specs
}
} else {
attr(object[[i]], "mcpar") <- specs
}
}
}
return(object)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.