3.2.1. Linear Granger Causality Test

Under framework (2), linear causality is typically tested via VAR. For two stationary series \(X_t\) and \(Y_t\), the null hypothesis that \(X\) does not Granger-cause \(Y\) is tested based on an \(F\)-statistic comparing restricted and unrestricted VAR models ¹,²,³. Although simple and intuitive, linear tests cannot capture non-linear dependencies.

3.2.2. Nonlinear Granger Causality Tests

Nonparametric and semiparametric tests were conducted to capture nonlinearities. The DP test ⁷ uses kernel density estimation to test the conditional independence in Eq. (1) directly, thereby overcoming the over-rejection issue of earlier methods. The LASSO-based test ⁸ approximates an unknown function via a Taylor expansion and uses LASSO for variable selection, thus enabling the identification of specific nonlinear patterns (e.g., quadratic terms). Quantile causality tests ⁹,¹⁰,¹¹ examine the predictability at specific quantiles and reveal asymmetric risk spillovers.

3.2.3. Tests Limitations of Existing Methods

Despite these advances, classical tests exhibit several limitations when applied to price discovery in futures and spot markets:

1. (1)

   Focus on conditional moments: Most tests compare conditional expectations or specific quantiles while disregarding the entire distribution. Consequently, volatility spillovers and tail dependencies are overlooked.

2. (2)

   In-sample fit vs. out-of-sample prediction: They rely on in-sample goodness-of-fit (e.g., RSS comparisons), which is not equivalent to out-of-sample predictive ability-the essence of Granger causality ¹².

3. (3)

   Non-integration with time-varying volatility: Nonlinear mean equations are commonly estimated separately from volatility dynamics, thus yielding incomplete conditional distribution forecasts.

4. (4)

   Implicit stability assumption: Causality is assumed to be constant over time, whereas financial markets undergo structural changes.

3.3.1. Overall Logic and Rolling Forecast Mechanism

We return to Granger’s original definition ¹² and compare the entire conditional distribution of futures returns with and without another market’s information. A rolling-window scheme is used to account for time variation. For each window, we estimate the following two models:

1. (1)

   Benchmark model: Uses only own-market lags;

2. (2)

   Extended model: Incorporates the lags (linear and quadratic) of another market.

One-step-ahead distribution forecasts are generated, and the out-of-sample sequences of PITs and logarithmic scores are evaluated.

3.3.2. Nonlinear Mean Equation: LASSO Variable Selection

Let \(Y_t\) be the return series (spot \(S_t\) or futures \(F_t\)). The conditional mean \(\mu_{t} = E(Y_{t}\mid I_{t-1})\) is modeled as a function of lagged linear and quadratic terms (up to lag \(p\) order). To avoid overfitting and identify relevant predictors, we apply LASSO to each rolling window as follows:

3.3.3. Time-Varying Volatility: EGARCH(1,1) and Skewed-\(t\) Distribution

To capture volatility clustering and asymmetry (the leverage effect), we use EGARCH(1,1) for the conditional variance as follows:

Stock-market volatility typically exhibits clustering, which justifies the use of GARCH-type models. However, when the actual data lack asymmetric effects, \(\gamma\) is insignificant. Consequently, in the Monte Carlo simulations presented in Section 4, the volatility process is specified as a GARCH model instead of an EGARCH model.

3.3.4. Distributional Forecast Evaluation: HL and AG Tests

After obtaining one-step-ahead forecast distributions from the benchmark and extended models over all the rolling windows, we apply two evaluation criteria:

1. (1)

   PIT and HL test

   For each forecast cumulative distribution function \(F_{t\mid t-1}\), we compute

   \begin{equation} \mathit{PIT}_{t} = F_{t\mid t-1}\left(y_{t}\right). \label{eq:5} \tag{5} \end{equation}

   If the forecast distribution is correctly specified, then \(\mathit{PIT}_t\) should be i.i.d. uniform in \([0,1]\). The HL test ¹³ provides a nonparametric omnibus test for this property. A nonsignificant HL statistic indicates that the forecast distribution is well specified.

