Papers

Stepped wedge design (SWD) trials typically assume that treatment effects are immediate following implementation, but this assumption is often unrealistic in pragmatic settings where effects evolve over exposure time. This impact is further complicated in multiple-intervention stepped wedge designs (M-SWDs), due to their complex crossover patterns. Existing work has addressed time-varying treatment effects for main effects at the analysis stage using nonparametric approaches; however, the implications for design-stage power and estimand interpretation in M-SWDs remain largely unaddressed. We develop a unified framework that incorporates parametric exposure-response functions into a modified fixed-effects design matrix, Z*, and derive corresponding generalized least squares (GLS) expressions for treatment effect estimands, variance, power, and bias under model misspecification. Analytic and simulation results demonstrate that time-varying effects can require multiple-fold increase in the number of clusters required to achieve the nominal target power. Misspecification can also induce large bias in treatment effect estimates, with particularly pronounced and non-uniform impacts on interaction effects. These findings demonstrate that treatment effect estimands in M-SWDs depend critically on assumptions about exposure-time dynamics and that failure to account for time-varying effects can lead to substantial miscalibration of power and biased inference, affecting main and interaction effects differently.

Learning causal graphs from interventional data is a challenging problem with broad applications. In molecular biology, for example, a central goal is to uncover gene regulatory networks from large-scale perturbation data. An ideal algorithm for this task should scale to thousands of nodes, incorporate interventions even when their targets are unknown, quantify uncertainty, and provide identifiability guarantees. However, existing approaches—e.g. approaches using score-based optimization or approximate Bayesian inference—often fail to meet all of these criteria. To address these limitations, we develop Amortized Bayesian Causal Discovery of Extended Factor Graphs (ABCDEFG). Our method guarantees exact acyclicity, scales to graphs with thousands of nodes, and naturally handles interventions even when their targets are unknown. Additionally, ABCDEFG estimates a posterior distribution whose maximum a posteriori estimate provably identifies the true causal graph up to an equivalence class. On simulated datasets, ABCDEFG achieves state-of-the-art accuracy, producing a well-calibrated posterior distribution while outperforming previous score-based and approximate Bayesian methods. Applied to large-scale single-cell perturbation data, ABCDEFG identifies both established and novel gene targets of growth factors.

Bayesian group sequential designs (GSDs) extend frequentist GSDs with interpretable decision-making and external evidence borrowing, but their use is limited by the computational burden of design-stage operating-characteristic evaluation. Conventional methods simulate virtual trials with Markov chain Monte Carlo or approximate analytical posterior updates at each interim look, making joint calibration of decision thresholds and design skeletons impractical on commodity hardware. Here we introduce a semi-simulation framework with two innovations. First, finite conjugate-mixture priors replace the posterior computation for each look (“per-look”) with closed-form conjugate updates and low-dimensional numerical integration for decision-rule tail probabilities. Second, a precomputation strategy caches per-look posterior tail probabilities from a single Monte Carlo pass at the union of all candidate analysis times, and each design in the calibration grid is evaluated against the same cache by a sub-second sweep, with no further simulation cost. The framework supports posterior-probability decision rules with multiple efficacy and futility criteria under either binding or non-binding futility, and derives closed-form per-look updates for binary, continuous, count and time-to-event endpoints, with benchmarking here focused on the binary endpoint. When applied to re-design the ADRENAL trial, with up to nine analyses, the framework reproduces the operating characteristics of BATSS and adaptr within Monte Carlo error while running, per GSD, approximately $7\times$ to $16\times$ faster than adaptr at a matched budget (several hundredfold at the million-trial calibration budget) and $3{,}700\times$ to $6{,}600\times$ faster than BATSS. This brings routine Bayesian GSD calibration within computational reach for confirmatory trials.

