Research
Estimating Bidirectional Causal Effects with Large Scale Online Kernel Learning
Overview Research area: Causal inference in econometrics combined with large-scale machine learning (stat.ML), specifically simultaneous equation models, heteroskedasticity-based identification, kerne
- arXiv
- 2511.05050
- Published
- 2025-11-07
- Authors
- Masahiro Tanaka
AI summary
Overview
Research area: Causal inference in econometrics combined with large-scale machine learning (stat.ML), specifically simultaneous equation models, heteroskedasticity-based identification, kernel methods, and online optimization.
Technical level: Advanced. The paper assumes familiarity with simultaneous equation models, quasi-maximum likelihood, instrumental variables, kernel approximations, and stochastic gradient optimization.
Scope: The paper proposes and evaluates a single framework that jointly estimates bidirectional causal effects and their nonlinear conditional variances using random Fourier features and online gradient descent.
What This Paper Is About
Most machine-learning-based causal inference estimates only one-directional effects, where one variable influences another but not the reverse. Many real systems, however, are mutually dependent — price and quantity, brain activity and behavior, policy and public response — so estimating each direction requires a model of two equations that determine one another. The paper's goal is to estimate both directions of such a relationship flexibly, at scale, and in a single joint model rather than two separate regressions.
Key Contributions
-
A joint estimator for bidirectional causal effects. The method integrates a quasi-maximum likelihood estimator for simultaneous equation models with large-scale online kernel learning, so both structural equations are estimated together in one step.
-
Flexible nonlinear modeling of both means and variances. Random Fourier feature approximations are used to model nonlinear conditional means and nonlinear conditional variances, addressing the problem that domain theory usually gives little guidance on how to specify the variance functions that identification depends on.
-
An efficient online learning algorithm. An adaptive online gradient descent procedure with adaptive gradient clipping and Adam-tuned step sizes allows parameter updates as data stream in, giving near-linear scaling with the number of observations and covariates.
-
Identification conditions and a simulation benchmark. The paper states local identification conditions and reports an extensive simulation study comparing the approach against two baselines across three data-generating processes, two sample sizes, two dimensionalities, and 1000 Monte Carlo replications, plus a sensitivity analysis and computation-time comparison.
Main Findings
-
The proposed method beat both baselines on bias and RMSE. Across the three data-generating processes and both sample sizes, the proposed SEM-Kernel method consistently produced lower bias and root mean squared error for γ₁ and γ₂ (true values fixed to γ₁ = −0.5 and γ₂ = 1.0).
-
The single-equation baseline was systematically biased. Single-Kernel, which ignores the simultaneous-equation structure, showed large, consistent bias — for example, at n = 20,000, DGP-1, d = 100, its γ₁ bias was 0.732 versus −0.007 for SEM-Kernel — confirming that neglecting the endogeneity between the two outcome variables yields inconsistent estimates even with a flexible kernel specification.
-
The polynomial baseline was numerically unstable. SEM-PAB, presented as the most flexible specification from the earlier literature, showed severe bias (for example, γ₁ bias of 1.212 at n = 5000 for DGP-1) and much larger errors across nearly all settings.
-
Results were strongest in low dimension. The paper reports that the improvement was most pronounced at d = 100, which it describes as demonstrating the scalability and robustness of the kernel representation in high-dimensional settings.
-
Performance gains were largest under complex heteroskedasticity. The authors state that gains were particularly strong for DGP-3, the data-generating process with the most complex functional forms, suggesting the random feature approximation captures local smoothness and heterogeneity in conditional variances.
-
A moderate number of Fourier components suffices. In the sensitivity analysis over m ∈ {100, 200, 500, 1000, 2000, 5000} on DGP-3, bias and RMSE dropped rapidly up to roughly m = 500 and then showed no notable improvement, while computation time grew approximately linearly with m. The paper concludes that m = 500–1000 is adequate.
-
Computation time was close to the single-equation baseline. Despite jointly estimating both equations and modeling heteroskedasticity, SEM-Kernel was only slightly slower than Single-Kernel: about 0.7 s at n = 5000 and 4.5 s at n = 20,000 according to the text, with the table reporting, for example, 1.33 s versus 1.10 s at n = 5000, d = 1000.
-
The polynomial baseline was more than an order of magnitude slower. SEM-PAB took 27.89 s and 112.15 s at n = 5000 and n = 20,000 respectively for DGP-1 with d = 1000, versus 1.32 s and 4.53 s for SEM-Kernel.
-
Local identification reduces to a simple condition. Under the paper's assumptions and its kernel specification, the true parameter vector is locally identified if and only if the two conditional-variance coefficient vectors are not equal (α₁ ≠ α₂).
-
Reported limitations. The authors state the method estimates only linear causal effects and assumes symmetry of the error terms.
Methodology in Plain English
The authors start from an existing idea in econometrics called heteroskedasticity-based identification. In ordinary instrumental-variable work, you identify causal slopes using outside shifts in the average level of a variable. Here, identification instead comes from how the variability of one variable changes while the other stays stable. If the noise variance of one equation moves with some observed variables while the other equation's does not, the slopes of both equations become identifiable — and both directions can be estimated at once.
Their model has two equations, each with one dependent variable, a term capturing the causal pull from the other dependent variable (the parameters of interest, γ₁ and γ₂), a flexible function of observed covariates, and an error term whose variance depends on covariates. Both the flexible mean functions and the variance functions are approximated with kernels, which the authors make practical using random Fourier features: instead of computing kernel values directly, they draw random frequencies from the kernel's spectral density and build a finite feature vector of sines and cosines. With a Gaussian kernel, the required spectral density is normal, and the bandwidth is chosen with a median heuristic on pairwise distances.
Because kernel inner products are now just inner products in a modest-dimensional Euclidean space, the parameters become ordinary coefficient vectors, and the whole model becomes a log quasi-likelihood that can be optimized one observation at a time. The authors derive the loss analytically and update parameters online with adaptive gradient clipping — the gradient is rescaled by a ratio of a running average of past gradient norms to the current gradient norm — combined with Adam for the step size. They also state assumptions and local identification conditions for the estimator.
The evaluation is a simulation study. Synthetic data follows three data-generating processes with varying complexity of mean and variance functions; covariates are drawn from a zero-mean multivariate normal distribution with a randomly generated correlation matrix from an inverse Wishart distribution rescaled to a correlation matrix. Errors are generated as chi-squared variables with 10 degrees of freedom, normalized to zero mean and unit variance, then multiplied by the square root of the conditional variance. The method is compared against a version that models the mean with the same kernel but fits each equation separately, and a version that uses polynomial approximations with beta function weights and a Box–Cox transformation. Each configuration is repeated 1000 times and judged by mean bias, standard deviation, and root mean squared error.
Why This Matters
Impact on research. The work connects an econometric identification strategy that is normally presented as a single-sample procedure with the online, kernel-based tooling of modern machine learning. It broadens machine-learning causal inference — which the authors argue has largely focused on unidirectional relationships — to systems where causality runs both ways, and it removes a practical obstacle in the earlier quasi-maximum-likelihood approach by making the conditional variance specification a learnable, kernel-based object rather than a hand-specified one.
Real-world applications (these are motivating examples named in the paper, not applications tested in it):
- Interactions between brain activity and behavior.
- Predator and prey population dynamics.
- Policy interventions and public responses.
- Employee morale and organizational performance.
The paper also frames the contribution as general to natural and social science, policy making, business, and industrial applications.
Industry relevance. The online, streaming formulation and near-linear scaling suit settings where data arrives continuously and covariates are numerous. The authors highlight the near-parity in computation time with the simpler single-equation baseline and the roughly linear growth of the proposed method's runtime with m as evidence that the approach is practical for large-scale, high-dimensional causal inference in real time.
Future Directions
-
Nonlinear causal effects. The paper explicitly names its restriction to linear causal effects as a limitation and identifies extension to nonlinear causal relationships as an important direction.
-
Strongly skewed and asymmetric data. The authors also cite the assumption of symmetry of the error terms as a limitation, and propose accommodating strongly skewed data as future work. This connects to their simulation design, which uses chi-squared errors with 10 degrees of freedom.
-
Choice and calibration of the feature dimension. Although the sensitivity analysis suggests m = 500–1000 is adequate and stable across sample sizes, whether a fixed moderate m remains sufficient for more complex or higher-dimensional settings than those tested is left open.
-
Behavior under other data-generating processes. The simulation covers three data-generating processes, so how the method behaves under alternative heteroskedastic structures, other kernel families, or kernels not built on a Gaussian assumption is not reported in the paper.
Target Audience
Researchers at the intersection of econometrics and machine learning who work on causal inference with endogeneity or simultaneity; statisticians and ML researchers interested in kernel methods and online learning for structural models; applied quantitative analysts in economics, policy, and industry who need to estimate two-way causal relationships on large or streaming datasets; and graduate students with prior exposure to instrumental variables and quasi-maximum likelihood who want to see identification theory paired with a scalable learning algorithm. Readers without that econometric background will find the identification arguments demanding, though the simulation study and computation-time comparison are readable on their own.
Authors’ abstract
In this study, a scalable online kernel learning framework is proposed for estimating bidirectional causal effects in systems characterized by mutual dependence and heteroskedasticity. Traditional causal inference often focuses on unidirectional effects, overlooking the common bidirectional relationships in real-world phenomena. Building on heteroskedasticity-based identification, the proposed method integrates a quasi-maximum likelihood estimator for simultaneous equation models with large scale online kernel learning. It employs random Fourier feature approximations to flexibly model nonlinear conditional means and variances, while an adaptive online gradient descent algorithm ensures computational efficiency for streaming and high-dimensional data. Results from extensive simulations demonstrate that the proposed method achieves superior accuracy and stability than single equation and polynomial approximation baselines, exhibiting lower bias and root mean squared error across various data-generating processes. These results confirm that the proposed approach effectively captures complex bidirectional causal effects with near-linear computational scaling. By combining econometric identification with modern machine learning techniques, the proposed framework offers a practical, scalable, and theoretically grounded solution for large scale causal inference in natural/social science, policy making, business, and industrial applications.