2. (2)

   Logarithmic score and AG test

   The logarithmic score is expressed as

   \begin{equation} \mathit{LS}_{t} = -\ln f_{t\mid t-1}\left(y_t\right), \label{eq:6} \tag{6} \end{equation}

   where \(f\) denotes the forecast density. A lower score indicates better predictive accuracy. For the two competing models, we define

   \begin{equation} d_{t} = \mathit{LS}_{t}^{\,(\mathit{bench})} - \mathit{LS}_t^{\,(\mathit{ext})}. \label{eq:7} \tag{7} \end{equation}

The AG test ¹⁴ verifies whether \(H_{0}\colon E(d_{t})=0\) (equal predictive accuracy). A significantly positive mean of \(d_{t}\) (the AG test rejects the null hypothesis) indicates that the extended model outperforms the benchmark.

Causality criterion: The extended model is said to Granger-cause \(Y\) (in the distributional forecasting sense) if (i) its HL test is insignificant (the model is correctly specified) and (ii) the AG test is significant with \(\bar{d}_t>0\) (better predictive accuracy).

3.3.5. Remark on Multi-Step-Ahead Forecasting

Granger’s original definition is based on a one-step-ahead prediction ¹². Multistep forecasts from EGARCH models tend to accumulate nonlinear errors, which can distort inferences. Therefore, this study focused on one-step-ahead forecasts.

5.2.1. Testing Approach

Following the standard Granger causality logic, we compared two sets of models (Table 2). Model 1 (S\(\sim\)S) uses only spot lags and Model 2 (S\(\sim\)S+F) incorporates futures lags.

Similarly, Model 3 (F\(\sim\)F) uses only futures lags and Model 4 (F\(\sim\)F+S) incorporates spot lags.

5.2.2. Linear Causality Test

We estimate a VAR(4) model (optimal lag by the SC criterion) for the full sample. Linear Granger causality results (Table 3) show that at the 5% level, we cannot reject “Futures does not Granger-cause spot” (\({p=0.0516}\)); however, we clearly rejected “Spot does not Granger-cause futures” (\({p=0.0302}\)). Thus, from a linear-mean perspective, spots lead futures, whereas futures’ lead over spots is only marginally significant.

5.2.3. Nonlinear-Causality Test

1. (1)

   The BDS test on VAR(4) residuals strongly rejects the i.i.d. (\(p<0.01\)), thus confirming a nonlinear structure.

2. (2)

   The DP nonparametric test ⁷ (last two columns of Table 3) reveals significant bidirectional nonlinear Granger causality between spot and futures markets at the 5% level.

3. (3)

   The LASSO-based test ⁸ (Taylor expansion to quadratic terms, Table 4) shows that adding information pertaining to other markets substantially improves the model fit (RSS decreases, LM statistics significant at 1%). The LASSO selects many quadratic terms (e.g., \(S_{t-1}^2\)), thus confirming nonlinearity.⁽¹⁾

5.2.4. Performance Comparison: Linear vs. Nonlinear Modeling

To clarify whether nonlinear modeling is superior to linear modeling, we compared the in-sample fit of linear (estimated via OLS) and nonlinear models (estimated via LASSO, including quadratic terms). An F-test-based RSS was used. The results are summarized in Table 5.

The comparison results show that, except for the spot-only model, the nonlinear models significantly outperformed the linear models (\(F\)-test \({p<0.01}\)). This justifies the use of the nonlinear framework.

5.2.5. Time-Varying Analysis of Causality

The rolling-window tests (window size 1143, step 1) generated 1845 subsamples. Table 6 shows the number of significant causalities at the 5% level: both linear and nonlinear tests showed frequent futures-to-spot, spot-to-futures, and bidirectional leadership, thus confirming time variation.

Figures 1–3 detail the evolution. Before July 2015 (2015 crash), futures significantly led spot while spot did not lead futures, which is attributable to spot liquidity dry-ups and trading halts. After October 2020 (COVID-19 period), spot’s linear lead over futures strengthened, which is attributable to capital-market reforms improving spot pricing efficiency.