Causal forests that estimate conditional average treatment effects by averaging honest leaf-level effects across trees are widely used in fixed-effects panel settings. We show that this averaging systematically attenuates the estimated heterogeneity: the raw prediction behaves like a + b*tau(x) with slope b < 1, so the spread of the CATEs is compressed toward the average effect, and the additive recentering used to report an unbiased average treatment effect does not fix it. Benchmarking against a similarity-weight generalized random forest on the same within-transformed signal, we find both estimators attenuate but the leaf-averaging construction attenuates materially more. We characterize how b moves with the design, worsening with lower signal-to-noise, smaller panels, and higher dimension; this diagnosis is our main contribution. As a remedy we adapt the best-linear-predictor calibration of Chernozhukov et al., estimating the de-attenuation slope out-of-bag so that it is self-contained within the observational panel and asymptotically inert under a homogeneous effect. In simulations the correction cuts CATE mean-squared error by 25-42% relative to the recentering default; on a standard county minimum-wage panel the attenuation is present but mild and the correction restores the imposed spread. We ship the method in the causalfe Python package.

Given a strongly convex function $u$, equip $R^d$ with a Riemannian metric given by the Hessian $\nabla^2 u$. This is a so-called Hessian manifold. Given a probability density $\mu$ one may run a Langevin diffusion intrinsic to the manifold with stationary distribution $\mu$. Such (Hessian) manifold-valued Langevin diffusions are called Mirror Langevin diffusions (MLD) which have recently become popular. One of the questions we explore is whether, given $\mu$, one can choose $u$ to get an exponential convergence to equilibrium for the MLD, especially if $\mu$ is not strongly log-concave. Our results are based on Lyapunov function methods and give sufficient conditions for a Poincaré or a log-Sobolev inequality to hold for the MLD. These, in turn, imply exponential convergence. We also introduce a Markov chain approximation to the MLD given by a two step Gibbs sampler with stationary distribution $\mu$. This Markov chain is a variant of the Sinkhorn Markov chain introduced in arXiv:2307.16421 that is conjectured to converge to a time-inhomogeneous generalization of the MLD. Under suitable assumptions, we prove that the Markov chain has a guaranteed convergence rate in $\chi^2$ that is consistent with the diffusion time scale. Our proofs are based on ideas from entropic optimal transport and strong data processing inequalities.

Learning the natural parameters $z \in \mathbb{R}^n$ of discrete distributions $\mu_z$ from independent samples constrained to a subset $S \subseteq \{0,1\}^n$ is a foundational challenge in high-dimensional statistics. Existing methods for efficiently estimating truncated Boolean product distributions, notably the work of [Fotakis et al' COLT'20, Algorithmica '22], require either strong local connectivity assumptions on $S$ – a property denoted fatness – or stringent anti-concentration assumptions and necessitate the total mass of the truncation set to be a constant with respect to $n$. Moreover, the results in [Fotakis et al' COLT'20, Algorithmica '22] suffer from sample complexities that scale as $\Omega(2^n)$ if the mass of $S$ is exponentially small in $n$. In this work, we circumvent these limitations by analyzing the geometry of $S$ under the measure $\mu_z$. We refine the existing parameter estimation guarantees under the fatness assumption, improving the prior sample complexity to $O( \log n / \epsilon^2)$ for $\ell_\infty$-recovery, matching the untruncated minimax rate. We further generalize fatness using the notion of influence utilized in the analysis of Boolean functions and provide sufficient conditions for efficient inference. Notably, unlike previous work, our method does not require sampling at arbitrary parameterizations of the model. Lastly, we establish a theoretical lower bound demonstrating the sample complexity exhibits an intrinsic exponential dependence on the width of the model and the minimum distance between elements in the set.

Hybrid mechanistic/data-driven models, which combine first-principles with learned components, are increasingly used in process engineering and scientific machine learning. Common hybrid modeling designs are specified primarily through their architectures and training losses, which offers a limited basis for a shared semantic interface to compare or verify them across domains, with comparatively little attention paid to epistemic uncertainty in the mechanistic part. We bridge hybrid modeling and neuro-symbolic (NeSy) AI by reconstructing these designs as instances of NeSy interface. The resulting translation, Hybrid-to-NeSy (H2N), places mechanistic knowledge on the language side, learned modules on the belief side, and validity domains together with constraints on the logic side. For each design, H2N then yields an explicit NeSy inference functional and a logic-belief decomposition. From this decomposition we derive two metrics: structural violation rate (SVR), measuring whether the learned belief respects the mechanistic structure; and belief dispersion (BD), measuring how concentrated the learned plausibility is, serving as a hybrid model's epistemic uncertainty in its mechanistic part. We instantiate H2N on a case study of a structured hybrid model for binary classification under label noise and show that models with higher SVR and BD exhibit greater variability in held-out accuracy. Under structural distribution shift, H2N further quantifies a model's uncertainty during extrapolations, whereas test accuracy reveals the same shift only post hoc.

