The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
83,722 characters
Vine Copula VAR: From Recursive Margins to Joint Forecast Inference
\title{Vine Copula VAR:\\ From Recursive Margins to Joint Forecast Inference\thanks{Previous versions of this paper have been circulated under the title ``Vine Copula VAR". We thank seminar participants at the 2026 Midwest Econometrics Group Meeting and the CUNY Graduate Center for their valuable comments. Yubo Tao acknowledges the financial support from the National Natural Science Foundation of China (Project No.~72303003) and a seed grant from the Asia-Pacific Academy of Economics and Management of the University of Macau.}}
\author[a]{Hunter Ng}
\affil[a]{\small\emph{Zicklin School of Business, Baruch College, City University of New York, USA}}
\author[b]{Yubo Tao}
\affil[b]{\small\emph{Faculty of Social Sciences, University of Macau, Macao SAR, China}}
\date{October 2026}
\maketitle
\begin{abstract}
\noindent
Joint-event forecasts often combine a dependence estimate based on past forecast errors with newly estimated marginal distributions. When each historical error retains the marginal fit available at its issue date, inference must account for an overlapping sequence of estimation errors. We derive their joint influence with the terminal forecast estimates in
a stable Vine Copula VAR with normal innovation margins and a fixed,
correctly specified Gaussian or positive Clayton vine. An intercept
identity and the stable VAR filter reduce the historical correction to
harmonically weighted innovation moments, while terminal slope uncertainty
remains. The resulting covariance estimator gives asymptotically valid
repeated-sample intervals for fixed one-sided event probabilities at the
realized forecast state. In the Gaussian submodel, retaining issued
transforms adds a positive semidefinite covariance term relative to
refitting margins on the same observations. Monte Carlo simulations show
that terminal-margin uncertainty is quantitatively more important than
this additional term and that logit intervals improve lower-tail coverage
in the designs studied. A real-time forecasting application to U.S.
macroeconomic releases shows how marginal estimation contributes to
uncertainty in predicted probabilities of joint contractions and identifies
limitations of the stationary marginal model.
\par\medskip
\noindent\textbf{Keywords}: Vine copulas; Recursive estimation; Generated regressors; Joint-event forecasts; Tail dependence.
\par\medskip
\noindent\textbf{JEL Classification}: C32; C53; C58.
\end{abstract}
\clearpage
\section{Introduction}\label{sec:introduction}
A forecaster assessing the probability of a joint contraction in
employment and production needs more than forecasts of each series and
their conditional covariance. The probability also depends on how
adverse innovations occur together. A Gaussian VAR ties innovation
dependence to a correlation matrix and imposes symmetry between joint
lower and upper tails. Even with normal marginal forecast errors, a
joint distribution can exhibit asymmetric dependence. A Vine Copula VAR
retains the familiar linear dynamics of conditional means while
allowing non-Gaussian contemporaneous dependence.
Vines provide a tractable way to specify this dependence. They build
a multivariate copula from bivariate copulas for pairs of variables
and conditional pairs \citep{aas_pair-copula_2009,nagler_kruger_min_2022}.
Different pairs can have different strengths and forms of dependence,
so joint downside risk need not be governed by a common symmetric
family. Gaussian pairs provide a symmetric benchmark, whereas positive
Clayton pairs allow lower-tail dependence \citep{czado_nagler_2022}.
Both families are covered by our regularity results. The pairwise
construction also supports sequential estimation, with later fits using
conditional transforms obtained from earlier pairs. These features make
the vine useful for translating marginal VAR forecasts into joint-event
probabilities.
We study inference for those probabilities when the dependence model
uses a record of marginal forecasts issued over time. Historical
transforms then retain different marginal fits, whereas the next
forecast uses the current fit. Their estimation errors overlap and
also enter the vine's conditional transforms. We derive their joint
influence and construct feasible confidence intervals for conditional
event probabilities.
In a joint-contraction application, the
training observations record whether employment and production were
simultaneously low relative to their preceding forecast distributions.
Each probability integral transform (PIT) is computed from a marginal
fit that excludes the outcome being transformed. Retaining it preserves
this timing convention. Recomputing the transform from the terminal fit
would let the outcome being assessed help determine its own reference
distribution. For the next, unobserved outcome, all current observations
are available to estimate the margins. Updating terminal margins and
retaining historical transforms therefore serve the same objective:
each outcome is standardized using information available before it
arrives.
This training protocol is motivated by the sequential assessment of
issued forecasts in prequential analysis \citep{dawid_1984} and by the
PIT approach to density-forecast evaluation
\citep{diebold_gunther_tay_1998}. We apply that logic to the marginal
inputs of copula estimation. The raw history remains available, so a
terminal refit on the same observations is feasible and provides an
efficiency benchmark. Under the maintained model, both procedures
target the same constant innovation copula and the same conditional
event probability.
Our specification uses a stable VAR with constant innovation scales
and a fixed vine for contemporaneous innovation dependence. Marginal
equations are estimated on expanding samples, the vine is fitted
sequentially to the retained
PITs, and the next event probability combines the vine estimate with the
terminal VAR estimates. The general need to propagate first-stage error
is familiar from two-step econometric inference
\citep{murphy_topel_1985} and copula estimation with fitted margins
\citep{joe_2005,patton_2012}. Here the first stage is an entire sequence
of overlapping marginal fits. An observation affects many subsequent
PITs, and the resulting errors also pass through the conditional
transforms used in later vine trees.
The main argument identifies which of these errors survive averaging.
An intercept identity and the stable VAR filter cancel the centered
historical lag-regressor contribution, leaving harmonic weights on
innovation location and scale moments. Earlier observations receive
larger weights because they enter more subsequent regressions. This
reduction applies to the historical scores. The next forecast is
evaluated at a single lag state, so every terminal slope remains in the
forecast gradient. We derive the joint covariance of the sequential
vine estimate, all terminal coefficients, and the terminal scales,
including their cross covariances. The resulting estimator is positive
semidefinite and uses quantities available from the forecast record.
Inference also requires care because the terminal state remains random
as the sample grows. We separate its recent innovations from the leading
estimation error and obtain a state-dependent Gaussian-mixture limit.
Studentization then yields a normal pivot. The intervals have
repeated-sample coverage for the conditional probability evaluated at
each sample's realized state. The theory is pointwise under a fixed
stationary law, with fixed dimension and lag order, iid innovation
vectors, normal coordinate margins, and a correctly specified simplified
vine. We verify the score and event conditions for Gaussian and positive
Clayton pairs, including their estimated ancestors. Fixed one-sided
events have positive limiting variance almost surely, including mixed
directions and subsets of variables. Bounded rectangles require an
event-specific nondegeneracy condition.
The Gaussian submodel permits a direct comparison between the archive
and a terminal refit on the same observations. Their limiting
terminal-margin and cross covariance blocks coincide. Retaining the
issued transforms adds a positive semidefinite term to the copula
block, which gives the variance cost of the archive rule for any covered
event. This comparison concerns a change in marginal standardization
on a common history. Replacing the history with revised observations
changes the data and is a separate comparison. We also derive a
covariance estimator that retains the exact finite training weights
and justify intervals on the log-odds scale. Both modifications are justified by the same first-order asymptotic theory.
Existing vine methods cover both conditional distributions and quantiles.
\citet{kraus_czado_2017} develop D-vine quantile regression, while
\citet{chang_joe_2019} compute conditional distributions under
regular-vine models with continuous and discrete variables. For
multivariate time series, \citet{smith_2015} uses a D-vine to represent
both nonlinear serial and cross-sectional dependence. Our question
concerns inference when the historical transforms retain
the marginal estimates available at their issue dates.
The closest results on inference for prediction are those of
\citet{nagler_kruger_min_2022}. Their stationary-vine framework uses
pair copulas to represent both serial and cross-sectional dependence,
extending copula time-series models such as
\citet{remillard_papageorgiou_soustra_2012}. Our vine describes
contemporaneous innovations after the VAR accounts for linear serial
dynamics. The inferential distinction concerns the marginal inputs.
Section~4.1 of
\citet{nagler_kruger_min_2022} forms historical transforms using a
common marginal estimate, and their Theorems~8--9 propagate marginal
and copula estimation error to prediction functionals. We calculate
the accumulated influence of dated marginal estimates and its joint
law with the terminal forecast parameters. The recursive weighting is
related to \citet{west_1996},
while the triangular score system builds on sequential vine estimation
\citep{hobaeck_haff_2013}. The sample-separation principle of
\citet{beutner_heinemann_smeekes_2021} addresses the additional issue
posed by a random conditioning state. We establish the required
separation for the recursive influence derived here.
Other copula time-series models address different sources of marginal
and dynamic uncertainty. \citet{chen_xiao_wang_2022} study filtered
nonstationarity, and \citet{chen_wang_xiao_yi_2025} develop a pseudo
sieve likelihood approach that jointly estimates the residual copula
parameter and invariant density. \citet{zhao_shi_zhang_2022} construct
conditional distributions through copula-linked univariate D-vines,
while \citet{tsionas_izzeldin_trapani_2022} allow time-varying VAR
parameters. Our fixed OLS and residual-variance rule is more restrictive
but makes the historical correction and its terminal cross covariances
explicit. Dimension remains fixed throughout. The common-fit,
increasing-dimension results of \citet{gauss_nagler_2026} concern a
different asymptotic regime.
Monte Carlo simulations assess interval coverage at 200--800
observations and use larger samples to examine the covariance
decomposition. For the lower-tail event, logit intervals give
93.8--96.1\% coverage in Gaussian vines and 93.3--95.2\% in Clayton
vines, compared with 85.8--92.6\% and 86.4--93.0\% for archive Wald
intervals. Some shortfalls remain, particularly with persistent dynamics,
and the smallest samples frequently activate the marginal safeguards.
Holding the fitted probability fixed while omitting covariance components
shows that terminal-margin uncertainty is much more consequential than
the additional archive term in these designs.
A real-time forecasting application examines joint contractions in U.S.
employment and economic activity at 450 historical origins. The dated
release record distinguishes changes in marginal estimates from revisions
to the underlying observations
\citep{croushore_stark_2001,koenig_dolmas_piger_2003}. Current marginal
estimation accounts for most of the adjustment to reported standard
errors. The archived vine has no consistent Brier-loss advantage over
the benchmarks, and the marginal diagnostics question the stationary,
constant-scale law even before the pandemic. The application therefore
shows how to quantify parameter uncertainty in a real-time forecast,
while identifying features that a richer marginal model would need to
address. Stochastic volatility and revision uncertainty are substantive
modeling issues in density forecasting
\citep{clark_2011,clements_galvao_2023}.
Section~\ref{sec:method} defines the forecast rule, and
Section~\ref{sec:theory} derives its joint limit and event intervals.
Sections~\ref{sec:simulation} and~\ref{sec:evaluation} present the
Monte Carlo simulations and real-time forecasting application. The main appendices contain
the central proofs. The online supplement provides supporting lemmas,
implementation details, and additional evidence.
\section{Model, forecast target, and estimation}\label{sec:method}\label{sec:model}
The construction assigns serial dynamics to the VAR and contemporaneous
innovation dependence to the vine. This separation follows the marginal
and copula approach to multivariate time series reviewed by
\citet{patton_2012}. The distinction needed for inference is temporal:
historical PITs use dated marginal estimates, whereas the event
probability uses the terminal estimates. We define both stages and the
information retained from them before stating the joint limit.
\subsection{A stationary VAR and a conditional event}
For fixed $n,p\ge1$, let $y_t=(y_{1t},\ldots,y_{nt})'$ follow the causal
stationary solution of
\begin{equation}\label{main:var-model}
y_t=a+\sum_{j=1}^p A_jy_{t-j}+D_\sigma Z_t,\qquad
D_\sigma=\text{diag}(\sigma_1,\ldots,\sigma_n),\quad \sigma_i>0.
\end{equation}
\begin{ass}[Stable VAR and innovations]\label{main:var-assumption}
The dimensions $n,p$ and parameters in \eqref{main:var-model} are fixed.
The lag companion matrix $\mathcal A$ and innovations satisfy
\[
\rho(\mathcal A)<1,\qquad Z_t\overset{\mathrm{iid}}{\sim}P_Z,
\qquad Z_{it}\sim N(0,1),\qquad R_Z=E(Z_tZ_t')\succ0.
\]
The process is its causal stationary solution, $\sigma_i>0$ for every
$i$, and the initial $p$ lag vectors are observed.
\end{ass}
Innovation vectors are independent across dates, while their coordinates
may be dependent. Thus Gaussian coordinate margins permit non-Gaussian
joint dependence, which the vine will specify. Write
$\mathcal F_t=\sigma(Z_s:s\le t)$. Stability makes the history of $y_t$
measurable with respect to this field. Define
\begin{equation}\label{eq:design}
L_t=(y_{t-1}',\ldots,y_{t-p}')',\qquad x_t=(1,L_t')',\qquad
y_{it}=x_t'\beta_i+\sigma_iZ_{it}.
\end{equation}
The vector $\beta_i$ contains an intercept and $np$ lag coefficients, so
its dimension is $d_x=1+np$. The stable VAR represents $L_t$ as a
geometrically decaying function of past innovations. Consequently, the
current innovation is independent of the regression design $x_t$, a
property used in the OLS influence calculation. The model's parameters
are constant. Its forecasts change with the observed lag state.
For a fixed deterministic threshold vector $c\in\mathbb R^n$, the
one-step forecast target is
\begin{equation}\label{main:conditional-target}
p_{T+1}(c)=P(y_{T+1}\le c\mid\mathcal F_T)
=C_{\theta_0}\left(\Phi\left\{
\frac{c_i-x_{T+1}'\beta_i}{\sigma_i}\right\}:i\le n\right).
\end{equation}
The inequalities hold coordinatewise. Although $c$ is fixed,
$p_{T+1}(c)$ is random across samples because the terminal state
$x_{T+1}$ is random. This distinction determines the coverage statement
in Section~\ref{main:forecast-section}. We first consider the lower
orthant in \eqref{main:conditional-target}, then extend inference to
fixed rectangles, mixed directions, and unrestricted coordinates. The
labor-market event in Section~\ref{sec:evaluation} uses that extension.
Estimating the threshold itself would introduce another influence term.
\subsection{Recursive margins and the archive}
At prefix $m$, estimate each marginal equation by guarded OLS and its
average squared residual:
\begin{align}
\widehat Q_m&=m^{-1}\sum_{s=1}^m x_sx_s',\notag\\
\widehat\beta_{i,m}
&=\begin{cases}
\widehat Q_m^{-1}m^{-1}\sum_{s=1}^m x_sy_{is},
&\lambda_{\min}(\widehat Q_m)\ge m^{-1/4},\\
0,&\text{otherwise},
\end{cases} \label{dyn:ols}\\
\widehat v_{i,m}
&=\max\left\{m^{-1/4},\,
m^{-1}\sum_{s=1}^m(y_{is}-x_s'\widehat\beta_{i,m})^2\right\},\notag\\
\widehat Z_{it}&=\frac{y_{it}-x_t'\widehat\beta_{i,t-1}}
{\sqrt{\widehat v_{i,t-1}}}.
\label{dyn:variance}
\end{align}
The Gram-matrix threshold and variance floor make the estimator defined
even for an ill-conditioned early regression. They become inactive
simultaneously on all retained prefixes with probability tending to one
(Lemma~\ref{dyn:first-stage}), but remain part of the procedure at
finite sample sizes. In particular, the residual $\widehat Z_{it}$
uses only estimates available through $t-1$.
The true and fitted conditional margins are
\begin{equation}\label{eq:gaussian-margin}
F^0_{i,t}(y)=\Phi\left(\frac{y-x_t'\beta_i}{\sigma_i}\right),\qquad
\widehat F_{i,t|t-1}(y)=
\Phi\left(\frac{y-x_t'\widehat\beta_{i,t-1}}
{\sqrt{\widehat v_{i,t-1}}}\right).
\end{equation}
Thus $y_t$ is assessed against a marginal distribution fitted before
it enters the regression. With $b_T=\lceil T^\eta\rceil$, $0<\eta<1$,
retain the archive
\begin{equation}\label{eq:stored-pits}
\widehat U_{it}=\widehat F_{i,t|t-1}(y_{it}),\qquad b_T<t\le T.
\end{equation}
Once recorded, a transform remains fixed. The current forecast margin
instead uses $\widehat\beta_{i,T}$ and $\widehat v_{i,T}$.
For each issue date, retain the marginal fit and the lag row, then
compute the PIT when the outcome arrives. The mathematical record
needed at origin $T$ is
\[
\mathcal A_T=\{(t,\widehat U_t,x_t):b_T<t\le T\},\qquad
\widehat Q_T,\quad
(\widehat\beta_{i,T},\widehat v_{i,T}:i\le n),\quad x_{T+1}.
\]
The PITs determine the score derivatives, while $x_t$ and
$\widehat Q_T$ determine the regression influences in
\eqref{main:joint-feasible}. Joint inference therefore requires more
than a PIT series. In numerical storage, retain the finite standardized
residuals $\widehat Z_t$ as well. Evaluating $\Phi$ can round an extreme
score to zero or one, after which inverse transformation cannot recover
it. Storing the scores preserves the information used by the estimator
without altering its mathematical definition.
The order of operations explains why the archive is retained while the
terminal margins are updated. On observing $y_t$, first evaluate it
using the marginal fit through $t-1$ and store $\widehat U_t$. Then
include $y_t$ in the marginal regression for the next forecast. At
origin $T$, the copula fit uses the resulting record of marginal
forecast surprises, while the forecast for $T+1$ uses the fit through
$T$. Thus every marginal training input was evaluated out of sample,
and the next forecast uses the full available estimation sample.
Replacing a stored PIT by a transform based on the fit through $T$
would change this training criterion even if the underlying observation
were held fixed.
This choice concerns the marginal inputs to estimation. The copula
parameters are estimated from the retained PITs, whose finite-sample
margins need not be exactly uniform because the dated marginal fits
contain estimation error. The population target remains the innovation
copula of the stationary model, and our inference accounts for that
generated error. A terminal refit targets the same copula using
in-sample marginal residuals. Section~\ref{main:dynamic-comparisons}
shows that refitting weakly reduces asymptotic variance in the Gaussian
submodel.
Retaining dated transforms implements the specified training protocol
and does not carry an efficiency claim.
For real-time data, let $\tau_m$ be the first date on which all
components of a declared reference-period vector are available, with
$\tau_m<\tau_{m+1}$. If $L_{i,m}^{[v]}$ denotes the level for period
$m$ available at date $v$, lock each growth transformation as
\begin{equation}\label{rli:locked-growth}
Y_{i,m}^{\rm L}=100\{\log L_{i,m}^{[\tau_m]}
-\log L_{i,m-1}^{[\tau_m]}\}.
\end{equation}
Both levels use the vintage available at $\tau_m$, so a growth rate
does not mix levels from different releases. Other declared
transformations complete the locked vector $Y_m^{\rm L}$. An early
component may have been revised before the final component appears.
The vector records the first common release rather than each component's
first estimate. Define the service history by
\begin{equation}\label{rli:service-field}
\mathcal G_m=\sigma(Y_s^{\rm L}:s\le m).
\end{equation}
Proposition~\ref{rli:transfer} applies the theory to these records if
the locked sequence itself satisfies \eqref{main:var-model} and the
innovation conditions. The recording rule does not establish that law.
The target conditions on $\mathcal G_T$. Using all calendar-time news
would require the next innovation to retain its maintained law given
that larger information set. Replacing the locked history with revised
data changes the observed process and lies outside the same-history
comparison. The stationary argument uses a two-sided history, while
implementation starts with observed initial lags.
Any fixed $\eta\in(0,1)$ is admissible in the proofs. The simulations use $\eta=0.7$ throughout. The growing prefix ensures
uniform control of retained marginal fits, while $b_T/T\to0$ preserves
the first-order sample normalization.
\subsection{Fixed vine and sequential fit}
For continuous predictive margins $F_{i,t}$ with densities $f_{i,t}$,
a valid copula density defines
\begin{equation}\label{eq:joint-model}
q_t(y\mid\mathcal F_{t-1})=
c_{\theta,V,K}\{F_{1,t}(y_1),\ldots,F_{n,t}(y_n)\}
\prod_{i=1}^n f_{i,t}(y_i).
\end{equation}
The regular
vine $V$, depth $K\in\{0,\ldots,n-1\}$, and pair families are fixed
before estimation. The simplified vine factorization is
\begin{equation}\label{eq:vine-product}
c_{\theta,V,K}(u)=\prod_{r=1}^K\prod_{e\in E_r(V)}
c_{e,\theta_e}(u_{a_e|D_e},u_{b_e|D_e}),\qquad |D_e|=r-1.
\end{equation}
An edge links $a_e$ and $b_e$ conditional on $D_e$. The arguments
are obtained by successive conditional-CDF maps
$h_{a|b;D}(v,w)=\partial_w C_{ab;D}(v,w)$, following the pair-copula
construction of \citet{aas_pair-copula_2009}. The simplifying restriction
makes the pair parameter invariant to the realized conditioning
variables. It is a dependence assumption, not a consequence of removing
conditional means and variances \citep{czado_nagler_2022}.
\begin{ass}[Fixed and correctly specified vine]\label{main:vine-assumption}
The copula of $Z_t$ is \eqref{eq:vine-product} at $\theta_0$.
The vine $V$, depth $K$, and pair families are fixed, and omitted
edges are independence copulas. For the retained free parameters,
\[
\Theta=\prod_{e=1}^E\Theta_e\subset\mathbb R^{q_\theta}
\text{ is compact},\qquad
\theta_0\in\operatorname{int}(\Theta),\qquad q_\theta<\infty.
\]
All these objects are fixed as $T\to\infty$.
Each free parameter belongs to one edge, and the product space permits
its variation independently of the other edge parameters. Pair densities
are positive on $(0,1)^2$ throughout their admissible parameter spaces.
\end{ass}
Correct specification includes the simplifying restriction and every
imposed independence edge. The product parameter space permits
stagewise estimation of the retained pairs, as in
\citet{hobaeck_haff_2013}. Parameters shared across edges would require
a different estimating system. Data-driven procedures such as
\citet{dissmann_selecting_2013} select the vine structure and pair
families. The present theory takes these choices, together with the
truncation depth, as fixed and does not account for selection uncertainty.
Positivity of each pair density makes the log objective
well defined in the open unit square. A uniform positive lower bound
is not imposed.
At origin $T$, enumerate retained edges in tree order and fit
\begin{equation}\label{eq:sequential-fit}
\widehat\theta_e\in\operatorname*{arg\,max}_{\theta_e\in\Theta_e}
\frac1{T-b_T}\sum_{t=b_T+1}^T
\log c_{e,\theta_e}
\{\widehat U_{a_e|D_e,t}(\widehat\theta_{<e}),
\widehat U_{b_e|D_e,t}(\widehat\theta_{<e})\}.
\end{equation}
Use fixed measurable tie-breaking. Optimization at edge $e$ holds the
estimated ancestors $\widehat\theta_{<e}$ fixed, but those estimates
determine the edge's conditional transforms. Their sampling error
therefore enters its score. The resulting system is triangular and
generally differs from the score of a jointly maximized likelihood.
Refitting at later origins updates estimates of the same constant
dependence law.
The fitted orthant probability combines this dependence estimate with
the terminal margins:
\begin{equation}\label{main:event-estimator}
\widehat p_{T+1}(c)=C_{\widehat\theta}\left(
\Phi\left\{\frac{c_i-x_{T+1}'\widehat\beta_{i,T}}
{\sqrt{\widehat v_{i,T}}}\right\}:i\le n\right).
\end{equation}
Theorem~\ref{main:event-inference} gives its interval. Section~\ref{main:implementation}
states the numerical accuracy needed for approximate optimization and
probability evaluation. Time-varying coefficients and latent volatility
require a separate analysis of the marginal estimation errors.
\section{Inference for conditional joint forecasts}\label{sec:theory}
The inferential problem has two parts. First, recursive marginal error
must be propagated through the sequential vine fit and joined to the
terminal VAR estimates. Second, that joint limit must be evaluated at
the random state used for the next forecast. The first part determines
the covariance matrix. The second justifies the event interval. We
develop them in that order, then use the same decomposition to study
terminal refitting, finite training weights, and logit intervals.
\subsection{Scores and regularity}
Let $V_e(\theta_{<e},z)$ be the two conditional PITs at edge $e$,
starting from coordinates $\Phi(z_i)$. Define the edge objective and
its own-parameter score by
\begin{equation}\label{main:edge-scores}
\ell_e(\theta_{\le e},z)=\log c_{e,\theta_e}\{V_e(\theta_{<e},z)\},
\qquad g_e(\theta,z)=\partial_{\theta_e}\ell_e(\theta_{\le e},z).
\end{equation}
The stacked score $g=(g_1',\ldots,g_E')'$ has dimension
$q_\theta=\dim(\theta)$. Ancestor dependence makes its Jacobian lower
block triangular.
\begin{ass}[Sequential identification and curvature]\label{main:identification}
For each retained edge, put
$Q_e(\vartheta)=E\ell_e(\theta_{<e,0},\vartheta,Z_t)$.
\textup{(i)} The population objective identifies the edge parameter:
\[
\operatorname*{arg\,max}_{\vartheta\in\Theta_e}Q_e(\vartheta)=\{\theta_{e,0}\}.
\]
\textup{(ii)} Its local curvature is nonsingular:
\[
I_e=-E\partial_{\theta_e}g_e(\theta_0,Z_t)\succ0.
\]
\end{ass}
\begin{ass}[Normal-score smoothness]\label{main:regularity}
Each $\ell_e$ is continuous on a neighborhood $\mathcal N$ of the
compact parameter product and has continuous mixed derivatives
satisfying, for finite $C,r$,
\[
\sup_{\theta\in\mathcal N}
\bigl\|\partial_\theta^\alpha\partial_z^\gamma
\ell_e(\theta_{\le e},z)\bigr\|
\le C(1+\|z\|^r),\qquad
|\alpha|\le3,\quad |\gamma|\le2,\quad |\alpha|+|\gamma|\le3.
\]
The constants may depend on the fixed model, but not on $T$ or $z$.
\end{ass}
Identification and curvature serve different purposes. Under the
remaining regularity conditions,
Assumption~\ref{main:identification}(i) identifies the limit of each
stagewise maximizer. Part (ii) makes the local score system invertible
around that limit. Neither global concavity nor uniqueness of every
stationary point is required. This distinction between a population
maximizer and a score root also matters for the numerical implementation
and appears in \citet[Section 4.3]{nagler_kruger_min_2022}.
The smoothness condition applies to the composite edge objective,
including every ancestor transformation. This is essential because an
earlier pair estimate changes the inputs to later pairs. Normal-score
coordinates make the condition tractable: derivatives that are difficult
to bound near a PIT boundary can instead be controlled by polynomials
in $z$. Normal coordinate margins give those polynomials all fixed
moments, even under non-Gaussian joint dependence. The family
verifications below follow the full recursion rather than checking
isolated pair densities.
A separate argument controls the recursive regressions. The stable VAR
filter and growing training prefix give uniform bounds over the retained
fits in Lemmas~\ref{dyn:moments}--\ref{dyn:first-stage}. Finite-lag
approximation of that filter also establishes the joint CLT in
Appendix~\ref{dji:section}, without imposing an additional mixing-rate
condition. These bounds depend on the fixed stationary law. They do
not provide uniform coverage near a unit root or a pair-parameter
boundary. Correct copula specification is needed for a further reason:
a pseudo-true limit, permitted in parts of
\citet{nagler_kruger_min_2022}, need not deliver the conditional event
probability in \eqref{main:conditional-target}.
For $a=(a_\mu',a_v')'$, put
$[T_a(z)]_i=(z_i-a_{\mu,i})/(1+a_{v,i})^{1/2}$. The score moments are
\begin{align}\label{main:score-matrices}
\varphi(z)&=(z_1,\ldots,z_n,z_1^2-1,\ldots,z_n^2-1)',\notag\\
J&=E\partial_\theta g(\theta_0,Z_t),\qquad
S=E\{g(\theta_0,Z_t)g(\theta_0,Z_t)'\},\notag\\
\Sigma_\varphi&=E\{\varphi(Z_t)\varphi(Z_t)'\},\qquad
\Gamma=E\{g(\theta_0,Z_t)\varphi(Z_t)'\},\notag\\
B_g&=E\,\partial_a g\{\theta_0,T_a(Z_t)\}\big|_{a=0}.
\end{align}
The columns of $B_g$ measure score sensitivity to standardized location
and relative variance errors. The diagonal blocks of $J$ are the
own-edge curvatures $-I_e$, while its lower blocks transmit estimation
error from ancestors. Thus the correction has both a marginal stage
and a sequence of dependence stages, consistent with the distinction
between marginal and copula estimation emphasized by
\citet{joe_2005,patton_2012}. Separate pairwise standard errors omit
the latter propagation.
Proposition~\ref{score:gaussian-family} verifies these conditions for
Gaussian vines: conditional normal scores are linear in the original
scores, making the edge objectives quadratic. For Clayton vines,
Online Appendix~\ref{tail:section} uses logarithms of conditional PITs
to propagate polynomial bounds through every ancestor, including
nonindependent higher-tree pairs. Free Clayton parameters stay in
compact subsets of $(0,\infty)$. Independence edges can be fixed in
advance.
The Clayton pair is
\[
C_\vartheta(u,v)=(u^{-\vartheta}+v^{-\vartheta}-1)^{-1/\vartheta},
\qquad \vartheta>0;
\]
its density and conditional transforms are obtained by differentiation
on $(0,1)^2$. Only the CDF is extended continuously to the closed square.
The additional conditions for forecast probabilities are stated in
Section~\ref{main:forecast-section}.
Selecting independence at a parameter boundary is a different problem.
\subsection{The joint dynamic limit}\label{main:dynamic-section}
Write $Q=E(x_tx_t')$ and $e_1=(1,0,\ldots,0)'$. Define the terminal
errors and their influence vector by
\begin{align}\label{main:terminal-influences}
d_{i,T}&=(\widehat\beta_{i,T}-\beta_i)/\sigma_i,\qquad
a^v_{i,T}=\widehat v_{i,T}/\sigma_i^2-1,\notag\\
\delta_T&=(d_{1,T}',\ldots,d_{n,T}',a^v_{1,T},\ldots,a^v_{n,T})',\notag\\
v_t&=(Z_{1t}^2-1,\ldots,Z_{nt}^2-1)',\qquad
w_{i,t}=Q^{-1}x_tZ_{it},\notag\\
e_t&=(w_{1,t}',\ldots,w_{n,t}',v_t')',\qquad
\Delta_T=((\widehat\theta-\theta_0)',\delta_T')'.
\end{align}
Coefficient blocks are stacked by response equation. The joint
covariance has the decomposition
\begin{align}
\xi_t&=\begin{pmatrix}
-J^{-1}\{g(\theta_0,Z_t)+B_g\varphi(Z_t)\}\\ e_t
\end{pmatrix},\notag\\
\mathcal W_{\rm VAR}
&=E(\xi_t\xi_t')+
\begin{pmatrix}J^{-1}B_g\Sigma_\varphi B_g'J^{-\prime}&0\\0&0\end{pmatrix}.
\label{main:covariance-decomposition}
\end{align}
The outer-product term combines marginal estimation with the copula
score and retains their cross covariances. The second term is the
additional variation generated by the dated marginal fits, and it acts
only on the copula block. Section~\ref{main:dynamic-comparisons} will
identify the first term as the terminal-refit covariance in the Gaussian
submodel. For the general joint limit, the decomposition is useful both
for interpretation and for constructing a positive semidefinite estimate.
To display its blocks, let
$R_Z=E(Z_tZ_t')$, $M_3=E(Z_tv_t')$, and $S_v=E(v_tv_t')$.
Partition $M=\Gamma+B_g\Sigma_\varphi=(M_\mu,M_v)$, with $n$ columns
in each part, and define
\begin{align}\label{main:joint-covariance}
W_\delta&=\begin{pmatrix}
R_Z\otimes Q^{-1}&M_3\otimes e_1\\
M_3'\otimes e_1'&S_v
\end{pmatrix},\notag\\
D_{\rm VAR}&=-J^{-1}
(M_{\mu,1}e_1',\ldots,M_{\mu,n}e_1',M_v),\notag\\
\mathcal W_{\rm VAR}&=
\begin{pmatrix}V_{\rm seq}&D_{\rm VAR}\\D_{\rm VAR}'&W_\delta\end{pmatrix},
\end{align}
where $M_{\mu,i}$ denotes column $i$ and
\begin{equation}\label{main:sequential-covariance}
V_{\rm seq}=J^{-1}
\{S+\Gamma B_g'+B_g\Gamma'+2B_g\Sigma_\varphi B_g'\}J^{-\prime}.
\end{equation}
The block $M_3$ links coefficient and scale uncertainty. Normal
coordinate margins alone do not make it zero: asymmetric dependence
can produce nonzero mixed third moments.
\begin{thm}[Joint inference after recursive VAR estimation]\label{main:dynamic}\label{dji:joint}
Under Assumptions~\ref{main:var-assumption}--\ref{main:regularity},
with the estimators and $b_T$ defined in Section~\ref{sec:method}, the sequential estimator is
consistent. With $N_T=T-b_T$ and
$h_{T,s}=\sum_{m=\max(b_T,s)}^{T-1}m^{-1}$, interpreted as zero for an
empty sum,
\begin{align}\label{main:sequential-expansion}
\sqrt T(\widehat\theta-\theta_0)
&=-\frac{J^{-1}}{\sqrt T}\sum_{s=1}^T
\{1\{s>b_T\}g(\theta_0,Z_s)+h_{T,s}B_g\varphi(Z_s)\}+o_p(1),
\\
\sqrt T\delta_T&=T^{-1/2}\sum_{s=1}^Te_s+o_p(1),\qquad
\sqrt T\Delta_T\Longrightarrow N(0,\mathcal W_{\rm VAR}).
\label{main:joint-limit}
\end{align}
The covariance estimator in \eqref{main:joint-feasible} is positive
semidefinite and consistent. The terminal fitted conditional margins
satisfy
\begin{equation}\label{main:dynamic-margin-rate}
\max_{i\le n}\sup_y
|\widehat F_{i,T+1|T}(y)-F^0_{i,T+1}(y)|=O_p(T^{-1/2}).
\end{equation}
\end{thm}
An innovation affects the dependence estimate both through its own
copula score and through the marginal fits used to standardize later
observations. Its contribution to a fit based on the first $m$
observations is proportional to $1/m$. Summing over the retained fits
that reuse innovation $s$ gives the harmonic weight $h_{T,s}$ in
\eqref{main:sequential-expansion}. Earlier innovations enter more of
these fits and receive greater weight. The normalized sum of squared
weights tends to two, whereas the corresponding average over retained
score dates tends to one. These limits explain the coefficients in
\eqref{main:sequential-covariance}. The recursive weighting mechanism
is familiar from \citet{west_1996}. Here it applies to the marginal
innovation moments and their joint covariance with the terminal VAR fit.
The reduction of the historical regression error to $\varphi(Z_s)$
depends on averaging over the lag designs. The intercept identity
$Qe_1=E(x_t)$ implies $E(x_t)'Q^{-1}x_s=1$, so the mean design leaves
only an innovation location term. Stability makes the normalized
contribution of the centered designs negligible, as established in
Lemma~\ref{dyn:interaction}. Location and scale errors therefore remain
in the historical correction, even though every prefix regression
estimates the full coefficient vector.
The terminal forecast is evaluated at one observed lag state, where
there is no such averaging. Errors in both the intercept and the slopes
affect its conditional mean, and their covariance retains the full
block $R_Z\otimes Q^{-1}$. The historical and terminal fits also share
innovations, which accounts for the cross block in
\eqref{main:joint-covariance}. The covariance comparisons in
Section~\ref{sec:sim-components} distinguish these contributions at a
common point forecast. Appendix~\ref{dji:section} gives the complete
argument, with the prefix and moment bounds in Online
Appendix~\ref{supp:joint-proofs}.
For implementation, use the stored normal scores, mathematically
$\widehat Z_{it}=\Phi^{-1}(\widehat U_{it})$, and put
$\widehat g_t=g(\widehat\theta,\widehat Z_t)$, and
$\widehat\varphi_t=\varphi(\widehat Z_t)$. Estimate $J$, $B_g$, and
$\Sigma_\varphi$ by the retained-sample averages of their defining
derivatives and outer products in \eqref{main:score-matrices}.
Let $\widehat e_t$ substitute $\widehat Q_T^{-1}$ and $\widehat Z_t$
in \eqref{main:terminal-influences}. A feasible joint covariance is
\begin{align}\label{main:joint-feasible}
\widehat\xi_t&=
\begin{pmatrix}-\widehat J^{-1}
(\widehat g_t+\widehat B_g\widehat\varphi_t)\\\widehat e_t\end{pmatrix},
\notag\\
\widehat{\mathcal W}_{\rm VAR}
&=\frac1{N_T}\sum_{t>b_T}\widehat\xi_t\widehat\xi_t'
+\begin{pmatrix}
\widehat J^{-1}\widehat B_g\widehat\Sigma_\varphi
\widehat B_g'\widehat J^{-\prime}&0\\0&0
\end{pmatrix}.
\end{align}
Use the identity matrix if a required inverse is absent. The probability
of this fallback tends to zero. Each term in
\eqref{main:joint-feasible} is positive semidefinite, so the construction
avoids a separate adjustment to enforce that property. The common outer
product already contains the copula--terminal cross covariances.
Only the copula block receives the additional historical correction.
\subsection{A confidence interval at the random forecast state}
\label{main:forecast-section}
For a lag state $r\in\mathbb R^{d_x}$, write
$z_i(r)=(c_i-r'\beta_i)/\sigma_i$ and $u_i(r)=\Phi\{z_i(r)\}$.
Let $C_i(\theta,u)=\partial C_\theta(u)/\partial u_i$.
The gradient in the normalized coordinates of $\Delta_T$ has blocks
\begin{align}\label{main:event-gradient}
d_\theta(r)&=\partial_\theta C_{\theta_0}\{u(r)\},\notag\\
d_{\beta_i}(r)&=-C_i(\theta_0,u(r))\phi\{z_i(r)\}\,r,\notag\\
d_{v_i}(r)&=-\tfrac12 C_i(\theta_0,u(r))z_i(r)\phi\{z_i(r)\}.
\end{align}
Stack these columns in the same order as $\Delta_T$ and put
$\tau^2(r)=d(r)'\mathcal W_{\rm VAR}d(r)$.
The coefficient gradient includes the entire state $r$.
Let $X$ have the stationary law of $x_t$.
\begin{ass}[Event smoothness and nondegeneracy]\label{main:event-regularity}
\textup{(i)} The map $(\theta,u)\mapsto C_\theta(u)$ is twice
continuously differentiable on $\mathcal N_0\times(0,1)^n$, for an
open neighborhood $\mathcal N_0$ of $\theta_0$.
\textup{(ii)} For the fixed event under consideration,
$P\{\tau^2(X)>0\}=1$.
\end{ass}
Part (i) permits a Taylor expansion of the event probability. Part (ii)
is needed for studentization and concerns the variance of this event
alone, so the full joint covariance may be singular. No positive lower
bound over all states is required. As in the prediction results of
\citet{nagler_kruger_min_2022}, estimation regularity and smoothness of
the prediction functional are separate requirements. The remaining
issue here is the joint behavior of estimation error and the terminal
state, resolved by Lemma~\ref{dfe:state-separation}.
\begin{thm}[Conditional event inference]\label{main:event-inference}\label{dfe:event-clt}
Under Theorem~\ref{main:dynamic} and
Assumption~\ref{main:event-regularity}(i), let
$G\sim N(0,\mathcal W_{\rm VAR})$ be independent of $X$. Then
\begin{equation}\label{main:event-limit}
\bigl(\sqrt T\{\widehat p_{T+1}(c)-p_{T+1}(c)\},x_{T+1}\bigr)
\Longrightarrow (d(X)'G,X).
\end{equation}
Let $\widehat d_T$ substitute the fitted parameters and terminal
standardized thresholds in \eqref{main:event-gradient}, and define
$\widehat\tau_T^2=\widehat d_T(x_{T+1})'
\widehat{\mathcal W}_{\rm VAR}\widehat d_T(x_{T+1})$.
Then $\widehat\tau_T^2-\tau^2(x_{T+1})\to_p0$.
If Assumption~\ref{main:event-regularity}(ii) also holds, then for every
fixed $0<\alpha<1$,
\begin{align}
\frac{\sqrt T\{\widehat p_{T+1}(c)-p_{T+1}(c)\}}{\widehat\tau_T}
&\Longrightarrow N(0,1),\label{main:pivot}\\
P\left\{p_{T+1}(c)\in
[\widehat p_{T+1}(c)\pm z_{1-\alpha/2}\widehat\tau_T/\sqrt T]
\cap[0,1]\right\}&\longrightarrow1-\alpha.
\label{main:interval}
\end{align}
Set the statistic to zero when its estimated variance is zero. That
event has vanishing probability under the positivity condition.
\end{thm}
The estimates and the terminal state use the same observations, but
they depend on different parts of the innovation history in the limit.
Estimation error accumulates over the full sample, and a slowly growing
block of the most recent innovations contributes negligibly to its
normalized influence. In a stable VAR, the terminal state can instead
be approximated by those recent innovations, since the effect of the
distant past decays geometrically. Removing that block from the
estimation influence and truncating the state to it leaves two
independent objects with asymptotically negligible approximation
errors. This separation explains the independence of $G$ and $X$ in
\eqref{main:event-limit}. It does not require independence of the
estimates and the forecast state in a finite sample.
For a given limiting state $X$, the leading forecast error is the
Gaussian projection $d(X)'G$, with variance $\tau^2(X)$. Different
states change the sensitivity of the forecast to parameter error,
producing a Gaussian mixture across samples. The feasible standard
error adapts to the realized state, so studentization removes this
random scale and yields \eqref{main:pivot}. Almost-sure positivity of
$\tau^2(X)$ suffices. States with arbitrarily small positive variance
can be handled by localization, without imposing a uniform lower bound.
Appendix~\ref{dfe:section} establishes the separation and studentization
arguments. The resulting coverage statement concerns a conditional
probability evaluated at a random realized state. It does not guarantee
coverage conditional on every training history or provide a prediction
set for the next binary outcome. Section~\ref{sec:sim-coverage} examines
this repeated-sample coverage in finite samples.
For the family verifications and rectangle extension, a sufficient
smoothness condition on the full log density
$\ell^*(\theta,z)=\sum_e\ell_e(\theta_{\le e},z)$ is
\begin{equation}\label{main:event-envelope}
\sup_{\theta\in\mathcal N_0}
\|\partial_\theta^\alpha\partial_z^\gamma\ell^*(\theta,z)\|
\le C(1+\|z\|^r),\qquad |\alpha|+|\gamma|\le3,
\end{equation}
with continuous derivatives on a neighborhood of a compact parameter
set containing $\theta_0$ in its interior. This includes pure third
derivatives in $z$, which Assumption~\ref{main:regularity} does not
require. Lemma~\ref{pred:cdf-smooth} shows that it implies
Assumption~\ref{main:event-regularity}(i). Gaussian log densities are
quadratic in $z$. Lemma~\ref{tail:propagation} verifies the bound for
the positive Clayton vines, including the ancestor transformations.
For a fixed rectangle $A=\prod_i(\ell_i,u_i]$, replace the CDF by
its probability over $A$ and use the corresponding gradient.
Corollary~\ref{dfe:rectangle-clt} proves the same interval, treating
infinite endpoints as inactive boundaries, under
\eqref{main:event-envelope} and the event-specific condition
$P\{\tau_A^2(X)>0\}=1$. This covers fixed mixed directions and target subsets without
deleting factors from a larger fitted vine. Thresholds estimated from
the same observations or moving into the tail require additional
influence or uniformity arguments.
\subsection{Dynamic refitting and finite training weights}\label{main:dynamic-comparisons}
The Gaussian VAR submodel provides a feasible benchmark for the archive
rule. Keep the observations and terminal marginal estimates fixed,
recompute every standardized residual using those terminal estimates,
and fit the same sequential Gaussian vine:
\begin{equation}
\widetilde Z_{it}=\frac{y_{it}-x_t'\widehat\beta_{i,T}}
{\sqrt{\widehat v_{i,T}}},\quad 1\le t\le T,
\qquad
\widetilde\theta_e\in\arg\max_{\theta_e\in\Theta_e}
\frac1T\sum_{t=1}^T
\ell_e(\widetilde\theta_{<e},\theta_e,\widetilde Z_t).
\label{dfs:estimator}
\end{equation}
All $T$ residuals enter this comparator, with the same parameter spaces
and measurable tie-breaking as in \eqref{eq:sequential-fit}.
Proposition~\ref{dfs:theorem} proves the joint influence
vector with copula block $-J^{-1}\{g_t+B_g\varphi_t\}$ and terminal
block $e_t$, together with the covariance ordering
\begin{equation}\label{main:dynamic-archive-cost}
\mathcal W_{\rm VAR}-\mathcal W_{\rm fs}
=\begin{pmatrix}J^{-1}B_g\Sigma_\varphi B_g'J^{-\prime}&0\\0&0\end{pmatrix}
\succeq0.
\end{equation}
Evaluate $\widehat g_t,\widehat\varphi_t,\widehat J,\widehat B_g$
at $(\widetilde\theta,\widetilde Z_t)$ using all $T$ observations,
and form $\widehat e_t$ from the same residuals and $\widehat Q_T^{-1}$.
The feasible full-sample covariance is
\begin{equation}
\widehat\xi_t^{\rm fs}=
\begin{pmatrix}-\widehat J^{-1}
(\widehat g_t+\widehat B_g\widehat\varphi_t)\\\widehat e_t\end{pmatrix},
\qquad
\widehat{\mathcal W}_{\rm fs}
=T^{-1}\sum_{t=1}^T\widehat\xi_t^{\rm fs}
\widehat\xi_t^{{\rm fs}\prime}.
\label{dfs:feasible}
\end{equation}
Use the identity fallback when a required inverse is absent.
The refit covariance retains every terminal slope and scale. Its
derivation combines the exact OLS residual cross-product identity with
the quadratic Gaussian edge scores, accounting for dependence between
the lag design and past innovations. Equation~\eqref{main:dynamic-archive-cost}
then compares limiting covariances at the same population parameters.
The extra archive term measures the variance cost of training the copula
on successive out-of-sample marginal forecasts while updating the
terminal margins.
The ordering need not extend to finite-sample standard errors evaluated at separately
fitted parameters. Nor does it compare estimators trained on different
data vintages or forecasting laws.
The archive's initial training period can occupy a substantial fraction
of a practical sample even though $b_T/T\to0$ asymptotically. We can
retain that fraction in the covariance of the leading array. Write
$b=b_T$, $N=T-b$, and
$H_T=\sum_{m=b}^{T-1}m^{-1}$. Retaining the exact $T/N$ normalization
in the leading score array replaces the copula-score covariance by
\begin{align}\label{main:finite-weights}
\Omega_T^{\rm fin}
&=k_{11}S+k_{12}(\Gamma B_g'+B_g\Gamma')+
k_{22}B_g\Sigma_\varphi B_g',\notag\\
k_{11}&=T/N,\qquad k_{12}=T(N-bH_T)/N^2,\qquad
k_{22}=T\{2N-(2b-1)H_T\}/N^2.
\end{align}
The copula--terminal cross coefficients remain one and the terminal
block is unchanged. To implement this covariance with one common
empirical moment matrix, set $a_{T,s}=(T/N)\mathbf1\{s>b\}$ and
$d_{T,s}=(T/N)h_{T,s}$ and form
\begin{align}
\widehat\zeta_t&=(\widehat g_t',\widehat\varphi_t',\widehat e_t')',&
\widehat C_T&=\frac1N\sum_{t>b}\widehat\zeta_t\widehat\zeta_t',\notag\\
\widehat L_{T,s}
&=\begin{pmatrix}
-a_{T,s}\widehat J^{-1}&
-d_{T,s}\widehat J^{-1}\widehat B_g&0\\
0&0&I_{n d_x+n}
\end{pmatrix},&
\widehat{\mathcal W}_T^{\rm fin}
&=\frac1T\sum_{s=1}^T
\widehat L_{T,s}\widehat C_T\widehat L_{T,s}'.
\label{fdw:feasible}
\end{align}
The score, nuisance and terminal quantities use the same retained
observations as \eqref{main:joint-feasible}. Use its identity fallback.
In computation, partition $\widehat C_T$ into its score, nuisance and
terminal blocks and apply the coefficients in \eqref{main:finite-weights}.
The displayed sum need not be evaluated date by date.
Proposition~\ref{fdw:proposition} proves that this estimator is positive
semidefinite and consistent and can replace
\eqref{main:joint-feasible} in the event interval. At $T=340$, the
simulation rule $b_T=\lceil T^{0.7}\rceil$ gives $b=60$ and coefficients
$1.2143$, $0.7611$, and $1.5298$, instead of $1,1,2$.
These coefficients are exact for the leading array. The estimator's
linearization remainder is still omitted, so the calculation supplies
neither an exact finite-sample variance nor a coverage improvement
theorem. Section~\ref{sec:sim-small} compares intervals based on these
weights with limiting-weight and terminal-refit intervals.
The terminal slopes also explain why one-sided event intervals are
nondegenerate. Centering the lag design isolates a coefficient influence
orthogonal to every function of the current innovation, including the
copula scores and scale influences. Its projection on the event gradient
is positive except at the mean lag state, which has probability zero.
Proposition~\ref{ond:proposition} therefore verifies the required
almost-sure positivity for Gaussian and positive Clayton archives.
Gaussian sign symmetry gives the stronger conclusion of positivity at
every finite state for both archive and refit
(Proposition~\ref{gnd:proposition}). General bounded rectangles retain
the event-specific condition.
\subsection{Probability boundaries and the cost of the archive}
Near zero or one, a symmetric probability-scale interval can have a
different finite-sample behavior from an interval constructed after
transformation. Applying the normal approximation to log odds gives
asymmetric probability endpoints and preserves the first-order coverage
result. Write $\widehat p_T$ for the fitted event probability
at $x_{T+1}$, and put
$\ell(p)=\log\{p/(1-p)\}$, $\Lambda(u)=(1+e^{-u})^{-1}$, and
$\widehat s_T=\widehat\tau_T/\sqrt T$. For an interior fitted probability,
the logit interval is
\begin{equation}\label{main:logit-interval}
I_{\ell,T}=\left[
\Lambda\left\{\ell(\widehat p_T)-
\frac{z_{1-\alpha/2}\widehat s_T}{\widehat p_T(1-\widehat p_T)}\right\},\quad
\Lambda\left\{\ell(\widehat p_T)+
\frac{z_{1-\alpha/2}\widehat s_T}{\widehat p_T(1-\widehat p_T)}\right\}
\right].
\end{equation}
Proposition~\ref{rli:logit} proves its repeated-sample coverage under the
corresponding event theorem, $0<p_A(X)<1$ and $\tau_A^2(X)>0$ almost
surely. Its proof localizes away from probability boundaries and small
projected variances, so no global lower bound over all lag states is
imposed. A deterministic vanishing clip defines the formula at numerical endpoints. Section~\ref{sec:sim-logit} examines finite-sample coverage.
The size of the archive correction depends on the event as well as the
forecast state. Partition the event gradient
into dependence and marginal components, and set
\[
q_A(r)=\|\Sigma_\varphi^{1/2}B_g'J^{-\prime}d_{A,\theta}(r)\|^2.
\]
Proposition~\ref{rli:cost} gives
$\tau_{\rm ar,A}^2(r)=\tau_{\rm fs,A}^2(r)+q_A(r)$ in the Gaussian
dynamic comparison. Thus $q_A(r)/\tau_{\rm ar,A}^2(r)$ measures the
share of event variance due specifically to preserving issued margins.
The same proposition quantifies limiting undercoverage when a refit
variance is used at the archived probability estimate. Holding that
estimate fixed also permits an exact calculation of when omitting the
archive term changes exclusion of a reference probability. These
comparisons isolate a covariance omission from changes in the fitted
forecast. The application reports both the archive share and the effect
of removing terminal-margin uncertainty.
\subsection{Numerical implementation}\label{main:implementation}
Implementation follows the same sequence as the theory. Fit the vine
to the dated transforms, then use their stored normal scores and lag
rows to form \eqref{main:joint-feasible}. Evaluate the event probability
and its gradient at $x_{T+1}$ using the terminal marginal estimates.
The resulting standard error gives either
\eqref{main:interval} or \eqref{main:logit-interval}.
\citet{stoeber_schepsmeier_2013} give algorithms for analytic scores
and observed information in regular-vine models. Here, the covariance
calculation also requires the derivatives with respect to the dated
marginal inputs and the terminal forecast parameters. Analytical and
numerical derivatives are both admissible under the accuracy conditions
below.
Propositions~\ref{op:numerical-score}--\ref{op:numerical-event} justify
a consistent local numerical solution whose average score residual is
$o_p(T^{-1/2})$, probability error is $o_p(T^{-1/2})$, and gradient
and covariance errors are $o_p(1)$. A small residual alone does not
select the correct stationary branch. Consistency or a vanishing
stagewise objective gap is also required. The bound in \eqref{op:finite-difference-bound} also justifies a forward
difference with step $s_T=T^{-1/4}$ when probability-evaluation error
at the required points is $o_p(T^{-1/2})$.
For event and gradient evaluation, these bounds hold uniformly on
each compact set of lag states and a neighborhood of the true
parameters. Stationarity then permits evaluation at the random
terminal state.
\section{Monte Carlo Simulations}\label{sec:simulation}
The simulations examine how the coverage of conditional joint-event
probability intervals changes with sample size and persistence.
We compare Wald and logit intervals over a common sample-size grid
and assess the role of marginal-estimation uncertainty. The coverage
target is the true conditional probability at the realized terminal
state of each simulated sample.
\subsection{Simulation design}\label{sec:sim-design}\label{replication:sim_dgp}
We generate observations from the VAR(1)
\[
Y_t=AY_{t-1}+Z_t,\qquad
A=0.35I_n+0.10P_n\quad\text{or}\quad A=0.70I_n+0.15P_n,
\]
where $P_n$ is the cyclic shift matrix. The two coefficient matrices
have spectral radii $0.45$ and $0.85$, respectively, and define the
baseline and persistent designs. Innovation vectors are independent
over time, with standard-normal coordinate margins and either a
Gaussian or a positive Clayton D-vine. This comparison allows for
symmetric dependence and for lower-tail dependence, which is relevant
to simultaneous adverse outcomes
\citep{joe_li_nikoloulopoulos_2010,nikoloulopoulos_joe_li_2012}.
The principal design has three variables. The pair parameters, in
tree order $(12,23,13\mid2)$, are $(0.4,0.35,0.2)$ for Gaussian pairs
and $(1,0.75,0.5)$ for Clayton pairs. The second tree is nonindependent
in both cases, so estimation error from the first tree enters the
last pair. We use $T\in\{200,250,340,500,800\}$ to examine shorter
estimation histories. This range brackets the 327--776 observations
available in the application. Samples of $T=2000$ and $8000$ provide
larger-sample benchmarks for the approximation under the same
stationary law. Additional
bivariate baseline results, with pair parameters $0.5$ and $1$ for
the Gaussian and Clayton families, are reported in the supplement.
Every specification uses 1,000 independent Monte Carlo replications.
We initialize the process at zero and discard 1,000 observations
before collecting the estimation sample.
Estimation follows Section~\ref{sec:method}. Although the true
intercepts are zero and innovation scales are one, all intercepts,
lag coefficients, and scales are estimated. Prefix estimates produce
the dated transforms, and the copula fit uses the transforms after
$b_T=\lceil T^{0.7}\rceil$. The terminal marginal fit uses all $T$
observations. We retain the Gram-matrix and variance safeguards in
every replication. The vine structure and pair families are correctly
specified and held fixed, as required by the maintained theory.
Four events allow the comparison to cover different parts of the
conditional distribution. Zero constrains every coordinate of
$Y_{T+1}$ to be at most zero, and Tail constrains every coordinate
to be at most $-1$. Mixed requires $Y_{1,T+1}\le0$ and
$Y_{2,T+1}>0$, leaving any other coordinate unrestricted.
Rectangle requires $-1<Y_{i,T+1}\le0.75$ for all $i$. In each
replication, deterministic integration evaluates the true and fitted
event probabilities at the same realized lag vector $Y_T$. The true
parameters enter only the calculation of the coverage target.
We construct nominal 95\% Wald and logit intervals. For a given
estimator and covariance estimate, the two intervals use the same
fitted probability and standard error. The archive estimator uses
either the limiting joint covariance in \eqref{main:joint-feasible}
or the exact finite training weights in
Proposition~\ref{fdw:proposition}. For Gaussian vines, we also use
the feasible terminal refit of Section~\ref{main:dynamic-comparisons},
which recomputes historical transforms with the current marginal
estimates and re-estimates dependence. All comparisons within a
specification use the same simulated samples. To isolate the role of
individual uncertainty components, we additionally construct intervals
that hold the archive point estimate fixed and omit selected terms
from its variance calculation. These omissions are diagnostics of the
joint covariance formula.
Coverage is the proportion of all 1,000 replications in which an
interval contains the true conditional event probability. Unavailable
intervals count as noncoverage, and mean interval length uses the
available intervals. A coverage estimate of 95\% has Monte Carlo
standard error 0.69 percentage points. Uncertainty for differences
between methods is calculated from the paired replication outcomes.
Online Appendix~\ref{app:simulation-details} gives the numerical
settings, replication seeds, availability counts, and supporting results.
\subsection{Simulation results}
\label{sec:sim-coverage}\label{sec:sim-small}\label{sec:sim-logit}
Table~\ref{tab:simulation-coverage} reports lower-tail coverage for
the seven sample sizes under baseline and persistent dynamics.
Wald coverage generally moves closer to the nominal 95\%
level as the sample grows, but substantial shortfalls remain at the
sample sizes relevant to the application. For the persistent Gaussian
design, coverage is 85.8\% at $T=200$, 90.9\% at $T=800$, and
91.8\% at $T=2000$. It reaches 95.1\% at $T=8000$.
\begin{table}[!t]
\centering\small
\caption{Lower-Tail Interval Coverage by Sample Size}
\label{tab:simulation-coverage}\label{tab:small-sample-coverage}
\setlength{\tabcolsep}{5pt}
\begin{tabularx}{\linewidth}{@{\extracolsep{\fill}}lcccc}
\toprule
& \multicolumn{2}{c}{Baseline dynamics} & \multicolumn{2}{c}{Persistent dynamics} \\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}
$T$ & Wald & Logit & Wald & Logit \\
\midrule
\multicolumn{5}{@{\extracolsep{\fill}}l}{\textit{Panel A: Gaussian vine}} \\
\addlinespace[2pt]
200 & $87.1\,(1.06)$ & $94.1\,(0.75)$ & $85.8\,(1.10)$ & $95.0\,(0.69)$ \\
250 & $87.3\,(1.05)$ & $96.1\,(0.61)$ & $86.4\,(1.08)$ & $94.1\,(0.75)$ \\
340 & $89.0\,(0.99)$ & $95.3\,(0.67)$ & $88.3\,(1.02)$ & $94.9\,(0.70)$ \\
500 & $90.6\,(0.92)$ & $94.7\,(0.71)$ & $87.5\,(1.05)$ & $93.8\,(0.76)$ \\
800 & $92.6\,(0.83)$ & $94.4\,(0.73)$ & $90.9\,(0.91)$ & $94.4\,(0.73)$ \\
2000 & $93.4\,(0.79)$ & $94.7\,(0.71)$ & $91.8\,(0.87)$ & $93.6\,(0.77)$ \\
8000 & $95.4\,(0.66)$ & $95.7\,(0.64)$ & $95.1\,(0.68)$ & $94.9\,(0.70)$ \\
\addlinespace[4pt]
\multicolumn{5}{@{\extracolsep{\fill}}l}{\textit{Panel B: Clayton vine}} \\
\addlinespace[2pt]
200 & $86.4\,(1.08)$ & $93.5\,(0.78)$ & $86.8\,(1.07)$ & $94.5\,(0.72)$ \\
250 & $89.5\,(0.97)$ & $93.9\,(0.76)$ & $88.4\,(1.01)$ & $95.2\,(0.68)$ \\
340 & $89.8\,(0.96)$ & $95.0\,(0.69)$ & $88.6\,(1.01)$ & $93.3\,(0.79)$ \\
500 & $91.9\,(0.86)$ & $94.8\,(0.70)$ & $91.9\,(0.86)$ & $94.8\,(0.70)$ \\
800 & $93.0\,(0.81)$ & $94.6\,(0.71)$ & $92.7\,(0.82)$ & $94.7\,(0.71)$ \\
2000 & $93.4\,(0.79)$ & $94.4\,(0.73)$ & $94.1\,(0.75)$ & $94.3\,(0.73)$ \\
8000 & $95.2\,(0.68)$ & $95.9\,(0.63)$ & $94.1\,(0.75)$ & $94.7\,(0.71)$ \\
\bottomrule
\end{tabularx}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize \textit{Note:}
Entries are coverage percentages for nominal 95\% intervals for
$\Pr(Y_{i,T+1}\le-1:i=1,2,3\mid Y_T)$, using the archive estimator
and limiting joint covariance. Parentheses give Monte Carlo standard errors
in percentage points. Baseline and persistent dynamics have spectral radii
$0.45$ and $0.85$. Each specification uses 1,000 replications, with Wald and
logit intervals computed from the same fitted probabilities and standard
errors. Unavailable intervals count as noncoverage. Samples of
$T=2000$ and $8000$ assess the large-sample approximation beyond the
327--776 observations in the application. Online
Appendix~\ref{app:simulation-details} reports the other events,
interval lengths, finite-weight and refit comparisons, and simulation details.
\end{minipage}
\end{table}
The $T=8000$ cases serve as benchmarks for the large-sample
approximation. They are far longer than the 327--776 observations
available in the real-time application. At this sample size,
lower-tail Wald coverage ranges from 94.1\% to 95.4\%, with Monte
Carlo standard errors of 0.66--0.75 percentage points. These results
support the asymptotic approximation in the designs considered when
the estimation history is sufficiently long. They do not establish
a sample-size threshold for accurate coverage, or imply comparable
accuracy at the shorter histories used in practice. Coverage need
not improve at every step of the grid. For example, persistent
Clayton Wald coverage is 94.1\% at both $T=2000$ and $T=8000$.
The logit interval improves lower-tail coverage at much shorter
histories. Over $T=200$--$800$, its coverage ranges from 93.3\%
to 96.1\%, compared with 85.8--93.0\% for the Wald interval.
The paired gains are 1.6--9.2 percentage points, with Monte Carlo
standard errors of 0.51--1.00 points. Each pair uses the same fitted
probability and standard error, so the difference comes from the
interval transformation. Coverage gains are accompanied by longer
intervals. The largest increase in mean length is 0.0032, in the
baseline Gaussian design at $T=200$, where mean length rises from
0.0402 to 0.0434. Online Table~\ref{tab:sim-short-comparators} reports
the finite-weight and Gaussian terminal-refit comparisons, which
give similar logit coverage over these shorter samples.
The remaining shortfalls are relevant to the use of the method.
Persistent Clayton logit coverage is 93.3\% at $T=340$, and
persistent Gaussian logit coverage is 93.6\% at $T=2000$.
The transformation also slightly lowers coverage for some other
events, as reported in Online Table~\ref{tab:small-other-events}.
The simulations therefore provide support for the limiting
approximation and for logit intervals in the reported lower-tail
designs, with finite-sample accuracy depending on the sample size,
persistence, and event. All safeguard activations are retained.
Both baseline Clayton replications with boundary fits at $T=200$
remain in the coverage denominator, with unavailable intervals
counted as noncoverage.
\label{sec:sim-components}
Online Table~\ref{sim:new:ablations} examines the contribution of
individual covariance terms. This diagnostic uses the three-variable
baseline designs at $T=8000$, holding the fitted probability fixed.
Removing all terminal-margin gradients lowers Wald coverage by
33.6--48.6 percentage points, with paired Monte Carlo standard errors
of 1.49--1.58 points. Removing only terminal lag-coefficient
gradients lowers coverage by 6.7--17.3 points, whereas removing the
additional historical term changes coverage by 0.1--0.9 points.
The importance of terminal slope uncertainty is consistent with
its role in the forecast gradient despite the cancellation of
centered historical lag terms. These comparisons explain the
components of the joint variance calculation at a long history.
The sample-size comparison in Table~\ref{tab:simulation-coverage}
shows where that calculation, together with the interval
approximation, provides accurate coverage in finite samples.
Online Table~\ref{sim:new:dynamic} reports Wald coverage for
the other events and bivariate designs. All simulations
maintain constant coefficients, constant scales, normal innovation
margins, and a correctly specified vine. The coverage evidence
therefore concerns that maintained model, leaving the marginal
misspecification seen in the application to a separate study.
\section{Real-time forecasting of joint contractions}\label{sec:evaluation}
We apply the procedure to forecasts of joint contractions in U.S.
employment and economic activity. At each historical origin, the rule
assigns probabilities to events in the next common release using the
dated observations available to it. This setting makes the distinction
between parameter updating and data revision consequential. Keeping the
recorded history fixed isolates the uncertainty from current marginal
estimation and the retained transforms. Refitting on revised histories
also changes the observations used to estimate the model. We compare
these forecasting rules and examine whether the stationary marginal
specification adequately describes their issued errors. All reported
forecasts are retrospective reconstructions from the information set
defined below.
\subsection{The release target and information set}
Real-time evaluation must specify both what the forecaster observes
and which release supplies the outcome
\citep{croushore_stark_2001,koenig_dolmas_piger_2003}. We use payroll
employment (\texttt{PAYEMS}), industrial production (\texttt{INDPRO}),
and unemployment (\texttt{UNRATE}), following the macroeconomic
application of \citet{tsionas_izzeldin_trapani_2022}. Together they
allow a labor contraction to be compared with a broader decline in
activity. Historical vintages are obtained from the ALFRED (Archival FRED)
database.\footnote{Federal Reserve Bank of St. Louis, ALFRED:
\url{https://alfred.stlouisfed.org/}.} These dated snapshots determine the
observations available at each origin and the release used to evaluate
the forecast. An event defined from that release can differ from the
event recorded in subsequently revised data. This distinction is
central to density forecasting under data uncertainty
\citep{clements_galvao_2023}, while \citet{clark_2011} emphasizes the
role of stochastic volatility in real-time density forecasts.
Let $\tau_m$ be the first ALFRED vintage date at which all three levels
for reference month $m$, and their preceding-month levels, are
available. Each change uses both levels as available at $\tau_m$:
\begin{align*}
E_m=100\log\frac{\mathrm{PAYEMS}_m^{[\tau_m]}}
{\mathrm{PAYEMS}_{m-1}^{[\tau_m]}},\qquad
I_m=100\log\frac{\mathrm{INDPRO}_m^{[\tau_m]}}
{\mathrm{INDPRO}_{m-1}^{[\tau_m]}},\\
U_m=\mathrm{UNRATE}_m^{[\tau_m]}-
\mathrm{UNRATE}_{m-1}^{[\tau_m]}.
\end{align*}
Lock the vector $Y_m^{\rm L}=(10E_m,I_m,-10U_m)'$ at this date.
It records the first \emph{common release}: a component published earlier
may already have been revised when the last one appears. Subsequent
regressions use these locked changes and lag rows. The factors of ten
put the scale-sensitive safeguards on useful numerical scales without
changing the event indicators.
The usable record begins in March 1960, after the initial unemployment
backfill, and ends in December 2024. Its 778 common-release dates are
strictly ordered. The first is April 15, 1960. The last is January 17,
2025. At $\tau_m$, the forecasting rule estimates margins using locked records
through $m$ and forms a forecast for $Y_{m+1}^{\rm L}$. Month $m+1$ has usually begun by then, giving the next-release forecast
a nowcasting interpretation. The forecasting rule uses the locked
history, a smaller information set than all public news at $\tau_m$.
The stationary theory applies only if this recorded sequence satisfies
its maintained law. Availability of historical vintages does not
establish that condition.
The primary evaluation covers July 1987--December 2019, giving 390
forecasts. January 2020--December 2024 supplies a separate 60-month
stress period. The three events, fixed before scoring, are
\begin{align*}
A_L&=\{E<0,\ U>0\},&
A_3&=\{E<0,\ I<0,\ U>0\},&
A_{EP}&=\{E<0,\ I<0\}.
\end{align*}
The recorded indicators use strict changes. The companion event
$A_{EP}$ avoids an outcome boundary on rounded unemployment, while
retaining unemployment as a predictor. It was fixed before scoring.
The continuous working model assigns no mass to an exact boundary.
The recorded unemployment changes do have ties. We therefore score
the observed binary events and examine a bracket for the corresponding
unrounded release event in the supplement. Revision sensitivity uses
September 30, 2026 as a separate outcome vintage.
\subsection{Forecast construction and predictive performance}
The main specification is an expanding VAR(1) with an intercept,
Gaussian coordinate margins, and a fixed D-vine ordered as employment,
production, and negative unemployment change. Each common release is
evaluated against the marginal forecast reconstructed
from the preceding locked history. Its PIT records how surprising that
release was before it entered the regression. We retain these PITs to
train the dependence model on simultaneous marginal forecast surprises,
and include the new release when updating the margins for the next
forecast. This choice is separate from locking the observations:
even with the same release values, recomputing historical PITs with
current estimates would change their reference distributions.
The first and last estimation samples contain 327 and 776 regression
observations.
The dependence fit uses the retained transforms after
$b_T=\lceil T^{0.7}\rceil$, and the terminal marginal fit uses all $T$
observations. Gaussian and positive Clayton vines are both estimated.
The Gaussian first-tree and independence restrictions use the same
margins. VAR(2) and unemployment boundaries of $0.05$ and $0.10$
percentage points are reported in the supplement.
The \emph{locked-history refit} recomputes all residuals with the
terminal marginal fit on the locked observations. It implements the
Gaussian same-history comparator, including its joint covariance
correction. The \emph{latest-vintage refit} instead re-estimates the
VAR and Gaussian dependence on the historical vintage available at
$\tau_m$, including the revised terminal lag. Both predict the next
locked release. Only the first holds the estimation data fixed, so
the covariance ordering applies to it alone.
Predictive performance is assessed by Brier loss, a proper scoring rule
for binary event probabilities \citep{gneiting_raftery_2007}, and by
aggregate calibration \citep{gneiting_balabdaoui_raftery_2007}.
Both assessments use the recorded event indicators. The replication
output also retains continuous joint log densities as diagnostics of
the working model. Those densities are not observation probabilities
for rounded unemployment, for which discrete predictive assessment
requires a different treatment \citep{czado_gneiting_held_2009}.
\begin{table}[htbp]
\centering\small
\caption{Next-Release Event Forecasts: Brier Loss}
\label{emp:alfred:scores}
\begin{tabularx}{\linewidth}{@{\extracolsep{\fill}}lrrrrrr}
\toprule
& \multicolumn{3}{c}{1987:07--2019:12} & \multicolumn{3}{c}{2020:01--2024:12} \\ \cmidrule{2-4}\cmidrule{5-7}
Model & Labor & Three & E--P & Labor & Three & E--P \\
\midrule
Archived Gaussian & 0.0882 & 0.0541 & 0.0819 & 0.1448 & 0.0833 & 0.0995 \\
Locked-history refit & 0.0872 & 0.0526 & 0.0807 & 0.1783 & 0.1146 & 0.1266 \\
Archived Clayton & 0.0871 & 0.0523 & 0.0806 & 0.1011 & 0.0610 & 0.0877 \\
Gaussian first tree & 0.0854 & 0.0526 & 0.0819 & 0.1040 & 0.0714 & 0.0995 \\
Independence & 0.0834 & 0.0476 & 0.0747 & 0.0996 & 0.0606 & 0.0869 \\
Latest-vintage refit & 0.0832 & 0.0509 & 0.0777 & 0.1761 & 0.1117 & 0.1193 \\
\bottomrule
\end{tabularx}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize \textit{Note:}
Smaller loss is better. All models forecast identical locked outcomes on 390 primary and 60 stress dates. E--P denotes employment--production contraction. Finite-weight and limiting archive procedures have identical point forecasts. These are descriptive loss comparisons. No significance ranking is asserted.
\end{minipage}
\end{table}
Table~\ref{emp:alfred:scores} gives little support for a predictive
advantage of the archive in the primary period.
For the primary labor event its Brier loss is $0.0882$, compared with
$0.0872$ for the locked-history refit, $0.0834$ for independence, and
$0.0832$ for the latest-vintage refit. For employment--production
contraction, independence again has lower loss: $0.0747$ against
$0.0819$ for the archive and $0.0777$ for the latest-vintage refit.
Clayton dependence lowers loss relative to the archived Gaussian fit
for all three primary events, yet independence has lower loss than
either archive. The forecasting exercise thus shows how to assess
estimation uncertainty and diagnose limitations of the marginal model.
\subsection{Sources of estimation uncertainty}
The variance decomposition measures parameter uncertainty within the
fitted continuous law. Its interpretation must account for two visible
departures from that law. Recorded unemployment changes have ties, and
the issued marginal errors question the assumed scales and serial
independence. Before 2020, standardized employment errors have standard
deviation $0.563$, while production errors have lag-one correlation
$-0.304$ and lag-one correlation of squared errors $0.404$
(Table~\ref{emp:alfred:diagnostics}). We report model-based standard errors and intervals as measures of sensitivity to parameter estimation under the maintained law. Their coverage under the observed release process is not established.
For the labor event, the median archived-Gaussian standard error is
$0.0178$, and the median logit interval width is $0.0700$. Retaining
finite harmonic weights changes these to $0.0180$ and $0.0708$.
The locked-history refit has median standard error $0.0191$.
The refit's larger median standard error is compatible with the
population covariance ordering. In finite samples, the two procedures
evaluate their covariance matrices and event gradients at different
parameter estimates.
\begin{table}[htbp]
\centering\small
\caption{Estimation Uncertainty about Release Probabilities}
\label{emp:alfred:intervals-main}
\begin{tabularx}{\linewidth}{@{\extracolsep{\fill}}lcc}
\toprule
Event and covariance & Median SE & Median logit width \\
\midrule
Labor: Archived Gaussian & 0.0178 & 0.0700 \\
Labor: Finite weights & 0.0180 & 0.0708 \\
Labor: Locked-history refit & 0.0191 & 0.0749 \\
Employment--production: Archived Gaussian & 0.0177 & 0.0695 \\
Employment--production: Finite weights & 0.0178 & 0.0697 \\
Employment--production: Locked-history refit & 0.0189 & 0.0740 \\
\bottomrule
\end{tabularx}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize \textit{Note:}
Primary period, VAR(1), nominal 95\% intervals. Every interval is numerically available. The intervals measure parameter uncertainty under the maintained model. Their coverage under the observed release process is not established. The supplement reports Wald widths and the remaining events and models.
\end{minipage}
\end{table}
The extra archive term accounts for a median $1.86\%$ of the labor
event variance and $2.89\%$ of the employment--production event
variance. Dropping all terminal-margin gradients, by contrast,
gives median ratios of the reduced standard errors to the complete
values of $37.7\%$ and $32.1\%$. Dropping only terminal slopes gives
median ratios of $93.0\%$ and $87.2\%$. Current marginal estimation and its covariance with the dependence
fit account for most of the adjustment in these real-time forecasts.
Holding the fitted probability and all common covariance blocks fixed
shows how much the archive term changes interval conclusions. Removing
that term changes whether a 95\% symmetric labor-event interval excludes
the reference probability $0.20$ at one of 390 primary dates. The
corresponding count is two dates for employment--production contraction.
These fixed reference probabilities illustrate sensitivity of interval
exclusion, without estimating a policy threshold. The supplement also
reports $0.10$, $0.30$, and $0.50$. The small counts are consistent with
the modest archive variance shares.
\subsection{Calibration, revisions, and the stress period}
The primary period contains 47 labor events, with frequency $0.1205$,
against an average archived Gaussian forecast of $0.1996$. The
employment--production event has frequency $0.1179$ and average
forecast $0.2053$. A fixed-period martingale bound assesses these
aggregate discrepancies without assuming independent months or a
correctly specified VAR.
For an event indicator $I_t$, let
$q_t=E(I_t\mid\mathcal H_{t-1})$, where $\mathcal H_t$ includes the
locked release records, public vintage snapshots available by the issue
date, and predetermined forecasting rules. Issued forecasts $p_{m,t}$
are predictable for the next release. Over a declared period of $n$
dates, the target is the realized-path average
$\overline q_n=n^{-1}\sum_tq_t$. The average calibration gap is
$\overline q_n-\overline p_{m,n}$. Proposition~\ref{cal:proposition}{}
gives simultaneous intervals
\[
\overline q_n\in[\overline I_n-r_n,\overline I_n+r_n]\cap[0,1],
\qquad r_n=\sqrt{\frac{\log(2K/\alpha)}{2n}}.
\]
Here $K=8$ comprises the strict labor, three-variable,
employment--production and weak labor events in each of the two
fixed periods. With $\alpha=0.05$, the primary-period radius is
$0.0860$. Subtracting a model's average forecast gives its calibration
gap interval without an additional multiplicity cost. These bounds
require no correctly specified VAR, but concern fixed-period averages
rather than any individual forecast origin.
\begin{table}[htbp]
\centering\small\setlength{\tabcolsep}{4pt}
\caption{Observed Releases and Aggregate Calibration}
\label{emp:alfred:calibration}
\begin{tabularx}{\linewidth}{@{\extracolsep{\fill}}lrrrr}
\toprule
Event--period & Occurrences & Mean forecast & Probability bound & Gap bound \\
\midrule
Primary: Labor & 47/390 & 0.200 & $[0.035,0.207]$ & $[-0.165,0.007]$ \\
Primary: Three-variable & 28/390 & 0.154 & $[0.000,0.158]$ & $[-0.154,0.004]$ \\
Primary: Employment--production & 46/390 & 0.205 & $[0.032,0.204]$ & $[-0.173,-0.001]$ \\
Stress: Labor & 2/60 & 0.332 & $[0.000,0.253]$ & $[-0.332,-0.079]$ \\
Stress: Three-variable & 2/60 & 0.213 & $[0.000,0.253]$ & $[-0.213,0.040]$ \\
Stress: Employment--production & 2/60 & 0.249 & $[0.000,0.253]$ & $[-0.249,0.003]$ \\
\bottomrule
\end{tabularx}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize \textit{Note:}
Forecasts are archived Gaussian VAR(1). The gap is the average true conditional probability minus the average issued forecast. Negative values indicate overprediction. Bounds are simultaneous at 95\% across eight event--period pairs, including weak labor rises used in the measurement bracket. They require adapted binary outcomes, not independent months or a correctly specified VAR. They are not monthly confidence intervals.
\end{minipage}
\end{table}
The simultaneous labor-gap interval includes zero. For
employment--production, the upper endpoint is $-0.0014$: aggregate
overprediction is detected, though the exclusion is narrow. In the
stress period, the labor event occurs twice and the archived Gaussian
forecast averages $0.3315$. Its average overprediction is also detected.
Changing the outcome vintage changes the economic question.
September 2026 data reclassify 21 primary labor dates, adding nine
events and removing twelve. They reclassify 25 employment--production
dates, even though the net event count changes only from 46 to 45.
This is why the release target must be fixed before comparing
forecasts. It is also why a latest-vintage series truncated at each
historical month cannot reproduce the information available at the
original forecast date.
Model limitations become particularly clear in the stress period.
The archived Gaussian labor interval has median width $0.989$, and one
date is numerically degenerate. Clayton estimates reach a parameter
boundary at 56 of 60 origins, where intervals based on the interior
theory are withheld. Large changes in the scale of issued errors and
persistence in their squares point to a marginal modeling problem.
The stochastic-volatility results of \citet{clark_2011} are relevant
here: a covariance correction for estimation error cannot make a
constant-scale law accommodate changing volatility.
Figure~\ref{emp:alfred:path} retains the full stress period.
\begin{figure}[!htb]
\centering
\includegraphics[width=0.92\textwidth]{release_forecasts.pdf}
\caption{Reconstructed event forecasts and model-based 95\% logit
intervals. Crosses mark recorded events. The vertical line starts the
stress period. Intervals quantify parameter uncertainty in the fitted
working law.}
\label{emp:alfred:path}
\end{figure}
The release record separates the uncertainty calculation from model
adequacy. Within the working law, current marginal estimation contributes
much more than the extra archive term. Across the observed releases,
the diagnostics reveal limitations of that law that the uncertainty
calculation does not resolve. Valid inference at an individual origin under an adapting
marginal model would require a new first-stage analysis.
\section{Conclusion}\label{sec:conclusion}
Inference for a joint-event forecast depends on both the marginal fits
that generated the dependence-training record and the fit used for the
next forecast. In the Vine Copula VAR, the historical regression
correction reduces to harmonically weighted innovation moments, whereas
the terminal forecast retains uncertainty in every coefficient. Their
joint influence yields a feasible covariance estimator and
repeated-sample intervals for conditional event probabilities at a
random forecast state. For fixed one-sided events, nondegeneracy follows
under the maintained Gaussian and positive Clayton specifications.
The Gaussian refit comparison assigns a precise variance cost to
retaining issued transforms on the same history. In the simulations
and real-time forecasting application, that cost is smaller than the contribution of
current marginal estimation. Logit intervals improve lower-tail
coverage in the reported designs, while exact training weights make
the finite-prefix calculation explicit. Both procedures have a
first-order justification. The larger-sample results support the
asymptotic approximation, but coverage shortfalls remain at the shorter
histories relevant for the application.
The real-time forecasts point to a substantive next step. Changes in marginal
scales and serial dependence in issued errors limit the stationary
benchmark, even when its estimation uncertainty is fully accounted for.
Allowing volatility to adapt requires the historical correction and
terminal influence to be derived for that marginal estimator. Discrete
release distributions and explicit revision models pose related
problems. In each case, inference must follow the forecasts actually
used to construct the historical record as well as the estimates used
at the current origin.
\clearpage