5.3.1. Rolling Forecast Setup

Before implementing the rolling-window distributional forecasting procedure, we performed a preliminary model selection on the full sample to determine the appropriate lag order, volatility specification, and error distribution.

1. ・

   Lag-order selection: A VAR model was fitted to the bivariate return series \((S_t,F_t)\). Based on the SC, an optimal lag order of four was selected. This indicates that including up to four lags in each market’s return provides the best in-sample fit.

2. ・

   Volatility and distribution selection: For each market, we estimated GARCH(1,1) and EGARCH(1,1) models under three distributional assumptions: normal, Student-\(t\), and skewed-\(t\). The models were compared using the AIC and BIC. EGARCH(1,1), which exhibited a skewed-\(t\) distribution, consistently yielded the lowest information criteria, thus confirming its ability to capture the leverage effects, fat tails, and skewness typical of financial returns.

Based on these preliminary results, we adopted the following rolling forecast setup:

1. ・

   Mean equation: The conditional mean was modeled as a function of lagged terms up to order four, including all first-order (linear) terms, second-order quadratic (squared) terms, and second-order cross-product terms (e.g., \(S_{t-i}\cdot S_{t-j}\), \(F_{t-i}\cdot F_{t-j}\), and \(S_{t-i}\cdot F_{t-j}\)). For the benchmark model (using only own-market information), the candidate set contained 14 variables (four linear lags, four squared lags, and six within-market cross-products). For the extended model (including other market’s information), the total number of candidate variables was 44. Variables were selected using the LASSO with a cross-validated penalty \(\lambda\).

2. ・

   Volatility equation: EGARCH(1,1) to capture volatility clustering and leverage effects.

3. ・

   Distributional assumption: Skewed-\(t\) distribution to accommodate fat tails and skewness.

4. ・

   Rolling forecast mechanism: The first estimation window was April 19, 2010 to December 31, 2014 (1143 trading days). Subsequently, the window was shifted forward one day at a time, thus generating 1845 one-step-ahead distribution forecasts for each model.

We compared the benchmark model (own-lags only) with the extended model (own lags \(+\) other market’s lags, including all linear, squared, and cross-product terms). If the extended model yields a statistically superior out-of-sample distributional forecast, then we conclude that the other market Granger-causes the target market distributionally.

5.3.2. Test Based on PIT

We applied the PIT to the 1845 one-step-ahead forecast distributions derived from four models: S\(\sim\)S, S\(\sim\)S+F, F\(\sim\)F, and F\(\sim\)F+S. Fig. 4 presents the histograms and autocorrelation functions of the resultant PIT sequences, which shows no clear visual evidence of severe deviation from uniformity or substantial autocorrelation. The HL test ¹³ serves as a rigorous diagnostic verification. Under the null hypothesis of the correct model specification, the test statistic \(W\) and moment statistic \(M(i,j)\) asymptotically reflect the standard normal distribution, with a 5% critical value of 1.645.

Table 7 presents the HL test results. For spot forecasting, the benchmark model S\(\sim\)S yielded a significant \(W(4.002)\), thus indicating misspecification. After incorporating futures information (S\(\sim\)S+F), \(W\) declined to 0.721 (insignificant), whereas the second-moment statistic \(M(2,2)\) improved from 2.531 to \(-0.419\). For futures forecasting, the benchmark F\(\sim\)F indicated a significant \(W(2.517)\); incorporating spot information (F\(\sim\)F+S) reduced \(W\) to 2.252, with clear improvement in \(M(1,1)\) (from 2.974 to 0.267). These results confirm that including another market’s information significantly improves the specification of the conditional return distribution, particularly its volatility component.

5.3.3. Test Based on Logarithmic Scoring Rule

