Paper deep dive
Explainable Deep Learning for Price-Trade Dynamics: From Black-Box Forecasts to Effective Parametric Models
Manuel Naviglio, Fabrizio Lillo
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 9/9/2026, 5:44:03 AM
Summary
This paper proposes a methodology to bridge the gap between black-box deep learning models and interpretable parametric models in financial market microstructure. By training a deep feed-forward neural network on high-frequency returns and signed volumes, the authors use Shapley-based explainability to uncover nonlinear dependencies. They identify that lagged signed volume has sign-preserving and saturating effects, while lagged returns act as state variables. These insights are used to construct a parsimonious SHAP-inspired nonlinear parametric model that outperforms linear VAR benchmarks and matches neural network performance while maintaining interpretability.
Entities (8)
Relation Signals (6)
SHAP → usedfor → Explainability
confidence 98% · Using Shapley-based explainability, we show that the dominant contributions are concentrated at the most recent lags.
SHAP-inspired Nonlinear Parametric Model → derivedfrom → Deep Feed-Forward Neural Network
confidence 95% · Building on these findings, we introduce a parsimonious SHAP-inspired nonlinear parametric model. It reproduces the main return-volume dependencies...
SHAP-inspired Nonlinear Parametric Model → outperforms → Vector Autoregression (VAR)
confidence 93% · It reproduces the main return–volume dependencies, outperforms the linear VAR benchmark, and achieves performance comparable to the neural network.
Deep Feed-Forward Neural Network → outperforms → Vector Autoregression (VAR)
confidence 92% · The neural network improves predictive performance, especially for returns, revealing nonlinear dependencies beyond the linear specification.
Signed Volume → exhibits → Nonlinear Price Impact
confidence 90% · Lagged signed volume generates sign-preserving and saturating effects, consistent with nonlinear price impact and order-flow persistence.
Returns → actsas → State Variables
confidence 89% · Lagged returns act as state variables: when the previous trade does not move the price, the model predicts continuation... whereas non-zero returns generate attenuation or reversal.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Understanding the joint dynamics of prices and trades is central to market microstructure, where returns and order flow interact through nonlinear and state-dependent mechanisms. Linear models are interpretable but may miss these effects, while deep neural networks improve forecasting at the cost of transparency. We use neural networks as tools for structural discovery rather than only for prediction. A deep feed-forward network is trained on high-frequency returns and signed volumes for large- and small-tick stocks and compared with a linear VAR benchmark. The neural network improves predictive performance, especially for returns, revealing nonlinear dependencies beyond the linear specification. Using Shapley-based explainability, we show that the dominant contributions are concentrated at the most recent lags. Model-implied responses are consistent with conditional averages reconstructed from the data. Unlike empirical averages, however, the neural-network decomposition isolates individual regressor contributions to the aggregate dependence. Lagged signed volume generates sign-preserving and saturating effects, consistent with nonlinear price impact and order-flow persistence. Lagged returns act as state variables: when the previous trade does not move the price, the model predicts continuation in the direction of past order flow, whereas non-zero returns generate attenuation or reversal. Building on these findings, we introduce a parsimonious SHAP-inspired nonlinear parametric model. It reproduces the main return-volume dependencies, outperforms the linear VAR benchmark, and achieves performance comparable to the neural network. A multi-lag extension captures residual longer-memory effects while preserving interpretability. Overall, explainability offers a route from black-box prediction to economically meaningful parametric models of price and trade dynamics.
Tags
Links
- Source: https://arxiv.org/abs/2609.06085v1
- Canonical: https://arxiv.org/abs/2609.06085v1
Trouble viewing inline? Open PDF directly →
Full Text
154,424 characters extracted from source content.
Expand or collapse full text
Explainable Deep Learning for Price–Trade Dynamics: From Black-Box Forecasts to Effective Parametric Models Manuel Naviglio Fabrizio Lillo Scuola Normale Superiore, Pisa, Italy †thanks: Email: manuel.naviglio@sns.it†thanks: Email: fabrizio.lillo@sns.it Abstract Understanding the joint dynamics of prices and trades is central to market microstructure, where returns and order flow interact through nonlinear and state-dependent mechanisms. Linear models are interpretable but may miss these effects, while deep neural networks improve forecasting at the cost of transparency. We use neural networks as tools for structural discovery rather than only for prediction. A deep feed-forward network is trained on high-frequency returns and signed volumes for large- and small-tick stocks and compared with a linear VAR benchmark. The neural network improves predictive performance, especially for returns, revealing nonlinear dependencies beyond the linear specification. Using Shapley-based explainability, we show that the dominant contributions are concentrated at the most recent lags. Model-implied responses are consistent with conditional averages reconstructed from the data. Unlike empirical averages, however, the neural-network decomposition isolates individual regressor contributions to the aggregate dependence. Lagged signed volume generates sign-preserving and saturating effects, consistent with nonlinear price impact and order-flow persistence. Lagged returns act as state variables: when the previous trade does not move the price, the model predicts continuation in the direction of past order flow, whereas non-zero returns generate attenuation or reversal. Building on these findings, we introduce a parsimonious SHAP-inspired nonlinear parametric model. It reproduces the main return–volume dependencies, outperforms the linear VAR benchmark, and achieves performance comparable to the neural network. A multi-lag extension captures residual longer-memory effects while preserving interpretability. Overall, explainability offers a route from black-box prediction to economically meaningful parametric models of price and trade dynamics. 1 Introduction The dynamics of trading volumes and prices are intrinsically linked through the microstructure of financial markets. In modern electronic markets, this interaction is mediated by the limit order book (LOB), where buy and sell orders are continuously submitted, canceled, and executed. Price changes arise from the interaction between incoming market orders and available liquidity, making order flow a primary driver of short-term price dynamics. This mechanism underlies several well-documented stylized facts, including persistent order flow, transient price impact, and short-term return reversals [1, 2, 3]. Modeling the joint dynamics of prices and trades is therefore central to market microstructure. On the one hand, accurate forecasts of returns and trading activity are essential for applications such as optimal execution, risk management, and algorithmic trading [4]. On the other hand, models that faithfully reproduce the statistical properties of the data are crucial for realistic market simulations, which are widely used to study liquidity provision and market impact. In particular, market impact—the price response induced by trading—plays a fundamental role in determining transaction costs, especially for metaorders executed over time [1, 5]. Capturing these mechanisms requires models that account for the complex and nonlinear interaction between order flow and price formation. A large body of literature has addressed this problem using linear parametric models. A prominent example is the vector autoregression (VAR) framework introduced by Hasbrouck [6], as well as its variants such as the Transient Impact Model [1]. These models describe price and trade dynamics through linear dependencies on past observations and have been widely adopted due to their interpretability and tractability. However, their ability to capture the full complexity of high-frequency financial data is inherently limited, as they impose restrictive assumptions on the functional form of the underlying relationships. More recently, machine learning methods have emerged as flexible alternatives capable of capturing nonlinear and high-dimensional dependencies directly from data. In particular, deep learning models applied to limit order book data have demonstrated strong predictive performance [7, 8]. Despite these advances, a key limitation remains: such models are typically treated as black boxes, making it difficult to understand which features drive their predictions and whether the extracted patterns correspond to economically meaningful mechanisms. In this work, we adopt a different perspective. Rather than using nonlinear models solely as forecasting devices, we employ them as tools for structural discovery. The idea is to exploit the flexibility of neural networks in a first stage, allowing the model to learn nonlinear dependencies directly from the data, and then to use explainability techniques to identify and interpret the functional forms that drive the predictions. By combining neural networks with Shapley-based explainability methods [9], we decompose the fitted predictions into feature-level contributions and reconstruct the nonlinear dependence of returns and signed volumes on their lagged values. The central hypothesis of this work is that, although nonlinear models are complex at the parametric level, the relationships they learn need not be arbitrary or uninterpretable. On the contrary, if the model is capturing stable features of the data-generating process, its predictions should be organized around structured and economically meaningful mechanisms. In this sense, explainability is used not merely as a diagnostic tool for understanding a black-box predictor, but as a way to extract candidate structural relations from the model. Our analysis shows that the neural network learns a set of robust nonlinear patterns that are consistent with well-known mechanisms of market microstructure. Lagged signed volume generates sign-preserving and saturating responses, which can be interpreted as nonlinear price impact in the return equation and order-flow persistence in the signed-volume equation. Lagged returns, instead, mainly act as state variables: they modulate the strength of the order-flow signal and generate short-term reversal or attenuation effects, depending on the magnitude of the previous price move. These relationships are stable across assets with some differences between large-tick and small-tick stocks, indicating that the model captures features that are sensitive to the underlying microstructure regime. We further validate these findings by comparing the model-implied predictive surfaces with analogous conditional averages computed directly from the data. The close qualitative agreement between the two provides an important consistency check: the nonlinear structures revealed by the neural network are not artifacts of the architecture or of the explainability method, but reflect empirical regularities already present in the data. Furthermore, the advantage of the explainability analysis is that it allows these regularities to be decomposed into the separate contributions of lagged returns, lagged signed volumes, and the remaining past information. This provides access to the internal structure of the nonlinear predictor and makes it possible to isolate the functional role of individual regressors. In particular, the Shapley decomposition shows that the lagged-volume contribution has a clean sign-like and saturating structure, with no substantial spurious dependence on lagged returns. The return contribution is instead more state-dependent and is strongly organized by the sign of previous order flow: when the previous trade does not move the price, the model tends to preserve the sign of lagged signed volume, whereas a non-zero previous return generates attenuation or reversal. This decomposition allows us to design simpler parametric models. In particular, we introduce a SHAP-inspired nonlinear parametric specification that incorporates the main functional forms identified by the neural network: a saturating dependence on lagged signed volume and a state-dependent correction driven by the magnitude of lagged returns and by the sign of previous order flow. A key result of the paper is that this reduced specification improves upon the linear VAR benchmark and achieves predictive performance comparable to that of the neural network, while being more stable with respect to the choice of lag order. The model therefore retains a large fraction of the predictive power of the flexible nonlinear architecture, but in a parsimonious and economically interpretable form. Overall, our results demonstrate that explainability provides a powerful framework to uncover and interpret nonlinear structure in financial data. It offers a principled way to move from black-box prediction to interpretable model construction: neural networks can be used as discovery devices, while Shapley decompositions and aggregated conditional response analysis can translate the learned relationships into economically meaningful mechanisms and parsimonious parametric specifications. In this way, explainability helps bridge flexible machine-learning models and traditional econometric modeling, and provides a way to assess whether the learned nonlinear relationships reflect genuine market mechanisms rather than spurious statistical patterns. The remainder of the paper is organized as follows. Section 2 introduces the modeling framework and presents both linear and nonlinear specifications. Section 3 describes the data and preprocessing steps. Section 4 compares the forecasting performance of the models. Section 5 develops the explainability analysis and extracts the nonlinear dependencies, decomposing the contributions of the most relevant regressors. Section 6 introduces and calibrates a SHAP-inspired parametric nonlinear model. Section 7 concludes. 2 Price and trades models A price and trade model describes the joint dynamics of prices and trading activity, capturing the relationships between order flow, price formation, and market impact. These models aim to characterize how trades influence prices and how market participants adjust their behavior in response to price movements. They are essential for forecasting market dynamics, optimizing execution strategies, and simulating realistic market conditions. In this section, we first provide a general definition of price and trade models, outlining their key components. We then present two classes of models considered in this study: parametric autoregressive linear models and nonparametric nonlinear models. Price and trade models in discrete time assume a recursive dependence on lagged variables and a cross-dependence between price and trade volume. They can be formulated in the following general form: pt=ℱ1(pt−1,pt−2,…,vt−1,vt−2,…),vt=ℱ2(pt−1,pt−2,…,vt−1,vt−2,…), split&p_t=F_1(p_t-1,p_t-2,…,v_t-1,v_t-2,…),\\ &v_t=F_2(p_t-1,p_t-2,…,v_t-1,v_t-2,…), split (1) where t∈ℕt represents the trade time, ptp_t is the midpoint price just before the trade, and vtv_t is the signed volume of the trade at time t11 1 As usual, vt>0v_t>0 (vt<0v_t<0) for buyer- (seller-) initiated trades. We will also refer to vtv_t as the order flow.. In this formulation, the functions ℱ1F_1 and ℱ2F_2 define the model, describing the relationship between the present values of ptp_t and vtv_t with their respective lagged values. The first equation expresses the dependence of the midpoint price on its past values and on the current and past trade volumes, while the second equation models the order flow as a function of past price and volume values. This structure captures both the time series properties of each variable and their mutual influence. To ensure stationarity, the first equation is rewritten in terms of midpoint price increments, Δpt=pt+1−pt p_t=p_t+1-p_t, or price returns. Stationarity is a fundamental property in time series modeling, as it ensures that statistical properties such as mean and variance do not change over time, making the model more reliable for both estimation and forecasting. Non-stationary models can lead to spurious correlations and unstable predictions, limiting their applicability in real-world financial settings. In this work, we consider logarithmic returns, i.e., rt=log(pt+1/pt)r_t= (p_t+1/p_t). Thus, the models are always expressed in the form: rt=1(rt−1,rt−2,…,vt−1,vt−2,…),vt=2(rt−1,rt−2,…,vt−1,vt−2,…). split&r_t=G_1(r_t-1,r_t-2,…,v_t-1,v_t-2,…),\\ &v_t=G_2(r_t-1,r_t-2,…,v_t-1,v_t-2,…). split (2) We focus here on lagged dependencies only. In contrast, several traditional market-impact models, such as propagator models and Hasbrouck-type specifications, also allow for a contemporaneous effect of order flow on returns. In our framework, the contemporaneous contribution of vtv_t is not included explicitly when predicting rtr_t, but is instead absorbed into the residual component, which may therefore retain part of the instantaneous dependence between returns and volumes. We return to this point later in the paper. This choice also keeps all time series aligned in length and preserves a consistent structure across different datasets. As before, the functions 1G_1 and 2G_2 define the model structure and can be specified using different approaches. When building 1G_1 and 2G_2, a key distinction lies between parametric and nonparametric models. In a parametric approach, the functional forms of 1G_1 and 2G_2 are explicitly defined based on theoretical principles or empirical observations. This requires specifying a limited set of parameters that characterize the model, which are then estimated from the data. Parametric models offer interpretability and computational efficiency but may struggle to capture complex, nonlinear dependencies in financial time series. In contrast, a nonparametric approach does not impose a predefined functional form. Instead, the model is inferred directly from historical data, allowing for greater flexibility in capturing intricate patterns and nonlinear relationships. Machine learning techniques, such as neural networks, exemplify this approach by learning the mapping between inputs and outputs without assuming a specific mathematical structure. While nonparametric models can achieve higher predictive accuracy, they may lack interpretability and require larger datasets for robust estimation. The choice between parametric and nonparametric models depends on the trade-off between interpretability and predictive power, as well as the availability of domain knowledge to guide the modeling process. In this work, we first compare a linear parametric benchmark, implemented as a VAR model, with a nonlinear nonparametric benchmark, implemented as a deep neural network. Since the latter achieves higher predictive accuracy, especially for returns, we use explainability techniques to identify the nonlinear dependencies learned by the network. The final objective is to translate these dependencies into an interpretable nonlinear parametric model with competitive predictive performance. 2.1 Parametric Linear models We consider a linear version of the general model in Eq. (2), based on the VAR framework defined as: rt=∑i=1pairt−i+∑i=1pbivt−i+utr,vt=∑i=1pcirt−i+∑i=1pdivt−i+utv, splitr_t= _i=1^pa_ir_t-i+ _i=1^pb_iv_t-i+u_t^r,\\ v_t= _i=1^pc_ir_t-i+ _i=1^pd_iv_t-i+u_t^v, split (3) where rtr_t is the midpoint price return, and vtv_t is the signed volume. The terms utru_t^r and utvu_t^v are white noise processes, which are generally correlated, thereby capturing the residual contemporaneous dependence between returns and order flow. The parameter p represents the lag, i.e., the number of past observations used for prediction. Although the two equations appear decoupled, they are linked through these residual terms, which encode information about the instantaneous interaction between price changes and trading activity. As we discuss later, this contributions can be decoupled recovering the Hasbrouck model [6] and, as a consequence, the impact models such as the Transient Impact Model [1]. Linear models of this type have been widely used to describe price formation and order flow dynamics because of their tractability and interpretability. They provide a natural framework to model transient price impact and order flow persistence, and have been extensively employed in empirical market microstructure studies. At the same time, their linear structure imposes strong restrictions on the form of the dependencies between returns and order flow. This is particularly relevant in the context of market impact, where empirical evidence points to nonlinear and state-dependent effects. More generally, models calibrated only on public trade and quote data may fail to reproduce the price trajectories observed during real metaorder executions. In particular, recent evidence shows that both linear and nonlinear models calibrated on public market data can generate almost linear impact profiles with limited post-trade reversion, in contrast with the concave execution profiles and subsequent reversion typically observed for real metaorders [10]. This motivates the use of more flexible nonlinear specifications, while also highlighting the need to interpret their predictions with care. 2.2 Non parametric Non Linear models Figure 1: Schematic representation of non parametric model. Input data are processed through multiple hidden layers, each applying a nonlinear transformation, to produce the final output. A single function fθf_θ maps the input data to both return and volume forecasts, allowing the network to learn joint representations and capture complex nonlinear dependencies between the two variables. In contrast to linear parametric models, nonparametric models do not impose a predefined functional form, but instead learn the relationship between inputs and outputs directly from the data. In this work, we adopt neural networks (NNs) to construct such models. Their use is motivated by the fact that feedforward neural networks provide flexible function approximators and, under suitable conditions, can approximate a broad class of nonlinear mappings [11, 12, 13]. Neural networks define nonlinear mappings obtained by composing parametrized transformations across multiple layers. Given an input sequence tX_t of length p, the network produces a prediction through a function fθ(t)f_θ(X_t), where θ denotes the collection of trainable parameters. As in the VAR setting, the parameter p determines the maximum lookback horizon and will therefore be referred to as the lag. The function fθf_θ is not specified ex ante, but is learned from the data by minimizing a loss function measuring the discrepancy between predictions and realized values, through gradient-based optimization [13]. Thus, in this case, we have the following relation between lagged quantities and predictions xt=fθ(t)+ut,x_t=f_θ(X_t)+u_t, (4) where, in complete analogy with the linear model, xt=(rt,vt)x_t=(r_t,v_t) are the outputs of the model, ut=(utr,utv)u_t=(u_t^r,u_t^v) are the correlated residuals, while t=(rt−1,…,rt−p,vt−1,…,vt−p)X_t=(r_t-1,...,r_t-p;v_t-1,...,v_t-p) are the lagged inputs. In analogy with the VAR specification, the lag index ℓ denotes the temporal distance from the predicted variable. Thus, the ℓ -th lag of the DNN corresponds to the input pair (rt−ℓ,vt−ℓr_t- ,v_t- ), with ℓ=1,…,p =1,…,p. Unlike in the VAR case, however, these lagged variables enter as coordinates of the nonlinear map fθf_θ, so that their effect on the prediction is not restricted to be linear or additive. Figure 1 illustrates the neural network framework used in this work. In contrast to the linear setting, where the dynamics of returns and volumes are described by separate equations and contemporaneous dependence is reflected in correlated residuals, here both components are captured jointly by a single nonlinear mapping. Several neural-network architectures can be used to parameterize the function fθf_θ, including convolutional and recurrent architectures designed to exploit local or long-range temporal dependencies in sequential data [14, 15, 13]. In this work, we use a deep feed-forward fully connected neural network with hidden layers as a flexible nonlinear benchmark. For a lag order p, the input vector is tX_t so that the input dimension is 2p2p. The network has two outputs, corresponding to the one-step-ahead predictions of rtr_t and vtv_t. The architecture consists of three hidden dense layers with 6464, 6464, and 3232 neurons, respectively. Each hidden layer uses a hyperbolic tangent activation function, while the output layer is linear, as appropriate for a regression task. All models are trained by minimizing the mean squared error (MSE) loss using the Adam optimizer [16]. This choice is motivated by two considerations. First, deep feed-forward neural networks with hidden layers provide flexible nonlinear function approximators [11, 12]. Second, compared with more complex sequence-specific architectures, this architecture defines a more transparent and easily manipulable mapping from lagged inputs to predictions. This makes it particularly suitable for the aggregated dependence and Shapley-based analyses developed below [17, 9]. The same architecture is used across stocks and lag specifications, so that differences in performance and explainability patterns are not driven by changes in the neural-network design. While we focus on this architecture for interpretability and reproducibility, the approach could in principle be extended to other nonlinear architectures. 2.3 Residual Contemporaneous Dependence Both the linear VAR model in Eq. (3) and the nonlinear models introduced above use lagged information only. In particular, the return equation does not include the contemporaneous signed volume vtv_t as a direct regressor. At first sight, this may seem restrictive from a market microstructure perspective, since trades are expected to have an immediate effect on prices. The point, however, is that contemporaneous order-flow effects are not ignored. Since the models estimate reduced-form conditional expectations based on past information, any instantaneous relation between unexpected order flow and unexpected price changes remains encoded in the residuals. Let tX_t denote the lagged information set used by the model. For a generic reduced-form predictor we write rt=fr(t)+utr,vt=fv(t)+utv, splitr_t&=f_r(X_t)+u_t^r,\\ v_t&=f_v(X_t)+u_t^v, split (5) where utru_t^r and utvu_t^v are the components of returns and signed volumes that cannot be predicted from lagged variables alone. If trades have an immediate effect on prices, these two residuals should remain contemporaneously dependent even after the lagged dynamics have been removed. This distinction is important because directly including vtv_t in the return equation may lead to a misleading interpretation. Indeed, vtv_t contains both a component predictable from past information and a genuinely contemporaneous innovation: vt=[vt∣t]+utv.v_t=E[v_t _t]+u_t^v. Therefore, a model of the form rt=G(t,vt)+ηtr_t=G(X_t,v_t)+ _t receives as input not only the unexpected order-flow shock utvu_t^v, but also the predictable part of order flow, which is already determined by lagged variables. In nonlinear models, these two effects are generally mixed inside the same function G, making it difficult to distinguish lagged predictability from instantaneous impact. For this reason, we first estimate the reduced-form model using only lagged information and then analyze the residual dependence. Following a Hasbrouck-type recursive interpretation, we separate the predictable component of the dynamics from the contemporaneous innovation. We assume that, in trade time, the innovation in order flow is contemporaneously prior to the quote revision. This motivates the triangular structural representation utv=ϵtv,utr=ϕ(utv)+ϵtr.u_t^v= _t^v, u_t^r=φ(u_t^v)+ _t^r. (6) Here, ϵtv _t^v is the structural order-flow shock, while ϵtr _t^r is the return innovation that remains after removing the contemporaneous impact of order flow. The function ϕφ therefore captures the instantaneous price impact of the unexpected component of signed volume. In the linear Hasbrouck model this reduces to ϕ(u)=b0uφ(u)=b_0u, whereas in our nonlinear setting we allow ϕφ to be a nonparametric function. This decomposition separates two objects. The reduced-form forecast, r^tlagged=fr(t), r_t^\,lagged=f_r(X_t), uses only information available up to time t−1t-1 and is therefore the object used in the forecasting comparisons. By contrast, once the contemporaneous signed volume is observed, one may form the structural conditional prediction r^tstruct=fr(t)+ϕ^(vt−fv(t)). r_t^\,struct=f_r(X_t)+ φ\! (v_t-f_v(X_t) ). (7) This second quantity is not a purely lagged forecast; it is conditional on the realized order-flow innovation. In the empirical analysis below, we therefore proceed in two steps. First, we evaluate the forecasting performance of the reduced-form models using lagged information only. Second, we estimate the residual contemporaneous projection defined in Equation (6). A more detailed discussion of the identification logic and of the difference between direct contemporaneous conditioning and residual decomposition is provided in Appendix A. 3 Data In this work, we use data obtained from LOBSTER22 2 lobsterdata.com, containing observations from NASDAQ-listed stocks. The dataset consists of two files: the message file and the order book file. The message file records all transactions that occur within the limit order book (LOB), while the order book file provides the state of the LOB immediately after each transaction. In our analysis, we focus exclusively on “Execution of a visible limit order” events, excluding hidden orders. When adopting such a selection criterion, it is necessary to define a consistent price measure. We define the price as the mid-price of the LOB immediately before the execution of a visible order. Thus, we compute it based on the LOB configuration just before the execution of the visible order. The LOB then evolves due to subsequent events until another visible order is executed, at which point we again record the state of the LOB just before the event to define the price at that time. This approach follows the TIM convention (see [10] for further details). Subsequently, we aggregated orders occurring at the same timestamp with the same sign. To mitigate potential distortions due to auction effects, we also removed the first and last 30 minutes of trading data. Finally, to reduce the influence of extreme observations and potential outliers, we applied a quantile filter by removing observations outside the [0.005,0.995][0.005,0.995] quantile range. Furthermore, the training of the DNN has been obtained by normalizing data. For all the models presented, the first 50% of the sample is used for training, while the remaining 50% is used for out-of-sample testing. In the analysis carried out in this work we consider June 2024 for 5 small tick stocks and 5 large tick stocks. As reported in Table 1, large tick stocks are characterized by bid–ask spreads that remain close to one tick for most of the time, whereas small tick stocks typically display spreads of several ticks (see [18] for details). ticker mean spread (ticks) large tick BAC 1.16 INTC 1.16 CSCO 1.18 CMCSA 1.19 PFE 1.20 small tick AMZN 1.70 AAPL 1.82 GILD 2.03 TSLA 3.69 NVDA 6.27 Table 1: Classification of the stocks considered in the analysis, together with their mean spread measured in ticks. 4 Forecasting Performance of the Models In this section, we evaluate the forecasting performance of the considered models, comparing parametric linear models with nonparametric nonlinear models. Performance is measured using the coefficient of determination R2R^2, which quantifies how well the predicted values explain the variance of the actual data. Lag choice. Before comparing the different model classes, we discuss the choice of the lag order used in the analysis. A detailed study of the dependence of forecasting performance on the lag length is reported in Appendix B. For the VAR benchmark, predictive performance initially improves as the lag order increases, but the out-of-sample R2R^2 rapidly reaches a stable level and does not display systematic gains at longer lags. This suggests that most of the useful predictive information captured by the linear model is concentrated in the recent past. The neural-network model leads to the same qualitative conclusion, although its lag-dependence is less regular. After the first lags, the test performance fluctuates around a comparable range of values rather than increasing monotonically with the lag order. As discussed in Appendix B, these fluctuations are consistent with the fact that changing the lag order also changes the input representation and the nonlinear mapping learned by the network. Overall, both analyses indicate that very long lag structures are not required for the forecasting task considered here. We therefore set p=20p=20 in the following analysis. This choice is sufficiently large to include the relevant short-memory dependencies observed in the data, while avoiding unnecessarily long input vectors and excessive parameter proliferation. ticker out()R^2_out(r) VAR out()R^2_out(r) DNN out()R^2_out(v) VAR out()R^2_out(v) DNN large tick BAC 0.0494 0.0591 0.0299 0.0215 INTC 0.0814 0.1217 0.0198 0.0288 CSCO 0.0764 0.0982 0.0328 0.0273 CMCSA 0.0604 0.0859 0.0350 0.0486 PFE 0.0356 0.0367 0.0183 0.0091 small tick AMZN 0.0246 0.0341 0.0685 0.0705 AAPL 0.0643 0.0951 0.1147 0.1372 GILD 0.0499 0.0429 0.1055 0.1173 TSLA 0.0433 0.0906 0.1115 0.1475 NVDA 0.0166 0.0454 0.0872 0.1068 Table 2: Out-of-sample R2R^2 for the VAR and DNN models at lag p=20p=20, separately for returns r and signed volumes v. Stocks are grouped according to the large-tick/small-tick classification used in Table 1. Performances. The empirical performances of the VAR and DNN models are summarized in Table 2, which reports the out-of-sample R2R^2 obtained at lag order p=20p=20 for both returns and signed volumes. The comparison shows that the DNN generally improves upon the VAR benchmark, although the strength of the improvement depends on the predicted variable and on the tick-size regime. For returns, the improvement is rather systematic: the DNN achieves a higher out-of-sample R2R^2 for almost all stocks in the sample. This indicates that the nonlinear model is able to extract predictive structure that is not captured by the linear VAR specification. For signed volumes, the evidence is more heterogeneous. The DNN improves upon the VAR for all small-tick stocks, whereas for large-tick stocks the two models are closer and the VAR remains competitive in some cases. Overall, this suggests that nonlinear dependencies contribute to forecasting performance beyond the linear benchmark, but that the incremental gain is not uniform across variables and market regimes. An interesting pattern also emerges when comparing the relative predictability of returns and signed volumes across tick-size regimes. For large-tick stocks, returns are more predictable than signed volumes, both under the VAR and under the DNN. Conversely, for small-tick stocks, signed volumes are more predictable than returns. This difference is consistent with the distinct microstructural nature of the two regimes. In large-tick assets, the spread is typically constrained close to one tick, so price changes are strongly affected by discreteness and queue depletion at the best quotes. Recent return–volume states can therefore contain relevant information about the direction of the next price move. Signed volumes, instead, may remain noisier, since they also reflect liquidity provision, queue-positioning decisions and trading motives not fully summarized by the variables included in the model. In small-tick assets, price changes are less mechanically constrained by the minimum tick size and may depend on a richer liquidity state across several price levels. Returns are therefore harder to forecast from the recent history of returns and signed volumes alone. Signed order flow, by contrast, typically displays stronger persistence: buyer-initiated trades tend to be followed by further buyer-initiated trades, and similarly for seller-initiated trades. This persistence makes signed volumes more predictable than returns in the small-tick regime. Overall, the out-of-sample results indicate that the DNN provides a useful improvement over the VAR benchmark, especially for return prediction and for signed-volume prediction in small-tick stocks. At the same time, the comparison between returns and signed volumes shows that predictive performance is itself informative about market microstructure: large-tick stocks display a stronger link between recent order-flow pressure and discrete price changes, while small-tick stocks reveal a clearer persistence in the direction of trading pressure. These results motivate the analysis of the next section. Since the DNN generally outperforms the linear VAR benchmark, we ask which additional dependencies are captured by the nonlinear model and how they contribute to the improvement in forecasting performance. The comparison of predictive accuracy therefore leads naturally to an explainability problem: understanding what the neural network has learned about the joint dynamics of returns and signed volumes. In the following, we use explainability tools to characterize the functional relationships learned by the DNN. In particular, we investigate whether the network mainly reproduces and refines the linear effects already present in the VAR, or whether it captures genuinely nonlinear dependencies that are absent from the linear benchmark. 5 Explaining the Neural Network As shown in the previous section, nonlinear models tend to improve predictive performance relative to their linear counterparts. This suggests that the data may contain dependencies that are not fully captured by a purely linear specification. At the same time, the flexibility of neural networks makes their predictions less immediately interpretable, since the functional form learned from the data is not specified a priori. For this reason, we now analyze the trained neural network itself. The goal is to understand which dependencies the model has learned and to assess whether the improvement over the VAR benchmark is associated with genuinely nonlinear effects in the joint dynamics of returns and signed volumes. The key point of the present section is that, once a nonlinear neural network is trained on the data, explainability methods allow us to decompose its prediction into the contributions of the individual input variables. This provides information that is not directly accessible from empirical conditional averages alone. The data can reveal that a nonlinear dependence exists between returns and signed volumes, but the trained model, through its attribution structure, allows us to separate the role of lagged signed volume, lagged returns, and the remaining past history in the construction of the prediction. In this sense, the neural network is used as a structural discovery device. If the fitted model captures stable patterns that are also consistent with direct empirical evidence, its explanations can be interpreted as a guide to the minimal nonlinear mechanisms present in the data. As we show below, this makes it possible to move from a flexible nonparametric predictor to a parsimonious parametric model whose functional form is designed to reproduce the patterns identified by the neural network. As defined above, the nonlinear predictor is a map f:ℝ2p→ℝ2,f(t)=(r^t,v^t),f:R^2p ^2, f(X_t)= ( r_t, v_t ), which returns the model-based forecasts of next-period return and signed volume. In what follows, we focus on each component individually and indicate them with fr(t)f_r(X_t) and fv(t)f_v(X_t), for the return and volume components, respectively. Thus, in this section we first extract a set of diagnostic observations that allow us to assess whether the model behaves consistently with the financial mechanisms we expect. We begin by computing Shapley values, which enable us to identify interpretable and economically meaningful patterns in the model predictions. We then construct empirical response surfaces relating the relevant variables, both from the data and from the trained model, and compare them directly. This comparison provides an additional benchmark to verify that the model reproduces the main structures present in the data. Once these checks indicate that the model is capturing the expected financial relationships, we can use the Shapley decomposition with greater confidence to isolate and interpret the marginal contribution of each input variable to the model output, while keeping in mind that these contributions should be understood as model-based explanations rather than causal effects. 5.1 The Shapley decomposition For a given scalar output fa(t)f_a(X_t), with a∈r,va∈\r,v\ the Shapley additive explanations (SHAP) [9] express the prediction as the sum of a baseline term and the contributions of all input features. Denoting by J the full set of input regressors, we write fa(t)=[fa(t)]+∑j∈ϕja(t).f_a(X_t)=E[f_a(X_t)]+ _j _j^a(X_t). (8) In our setting, the input vector consists of lagged returns and lagged signed volumes, so that the decomposition can be written explicitly as fa(t)=[fa(t)]+∑ℓ=1pϕrt−ℓa(t)+∑ℓ=1pϕvt−ℓa(t).f_a(X_t)=E[f_a(X_t)]+ _ =1^p _r_t- ^a(X_t)+ _ =1^p _v_t- ^a(X_t). (9) Here: • [fa(t)]E[f_a(X_t)] is the baseline prediction, i.e. the average output of the model over the reference sample; • ϕrt−ℓa(t) _r_t- ^a(X_t) is the Shapley contribution of the lagged return rt−ℓr_t- on the predicted variable a; • ϕvt−ℓa(t) _v_t- ^a(X_t) is the Shapley contribution of the lagged signed volume vt−ℓv_t- on the predicted variable a. This representation makes it possible to quantify the marginal contribution of each lagged input separately. In particular, it allows us to distinguish the role of recent lags, such as rt−1r_t-1 and vt−1v_t-1, from the aggregate effect of more distant past values. For a scalar output faf_a, given the decomposition (9), the Shapley value of feature j at point xtx_t can be written as ϕja(xt)=∑S⊆∖j|S|!(||−|S|−1)!||![va(S∪j,xt)−va(S,xt)], _j^a(x_t)= _S \j\ |S|!\, (|J|-|S|-1 )!|J|! [v_a(S∪\j\;x_t)-v_a(S;x_t) ], (10) where va(S,xt)v_a(S;x_t) denotes the value of the model prediction on the variable a when only the subset of features S is observed. In practice, this quantity is defined by marginalizing or replacing the missing features with respect to a reference distribution, usually represented by a background sample. Thus, Eq. (10) expresses the contribution of feature j as the average change in the model prediction obtained when j is added to all possible subsets of the remaining features. It is useful to keep in mind what this decomposition implies in the linear benchmark defined in Eq. (3). Under the background formulation adopted here, since the VAR predictor is linear in the lagged regressors, the Shapley contribution of each input reduces to its linear coefficient times the deviation of that input from its background mean. For instance, for the return equation one obtains ϕrt−ir(t)=ai(rt−i−[rt−i]),ϕvt−ir(t)=bi(vt−i−[vt−i]), _r_t-i^r(X_t)=a_i (r_t-i-E[r_t-i] ), _v_t-i^r(X_t)=b_i (v_t-i-E[v_t-i] ), (11) while for the signed-volume equation ϕrt−iv(t)=ci(rt−i−[rt−i]),ϕvt−iv(t)=di(vt−i−[vt−i]). _r_t-i^v(X_t)=c_i (r_t-i-E[r_t-i] ), _v_t-i^v(X_t)=d_i (v_t-i-E[v_t-i] ). (12) Thus, in the VAR case, the Shapley decomposition recovers the usual linear contribution of each lagged regressor, centered with respect to the reference sample. This observation provides a useful benchmark for the nonlinear case. For the neural network, the same decomposition summarizes a nonparametric prediction in terms of feature-wise contributions, and deviations from the linear pattern above can reveal dependencies that are not captured by the VAR specification. The exact computation of Eq. (10) is infeasible in our setting, since it requires averaging over all possible subsets of the input variables. With 2p2p input features, corresponding to p lagged returns and p lagged signed volumes, the number of coalitions grows exponentially with the lag length. For this reason, we compute Shapley-based attributions using the DeepExplainer method implemented in the SHAP library [9]. This method provides an efficient approximation of Shapley values for neural networks by combining ideas from DeepLIFT-style backpropagation rules with an expectation over a background reference sample [19, 9]. In this way, the contribution of each feature is evaluated relative to the distribution of background observations, rather than by explicitly enumerating all feature coalitions. Therefore, throughout the paper we refer to these quantities as SHAP values. The resulting attributions preserve the additive interpretation in Eq. (8) and allow us to decompose the fitted nonlinear prediction into feature-level contributions. 5.2 Shapley values In the following we want to understand whether the improvement observed in the previous section is associated with economically meaningful dependencies in the data. In other words, we would like to verify that the model is not only exploiting its greater flexibility to improve forecasting accuracy, but is also learning patterns that are consistent with plausible financial and microstructural mechanisms. A first diagnostic in this direction is provided by the raw Shapley values. Since Shapley values decompose each prediction into feature-wise contributions, they allow us to inspect which lagged returns and signed volumes are used by the network, and with which sign and intensity. Systematic patterns in these contributions can therefore provide preliminary evidence that the model is extracting stable structure from the return–volume dynamics, rather than relying on purely noisy fluctuations. In Figure 2, we show the Shapley values for a representative stock, namely AAPL. The results presented below are qualitatively robust across all considered tickers, both large- and small-tick, and AAPL can therefore be viewed as a benchmark case. We consider four configurations, corresponding to the effect of lagged signed volumes and lagged returns on both predicted returns and predicted volumes. Each panel reports the distribution of Shapley values across observations for each lag ℓ=1,…,p =1,…,p, with color indicating the realized value of the corresponding feature: blue for low values and red for high values. (a) Effect of lagged signed volume vt−ℓv_t- on predicted volume vtv_t. Positive (negative) volumes generate positive (negative) contributions, reflecting persistence of order flow. (b) Effect of lagged signed volume vt−ℓv_t- on predicted return rtr_t. Positive order flow induces positive contributions, consistent with price impact. (c) Effect of lagged return rt−ℓr_t- on predicted volume vtv_t. Large price moves increase expected trading activity, consistent with the volatility–volume relationship. (d) Effect of lagged return rt−ℓr_t- on predicted return rtr_t. Positive past returns generate negative contributions, revealing short-term mean reversion. Figure 2: Shapley decomposition for AAPL. Colors indicate feature values (blue: low, red: high). A first important feature is the strong concentration of Shapley magnitude at the first lag. Across all panels, the largest contributions are associated with ℓ=1 =1, while the explanatory power rapidly decays for larger lags. This suggests that the neural network mainly relies on the most recent return–volume pair when forming its predictions. The sign structure of the Shapley values is also economically meaningful. Lagged signed volumes generate sign-preserving contributions: positive order flow pushes both the predicted volume and the predicted return upward, while negative order flow pushes them downward. This captures, respectively, order-flow persistence in the volume equation and directional price impact in the return equation. Lagged returns display a different pattern. In the return prediction, positive past returns tend to generate negative contributions and negative past returns positive contributions, revealing a short-term reversal effect. In the volume prediction, the contribution of lagged returns is more symmetric and mainly related to the magnitude of past price changes, consistently with the link between volatility and trading activity. Overall, the SHAP profiles indicate that the neural network has learned interpretable short-horizon mechanisms: order-flow persistence, directional price impact, and return reversal. A more detailed analysis of the distribution of Shapley importance across lags is reported in Appendix C. As already suggested by Figure 2, Shapley importance is strongly concentrated at the first lag across all prediction blocks and both tick-size regimes. This motivates the focus on the first-lag coordinates (rt−1,vt−1)(r_t-1,v_t-1) in the predictive-surface analysis below. 5.3 Model-Implied Predictive Surface As justified by the previous observations, we focus on the most recent return–volume observations. Importantly, we do not re-estimate the model using only the first lag. The predictor is still fitted using the full p-lag input vector tX_t; we only examine its dependence on the most recent pair (rt−1,vt−1)(r_t-1,v_t-1), while the remaining lagged coordinates retain the values observed in each empirical realization. This choice allows us to isolate and visualize the dominant short-term dependence without discarding the additional information contained in deeper lags. Re-fitting a model using only t−1t-1 would instead define a different forecasting specification and would therefore confound the analysis of the fitted nonlinear response with a change in the information set available to the model. At the pointwise level, this is simply the map t↦f(t),X_t f(X_t), viewed as a function of the pair (rt−1,vt−1)(r_t-1,v_t-1), where f denotes the fitted model predictor. In practice, when we plot the cloud of points (rt−1,vt−1,f(t)), (r_t-1,\,v_t-1,\,f(X_t) ), we are visualizing the model response on the empirical support of the data. This representation is useful because it shows directly how the fitted nonlinear predictor behaves in regions of the state space that are actually visited by the market. Moreover, this object has a natural empirical counterpart obtained by replacing the model prediction f(t)f(X_t) with the realized target variable. Comparing the model-implied response with the corresponding empirical cloud therefore allows us to assess whether the dependencies reconstructed by the neural network are genuinely present in the data, or whether they are artifacts introduced by the flexibility of the model. To obtain more readable objects, we consider bin-averaged versions of the predictive surface. The construction, however, differs depending on whether the horizontal axis is given by vt−1v_t-1 or by rt−1r_t-1. When the response is represented as a function of vt−1v_t-1, we partition the support of vt−1v_t-1 into uniform bins, denoted by Buvu\B_u^v\_u. We then compute conditional averages of the model prediction within each bin, further splitting the sample according to the regime of rt−1r_t-1, namely rt−1>0r_t-1>0, rt−1<0r_t-1<0, and rt−1=0r_t-1=0. The corresponding model-implied response is therefore m^a|v(u,kr):=[fa(t)∣vt−1∈Buv,sign(rt−1)=kr], m_a|v(u,k_r):=E\! [f_a(X_t) v_t-1∈ B_u^v,\ sign(r_t-1)=k_r ], (13) where kr∈−1,0,+1k_r∈\-1,0,+1\ indexes the sign regime of rt−1r_t-1 and a∈r,va∈\r,v\. When the response is represented instead as a function of rt−1r_t-1, uniform binning is less appropriate, since log-returns are concentrated around a finite set of local clusters induced by price discreteness. For this reason, we identify the relevant regions of the support of rt−1r_t-1 through a Gaussian kernel density estimator, and define clusters by locating the local maxima of the estimated density and assigning each observation to the interval delimited by the midpoints between adjacent peaks. Denoting these data-driven clusters by Buru\B_u^r\_u, we define m^a|r(u,kv):=[fa(t)∣rt−1∈Bur,sign(vt−1)=kv], m_a|r(u,k_v):=E\! [f_a(X_t) r_t-1∈ B_u^r,\ sign(v_t-1)=k_v ], (14) where kv∈−1,+1k_v∈\-1,+1\ denotes the sign of vt−1v_t-1. Clearly, volumes are non vanishing and cannot have zero sign. Similarly, we can define the same quantities on the Shapley functions sva(u,kr) s_v^a(u,k_r) :=[ϕvt−1a(t)∣vt−1∈Buv,sign(rt−1)=kr], :=E\! [ _v_t-1^a(X_t) v_t-1∈ B_u^v,\ sign(r_t-1)=k_r ], (15) sra(u,kv) s_r^a(u,k_v) :=[ϕrt−1a(t)∣rt−1∈Bur,sign(vt−1)=kv], :=E\! [ _r_t-1^a(X_t) r_t-1∈ B_u^r,\ sign(v_t-1)=k_v ], sr∗,a(u,kr) s_r^*,a(u,k_r) :=[ϕrt−1a(t)∣vt−1∈Bvu,sign(rt−1)=kr], :=E\! [ _r_t-1^a(X_t) v_t-1∈ B^v_u,\ sign(r_t-1)=k_r ], sv∗,a(u,kv) s_v^*,a(u,k_v) :=[ϕvt−1a(t)∣rt−1∈Bru,sign(vt−1)=kv], :=E\! [ _v_t-1^a(X_t) r_t-1∈ B^r_u,\ sign(v_t-1)=k_v ], srest,va(u,kr) s_rest,v^a(u,k_r) :=[ϕresta(t)∣vt−1∈Bvu,sign(rt−1)=kr], :=E\! [ _rest^a(X_t) v_t-1∈ B^v_u,\ sign(r_t-1)=k_r ], srest,ra(u,kv) s_rest,r^a(u,k_v) :=[ϕresta(t)∣rt−1∈Bru,sign(vt−1)=kv]. :=E\! [ _rest^a(X_t) r_t-1∈ B^r_u,\ sign(v_t-1)=k_v ]. Here, we have defined ϕresta=∑ℓ=2p[ϕrt−ℓa(t)+ϕvt−ℓa(t)]. _rest^a= _ =2^p [ _r_t- ^a(X_t)+ _v_t- ^a(X_t) ]. (16) Using these definitions, we can recover Equation (8) in conditional form as m^a|v(u,kr)=[fa(t)]+sva(u,kr)+sr∗,a(u,kr)+srest,va(u,kr),m^a|r(u,kv)=[fa(t)]+sra(u,kv)+sv∗,a(u,kv)+srest,ra(u,kv). split m_a|v(u,k_r)&=E[f_a(X_t)]+s_v^a(u,k_r)+s_r^*,a(u,k_r)+s_rest,v^a(u,k_r),\\ m_a|r(u,k_v)&=E[f_a(X_t)]+s_r^a(u,k_v)+s_v^*,a(u,k_v)+s_rest,r^a(u,k_v). split (17) The quantities defined above provide a Shapley-based decomposition of the binned predictive surfaces. The non-starred terms, sva(u,kr)andsra(u,kv),s_v^a(u,k_r) s_r^a(u,k_v), measure the average Shapley contribution of the variable displayed on the horizontal axis of the corresponding aggregated conditional response surfaces plot. More precisely, sva(u,kr)s_v^a(u,k_r) captures the contribution of vt−1v_t-1 when the model response is studied as a function of vt−1v_t-1, conditional on the regime of rt−1r_t-1. Analogously, sra(u,kv)s_r^a(u,k_v) captures the contribution of rt−1r_t-1 when the model response is studied as a function of rt−1r_t-1, conditional on the sign of vt−1v_t-1. The starred quantities, sr∗,a(u,kr)andsv∗,a(u,kv),s_r^*,a(u,k_r) s_v^*,a(u,k_v), measure the average Shapley contribution of the conditioning variable, rather than of the variable shown on the horizontal axis. Therefore, they are useful to understand how the separation between the different conditional curves is modulated by the second covariate. For example, when the response is plotted as a function of vt−1v_t-1, the quantity sr∗,a(u,kr)s_r^*,a(u,k_r) measures the contribution of rt−1r_t-1, which is the variable defining the conditioning regimes. Similarly, when the response is plotted as a function of rt−1r_t-1, the quantity sv∗,a(u,kv)s_v^*,a(u,k_v) measures the contribution of vt−1v_t-1, which determines the conditioning regimes. For the linear VAR benchmark, these quantities have a simple expected structure. Since the fitted response is affine in the lagged variables, the non-starred profiles are expected to be linear functions of the variable displayed on the horizontal axis. By contrast, the starred profiles measure the contribution of the conditioning variable and should therefore be approximately flat within each conditioning regime, up to residual effects induced by empirical correlations between rt−1r_t-1 and vt−1v_t-1 and by the binning procedure. 5.3.1 Aggregated conditional surfaces plots results In this section, we first consider the aggregated surfaces defined on the left-hand side of Equations (17). These surfaces provide a useful diagnostic tool to assess whether the model learns relationships that are meaningful from a financial point of view. Moreover, as anticipated, they can be directly benchmarked against the data whose results are reported in Appendix D. In figures 3 and 4 we report the resulting plots for the left hand side of Equations (17), namely for Eq.s (13)-(14) for the small tick (ST) and large tick (LT) buckets, starting from ticker-specific grouped aggregated conditional response surfaces computed individually for each stock. These results are obtained starting from the four grouped aggregated conditional response surfaces computed for each stock: the effect of rt−1r_t-1 on vtv_t conditional on the sign of vt−1v_t-1, the effect of rt−1r_t-1 on rtr_t conditional on the sign of vt−1v_t-1, the effect of vt−1v_t-1 on vtv_t conditional on the regime of rt−1r_t-1, and the effect of vt−1v_t-1 on rtr_t conditional on the regime of rt−1r_t-1. Then an average across stocks is performed. The histograms reported in transparency represent the number of observations in each bin, thus providing a measure of empirical support. In Appendix D, we report a more detailed discussion of the procedure used to obtain these results. Figure 3: Results for Equations (13) and (14) for small tick (ST) stocks. Figure 4: Results for Equations (13) and (14) for large tick (LT) stocks. From an interpretative viewpoint, the figures reveal that the model has learned clear nonlinear and state-dependent relationships between returns and signed volumes. These effects are more pronounced for ST stocks. The most direct structure appears in the panels where the horizontal axis is vt−1v_t-1. Both the predicted return, m^r|v m_r|v, and the predicted signed volume, m^v|v m_v|v, display a large discontinuity around the origin. For returns, the sign-like response to vt−1v_t-1 can be interpreted as a nonlinear price-impact relation: buy-initiated volume predicts positive price changes, while sell-initiated volume predicts negative price changes. The response is not purely linear, but tends to flatten for large absolute values of vt−1v_t-1, consistently with empirical evidence showing that price impact is nonlinear and often concave in traded volume [20, 21]. For signed volumes, the same structure reflects the persistence of order flow: positive lagged signed volume is followed by positive predicted signed volume, and negative lagged signed volume by negative predicted signed volume. In this sense, the model recovers the nonlinear analogue of the order-flow autocorrelation structure documented in the empirical literature [22, 23]. It is particularly interesting to compare the curves obtained by conditioning on the regime of rt−1r_t-1, namely rt−1<0r_t-1<0, rt−1=0r_t-1=0, and rt−1>0r_t-1>0. The three regimes generate visibly distinct predictive responses, indicating that the effect of lagged signed volume depends on the recent price state. Especially in the ST case, a sharp transition around vt−1=0v_t-1=0, although less pronounced than in the rt−1=0r_t-1=0 regime, remains clearly visible for both rt−1<0r_t-1<0 and rt−1>0r_t-1>0. This suggests that, in ST stocks, the model captures a sharper change between sell- and buy-pressure regimes, whereas in LT stocks the response to lagged signed volume is smoother once prices have just moved. The panels where the horizontal axis is rt−1r_t-1 convey a different message. In both m^r|r m_r|r and m^v|r m_v|r, the dominant feature is not a simple monotonic dependence on past returns, but rather a strong separation between the curves associated with vt−1>0v_t-1>0 and vt−1<0v_t-1<0. This indicates that, even when the model response is projected onto lagged returns, the sign of lagged order flow remains the main state variable organizing the prediction. More specifically, conditional on vt−1>0v_t-1>0, the model tends to predict positive values for both rtr_t and vtv_t over most of the support of rt−1r_t-1, whereas conditional on vt−1<0v_t-1<0, the predicted values remain negative. Thus, both positive and negative lagged returns mainly affect the intensity of the response, while the direction of the prediction is primarily governed by the sign of order flow. This suggests that, at this horizon, rt−1r_t-1 does not act as an autonomous directional predictor. Rather, it works as a state variable modulating the strength of the continuation of order-flow pressure. A positive lagged signed volume reflects recent buy pressure and therefore tends to predict both positive future signed volume and positive short-term price changes; analogously, a negative lagged signed volume reflects sell pressure and leads to negative predicted responses. This behavior is consistent with evidence that short-horizon price changes are strongly related to order-flow imbalance and to the state of liquidity at the best quotes [24]. In large tick assets, where price changes are constrained by the tick size, part of the information contained in trading pressure may therefore be expressed through order-flow and queue dynamics rather than through smooth price adjustment. This interpretation is also consistent with the distinction between large- and small-tick regimes discussed in Eisler, Bouchaud and Kockelkoren [18]. Thus, rather than determining the sign of the prediction, rt−1r_t-1 primarily controls its intensity. The directional effect is strongest when rt−1r_t-1 is close to zero. As |rt−1||r_t-1| increases, the predicted response moves back toward zero. This suggests that large previous price changes weaken the continuation of the order-flow signal. Economically, this can be interpreted as an attenuation or reversal mechanism. If the previous return is small, the order-flow pressure may not yet have been fully incorporated into prices, so the model reinforces the prevailing buy/sell-pressure regime. If the previous return is large in absolute value, part of the order-flow pressure may already have been absorbed by the price move, or liquidity provision may induce a short-term reversal. In that case, the model reduces the expected continuation of the same directional signal. Note that small tick stocks exhibit much stronger nonlinearities and sharper transitions, while large tick stocks display smoother and more constrained dynamics, consistent with the role of tick size in limiting price adjustments. In Appendix D, we report the results obtained by applying the same aggregation procedure directly to the data. The only difference is that, instead of averaging the model prediction, we average the realized target variable. Overall, the comparison shows that the model reproduces fairly well, at least qualitatively, the main dependencies observed in the empirical conditional surfaces. The agreement is not perfect, as expected. In particular, especially for large-tick stocks, the model-implied profiles m^a|r m_a|r display a more pronounced asymmetry with respect to the vertical axis than the corresponding empirical surfaces. This suggests that the neural network may overemphasize some state-dependent asymmetries or decay patterns in this direction. These differences are natural, since the model provides a fitted and smoothed representation of empirical conditional averages, which are themselves affected by finite-sample fluctuations and uneven empirical support across bins. Nevertheless, the comparison remains an important consistency check: the main nonlinear structures learned by the model are also visible when the same procedure is applied directly to the data. Finally, in Appendix E we also report the corresponding aggregated predictive surfaces obtained from the VAR benchmark. As expected, the comparison reveals clear differences with respect to both the empirical surfaces and the neural-network ones. In the VAR case, the aggregated dependencies are essentially linear, with only mild deviations induced by the empirical correlations between vt−1v_t-1 and rt−1r_t-1. This shows that such correlations, once processed through a linear model, are not sufficient to reproduce the nonlinear aggregate dependencies observed in the data. These dependencies must therefore be explicitly captured by the predictive model, as happens in the nonlinear neural-network specification. Figure 5: Average binned Shapley contributions for small-tick (ST) stocks for the volume’s term. Figure 6: Average binned Shapley contributions for ST stocks for the return’s term. 5.3.2 Shap values results We now turn to the Shapley-based decomposition introduced in Equation (17). The goal of this section is to move from the analysis of the aggregate predictive response of the model to the interpretation of the isolated contributions of the individual regressors. In the previous sections, we showed that the fitted model is able to reproduce meaningful dependence patterns between lagged returns and lagged signed volumes. This was assessed both by looking at the raw Shapley values and also by looking at the model-based response surfaces compared with the corresponding empirical surfaces extracted from the data. These consistency checks suggest that the nonlinear structures learned by the model are not merely artifacts of the estimation procedure, but reflect relevant patterns already present in the data. We can therefore use the Shapley decomposition with greater confidence to study how each input variable contributes to the model prediction. In particular, Equation (17) separates two types of quantities. The non-starred components, svas_v^a and sras_r^a, measure the average Shapley contribution of the variable along which the profile is constructed: svas_v^a for profiles in vt−1v_t-1, and sras_r^a for profiles in rt−1r_t-1. They therefore describe how the contribution assigned to a variable changes as that same variable varies. The starred components, sv∗,as_v^*,a and sr∗,as_r^*,a, instead project the contribution of one variable along the direction of the other. Thus, sv∗,as_v^*,a measures how the contribution of vt−1v_t-1 varies across rt−1r_t-1, while sr∗,as_r^*,a measures how the contribution of rt−1r_t-1 varies across vt−1v_t-1. These profiles are therefore useful to detect residual cross-dependencies or interaction effects between the two regressors. This distinction is useful also as a benchmark against the linear VAR case. In a purely linear model without interaction terms, the contribution of each regressor depends only on that regressor itself. Consequently, the non-starred profiles reproduce the corresponding linear dependence, up to an additive constant, while the starred profiles are flat. A non-flat starred component therefore signals that, in the nonlinear model, the contribution assigned to one variable is modulated by the value of the other variable. This provides a way to identify nonlinear cross-dependencies learned by the model and to distinguish them from purely additive linear effects. In the following, we first analyse the volume component by considering the profiles svas_v^a and sv∗,as_v^*,a. We then apply the same analysis to the return component, studying sras_r^a and sr∗,as_r^*,a. For reasons of space, in this section we report only the results for small-tick stocks. The corresponding results for large-tick stocks, which are qualitatively similar, are provided in Appendix F. Let us first consider the results in Figure 5 obtained for the volume term. We observe that, while in the aggregated case in the right panels of Figures 3 and 4 the curves conditional on rt−1>0r_t-1>0, rt−1<0r_t-1<0, and rt−1=0r_t-1=0 are visibly separated, in this case the three curves in the left panels become much closer to each other. This is confirmed also by the right panels showing the spurious dependence of the volume term as a function of rt−1r_t-1. This dependence is, to a first approximation, quite flat, meaning that there is no relevant dependence on rt−1r_t-1 in the contribution associated with vt−1v_t-1. Thus, the separation previously observed in the aggregated profiles as a function of vt−1v_t-1, which is not fully accounted for by the direct vt−1v_t-1 contribution, is likely captured by a residual contribution in the rt−1r_t-1 term. Let us now consider the profiles constructed along the rt−1r_t-1 direction. In the non-starred component srvs_r^v, we observe a change in the role played by past returns depending on their magnitude. For small price changes, the contribution of rt−1r_t-1 tends to preserve the sign of the subsequent signed volume. By contrast, for larger absolute values of rt−1r_t-1, the return contribution changes sign, generating a mean-reverting effect on future order flow. This suggests that, for small-tick stocks, zero or very small price changes are associated with persistence in trading pressure, whereas larger price changes induce a reversal contribution to vtv_t. Furthermore, although the separation between the two curves is reduced with respect to the aggregated profiles, it is still present in the Shapley profiles when conditioning on the sign of vt−1v_t-1. This suggests that the return contribution also contains a residual component associated with lagged signed volume. The starred components sr∗,as_r^*,a provide complementary information on the separation observed in the full predictive surfaces. In fact, these two panels reveal a clear and interpretable residual dependence of the return contribution ϕrv _r^v on vt−1v_t-1. When the previous signed volume does not move the price, namely when rt−1r_t-1 is close to zero, the model attributes to the return term a contribution that reinforces the persistence of order flow and, as a consequence, maintains the sign of the previous trade. When, instead, the previous order flow is accompanied by a non-zero price change, the return contribution becomes mean-reverting and tends to reverse the sign of future signed volume and return. In Appendix F, we report the corresponding results for LT stocks. Although some differences are present, the qualitative picture remains similar. This is particularly true for the starred components, which again show persistence of the order-flow sign when there is no price movement, and a sign reversal when the lagged return is non-zero. The same mechanism is also visible in the non-starred profiles: when rt−1=0r_t-1=0, the return contribution remains consistent with the sign of lagged signed volume, whereas for non-zero returns it tends to reverse it. At the same time, clear differences between ST and LT stocks emerge, especially in the non-starred components. For LT stocks, the profiles display more pronounced asymmetries and a more evident mean-reverting behavior. Moreover, their shape is closer to a linear, almost straight-line dependence than in the ST case. Overall, this Shapley-based decomposition provides a more detailed interpretation of the aggregate dependence plots. While the latter can also be constructed directly from the data, they do not separate the different contributions entering the fitted response. The Shapley decomposition, instead, allows us to decompose the aggregate predictive patterns into the individual contributions associated with the different regressors. Before constructing the reduced parametric model, we first extract the residual contemporaneous component associated with the unexpected part of signed volume. Since the models considered so far use only lagged information, the instantaneous contribution of vtv_t remains encoded in the residual dependence between return and signed-volume innovations. As shown below, this residual decomposition reveals interpretable nonlinear patterns, similar to those identified in the lagged predictive structure. 5.4 The structural shocks We now use the fitted nonlinear model to extract structural shocks defined in Equation (6) from the reduced-form residuals. In this way, we aim to extract the contemporaneous effect of order flow on returns, not included in the previous model, which remains encoded in the residual dependence between return and signed-volume innovations. Following the triangular structural interpretation in Equation (6), the resulting structural shocks are estimated as ϵ^tv=utv,ϵ^tr=utr−ϕ^(utv). ε_t^v=u_t^v, ε_t^r=u_t^r- φ(u_t^v). (18) We estimate ϕφ using a nonparametric binned estimator, ϕbin(uj)=[utr∣utv∈Bj], _bin(u_j)=E [u_t^r u_t^v∈ B_j ], (19) which provides an agnostic estimate of the conditional mean ϕ(u)=[utr∣utv=u].φ(u)=E[u_t^r u_t^v=u]. This allows us to capture possible nonlinearities in the residual contemporaneous impact. Moreover this choice implies [ϵtr∣utv]=0,E[ _t^r u_t^v]=0, and therefore, in population, Cov(ϵtr,utv)=0.Cov( _t^r,u_t^v)=0. To quantify the relevance of the residual impact map, we compute the fraction of variation in the return residual accounted for by the contemporaneous order-flow projection, Rϕ2=1−∑t(utr−ϕ^(utv))2∑t(utr−ur¯)2.R_φ^2=1- _t (u_t^r- φ(u_t^v) )^2 _t (u_t^r- u^r )^2. (20) Here, utru_t^r is the component of the return that remains unexplained by lagged returns and lagged signed volumes, while ϕ^(utv) φ(u_t^v) is the component of this residual return accounted for by the contemporaneous innovation in signed volume. Thus, Rϕ2R_φ^2 measures how much of the lag-unexplained return variation is captured by the residual contemporaneous order-flow relation. If Rlag2R_lag^2 denotes the fraction of total return variation explained by the lagged predictor alone, then for zero-mean return residuals, ur¯=0, u^r=0, the corresponding gain in total explanatory power is ΔRtotal2=(1−Rlag2)Rϕ2. R_total^2= (1-R_lag^2 )R_φ^2. (21) This zero-mean condition holds in population when the reduced-form predictor coincides with the conditional mean, and is expected to hold approximately in finite samples. Hence, Rϕ2R_φ^2 should be interpreted as the fraction of the remaining return variation explained after the lagged predictive component has been removed. The analysis is performed separately for LT and ST stocks. This is consistent with the evidence presented above, which shows that the two buckets differ in their predictive surfaces, Shapley profiles, and relative predictability of returns and signed volumes. Since ϕφ captures the instantaneous impact of unexpected order flow on unexpected returns, its shape and explanatory power may also depend on the tick-size regime. Empirical results. Table 3 reports the results of the binned structural decomposition for each ticker, together with equal-weight averages within the LT and ST buckets. Table 3: Structural shock decomposition by ticker. Bucket Ticker Rϕ2R_φ^2 LT BAC 0.163 LT CMCSA 0.117 LT CSCO 0.148 LT INTC 0.178 LT PFE 0.205 LT Average 0.162 ST AAPL 0.143 ST AMZN 0.266 ST GILD 0.222 ST NVDA 0.146 ST TSLA 0.116 ST Average 0.179 We observe that the explanatory power of the residual impact function is economically meaningful. The average value of Rϕ2R_φ^2 is 0.1620.162 for LT stocks and 0.1790.179 for ST stocks. Thus, around 16%16\%–18%18\% of the variance of the reduced-form return innovation can be explained by the contemporaneous innovation in signed volume. This fraction is sizeable, especially because the lagged nonlinear model has already removed the predictable component based on past returns and volumes. The result therefore confirms that a relevant part of price formation occurs through an instantaneous order-flow channel. At the ticker level, Rϕ2R_φ^2 ranges from 0.1170.117 to 0.2050.205 in the LT bucket and from 0.1160.116 to 0.2660.266 in the ST bucket, showing that the residual contemporaneous impact is present across all assets, although with heterogeneous intensity. The shape of the estimated residual impact functions is shown in Figure 7. In both LT and ST stocks, the function ϕ φ is increasing and sign-preserving: negative order-flow innovations generate negative return innovations, while positive order-flow innovations generate positive return innovations. The response is not purely linear. It is steeper around the origin and tends to flatten in the tails, suggesting a nonlinear and saturating instantaneous impact. This behavior is consistent with the interpretation of ϕφ as a residual price-impact function: after removing the lagged predictable component, unexpected buy pressure pushes the return innovation upward, while unexpected sell pressure pushes it downward. The two buckets display similar qualitative shapes, but with different scales. The LT residual impact curve is defined over a wider range of order-flow innovations and reaches a larger return-residual magnitude. The ST curve is more concentrated around the origin and displays a smoother saturation pattern. This suggests that the contemporaneous residual impact channel is present in both regimes, but its scale and effective support are affected by tick-size constraints. Figure 7: Average binned residual impact function ϕ^(utv) φ(u_t^v) for LT stocks (left) and ST stocks (right). The function maps order-flow innovations into return innovations after removing the lagged predictable component. 6 From Explanations to a Parametric Model The analysis developed in the previous sections suggests that the nonlinear neural network learns economically interpretable dependencies. In particular, the aggregated conditional response surfaces describe the model-implied response in selected regions of the state space and are qualitatively validated by the corresponding empirical surfaces reconstructed directly from the data. The Shapley decomposition provides an additional layer of information by separating this aggregate response into the contributions of the individual input variables. We use these insights to introduce a simpler, calibratable, and interpretable parametric model. The idea is therefore to construct a reduced-form model based on the explanations developed in the previous section. In this way, the model keeps the interpretability and tractability of a VAR-type specification, but incorporates the nonlinear mechanisms revealed by the neural network. 6.1 One-lag model construction The SHAP analysis showed that a large part of the nonlinear structure learned by the neural network is concentrated in the most recent return–volume coordinates, (rt−1,vt−1)(r_t-1,v_t-1). This does not imply that higher-order lags are negligible: they still contain useful predictive information, in particular through the persistence of order flow. Rather, the evidence suggests that the strongest nonlinear effects are already visible at the first lag. For this reason, we first ask how much explanatory power is lost when replacing a linear model with several lags by a nonlinear Markovian specification based only on the most recent return and signed-volume observations. In this way, the analysis compares two different sources of predictability: the use of a richer lag structure in a linear model and the use of nonlinear interactions in a low-dimensional, one-lag representation. Volume terms. We first consider the volume-related contributions. For small-tick stocks, these profiles are reported in Figure 5, while the corresponding large-tick profiles are reported in Appendix Figure 18. The figure shows the profiles svrs_v^r and svvs_v^v, which describe the average Shapley contribution of the lagged signed-volume term as a function of the lagged signed volume itself, respectively for the return and signed-volume predictions. We also report the corresponding starred profiles, sv∗,rs_v^*,r and sv∗,vs_v^*,v. These quantities measure how the remaining Shapley contribution varies when projected on the same signed-volume coordinate, and are therefore useful to detect possible spurious interaction effects attributed to the other variables. As widely discussed, this term appears to have no spurious dependencies on rt−1r_t-1 but only a sign-like dependence on the signed volume. Given this, we can guess the term Φva _v^a, where a∈r,va∈\r,v\, to be an odd, sign-preserving, saturating function: Φva(v)=Avatanh(βvav). _v^a(v)=A_v^a ( _v^av). (22) This preserves the sign-like directional response observed in the projections of Shapley term ϕva _v^a, while avoiding a discontinuous jump at zero. The steepness parameter allows the same specification to capture smoother LT responses and sharper ST responses. This term appears in the same way in the return equation, a=ra=r, and in the volume equation a=va=v. Return terms. We now consider the return-related contributions, reported in Figure 6 for small-tick stocks. The analogous profiles for large-tick stocks are reported in Appendix Figure 18. The interpretation is more delicate than for the volume terms, because a substantial part of the structure associated with the return contribution is organized by lagged signed volume. A relevant structure is present in the starred components. These profiles show how the remaining Shapley contribution varies when projected on the signed-volume coordinate. They indicate that the sign of vt−1v_t-1 plays a central role in determining the direction of the next-step contribution. In particular, when the previous trade did not induce a price change, rt−1=0r_t-1=0, the contribution tends to preserve the sign of vt−1v_t-1: positive lagged signed volume produces a positive contribution, while negative lagged signed volume produces a negative contribution, both for the next return and for the next signed volume. Conversely, when a price change occurred, rt−1≠0r_t-1≠ 0, the contribution tends to generate an opposite pressure, favouring a sign reversal in the next-step prediction. This means that the return-related term contains information about the joint return–volume state, rather than a purely return-driven effect. The same mechanism is visible directly in the non-starred return profiles. Around rt−1=0r_t-1=0, the sign of the contribution is mainly determined by the sign of vt−1v_t-1: the two conditional curves lie on opposite sides of zero, consistently with the sign-preserving behavior observed in the starred components. Thus, even when the profile is plotted as a function of rt−1r_t-1, the dominant separation between curves is induced by the lagged signed-volume state. For non-zero lagged returns, the dependence on rt−1r_t-1 becomes more subtle and differs between tick-size regimes. In the small-tick case, shown in Figure 6, the dependence is approximately symmetric with respect to the sign of rt−1r_t-1. As |rt−1||r_t-1| increases, the sign-preserving contribution associated with vt−1v_t-1 is progressively weakened and may become mean-reverting. This motivates the following parsimonious reduced-form specification: Φra(r,v)=q[Ara−Bratanh(λra|r|)],q=sign(vt−1). _r^a(r,v)=q [A_r^a-B_r^a ( _r^a|r| ) ], q=sign(v_t-1). (23) This functional form captures the dominant symmetric component of the empirical profiles: for small lagged returns, the contribution is mainly aligned with the sign of the previous signed volume, while for larger |rt−1||r_t-1| this contribution is attenuated through the saturating hyperbolic tangent term. The large-tick profiles, reported in Appendix Figure 18, suggest that this symmetry is only approximate. In particular, away from rt−1=0r_t-1=0, the dependence on positive and negative lagged returns appears partially asymmetric, indicating a possible symmetry-breaking correction. Such a correction could be introduced, for instance, through an additional term depending on the signed return rt−1r_t-1, or more generally on the relative sign variable qrt−1qr_t-1. We do not include this correction in the baseline specification, since our goal is to construct a parsimonious reduced model that captures the dominant and most stable empirical pattern across stocks. So we keep Equation (23) as a first-order reduced form for the symmetric component of the return contribution, while the asymmetric correction is left as a natural extension to be tested separately. Note that taking Bra>AraB_r^a>A_r^a, the function Φra _r^a changes sign at a finite threshold in |r||r|. Indeed, the term in brackets, Ara−Bratanh(λra|r|),A_r^a-B_r^a ( _r^a|r|), is positive for sufficiently small |r||r|, since tanh(λra|r|)≃0 ( _r^a|r|) 0 when r is close to zero. In this region, Φra(r,v)≃qAra, _r^a(r,v) qA_r^a, so the contribution has the same sign as the lagged signed volume. Therefore, when vt−1>0v_t-1>0, the term contributes positively to the prediction of rtr_t, while when vt−1<0v_t-1<0, it contributes negatively. In this regime, the previous return is too small to indicate that the order-flow pressure has already been incorporated into prices, and the model therefore reinforces the directional signal carried by vt−1v_t-1. As |r||r| increases, the quantity tanh(λra|r|) ( _r^a|r|) increases from zero to one. If Bra>AraB_r^a>A_r^a, there exists a threshold |r|⋆=1λraarctanh(AraBra)|r| = 1 _r^aarctanh\! ( A_r^aB_r^a ) such that Ara−Bratanh(λra|r|)>0for |r|<|r|⋆,A_r^a-B_r^a ( _r^a|r|)>0 |r|<|r| , while Ara−Bratanh(λra|r|)<0for |r|>|r|⋆.A_r^a-B_r^a ( _r^a|r|)<0 |r|>|r| . Beyond this threshold, the sign of Φra _r^a becomes opposite to the sign of q=sign(vt−1)q=sign(v_t-1). Hence, for large absolute returns, the term no longer reinforces the order-flow regime, but instead attenuates or reverses it. This behavior is consistent with the combined evidence from the aggregated surfaces and the SHAP decomposition: the return-state contribution reinforces the order-flow regime when the previous price move is small, while it attenuates or reverses it when the previous return is large in magnitude. Summary of the parametric model. Given these specifications we can, thus, summarize the two equations as rt=αr+Φvr(v)+Φr(r,v)=αr+Avrtanh(βvrvt−1)+sign(vt−1)[Arr−Brrtanh(λr|rt−1|)]+utr,vt=αv+Φv(v)+Φrv(r,v)=αv+Avvtanh(βvvt−1)+sign(vt−1)[Arv−Brvtanh(λrv|rt−1|)]+utv. splitr_t&= _r+ _v^r(v)+ _r^r(r,v)= _r+A_v^r ( _v^rv_t-1)+sign(v_t-1) [A_r^r-B_r^r ( _r^r|r_t-1|) ]+u_t^r,\\ v_t&= _v+ _v^v(v)+ _r^v(r,v)= _v+A_v^v ( _v^vv_t-1)+sign(v_t-1) [A_r^v-B_r^v ( _r^v|r_t-1|) ]+u_t^v. split (24) The first nonlinear term in each equation captures the direct contribution of lagged signed volume: directional order-flow impact in the return equation and order-flow persistence in the volume equation. The second term captures the correction induced by the previous return state. For small |rt−1||r_t-1|, this term reinforces the sign of vt−1v_t-1; for sufficiently large |rt−1||r_t-1|, it attenuates or reverses the order-flow signal. Thus, the model provides a compact parametric representation of the main mechanisms isolated by the aggregated surfaces and SHAP analyses. Recovering the SHAP profiles. The projected quantities in Equation (15) and the relations in Equation (17) can be recovered by conditioning Equations (24) in the directions rt−1,vt−1r_t-1,v_t-1. If we compute the response of the reduced model equation as a function of vt−1v_t-1, conditioning on the sign kr=sign(rt−1)k_r=sign(r_t-1), we obtain approximately ma|vpar(u,k)≃αa+Φva(vu)+[Φra(r,v)∣vt−1∈Buv,sign(rt−1)=kr]=αa+sv,para+sr,par∗,a.m_a|v^par(u,k) _a+ _v^a(v_u)+E\! [ _r^a(r,v) v_t-1∈ B_u^v,\ sign(r_t-1)=k_r ]= _a+s_v,par^a+s_r,par^*,a. (25) In fact, in this decomposition, Φva(vu) _v^a(v_u) reproduces the observed non-starred profile svas_v^a, while the conditional expectation of Φra _r^a reproduces the regime-dependent shift associated with the observed sr∗,as_r^*,a. Here, we observe that sr,par∗,as_r,par^*,a reproduces the sr∗,as_r^*,a shown in the right panels of Figures 6 and 19. In fact, for |r||r| sufficiently large, we have that sr,par∗,a≃−Csign(vt−1),C=Bratanh(λra|r|)−Ara>0,s_r,par^*,a -Csign(v_t-1), C=B_r^a ( _r^a|r|)-A_r^a>0, which reproduces the blue and red curves in the top left panels in Figures 18 and 19. Differently, for |r|∼0|r| 0 we have that sr,par∗,a≃Arasign(vt−1),Ara>0,s_r,par^*,a A_r^asign(v_t-1), A_r^a>0, which reproduces the green curve. Hence, the separation between the three curves in m^a|v m_a|v is not attributed to different direct volume effects, but to the return-state correction induced by Φra _r^a, exactly as observed in the explainability of the model. Conversely, if we compute the response as a function of rt−1r_t-1, conditioning on q=sign(vt−1)q=sign(v_t-1), we obtain ma|rpar(u,q)≃αa+[Φva(vt−1)∣rt−1∈Bur,sign(vt−1)=q]+Φra(ru)=αa+sv,par∗,a+sr,para.m_a|r^par(u,q) _a+E\! [ _v^a(v_t-1) r_t-1∈ B_u^r,\ sign(v_t-1)=q ]+ _r^a(r_u)= _a+s_v,par^*,a+s_r,par^a. (26) Also in this case, sv,par∗,as_v,par^*,a reproduces what observed in the right panels of Figures 5 and 18, being just a constant as a function of rt−1r_t-1 but whose sign depends on the sign of vt−1v_t-1. Hence, again, the expectation of Φva _v^a generates the vertical separation between the vt−1>0v_t-1>0 and vt−1<0v_t-1<0 curves, corresponding to the starred contribution sv∗,as_v^*,a, while Φra(ru) _r^a(r_u) describes the shape of the response along the rt−1r_t-1 axis, corresponding to the non-starred profile sras_r^a. In the next section, we calibrate the reduced model on data and test whether the proposed specifications for the return and volume terms are sufficient to approximately recover the empirical map m m. To this end, we compare the results of the reduced parametric model with the neural-network profiles discussed in Section 5.3.1, and more specifically with the corresponding empirical profiles computed directly from the data and reported in Appendix D. 6.2 Calibration of the SHAP-Inspired Parametric Model We estimate the parameters separately for the two equations by nonlinear least squares, minimizing the mean squared prediction error on the training sample. To ensure numerical stability and preserve the qualitative structure of the model, the positivity constraints are imposed by reparametrization. In particular, the scale and steepness parameters are constrained to be positive, while the condition Brr>Arr,Brv>ArvB_r^r>A_r^r, B_r^v>A_r^v can be imposed through Brr=Arr+softplus(δr),Brv=Arv+softplus(δrv),B_r^r=A_r^r+softplus( _r^r), B_r^v=A_r^v+softplus( _r^v), where softplus(x)=log(1+ex)softplus(x)= (1+e^x) is a smooth, strictly positive function, ensuring Brr>ArrB_r^r>A_r^r and Brv>ArvB_r^v>A_r^v. This condition guarantees that the correction terms can change sign at finite thresholds in |rt−1||r_t-1|. For instance, in the return equation the threshold is |r|⋆=1λrarctanh(ArrBrr),|r| = 1 _r^rarctanh\! ( A_r^rB_r^r ), provided Brr>ArrB_r^r>A_r^r. Hence, below this threshold the correction reinforces the signed-volume regime, whereas above it the correction acts in the opposite direction. The model is calibrated on the same preprocessed data used for the VAR and neural-network analysis. (a) r^t r_t as a function of rt−1r_t-1. (b) r^t r_t as a function of vt−1v_t-1. (c) v^t v_t as a function of rt−1r_t-1. (d) v^t v_t as a function of vt−1v_t-1. Figure 8: Fitted one-lag parametric response functions for large-tick stocks. The top row reports the response of r^t r_t to rt−1r_t-1 and vt−1v_t-1, while the bottom row reports the response of v^t v_t to rt−1r_t-1 and vt−1v_t-1. (a) r^t r_t as a function of rt−1r_t-1. (b) r^t r_t as a function of vt−1v_t-1. (c) v^t v_t as a function of rt−1r_t-1. (d) v^t v_t as a function of vt−1v_t-1. Figure 9: Fitted one-lag parametric response functions for small-tick stocks. The top row reports the response of r^t r_t to rt−1r_t-1 and vt−1v_t-1, while the bottom row reports the response of v^t v_t to rt−1r_t-1 and vt−1v_t-1. The fitted response functions are reported in Figures 8–9. Overall, the one-lag parametric model is able to reproduce, at least qualitatively, the main dependencies identified through the neural-network analysis. More importantly, the same qualitative structures are consistent with the aggregate dependencies reconstructed directly from the data and reported in Appendix D, Figures 14–15. The agreement is not exact, especially when the parametric surfaces are compared with the neural-network predictive surfaces. Some differences are visible also with respect to the empirical averages, and are mainly related to asymmetries in the response functions. These discrepancies may have different origins. They may reflect residual learning errors or smoothing effects in the neural-network prediction, finite-sample fluctuations in the empirical conditional averages, or higher-order corrections that break the symmetric structure imposed in the baseline reduced model. In this first specification, we do not include such asymmetric corrections, in order to preserve a parsimonious and interpretable form; they represent a natural direction for future extensions. The important point is that, starting from the isolated Shapley contributions of the individual regressors, the reduced model is able to recover aggregate patterns that are comparable with those observed in the data. Thus, despite the expected differences due to its reduced complexity, the parametric specification captures the dominant qualitative structure of the return–volume dependence in both large-tick and small-tick stocks. Table 4 reports the out-of-sample R2R^2 obtained by the one-lag nonlinear parametric model for returns and signed volumes, together with the corresponding results of the VAR and DNN models estimated with p=20p=20 lags. The comparison is particularly informative because the reduced model uses only the most recent return–order-flow pair, (rt−1,vt−1)(r_t-1,v_t-1), whereas the VAR and the DNN exploit the full lagged state up to twenty lags. The main finding is that the Shapley-motivated one-lag nonlinear specification already captures most of the predictive structure in the data. For returns, the reduced model systematically improves on the linear VAR benchmark across all stocks. The gains are especially pronounced for large-tick stocks: for example, the out-of-sample R2R^2 increases from 0.04940.0494 to 0.09340.0934 for BAC, from 0.07640.0764 to 0.11560.1156 for CSCO, and from 0.03560.0356 to 0.08450.0845 for PFE. These improvements indicate that the relevant predictability is not only a lagged linear effect, but is largely associated with the nonlinear interaction between the most recent return and signed volume. The comparison with the DNN is even more revealing. Despite its much lower dimensionality, the reduced model often matches or outperforms the neural network. This is the case for most large-tick stocks in the return equation and for several small-tick stocks as well, such as AMZN, GILD, and NVDA. The DNN retains an advantage in some small-tick cases, most notably AAPL and TSLA for returns, suggesting that additional nonlinear or longer-memory effects may still be relevant for these stocks. Nevertheless, the fact that a transparent two-variable parametric specification is competitive with, and frequently superior to, a twenty-lag neural network is one of the central empirical results of the paper. A similar pattern appears for signed volumes. The reduced model improves on the VAR benchmark for every stock and is particularly effective for large-tick stocks, where it also dominates the DNN in all cases. For small-tick stocks the DNN remains stronger in some cases, notably GILD and NVDA, while the two models perform identically for TSLA. Overall, Table 4 shows that most of the out-of-sample predictability can be recovered from a parsimonious nonlinear description of the latest return–order-flow state. This provides direct empirical support for the reduced-form structure identified by the Shapley analysis and motivates the extension, considered next, in which the same nonlinear building blocks are enriched with multiple lagged inputs. Ticker Rout2(r)R^2_out(r) Rout2(v)R^2_out(v) VAR DNN Reduced VAR DNN Reduced large tick BAC 0.0494 0.0591 0.0934 0.0299 0.0215 0.0623 INTC 0.0814 0.1217 0.1065 0.0198 0.0288 0.0644 CSCO 0.0764 0.0982 0.1156 0.0328 0.0273 0.0735 CMCSA 0.0604 0.0859 0.1020 0.0350 0.0486 0.0820 PFE 0.0356 0.0367 0.0845 0.0183 0.0091 0.0478 small tick AMZN 0.0246 0.0341 0.0528 0.0685 0.0705 0.0892 AAPL 0.0643 0.0951 0.0858 0.1147 0.1372 0.1421 GILD 0.0499 0.0429 0.0723 0.1055 0.1173 0.1165 TSLA 0.0433 0.0906 0.0826 0.1115 0.1475 0.1475 NVDA 0.0166 0.0454 0.0498 0.0872 0.1068 0.0967 Table 4: Out-of-sample R2R^2 for the VAR, DNN, and reduced parametric models, separately for returns r and signed volumes (v). The VAR and DNN models are estimated at lag p=20p=20, while the reduced parametric model is the one-lag specification obtained from the Shapley-based reduced-form analysis. Stocks are grouped according to the large-tick/small-tick classification used in Table 1. 6.3 Multi-lag extension Although the one-lag specification is the most direct translation of the SHAP analysis, it is natural to ask whether the same reduced-form structure can be extended to longer histories. A fully flexible multi-lag version would assign independent nonlinear parameters to each lag. For example, the return equation would contain terms of the form Arv,ℓtanh(βrv,ℓvt−ℓ)+sign(vt−ℓ)[Arr,ℓ−Brr,ℓtanh(λrr,ℓ|rt−ℓ|)],ℓ=1,…,p.A_rv, ( _rv, v_t- )+sign(v_t- ) [A_r, -B_r, ( _r, |r_t- |) ], =1,…,p. This specification is the most general extension of the one-lag model, but it also introduces a large number of nonlinear parameters. As a consequence, it risks fitting lag-specific noise rather than persistent predictive structure, especially in light of the SHAP lag-importance analysis, which shows that most of the predictive contribution is concentrated at the first lags and then decays smoothly. We therefore impose a structured multi-lag restriction motivated by the empirical patterns documented above. The key assumption is that the qualitative nonlinear shape identified at the first lag remains the same across lags, while the strength of the corresponding contribution decreases with lag distance. Thus, instead of estimating an independent nonlinear function for each lag, we share the nonlinear shape parameters across lags and let lag-specific intensities follow a geometric decay. This preserves the economic structure suggested by the SHAP analysis, while allowing the model to incorporate longer-memory effects without introducing a separate set of nonlinear parameters for every lag. The return equation becomes rt=αr+Arv∑ℓ=1pρrvℓ−1tanh(βrvvt−ℓ)+∑ℓ=1pρrrℓ−1sign(vt−ℓ)[Arr−Brrtanh(λrr|rt−ℓ|)]+εtr, splitr_t= _r&+A_rv _ =1^p _rv -1 ( _rvv_t- )\\ &+ _ =1^p _r -1sign(v_t- ) [A_r-B_r ( _r|r_t- |) ]+ _t^r, split (27) with 0<ρrv<10< _rv<1 and 0<ρrr<10< _r<1. Similarly, the signed-volume equation is vt=αv+Avv∑ℓ=1pρvvℓ−1tanh(βvvvt−ℓ)+∑ℓ=1pρvrℓ−1sign(vt−ℓ)[Avr−Bvrtanh(λvr|rt−ℓ|)]+εtv, splitv_t= _v&+A_v _ =1^p _v -1 ( _vv_t- )\\ &+ _ =1^p _vr -1sign(v_t- ) [A_vr-B_vr ( _vr|r_t- |) ]+ _t^v, split (28) with 0<ρvv<10< _v<1 and 0<ρvr<10< _vr<1. This specification separates the nonlinear shape of the dependence from its temporal decay. The parameters A, β, B, and λ determine the common nonlinear response, while the ρ’s determine how quickly the corresponding contribution fades with the lag. In this way, the model remains parsimonious even for large p: the number of nonlinear parameters does not grow linearly with the lag length. At the same time, it incorporates longer memory in a way that is consistent with the SHAP evidence of decreasing lag importance. Table 5 reports the out-of-sample R2R^2 obtained by the VAR, the deep neural network, the one-lag reduced parametric model, and the multi-lag parametric specification with exponential decay and p=20p=20. The table shows that the reduced parametric models provide a clear improvement over the linear VAR benchmark. This holds for both returns and signed volumes and across both large-tick and small-tick stocks. For returns, the one-lag reduced model already captures most of the predictive gain. The multi-lag extension produces only moderate additional improvements: the R2R^2 increases slightly for most stocks, while remaining essentially unchanged for some small-tick names such as TSLA and NVDA. This pattern is consistent with the Shapley analysis, which shows that the dominant nonlinear dependence is concentrated at the first lag. Additional lags therefore contain some residual predictive information, but their contribution to return forecasting is limited. Importantly, the reduced parametric specifications are often comparable to the deep neural network and outperform it for several stocks, especially among large-tick names such as BAC, CSCO, CMCSA, and PFE. The gains from the multi-lag extension are more pronounced for signed volumes. While the one-lag reduced model already improves substantially over the VAR benchmark, the multi-lag specification further increases the out-of-sample R2R^2, particularly for small-tick stocks. This suggests that signed-volume dynamics contain a more persistent component than returns, so that longer histories provide useful additional information beyond the first lag. Consistently, for the signed-volume equation, the multi-lag parametric model achieves higher out-of-sample R2R^2 than the deep neural-network benchmark for all stocks in the table. Overall, these results show that the explainability-guided parametric specification retains the main nonlinear mechanisms learned by the neural network while offering a more structured and interpretable representation. The one-lag version already captures most of the nonlinear predictability of returns, whereas the multi-lag extension is especially useful for signed-volume forecasting, where persistence across lags plays a more important role. Ticker Rout2(r)R^2_out(r) Rout2(v)R^2_out(v) VAR DNN Reduced Multi-lag VAR DNN Reduced Multi-lag large tick BAC 0.0494 0.0591 0.0934 0.0945 0.0299 0.0215 0.0623 0.0623 INTC 0.0814 0.1217 0.1065 0.1073 0.0198 0.0288 0.0644 0.0651 CSCO 0.0764 0.0982 0.1156 0.1168 0.0328 0.0273 0.0735 0.0745 CMCSA 0.0604 0.0859 0.1020 0.1031 0.0350 0.0486 0.0820 0.0828 PFE 0.0356 0.0367 0.0845 0.0855 0.0183 0.0091 0.0478 0.0479 small tick AMZN 0.0246 0.0341 0.0528 0.0529 0.0685 0.0705 0.0892 0.1047 AAPL 0.0643 0.0951 0.0858 0.0866 0.1147 0.1372 0.1421 0.1552 GILD 0.0499 0.0429 0.0723 0.0731 0.1055 0.1173 0.1165 0.1248 TSLA 0.0433 0.0906 0.0826 0.0826 0.1115 0.1475 0.1475 0.1589 NVDA 0.0166 0.0454 0.0498 0.0498 0.0872 0.1068 0.0967 0.1141 Table 5: Out-of-sample R2R^2 for the VAR, DNN, one-lag reduced parametric model, and multi-lag parametric model, separately for returns r and signed volumes v. The VAR and DNN models are estimated with p=20p=20 lags. The reduced parametric model is the one-lag specification obtained from the Shapley-based reduced-form analysis, while the multi-lag parametric model uses p=20p=20 lags with a shared nonlinear shape and exponential lag decay. Stocks are grouped according to the large-tick/small-tick classification used in Table 1. Additional evidence on the role of the lag length is reported in Appendix G. There we study the evolution of the in-sample and out-of-sample R2R^2 of the multi-lag parametric model as a function of the maximum lag ℓ . The results show that predictive performance increases only mildly after the first few lags and then rapidly saturates, confirming that the dominant nonlinear dependence is short-lived. This supports the use of the exponential-decay parametrization, which incorporates residual longer-lag information without introducing a separate set of nonlinear parameters for each lag. 7 Conclusions In this paper, we studied the joint high-frequency dynamics of returns and signed volumes through the combination of nonlinear forecasting models, explainability techniques, and interpretable parametric modeling. The starting point of the analysis was the comparison between a linear VAR benchmark and a deep feed-forward neural network trained on five large-tick and five small-tick stocks. The neural network generally improves upon the VAR benchmark, especially for return prediction, suggesting that relevant nonlinear dependencies are present in the data and are not fully captured by a linear specification. Before interpreting the nonlinear model, we analyzed the dependence of forecasting performance on the lag order. For the VAR benchmark, the out-of-sample performance increases at short lags and then rapidly saturates. The neural network exhibits a less regular dependence on the lag length, with fluctuations in out-of-sample performance, but the overall evidence again points to a short-memory structure. We therefore use p=20p=20 as a conservative lag choice: it is sufficiently large to include the relevant short-horizon dependencies, while avoiding an excessive expansion of the input space. We then used SHAP-based explainability to understand whether the nonlinear model learns economically meaningful relationships between the regressors. The Shapley lag-importance analysis shows that, across both tick-size regimes, the dominant contributions are concentrated at the most recent lags. This motivates the construction of projected predictive surfaces conditional on the first return–volume pair, (rt−1,vt−1)(r_t-1,v_t-1). These surfaces allow us to study the aggregate dependence of the neural-network predictions on the most relevant lagged variables. A first important result is that the model-implied predictive surfaces are qualitatively consistent with the corresponding empirical conditional averages computed directly from the data. This provides a crucial validation step: the nonlinear structures learned by the neural network are not merely artifacts of the architecture or of the training procedure, but reflect regularities that are already present in the data. Once this consistency is established, the Shapley decomposition can be used to go beyond aggregate conditional averages and separate the predictive response into the isolated contributions of individual regressors. The isolated Shapley components reveal a clear and interpretable structure. The contribution of lagged signed volume is strongly nonlinear, sign-preserving and saturating. Positive lagged signed volume contributes positively to future returns and signed volumes, while negative lagged signed volume contributes negatively. This term does not display a substantial spurious dependence on the lagged return, indicating that the direct volume contribution is well described as a sign-like saturating function of vt−1v_t-1. The return-related contribution is more complex and state-dependent. It depends on the magnitude of the lagged return, but it is also strongly organized by the sign of the previous signed volume. In particular, when the previous trade does not move the price, rt−1=0r_t-1=0, the model tends to predict continuation: both the next return and the next signed volume are shifted in the same direction as vt−1v_t-1. Conversely, when the previous trade is associated with a non-zero price change, the model tends to generate an opposite pressure: the next return and signed volume are shifted against the sign of vt−1v_t-1. This mechanism can be interpreted as a transition from order-flow persistence when price adjustment has not yet occurred to attenuation or reversal when the previous order-flow pressure has already produced a price move. The same qualitative structure is observed for both small-tick and large-tick stocks, although with expected differences in intensity, smoothness, and asymmetry. We also analyzed the contemporaneous residual dependence between return and signed-volume innovations. Since the reduced-form models use lagged information only, the instantaneous impact of unexpected order flow remains encoded in the residuals. By estimating a nonparametric structural projection utr=ϕ(utv)+ϵtru_t^r=φ(u_t^v)+ _t^r, we find that the residual impact function is again increasing, sign-preserving, and saturating, for both large-tick and small-tick stocks. This confirms that, after removing lagged nonlinear predictability, a sizeable instantaneous order-flow channel remains in the return residuals. Building on these observations, we introduced a SHAP-inspired nonlinear parametric model. The model contains the main dependencies extracted from the explainability analysis: a saturating sign-like term in lagged signed volume and a state-dependent correction driven by the magnitude of lagged returns and by the sign of previous order flow. The one-lag version of the model provides a compact Markovian representation of the dominant nonlinear mechanisms identified by the neural network. Despite its simplicity, it reconstructs the main qualitative features of the model-implied and empirical predictive surfaces. A central empirical result of the paper is that this reduced parametric model outperforms the linear VAR benchmark and achieves predictive performance comparable to that of the deep neural network, and in several cases even higher. This shows that much of the forecasting advantage of the neural network does not come from an opaque high-dimensional structure, but from a small set of stable and interpretable nonlinear mechanisms. In this sense, the neural network acts as a discovery device, while the explainability analysis translates its learned dependencies into a parsimonious econometric specification. Finally, we extended the reduced model to multiple lags by assuming that the same nonlinear shapes apply across lagged observations, with intensities that decay exponentially with lag distance. This multi-lag specification preserves the interpretability of the one-lag model while incorporating residual longer-memory information. Its performance improves relative to the one-lag version, with the clearest gains appearing in the signed-volume equation, but the gain saturates rapidly as the maximum lag increases. This is consistent with both the VAR lag analysis and the SHAP lag-importance profiles: longer histories contain some additional predictive information, but the dominant nonlinear contribution is concentrated at short horizons. Moreover, the structured parametric lag decay provides a more stable dependence on the lag length than the unconstrained neural-network benchmark. Overall, the paper proposes a constructive use of explainability in financial modeling. Starting from a flexible nonlinear model calibrated on data, we use SHAP values and aggregated conditional response surfaces analysis to extract isolated, economically meaningful dependencies among the regressors. These dependencies are then translated into interpretable parametric models that retain most of the predictive power of the neural network while offering greater transparency and stability. The approach has two complementary benefits. First, it provides evidence that the neural network is learning meaningful financial mechanisms rather than spurious statistical patterns. Second, it turns these mechanisms into explicit parametric forms that can be calibrated, interpreted, compared with standard benchmarks, and potentially used in simulation or structural market-microstructure models. Several directions remain open. First, the parametric model could be extended to include richer state variables from the limit order book, such as queue imbalance, depth, spread regimes, or measures of liquidity at the best quotes. Second, the framework could be tested on longer time periods, additional assets, and different market regimes to assess the stability of the extracted nonlinear mechanisms. These extensions would further clarify how machine-learning explainability can be used to build structural models of price formation in modern electronic markets. Acknowledgements This publication has received funding from the European Union – NextGenerationEU – PNRR – Project: “SoBigData.it – Strengthening the Italian RI for Social Mining and Big Data Analytics” – Prot. IR0000013. The research was carried out also with the contribution of research groups linked to the Horizon 2020 Program under the scheme “INFRAIA-01-2018-2019 – Integrating Activities for Advanced Communities”, Grant Agreement n.871042, “SoBigData++: European Integrated Infrastructure for Social Mining and Big Data Analytics”, to the HORIZON Europe Program under the scheme “INFRA-2021-DEV-02- SoBigData RI Preparatory Phase Project” Grant Agreement n. 101079043 — SoBigData RI P. References [1] J. Bouchaud, Y. Gefen, M. Potters, and M. Wyart (2004) Fluctuations and response in financial markets: the subtle nature of random price changes. Quantitative Finance 4 (2), p. 176–190. Cited by: §1, §1, §1, §2.1. [2] F. Lillo and J. D. Farmer (2004) Long memory in the efficient market. Studies in Nonlinear Dynamics & Econometrics 8 (3). Cited by: §1. [3] B. Tóth, Y. Lemperiere, C. Deremble, J. de Lataillade, J. Kockelkoren, and J. Bouchaud (2011) Anomalous price impact and the critical nature of liquidity in financial markets. Physical Review X 1 (2), p. 021006. Cited by: §1. [4] R. Almgren and N. Chriss (2001) Optimal execution of portfolio transactions. Journal of Risk 3, p. 5–40. Cited by: §1. [5] A. Bucci, I. Mastromatteo, M. Benzaquen, and J. Bouchaud (2019) Are trading invariants really invariant? trading costs matter. Quantitative Finance. Cited by: §1. [6] J. Hasbrouck (1991) Measuring the information content of stock trades. Journal of Finance 46 (1), p. 179–207. Cited by: Appendix A, §1, §2.1. [7] J. Sirignano and R. Cont (2019) Deep learning for limit order books. Quantitative Finance 19 (4), p. 549–570. Cited by: §1. [8] Z. Zhang, S. Zohren, and S. Roberts (2019) DeepLOB: deep convolutional neural networks for limit order books. IEEE Transactions on Signal Processing 67 (11), p. 3001–3012. Cited by: §1. [9] S. M. Lundberg and S. Lee (2017) A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems 30. Cited by: §1, §2.2, §5.1, §5.1. [10] M. Naviglio, G. Bormetti, F. Campigli, G. Rodikov, and F. Lillo (2026) Why is the estimation of metaorder impact with public market data so challenging?. Quantitative Finance 26 (4), p. 505–523. External Links: Document, Link, https://doi.org/10.1080/14697688.2026.2628221 Cited by: §2.1, §3. [11] G. Cybenko (1989) Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems 2 (4), p. 303–314. External Links: Document Cited by: §2.2, §2.2. [12] K. Hornik, M. Stinchcombe, and H. White (1989) Multilayer feedforward networks are universal approximators. Neural Networks 2 (5), p. 359–366. External Links: Document Cited by: §2.2, §2.2. [13] I. Goodfellow, Y. Bengio, and A. Courville (2016) Deep learning. MIT Press. External Links: Link Cited by: §2.2, §2.2, §2.2. [14] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner (1998) Gradient-based learning applied to document recognition. Proceedings of the IEEE 86 (11), p. 2278–2324. External Links: Document Cited by: §2.2. [15] S. Hochreiter and J. Schmidhuber (1997) Long short-term memory. Neural Computation 9 (8), p. 1735–1780. External Links: Document Cited by: §2.2. [16] D. P. Kingma and J. Ba (2014) Adam: a method for stochastic optimization. External Links: 1412.6980 Cited by: §2.2. [17] J. H. Friedman (2001) Greedy function approximation: a gradient boosting machine. The Annals of Statistics 29 (5), p. 1189–1232. External Links: Document Cited by: §2.2. [18] Z. Eisler, J. Bouchaud, and J. Kockelkoren (2012) The price impact of order book events: market orders, limit orders and cancellations. Quantitative Finance 12 (9), p. 1395–1419. External Links: Document, 0904.0900 Cited by: §3, §5.3.1. [19] A. Shrikumar, P. Greenside, and A. Kundaje (2017) Learning important features through propagating activation differences. In Proceedings of the 34th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 70, p. 3145–3153. Cited by: §5.1. [20] F. Lillo, J. D. Farmer, and R. N. Mantegna (2003) Master curve for price-impact function. Nature 421 (6919), p. 129–130. External Links: Document Cited by: §5.3.1. [21] J. Bouchaud, J. D. Farmer, and F. Lillo (2009) How markets slowly digest changes in supply and demand. Handbook of Financial Markets: Dynamics and Evolution, p. 57–160. External Links: Document Cited by: §5.3.1. [22] J. Bouchaud, Y. Gefen, M. Potters, and M. Wyart (2004) Fluctuations and response in financial markets: the subtle nature of ‘random’ price changes. Quantitative Finance 4 (2), p. 176–190. External Links: Document Cited by: §5.3.1. [23] B. Tóth, I. Palit, F. Lillo, and J. D. Farmer (2015) Why is equity order flow so persistent?. Journal of Economic Dynamics and Control 51, p. 218–239. External Links: Document Cited by: §5.3.1. [24] R. Cont, A. Kukanov, and S. Stoikov (2014) The price impact of order book events. Journal of Financial Econometrics 12 (1), p. 47–88. External Links: Document Cited by: §5.3.1. [25] L. Kilian and H. Lütkepohl (2017) Nonlinear structural VAR models. In Structural Vector Autoregressive Analysis, p. 609–658. External Links: Document Cited by: Appendix A, Appendix A. [26] F. Cordoni, N. Doremus, and A. Moneta (2024) Identification of vector autoregressive models with nonlinear contemporaneous structure. Journal of Economic Dynamics and Control 161, p. 104846. External Links: Document Cited by: Appendix A. Appendix A Residual Contemporaneous Dependence and Structural Identification In this Appendix we provide additional details on the residual structural interpretation used in the main text. The key point is that the reduced-form models considered in the paper use lagged information only. Therefore, contemporaneous effects are not included as ordinary regressors, but appear as dependence between reduced-form residuals. The reduced-form model can be written as rt=fr(t)+utr,vt=fv(t)+utv, splitr_t&=f_r(X_t)+u_t^r,\\ v_t&=f_v(X_t)+u_t^v, split (29) where tX_t denotes the lagged information set, while utru_t^r and utvu_t^v are the reduced-form innovations. The Hasbrouck-type structural interpretation is obtained by imposing a recursive contemporaneous ordering. In trade time, the signed volume is observed at the event and can be regarded as causally prior to the associated quote revision. This motivates vt=fv(t)+ϵtv,rt=fr(t)+ψ(ϵtv)+ϵtr. splitv_t&=f_v(X_t)+ _t^v,\\ r_t&=f_r(X_t)+ψ( _t^v)+ _t^r. split (30) In the linear Hasbrouck model, ψ(ϵtv)=b0ϵtvψ( _t^v)=b_0 _t^v. Hence, utv=ϵtv,utr=b0utv+ϵtr.u_t^v= _t^v, u_t^r=b_0u_t^v+ _t^r. (31) Thus, the contemporaneous impact coefficient is recovered from the residual dependence between utru_t^r and utvu_t^v, rather than by inserting vtv_t directly into the reduced-form forecasting equation. This is the usual logic of structural VAR identification: reduced-form residuals are correlated because they are mixtures of underlying structural shocks, and additional restrictions are required to recover those shocks [6, 25]. Directly fitting the return equation by inserting vtv_t as an ordinary contemporaneous regressor is generally not the cleanest approach when the goal is to separate lagged predictability from instantaneous impact. In the reduced-form formulation, the model estimates fr(t)≃[rt∣t],fv(t)≃[vt∣t],f_r(X_t) [r_t _t], f_v(X_t) [v_t _t], and the residuals are the corresponding innovations, utr=rt−[rt∣t],utv=vt−[vt∣t].u_t^r=r_t-E[r_t _t], u_t^v=v_t-E[v_t _t]. Thus, utvu_t^v represents the unexpected component of order flow, namely the part of vtv_t that is not predictable from past returns and volumes. By contrast, if one fits rt=G(t,vt)+ηt,r_t=G(X_t,v_t)+ _t, then the model receives as input vt=[vt∣t]+utv.v_t=E[v_t _t]+u_t^v. Hence, the return equation receives both the predictable component of order flow, already determined by lagged information, and the genuinely contemporaneous innovation utvu_t^v. These two effects are mixed inside the same function G. This can lead to an incorrect interpretation of the learned dependence: what appears to be an instantaneous effect of vtv_t may partly come from its predictable lagged component, while the true contemporaneous innovation is not isolated. In a linear model this mixing is less problematic, because the predictable component can be reabsorbed into the reduced-form lagged predictor. Suppose, for example, that the structural model is vt=fv(t)+ϵtv,rt=gr(t)+b0vt+ϵtr. splitv_t&=f_v(X_t)+ _t^v,\\ r_t&=g_r(X_t)+b_0v_t+ _t^r. split (32) Substituting the first equation into the second gives rt=gr(t)+b0fv(t)+b0ϵtv+ϵtr.r_t=g_r(X_t)+b_0f_v(X_t)+b_0 _t^v+ _t^r. (33) The predictable component fv(t)f_v(X_t) is absorbed into the lagged reduced-form predictor, while the innovation ϵtv=utv _t^v=u_t^v remains in the residual component. Therefore, in the linear case, the direct contemporaneous formulation and the reduced-form residual formulation are tightly connected. This equivalence is much less transparent in nonlinear models. If one directly fits rt=G(t,vt)+ηt,r_t=G\! (X_t,v_t )+ _t, (34) then, since vt=fv(t)+utvv_t=f_v(X_t)+u_t^v, the nonlinear function G acts jointly on the predictable component of order flow and on its contemporaneous innovation: G(t,fv(t)+utv).G\! (X_t,f_v(X_t)+u_t^v ). In general, this expression cannot be decomposed uniquely into a purely lagged term plus a purely contemporaneous innovation term. As a result, direct conditioning on vtv_t may mix lagged predictable order-flow effects with instantaneous impact. This is precisely the effect that the residual decomposition avoids. For this reason, especially in nonlinear models, we estimate first the reduced-form predictive model using only lagged information, xt=fθ(t)+ut,xt=(rt,vt),x_t=f_θ(X_t)+u_t, x_t=(r_t,v_t), (35) and then analyze the residual dependence. Even when fθf_θ is a single vector-valued nonlinear map, the residual vector can be written componentwise as ut=(utrutv)=(rt−fr(t)vt−fv(t)).u_t= pmatrixu_t^r\\ u_t^v pmatrix= pmatrixr_t-f_r(X_t)\\ v_t-f_v(X_t) pmatrix. A residual dependence between utru_t^r and utvu_t^v indicates that, after removing the lagged predictable dynamics, an instantaneous relation remains between unexpected order flow and unexpected price changes. In the most general nonlinear structural representation, one could write the reduced-form residuals as nonlinear functions of structural shocks, ut=Ht(ϵtr,ϵtv).u_t=H_t( _t^r, _t^v). (36) Recovering the shocks then requires identifying restrictions on HtH_t. This is precisely the point emphasized by nonlinear structural VAR models: when contemporaneous relations are nonlinear, the mapping from reduced-form innovations to structural shocks is no longer characterized by a constant impact matrix, and identification requires additional assumptions [25, 26]. In our market microstructure setting, we impose the same recursive economic ordering used in the Hasbrouck framework: order-flow innovations are allowed to affect return innovations contemporaneously, but not vice versa within the same trade-time event. This leads to utv=ϵtv,utr=ϕ(utv)+ϵtr.u_t^v= _t^v, u_t^r=φ(u_t^v)+ _t^r. (37) The function ϕφ captures the contemporaneous impact of the unexpected component of signed volume on the unexpected component of returns. This decomposition also clarifies the distinction between reduced-form forecasting and structural conditioning. If one is interested in pure forecasting before the event at time t, the relevant prediction is r^tlagged=fr(t), r_t^\,lagged=f_r(X_t), which uses only information available up to t−1t-1. This is the object used when comparing the predictive performance of the VAR, the nonlinear parametric model, and the neural network. By contrast, if the signed volume vtv_t is observed before the associated quote revision, one can form r^tstruct=fr(t)+ϕ^(vt−fv(t)), r_t^\,struct=f_r(X_t)+ φ\! (v_t-f_v(X_t) ), which is a prediction conditional on the contemporaneous realization of order flow, rather than a purely lagged forecast. Appendix B Lag-order dependence of forecasting performance Figure 10: Predictive performance of VAR model of Eq. (3) as a function of the lag order p for BAC. The panels report the in-sample and out-of-sample (R2R^2) as functions of the lag order. BAC is shown as a representative example, since the same qualitative behavior is observed across the other stocks in the dataset. In this appendix we provide a more detailed discussion of the dependence of the forecasting performance on the lag order p. The purpose of this analysis is not to identify an asset-specific optimal lag length, but rather to verify whether increasing the memory of the model leads to systematic predictive gains. This serves as a robustness check for the choice p=20p=20 adopted in the main text. For each value of p, the models are trained using the lagged state vector Xt(p)=(rt−1,vt−1,rt−2,vt−2,…,rt−p,vt−p),X_t^(p)= (r_t-1,v_t-1,r_t-2,v_t-2,…,r_t-p,v_t-p ), (38) where rtr_t denotes returns and vtv_t denotes signed volumes. We then evaluate the in-sample and out-of-sample R2R^2 for the prediction of the target variables. This allows us to compare how the amount of past information included in the model affects predictive performance. Figure 10 reports the lag-dependence for the VAR benchmark in the case of the BAC stock. The behavior is rather regular. At small lag orders, increasing p substantially improves both in-sample and out-of-sample performance, since the model is progressively allowed to exploit short-term linear dependencies in the joint dynamics of returns and signed volumes. After the first few lags, however, the out-of-sample R2R^2 reaches a relatively stable level and further increases in p do not lead to systematic improvements. In some cases, the out-of-sample performance even slightly decreases at larger lag orders. This behavior is natural for a linear VAR model. Increasing the lag order amounts to adding additional linear regressors to the specification. Once the dominant short-lag dependencies have been included, the marginal predictive content of further lagged variables becomes small. At the same time, the number of parameters to estimate increases with p. Longer lag structures may therefore introduce additional estimation noise without providing a compensating increase in useful predictive information. The resulting pattern is an initial improvement followed by a saturation of the out-of-sample performance. A similar qualitative behavior is observed for the other stocks in the sample. The BAC results are therefore reported as a representative illustration of the general lag-dependence pattern, rather than as an asset-specific finding. The neural-network model displays a related, but less regular, behavior, as shown in Figure 11. In this case, increasing the lag order does not produce a smooth monotonic improvement of the test R2R^2. After the first lags, the performance fluctuates around a comparable range of values, with local increases and decreases depending on the specific lag length. These oscillations are more pronounced than in the VAR case and should not be interpreted as evidence that particular lag orders are structurally privileged. Figure 11: Lag-dependent predictive performance for AAPL using the neural network model. The left panel reports the in-sample and out-of-sample (R2R^2) for returns as a function of the lag order, while the right panel shows the corresponding results for signed volumes. AAPL is shown as a representative example, as similar qualitative patterns are observed across all stocks. Figure 12: Normalized lag importance for large-tick stocks (LT). Each curve represents the average fraction of total Shapley magnitude explained by lag ℓ . Figure 13: Normalized lag importance for small-tick stocks (ST). Each curve represents the average fraction of total Shapley magnitude explained by lag ℓ . The reason is that changing p in the neural-network setting does more than simply adding a few additional coefficients. For each lag order, the model is trained on a different input representation Xt(p)X_t^(p). Therefore, the network does not only receive a longer history, but also explores a different nonlinear mapping from the lagged state variables to the target variable. Additional distant lags may contain weak residual information, especially for signed volumes, but they can also act as noisy inputs. Because the model is flexible, such weak signals can interact nonlinearly with the most recent variables and with finite-sample fluctuations. This can generate a noisy dependence of the test performance on p, even when there is no genuine systematic gain from using longer histories. This distinction also explains why the VAR curves are typically smoother than the neural-network curves. In the linear case, the effect of increasing p is more controlled: the model class changes by adding linear terms. In the nonlinear case, instead, the effective model explored by training can vary more substantially across lag orders, because the network can combine the additional inputs through nonlinear interactions. Hence, a larger lag order does not necessarily imply a better forecast, even if the training set contains more explanatory variables. Overall, the lag-dependence analysis supports two conclusions. First, the relevant predictive memory of the system is concentrated at relatively short horizons: after the initial lags, increasing p does not provide robust out-of-sample improvements. Second, the noisy oscillations observed for the neural network are consistent with the interaction between weak long-lag signals, nonlinear model flexibility and finite-sample estimation effects, rather than with the existence of sharply optimal lag lengths. On the basis of these results, we use p=20p=20 in the main analysis. This choice is conservative: it is large enough to include the short-memory effects that drive the initial improvement in predictive performance, but not so large as to excessively increase the input dimension and the number of parameters. Moreover, the absence of systematic gains at larger lag orders indicates that the main empirical conclusions are not driven by a fine tuning of the memory length. Appendix C Shapley Lag-Importance Profiles To quantify how the explanatory contribution of the neural network is distributed across lags, we construct an aggregated measure of lag importance based on mean absolute Shapley values. For each ticker i and each prediction block, for example r→rtr→ r_t, v→rtv→ r_t, r→vtr→ v_t, or v→vtv→ v_t, we compute Ii(ℓ)=[|ϕjℓf(xt)|],I_i( )=E [| _j_ ^f(x_t)| ], (39) where jℓj_ denotes the input feature corresponding to lag ℓ , either rt−ℓr_t- or vt−ℓv_t- , and the expectation is taken over all observations for that ticker. Since raw Shapley magnitudes are not directly comparable across tickers and outputs, we normalize them within each ticker and block: I~i(ℓ)=Ii(ℓ)∑k=1pIi(k). I_i( )= I_i( ) _k=1^pI_i(k). (40) By construction, ∑ℓ=1pI~i(ℓ)=1, _ =1^p I_i( )=1, so that I~i(ℓ) I_i( ) can be interpreted as the fraction of the total Shapley magnitude attributed to lag ℓ . We then average these profiles across tickers, separately for large-tick and small-tick stocks: I¯(ℓ)=1N∑i=1NI~i(ℓ). I( )= 1N _i=1^N I_i( ). (41) Figures 12 and 13 report the resulting profiles for LT and ST stocks. The first lag dominates across all prediction blocks and both tick-size regimes. This suggests that most of the nonlinear structure learned by the neural network is concentrated in the most recent return–volume pair (rt−1,vt−1)(r_t-1,v_t-1). At the same time, the decay is not identical across channels. Lagged signed volumes generally display a more persistent profile than lagged returns, especially in the v→vtv→ v_t block, consistently with the persistence of order flow. Lagged returns instead tend to decay more rapidly, indicating that their predictive contribution is mainly concentrated at very short horizons. The comparison between LT and ST stocks also shows some differences. In LT stocks, volume lags retain a relatively more gradual contribution in the v→rtv→ r_t and v→vtv→ v_t blocks. In ST stocks, the effect of lagged volumes on returns is more concentrated at the first lag, while lagged returns contribute somewhat more persistently to future signed volumes. This suggests that the interaction between past price changes and subsequent trading activity is richer in small-tick assets. Overall, these results justify focusing the main predictive-surface analysis on the first lag, while also motivating the multi-lag parametric extension with decaying lag weights introduced in the main text. Appendix D Empirical surfaces and validation The binning procedure depends on the variable placed on the horizontal axis. When the response is plotted as a function of vt−1v_t-1, the support of signed volumes is partitioned into uniform bins within each ticker. Since these bins are defined over ticker-specific ranges, they are not directly comparable across assets. When the response is plotted as a function of rt−1r_t-1, uniform binning is less appropriate, as returns are concentrated around a finite number of clusters due to price discreteness. In this case, a Gaussian kernel density estimator is used to identify local maxima of the empirical distribution, and bins are defined through the intervals between consecutive peaks. This results in a data-driven clustering of returns. Each ticker-specific curve is then interpolated linearly on a common grid, separately for each conditioning regime. At each grid point, the cross-sectional average of the model prediction is computed across tickers. To avoid distortions from sparsely populated regions, only bins with at least 150 average observations are retained in the final plots, while the histograms are kept complete. The shaded bands represent cross-sectional uncertainty and are constructed as ± one standard error of the mean across tickers. These bands capture the dispersion of the estimated relationship across assets and should not be interpreted as confidence intervals at the observation level. The results of the previous analysis suggest that the predictive relationships captured by the model are strongly nonlinear and state-dependent. In this section, we want to validate these results by tracing them back to empirical regularities directly observed in the data. Of course, when relying only on the observed data, we can recover only the empirical analogues of the aggregated conditional response surfaces defined in Eqs. (13) and (14). What is lost is the additive attribution of the prediction across the different input variables, which is provided by the Shapley decomposition. In this sense, a well-calibrated model that captures genuine regularities in the data can be used not only for prediction, but also as an interpretative device: it allows us to decompose the observed predictive structure and to assess the relative contribution of returns, signed volumes, and the remaining lagged variables. The empirical results are constructed starting from the time series of returns rtr_t and signed volumes vtv_t, considered separately for large-tick (LT) and small-tick (ST) stocks. We then study bivariate relationships between a response variable at time t and a lagged explanatory variable at time t−1t-1, considering the four configurations: (rt vs rt−1),(rt vs vt−1),(vt vs rt−1),(vt vs vt−1).(r_t vs r_t-1), (r_t vs v_t-1), (v_t vs r_t-1), (v_t vs v_t-1). In addition to unconditional relationships, we also consider conditional ones, splitting the sample according to the sign of the second lagged variable. In particular: • when conditioning on rt−1r_t-1, we consider the regimes rt−1<0r_t-1<0, rt−1=0r_t-1=0, and rt−1>0r_t-1>0; • when conditioning on vt−1v_t-1, we consider the regimes vt−1<0v_t-1<0 and vt−1>0v_t-1>0. For each ticker, the explanatory variable is discretized into bins and, within each bin, we compute the empirical mean of the response variable. The resulting curves represent bin-wise averages of the data. To obtain bucket-level results (LT and ST), the ticker-specific curves are interpolated onto a common grid and averaged cross-sectionally, using the same procedure adopted in the previous sections. The final plots display the average curve, a confidence band based on the standard error, and histograms representing the average number of observations per bin. Figure 14: Empirical nonlinear relationships for large-tick (LT) stocks. The empirical relationships for large-tick (LT) and small-tick (ST) stocks are shown in Figures 14 and 15. Overall, the agreement between the model-implied aggregated conditional response surfaces and the empirical ones reconstructed directly from the data provides an important sanity check. It indicates that the nonlinear model is capturing dependence patterns that are genuinely present in the data, rather than introducing spurious nonlinearities through its flexibility. This validation is important for the interpretation developed in the previous sections: once the model is shown to reproduce the empirical predictive structure, we are more confident that the Shapley decomposition offers a meaningful way to separate the contributions of individual lagged variables and to identify which components drive the observed relationships. Figure 15: Empirical nonlinear relationships for small-tick (ST) stocks. Appendix E VAR-implied predictive surfaces In this appendix we repeat the analysis of the model-implied predictive surfaces using the linear VAR benchmark instead of the neural-network model. The goal is to provide a direct comparison between the nonlinear response functions learned by the DNN and the corresponding structures implied by a purely linear specification. Starting from the VAR model in Eq. (3), the one-step-ahead predictions for returns and signed volumes are linear functions of the lagged return–volume state. Therefore, when the prediction is viewed as a function of (rt−1,vt−1)(r_t-1,v_t-1), the resulting predictive surface is affine in the current lagged variables, up to the contribution of the remaining lags. In particular, the dependence on vt−1v_t-1 is controlled by the first-lag volume coefficients, whereas the dependence on rt−1r_t-1 is controlled by the first-lag return coefficients. This makes the VAR surfaces a natural benchmark for assessing which patterns can already be explained by linear lagged dependencies and which ones require the additional flexibility of the neural network. We construct the binned VAR surfaces using the same procedure adopted in the main text for the DNN. Namely, we compute conditional averages of the VAR predictions within bins of either vt−1v_t-1 or rt−1r_t-1, further conditioning on the sign regime of the other lagged variable. The same binning and clustering procedure is used, so that the comparison with the DNN surfaces is not affected by differences in empirical support. The resulting surfaces are reported in Figures 16 and 17, respectively for small-tick and large-tick stocks. Figure 16: VAR-implied binned predictive surfaces for small-tick (ST) stocks. The panels are organized as in Figure 3: predicted returns and signed volumes are represented as functions of either rt−1r_t-1 or vt−1v_t-1, conditional on the sign regime of the other lagged variable. Figure 17: VAR-implied binned predictive surfaces for large-tick (LT) stocks. The panels are organized as in Figure 4: predicted returns and signed volumes are represented as functions of either rt−1r_t-1 or vt−1v_t-1, conditional on the sign regime of the other lagged variable. Figures 16 and 17 show that the VAR captures the main first-order directional effects in the data. In particular, when the response is represented as a function of vt−1v_t-1, the predicted signed volume is approximately monotone increasing: positive lagged signed volume predicts positive future signed volume, while negative lagged signed volume predicts negative future signed volume. This is the linear counterpart of order-flow persistence. Similarly, the predicted return tends to increase with vt−1v_t-1, consistently with the idea that buy-initiated volume is associated with positive short-horizon price changes and sell-initiated volume with negative ones. However, because the VAR is linear, these relations are necessarily smoother and more constrained than those obtained with the neural network. The sharp transitions around vt−1=0v_t-1=0, the flattening at large absolute values of signed volume, and the more pronounced state dependence observed in the DNN surfaces cannot be fully reproduced by the linear benchmark. In the VAR case, the dependence on vt−1v_t-1 is essentially governed by a single linear slope, and the differences between the curves associated with rt−1<0r_t-1<0, rt−1=0r_t-1=0, and rt−1>0r_t-1>0 mainly appear as level shifts induced by the conditioning variable and by correlations with the remaining regressors. A similar interpretation applies to the panels where the horizontal axis is rt−1r_t-1. The VAR produces approximately linear responses in lagged returns, with the curves separated according to the sign of vt−1v_t-1. This confirms that lagged signed volume remains an important state variable also in the linear model. Nevertheless, the role of rt−1r_t-1 is constrained by the linear structure: it can only affect the prediction proportionally. As a consequence, the VAR cannot reproduce the more flexible attenuation patterns observed in the neural-network surfaces, where large absolute lagged returns tend to weaken the continuation of the order-flow signal. The comparison across tick-size regimes is also informative. For small-tick stocks, shown in Figure 16, the VAR recovers the average direction of the main relationships, but it provides a substantially smoother representation than the DNN. This is especially visible in the response to lagged signed volume, where the neural network displays sharper nonlinearities and stronger state dependence. For large-tick stocks, shown in Figure 17, the difference between the VAR and the DNN is less pronounced. This is consistent with the idea that large-tick dynamics are more strongly constrained by price discreteness, so that a linear benchmark already captures a relevant part of the directional dependence between recent order-flow pressure and price changes. Overall, the VAR predictive surfaces confirm that part of the return–volume dynamics can be described by linear lagged dependencies. In particular, the linear model captures the sign of price impact and the persistence of signed order flow. At the same time, the comparison with the neural-network surfaces shows that the DNN is not merely reproducing the VAR relation with additional noise. Rather, it learns a richer response characterized by sharper transitions, saturation effects, and state-dependent modulation of the predictive signal. This provides a useful benchmark for the explainability analysis in the main text: deviations from the VAR surfaces can be interpreted as evidence of nonlinear dependencies learned by the neural network beyond the linear specification. Appendix F Large tick Shapley Components In this appendix, we report the Shapley components for large-tick stocks, following the same structure used for the small-tick stocks in the main text. The corresponding profiles are shown in Figures 18 and 19. Overall, the results for large-tick stocks are qualitatively similar to those obtained for small-tick stocks, although some differences are relevant and deserve to be discussed. We start from the term associated with the contribution of vt−1v_t-1 to the next-step prediction. The behavior is essentially analogous to the one observed for small-tick stocks. In particular, the volume term displays a sign-like dependence on vt−1v_t-1: positive lagged signed volume gives a positive contribution, while negative lagged signed volume gives a negative contribution. Moreover, the conditional curves corresponding to non-zero values of rt−1r_t-1 are, within the error bars, almost collapsed onto the curve associated with rt−1=0r_t-1=0. Therefore, also in the large-tick case, we do not observe an additional dependence on rt−1r_t-1 in the volume contribution. This conclusion is confirmed by the starred terms, shown in the right panels of Figures 18, which are essentially flat as functions of rt−1r_t-1. Hence, the dependence on lagged signed volume appears to be genuinely captured by the volume Shapley term, without a relevant spurious component attributed to the return term. We then consider the return-related contribution. Also in this case, as for small-tick stocks, a relevant part of the structure is mediated by the lagged signed volume. Looking at the starred terms in the right panels of Figures 19, we recover the same qualitative mechanism observed in the ST bucket. When the previous trade did not induce a price change, namely rt−1=0r_t-1=0, the return term contributes to the next-step return and volume by preserving the sign of the previous signed volume. By contrast, when a price change occurred, rt−1≠0r_t-1≠ 0, the average contribution tends to act in the opposite direction, generating a pressure toward a sign reversal both in the next signed volume and in the next return. This behavior is also visible in the non-starred return profiles. For rt−1=0r_t-1=0, the contribution of the return term tends to preserve the sign of vt−1v_t-1, as in the small-tick case. However, differently from the ST bucket, the LT profiles display a more pronounced asymmetry with respect to rt−1r_t-1. For rt−1>0r_t-1>0, the return contribution to the next-step volume and return tends to have the opposite sign, while for rt−1<0r_t-1<0 the opposite pattern is observed. This suggests a stronger mean-reverting component in prices for large-tick stocks. A similar effect was already visible, although in a weaker form, for small-tick stocks in the profile srrs_r^r, reported in the top-left panel of Figure 6. Figure 18: Average Shapley profiles for large tick (LT) stocks for the volume term. Top panels refer to the return prediction r^t r_t, while bottom panels refer to the signed-volume prediction v^t v_t. Left panels show the contribution of the conditioning variable vt−1v_t-1 when the horizontal axis is vt−1v_t-1. Right panels show the contribution of the conditioning variable vt−1v_t-1 when the horizontal axis is rt−1r_t-1. Shaded bands denote ± one cross-sectional standard error, while histograms report the average empirical support. Figure 19: Average Shapley profiles for large tick (LT) stocks for the return term. Panels are organized as in Figure 18. Appendix G Lag dependence for the parametric multi-lag model (a) Returns. (b) Signed volumes. Figure 20: Evolution of the in-sample and out-of-sample R2R^2 as a function of the maximum lag p for the GILD ticker. The improvement after the first lag is limited and the predictive performance rapidly saturates, consistently with a short-memory nonlinear dependence structure. Figure 20 reports the evolution of the R2R^2 as a function of the maximum lag p for the representative case of GILD. For both returns and signed volumes, the largest improvement occurs when moving beyond the one-lag specification, while the curves flatten as p increases. This behaviour is consistent with the evidence obtained from the linear VAR benchmark and with the SHAP lag-importance analysis: the process has a short dominant nonlinear memory, whereas more distant lags contribute only marginally. The exponential-decay parametrization therefore provides a useful compromise between flexibility and parsimony, since it incorporates the residual information contained in longer histories without assigning independent nonlinear parameters to each lag. Appendix H Structural shocks extraction for the parametric multi-lag model Figure 21 reports the estimated functions ϕ^(utv) φ(u_t^v), averaged separately over large-tick and small-tick stocks. In both cases, we observe similar results as the neural network model. The scale of the relation is larger for large-tick stocks, consistently with the stronger discreteness of price changes and the queue-driven nature of price formation in the large-tick regime. Table 6 reports the residual correlations before and after the structural correction. Before the correction, the correlation between utru_t^r and utvu_t^v is sizeable, with average values around 0.290.29 for large-tick stocks and 0.320.32 for small-tick stocks. After subtracting the fitted binned projection ϕ^(utv) φ(u_t^v), the correlation is reduced to approximately 10−210^-2 in both buckets. The same pattern is observed for Spearman correlations, which fall from about 0.370.37 to values close to zero for large-tick stocks, and from about 0.420.42 to slightly negative values for small-tick stocks. This indicates that the binned structural projection removes most of the contemporaneous signed dependence between return and signed-volume residuals. The last two columns report additional diagnostics based on nonlinear transformations of the residuals. The correlation with squared residuals remains very small on average, while the correlation with absolute residuals is larger but still moderate, especially for large-tick stocks. Hence, the structural correction removes the main directional dependence, although a residual volatility-level dependence may remain. (a) Large tick. (b) Small tick. Figure 21: Estimated structural projection ϕ^(utv) φ(u_t^v) for the multi-lag parametric model. The function is estimated through binned conditional averages and then averaged separately over large-tick and small-tick stocks. Table 6: Residual structural decomposition for the multi-lag parametric model using the binned estimator of ϕφ. The last row of each bucket reports equal-weight averages across tickers. Bucket Ticker Corr. before Corr. after Spearman before Spearman after Corr. sq. Corr. abs. LT BAC 0.301 0.0121 0.361 -0.0214 0.0034 0.111 LT CMCSA 0.250 0.0083 0.294 -0.0347 -0.0123 0.049 LT CSCO 0.292 0.0109 0.348 -0.0221 0.0075 0.107 LT INTC 0.271 0.0119 0.323 0.0027 0.0175 0.126 LT PFE 0.331 0.0135 0.423 0.0136 -0.0152 0.086 LT Average 0.289 0.0113 0.369 -0.0108 0.0002 0.096 ST AAPL 0.309 0.0078 0.352 -0.0516 0.0064 0.051 ST AMZN 0.382 0.0118 0.477 -0.0164 -0.0196 0.001 ST GILD 0.340 0.0077 0.422 -0.0325 -0.0143 0.020 ST NVDA 0.281 0.0149 0.385 -0.0223 0.0018 0.045 ST TSLA 0.271 0.0100 0.345 -0.0425 0.0054 0.039 ST Average 0.317 0.0104 0.417 -0.0322 -0.0041 0.031