AC optimal power flow determines the minimum-cost generation dispatch under nonlinear power balance constraints and is solved thousands of times daily in electricity market operations. Learning a direct mapping from load conditions to OPF solutions can accelerate this computation, yet with deepening renewable penetration, a single optimal dispatch is no longer sufficient. Operators require a characterization of the distribution of feasible near-optimal solutions for risk quantification, sensitivity analysis, and multi-objective trade-off assessment. Supervised neural networks provide fast point predictions but cannot capture this conditional distribution. Diffusion-based generative models can sample diverse solutions in principle, yet existing methods operating in the raw state space exhibit degraded solution quality and fail to scale beyond medium-sized systems. We identify the root cause as the conflation of two distinct tasks within a single model. Compressing the high-dimensional OPF solution manifold is one task, and learning the conditional mapping from loads to that manifold is another. This paper presents FMOPF, a framework that resolves this conflation by decoupling compression from generation through latent flow matching and by explicitly modeling load-state coupling through a Constraint-Aware Interaction Prior Network. Experiments on four IEEE test systems demonstrate that FMOPF provides the most effective Newton-Raphson warm starts, achieves the lowest tail risk among generative methods, and is the first such method to scale to systems with several hundred buses while preserving full feasibility. Ablation studies confirm that the latent generation pipeline is a necessary condition for physical feasibility and that the interaction prior functions as a late-stage tail-risk controller.

In this work, we explore how the inference time of a Transformer Neural Network can be efficiently optimized with applications to real-time anomaly detection in financial time series. The financial time series are price series such as asset prices. Unfortunately, the data is often with errors or outliers that make the downstream data processing tasks useless, unstable or even harmful. Moreover, the amount of financial time-series data has been significantly increasing. Hence, there is a need for better data-cleaning methods in terms of accuracy and in terms of processing speed. Transformers as a neural network architecture have achieved superior performances in many tasks such as Natural Language Processing and Computer Vision. Time series modelling and especially anomaly detection tasks can benefit from the features of transformers architecture in multiple ways, including the capacity to capture long-range dependencies and interactions. Increasingly powerful hardware, such as field-programmable gate arrays (FPGAs), have seen increasing usage in recent years due to their reconfigurability and high performance. They can be efficiently utilized to speed up the computations of the Transformer architecture. We explore different Transformer architectures for time series modelling and how they can be efficiently implemented on an FPGA board (PYNQ-Z2). In particular, we examine the application of Transformers to detect anomalies in time series and we show how they can be efficiently implemented on an FPGA board to minimize latency. The code is available at https://github.com/thxi/icl_thesis

Automating theoretical research requires generating candidate results and evaluating them reliably. Models keep getting better at the first, while the second remains hard. A common approach asks one large language model (LLM) to review what another produced, yet such reviewers are empirically unreliable: they may accept fabricated papers and catch the fabrication at close to chance rates~\citep{badscientist2025}. We present CausalSmith, a framework for automated theoretical research in causal inference built on the Lean proof assistant, where a proof is checked by a program rather than read by a referee. CausalSmith rests on Causalean, a foundational Lean library for causal inference holding 8,179 machine-checked definitions and theorems, developed with language-model assistance under human design and review. Around it, we build a self-improving agentic pipeline that selects research topics, proposes results, formalizes statements, constructs proofs, and presents the resulting artifacts for human inspection. Moreover, the pipeline pairs Lean verification with a statement audit that compares each formal theorem against the informal claim behind it. We evaluate the system using artifacts produced by completed autonomous research runs. The source code, formal library, and run records are available at https://github.com/Jiyuan-Tan/CausalSmith.