As a complementary evaluation criterion, we employed the logarithmic scoring rule and the associated AG test ¹⁴. Based on the definitions in Section 3.3.4, the logarithmic score \(\mathit{LS}_t=\ln f_{t\mid t-1}(y_t)\) measures the accuracy of each density forecast, with lower values indicating better performance. For each pair of benchmark and extended models, we computed the score difference \(d_t=\mathit{LS}_t^{\,(\mathit{bench})}-\mathit{LS}_t^{\,(\mathit{ext})}\). The AG test examines whether the mean of \(d_t\) is significantly positive, which would imply that the extended model yields superior distributional forecasts.

Table 8 presents the average logarithmic scores over all rolling windows for the empirical period, in addition to the AG test \(p\)-values. The salient findings are as follows:

1. ・

   Spot market forecasting: The extended model S\(\sim\)S+F attains an average score of \(-1.625\), which is considerably higher (less negative) than \(-1.668\) of the benchmark S\(\sim\)S. The AG test \(p\)-value is \(1.660\times 10^{-16}\), which implies that the null of equal predictive accuracy should be rejected. Thus, adding futures information significantly improves the spot return distributional forecasts.

2. ・

   Futures market forecasting: Similarly, the extended model F\(\sim\)F+S yields an average score of \(-1.700\), thereby outperforming the benchmark F\(\sim\)F (\(-1.744\)). The AG test \(p\)-value is \(4.173\times 10^{-23}\), which is highly significant. Thus, spot information provides substantial incremental value for predicting the conditional distribution of futures returns.

These results are consistent with the HL test findings described in Section 5.3.2. Collectively, they provide robust evidence of bidirectional Granger causality in the distributional sense: historical information from one market not only aids one in specifying the correct conditional distribution of the other market (HL test) but also delivers economically meaningful improvements in out-of-sample predictive accuracy (AG test). The magnitude of the log-score gains (approximately 0.04 for both markets) is non-negligible, particularly in the context of financial risk management, where accurate density forecasts are critical for tail-risk assessment and portfolio optimization.

5.4.1. Variable Importance Measurement

To identify the lagged variables that contribute persistently to distributional forecasting, this study employed two complementary indicators:

1. (1)

   Selection frequency (\(\boldsymbol{\mathit{SI}}_{\boldsymbol{i}}\)): The proportion of rolling windows in which variable \(\textrm{i}\) is selected by the LASSO estimation, thus reflecting its stability.

   \begin{equation} \mathit{SI}_{i} = \dfrac{1}{T}\sum_{t=1}^{T} I\left(\hat{\beta}_{i,t}\ne 0\right), \label{eq:14} \tag{14} \end{equation}

   where \(T=1845\) denotes the total number of rolled windows. A higher \(\mathit{SI}_i\) value indicates that the variable is an effective predictor by the model across most market states.

2. (2)

   Conditional coefficient mean \(\boldsymbol{\vert\bar{\beta}\vert_{i}^{\textit{cond}}}\): The average absolute value of the standardized regression coefficient of variable \(i\) over the windows in which it is selected, thus reflecting its average marginal contribution strength when active.

   \begin{equation} \left\vert\bar{\beta}\right\vert_i^{\mathit{cond}} = \dfrac{1}{T_i}\sum_{\hat{\beta}_{{i,t}\ne 0}} \left\vert\hat{\beta}_{i,t}\right\vert, \label{eq:15} \tag{15} \end{equation}

   where \(T_i\) is the number of times variable \(i\) is selected. Because all variables are standardized before LASSO estimation, the coefficients of the different variables are comparable.

The two indicators capture the importance from different dimensions: a high \(\mathit{SI}_i\) indicates persistent and stable predictive information, and a high \(\vert\bar{\beta}\vert_i^{\mathit{cond}}\) indicates strong influence when selected. Combining both enables a more comprehensive identification of core transmission paths. Herein, variables with \(\mathit{SI}_i>0.8\) are defined as “important variables,” and their conditional coefficient means are reported.

5.4.2. Identification Results of Important Lagged Variables

