Corrected QML, Monte Carlo inference, and empirical estimation
This document provides a reproducible implementation of the multivariate spatial and spatiotemporal ARCH (vec-spARCH) model introduced in Otto (2024), together with the corrected Gaussian quasi-log-likelihood and the additional inferential analysis reported in the accompanying corrigendum.
The main correction concerns the normalisation of the Gaussian quasi-log-likelihood. In the original article, both parameter-dependent terms, namely the log-determinant term and the quadratic residual term, were divided by the combined dimension \(np\). This does not change the maximiser of the criterion, because the parameter-dependent part is multiplied by a common positive constant. It does, however, rescale the curvature of the criterion and therefore affects inverse-Hessian covariance estimates and the corresponding standard errors.
The corrected implementation below separates three questions:
The document additionally reproduces the Monte Carlo design of the original article, computes Monte Carlo standard deviations, standard-error calibration ratios, and 95% coverage probabilities, and provides code for the empirical estimator with variable-specific intercepts.
The corrected likelihood normalisation leaves the QML point estimates essentially unchanged, but substantially reduces the Hessian-based standard errors. In the Berlin real-estate application this strengthens the inferential conclusions: temporal dependence remains stronger than spatial dependence, but the corrected inference provides evidence for instantaneous spatial cross-variable interactions in addition to temporally lagged cross-variable interactions.
The Monte Carlo study shows that the corrected Hessian-based standard errors are generally of the right order of magnitude relative to the empirical Monte Carlo standard deviation. The largest deviations occur for some zero temporal cross-effects and for the intercept and selected spatial coefficients in the purely spatial Model B. Sandwich-based standard errors behave very similarly to the inverse-Hessian standard errors.
A special feature of Model B is that \(\Pi_0=0\). With a spatially constant intercept and a row-standardised spatial weights matrix, the unconditional mean belongs to the constant spatial subspace. The mean alone therefore does not separately identify the intercept and \(\Psi\); their separation additionally relies on the spatial covariance structure and the Jacobian term of the QML criterion. This helps explain the stronger finite-sample interaction between the corresponding estimates in that design.
Let
\[ \ddot{\mathbf Y}_t = \operatorname{vec}\!\left(\mathbf Y_t^{(\ln,2)}\right), \qquad \mathbf S_{np} = \mathbf I_{np}-\mathbf\Psi^\prime\otimes\mathbf W, \]
and define
\[ \mathbf r_t(\boldsymbol\vartheta) = \mathbf S_{np}\ddot{\mathbf Y}_t - \operatorname{vec}(\widetilde{\mathbf A}) - (\mathbf\Pi^\prime\otimes\mathbf I_n) \ddot{\mathbf Y}_{t-1}. \]
For \(T\) conditional transitions, the correctly scaled Gaussian quasi-log-likelihood is
\[ \ell_{nT}^{\mathrm{cor}}(\boldsymbol\vartheta\mid\mathbf Y_0) = -\frac{Tnp}{2}\log(2\pi) -\frac{Tnp}{2}\log(\sigma_u^2) +T\log|\mathbf S_{np}| -\frac{1}{2\sigma_u^2}\sum_{t=1}^T \mathbf r_t(\boldsymbol\vartheta)^\prime\mathbf r_t(\boldsymbol\vartheta). \]
In the implementation below, the first observed time point is conditioned upon, so an array containing TT observations contributes TT - 1 conditional likelihood terms.
If the originally normalised criterion is denoted by \(\ell^{\mathrm{old}}_{nT}\), then
\[ \ell^{\mathrm{old}}_{nT}(\boldsymbol\vartheta) = \frac{1}{np}\ell^{\mathrm{cor}}_{nT}(\boldsymbol\vartheta)+C, \]
where \(C\) is independent of \(\boldsymbol\vartheta\). Consequently,
\[ \arg\max_{\boldsymbol\vartheta}\ell^{\mathrm{old}}_{nT} = \arg\max_{\boldsymbol\vartheta}\ell^{\mathrm{cor}}_{nT}. \]
For the Hessians,
\[ \nabla^2\ell^{\mathrm{old}}_{nT} = \frac{1}{np}\nabla^2\ell^{\mathrm{cor}}_{nT}, \]
and hence
\[ \widehat{\operatorname{Var}}_{\mathrm{cor}}(\widehat{\boldsymbol\vartheta}) = \frac{1}{np}\widehat{\operatorname{Var}}_{\mathrm{old}}(\widehat{\boldsymbol\vartheta}), \qquad \operatorname{se}_{\mathrm{cor}} = \frac{1}{\sqrt{np}}\operatorname{se}_{\mathrm{old}}. \]
For standard normal innovations,
\[ E(\log\varepsilon^2)=\psi(1/2)+\log 2=-\gamma-\log 2, \]
and
\[ \sigma_u^2=\operatorname{Var}(\log\varepsilon^2)=\psi_1(1/2)=\frac{\pi^2}{2}. \]
The historical simulation code used raw \(t_3\) innovations. Such innovations have variance 3, whereas the model formulation normalises the innovation variance to one. This document therefore exposes the choice explicitly through the parameter standardise_t.
standardise_t: false reproduces the historical simulation convention, with \(E(\log\varepsilon^2)=\log 3-2\).standardise_t: true uses \(\varepsilon=T_3/\sqrt{3}\), so that \(\operatorname{Var}(\varepsilon)=1\) and \(E(\log\varepsilon^2)=-2\).The variance of \(\log\varepsilon^2\) is unchanged by this rescaling:
\[ \operatorname{Var}(\log\varepsilon^2)=\pi^2-4. \]
For exact reproduction of the historical Monte Carlo output, keep the default standardise_t: false. For a simulation aligned exactly with the unit-variance innovation convention of the model, set it to true.
required_packages <- c(
"spdep", "Rsolnp", "numDeriv",
"ggplot2", "dplyr", "tidyr", "parallel"
)
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
stop(
"Please install the following packages before rendering: ",
paste(missing_packages, collapse = ", ")
)
}
library(spdep)
library(Rsolnp)
library(numDeriv)
library(ggplot2)
library(dplyr)
library(tidyr)
library(parallel)logsq <- function(x) log(x^2)
logsq_back <- function(x, sgn) sgn * sqrt(exp(x))
vec <- function(x) as.vector(x)
innovation_spec <- function(
errortype = c("norm", "t"),
standardise_t = params$standardise_t
) {
errortype <- match.arg(errortype)
if (errortype == "norm") {
return(list(
generator = function(n, p) array(rnorm(n * p), dim = c(n, p)),
E_logsq = digamma(0.5) + log(2),
sig_u2 = trigamma(0.5)
))
}
if (isTRUE(standardise_t)) {
generator <- function(n, p) {
array(rt(n * p, df = 3) / sqrt(3), dim = c(n, p))
}
E_logsq <- -2
} else {
generator <- function(n, p) {
array(rt(n * p, df = 3), dim = c(n, p))
}
E_logsq <- log(3) - 2
}
list(
generator = generator,
E_logsq = E_logsq,
sig_u2 = pi^2 - 4
)
}The default burn-in is 20 iterations to reproduce the historical simulation design. For new applications, a longer burn-in can be supplied through burn_in.
simulate_vec_spARCH <- function(
n, t, p, W, parameters,
burn_in = 20,
standardise_t = params$standardise_t
) {
if (!all(dim(W) == c(n, n))) stop("Dimension of W is incompatible with n.")
if (is.null(parameters$a)) stop("parameters$a must be supplied.")
if (!all(dim(parameters$Psi) == c(p, p))) stop("parameters$Psi must be p x p.")
if (!all(dim(parameters$Pi) == c(p, p))) stop("parameters$Pi must be p x p.")
spec <- innovation_spec(parameters$errortype, standardise_t)
err_distr <- spec$generator
E_logsq_errors <- spec$E_logsq
Psi <- parameters$Psi
Pi <- parameters$Pi
A <- matrix(parameters$a, nrow = n, ncol = p)
A_tilde <- A + E_logsq_errors
S <- diag(n * p) - kronecker(t(Psi), W)
S_inv <- solve(S)
eps <- err_distr(n, p)
Y_lagged <- eps
simulate_one <- function(Y_lagged) {
eps <- err_distr(n, p)
u <- logsq(eps) - E_logsq_errors
vec_log_y2 <- S_inv %*% (
vec(A_tilde) +
kronecker(t(Pi), diag(n)) %*% vec(logsq(Y_lagged)) +
vec(u)
)
y <- logsq_back(vec_log_y2, sign(vec(eps)))
Y_new <- matrix(y, nrow = n, ncol = p)
log_h <- matrix(vec_log_y2, nrow = n, ncol = p) - logsq(eps)
H_new <- exp(log_h)
list(Y = Y_new, H = H_new)
}
for (b in seq_len(burn_in)) {
step <- simulate_one(Y_lagged)
Y_lagged <- step$Y
}
Y <- array(NA_real_, dim = c(n, p, t))
H <- array(NA_real_, dim = c(n, p, t))
for (tt in seq_len(t)) {
step <- simulate_one(Y_lagged)
Y[, , tt] <- step$Y
H[, , tt] <- step$H
Y_lagged <- step$Y
}
list(
Y = Y,
H = H,
parameters = parameters,
E_logsq_errors = E_logsq_errors,
sig_u2 = spec$sig_u2
)
}The score function supports both a single common intercept, as used in the Monte Carlo study, and \(p\) variable-specific intercepts, as used in the empirical application. The returned matrix has one row for each conditional transition and is used to construct the outer-product matrix \(J\) in the sandwich covariance estimator.
qml_score_time <- function(
par, Y, W, sig_u2,
intercept = c("scalar", "vector")
) {
intercept <- match.arg(intercept)
n <- dim(Y)[1]
p <- dim(Y)[2]
TT <- dim(Y)[3]
np <- n * p
q <- if (intercept == "scalar") 1L else p
idx_a <- seq_len(q)
idx_psi <- q + seq_len(p^2)
idx_pi <- q + p^2 + seq_len(p^2)
a_tilde <- par[idx_a]
Psi <- matrix(par[idx_psi], p, p)
Pi <- matrix(par[idx_pi], p, p)
if (intercept == "scalar") {
A_tilde <- matrix(a_tilde, nrow = n, ncol = p)
X_a <- matrix(1, nrow = np, ncol = 1)
} else {
A_tilde <- matrix(rep(a_tilde, each = n), nrow = n, ncol = p)
X_a <- kronecker(diag(p), matrix(1, nrow = n, ncol = 1))
}
S <- diag(np) - kronecker(t(Psi), W)
S_inv <- solve(S)
K_list <- vector("list", p^2)
tr_term <- numeric(p^2)
counter <- 1L
# Column-major order: psi_11, psi_21, ..., psi_p1, psi_12, ...
for (j in seq_len(p)) {
for (i in seq_len(p)) {
Eji <- matrix(0, p, p)
Eji[j, i] <- 1
K <- kronecker(Eji, W)
K_list[[counter]] <- K
tr_term[counter] <- sum(S_inv * t(K))
counter <- counter + 1L
}
}
scores <- matrix(NA_real_, nrow = TT - 1, ncol = q + 2 * p^2)
for (tt in 2:TT) {
z_t <- vec(logsq(Y[, , tt]))
Y_lag <- logsq(Y[, , tt - 1])
X_pi <- kronecker(diag(p), Y_lag)
r_t <- S %*% z_t - vec(A_tilde) - X_pi %*% vec(Pi)
score <- numeric(q + 2 * p^2)
score[idx_a] <- -as.vector(crossprod(X_a, r_t)) / sig_u2
for (k in seq_len(p^2)) {
Kz <- K_list[[k]] %*% z_t
score[idx_psi[k]] <-
tr_term[k] - drop(crossprod(r_t, Kz)) / sig_u2
}
score[idx_pi] <- -as.vector(crossprod(X_pi, r_t)) / sig_u2
scores[tt - 1, ] <- score
}
scores
}The same estimator is used for the Monte Carlo study and the empirical application. Set intercept = "scalar" for the simulation design and intercept = "vector" for variable-specific empirical intercepts. The optimiser works with \(\widetilde A\); the returned theta transforms the intercept back to \(A\) using \(E(\log\varepsilon^2)\).
qml_vec_spARCH <- function(
Y, W, parameters,
intercept = c("scalar", "vector"),
start = NULL,
standardise_t = params$standardise_t,
trace = FALSE
) {
intercept <- match.arg(intercept)
dimY <- dim(Y)
n <- dimY[1]
p <- dimY[2]
TT <- dimY[3]
np <- n * p
if (!all(dim(W) == c(n, n))) stop("Dimension of W is incompatible with Y.")
spec <- innovation_spec(parameters$errortype, standardise_t)
E_logsq_errors <- spec$E_logsq
sig_u2 <- spec$sig_u2
q <- if (intercept == "scalar") 1L else p
unpack_parameters <- function(pars) {
idx_a <- seq_len(q)
idx_psi <- q + seq_len(p^2)
idx_pi <- q + p^2 + seq_len(p^2)
a_tilde <- pars[idx_a]
if (intercept == "scalar") {
A_tilde <- matrix(a_tilde, nrow = n, ncol = p)
} else {
A_tilde <- matrix(rep(a_tilde, each = n), nrow = n, ncol = p)
}
list(
idx_a = idx_a,
idx_psi = idx_psi,
idx_pi = idx_pi,
a_tilde = a_tilde,
A_tilde = A_tilde,
Psi = matrix(pars[idx_psi], p, p),
Pi = matrix(pars[idx_pi], p, p)
)
}
negative_qml <- function(pars, Y, W, sig_u2) {
par_list <- unpack_parameters(pars)
S <- diag(np) - kronecker(t(par_list$Psi), W)
detS <- determinant(S, logarithm = TRUE)
log_det_S <- as.numeric(detS$modulus)
if (!is.finite(log_det_S)) return(1e12)
sum_res2 <- 0
for (tt in 2:TT) {
r_t <- S %*% vec(logsq(Y[, , tt])) -
vec(par_list$A_tilde) -
kronecker(diag(p), logsq(Y[, , tt - 1])) %*% vec(par_list$Pi)
if (any(!is.finite(r_t))) return(1e12)
sum_res2 <- sum_res2 + sum(r_t^2)
}
n_transitions <- TT - 1
ll <-
-n_transitions * np / 2 * log(2 * pi) -
n_transitions * np / 2 * log(sig_u2) +
n_transitions * log_det_S -
sum_res2 / (2 * sig_u2)
-ll
}
residuals_fun <- function(pars) {
par_list <- unpack_parameters(pars)
S <- diag(np) - kronecker(t(par_list$Psi), W)
U <- matrix(NA_real_, nrow = np, ncol = TT - 1)
for (tt in 2:TT) {
U[, tt - 1] <- S %*% vec(logsq(Y[, , tt])) -
vec(par_list$A_tilde) -
kronecker(diag(p), logsq(Y[, , tt - 1])) %*% vec(par_list$Pi)
}
U
}
if (is.null(start)) {
start_a <- if (intercept == "scalar") 1 else rep(1, p)
if (p == 2) {
start_Psi <- matrix(c(0.5, 0.1, 0.1, 0.5), p, p)
start_Pi <- matrix(c(0.3, 0, 0, 0.3), p, p)
} else {
start_Psi <- diag(0.2, p)
start_Pi <- diag(0.2, p)
}
start <- c(start_a, vec(start_Psi), vec(start_Pi))
}
out <- solnp(
pars = start,
fun = negative_qml,
Y = Y,
W = W,
sig_u2 = sig_u2,
control = list(trace = trace)
)
par_opt <- out$pars
par_list <- unpack_parameters(par_opt)
H <- numDeriv::hessian(
func = negative_qml,
x = par_opt,
Y = Y,
W = W,
sig_u2 = sig_u2
)
H <- (H + t(H)) / 2
V_hessian <- tryCatch(
solve(H),
error = function(e) matrix(NA_real_, nrow(H), ncol(H))
)
stderr_hessian <- rep(NA_real_, nrow(H))
if (all(is.finite(V_hessian))) {
diag_v <- diag(V_hessian)
keep <- diag_v > 0
stderr_hessian[keep] <- sqrt(diag_v[keep])
}
score_t <- qml_score_time(
par = par_opt,
Y = Y,
W = W,
sig_u2 = sig_u2,
intercept = intercept
)
J <- crossprod(score_t)
V_sandwich <- V_hessian %*% J %*% V_hessian
V_sandwich <- (V_sandwich + t(V_sandwich)) / 2
stderr_sandwich <- rep(NA_real_, nrow(V_sandwich))
if (all(is.finite(V_sandwich))) {
diag_v <- diag(V_sandwich)
keep <- diag_v > 0
stderr_sandwich[keep] <- sqrt(diag_v[keep])
}
a_est <- par_list$a_tilde - E_logsq_errors
if (intercept == "scalar") {
A_est <- matrix(a_est, nrow = n, ncol = p)
} else {
A_est <- matrix(rep(a_est, each = n), nrow = n, ncol = p)
}
theta <- c(a_est, vec(par_list$Psi), vec(par_list$Pi))
theta_tilde <- c(par_list$a_tilde, vec(par_list$Psi), vec(par_list$Pi))
res <- residuals_fun(par_opt)
list(
theta = theta,
theta_tilde = theta_tilde,
par_opt = par_opt,
a_est = a_est,
a_tilde_est = par_list$a_tilde,
A_est = A_est,
A_tilde_est = par_list$A_tilde,
Psi_est = par_list$Psi,
Pi_est = par_list$Pi,
H = H,
vcov_hessian = V_hessian,
stderr_hessian = stderr_hessian,
score_t = score_t,
J = J,
vcov_sandwich = V_sandwich,
stderr_sandwich = stderr_sandwich,
ll = -tail(out$values, 1),
residuals = res,
var_res = mean(res^2, na.rm = TRUE),
E_logsq_errors = E_logsq_errors,
sig_u2 = sig_u2,
intercept = intercept
)
}set.seed(1)
p <- 2
d <- 5
n <- d^2
TT <- 50
W <- nb2mat(
cell2nb(d, d, type = "queen"),
style = "W"
)
parameters <- list(
a = 1,
Psi = matrix(
c(0.5, 0.1,
0.1, 0.5),
p, p
),
Pi = matrix(
c(0.3, 0,
0, 0.3),
p, p
),
errortype = "norm"
)
# Simulate data
sim <- simulate_vec_spARCH(
n = n,
t = TT,
p = p,
W = W,
parameters = parameters
)
# Estimate model
est <- qml_vec_spARCH(
Y = sim$Y,
W = W,
parameters = parameters,
intercept = "scalar"
)
# ------------------------------------------------------------
# Compare true and estimated parameters
# ------------------------------------------------------------
true_parameters <- c(
parameters$a,
vec(parameters$Psi),
vec(parameters$Pi)
)
parameter_names <- c(
"a",
"psi_11", "psi_21", "psi_12", "psi_22",
"pi_11", "pi_21", "pi_12", "pi_22"
)
results <- data.frame(
Parameter = parameter_names,
True = true_parameters,
Estimate = est$theta,
SE_Hessian = est$stderr_hessian,
SE_Sandwich = est$stderr_sandwich
)
results[, -1] <- round(results[, -1], 3)
results Parameter True Estimate SE_Hessian SE_Sandwich
1 a 1.0 1.002 0.073 0.057
2 psi_11 0.5 0.454 0.033 0.023
3 psi_21 0.1 0.073 0.041 0.038
4 psi_12 0.1 0.118 0.046 0.050
5 psi_22 0.5 0.472 0.030 0.032
6 pi_11 0.3 0.310 0.024 0.021
7 pi_21 0.0 0.042 0.024 0.022
8 pi_12 0.0 -0.022 0.026 0.024
9 pi_22 0.3 0.314 0.023 0.026The original simulation study considers three bivariate models:
\[ \text{Model A:}\qquad \Psi_0= \begin{pmatrix}0.5&0.1\\0.1&0.5\end{pmatrix}, \qquad \Pi_0= \begin{pmatrix}0.3&0\\0&0.3\end{pmatrix}, \]
\[ \text{Model B:}\qquad \Psi_0= \begin{pmatrix}0.5&0.1\\0.1&0.5\end{pmatrix}, \qquad \Pi_0=\mathbf 0, \]
and
\[ \text{Model C:}\qquad \Psi_0= \begin{pmatrix}0.2&0.4\\0.4&0.2\end{pmatrix}, \qquad \Pi_0= \begin{pmatrix}0.3&0\\0&0.3\end{pmatrix}. \]
For all three models, \(a_0=1\). The sample-size combinations are
\[ (n,T)\in\{(25,30),(49,100),(100,200)\}, \]
with row-standardised Queen contiguity weights and both normal and \(t_3\) innovations.
n_set <- c(25, 49, 100)
t_set <- c(30, 100, 200)
p_set <- 2
make_design <- function(model, n, error) {
p <- p_set
parameters <- list(a = 1)
if (model == "M1") {
parameters$Psi <- matrix(c(0.5, 0.1, 0.1, 0.5), p, p)
parameters$Pi <- matrix(c(0.3, 0, 0, 0.3), p, p)
}
if (model == "M2") {
parameters$Psi <- matrix(c(0.5, 0.1, 0.1, 0.5), p, p)
parameters$Pi <- matrix(0, p, p)
}
if (model == "M3") {
parameters$Psi <- matrix(c(0.2, 0.4, 0.4, 0.2), p, p)
parameters$Pi <- matrix(c(0.3, 0, 0, 0.3), p, p)
}
parameters$errortype <- error
W <- nb2mat(
cell2nb(sqrt(n), sqrt(n), type = "queen"),
style = "W"
)
list(
parameters = parameters,
W = W,
dgp = c(parameters$a, vec(parameters$Psi), vec(parameters$Pi))
)
}The full simulation is not run automatically when rendering the web page. Run the following chunk interactively, or change its chunk option to eval: true, to generate the 1000-replication result file.
run_replication <- function(
i, n, TT, p, W, parameters, dgp,
seed_base, standardise_t
) {
set.seed(seed_base + i)
tryCatch({
sim <- simulate_vec_spARCH(
n = n,
t = TT,
p = p,
W = W,
parameters = parameters,
standardise_t = standardise_t
)
est <- qml_vec_spARCH(
Y = sim$Y,
W = W,
parameters = parameters,
intercept = "scalar",
standardise_t = standardise_t
)
list(
success = TRUE,
replication = i,
output = est,
n = n,
t = TT,
p = p,
dgp = dgp,
error = parameters$errortype
)
}, error = function(e) {
list(
success = FALSE,
replication = i,
message = conditionMessage(e),
n = n,
t = TT,
p = p,
dgp = dgp,
error = parameters$errortype
)
})
}if (isTRUE(params$run_full_mc)) {
m <- params$mc_replications
available_cores <- max(1, parallel::detectCores() - 2)
n_cores <- min(params$cores, available_cores)
cat("Using", n_cores, "cores\n")
cl <- makeCluster(n_cores)
clusterEvalQ(cl, {
library(spdep)
library(Rsolnp)
library(numDeriv)
NULL
})
clusterExport(
cl,
varlist = c(
"logsq", "logsq_back", "vec",
"innovation_spec", "simulate_vec_spARCH",
"qml_score_time", "qml_vec_spARCH",
"run_replication"
),
envir = environment()
)
M1 <- M2 <- M3 <- vector("list", length = 6)
counter <- 1L
for (size in seq_along(n_set)) {
n <- n_set[size]
TT <- t_set[size]
p <- p_set
for (error_distrib in c("norm", "t")) {
for (model in c("M1", "M2", "M3")) {
design <- make_design(model = model, n = n, error = error_distrib)
cat(
"\nRunning", model,
"| n =", n,
"| T =", TT,
"| error =", error_distrib,
"| replications =", m, "\n"
)
result_list <- parLapplyLB(
cl,
X = seq_len(m),
fun = run_replication,
n = n,
TT = TT,
p = p,
W = design$W,
parameters = design$parameters,
dgp = design$dgp,
seed_base = 100000 * counter,
standardise_t = params$standardise_t
)
if (model == "M1") M1[[counter]] <- result_list
if (model == "M2") M2[[counter]] <- result_list
if (model == "M3") M3[[counter]] <- result_list
save(M1, M2, M3, file = "Output_vecspARCH_MC_temp.rda")
}
counter <- counter + 1L
}
}
stopCluster(cl)
save(M1, M2, M3, file = "Output_vecspARCH_MC_1000.rda")
} else {
message("Full Monte Carlo simulation skipped. Set run_full_mc: true to run it.")
}For a website render, run the full simulation once, save Output_vecspARCH_MC_1000.rda in the same directory, and let the page load the stored results.
result_file <- "Output_vecspARCH_MC_1000.rda"
if (file.exists(result_file)) load(result_file)
have_mc_results <- all(
c("M1", "M2", "M3") %in% ls(envir = environment())
)par_names <- c(
"a",
"psi_11", "psi_21", "psi_12", "psi_22",
"pi_11", "pi_21", "pi_12", "pi_22"
)
successful_results <- function(results) {
Filter(function(x) isTRUE(x$success), results)
}
summarise_MC <- function(results, model, error, n, TT, par_names) {
results <- successful_results(results)
if (length(results) == 0) stop("No successful Monte Carlo replications.")
theta <- do.call(rbind, lapply(results, function(x) x$output$theta))
se_hessian <- do.call(rbind, lapply(results, function(x) x$output$stderr_hessian))
se_sandwich <- do.call(rbind, lapply(results, function(x) x$output$stderr_sandwich))
dgp <- results[[1]]$dgp
p <- results[[1]]$p
error_mat <- sweep(theta, 2, dgp, "-")
z <- qnorm(0.975)
coverage_hessian <- colMeans(
abs(error_mat) <= z * se_hessian,
na.rm = TRUE
)
coverage_sandwich <- colMeans(
abs(error_mat) <= z * se_sandwich,
na.rm = TRUE
)
se_old <- sqrt(n * p) * se_hessian
coverage_old <- colMeans(
abs(error_mat) <= z * se_old,
na.rm = TRUE
)
MCSD <- apply(theta, 2, sd, na.rm = TRUE)
mean_se_hessian <- colMeans(se_hessian, na.rm = TRUE)
mean_se_sandwich <- colMeans(se_sandwich, na.rm = TRUE)
data.frame(
model = model,
error = error,
n = n,
T = TT,
parameter = par_names,
true = dgp,
bias = colMeans(error_mat, na.rm = TRUE),
rmse = sqrt(colMeans(error_mat^2, na.rm = TRUE)),
MCSD = MCSD,
mean_SE_H = mean_se_hessian,
SE_ratio_H = mean_se_hessian / MCSD,
coverage_H = coverage_hessian,
mean_SE_S = mean_se_sandwich,
SE_ratio_S = mean_se_sandwich / MCSD,
coverage_S = coverage_sandwich,
coverage_old = coverage_old,
valid_H = colSums(is.finite(se_hessian)),
valid_S = colSums(is.finite(se_sandwich))
)
}if (have_mc_results) {
MC_summary <- data.frame()
for (model in c("M1", "M2", "M3")) {
object <- get(model)
for (err in 1:2) {
for (size in seq_along(n_set)) {
ind_model <- 2 * (size - 1) + err
error_name <- if (err == 1) "Normal" else "t3"
tmp <- summarise_MC(
results = object[[ind_model]],
model = model,
error = error_name,
n = n_set[size],
TT = t_set[size],
par_names = par_names
)
MC_summary <- bind_rows(MC_summary, tmp)
}
}
}
MC_summary <- MC_summary %>%
mutate(
model = factor(
model,
levels = c("M1", "M2", "M3"),
labels = c("Model A", "Model B", "Model C")
),
T_fac = factor(T),
bias_ratio_H = bias / mean_SE_H,
bias_ratio_S = bias / mean_SE_S
)
}if (have_mc_results) {
MC_summary %>%
select(
model, error, n, T, parameter,
bias, rmse, MCSD,
mean_SE_H, SE_ratio_H, coverage_H,
mean_SE_S, SE_ratio_S, coverage_S,
coverage_old
) %>%
mutate(across(where(is.numeric), ~ round(.x, 4)))
} model error n T parameter bias rmse MCSD mean_SE_H SE_ratio_H
1 Model A Normal 25 30 a -0.0642 0.1508 0.1371 0.1165 0.8497
2 Model A Normal 25 30 psi_11 -0.0081 0.0420 0.0414 0.0417 1.0066
3 Model A Normal 25 30 psi_21 -0.0085 0.0579 0.0575 0.0579 1.0074
4 Model A Normal 25 30 psi_12 -0.0015 0.0626 0.0629 0.0585 0.9300
5 Model A Normal 25 30 psi_22 -0.0114 0.0439 0.0426 0.0414 0.9717
6 Model A Normal 25 30 pi_11 -0.0070 0.0310 0.0304 0.0317 1.0445
7 Model A Normal 25 30 pi_21 0.0023 0.0290 0.0290 0.0333 1.1483
8 Model A Normal 25 30 pi_12 -0.0016 0.0349 0.0350 0.0338 0.9647
9 Model A Normal 25 30 pi_22 -0.0062 0.0349 0.0345 0.0315 0.9131
10 Model A Normal 49 100 a -0.0165 0.0465 0.0437 0.0448 1.0261
11 Model A Normal 49 100 psi_11 -0.0008 0.0178 0.0179 0.0166 0.9294
12 Model A Normal 49 100 psi_21 -0.0012 0.0230 0.0231 0.0227 0.9829
13 Model A Normal 49 100 psi_12 -0.0007 0.0229 0.0230 0.0226 0.9856
14 Model A Normal 49 100 psi_22 -0.0034 0.0156 0.0152 0.0166 1.0902
15 Model A Normal 49 100 pi_11 -0.0034 0.0130 0.0126 0.0122 0.9653
16 Model A Normal 49 100 pi_21 0.0007 0.0124 0.0124 0.0130 1.0430
17 Model A Normal 49 100 pi_12 -0.0001 0.0095 0.0096 0.0129 1.3520
18 Model A Normal 49 100 pi_22 -0.0016 0.0122 0.0121 0.0122 1.0101
19 Model A Normal 100 200 a 0.0008 0.0268 0.0269 0.0227 0.8428
20 Model A Normal 100 200 psi_11 -0.0008 0.0107 0.0107 0.0086 0.8079
21 Model A Normal 100 200 psi_21 0.0011 0.0129 0.0129 0.0116 0.9034
22 Model A Normal 100 200 psi_12 0.0029 0.0118 0.0115 0.0116 1.0091
23 Model A Normal 100 200 psi_22 -0.0011 0.0098 0.0098 0.0086 0.8842
24 Model A Normal 100 200 pi_11 0.0001 0.0067 0.0067 0.0061 0.9104
25 Model A Normal 100 200 pi_21 0.0008 0.0052 0.0052 0.0064 1.2485
26 Model A Normal 100 200 pi_12 -0.0001 0.0052 0.0053 0.0064 1.2250
27 Model A Normal 100 200 pi_22 -0.0001 0.0078 0.0079 0.0061 0.7758
28 Model A t3 25 30 a 0.0308 0.1063 0.1022 0.0788 0.7714
29 Model A t3 25 30 psi_11 -0.0098 0.0427 0.0418 0.0419 1.0014
30 Model A t3 25 30 psi_21 -0.0095 0.0737 0.0735 0.0683 0.9302
31 Model A t3 25 30 psi_12 0.0041 0.0771 0.0774 0.0688 0.8893
32 Model A t3 25 30 psi_22 -0.0140 0.0486 0.0468 0.0416 0.8886
33 Model A t3 25 30 pi_11 -0.0048 0.0310 0.0308 0.0317 1.0291
34 Model A t3 25 30 pi_21 -0.0027 0.0291 0.0291 0.0336 1.1517
35 Model A t3 25 30 pi_12 -0.0046 0.0283 0.0281 0.0337 1.2003
36 Model A t3 25 30 pi_22 -0.0040 0.0348 0.0348 0.0316 0.9088
37 Model A t3 49 100 a 0.0016 0.0322 0.0323 0.0285 0.8824
38 Model A t3 49 100 psi_11 -0.0025 0.0179 0.0178 0.0167 0.9366
39 Model A t3 49 100 psi_21 0.0014 0.0262 0.0263 0.0273 1.0378
40 Model A t3 49 100 psi_12 -0.0012 0.0254 0.0255 0.0275 1.0776
41 Model A t3 49 100 psi_22 -0.0012 0.0161 0.0162 0.0167 1.0317
42 Model A t3 49 100 pi_11 -0.0028 0.0136 0.0134 0.0122 0.9167
43 Model A t3 49 100 pi_21 0.0000 0.0082 0.0082 0.0131 1.5903
44 Model A t3 49 100 pi_12 0.0001 0.0071 0.0072 0.0130 1.8165
45 Model A t3 49 100 pi_22 -0.0007 0.0125 0.0125 0.0122 0.9798
46 Model A t3 100 200 a -0.0012 0.0251 0.0252 0.0142 0.5627
47 Model A t3 100 200 psi_11 0.0001 0.0090 0.0091 0.0086 0.9518
48 Model A t3 100 200 psi_21 -0.0011 0.0135 0.0136 0.0141 1.0388
49 Model A t3 100 200 psi_12 0.0018 0.0129 0.0128 0.0140 1.0927
50 Model A t3 100 200 psi_22 -0.0017 0.0086 0.0085 0.0086 1.0152
51 Model A t3 100 200 pi_11 0.0008 0.0085 0.0085 0.0061 0.7192
52 Model A t3 100 200 pi_21 0.0000 0.0028 0.0028 0.0065 2.2864
53 Model A t3 100 200 pi_12 -0.0001 0.0030 0.0030 0.0065 2.1280
54 Model A t3 100 200 pi_22 -0.0008 0.0086 0.0086 0.0061 0.7112
55 Model B Normal 25 30 a 0.0165 0.1165 0.1159 0.0674 0.5812
56 Model B Normal 25 30 psi_11 0.0174 0.0817 0.0802 0.0478 0.5963
57 Model B Normal 25 30 psi_21 0.0114 0.1121 0.1121 0.0893 0.7972
58 Model B Normal 25 30 psi_12 0.0175 0.1285 0.1279 0.0902 0.7051
59 Model B Normal 25 30 psi_22 0.0142 0.0805 0.0796 0.0468 0.5872
60 Model B Normal 25 30 pi_11 -0.0049 0.0295 0.0293 0.0335 1.1442
61 Model B Normal 25 30 pi_21 0.0031 0.0250 0.0250 0.0332 1.3307
62 Model B Normal 25 30 pi_12 -0.0005 0.0331 0.0333 0.0334 1.0043
63 Model B Normal 25 30 pi_22 -0.0022 0.0368 0.0369 0.0333 0.9013
64 Model B Normal 49 100 a 0.0297 0.1124 0.1089 0.0260 0.2387
65 Model B Normal 49 100 psi_11 0.0227 0.0818 0.0790 0.0190 0.2401
66 Model B Normal 49 100 psi_21 0.0249 0.0960 0.0932 0.0369 0.3962
67 Model B Normal 49 100 psi_12 0.0167 0.0912 0.0901 0.0369 0.4097
68 Model B Normal 49 100 psi_22 0.0208 0.0808 0.0785 0.0191 0.2430
69 Model B Normal 49 100 pi_11 -0.0011 0.0154 0.0154 0.0129 0.8366
70 Model B Normal 49 100 pi_21 0.0005 0.0089 0.0090 0.0129 1.4406
71 Model B Normal 49 100 pi_12 0.0004 0.0078 0.0078 0.0129 1.6524
72 Model B Normal 49 100 pi_22 0.0004 0.0142 0.0143 0.0130 0.9072
73 Model B Normal 100 200 a 0.0094 0.0790 0.0788 0.0129 0.1641
74 Model B Normal 100 200 psi_11 0.0036 0.0276 0.0275 0.0095 0.3451
75 Model B Normal 100 200 psi_21 0.0065 0.0383 0.0379 0.0200 0.5268
76 Model B Normal 100 200 psi_12 0.0052 0.0379 0.0377 0.0200 0.5299
77 Model B Normal 100 200 psi_22 0.0048 0.0263 0.0260 0.0094 0.3601
78 Model B Normal 100 200 pi_11 0.0004 0.0115 0.0116 0.0064 0.5567
79 Model B Normal 100 200 pi_21 0.0003 0.0038 0.0038 0.0064 1.7029
80 Model B Normal 100 200 pi_12 0.0004 0.0039 0.0039 0.0064 1.6517
81 Model B Normal 100 200 pi_22 0.0001 0.0116 0.0117 0.0064 0.5516
82 Model B t3 25 30 a 0.0042 0.0754 0.0756 0.0669 0.8842
83 Model B t3 25 30 psi_11 0.0002 0.0457 0.0459 0.0478 1.0411
84 Model B t3 25 30 psi_21 -0.0152 0.1272 0.1269 0.1162 0.9151
85 Model B t3 25 30 psi_12 0.0163 0.1217 0.1212 0.1171 0.9661
86 Model B t3 25 30 psi_22 -0.0047 0.0466 0.0466 0.0479 1.0269
87 Model B t3 25 30 pi_11 -0.0042 0.0311 0.0310 0.0337 1.0890
88 Model B t3 25 30 pi_21 -0.0005 0.0317 0.0318 0.0336 1.0550
89 Model B t3 25 30 pi_12 -0.0042 0.0338 0.0337 0.0337 0.9991
90 Model B t3 25 30 pi_22 -0.0018 0.0362 0.0363 0.0336 0.9246
91 Model B t3 49 100 a 0.0023 0.0293 0.0294 0.0252 0.8585
92 Model B t3 49 100 psi_11 0.0004 0.0196 0.0197 0.0188 0.9531
93 Model B t3 49 100 psi_21 -0.0019 0.0482 0.0484 0.0541 1.1187
94 Model B t3 49 100 psi_12 0.0025 0.0469 0.0471 0.0541 1.1475
95 Model B t3 49 100 psi_22 0.0002 0.0185 0.0186 0.0186 1.0034
96 Model B t3 49 100 pi_11 -0.0028 0.0148 0.0146 0.0130 0.8890
97 Model B t3 49 100 pi_21 0.0001 0.0079 0.0080 0.0130 1.6267
98 Model B t3 49 100 pi_12 0.0000 0.0060 0.0061 0.0130 2.1446
99 Model B t3 49 100 pi_22 -0.0006 0.0132 0.0132 0.0130 0.9821
100 Model B t3 100 200 a 0.0111 0.0255 0.0230 0.0124 0.5390
101 Model B t3 100 200 psi_11 0.0006 0.0090 0.0090 0.0097 1.0808
102 Model B t3 100 200 psi_21 0.0011 0.0260 0.0261 0.0304 1.1629
103 Model B t3 100 200 psi_12 0.0008 0.0270 0.0272 0.0303 1.1175
104 Model B t3 100 200 psi_22 -0.0004 0.0085 0.0085 0.0097 1.1455
105 Model B t3 100 200 pi_11 0.0024 0.0078 0.0075 0.0065 0.8593
106 Model B t3 100 200 pi_21 -0.0002 0.0023 0.0023 0.0065 2.7530
107 Model B t3 100 200 pi_12 0.0006 0.0029 0.0029 0.0065 2.2443
108 Model B t3 100 200 pi_22 0.0005 0.0074 0.0074 0.0065 0.8681
109 Model C Normal 25 30 a -0.0598 0.1432 0.1308 0.1154 0.8826
110 Model C Normal 25 30 psi_11 0.0011 0.0575 0.0578 0.0574 0.9935
111 Model C Normal 25 30 psi_21 -0.0153 0.0755 0.0743 0.0729 0.9809
112 Model C Normal 25 30 psi_12 -0.0001 0.0790 0.0794 0.0734 0.9250
113 Model C Normal 25 30 psi_22 -0.0097 0.0574 0.0568 0.0571 1.0053
114 Model C Normal 25 30 pi_11 -0.0053 0.0316 0.0314 0.0336 1.0732
115 Model C Normal 25 30 pi_21 -0.0007 0.0274 0.0276 0.0327 1.1872
116 Model C Normal 25 30 pi_12 -0.0029 0.0352 0.0352 0.0331 0.9383
117 Model C Normal 25 30 pi_22 -0.0058 0.0368 0.0365 0.0334 0.9156
118 Model C Normal 49 100 a -0.0157 0.0465 0.0440 0.0448 1.0176
119 Model C Normal 49 100 psi_11 -0.0017 0.0261 0.0262 0.0230 0.8774
120 Model C Normal 49 100 psi_21 -0.0001 0.0326 0.0328 0.0288 0.8791
121 Model C Normal 49 100 psi_12 -0.0003 0.0300 0.0302 0.0288 0.9555
122 Model C Normal 49 100 psi_22 -0.0047 0.0209 0.0205 0.0230 1.1222
123 Model C Normal 49 100 pi_11 -0.0033 0.0133 0.0130 0.0129 0.9969
124 Model C Normal 49 100 pi_21 0.0003 0.0121 0.0121 0.0127 1.0482
125 Model C Normal 49 100 pi_12 0.0013 0.0113 0.0113 0.0127 1.1213
126 Model C Normal 49 100 pi_22 -0.0020 0.0133 0.0132 0.0130 0.9840
127 Model C Normal 100 200 a 0.0016 0.0253 0.0254 0.0227 0.8942
128 Model C Normal 100 200 psi_11 -0.0001 0.0135 0.0135 0.0119 0.8779
129 Model C Normal 100 200 psi_21 -0.0003 0.0181 0.0182 0.0148 0.8117
130 Model C Normal 100 200 psi_12 0.0006 0.0165 0.0166 0.0148 0.8929
131 Model C Normal 100 200 psi_22 -0.0007 0.0135 0.0136 0.0119 0.8761
132 Model C Normal 100 200 pi_11 -0.0009 0.0070 0.0070 0.0064 0.9180
133 Model C Normal 100 200 pi_21 0.0000 0.0050 0.0050 0.0063 1.2678
134 Model C Normal 100 200 pi_12 0.0000 0.0047 0.0047 0.0063 1.3380
135 Model C Normal 100 200 pi_22 -0.0009 0.0082 0.0082 0.0064 0.7871
136 Model C t3 25 30 a 0.0238 0.1009 0.0985 0.0775 0.7862
137 Model C t3 25 30 psi_11 -0.0046 0.0542 0.0542 0.0577 1.0644
138 Model C t3 25 30 psi_21 -0.0148 0.0791 0.0781 0.0825 1.0574
139 Model C t3 25 30 psi_12 0.0095 0.0983 0.0983 0.0831 0.8458
140 Model C t3 25 30 psi_22 -0.0167 0.0690 0.0673 0.0573 0.8513
141 Model C t3 25 30 pi_11 -0.0046 0.0307 0.0305 0.0336 1.1021
142 Model C t3 25 30 pi_21 -0.0017 0.0305 0.0307 0.0328 1.0703
143 Model C t3 25 30 pi_12 -0.0044 0.0300 0.0299 0.0330 1.1048
144 Model C t3 25 30 pi_22 -0.0054 0.0374 0.0372 0.0335 0.9001
145 Model C t3 49 100 a 0.0028 0.0311 0.0311 0.0284 0.9132
146 Model C t3 49 100 psi_11 -0.0017 0.0237 0.0237 0.0232 0.9793
147 Model C t3 49 100 psi_21 0.0001 0.0318 0.0320 0.0333 1.0404
148 Model C t3 49 100 psi_12 0.0007 0.0294 0.0295 0.0334 1.1299
149 Model C t3 49 100 psi_22 -0.0013 0.0208 0.0209 0.0232 1.1132
150 Model C t3 49 100 pi_11 -0.0029 0.0144 0.0142 0.0130 0.9128
151 Model C t3 49 100 pi_21 -0.0002 0.0063 0.0064 0.0127 2.0019
152 Model C t3 49 100 pi_12 -0.0014 0.0069 0.0068 0.0127 1.8623
153 Model C t3 49 100 pi_22 -0.0009 0.0120 0.0120 0.0130 1.0788
154 Model C t3 100 200 a 0.0015 0.0233 0.0234 0.0142 0.6060
155 Model C t3 100 200 psi_11 -0.0003 0.0146 0.0147 0.0120 0.8146
156 Model C t3 100 200 psi_21 -0.0021 0.0187 0.0186 0.0171 0.9193
157 Model C t3 100 200 psi_12 0.0004 0.0197 0.0197 0.0171 0.8659
158 Model C t3 100 200 psi_22 -0.0028 0.0142 0.0140 0.0120 0.8524
159 Model C t3 100 200 pi_11 -0.0012 0.0089 0.0089 0.0064 0.7272
160 Model C t3 100 200 pi_21 0.0003 0.0034 0.0034 0.0063 1.8814
161 Model C t3 100 200 pi_12 0.0007 0.0034 0.0033 0.0063 1.9146
162 Model C t3 100 200 pi_22 -0.0026 0.0090 0.0086 0.0064 0.7481
coverage_H mean_SE_S SE_ratio_S coverage_S coverage_old
1 0.9100 0.1108 0.8084 0.88 1.0000
2 0.9600 0.0403 0.9730 0.93 1.0000
3 0.9800 0.0571 0.9932 0.97 1.0000
4 0.9400 0.0586 0.9322 0.94 1.0000
5 0.9300 0.0396 0.9282 0.92 1.0000
6 0.9700 0.0312 1.0262 0.95 1.0000
7 0.9400 0.0314 1.0811 0.95 1.0000
8 0.9400 0.0326 0.9306 0.94 1.0000
9 0.9400 0.0301 0.8709 0.93 1.0000
10 0.9200 0.0425 0.9715 0.93 1.0000
11 0.9100 0.0169 0.9460 0.91 1.0000
12 0.9600 0.0223 0.9687 0.95 1.0000
13 0.9500 0.0224 0.9749 0.95 1.0000
14 0.9700 0.0169 1.1050 0.96 1.0000
15 0.9400 0.0121 0.9602 0.93 1.0000
16 0.9600 0.0130 1.0413 0.94 1.0000
17 0.9900 0.0128 1.3368 0.97 1.0000
18 0.9500 0.0122 1.0104 0.92 1.0000
19 0.9400 0.0217 0.8066 0.93 1.0000
20 0.9100 0.0089 0.8292 0.91 1.0000
21 0.9100 0.0116 0.8993 0.91 1.0000
22 0.9300 0.0115 0.9987 0.93 1.0000
23 0.9300 0.0087 0.8931 0.92 1.0000
24 0.9500 0.0061 0.9082 0.95 1.0000
25 0.9600 0.0064 1.2450 0.96 1.0000
26 0.9800 0.0065 1.2283 0.98 1.0000
27 0.8900 0.0061 0.7784 0.89 1.0000
28 0.8900 0.0811 0.7933 0.86 1.0000
29 0.9200 0.0410 0.9819 0.93 1.0000
30 0.9300 0.0681 0.9268 0.92 1.0000
31 0.9300 0.0676 0.8734 0.94 1.0000
32 0.9200 0.0404 0.8638 0.91 1.0000
33 0.9600 0.0308 1.0005 0.93 1.0000
34 0.9600 0.0322 1.1046 0.96 1.0000
35 0.9400 0.0332 1.1802 0.94 1.0000
36 0.9300 0.0302 0.8697 0.91 1.0000
37 0.9000 0.0293 0.9052 0.91 1.0000
38 0.9000 0.0169 0.9458 0.91 1.0000
39 0.9800 0.0273 1.0360 0.97 1.0000
40 0.9700 0.0271 1.0624 0.95 1.0000
41 0.9400 0.0168 1.0393 0.94 1.0000
42 0.9100 0.0122 0.9122 0.91 1.0000
43 0.9800 0.0129 1.5711 0.98 1.0000
44 0.9800 0.0130 1.8109 0.98 1.0000
45 0.9400 0.0122 0.9767 0.92 1.0000
46 0.7600 0.0147 0.5833 0.76 1.0000
47 0.9300 0.0087 0.9636 0.94 1.0000
48 0.9600 0.0139 1.0265 0.96 1.0000
49 0.9500 0.0139 1.0794 0.95 1.0000
50 0.9500 0.0087 1.0200 0.96 1.0000
51 0.8400 0.0061 0.7208 0.84 1.0000
52 0.9900 0.0064 2.2654 0.99 1.0000
53 0.9900 0.0064 2.1143 0.99 1.0000
54 0.8900 0.0061 0.7079 0.89 1.0000
55 0.9100 0.0630 0.5432 0.90 1.0000
56 0.9300 0.0487 0.6071 0.89 1.0000
57 0.8600 0.1001 0.8930 0.88 1.0000
58 0.8400 0.1015 0.7937 0.87 1.0000
59 0.9000 0.0465 0.5835 0.89 1.0000
60 0.9800 0.0323 1.1036 0.98 1.0000
61 0.9700 0.0311 1.2450 0.97 1.0000
62 0.9300 0.0321 0.9641 0.94 1.0000
63 0.9400 0.0319 0.8655 0.87 1.0000
64 0.8500 0.0243 0.2232 0.84 0.9800
65 0.8485 0.0219 0.2773 0.87 1.0000
66 0.8384 0.0463 0.4967 0.87 0.9798
67 0.8586 0.0461 0.5111 0.88 0.9798
68 0.8700 0.0223 0.2840 0.89 1.0000
69 0.9000 0.0129 0.8372 0.91 1.0000
70 0.9800 0.0128 1.4257 0.97 1.0000
71 0.9700 0.0128 1.6427 0.97 1.0000
72 0.9400 0.0130 0.9131 0.93 1.0000
73 0.8800 0.0125 0.1590 0.86 0.9900
74 0.8800 0.0100 0.3658 0.88 0.9900
75 0.8788 0.0223 0.5888 0.86 1.0000
76 0.8788 0.0223 0.5898 0.93 1.0000
77 0.8900 0.0095 0.3640 0.90 0.9900
78 0.9000 0.0065 0.5598 0.89 1.0000
79 0.9800 0.0064 1.7033 0.98 1.0000
80 0.9800 0.0065 1.6589 0.99 1.0000
81 0.8700 0.0065 0.5578 0.88 1.0000
82 0.8900 0.0706 0.9335 0.92 1.0000
83 0.9400 0.0478 1.0402 0.95 1.0000
84 0.9100 0.1228 0.9676 0.91 1.0000
85 0.9100 0.1237 1.0203 0.91 1.0000
86 0.9700 0.0483 1.0355 0.94 1.0000
87 0.9600 0.0323 1.0424 0.92 1.0000
88 0.9300 0.0318 1.0000 0.93 1.0000
89 0.9400 0.0332 0.9861 0.93 1.0000
90 0.9200 0.0326 0.8986 0.93 1.0000
91 0.9100 0.0258 0.8798 0.91 1.0000
92 0.9200 0.0191 0.9677 0.92 1.0000
93 0.9800 0.0544 1.1243 0.95 1.0000
94 0.9900 0.0544 1.1553 0.96 1.0000
95 0.9400 0.0189 1.0168 0.96 1.0000
96 0.8900 0.0128 0.8785 0.88 1.0000
97 0.9800 0.0128 1.6000 0.97 1.0000
98 0.9800 0.0129 2.1319 0.99 1.0000
99 0.9400 0.0130 0.9811 0.93 1.0000
100 0.7500 0.0127 0.5533 0.77 1.0000
101 0.9800 0.0099 1.1021 1.00 1.0000
102 0.9900 0.0300 1.1505 0.97 1.0000
103 0.9700 0.0301 1.1080 0.98 1.0000
104 0.9800 0.0099 1.1592 0.98 1.0000
105 0.8800 0.0065 0.8588 0.89 1.0000
106 0.9900 0.0064 2.7301 0.98 1.0000
107 0.9900 0.0064 2.2426 0.98 1.0000
108 0.9300 0.0064 0.8620 0.90 1.0000
109 0.9300 0.1099 0.8399 0.88 1.0000
110 0.9400 0.0552 0.9551 0.92 1.0000
111 0.9600 0.0771 1.0383 0.96 1.0000
112 0.9200 0.0793 0.9994 0.96 1.0000
113 0.9200 0.0560 0.9858 0.92 1.0000
114 0.9600 0.0331 1.0566 0.97 1.0000
115 0.9600 0.0310 1.1248 0.96 1.0000
116 0.9100 0.0315 0.8940 0.92 1.0000
117 0.9300 0.0325 0.8911 0.91 1.0000
118 0.9300 0.0424 0.9640 0.92 1.0000
119 0.9300 0.0234 0.8934 0.93 1.0000
120 0.9200 0.0313 0.9538 0.94 1.0000
121 0.9400 0.0311 1.0299 0.94 1.0000
122 0.9700 0.0230 1.1222 0.97 1.0000
123 0.9400 0.0130 1.0007 0.94 1.0000
124 0.9200 0.0127 1.0504 0.93 1.0000
125 0.9500 0.0127 1.1184 0.97 1.0000
126 0.9300 0.0131 0.9958 0.92 1.0000
127 0.9100 0.0217 0.8570 0.89 1.0000
128 0.8800 0.0121 0.8952 0.90 1.0000
129 0.8500 0.0160 0.8773 0.92 1.0000
130 0.9200 0.0158 0.9534 0.92 1.0000
131 0.9300 0.0118 0.8720 0.94 1.0000
132 0.9500 0.0065 0.9235 0.95 1.0000
133 0.9600 0.0063 1.2692 0.96 1.0000
134 0.9800 0.0063 1.3417 0.99 1.0000
135 0.8800 0.0065 0.7943 0.90 1.0000
136 0.8800 0.0791 0.8023 0.89 1.0000
137 0.9500 0.0552 1.0174 0.93 1.0000
138 0.9700 0.0797 1.0215 0.91 1.0000
139 0.9200 0.0806 0.8205 0.87 1.0000
140 0.9100 0.0552 0.8202 0.90 1.0000
141 0.9500 0.0326 1.0694 0.92 1.0000
142 0.9600 0.0313 1.0219 0.95 1.0000
143 0.9700 0.0322 1.0786 0.92 1.0000
144 0.9100 0.0321 0.8641 0.89 1.0000
145 0.9200 0.0291 0.9351 0.92 1.0000
146 0.9500 0.0229 0.9640 0.94 1.0000
147 0.9700 0.0324 1.0130 0.96 1.0000
148 0.9800 0.0322 1.0920 0.98 1.0000
149 0.9500 0.0229 1.0991 0.95 1.0000
150 0.9000 0.0130 0.9147 0.91 1.0000
151 1.0000 0.0126 1.9843 1.00 1.0000
152 1.0000 0.0127 1.8585 0.99 1.0000
153 0.9500 0.0129 1.0761 0.94 1.0000
154 0.8000 0.0148 0.6310 0.80 1.0000
155 0.9000 0.0119 0.8123 0.90 1.0000
156 0.9400 0.0165 0.8857 0.94 1.0000
157 0.9000 0.0165 0.8367 0.90 1.0000
158 0.9000 0.0119 0.8476 0.90 1.0000
159 0.8800 0.0064 0.7283 0.89 1.0000
160 0.9800 0.0063 1.8742 0.98 1.0000
161 0.9900 0.0063 1.9109 0.99 1.0000
162 0.9000 0.0064 0.7483 0.91 1.0000model_cols <- c(
"Model A" = "#D55E00",
"Model B" = "#009E73",
"Model C" = "#0072B2"
)
parameter_labeller <- as_labeller(
c(
"a" = "a",
"pi_11" = "pi[11]",
"pi_12" = "pi[12]",
"pi_21" = "pi[21]",
"pi_22" = "pi[22]",
"psi_11" = "psi[11]",
"psi_12" = "psi[12]",
"psi_21" = "psi[21]",
"psi_22" = "psi[22]"
),
label_parsed
)
error_labeller <- as_labeller(
c("Normal" = "Normal", "t3" = "t[3]"),
label_parsed
)
plot_MC_measure <- function(data, measure, ylab, reference = NULL, ylim = NULL) {
p <- ggplot(
data,
aes(
x = n,
y = .data[[measure]],
colour = model,
group = model
)
) +
geom_line(linewidth = 0.8, alpha = 0.95) +
geom_point(aes(shape = T_fac), size = 2.5, stroke = 0.7) +
facet_grid(
error ~ parameter,
scales = "free_y",
labeller = labeller(
error = error_labeller,
parameter = parameter_labeller
)
) +
scale_colour_manual(values = model_cols) +
scale_shape_manual(values = c("30" = 16, "100" = 17, "200" = 15)) +
labs(
x = "Number of spatial locations",
y = ylab,
colour = "Model",
shape = "Number of time points"
) +
theme_minimal(base_size = 11) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
panel.grid.major.x = element_line(colour = "grey88", linewidth = 0.35),
panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.35),
strip.background = element_rect(fill = "grey95", colour = NA),
strip.text = element_text(face = "bold", size = 9)
)
if (!is.null(reference)) {
p <- p + geom_hline(
yintercept = reference,
linetype = "dashed",
colour = "grey35",
linewidth = 0.55
)
}
if (!is.null(ylim)) p <- p + coord_cartesian(ylim = ylim)
p
}if (have_mc_results) {
p_rmse <- plot_MC_measure(
MC_summary,
measure = "rmse",
ylab = "RMSE"
)
p_rmse
}
if (have_mc_results) {
p_ratio_h <- plot_MC_measure(
MC_summary,
measure = "SE_ratio_H",
ylab = "Mean SE / Monte Carlo SD",
reference = 1
)
p_ratio_h
}
if (have_mc_results) {
p_cov_h <- plot_MC_measure(
MC_summary,
measure = "coverage_H",
ylab = "Empirical coverage probability",
reference = 0.95,
ylim = c(0.75, 1)
)
p_cov_h
}
ggsave(
"Results_MC_RMSE.pdf",
p_rmse,
width = 14,
height = 7.5,
units = "in",
device = cairo_pdf
)
ggsave(
"Results_MC_SE_ratio.pdf",
p_ratio_h,
width = 14,
height = 7.5,
units = "in",
device = cairo_pdf
)
ggsave(
"Results_MC_Coverage.pdf",
p_cov_h,
width = 14,
height = 7.5,
units = "in",
device = cairo_pdf
)For Model B, the absence of temporal dependence makes the separation of the constant intercept and spatial autoregressive parameters more reliant on the spatial covariance structure. The following diagnostic reproduces the empirical correlation matrix among the estimated intercept and the four entries of \(\Psi\) for the largest normal-error design.
if (have_mc_results) {
theta_B <- do.call(
rbind,
lapply(
successful_results(M2[[5]]),
function(x) x$output$theta
)
)
round(cor(theta_B[, 1:5]), 3)
} [,1] [,2] [,3] [,4] [,5]
[1,] 1.000 0.770 0.483 0.521 0.734
[2,] 0.770 1.000 0.131 0.745 0.617
[3,] 0.483 0.131 1.000 -0.317 0.679
[4,] 0.521 0.745 -0.317 1.000 0.249
[5,] 0.734 0.617 0.679 0.249 1.000The Berlin application uses \(p=3\) variable-specific intercepts, so the estimator is called with intercept = "vector". The current document does not include the original Berlin data. Once the empirical response array Y_emp and the corresponding row-standardised spatial weights matrix W_emp are available, the corrected estimator is run as follows.
parameters_emp <- list(errortype = "norm")
est_emp <- qml_vec_spARCH(
Y = Y_emp,
W = W_emp,
parameters = parameters_emp,
intercept = "vector"
)
p <- dim(Y_emp)[2]
round(est_emp$a_tilde_est, 3)
round(
matrix(est_emp$theta_tilde[(p + 1):(p^2 + p)], p, p),
3
)
round(
matrix(est_emp$theta_tilde[(p^2 + p + 1):(2 * p^2 + p)], p, p),
3
)
round(est_emp$stderr_hessian, 3)
round(est_emp$stderr_sandwich, 3)The corrigendum reports the estimated centred intercept \(\widetilde A\) together with \(\Psi\) and \(\Pi\). The corresponding estimates are approximately
\[ \widehat{\widetilde A} = (-4.731,\;0.226,\;-2.544)^\prime, \]
\[ \widehat\Psi = \begin{pmatrix} 0.108 & 0.019 & -0.054\\ 0.014 & 0.140 & 0.000\\ -0.084 & 0.013 & 0.118 \end{pmatrix}, \]
and
\[ \widehat\Pi = \begin{pmatrix} 0.581 & 0.128 & -0.014\\ 0.082 & 0.552 & 0.026\\ -0.027 & 0.080 & 0.610 \end{pmatrix}. \]
The corrected Hessian-based standard errors are much smaller than those reported under the original likelihood normalisation. All displayed dependence coefficients are statistically significant at the 5% level except the spatial cross-effect whose estimate is essentially zero.
Let
\[ \mathbf B = \mathbf S_{np}^{-1}(\mathbf\Pi^\prime\otimes\mathbf I_n), \qquad \mathbf c = \mathbf S_{np}^{-1}\operatorname{vec}(\widetilde{\mathbf A}), \qquad \mathbf e_t = \mathbf S_{np}^{-1}\operatorname{vec}(\mathbf U_t). \]
Then
\[ \ddot{\mathbf Y}_t = \mathbf c+\mathbf B\ddot{\mathbf Y}_{t-1}+\mathbf e_t, \]
and, if \(\rho(\mathbf B)<1\),
\[ \ddot{\mathbf Y}_t = (\mathbf I_{np}-\mathbf B)^{-1}\mathbf c + \sum_{k=0}^{\infty}\mathbf B^k\mathbf e_{t-k}. \]
A convenient sufficient condition follows from any submultiplicative matrix norm:
\[ \rho\!\left(\mathbf S_{np}^{-1}(\mathbf\Pi^\prime\otimes\mathbf I_n)\right) \le \|\mathbf S_{np}^{-1}\| \|\mathbf\Pi^\prime\otimes\mathbf I_n\| <1. \]
For \(\Pi\), a Gershgorin-based sufficient condition is
\[ \max_i\left\{|\pi_{ii}|+\sum_{j\ne i}|\pi_{ij}|\right\}<1. \]
Identification additionally presumes a non-degenerate spatial structure, so that distinct admissible values of \(\Psi\) induce distinct spatial covariance structures. Degenerate cases such as \(W=0\), for which \(\Psi\) is not identifiable, are excluded.
To reproduce the complete results:
mc_replications: 1000 and the desired value of standardise_t.Output_vecspARCH_MC_1000.rda in the same directory as this QMD.sessionInfo()R version 4.6.1 (2026-06-24)
Platform: aarch64-apple-darwin23
Running under: macOS Tahoe 26.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] tidyr_1.3.2 dplyr_1.2.1 ggplot2_4.0.3
[4] numDeriv_2016.8-1.1 Rsolnp_2.0.1 spdep_1.4-2
[7] sf_1.1-2 spData_2.3.5
loaded via a namespace (and not attached):
[1] s2_1.1.11 future_1.75.0 generics_0.1.4
[4] class_7.3-24 KernSmooth_2.23-27 lattice_0.22-9
[7] listenv_1.0.0 digest_0.6.39 magrittr_2.0.5
[10] evaluate_1.0.5 grid_4.6.1 RColorBrewer_1.1-3
[13] fastmap_1.2.0 jsonlite_2.0.0 e1071_1.7-17
[16] DBI_1.3.0 purrr_1.2.2 scales_1.4.0
[19] truncnorm_1.0-9 codetools_0.2-20 cli_3.6.6
[22] rlang_1.3.0 units_1.0-1 parallelly_1.48.0
[25] future.apply_1.20.2 withr_3.0.3 yaml_2.3.12
[28] otel_0.2.0 tools_4.6.1 deldir_2.0-4
[31] boot_1.3-32 globals_0.19.1 vctrs_0.7.3
[34] R6_2.6.1 proxy_0.4-29 lifecycle_1.0.5
[37] classInt_0.4-11 htmlwidgets_1.6.4 pkgconfig_2.0.3
[40] pillar_1.11.1 gtable_0.3.6 glue_1.8.1
[43] Rcpp_1.1.2 tidyselect_1.2.1 xfun_0.60
[46] tibble_3.3.1 knitr_1.51 farver_2.1.2
[49] htmltools_0.5.9 labeling_0.4.3 rmarkdown_2.31
[52] wk_0.9.5 compiler_4.6.1 S7_0.2.2
[55] sp_2.2-3 Otto, P. (2024). A multivariate spatial and spatiotemporal ARCH model. Spatial Statistics, 60, 100823.