For decades, the bootstrap has been a default tool for statistical inference because of its broad applicability and minimal analytic requirements. Although its validity is well understood for smooth parametric estimators, its theoretical properties for many modern semiparametric and machine-learning estimators remain largely unstudied. Nevertheless, bootstrap procedures are often used routinely in such settings, even when their validity is unknown and their computational cost is substantial. We develop the $V$-fold jackknife as a computationally efficient and theoretically justified alternative for semiparametric inference. It requires only $V$ leave-fold-out refits and uses the empirical dispersion of jackknife pseudo-values to quantify uncertainty, without deriving or evaluating an influence function. For regular asymptotically linear estimators of pathwise differentiable parameters, we show that, for fixed $V$, the Studentized $V$-fold jackknife statistic converges to a $t$-distribution with $V-1$ degrees of freedom, giving valid confidence intervals even though the jackknife variance estimator does not converge in probability. When $V\to\infty$, we establish consistency of the variance estimator at rate $V^{-1/2}$, allowing $V$ to diverge slowly, for example at rate $\log n$. We also develop simultaneous confidence bands based on the correct componentwise-Studentized limiting distribution. Finally, we extend the theory to generalized asymptotically linear estimators with diverging influence-function variance and slower-than-$\sqrt n$ convergence; scale invariance of Studentization eliminates the need to know the effective convergence rate. Simulations on the average treatment effect, Kaplan–Meier survival curve, and highly adaptive lasso dose-response curves confirm reliable inference, including where influence-function-based standard errors are anti-conservative or unstable.

In overparameterized linear regression, many weak spectral directions act like a ridge penalty on the signal-bearing spectrum; negative ridge is the natural correction, pushing filters above one. The stable negative-ridge endpoint, however, is structurally limited: its pole must stay below the smallest nonzero empirical eigenvalue, and it anti-shrinks smaller eigenvalues more than larger ones. Early-stopped negative-shifted gradient descent escapes this constraint. Its filter is smooth at the would-be pole and mixed-sign-capable: above-ridgeless directions form a leading prefix, with lower directions shrunk or exposure-controlled while stopping sets the crossover. In a Gaussian spike-plus-flat model we discover a Marchenko-Pastur barrier: the shift that cancels the implicit penalty lies a bulk width above the smallest empirical eigenvalue, and the stopped path improves on every admissible endpoint by a polynomial factor in risk under explicit conditions. Our main theorem permits a general high-effective-rank tail: its trace sets the implicit floor, its squared spectrum controls exposure, and the floor-critical path recovers all head scales at once, beyond positive shrinkage and, once scales separate, every uniform rescaling of ridgeless. Handling the noncontractive shifted dynamics is the central technical challenge; localized Duhamel integrals control them. A finite-grid hold-out inequality transfers the separations to the validation-selected algorithm.

This study evaluates the out-of-sample forecasting ability of six model types: ARIMA, VAR, Random Forest, XGBoost, LSTM, and GRU, for monthly Canadian inflation from January 2012 to April 2026 (n = 172). The evaluation employs expanding-window walk-forward validation across 1-, 3-, 6-, and 12-month horizons. Results reveal a horizon-dependent shift: ARIMA significantly outperforms all machine learning and deep learning models at the one-month horizon (Diebold-Mariano p < 0.05). However, Random Forest and XGBoost become notably superior at six and twelve months, reducing RMSE by approximately 30-75 percent compared to ARIMA and VAR. LSTM and GRU perform well only at the shortest horizon, likely due to overfitting given the limited data. Analyzing four macroeconomic sub-periods shows that no single model consistently dominates. SHAP analysis of the top-performing XGBoost model indicates that lagged inflation is more influential than unemployment, which only becomes significantly impactful during the pandemic tail. The findings clarify when machine learning methods can surpass traditional benchmarks.

Differential privacy is the standard for ensuring data privacy and is widely used in major data publications, including reporting results from the U.S. decennial census. Common implementation of differential privacy uses independent Gaussian or Laplace noise addition to the database. However, there could be aggregate (linear) queries to the database that are excluded from the privacy budget, for example, state totals that can not be perturbed due to constitutional mandates. Any implementation of a differential privacy is required to honor these constraints, also referred to as invariants. Under aggregation constraints, the noise vector is no longer independent and the traditional differential privacy guarantees have to be re-evaluated. We propose a high entropy differential privacy implementation that maintains the aggregation invariants with probability one or exponentially close to one and derive the privacy guarantees for the implementation under the invariants. The theoretical proof covers a partial solution to an open question about the null space of correlation matrices. Moreover, the methodology has general use in the context of sampling from normal mixture models under linear equality constraints.