Based on the method above, the important variables identified for the four forecasting models during the empirical period are listed in Table 9. All the selected variables have selection frequencies above 0.84, with most exceeding 0.93, thus indicating that these lagged terms provide persistently stable contributions to distributional forecasting.

1. (1)

   Spot-only model (S\(\boldsymbol{\sim}\)S): Two variables are important: the linear term \(S_{t-4}\) (\(\mathit{SI}=0.9306\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0671\)) and the nonlinear product term \(S_{t-2}\cdot S_{t-3}\) (\(\mathit{SI}=0.8455\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0306\)). The relatively modest coefficients suggest that spot returns exhibit mild linear autocorrelation at lag 4 and a weak but persistent nonlinear interaction between lags 2 and 3.

2. (2)

   Spot model with futures information (S\(\boldsymbol{\sim}\)S\(\boldsymbol{+}\)F): Futures lags \(F_{t-3}\) (\(\mathit{SI}=0.9447\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0798\)) and \(F_{t-4}\) (\(\mathit{SI}=0.9144\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0509\)) enter the important set, while the squared spot term \(S_{t-1}^2\) (\(\mathit{SI}=0.9783\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0512\)) remains highly relevant. This confirms that futures lead spots with a lag of three to four periods and that nonlinear volatility persistence in the spot market remains relevant even after futures information is controlled.

3. (3)

   Futures-only model (F\(\boldsymbol{\sim}\)F): The important variables are \(F_{t-3}\) (\(\mathit{SI}=0.9447\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0551\)), \(F_{t-4}\) (\(\mathit{SI}=0.9263\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0678\)) and \(F_{t-1}^2\) (\(\mathit{SI}=0.9534\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0145\)). The small coefficient of the quadratic term indicates that the futures market dynamics are predominantly linear.

4. (4)

   Futures model with spot information (F\(\boldsymbol{\sim}\)F\(\boldsymbol{+}\)S): Five variables are important: futures lags \(F_{t-3}\) (\(\mathit{SI}=0.9745\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.1535\)) and \(F_{t-4}\) (\(\mathit{SI}=0.9680\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0554\)); the squared spot term \(S_{t-1}^2\) (\(\mathit{SI}=0.9474\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0966\)); the cross-product \(S_{t-1}\cdot F_{t-4}\) (\(\mathit{SI}=0.9344\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.0721\)); and the nonlinear product \(S_{t-3}\cdot S_{t-4}\) (\(\mathit{SI}=0.9631\), \(\vert\bar{\beta}\vert^{\mathit{cond}}=0.1286\)). These results reveal that spot-market information exerts a strong and fast (lag-1) influence on futures returns through both linear (\(S_{t-1}\) implied by the cross-product and the large coefficient on \(F_{t-3}\) may capture indirect effects as well) and nonlinear channels. The high conditional coefficients for \(F_{t-3}\) (0.1535) and \(S_{t-3}\cdot S_{t-4}\) (0.1286) indicate that the spot-to-futures transmission is generally larger in magnitude than the reverse direction.

Collectively, the evidence supports the study’s two main conclusions: (i) A robust bidirectional nonlinear time-varying price-discovery process exists, and (ii) cross-market leadership is dominated by linear spillovers (e.g., \(F_{t-3}\) and \(F_{t-4}\)), while nonlinear effects (squared and product terms) are primarily internalized within each market’s own volatility structure.

5.4.3. Economic Implications and Strategy Insights

The identified transmission paths reveal a “dual-track asynchronous” price-discovery pattern: futures affect spot with a lag of three to four periods, whereas spots affect futures with a rapid one-period impact and sustained feedback at lags 3–4. Moreover, the spot-to-futures coefficients are generally larger (up to 0.1286) than the futures-to-spot coefficients (up to 0.0798), thus indicating that the spot market is the dominant pricing anchor.

