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.
54,579 characters
\runningtitle{Causal Lag Structure Discovery in Confounded Time Series}
\runningauthor{Tan, Zanoria, Chen, Lyu, Cucuringu}
\twocolumn[
\aistatstitle{Causal Lag Structure Discovery in Confounded\\Time Series via Orthogonalized Adaptive Estimation}
\aistatsauthor{Hong Kiat Tan$^{1,*}$ \And Isaac-Neil Zanoria$^{2}$ \And James Chen$^{1}$ \And Haoyang Lyu$^{1}$ \And Mihai Cucuringu$^{1,3,4}$}
\aistatsaddress{$^{1}$Department of Mathematics, University of California, Los Angeles\\ $^{2}$Department of Electrical and Computer Engineering, University of California, Los Angeles\\ $^{3}$Oxford-Man Institute of Quantitative Finance, University of Oxford\\ $^{4}$Department of Statistics, University of Oxford}
]
{\footnotetext{$^{*}$Correspondence to: \texttt{[email removed]}}}
\begin{abstract}
Finding which variables cause which others in multivariate time series, and at what lags, is central to science and policy, yet existing methods force a choice between flexible confounder adjustment, data-driven lag selection, and inference that controls the false discovery rate (FDR). ORACLE-VARX does all three in one pipeline. First, double/debiased machine learning (DML) removes nonlinear confounder effects from the outcomes and the lagged series. Second, adaptive causal lag estimation (ACLE) picks the lag order at each time step by sequential significance tests, tracking regime changes. Third, entry-wise $z$-tests with Benjamini--Hochberg correction select directed edges at a target FDR. We prove that in each rolling window, the debiased coefficients are asymptotically normal around a window-averaged target, so their $z$-tests are asymptotically valid. On a synthetic benchmark with time-varying structure and nonlinear confounding, ORACLE-VARX (LightGBM) tracks the true lag order best (RMSE $0.96$ vs $1.1$--$1.5$), has edge FDR $0.047$, close to PCMCI ($0.045$) and below VAR ($0.129$) and VAR-LiNGAM ($0.187$), and forecasts better than all three. On nine U.S. sector ETFs with macroeconomic confounders, it yields interpretable causal graphs whose lag order rises in high-volatility regimes.
\end{abstract}
\begin{center}
\faGithub\ \,\href{https://github.com/HK-Tan/ORACLE-VARX}{\texttt{HK-Tan/ORACLE-VARX}}
\end{center}
\section{Introduction}\label{sec:intro}
Causal discovery, learning which variables directly influence
which others from observational data, is a cornerstone of
scientific reasoning~\citep{pearl2009causality}. For time series, one must determine not only
\emph{which} variables are causally linked, but also \emph{the lags at which}
those effects operate. Such lag-annotated graphs arise across
domains: lead--lag relationships among financial assets signal
information flow and systemic
risk~\citep{lo1990contrarian,billio2012econometric}; climate
teleconnections link weather patterns across distant
regions~\citep{wallace1981teleconnections}; and neural connectivity
studies seek directed interactions at millisecond
resolution~\citep{brovelli2004beta,seth2015granger}. In each setting, recovering the causal lag structure is the first
step toward reliable prediction, intervention, and policy
evaluation.
Time series pose challenges that generic causal
discovery methods do not address. Autocorrelation inflates the
effective number of tests, complicating inference.
The true lag order may be unknown and time-varying. Moreover, exogenous confounders often enter through nonlinear channels. Classical Granger causality and VAR methods~\citep{granger1969investigating,sims1980macroeconomics,lutkepohl2005new}
assume a single, fixed lag order and either ignore confounders or
include them linearly. Temporal structure-learning methods---PCMCI~\citep{runge2019detecting},
DYNOTEARS~\citep{pamfil2020dynotears}, CD-NOTS~\citep{sadeghi2025causal},
VAR-LiNGAM~\citep{hyvarinen2010estimation}---each address some of these
issues, but none handles all of (a)~nonlinear confounding,
(b)~adaptive lag order selection, and (c)~FDR-controlled
edge discovery with entry-wise confidence intervals.
As a concrete illustration, consider a system of three endogenous time
series $Y^{(1)}_t, Y^{(2)}_t, Y^{(3)}_t$ driven by a lag-1 cycle
$Y^{(3)} \!\to\! Y^{(1)} \!\to\! Y^{(2)} \!\to\! Y^{(3)}$ (each edge with
coefficient $0.60$), where additional
self-effects at lags~2 and~3 decay linearly to zero over time,
producing three successive regimes (6~edges, then~5, then only the
lag-1 cycle). Three exogenous confounders enter through distinct
nonlinear functions (tanh, quadratic, sine), creating correlated noise
that cannot be removed by linear regression. The challenge is: given
observed trajectories of $\mathbf{Y}_t$ and $\mathbf{W}_t$, can we simultaneously
discover \emph{(i)}~which directed edges exist, \emph{(ii)}~at which
lags they operate, and
(iii) provide valid statistical inference that controls false discoveries?
Figure~\ref{fig:running-dag} depicts this data-generating process as a time-unrolled
acyclic digraph.
\begin{figure}[t]
\centering
\begin{tikzpicture}[
scale=0.72,
every node/.style={font=\small},
varnode/.style={circle, draw, thick, minimum size=6.5mm, inner sep=0pt},
confbox/.style={rectangle, draw, thick, rounded corners=3pt,
fill=gray!12, minimum width=12mm, minimum height=7mm,
inner sep=2pt, font=\small},
lag1/.style={-{Stealth[length=2.2mm]}, thick, blue!70!black},
lag2/.style={-{Stealth[length=2.2mm]}, thick, dashed, red!70!black},
lag3/.style={-{Stealth[length=2.2mm]}, thick, densely dotted,
ForestGreen!90!black},
conf/.style={-{Stealth[length=2mm]}, thick, gray!60, densely dashed},
]
\foreach \col/\lab [count=\c] in {0/{t\!-\!3}, 2.6/{t\!-\!2}, 5.2/{t\!-\!1}, 7.8/{t}} {
\node[anchor=south, font=\footnotesize\itshape] at (\col, 3.55) {$\lab$};
}
\foreach \col [count=\c] in {0, 2.6, 5.2, 7.8} {
\node[varnode] (Y1-\c) at (\col, 2.7) {$Y^{(1)}$};
\node[varnode] (Y2-\c) at (\col, 1.35) {$Y^{(2)}$};
\node[varnode] (Y3-\c) at (\col, 0) {$Y^{(3)}$};
}
\foreach \s/\e in {1/2, 2/3, 3/4} {
\draw[lag1] (Y3-\s) -- (Y1-\e)
node[midway, above, sloped, font=\scriptsize, blue!70!black] {0.60};
}
\foreach \s/\e in {1/2, 2/3, 3/4} {
\draw[lag1] (Y1-\s) -- (Y2-\e)
node[midway, above, sloped, font=\scriptsize, blue!70!black] {0.60};
}
\foreach \s/\e in {1/2, 2/3, 3/4} {
\draw[lag1] (Y2-\s) -- (Y3-\e)
node[midway, below, sloped, font=\scriptsize, blue!70!black] {0.60};
}
\draw[lag2] (Y1-1) to[bend left=42] node[above, font=\scriptsize,
red!70!black, yshift=4pt, fill=white, inner sep=1pt] {$a(t)$} (Y1-3);
\draw[lag2] (Y1-2) to[bend left=42] node[above, font=\scriptsize,
red!70!black, yshift=4pt, fill=white, inner sep=1pt] {$a(t)$} (Y1-4);
\draw[lag2] (Y2-1) to[bend left=22] node[pos=0.35, below, font=\scriptsize,
red!70!black, yshift=-5pt, xshift = -5pt, fill=white, inner sep=1pt] {$a(t)$} (Y2-3);
\draw[lag2] (Y2-2) to[bend left=22] node[pos=0.35, below, font=\scriptsize,
red!70!black, yshift=-5pt, xshift = -5pt, fill=white, inner sep=1pt] {$a(t)$} (Y2-4);
\draw[lag3] (Y3-1) to[bend right=22] node[below, font=\scriptsize,
ForestGreen!90!black, yshift=-6pt, fill=white, inner sep=1pt] {$b(t)$} (Y3-4);
\node[confbox] (W) at (3.9, -2.35) {$\mathbf{W}_t$};
\foreach \c in {1,2,3,4} {
\draw[conf] (W) -- (Y1-\c);
\draw[conf] (W) -- (Y2-\c);
\draw[conf] (W) -- (Y3-\c);
}
\matrix[anchor=north west, draw, rounded corners=2pt, inner sep=3pt,
column sep=4pt, row sep=1pt, font=\scriptsize,
nodes={anchor=west}]
at (5, -1.2) {
\draw[lag1] (0,0) -- (0.5,0); & \node{Lag-1 cycle (0.60)}; \\
\draw[lag2] (0,0) -- (0.5,0); & \node{Lag-2 self, $a(t)\!\downarrow\!0$}; \\
\draw[lag3] (0,0) -- (0.5,0); & \node{Lag-3 self, $b(t)\!\downarrow\!0$}; \\
\draw[conf] (0,0) -- (0.5,0); & \node{Nonlinear confounders}; \\
};
\end{tikzpicture}
\vspace{-2mm}
\caption{An example of a time-unrolled causal graph of a 3-variable system.
Solid blue arrows show the lag-1 cycle (coefficient~$0.60$).
Dashed red and dotted green arrows show lag-2 and lag-3
self-effects $a(t)$ and $b(t)$ that decay to zero
(Section~\ref{sec:synthetic}).
Gray dashed arrows indicate nonlinear confounder effects.}
\label{fig:running-dag}
\end{figure}
We propose ORACLE-VARX (Orthogonalized Regression for Adaptive Causal
Lag Estimation in VARX), a framework that addresses all three challenges
within a unified inferential framework. Our contributions are:
\vspace{-2mm}
\begin{enumerate}[itemsep=0pt]
\item We introduce \textbf{DML orthogonalization for VAR coefficients},
applying double/debiased machine
learning~\citep{chernozhukov2018double} to residualize both outcomes
and treatments against nonlinear confounder effects via flexible
first-stage learners and time-series cross-fitting, yielding
$\sqrt{T}$-consistent and asymptotically normal coefficient
estimates (Section~\ref{sec:dml-var}).
\item We propose \textbf{adaptive causal lag estimation (ACLE)},
a two-stage procedure that combines sequential significance testing
across a grid of significance levels with rolling validation-based
tuning to select the lag order $\hat{p}_t$ at each time step,
tracking regime changes in the underlying causal structure without
requiring $p$ to be fixed \emph{a priori}
(Section~\ref{sec:acle}).
\item We perform \textbf{FDR-controlled edge discovery} via
entry-wise $z$-tests on the debiased coefficients, corrected by
the Benjamini--Hochberg procedure, yielding directed causal graphs
with false discovery rate control at level~$q$.
\item We establish \textbf{formal guarantees}, proving
local asymptotic normality of the debiased coefficients in each
rolling window, around a window-averaged target and with valid
standard errors (Theorem~\ref{thm:dml-var}), extending the DML framework
of~\citet{chernozhukov2018double} to the VARX setting, and
connecting zero debiased coefficients to Granger non-causality
(Proposition~\ref{prop:granger}).
\item We provide \textbf{empirical validation} on both synthetic
and real data. On a synthetic benchmark with time-varying structure
and nonlinear confounding, ORACLE-VARX closely tracks the true lag
order with near-unit RMSE across five first-stage learners. On nine U.S.\ sector ETFs with macroeconomic confounders, the ACLE adaptive lag procedure alone achieves an annualized Sharpe ratio of approximately
1.3 on open-to-close returns, a setting where broad market indices
earn near-zero
returns~\citep{haghani2022nightmoves}, and the resulting causal
graphs yield interpretable lead--lag structure whose adaptive lag
order rises in high-volatility regimes (Sections~\ref{sec:synthetic}--\ref{sec:etf}).
\end{enumerate}
\vspace{-2mm}
\section{Related Work}\label{sec:related}
\vspace{-2mm}
\paragraph{Granger causality and VAR models.}
Classical Granger causality~\citep{granger1969investigating} tests
whether past values of one series improve the prediction of another
beyond the target's own history, and remains the most widely used
framework for analyzing lead--lag effects in multivariate time series.
\citet{sims1980macroeconomics} popularized vector autoregressive (VAR)
models as a flexible, data-driven tool for macroeconomic analysis, and
\citet{lutkepohl2005new} codified the statistical theory
for estimation, testing, and lag selection---mostly under linear dynamics
and a fixed, globally chosen lag order.
Sparse VAR extensions handle high-dimensional settings:
\citet{basu2015regularized} analyzed $\ell_1$-penalized estimation
with non-asymptotic guarantees under stability conditions, and
\citet{davis2016sparse} fit sparse VARs in two stages via partial
spectral coherence.
\citet{uematsu2025discovering} test each lag coefficient of a large
VAR with debiased-LASSO statistics and control the FDR across edges, but
in a linear, stationary VAR of known order.
All of these methods assume that confounders enter linearly or are
absent, and none adapts the lag order over time.
\vspace{-1mm}
\paragraph{Causal discovery for time series.}
PCMCI~\citep{runge2019detecting} combines PC-style parent
selection with conditional independence tests (optionally
nonlinear) that handle autocorrelation, but it neither adapts the active
lags over time nor gives lag coefficients with confidence intervals
(Table~\ref{tab:toy-results} compares it and VAR-LiNGAM with ours).
VAR-LiNGAM~\citep{hyvarinen2010estimation} uses non-Gaussianity to
identify contemporaneous and lagged structure, but assumes linear effects
and no hidden confounders.
TiMINo~\citep{peters2013causal} establishes identifiability under
independent noise and nonlinear relationships, but its algorithm
rules out feedback loops between series.
DYNOTEARS~\citep{pamfil2020dynotears} learns structure by
continuous optimization but gives no inference on individual edges.
CD-NOTS~\citep{sadeghi2025causal} handles nonstationarity by detecting
mechanism changes, but assumes every hidden confounder is a smooth
function of time.
Recent ML-based Granger and conditional independence tests for time
series~\citep{angelis2025doubly,wiecksosa2025dynamic} adjust for confounders
flexibly, and \citet{shi2026conditional} also search over the lag depth,
one pair of series at a time.
These tests detect links but give no lag coefficients with
confidence intervals.
\vspace{-3mm}
\paragraph{Double/debiased machine learning.}
\citet{robinson1988root} showed that in a partially linear model,
nonparametric nuisance fits (regressions on the confounders) followed by
OLS give $\sqrt{T}$-consistent estimates of the linear part.
\citet{chernozhukov2018double} generalized this into DML: Neyman
orthogonal scores and cross-fitting allow ML first stages that converge fast
enough, while keeping asymptotically normal inference.
Recent works extend DML to dependent data~\citep{ciganovic2026double,yu2026semiparametric}.
For high-dimensional VARs, \citet{hecq2023granger} and
\citet{adamek2023lasso} give LASSO-based inference with a fixed lag order.
Closest to our work, \citet{ballinari2024irfs} prove DML
consistency for impulse response functions under $\alpha$-mixing,
showing that blocked cross-fitting preserves asymptotic validity
(see Remark~\ref{rem:buffer}).
We apply their cross-fitting theory to VAR
coefficient matrices and add ACLE and BH-FDR control.
\citet{yu2026semiparametric} allow vector outcomes and treatments;
with the lags as treatments and instruments, their estimator reduces to our single-window
one. They assume time-invariant nuisance functions;
Theorem~\ref{thm:dml-var} covers drifting ones for our window average.
ORACLE-VARX is, to our knowledge, the first to combine
nonparametric DML debiasing, adaptive lag discovery, and FDR-controlled edge selection, with a confidence interval per lag coefficient.
\vspace{-2mm}
\section{Methodology}\label{sec:method}
\vspace{-2mm}
\subsection{The VARX Model and Its Partially Linear Extension}\label{sec:setup}
\vspace{-1mm}
Assume that we have $n$ endogenous time series with
$\mathbf{Y}_t = (Y^{(1)}_t,\ldots,Y^{(n)}_t)^\top \in \mathbb{R}^n$ for
$t=1,\ldots,T$, together with $d$-dimensional exogenous confounders
$\mathbf{W}_t = (W_{1,t},\ldots,W_{d,t})^\top \in \mathbb{R}^d$.
\paragraph{Linear VARX.}
The standard linear VARX($p$) model is
\begin{equation}\label{eq:varx}
\mathbf{Y}_t
= \sum_{k=1}^{p} \mathbf{A}_k\,\mathbf{Y}_{t-k}
+ \sum_{k=1}^{p} \mathbf{B}_k\,\mathbf{W}_{t-k}
+ \boldsymbol{\varepsilon}_t,
\end{equation}
where $\mathbf{A}_k \in \mathbb{R}^{n\times n}$ are lag-$k$ coefficient matrices
encoding causal effects among the endogenous variables,
$\mathbf{B}_k \in \mathbb{R}^{n\times d}$ capture linear confounder effects,
and $\boldsymbol{\varepsilon}_t$ is a zero-mean innovation.
For simplicity, we use the same lag order~$p$ for both endogenous and exogenous terms.
\paragraph{Partially linear VARX.}
To allow flexible confounder adjustment, we generalize to a
\emph{partially linear VARX}:
\begin{equation}\label{eq:plvar}
\mathbf{Y}_t
= \sum_{k=1}^{p^*} \mathbf{A}_k\,\mathbf{Y}_{t-k}
+ g(\mathbf{W}_{t-1},\ldots,\mathbf{W}_{t-p})
+ \boldsymbol{\varepsilon}_t,
\end{equation}
where $g:\mathbb{R}^{d \cdot p}\!\to\mathbb{R}^n$ is an unknown, potentially nonlinear
nuisance function,
$p^*$ is the true (unknown) maximum causal lag order, $p\ge p^*$ is the lag order used in the model, and
$\mathbb{E}[\boldsymbol{\varepsilon}_t \mid \mathbf{Y}_{t-1},\ldots,\mathbf{W}_{t-p}]=\mathbf{0}$.
The effect of lagged $\mathbf{Y}$ on current $\mathbf{Y}$ remains
linear---yielding interpretable, testable coefficients---while the
confounder effect is left fully flexible.
This is the model underlying OR-VARX and ORACLE-VARX, and we can recover
the standard linear VARX~\eqref{eq:varx} when $g$ is linear.
\begin{definition}[Causal Lag Graph]\label{def:lag-graph}
The \emph{causal lag graph}
$\mathcal{G}=(V,E)$ has vertex set
$V=\{Y^{(1)},\ldots,Y^{(n)}\}$
and directed, lag-annotated edge set
$E=\bigl\{(Y^{(j)}\!\to\! Y^{(i)},\,k): \mathbf{A}_k[i,j]\neq 0\bigr\}$.
An edge $(Y^{(j)}\!\to\! Y^{(i)},\,k)$ means $Y^{(j)}$ at lag~$k$
Granger-causes $Y^{(i)}$.
\end{definition}
\noindent
\textbf{Goal.}\;
Recover $\mathcal{G}$ from observational data---that is, determine
which entries $\mathbf{A}_k[i,j]$ are nonzero, at which lags $k$, and for
which directed pairs $(j,i)$ while providing valid statistical inference for each edge, while making minimal assumptions on the nuisance function~$g$.
\subsection{DML for VAR: Orthogonalization Procedure}\label{sec:dml-var}
Following the double/debiased machine learning (DML) framework
of~\citet{chernozhukov2018double}, we cast the partially linear
VARX~\eqref{eq:plvar} into a treatment-effect setup.
\paragraph{DML formulation.}
Define the \emph{stacked treatment} and \emph{stacked confounders}:
\begin{equation}\label{eq:stacked}
\mathbf{T}_t \;:=\;
\begin{pmatrix} \mathbf{Y}_{t-1} \\ \vdots \\ \mathbf{Y}_{t-p} \end{pmatrix}
\!\in \mathbb{R}^{np},
\quad
\mathbf{C}_t \;:=\;
\begin{pmatrix} \mathbf{W}_{t-1} \\ \vdots \\ \mathbf{W}_{t-p} \end{pmatrix}
\!\in \mathbb{R}^{dp}.
\end{equation}
Let $\boldsymbol{\Theta} \in \mathbb{R}^{np\times n}$ denote the stacked coefficient matrix
formed by vertically stacking $\mathbf{A}_1^\top,\ldots,\mathbf{A}_p^\top$, so that $\boldsymbol{\Theta}^\top\mathbf{T}_t=\sum_{k}\mathbf{A}_k\mathbf{Y}_{t-k}$.
Then~\eqref{eq:plvar} becomes the partially linear model
\begin{align}
\mathbf{Y}_t &= \boldsymbol{\Theta}^\top \mathbf{T}_t + g(\mathbf{C}_t) + \mathbf{U}_t, \label{eq:dml-Y}\\
\mathbf{T}_t &= g_T(\mathbf{C}_t) + \mathbf{V}_t, \label{eq:dml-T}
\end{align}
where $g:\mathbb{R}^{dp}\!\to\mathbb{R}^n$ and $g_T:\mathbb{R}^{dp}\!\to\mathbb{R}^{np}$ are
nuisance functions,
$\mathbf{U}_t$ is the structural innovation with $\mathbb{E}[\mathbf{U}_t\mid\mathbf{C}_t,\mathbf{T}_t]=\mathbf{0}$,
and $\mathbf{V}_t$ is the treatment residual with $\mathbb{E}[\mathbf{V}_t\mid\mathbf{C}_t]=\mathbf{0}$.
The residualization below uses the outcome's prediction from the
confounders, $g_Y(\mathbf{C}_t):=\mathbb{E}[\mathbf{Y}_t\mid\mathbf{C}_t]=\boldsymbol{\Theta}^\top g_T(\mathbf{C}_t)+g(\mathbf{C}_t)$.
The parameter of interest is $\boldsymbol{\Theta}$ (equivalently, the individual
$\mathbf{A}_k$ matrices), and the nuisance functions $g_Y$, $g_T$ are treated as
infinite-dimensional parameters that must be estimated.
\paragraph{Three-step DML procedure.}
\begin{enumerate}[leftmargin=*,itemsep=2pt]
\item \textbf{Cross-fitted nuisance estimation.}
Using flexible machine learning methods (extra trees, gradient boosting,
random forests, or a foundational model), estimate $\hat{g}_Y$ and $\hat{g}_T$ via
\emph{time-series cross-fitting}.
\item \textbf{Residualization.}
Form the residualized outcome and treatment:
\begin{equation}\label{eq:residuals}
\tilde{\mathbf{Y}}_t := \mathbf{Y}_t - \hat{g}_Y(\mathbf{C}_t), \quad
\tilde{\mathbf{T}}_t := \mathbf{T}_t - \hat{g}_T(\mathbf{C}_t).
\end{equation}
\item \textbf{Orthogonalized OLS.}
Estimate the debiased coefficient matrix by ordinary least squares on
the residualized data:
\begin{equation}\label{eq:theta-hat}
\hat{\boldsymbol{\Theta}}
= \Bigl(\textstyle\sum_{t}\tilde{\mathbf{T}}_t\,\tilde{\mathbf{T}}_t^\top\Bigr)^{-1}
\Bigl(\textstyle\sum_{t}\tilde{\mathbf{T}}_t\,\tilde{\mathbf{Y}}_t^\top\Bigr).
\end{equation}
\end{enumerate}
The estimator~\eqref{eq:theta-hat} satisfies Neyman orthogonality
(Section~\ref{sec:theory}), so first-stage estimation errors in
$\hat{g}_Y$ and $\hat{g}_T$ affect $\hat{\boldsymbol{\Theta}}$ only at second order.
This permits flexible ML learners in the first stage without
sacrificing inferential validity.
\paragraph{Time-series cross-fitting.}
Note that standard $K$-fold cross-validation violates temporal ordering and induces
dependence between training and test sets. Hence,
we instead use a \emph{blocked rolling window split}:
fold~$k$ trains on the interval
$[s_k,\;s_k+n_{\mathrm{train}})$ and tests on
$[s_k+n_{\mathrm{train}},\;s_k+n_{\mathrm{train}}+n_{\mathrm{test}})$,
with $s_{k+1}=s_k+n_{\mathrm{test}}$ so that folds tile without overlap.
Nuisance residuals $\tilde{\mathbf{Y}}_t$ and $\tilde{\mathbf{T}}_t$
are computed and aggregated only on the test folds, for which
the second-stage OLS~\eqref{eq:theta-hat} then uses only these out-of-sample
residuals accordingly.
We \emph{memoize} the first stage: because the blocks lie on one grid shared by all windows, each block is fit once and reused by every window that overlaps it. This needs over $200\times$ fewer first-stage fits than recomputing each window, and a full synthetic run (all 2,595 windows, one tree-based learner) takes 6--17 minutes on a 10-core Apple M5 laptop CPU, with no GPU (Appendix~\ref{app:compute}).
\paragraph{Standard errors and inference.}
After obtaining $\hat{\boldsymbol{\Theta}}$, we compute standard errors using the
classical OLS variance formula on the residualized data:
\begin{equation}\label{eq:ols-se}
\widehat{\operatorname{Var}}(\hat{\boldsymbol{\Theta}}_{\cdot i})
= \hat{\sigma}_i^2\,
\Bigl(\textstyle\sum_{t}\tilde{\mathbf{T}}_t\,\tilde{\mathbf{T}}_t^\top\Bigr)^{-1},
\end{equation}
for each equation $i$, where $\hat{\boldsymbol{\Theta}}_{\cdot i}$ is the $i$-th column of $\hat{\boldsymbol{\Theta}}$,
$\hat{\sigma}_i^2 = (T_{\mathrm{eff}}-np)^{-1}\sum_t
(\tilde{Y}_{t,i} - \hat{\boldsymbol{\Theta}}_{\cdot i}^\top\tilde{\mathbf{T}}_t)^2$
and $T_{\mathrm{eff}}$ is the number of out-of-sample residualized
observations. Under conditionally homoskedastic innovations this is the
correct asymptotic variance, with no HAC correction
(Theorem~\ref{thm:dml-var}).
The diagonal of $\widehat{\operatorname{Var}}(\hat{\boldsymbol{\Theta}})$ then provides entry-wise standard
errors, enabling $z$-tests for the hypothesis
$H_0:\mathbf{A}_k[i,j]=0$ on every potential edge in the causal graph.
\begin{remark}[Buffer Zones]\label{rem:buffer}
Time-series DML theory often inserts buffer zones of $k_T$ observations
(e.g.\ $k_T=\lfloor T/K\rfloor$) between training and test folds
\citep{ballinari2024irfs}.
We use adjacent windows ($k_T=0$). The nuisance functions themselves drift
over time, so a buffer would make every nuisance fit staler, and it would drop
$k_T$ observations from every training fold. Our theory
(Section~\ref{sec:theory}) does not use a buffer: it assumes directly that
the cross-fitted nuisance errors are negligible (condition (N),
Appendix~\ref{app:assumptions}).
\end{remark}
\vspace{0mm}
\subsection{Lag Order Selection}\label{sec:acle}
\vspace{0mm}
The lag order $p$ controls the number of coefficient matrices
$\mathbf{A}_1,\ldots,\mathbf{A}_p$ to estimate.
Fixing $p$ \emph{a priori} is problematic: setting it too low misses
genuine edges at higher lags, while setting it too high inflates the
parameter count and increases false positives.
The true lag order may also shift over time.
We consider two data-driven approaches.
\vspace{0mm}
\paragraph{RMSE-based lag selection.}
For methods without significance testing (VAR, VARX, OR-VARX),
we iterate over candidate lags $p \in \{1,\ldots,p_{\max}\}$,
compute rolling forecasts $\hat{\mathbf{Y}}_s^{(p)}$ for each,
and select the lag minimizing root mean squared forecast error
on a trailing validation window of $v$ days.
Concretely, define the rolling validation RMSE for any candidate
lag~$p$ as
\begin{equation}\label{eq:rmse-select}
\mathrm{RMSE}_v(p;\,t)
= \sqrt{\frac{1}{v}\sum_{s=t-v}^{t-1}
\bigl\|\mathbf{Y}_s - \hat{\mathbf{Y}}_s^{(p)}\bigr\|^2},
\end{equation}
and set $\hat{p}_t = \operatorname*{arg\,min}_{p \in \{1,\ldots,p_{\max}\}}
\mathrm{RMSE}_v(p;\,t)$.
This is simple and assumption-free but ignores the statistical evidence for individual edges.
\vspace{0mm}
\paragraph{Adaptive Causal Lag Estimation (ACLE).}
ACLE uses statistical significance to
decide when to stop adding lags, so it tracks causal signal
strength, not only forecast accuracy.
It works with both OLS-VAR and DML coefficients.
\vspace{0mm}
\paragraph{Stage~1: Significance-based lag selection.}
For each rolling window centred at time~$t$, we evaluate candidate
lags $p \in \{2,\ldots,p_{\max}\}$ via sequential hypothesis testing.
Figure~\ref{fig:acle-dag} illustrates the procedure.
\begin{figure}[t]
\centering
\begin{tikzpicture}[
scale=0.45,
circ/.style={circle, draw, thick, minimum size=0.8cm, align=center, font=\small},
wcircle/.style={circle, draw, thick, minimum size=0.8cm, align=center, font=\small},
rect/.style={rectangle, draw, thick, minimum size=0.8cm, rounded corners, align=center, font=\small},
arrow/.style={->, thick},
dashed_arrow/.style={->, thick, densely dashed}
]
\begin{scope}
\node[font=\bfseries] at (3, 7.5) {Lags $1$ to $p{-}1$};
\node[circ] (T_pre) at (0,0) {$T^{(p-1)}$};
\node[wcircle] (W_pre) at (3,5) {$W$};
\node[circ] (Y_pre) at (6,0) {$\mathbf{Y}_t$};
\draw[arrow] (W_pre) -- (T_pre);
\draw[arrow] (W_pre) -- (Y_pre);
\draw[arrow] (T_pre) -- (Y_pre);
\node[] at (3,0.6) {$\mathbf{A}_{1:p-1}$};
\end{scope}
\begin{scope}[xshift=10cm]
\node[font=\bfseries] at (3, 7.5) {Testing lag $p$};
\node[circ] (T) at (0,0) {$T^{(p-1)}$};
\node[wcircle] (W) at (3,5) {$W$};
\node[circ] (Y) at (6,0) {$\mathbf{Y}_t$};
\node[rect] (R) at (0,-4) {$\mathbf{Y}_{t-p}$};
\draw[arrow] (W) -- (T);
\draw[arrow] (W) -- (Y);
\draw[arrow] (T) -- (Y);
\node[] at (3,0.6) {$\mathbf{A}_{1:p-1}$};
\draw[arrow, densely dotted] (R) -- (T);
\draw[dashed_arrow] (R) -- (Y);
\node[] at (4,-2.5) {$\mathbf{A}_p$};
\end{scope}
\end{tikzpicture}
\caption{Sequential lag testing in ACLE.
Left: the model with lags $1,\ldots,p{-}1$.
Right: testing $H_0\colon \mathbf{A}_p = 0$ on the lag-$p$ block.
If rejected at level $\alpha$, lag $p$ is retained and we proceed to
$p{+}1$; otherwise we stop and set
$p_\alpha = p{-}1$.}
\label{fig:acle-dag}
\end{figure}
\vspace{0mm}
For each candidate lag $p \in \{2,\ldots,p_{\max}\}$, we test
\[
H_0\colon \mathbf{A}_p[i,j] = 0 \quad \text{for all } (i,j) \text{ pairs simultaneously.}
\]
Entry-wise $z$-statistics
$z_{i,j} = \hat{\mathbf{A}}_p[i,j] \,/\, \hat s(\hat{\mathbf{A}}_p[i,j])$
are computed from the $p_{\max}$ fit, and Benjamini--Hochberg FDR correction~\citep{benjamini1995controlling} at level~$\alpha$
determines whether any entry of $\mathbf{A}_p$ is significant.
If none is, we stop and set
$p_\alpha = p - 1$.
This yields a candidate lag $p_\alpha$ for each $\alpha$. Because $\alpha$ governs the sensitivity of
lag detection, we evaluate a grid
$\alpha \in \{\alpha_1,\ldots,\alpha_M\}$, producing candidate lags $\{p_{\alpha_1},\ldots,p_{\alpha_M}\}$.
\vspace{0mm}
\paragraph{Stage~2: Validation-based $\alpha$-selection.}
To select the optimal significance level~$\alpha^*_t$ without an
\emph{a priori} choice, we apply the same rolling validation
criterion~\eqref{eq:rmse-select}, now over the $\alpha$-grid
rather than over raw lag candidates:
\begin{equation}\label{eq:alpha-select}
\alpha^*_t
= \operatorname*{arg\,min}_{\alpha \in \{\alpha_1,\ldots,\alpha_M\}}
\mathrm{RMSE}_v(p_\alpha;\,t),
\end{equation}
where $p_\alpha$ is the lag selected by Stage~1 at level~$\alpha$.
The final adaptively selected lag is
$\hat{p}_t = p_{\alpha^*_t}$.
The selected lag $\hat{p}_t$ may vary as the
rolling window advances: when higher-order self-effects decay to zero,
ACLE tracks the change by reducing $\hat{p}_t$.
\vspace{0mm}
\subsection{ORACLE-VARX: Full Pipeline}\label{sec:oracle-varx}
\vspace{0mm}
The complete ORACLE-VARX pipeline combines the three components above
into a unified framework for causal lag discovery.
For each rolling window at time~$t$:
\begin{enumerate}[leftmargin=*,itemsep=1pt]
\item \textbf{Input:}\;
Time series $\{\mathbf{Y}_s,\mathbf{W}_s\}$ in the current window, maximum lag
$p_{\max}$, first-stage learner~$\mathcal{L}$.
\item \textbf{Nuisance estimation:}\;
For each candidate $p\in\{1,\ldots,p_{\max}\}$, run DML with rolling
cross-fitting (Section~\ref{sec:dml-var}) to estimate nuisance
functions $\hat{g}_Y^{(p)}$ and $\hat{g}_T^{(p)}$, then form residuals
$\tilde{\mathbf{Y}}_t^{(p)}$ and $\tilde{\mathbf{T}}_t^{(p)}$.
\item \textbf{Orthogonalized regression:}\;
Compute $\hat{\boldsymbol{\Theta}}^{(p)}$ and its standard errors
via~\eqref{eq:theta-hat}--\eqref{eq:ols-se} for each candidate~$p$.
\item \textbf{Lag selection:}\;
For non-adaptive methods, select
$\hat{p}_t$ via RMSE minimization~\eqref{eq:rmse-select}.
For ACLE methods, run sequential significance testing on the $p_{\max}$ fit across the
$\alpha$-grid (Stage~1), then select $\alpha^*_t$ via rolling
validation~\eqref{eq:alpha-select} (Stage~2) to obtain
$\hat{p}_t = p_{\alpha^*_t}$.
\item \textbf{BH-FDR edge discovery:}\;
From the $p_{\max}$ fit $\hat{\boldsymbol{\Theta}}^{(p_{\max})}$, compute $z$-scores
$z_{k,i,j} = \hat{\mathbf{A}}_k[i,j]\,/\,\hat s(\hat{\mathbf{A}}_k[i,j])$
for all $k\in\{1,\ldots,\hat{p}_t\}$ and $(i,j)$ pairs.
Apply the Benjamini--Hochberg procedure~\citep{benjamini1995controlling} at level $q$ across
all $n^2\hat{p}_t$ tests.
\item \textbf{Output:}\;
Causal lag graph $\hat{\mathcal{G}}_t$ with lag-annotated edges,
and day-$t$ forecasts $\hat{\mathbf{Y}}_t$.
\end{enumerate}
The full pseudocode is given in Algorithm~\ref{alg:oracle} (Appendix~\ref{app:algorithm}).
\vspace{0mm}
\subsection{Method Hierarchy}\label{sec:hierarchy}
Table~\ref{tab:hierarchy} displays the nested hierarchy of methods
evaluated in this paper. The experiments in Section~\ref{sec:experiments} evaluate every row of
Table~\ref{tab:hierarchy}, allowing us to attribute performance gains
to each individual component.
\begin{table}[h!]
\centering
\caption{Hierarchy of methods. Each row adds one capability
relative to the row above. ``RMSE'' lag selection picks the lag order
minimizing rolling validation RMSE~\eqref{eq:rmse-select} over
$\{1,\ldots,p_{\max}\}$;
``ACLE'' uses significance-based sequential testing with
validation-based $\alpha$-tuning
(Section~\ref{sec:acle}).}\label{tab:hierarchy}
\small
\begin{tabular}{@{}lccc@{}}
\toprule
\textbf{Method}
& \textbf{Lag Sel.}
& \textbf{Conf.\ Adj.}
& \textbf{DML} \\
\midrule
VAR
& RMSE
& None
& -- \\
VARX
& RMSE
& Linear
& -- \\
ACLE-VAR
& ACLE
& None
& -- \\
ACLE-VARX
& ACLE
& Linear
& -- \\
OR-VARX
& RMSE
& Nonlinear
& \checkmark \\
ORACLE-VARX
& ACLE
& Nonlinear
& \checkmark \\
\bottomrule
\end{tabular}
\end{table}
\vspace{0mm}
\section{Theoretical Analysis}\label{sec:theory}
The coefficients in our setting drift over time, so we study a single
rolling window. We use the stacked form of Section~\ref{sec:dml-var} with
$p=p_{\max}$ lags and let the coefficients vary slowly:
$\boldsymbol{\Theta}_t=\boldsymbol{\Theta}(t/T)$ for a function
$\boldsymbol{\Theta}:[0,1]\to\mathbb{R}^{np\times n}$, so that
$\mathbf{Y}_t=\boldsymbol{\Theta}_t^\top\mathbf{T}_t+g_t(\mathbf{C}_t)+\boldsymbol{\varepsilon}_t$.
Fix $u_0\in(0,1]$ and the window
$\mathcal W=\{t_0-T_w+1,\ldots,t_0\}$, $t_0=\lfloor u_0T\rfloor$, on which
the second-stage OLS~\eqref{eq:theta-hat} is run.
\paragraph{What a window estimates.}
Let $\mathbf{V}_t=\mathbf{T}_t-\mathbb{E}[\mathbf{T}_t\mid\mathbf{C}_t]$ be the part of the lags that the
confounders cannot predict, and $J_t=\mathbb{E}[\mathbf{V}_t\mathbf{V}_t^\top]$. The window
estimate targets
\begin{equation}\label{eq:target}
\bar\boldsymbol{\Theta}=\Bigl(\sum_{t\in\mathcal W}J_t\Bigr)^{-1}
\sum_{t\in\mathcal W}J_t\,\boldsymbol{\Theta}_t,
\end{equation}
a weighted average of the coefficients over the window. These weights make
the drift cancel exactly in the estimation error, since
$\sum_{t\in\mathcal W}J_t(\boldsymbol{\Theta}_t-\bar\boldsymbol{\Theta})=0$, and $\bar\boldsymbol{\Theta}$
lies within $O(T_w/T)$ of $\boldsymbol{\Theta}(u_0)$ (Lemma~\ref{lem:A-target}).
\begin{assumption}\label{asm:main}
We assume the following conditions:
\begin{itemize}[wide=0pt,labelsep=0.4em,itemsep=0pt,topsep=-\parskip,parsep=0pt,partopsep=0pt]
\item[\textbf{(M)}] $\mathbb{E}[\boldsymbol{\varepsilon}_t\mid\mathcal F_{t-1}]=\mathbf{0}$, where
$\mathcal F_{t-1}=\sigma(\mathbf{Y}_s,\mathbf{W}_s:s<t)$: given the observed past, the
innovation is unpredictable (all confounders are observed and
$p_{\max}\ge p^*$).
\item[\textbf{(S)}] $\boldsymbol{\Theta}$ is Lipschitz and $\lambda_{\min}(J_t)\ge c_J>0$:
coefficients drift slowly and confounder-adjusted lags are not collinear.
\item[\textbf{(W)}] $T_w\to\infty$ and $T_w=O(T^{2/3})$.
\item[\textbf{(R)}] The regularity conditions (D), (H), (N) of
Appendix~\ref{app:assumptions} hold: bounded moments and short memory,
conditionally homoskedastic innovations, and the usual DML requirement that
the cross-fitted nuisance errors are small (a product rate) and
uncorrelated with $\boldsymbol{\varepsilon}_t$ and $\mathbf{V}_t$.
\end{itemize}
\end{assumption}
\begin{theorem}[Local asymptotic normality]\label{thm:dml-var}
Under Assumption~\ref{asm:main}, for every entry $(r,i)$ of $\boldsymbol{\Theta}$,
with $\hat s_{ri}$ the standard error from~\eqref{eq:ols-se},
\[
\frac{\hat{\boldsymbol{\Theta}}_{ri}-\bar\boldsymbol{\Theta}_{ri}}{\hat s_{ri}}\xrightarrow{d}\mathcal{N}(0,1).
\]
If moreover $T_w=o(T^{2/3})$, then $\bar\boldsymbol{\Theta}_{ri}$ may be replaced by
$\boldsymbol{\Theta}_{ri}(u_0)$.
\end{theorem}
\paragraph{Proof idea.}
We first write the error exactly as
$\sqrt{T_w}(\hat{\boldsymbol{\Theta}}-\bar\boldsymbol{\Theta})=\hat J^{-1}(S_0+R)$, where
$\hat J=T_w^{-1}\sum_{t\in\mathcal W}\tilde\mathbf{T}_t\tilde\mathbf{T}_t^\top$
is close to the average of the $J_t$. The leading term
$S_0=T_w^{-1/2}\sum_{t\in\mathcal W}\mathbf{V}_t\boldsymbol{\varepsilon}_t^\top$ is a sum of
martingale differences, so a martingale central limit theorem applies with
no autocovariance (HAC) terms. The remainder $R$ collects the drift, whose
mean cancels by~\eqref{eq:target}, and the nuisance errors, which enter
only through products or terms uncorrelated with the innovations (Neyman
orthogonality), so they vanish under~(N) (Appendix~\ref{app:thm1-proof}).
Last, $T_w\hat s_{ri}^2$ is close to the limiting variance, since the residual
variance is consistent and the noise variance drifts slowly
(Appendix~\ref{app:se}).
\paragraph{Lag selection.}
ACLE is a heuristic; we do not prove that it selects the true lag.
Theorem~\ref{thm:dml-var} explains why it tracks the lag order: if
$\mathbf{A}_k=0$ on the window for all lags above some $p^*$, their
$z$-statistics stay of order one, while every lag block with
$\mathbf{A}_k(u_0)\neq0$ has an entry whose $\lvert z\rvert$ grows like
$\sqrt{T_w}$. With a fixed Stage~1 level $\alpha$, a zero lag is still
added with probability bounded away from zero, so $\alpha$ must shrink
with $T_w$; if it shrinks slowly enough that its cutoff grows more slowly
than $\sqrt{T_w}$, true lags are still detected. ACLE does not fix
$\alpha$: Stage~2 picks it from a grid by validation, letting the data set
the level. A formal consistency result covering the BH step and this
data-driven $\alpha$ is left to future work.
\paragraph{From coefficients to causal edges.}
Theorem~\ref{thm:dml-var} gives valid $p$-values for each entry of
$\bar\boldsymbol{\Theta}$ (of $\mathbf{A}_k(u_0)$ if $T_w=o(T^{2/3})$); BH's FDR guarantee
also needs independent or positively dependent (PRDS)
$p$-values~\citep{benjamini1995controlling,benjamini2001control}.
Under causal sufficiency, meaning every common cause of the series is
observed in $\mathbf{W}_t$, a zero entry of $\mathbf{A}_k$ is equivalent to Granger
non-causality at lag $k$ (Proposition~\ref{prop:granger}, given a non-degeneracy condition). A hidden
variable breaks this, because its effect is absorbed into the lag
coefficients (Figure~\ref{fig:causal-suff}): a hidden common cause of two
series creates a spurious edge, and a hidden mediator makes an indirect link look direct.
The coefficients then estimate the best linear lead--lag predictor from
the observed variables, not direct causal effects, and their $p$-values
lose the guarantee of Theorem~\ref{thm:dml-var}. This happens for VAR and ACLE-VAR, which ignore the confounders
in our benchmark, and for the ETFs, where many common causes are
unobserved (Section~\ref{sec:etf}).
\vspace{-1mm}
\section{Experiments}\label{sec:experiments}
We evaluate the full method hierarchy from Table~\ref{tab:hierarchy} on two complementary benchmarks: a synthetic
time-varying system with known ground truth (Section~\ref{sec:synthetic}), and
a financial application to ETF markets with macroeconomic confounders
(Section~\ref{sec:etf}). The synthetic benchmark is our primary experiment,
providing controlled evaluation of lag recovery, coefficient estimation, and
false positive control. The financial application demonstrates that the
framework discovers interpretable causal structure in real data.
\vspace{-1mm}
\subsection{Synthetic Benchmark}\label{sec:synthetic}
\vspace{-1mm}
\paragraph{Data-generating process.}
We construct a three-variable ($n = 3$) partially linear VARX process whose
causal structure is depicted in Figure~\ref{fig:running-dag}. The lag-1
coefficient matrix is
\[
\mathbf{A}_1 = \begin{pmatrix} 0 & 0 & 0.60 \\ 0.60 & 0 & 0 \\ 0 & 0.60 & 0 \end{pmatrix},
\]
and time-varying self-effects at lags 2 and 3 decay linearly to zero:
$\mathbf{A}_2[1,1](t) = \mathbf{A}_2[2,2](t) = 0.30 \cdot \max(0,\, 1 - t/2000)$ and
$\mathbf{A}_3[3,3](t) = 0.30 \cdot \max(0,\, 1 - t/1000)$, creating three regimes
($p^*(t) = 3, 2, 1$). Three independent AR(1) confounders
($\rho_W = 0.5$) enter through nonlinear nuisance functions:
\begin{align}
g_1(\mathbf{W}) &= \tfrac{\tanh(2W_1) + W_2^2 + \sin(\pi W_3)}{2.2}, \notag\\
g_2(\mathbf{W}) &= \tfrac{\mathrm{ReLU}(W_1 - 0.5) + \sin(\pi W_2) + W_3^2}{2.1}, \label{eq:nuisance-fns}\\
g_3(\mathbf{W}) &= \tfrac{W_1 W_2 + \tanh(W_3) + \cos(\pi W_1)}{1.7}, \notag
\end{align}
each entering the $i$-th structural equation as $\nu\lambda \cdot g_i(\mathbf{W}_{t-1})$.
Key parameters: $T = 3000$, noise scale $\nu = 0.1$, confounder strength
$\lambda = 1.0$.
\vspace{0mm}
\paragraph{Observability conditions.}
All three confounders always drive the data; the observability level
sets how many of them a method is given: \textbf{All} (3 confounders observed), \textbf{Partial-2}
(2 observed), \textbf{Partial-1} (1 observed), and \textbf{None}
(no confounders; of our methods, only VAR and ACLE-VAR apply). Below All,
confounders that drive the data are left out of the model, which violates
causal sufficiency (Section~\ref{sec:theory}), so the FDR guarantee for
causal edges no longer applies. The main text therefore reports All, and
keeps VAR and ACLE-VAR (None) only as the reference for ignoring the
confounders altogether; Appendix~\ref{app:partial-obs} reports Partial-2
and Partial-1.
\vspace{0mm}
\paragraph{Methods.}
We evaluate all six method variants from Table~\ref{tab:hierarchy},
plus PCMCI (linear test, with an OLS refit) and VAR-LiNGAM, refit in every
window (Appendix~\ref{app:cd-baselines}).
For DML methods, we compare five first-stage learners: Extra Trees,
Random Forest, LightGBM, XGBoost, and TabPFN-2.5~\citep{hollmann2025tabpfn,grinsztajn2025tabpfn}.
TabPFN is a transformer pretrained on millions of synthetic tabular
datasets that performs inference via in-context learning,
achieving strong performance on small tabular datasets without
hyperparameter tuning, which is the regime of our DML cross-fitting
folds. We set $p_{\max} = 5$, a value larger than the true maximum lag of the data generating procedure.
\vspace{-1mm}
\paragraph{Metrics.}
We evaluate five complementary metrics:
\begin{enumerate}
\item \textbf{Lag RMSE}: $\sqrt{T_{\mathrm{out}}^{-1} \sum_{t}
(\hat{p}_t - p^*(t))^2}$, measuring how well the method tracks the
true time-varying lag order. For RMSE-based methods, $\hat{p}_t$ is the
validation-selected lag.
\item \textbf{Coefficient MAE} (nonzero edges): average
$|\hat{\mathbf{A}}_k[i,j] - \mathbf{A}_k^*[i,j](t)|$ over true nonzero entries,
measuring signal recovery.
\item \textbf{Forecast MAE}: out-of-sample prediction error.
\item \textbf{FDR} and \textbf{Power}: the share of edges in
$\hat{\mathcal{G}}_t$ (BH at level $q=0.05$ over all $n^2\hat{p}_t$
entries) that are false, and the share of true edges found. An edge is
true if its window-averaged coefficient is nonzero, since this
average is what the estimates target (Theorem~\ref{thm:dml-var}).
\end{enumerate}
\vspace{-1mm}
\begin{table}[t]
\centering
\caption{Synthetic benchmark results (All confounders observed unless
noted). ET: Extra Trees; LGBM: LightGBM; TP: TabPFN. FDR: BH at
$q=0.05$. Bold: best per column; lower is better except for Power.}\label{tab:toy-results}
\small
\setlength{\tabcolsep}{1.7pt}
\begin{tabular}{@{}lccccc@{}}
\toprule
\textbf{Method} &
\textbf{Lag} &
\textbf{Coeff.} &
\textbf{Fcst.} &
\textbf{FDR} &
\textbf{Power} \\
&
\textbf{RMSE} &
\textbf{MAE} &
\textbf{MAE} &
&
$\uparrow$ \\
\midrule
\multicolumn{6}{@{}l}{\emph{No confounders}} \\
VAR
& 1.475 & 0.060 & 0.110 & 0.129 & 0.655 \\
ACLE-VAR
& 1.006 & 0.066 & 0.109 & 0.147 & 0.661 \\
\addlinespace
\multicolumn{6}{@{}l}{\emph{All confounders, no DML}} \\
VARX
& 1.239 & 0.063 & 0.109 & 0.076 & 0.650 \\
ACLE-VARX
& 1.099 & 0.069 & 0.108 & 0.080 & 0.651 \\
\addlinespace
\multicolumn{6}{@{}l}{\emph{All confounders, causal discovery}} \\
PCMCI + OLS
& 1.214 & 0.065 & 0.108 & \textbf{0.045} & \textbf{0.685} \\
VAR-LiNGAM
& 1.132 & 0.070 & 0.107 & 0.187 & 0.641 \\
\addlinespace
\multicolumn{6}{@{}l}{\emph{All confounders, DML}} \\
OR-VARX, ET
& 1.212 & 0.063 & 0.095 & 0.094 & 0.649 \\
OR-VARX, LGBM
& 1.288 & 0.067 & 0.102 & \textbf{0.045} & 0.666 \\
OR-VARX, TP
& 1.159 & 0.067 & \textbf{0.091} & 0.056 & 0.660 \\
\addlinespace
ORACLE-VARX, ET
& 1.073 & 0.063 & 0.094 & 0.101 & 0.652 \\
ORACLE-VARX, LGBM
& \textbf{0.962} & \textbf{0.058} & 0.101 & 0.047 & 0.673 \\
ORACLE-VARX, TP
& 0.981 & 0.060 & \textbf{0.091} & 0.066 & 0.673 \\
\bottomrule
\end{tabular}
\end{table}
\paragraph{Results.}
\emph{Finding 1: ACLE consistently improves lag tracking.}
Adaptive lag selection reduces lag RMSE regardless of the
base method (Table~\ref{tab:toy-results}). The improvement is 32\% for VAR, 11\% for VARX, and 11\% for OR-VARX with Extra Trees.
It holds for every learner and observability level
(Table~\ref{tab:full-results}; Figures~\ref{fig:edge-traj}--\ref{fig:lag-recovery}).
\emph{Finding 2: Adjusting for confounders lowers the FDR, by an amount that depends on the learner.}
Ignoring the confounders, BH edge discovery overshoots its target of
$0.05$: FDR is 0.129 for VAR. Most false edges are lag-1 self-loops,
which the autocorrelated confounders induce: each series appears to
predict itself.
Adjusting for the confounders in VARX and DML cuts FDR to 0.045--0.080 at similar
power (Table~\ref{tab:toy-results}), except with ET (0.094--0.101). Our theory needs the first-stage fits
to be accurate enough (condition~(N), Appendix~\ref{app:assumptions}). We
suspect that ET, with our hyperparameters, falls short of this: the
confounding it leaves in the residuals inflates the edge tests.
\emph{Finding 3 (tentative): TabPFN is a strong first-stage learner without tuning.}
Among the five DML learners at full observability, ORACLE-VARX with TabPFN
has the best forecast MAE (0.091), the best edge power together with LightGBM (both 0.673),
the second-lowest FDR (0.066, after LightGBM's 0.047; the other three
learners are at 0.10--0.12), and the
second-best lag RMSE (0.981, after LightGBM's 0.962;
Table~\ref{tab:full-results}). This validates TabPFN as an effective
``plug-and-play'' nuisance function estimator for DML, requiring no
hyperparameter tuning.
\emph{Finding 4: PCMCI finds edges as well as ours, but tracks lags and forecasts worse.}
PCMCI has the lowest FDR (0.045) and highest power (0.685), close to
ORACLE-VARX with LGBM (0.047, 0.673). ORACLE-VARX (LGBM, TP) is better
on lag RMSE (0.96--0.98 vs 1.21), coefficient MAE (0.058--0.060 vs
0.065) and forecast MAE (0.091--0.101 vs 0.108). This PCMCI uses a
linear test. We did not run its nonlinear version, which suits nonlinear
confounding like ours: it would take about 30 hours per configuration on
9 cores, against about 1 minute for the linear one
(Appendix~\ref{app:compute}).
\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{plots/toy_oraclevarx_extra_trees_all_edge_trajectories.png}
\caption{Estimated vs.\ true edge coefficients over time for
ORACLE-VARX (Extra Trees, all confounders). Each panel shows one
edge; the black line is the true coefficient $\mathbf{A}_k^*[i,j](t)$ and
the colored line is the rolling estimate $\hat{\mathbf{A}}_k[i,j](t)$.
Vertical dashed lines mark regime boundaries ($t = 1000$, $t = 2000$).}\label{fig:edge-traj}
\end{figure}
\begin{figure}[t]
\centering
\begin{minipage}[t]{0.48\linewidth}
\centering
\includegraphics[width=\linewidth]{plots/toy_aclevar_none_lag_analysis.png}
\end{minipage}\hfill
\begin{minipage}[t]{0.48\linewidth}
\centering
\includegraphics[width=\linewidth]{plots/toy_oraclevarx_extra_trees_all_lag_analysis.png}
\end{minipage}
\caption{Lag tracking: estimated $\hat{p}_t$ (blue) vs.\ true $p^*(t)$
(black step function) over time. \emph{Left}: ACLE-VAR (no confounders).
\emph{Right}: ORACLE-VARX (Extra Trees, all confounders). The lag RMSE of every method is in
Table~\ref{tab:full-results}.}\label{fig:lag-recovery}
\end{figure}
\vspace{-2mm}
\subsection{Financial ETF Application}\label{sec:etf}
\vspace{-2mm}
\paragraph{Data.}
We apply ORACLE-VARX to daily open-to-close returns of nine U.S.\ sector ETFs
(XLY, XLP, XLE, XLF, XLV, XLI, XLB, XLK, XLU), covering January 2000
to December 2020 ($T \approx 5{,}200$ trading days). Ten macroeconomic
confounders serve as exogenous variables: VIX, the Federal Funds Rate,
5-Year inflation expectations, WTI crude oil, the Economic Policy
Uncertainty index, BBB corporate spreads, 10-Year real yield,
the broad U.S.\ Dollar Index, the emerging market Dollar Index, and
gold volatility. The model is estimated with a rolling window of
504~days ($\approx 2$~years) and $p_{\max} = 10$; see
Table~\ref{tab:config} in Appendix~\ref{app:implementation} for full details.
\vspace{-2mm}
\paragraph{Causal structure discovery.}
Unlike the synthetic benchmark, no ground truth is available for the ETF
application. We therefore focus on qualitative properties of the discovered
causal structure.
\emph{The detected lag order tracks market regimes.}
Figure~\ref{fig:etf-lag} plots the selected lag order $\hat{p}_t$
against the 21-day realized volatility of the S\&P~500 (SPY), the
standard deviation of its daily returns over the past 21 days. In calm markets (e.g., 2004--2006, 2013--2014),
$\hat{p}_t \approx 1$: lead--lag relationships are weak and
short-lived. During the 2008 financial crisis and the COVID-19
shock in early 2020, $\hat{p}_t$ spikes to 3--5: the lag structure
elongates as shocks propagate through the financial system, consistent with the
well-documented increase in cross-market equity correlations in bear
markets~\citep{longin2001extreme}.
\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{plots/real_oraclevarx_all10_lag_analysis_with_realized_volatility.png}
\caption{ETF application (ORACLE-VARX, Extra Trees, \texttt{all10}
confounders). Blue step function: adaptively selected lag order
$\hat{p}_t$ over time (2004--2020). Orange dashed line: 21-day
realized volatility of the S\&P~500 (right axis). The lag order spikes
during the 2008 financial crisis and COVID-2020, rising with the
volatility.}\label{fig:etf-lag}
\end{figure}
\vspace{-1mm}
\paragraph{Portfolio performance.}
Table~\ref{tab:etf-sharpe} reports Sharpe ratios of a long-short
portfolio weighted by the forecasts (Appendix~\ref{app:etf-trading}).
ACLE generally raises the Sharpe ratio, and the best portfolio comes from
ACLE-VAR, with neither DML nor confounders. On these data, causal
sufficiency is badly violated: each sector ETF pools many stocks driven by
many common factors, most of them unobserved, and, as in Figure~\ref{fig:causal-suff}b, they create lag links that are
not causal but still predict. ACLE tunes its lag order on validation
forecasts, so it keeps these links with a short lag order. DML, in turn,
no longer yields causal coefficients once sufficiency fails, and its
first stage adds parameters to fit from 504 noisy days.
\begin{table}[t]
\centering
\caption{ETF Sharpe ratios (weighted strategy, 2004--2020) by
confounder set; VAR and ACLE-VAR use none.}\label{tab:etf-sharpe}
\small
\begin{tabular}{@{}lccc@{}}
\toprule
\textbf{Method} & \texttt{vix} & \texttt{macro5} & \texttt{all10} \\
\midrule
VAR & \multicolumn{3}{c}{1.08} \\
ACLE-VAR & \multicolumn{3}{c}{\textbf{1.32}} \\
\addlinespace
VARX & 1.11 & 0.43 & 0.75 \\
ACLE-VARX & 1.07 & 0.81 & 0.86 \\
\addlinespace
OR-VARX, ET & 0.71 & 0.71 & 0.67 \\
ORACLE-VARX, ET & 0.83 & 0.85 & 0.68 \\
ORACLE-VARX, TP & 0.70 & 0.74 & 0.36 \\
\bottomrule
\end{tabular}
\end{table}
\vspace{-1mm}
\Needspace*{6\baselineskip}
\section{Conclusion}\label{sec:conclusion}
\vspace{-1mm}
We introduced ORACLE-VARX, a framework for causal lag discovery in
multivariate time series whose confounders act nonlinearly and whose lag
structure drifts over time. It combines (i) double machine learning, which
removes the confounders' nonlinear effects before the VAR coefficients are
estimated, (ii) ACLE, which picks the lag order in each rolling window by
sequential significance testing tuned on validation data, and (iii)
Benjamini--Hochberg edge selection. We proved that in each rolling window
the debiased coefficients are asymptotically normal around a
window-averaged target, with valid standard errors, even though the
coefficients drift; under causal sufficiency, a zero coefficient means
Granger non-causality. We validated it on synthetic data with known
ground truth and on nine U.S.\ sector ETFs.
\vspace{-2mm}
\section{Future Work}\label{sec:futurework}
Three directions stand out. First, a formal consistency result for
ACLE, covering the BH step and the validation-chosen level
(Section~\ref{sec:theory}). Second, weakening causal sufficiency: hidden
confounders turn lag coefficients into predictive rather than causal
quantities (Figure~\ref{fig:causal-suff}), so richer, high-dimensional
confounder sets, or tests that flag hidden confounding, would let more of
the discovered edges be read causally. Third, extending the model to a structural VAR (SVAR) with
contemporaneous effects~\citep{lutkepohl2005new} could distinguish direct
lagged effects from effects transmitted through same-time links, at the
cost of extra identifying assumptions.
\clearpage
\subsection*{AI use statement}
In this work, we used generative AI tools to state the main theorem and its assumptions precisely and to write its proof (Appendix~\ref{app:proofs}), to implement the method and the experiments, and to give feedback on the experimental design and its engineering, such as reusing repeated model fits.
We have not used generative AI tools to propose the research question or the method, to generate the synthetic data (the authors wrote the functions that simulate it), or to clean the financial data (done before this project used coding agents), and translation and qualitative data analysis are not applicable to this work.
We asked AI tools to comment on the results, but their comments listed details without picking out the main story, so the interpretations in the paper are our own.
Additionally, we used generative AI tools to draft parts of the paper, which the authors then heavily rewrote, to create Figures~\ref{fig:running-dag}, \ref{fig:acle-dag} and~\ref{fig:causal-suff}, to edit the text for readability, and to find related work and check every reference against its source.
We have reviewed all AI-assisted work.
One author, trained in mathematical analysis, checked every step and bound of the proofs by hand; Appendix~\ref{app:ai-use} describes this workflow with an example. The authors read the main mathematical and statistical functions of the code and checked that the outputs are correct or as expected.
We take responsibility for the final content of this work, including text, claims, or artifacts produced with the aid of generative AI.
\bibliography{references}
\clearpage