Piecewise Deterministic Monte Carlo (PDMC) algorithms utilize continuous time Markov processes to generate samples from continuous distributions, and provide a modern alternative to discrete time Markov chain Monte Carlo algorithms. In this work we provide a survey of recent results on scaling limit arguments to understand the behaviour and efficiency of two often-used Piecewise Deterministic Monte Carlo algorithms: the Bouncy Particle Sampler and the Zig-Zag Process. In particular we discuss a Functional Central Limit Theorem, scaling in the high-dimensional regime, scaling under anisotropy, and scaling in the big data regime. This work is intended as part of the proceedings of the 2024 Isaac Newton Institute Programme "Stochastic systems for anomalous diffusion".

Data analysis agents are becoming increasingly common tools for applied and scientific research. Yet, for highly specialized tasks such as Observational Causal Inference (OCI), human oversight remains necessary to ensure the validity of results. We introduce `oci-agent`, an open-source Python package that implements a human-in-the-loop agentic workflow for observational causal inference. `oci-agent` is designed to automate vital but laborious aspects of applied causal inference, such as covariate balance checking, propensity score trimming, and sensitivity analysis, so that humans can focus on more nuanced tasks, such as framing questions, scrutinizing assumptions, and evaluating diagnostics and results. We initially open-sourced `oci-agent` in June 2026 with support for doubly robust learning of the average treatment effect of a single binary treatment. Since then, we have added support for heterogeneous treatment effect estimation and for multiple continuous treatments via partially linear models. In this paper, we describe the principles behind `oci-agent` and offer internal Netflix case studies and evaluations on public data of its capabilities. Across numerous evaluations, `oci-agent` outperforms less structured baselines while remaining competitive with hand-tuned benchmarks. `oci-agent` is used extensively for causal inference at Netflix and has orchestrated more than 100 analyses per month since its release in June.

As data collecting technologies advance, data structures are getting more and more complex, from single vectors to multi-dimensional tensors. This article is motivated by a variable selection problem to detect important genes from an ultrahigh dimensional pool that are associated with human facial shape variations. We propose a data-driven trimmed feature screening method based on a tensor ridge regression model (TrimTenRidge) through setting thresholds on the tensor coefficients to perform a feature screening procedure. Unlike existing approaches, the TrimTenRidge does not require any sparse structures. In addition, it not only detects important predictors but also locates specific regions/components of the tensor response that are associated with each of the selected predictors. We prove the theoretical selection consistency and also assess its empirical performance through various simulation settings. The approach copes with ultra-high dimensional predictors and tensor responses simultaneously and contributes to the literature from theoretical, methodological, and five applicational aspects. We further apply the TrimTenRidge approach to genome-wide human facial shape data, from which the entire facial shapes form a $2,342\times 7,160\times 3$ tensor, and we successfully detect several novel genetic loci and also confirm some existing findings that are associated to facial shape.

This work addresses the generation of theoretical correlation matrices with prescribed sparsity patterns associated to graph structures. We propose a novel convex optimization framework in which an initial matrix is projected onto an elliptope under a positive semidefiniteness constraint. Several numerical schemes are implemented and compared. The problem falls within the broader class of matrix completion, where off-diagonal entries corresponding to absent edges are fixed to zero and diagonal entries are fixed to one. Beyond this structural constraint, the approach offers greater flexibility than existing methods by allowing control over the mean of the off-diagonal entry distribution, enabling the generation of correlation matrices that better reflect realistic data. This procedure is not designed to yield a uniform distribution over the feasible set; rather, it provides a principled and tunable way to construct correlation matrices suitable for benchmarking statistical methods for graphical model inference. Theoretical guarantees on the existence of solutions are established, both in the general setting and under the additional mean constraint. Simulation studies illustrate the properties of the generated matrices with respect to graph structure. The methodology is applied to two real-world datasets from neuroscience and finance, and a comparison with GAN-based correlation matrix generation is provided.

