setwd("C:/Users/strohsal/Dropbox/Arbeit/Luetkepohl/2025_Non_normality_and_instruments/5_JAE/3_Second_Revision/4_Replication_Package/2_This_papers_replication/Code")
setwd("C:/Users/strohsal/Dropbox/Arbeit/Luetkepohl/2025_Non_normality_and_instruments/5_JAE/3_Second_Revision/4_Replication_Package/2_This_papers_replication/Code")
rm(list = ls()); set.seed(42)
# WARNING: Full run (1,000 outer x 100 inner draws) may take 10 days or more.
pkgs <- c("vars", "svars", "zoo")
missing <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if (length(missing))
stop("Please install R package(s): ", paste(missing, collapse = ", "))
if (!requireNamespace("rhdf5", quietly = TRUE))
stop("Please install 'rhdf5' via Bioconductor: BiocManager::install('rhdf5', update = FALSE)")
suppressPackageStartupMessages({
library(vars); library(svars); library(zoo); library(rhdf5)
})
project_dir <- dirname(getwd())
in_dir  <- file.path(project_dir, "Data", "kaenzig_exports")
out_dir <- file.path(project_dir, "Data", "thispaper_exports")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
matfile <- file.path(in_dir, "kaenzig_nested_bootstrap_exports.mat")
checkpoint <- file.path(out_dir, "thispaper_wald_checkpoint.rds")
outfile <- file.path(out_dir, "thispaper_wald_test_bK_equals_bNG_nested.csv")
# Baseline nG estimate
D <- read.csv(file.path(in_dir, "kaenzig_baseline_data.csv"),
na.strings = c("NaN", "NA", ""))
ym <- as.yearmon(gsub("M", "-", D$time), "%Y-%m")
Y <- D[, 2:7]
Y[] <- lapply(Y, as.numeric)
ok <- complete.cases(Y)
Y <- Y[ok, , drop = FALSE]
ym <- ym[ok]
Y_ts <- ts(as.matrix(Y),
start = c(as.integer(format(ym[1], "%Y")),
as.integer(format(ym[1], "%m"))), frequency = 12)
news_full <- read.csv(file.path(in_dir, "kaenzig_baseline_structural_shock.csv"))$shock
plain.var0 <- VAR(Y_ts, p = 12, type = "const")
u0 <- resid(plain.var0)
T0 <- nrow(u0)
news <- tail(news_full, T0)
B0 <- id.ngml(plain.var0)$B
eps0 <- t(solve(B0, t(u0)))
oil_j0 <- which.max(abs(apply(eps0, 2, function(x) cor(x, news))))
b_ngml0 <- as.numeric(B0[, oil_j0])
if (!is.na(b_ngml0[1]) && b_ngml0[1] < 0) b_ngml0 <- -b_ngml0
# Bootstrap metadata
p_mat    <- as.integer(h5read(matfile, "p"))
nvar_mat <- as.integer(h5read(matfile, "nvar"))
nsim     <- as.integer(h5read(matfile, "nsim"))
ninner   <- as.integer(h5read(matfile, "ninner"))
nobs_boot <- T0 + p_mat
match_column <- function(Bstar, bref) {
sims <- rep(NA_real_, ncol(Bstar))
for (j in seq_len(ncol(Bstar))) {
bj <- as.numeric(Bstar[, j])
sims[j] <- abs(sum(bj * bref) / sqrt(sum(bj^2) * sum(bref^2)))
}
which.max(sims)
}
estimate_b <- function(Ymat, bref) {
Ymat <- drop(Ymat)
Ymat <- as.matrix(Ymat)
if (ncol(Ymat) != length(bref) && nrow(Ymat) == length(bref)) Ymat <- t(Ymat)
storage.mode(Ymat) <- "double"
Bstar <- id.ngml(VAR(ts(Ymat, frequency = 12), p = p_mat, type = "const"))$B
j <- match_column(Bstar, bref)
b <- as.numeric(Bstar[, j])
if (sum(b * bref) < 0) b <- -b
b
}
save_checkpoint <- function() saveRDS(list(
bootb_ngml = bootb_ngml, bootb_ngml_inner = bootb_ngml_inner,
success = success, success_inner = success_inner,
done = done, done_inner = done_inner,
nsim = nsim, ninner = ninner, rng_state = .Random.seed
), checkpoint)
# Initialize/resume
if (file.exists(checkpoint)) {
z <- readRDS(checkpoint)
if (z$nsim != nsim || z$ninner != ninner)
stop("Checkpoint dimensions do not match the bootstrap file. Delete the checkpoint and restart.")
bootb_ngml <- z$bootb_ngml
bootb_ngml_inner <- z$bootb_ngml_inner
success <- z$success
success_inner <- z$success_inner
done <- z$done
done_inner <- z$done_inner
.Random.seed <- z$rng_state
cat("Resuming saved computation.\n")
} else {
bootb_ngml <- matrix(NA_real_, length(b_ngml0), nsim)
bootb_ngml_inner <- array(NA_real_, c(length(b_ngml0), ninner, nsim))
success <- rep(FALSE, nsim)
success_inner <- matrix(FALSE, ninner, nsim)
done <- rep(FALSE, nsim)
done_inner <- matrix(FALSE, ninner, nsim)
}
# Outer draws
for (j in seq_len(nsim)) {
if (done[j]) next
cat("Outer", j, "of", nsim, "\n")
Yboot <- h5read(matfile, "bootDatas", index = list(1:nobs_boot, 1:nvar_mat, j))
z <- try(estimate_b(Yboot, b_ngml0), silent = TRUE)
if (!inherits(z, "try-error")) {
bootb_ngml[, j] <- z
success[j] <- TRUE
}
done[j] <- TRUE
save_checkpoint()
}
rm(list = ls()); set.seed(42)
# WARNING: Full run (1,000 outer x 100 inner draws) may take 10 days or more.
pkgs <- c("vars", "svars", "zoo")
missing <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if (length(missing))
stop("Please install R package(s): ", paste(missing, collapse = ", "))
if (!requireNamespace("rhdf5", quietly = TRUE))
stop("Please install 'rhdf5' via Bioconductor: BiocManager::install('rhdf5', update = FALSE)")
suppressPackageStartupMessages({
library(vars); library(svars); library(zoo); library(rhdf5)
})
project_dir <- dirname(getwd())
in_dir  <- file.path(project_dir, "Data", "kaenzig_exports")
out_dir <- file.path(project_dir, "Data", "thispaper_exports")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
matfile <- file.path(in_dir, "kaenzig_nested_bootstrap_exports.mat")
checkpoint <- file.path(out_dir, "thispaper_wald_checkpoint.rds")
outfile <- file.path(out_dir, "thispaper_wald_test_bK_equals_bNG_nested.csv")
# Baseline nG estimate
D <- read.csv(file.path(in_dir, "kaenzig_baseline_data.csv"),
na.strings = c("NaN", "NA", ""))
ym <- as.yearmon(gsub("M", "-", D$time), "%Y-%m")
Y <- D[, 2:7]
Y[] <- lapply(Y, as.numeric)
ok <- complete.cases(Y)
Y <- Y[ok, , drop = FALSE]
ym <- ym[ok]
Y_ts <- ts(as.matrix(Y),
start = c(as.integer(format(ym[1], "%Y")),
as.integer(format(ym[1], "%m"))), frequency = 12)
news_full <- read.csv(file.path(in_dir, "kaenzig_baseline_structural_shock.csv"))$shock
plain.var0 <- VAR(Y_ts, p = 12, type = "const")
u0 <- resid(plain.var0)
T0 <- nrow(u0)
news <- tail(news_full, T0)
B0 <- id.ngml(plain.var0)$B
eps0 <- t(solve(B0, t(u0)))
oil_j0 <- which.max(abs(apply(eps0, 2, function(x) cor(x, news))))
b_ngml0 <- as.numeric(B0[, oil_j0])
if (!is.na(b_ngml0[1]) && b_ngml0[1] < 0) b_ngml0 <- -b_ngml0
# Bootstrap metadata
p_mat    <- as.integer(h5read(matfile, "p"))
nvar_mat <- as.integer(h5read(matfile, "nvar"))
nsim     <- as.integer(h5read(matfile, "nsim"))
ninner   <- as.integer(h5read(matfile, "ninner"))
nobs_boot <- T0 + p_mat
match_column <- function(Bstar, bref) {
sims <- rep(NA_real_, ncol(Bstar))
for (j in seq_len(ncol(Bstar))) {
bj <- as.numeric(Bstar[, j])
sims[j] <- abs(sum(bj * bref) / sqrt(sum(bj^2) * sum(bref^2)))
}
which.max(sims)
}
estimate_b <- function(Ymat, bref) {
Ymat <- drop(Ymat)
Ymat <- as.matrix(Ymat)
if (ncol(Ymat) != length(bref) && nrow(Ymat) == length(bref)) Ymat <- t(Ymat)
storage.mode(Ymat) <- "double"
Bstar <- id.ngml(VAR(ts(Ymat, frequency = 12), p = p_mat, type = "const"))$B
j <- match_column(Bstar, bref)
b <- as.numeric(Bstar[, j])
if (sum(b * bref) < 0) b <- -b
b
}
save_checkpoint <- function() saveRDS(list(
bootb_ngml = bootb_ngml, bootb_ngml_inner = bootb_ngml_inner,
success = success, success_inner = success_inner,
done = done, done_inner = done_inner,
nsim = nsim, ninner = ninner, rng_state = .Random.seed
), checkpoint)
# Initialize/resume
if (file.exists(checkpoint)) {
z <- readRDS(checkpoint)
if (z$nsim != nsim || z$ninner != ninner)
stop("Checkpoint dimensions do not match the bootstrap file. Delete the checkpoint and restart.")
bootb_ngml <- z$bootb_ngml
bootb_ngml_inner <- z$bootb_ngml_inner
success <- z$success
success_inner <- z$success_inner
done <- z$done
done_inner <- z$done_inner
.Random.seed <- z$rng_state
cat("Resuming saved computation.\n")
} else {
bootb_ngml <- matrix(NA_real_, length(b_ngml0), nsim)
bootb_ngml_inner <- array(NA_real_, c(length(b_ngml0), ninner, nsim))
success <- rep(FALSE, nsim)
success_inner <- matrix(FALSE, ninner, nsim)
done <- rep(FALSE, nsim)
done_inner <- matrix(FALSE, ninner, nsim)
}
# Outer draws
for (j in seq_len(nsim)) {
if (done[j]) next
cat("Outer", j, "of", nsim, "\n")
Yboot <- h5read(matfile, "bootDatas", index = list(1:nobs_boot, 1:nvar_mat, j))
z <- try(estimate_b(Yboot, b_ngml0), silent = TRUE)
if (!inherits(z, "try-error")) {
bootb_ngml[, j] <- z
success[j] <- TRUE
}
done_inner[jj, j] <- TRUE
if (jj %% 10 == 0 || jj == ninner) save_checkpoint()
}
if (jj %% 10 == 0 || jj == ninner) save_checkpoint()
rm(list = ls()); set.seed(42)
# WARNING: Full run (1,000 outer x 100 inner draws) may take 10 days or more.
pkgs <- c("vars", "svars", "zoo")
missing <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if (length(missing))
stop("Please install R package(s): ", paste(missing, collapse = ", "))
if (!requireNamespace("rhdf5", quietly = TRUE))
stop("Please install 'rhdf5' via Bioconductor: BiocManager::install('rhdf5', update = FALSE)")
suppressPackageStartupMessages({
library(vars); library(svars); library(zoo); library(rhdf5)
})
project_dir <- dirname(getwd())
in_dir  <- file.path(project_dir, "Data", "kaenzig_exports")
out_dir <- file.path(project_dir, "Data", "thispaper_exports")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
matfile <- file.path(in_dir, "kaenzig_nested_bootstrap_exports.mat")
checkpoint <- file.path(out_dir, "thispaper_wald_checkpoint.rds")
outfile <- file.path(out_dir, "thispaper_wald_test_bK_equals_bNG_nested.csv")
# Baseline nG estimate
D <- read.csv(file.path(in_dir, "kaenzig_baseline_data.csv"),
na.strings = c("NaN", "NA", ""))
ym <- as.yearmon(gsub("M", "-", D$time), "%Y-%m")
Y <- D[, 2:7]
Y[] <- lapply(Y, as.numeric)
ok <- complete.cases(Y)
Y <- Y[ok, , drop = FALSE]
ym <- ym[ok]
Y_ts <- ts(as.matrix(Y),
start = c(as.integer(format(ym[1], "%Y")),
as.integer(format(ym[1], "%m"))), frequency = 12)
news_full <- read.csv(file.path(in_dir, "kaenzig_baseline_structural_shock.csv"))$shock
plain.var0 <- VAR(Y_ts, p = 12, type = "const")
u0 <- resid(plain.var0)
T0 <- nrow(u0)
news <- tail(news_full, T0)
B0 <- id.ngml(plain.var0)$B
eps0 <- t(solve(B0, t(u0)))
oil_j0 <- which.max(abs(apply(eps0, 2, function(x) cor(x, news))))
b_ngml0 <- as.numeric(B0[, oil_j0])
if (!is.na(b_ngml0[1]) && b_ngml0[1] < 0) b_ngml0 <- -b_ngml0
# Bootstrap metadata
p_mat    <- as.integer(h5read(matfile, "p"))
nvar_mat <- as.integer(h5read(matfile, "nvar"))
nsim     <- as.integer(h5read(matfile, "nsim"))
ninner   <- as.integer(h5read(matfile, "ninner"))
nobs_boot <- T0 + p_mat
match_column <- function(Bstar, bref) {
sims <- rep(NA_real_, ncol(Bstar))
for (j in seq_len(ncol(Bstar))) {
bj <- as.numeric(Bstar[, j])
sims[j] <- abs(sum(bj * bref) / sqrt(sum(bj^2) * sum(bref^2)))
}
which.max(sims)
}
estimate_b <- function(Ymat, bref) {
Ymat <- drop(Ymat)
Ymat <- as.matrix(Ymat)
if (ncol(Ymat) != length(bref) && nrow(Ymat) == length(bref)) Ymat <- t(Ymat)
storage.mode(Ymat) <- "double"
Bstar <- id.ngml(VAR(ts(Ymat, frequency = 12), p = p_mat, type = "const"))$B
j <- match_column(Bstar, bref)
b <- as.numeric(Bstar[, j])
if (sum(b * bref) < 0) b <- -b
b
}
save_checkpoint <- function() saveRDS(list(
bootb_ngml = bootb_ngml, bootb_ngml_inner = bootb_ngml_inner,
success = success, success_inner = success_inner,
done = done, done_inner = done_inner,
nsim = nsim, ninner = ninner, rng_state = .Random.seed
), checkpoint)
# Initialize/resume
if (file.exists(checkpoint)) {
z <- readRDS(checkpoint)
if (z$nsim != nsim || z$ninner != ninner)
stop("Checkpoint dimensions do not match the bootstrap file. Delete the checkpoint and restart.")
bootb_ngml <- z$bootb_ngml
bootb_ngml_inner <- z$bootb_ngml_inner
success <- z$success
success_inner <- z$success_inner
done <- z$done
done_inner <- z$done_inner
.Random.seed <- z$rng_state
cat("Resuming saved computation.\n")
} else {
bootb_ngml <- matrix(NA_real_, length(b_ngml0), nsim)
bootb_ngml_inner <- array(NA_real_, c(length(b_ngml0), ninner, nsim))
success <- rep(FALSE, nsim)
success_inner <- matrix(FALSE, ninner, nsim)
done <- rep(FALSE, nsim)
done_inner <- matrix(FALSE, ninner, nsim)
}
# Outer draws
for (j in seq_len(nsim)) {
if (done[j]) next
cat("Outer", j, "of", nsim, "\n")
Yboot <- h5read(matfile, "bootDatas", index = list(1:nobs_boot, 1:nvar_mat, j))
z <- try(estimate_b(Yboot, b_ngml0), silent = TRUE)
if (!inherits(z, "try-error")) {
bootb_ngml[, j] <- z
success[j] <- TRUE
}
done[j] <- TRUE
save_checkpoint()
}