Why do futures lead spot with a lag of three to four days? This seemingly counterintuitive finding can be explained by three interaction mechanisms. First, daily price limits in the Chinese A-share market (\(\pm 10\%\)) prevent spot prices from fully absorbing extreme information on the same day. When futures prices shift significantly (e.g., by –3%), the spot market may reach the limit, thus causing the price adjustment to continue over the following three to four trading days ¹⁸. Second, the retail-dominated investor structure (individual investors account for a large proportion of trading volume) results in slower information diffusion. Professional traders react to futures signals first, followed by retail investors, who require time for confirmation and herding, thus creating a multiday cascade ¹⁹. Third, arbitrage constraints. Short selling in the spot market is costly and operationally complex, thus limiting the ability of arbitrageurs to immediately eliminate mispricing and allowing the price discrepancy to persist for several days ²⁰. These institutional and behavioral factors simultaneously result in an observable 3–4 day lag, which is not a failure of market efficiency but rather a reflection of the market microstructure in an emerging market.

Practical implications. For dynamic hedging, portfolio managers should focus on futures returns and spot-market volatility at lags 3–4 and 1, respectively, thus expanding beyond conventional same-day or lag-1 models. The dominance of linear mean spillovers in cross-market transmission suggests that relatively simple linear models may be sufficient for constructing real-time cross-market risk monitoring systems. Finally, the persistent presence of the squared spot term \(S_{t-1}^2\) in both spot and futures models indicates that volatility forecasting should incorporate own-market nonlinear dynamics, even when other-market information is included.

5.5.1. Changing Sample Period (Excluding 2010–2013)

To assess whether the main conclusions are affected by the sample period, we repeated the analysis using a subsample that excluded the early market phase (2010–2013). The rationale is that the CSI 300 index futures were launched in April 2010. During the first few years, the market might have been affected by low liquidity, an immature investor base, and ongoing mechanism adjustments, which might have temporarily affected the price-discovery process. Focusing on the more mature period from 2014 onwards provides a clearer environment for examining the intrinsic futures–spot leadership relationship.

The HL test results for this subsample (Table 10) show that the benchmark models (S\(\sim\)S, F\(\sim\)F) remain misspecified (significant \(W\) statistics), whereas the extended models (S\(\sim\)S+F, F\(\sim\)F+S) yield insignificant \(W\) statistics, thereby indicating correctly specified forecast distributions. The AG test (Table 8, right column) confirms that the extended models significantly outperform the benchmarks (\({p<0.01}\)). Moreover, the important lagged variables identified in the robustness sample (Table 11) are nearly identical to those in the full sample, with \(F_{t-3}\), \(F_{t-4}\), \(S_{t-1}\), \(S_{t-1}^2\) remaining as the core drivers. Thus, the bidirectional distributional causality and specific transmission paths are not artifacts of the sample period.

5.5.2. Additional Remarks

Several other potential concerns are addressed through theoretical arguments instead of additional computations:

1. ・

   Multi-step-ahead forecasting: Granger causality was originally defined for one-step-ahead prediction ¹². Multistep forecasts from the EGARCH models accumulate nonlinear errors, which can distort inferences. Moreover, our analysis captures medium-term effects (3–4 day lags) through the lag structure in the mean equation (see Table 9). Therefore, we focus on the one-step-ahead forecasts, which are consistent with the theoretical definition.

2. ・

   Different lag orders (e.g., six or eight): The lag order of four was selected by the SC on the full sample. Nevertheless, the LASSO automatically excludes irrelevant longer lags if they do not contribute to predictive accuracy. Thus, increasing the maximum lag order would not alter the selection of the important variables identified in Table 9.

3. ・

   Alternative volatility models or error distributions: The preliminary model selection in Section 5.3.1 identified EGARCH with skewed-\(t\) as the optimal specification based on the AIC and BIC. Furthermore, the Monte Carlo simulations in Section 4 show that the proposed test maintains desirable finite-sample properties even under the more restrictive GARCH-normal setting. Consequently, deviations from the optimal specifications are unlikely to overturn our conclusions. Hence, we did not conduct additional robustness verifications beyond the sample-period change. The results reported in Sections 5.2–5.4 are thus reliable and robust to a range of reasonable alternatives.