We consider the problem of learning from a single finite trajectory of an ergodic stochastic dynamical system. More precisely, we study discrete-time autonomous stochastic systems defining time-homogeneous Markov processes. We first focus on estimating the optimal one-step prediction function by nonlinear least squares, and derive high-probability guarantees measured with respect to the invariant measure of the process. These results make explicit how the non-independent and non-identically distributed nature of trajectory data modifies the classical statistical learning analysis. We then extend the framework to higher-order systems and finite-state spaces. Finally, we show that the same least squares and concentration arguments naturally extend to learning Koopman operators. Our approach combines tools from statistical learning theory and quantitative ergodic theory for Markov chains. It relies, in particular, on a concentration inequality for Hilbert-space-valued additive functionals of uniformly geometrically ergodic Markov chains.

Curvature notions on graphs, particularly Ollivier-Ricci and Forman, have emerged as powerful tools for addressing fundamental issues in Graph Neural Networks (GNNs) such as oversmoothing and oversquashing, but rely almost exclusively on local edge-level comparisons and therefore fail to certify how information actually propagates over long distances. We introduce Entropic Curvature, a global, transport-based curvature obtained by extending the Lott-Sturm-Villani framework to graphs through the displacement convexity of entropy along Wasserstein geodesics. We define a tractable Weak Entropic Curvature proxy that lower-bounds the global entropic curvature, and from it derive (i) a Poincare-type inequality controlling oversmoothing, (ii) a transport-entropy generalization bound, and (iii) an expansion paradox proving that sparsity, strong spectral expansion, and positive entropic curvature cannot coexist in large graphs, unifying oversmoothing and oversquashing as opposite ends of a single curvature spectrum. We translate the theory into three practical mechanisms, the E-Gate aggregator, the ENT structural encoding, and Midpoint-Completion Rewiring (MCR), and benchmark them against SDRF, FoSR, BORF, LCP, and Graph Ricci Flow on six node-classification benchmarks, and graph-classification.

Spatial intensity maps are routinely used to summarize point patterns on geographically constrained regions, such as islands, coastlines, watersheds, ecological reserves, and administrative areas with physical barriers. In these settings, the domain is not a nuisance feature– it determines where probability mass may be assigned and which locations should be smoothed together. Standard kernel density estimators can place mass outside the study region and smooth according to Euclidean distance, while diffusion-based estimators respect the geometry but are costly to recompute when the data subset, bandwidth, or evaluation grid changes. Motivated by repeated hotspot mapping of theft and larceny incidents on Oahu, Hawaii, we propose the projected diffusion kernel density estimator (PDKDE). PDKDE expands the Neumann diffusion kernel in a truncated Laplacian eigenbasis, so that the geometry of a fixed domain is computed once and subsequent density estimates are obtained through explicit spectral coefficients. The resulting estimator preserves the boundary-aware and geometry-respecting behavior of diffusion smoothing while making repeated estimation and least-squares cross-validation computationally practical. For the exact projected estimator, we prove MISE consistency and pointwise consistency up to the boundary. Controlled simulations show that PDKDE reduces boundary and barrier artifacts relative to Euclidean kernels and is substantially faster than the geometry-aware comparators considered after the one-time domain computation. In the Oahu application, the method produces coastline-constrained descriptive maps across time windows in seconds, illustrating its intended fixed-domain, changing-data use case.

Ill-posed inverse problems are encountered in numerous applications, possibly characterized by a highly non-linear forward model, both additive and multiplicative sources of noise, and censored data. In the absence of ground truth, uncertainty quantification is crucial to assess estimation reliability. This motivates the use of a Bayesian model and stochastic inference methods such as Markov Chain Monte Carlo algorithms. Problems combining all these challenges often lead to a complex and potentially multimodal posterior distribution, difficult to handle in practice. Approximate approaches have been proposed in the literature by either neglecting a source of noise or by using a tractable approximation of the likelihood function. These approaches either lead to an inaccurate model, or may require a complex calibration of the approximate likelihood. This paper proposes to tackle such problems with a general hierarchical Bayesian model and an efficient MCMC algorithm. The proposed formulation bypasses the need for calibrating the hyperparameters of an approximate model and is more versatile. The proposed method is assessed on a challenging scenario encountered in astronomy using synthetic data, in a variety of noise and censoring configurations. Comparisons are conducted against two baselines and a state-of-the-art method applicable in this context. The proposed approach is general and yields state-of-the-art results in terms of point-wise estimates and computing costs, with superior predictive performance. Results in the supplementary material further complete this comprehensive and rigorous model study. This work can serve as a guide for practitioners to select the best likelihood model according to their specific application.

Estimating contemporaneous bidirectional interactions from observational data is difficult because each outcome is endogenous to the other, while flexible regressions may capture only reduced-form dependence. This paper proposes SEM-DNN, a heteroscedastic neural simultaneous-equation estimator that learns reciprocal structural interactions without external instruments. Identification exploits conditional covariance diagonalization: when structural shocks have zero conditional means, are conditionally uncorrelated given predetermined covariates, and exhibit nonproportional conditional variances, only the true interaction coefficients diagonalize the conditional residual covariance across the feature space. The method jointly approximates nonlinear structural mean functions and feature-dependent variances using a diagonal Gaussian quasi-likelihood that incorporates the simultaneous-system Jacobian. We establish unique identification and positive-definite local curvature of the profiled population criterion and show that, under neural-profile compatibility conditions, the implemented neural criterion inherits this curvature despite nonunique network parameterizations. The coefficients admit a causal interpretation when the structural equations represent autonomous mechanisms that remain invariant under the relevant interventions. Monte Carlo experiments with nonlinear, high-dimensional nuisance functions and non-Gaussian shocks show that SEM-DNN recovers structural effects more reliably than parametric, kernel-based, and separate-equation neural alternatives as information increases, although at greater computational cost. An application to ready-to-eat cereal scanner data illustrates how the method can study contemporaneous price-sales feedback and assess identification strength, residual diagonalization, variance calibration, and optimization sensitivity.

Day-ahead electricity prices are forecast by many competing models and at several temporal granularities, from hourly prices to block and baseload products. The resulting forecasts suffer from two distinct problems: forecasts produced at different granularities are incoherent, as the aggregates do not match the averages of their components, and no single model is the most accurate in every market, period and level. We develop a temporal coherent combination approach that addresses both problems in one step, pooling the forecasts that the competing experts produce at all the levels of a temporal hierarchy into a single forecast that satisfies the aggregation constraints and has minimum error variance among the linear unbiased combinations. Since the optimal solution depends on a high-dimensional error covariance matrix, we compare different estimators obtained by crossing correlation structures with linear and nonlinear shrinkage approaches. Using the base forecasts of four model classes for the German and Spanish day-ahead markets, the combined forecasts significantly outperform both the base forecasts and the best reconciled expert at nearly every temporal level.

Forecasting multiple time-series with high-dimensional covariates presents a core challenge: unifying common temporal patterns while retaining meaningful series-specific information. We introduce Hopformer (Homogeneity-Pursuit Transformer), a two-stage framework that addresses this challenge. In the first stage, we perform a Sparsity Pattern Aggregation (SPA) scheme extracting a common low-variance trend that incorporates the covariates. This acts as a homogenization layer. In the second stage, a LoRA-fine-tuned Transformer models the remaining complex dependencies in the residual. Our method is theoretically grounded. We prove that SPA achieves a near-optimal bias-variance trade-off via an oracle inequality. We also provide generalization bounds for the second stage under dependent time series data. Hopformer sets a new state of the art, improving MASE by an average of 6.56% across synthetic and real-world forecasting benchmarks.

Anomaly detection is inherently characterised by severe class imbalance, making the interpretation of evaluation metrics challenging. Although metrics such as AUROC, AUPR, F1-score, and MCC are widely used, their values convey different meanings depending on the anomaly ratio. In this work, we analyse the behaviour of those four common anomaly detection metrics under varying levels of imbalance. We focus on the study of metric landscapes, visualisations that relate metric values to true positive and true negative rates, providing an intuitive view of metric preferences and stability. Our analysis offers practical guidance for interpreting and comparing anomaly detection results across datasets with different imbalance ratios.

Pearson correlation is the default measure of association in most statistical software, yet it is only appropriate for pairs of continuous variables with a linear relationship. When variables are binary, ordinal, or categorical, specialized methods (e.g., point-biserial, polychoric, tetrachoric, and Cramér's~$V$) may be more appropriate, but practitioners rarely know which to select. The smartcor package for R (and its companion pysmartcor package for Python) automatically detects variable types, selects the statistically appropriate correlation method for each pair, and explains its reasoning. The package supports 14~correlation and association methods covering all 10~variable-type pair combinations, and distinguishes true correlation (for ordinal and continuous pairs) from statistical association (for nominal categorical pairs). Monte~Carlo simulations validate the selection logic, and a case study with General Social Survey data demonstrates substantive differences between naive and type-aware correlation analysis.

Remaining useful life (RUL) prediction and failure-mode classification are central tasks in predictive maintenance. Many data-driven pipelines use fixed-window supervised learning with complete terminal labels; such routes do not naturally encode the temporal recursion linking successive degradation-state predictions when observations are partial or unit identities are unavailable. We formulate prognostics as vector General Value Function (GVF) prediction on an absorbing degradation process, treating RUL and failure-mode probabilities as temporally consistent targets rather than independent window-level labels, and estimate them with a multi-step temporal-difference estimator, TD($n,\lambda$). Supporting theory identifies the Bellman fixed point of the vector GVFs, characterizes the linear projected-TD limit and its relation to complete-return Monte Carlo regression under realizability, and explains when bootstrapped TD targets are less variable than Monte Carlo returns. On an event-triggered multimode simulation and NASA C-MAPSS label-scarce stitch data, TD improves RUL and failure-mode prediction relative to a supervised same-backbone Monte Carlo control, especially under scarce complete labels. Practically, fragmented, identity-free degradation records can contribute local Bellman transitions instead of being discarded until complete run-to-failure labels are available.

A data-driven inverse optimization problem (DDIOP) is the problem of estimating the objective-function parameters (weights) that explain observed optimal-solution data, and it arises in many applications, including integer linear programming (ILP). It is known that, by applying gradient-based optimization methods to the suboptimality loss, the inverse optimization of ILPs can be solved exactly within finitely many oracle iterations, and that the required number of iterations is bounded as $T=O(1/\gamma(\ell_{\mathrm{sub}})^2)$ in terms of a problem-dependent geometric constant $\gamma(\ell_{\mathrm{sub}})$. However, no means of bounding $\gamma(\ell_{\mathrm{sub}})$ from below as a function of the problem size has been available, and hence the number of iterations could not be given as an explicit function of the problem size. We therefore give, when the forward problem is an integer linear program (ILP), the number of iterations sufficient for projected subgradient descent applied to the suboptimality loss to achieve exact consistency with the observed data, as a fully explicit function of the number of samples, the dimension of the features, the ranges of the features, and the structure of the constraint coefficient matrix, up to polynomial factors in the basic constants (the diameter of the weight set, the step-size parameter, and the Lipschitz constant of the suboptimality loss).

Modeling shared and subject-specific structure in multisubject spatiotemporal data remains challenging, particularly in neuroimaging, where both spatial and temporal patterns exhibit rich variability across subjects. Existing matrix and tensor decompositions provide interpretable factorizations, but rely on fixed multilinear structures or coupling schemes that may limit their flexibility in capturing complex variability. In this work, we introduce a spatiotemporal variational tensor decomposition (ST-VTD) framework that combines a tensor factorization generative model with structured priors to jointly represent spatial maps and temporal dynamics. Spatial factors are regularized to promote a low-rank structure inspired by the LL1 decomposition, while temporal factors are modeled using a learned Long short-term memory (LSTM)-based prior, enabling flexible and adaptive dynamics. Posterior inference is performed using an amortized variational formulation by unrolling iterations of an optimization algorithm, leading to an interpretable and parameter-efficient architecture. The proposed inference framework employs a warm-start strategy based on group independent component analysis, which we found to improve optimization performance. Experiments on a realistic synthetic functional MRI (fMRI) dataset demonstrate that the proposed approach significantly improves latent factor recovery compared with representative classical and probabilistic decomposition benchmarks.