% ============================================================================
\section{Methodology}
\label{sec:methodology}
% ============================================================================

\subsection{Realized Volatility}
\label{sec:rv}

For instrument $i$ on calendar day $t$, let $P_{i,t,j}$ denote the closing price of the $j$-th 5-minute bar, $j = 1, \ldots, M_{i,t}$, where $M_{i,t}$ is the number of bars available on that day. The 5-minute log return and the daily realized variance are
\begin{equation}
    r_{i,t,j} = \ln P_{i,t,j} - \ln P_{i,t,j-1},
    \qquad
    \RV_{i,t} = \sum_{j=1}^{M_{i,t}} r_{i,t,j}^2 .
    \label{eq:rv}
\end{equation}
If the log price is a continuous semimartingale with spot variance $\sigma_{i,t}^2(s)$ and the sampling interval shrinks, $\RV_{i,t}$ converges in probability to the integrated variance $\int_{t-1}^{t} \sigma_{i,t}^2(s)\,ds$ \citep{Andersen2003,BarndorffNielsen2002}; in the presence of jumps it converges to the integrated variance plus the sum of squared jumps \citep{Andersen2007}. At the 5-minute frequency the bias induced by microstructure noise is small for liquid instruments \citep{HansenLunde2006,Liu2015}. We report realized volatility in annualized form,
\begin{equation}
    \sigma_{i,t} = \sqrt{252 \cdot \RV_{i,t}},
    \label{eq:rv_annualized}
\end{equation}
and require $M_{i,t} \geq 20$ bars for a day to be retained. Because $\sigma_{i,t}$ is strongly right-skewed and heavy-tailed whereas its logarithm is close to Gaussian \citep{Andersen2001,Andersen2003}, the connectedness analysis is conducted on
\begin{equation}
    y_{i,t} = \ln \sigma_{i,t},
    \label{eq:logrv}
\end{equation}
which is equal to $\tfrac12 \ln \RV_{i,t}$ up to a constant. The transformation matters for the VAR because the generalized variance decomposition weights observations by their squared innovations: in levels, a handful of extreme days dominates the residual covariance matrix, and the decomposition then describes those days rather than the typical transmission mechanism. Results in levels are reported as a robustness check. Stationarity of $y_{i,t}$ is assessed with the augmented Dickey--Fuller test \citep{DickeyFuller1979} and the KPSS test \citep{Kwiatkowski1992}.

\subsection{Vector Autoregression}
\label{sec:var}

Let $\mathbf{y}_t = (y_{1,t}, \ldots, y_{N,t})'$ collect the $N = 10$ log realized volatilities on trading day $t$. The dynamics are approximated by a covariance-stationary VAR($p$),
\begin{equation}
    \mathbf{y}_t = \boldsymbol{\mu} + \sum_{\ell=1}^{p} \boldsymbol{\Phi}_\ell \mathbf{y}_{t-\ell} + \boldsymbol{\varepsilon}_t,
    \qquad \mathbb{E}[\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_t'] = \boldsymbol{\Sigma},
    \label{eq:var}
\end{equation}
estimated equation by equation by least squares \citep{Lutkepohl2005}. The lag order is chosen by comparing the Akaike, Schwarz and Hannan--Quinn criteria over $p \in \{1, \ldots, 10\}$ (\Cref{tab:lag_selection} in the appendix). The criteria select $p = 5$, $3$ and $4$ respectively; we adopt the Hannan--Quinn choice, $p = 4$, which coincides with the original specification of the paper, and report $p \in \{1, 2, 3, 6\}$ as robustness checks. The moving-average representation of \eqref{eq:var} is $\mathbf{y}_t = \boldsymbol{\mu}^{*} + \sum_{h=0}^{\infty} \boldsymbol{\Psi}_h \boldsymbol{\varepsilon}_{t-h}$ with
\begin{equation}
    \boldsymbol{\Psi}_0 = \mathbf{I}_N,
    \qquad
    \boldsymbol{\Psi}_h = \sum_{\ell=1}^{\min(h,p)} \boldsymbol{\Phi}_\ell \boldsymbol{\Psi}_{h-\ell}, \quad h \geq 1 .
    \label{eq:psi_recursion}
\end{equation}

\subsection{Generalized Forecast-Error Variance Decomposition}
\label{sec:gfevd}

Following \citet{DieboldYilmaz2012}, the $H$-step-ahead generalized forecast-error variance decomposition of \citet{Koop1996} and \citet{Pesaran1998} is
\begin{equation}
    \theta_{ij}(H) = \frac{\sigma_{jj}^{-1} \sum_{h=0}^{H-1} \left(\mathbf{e}_i' \boldsymbol{\Psi}_h \boldsymbol{\Sigma} \mathbf{e}_j\right)^2}{\sum_{h=0}^{H-1} \mathbf{e}_i' \boldsymbol{\Psi}_h \boldsymbol{\Sigma} \boldsymbol{\Psi}_h' \mathbf{e}_i},
    \label{eq:gfevd}
\end{equation}
where $\sigma_{jj}$ is the $j$-th diagonal element of $\boldsymbol{\Sigma}$ and $\mathbf{e}_i$ is the $i$-th unit vector. The numerator is the contribution of a one-standard-deviation shock to variable $j$---allowing for its contemporaneous correlation with the other shocks---to the $H$-step forecast-error variance of variable $i$, and the denominator is that forecast-error variance. Because the generalized shocks are not orthogonal, the rows of $[\theta_{ij}(H)]$ do not sum to one, and we normalize
\begin{equation}
    \tilde{\theta}_{ij}(H) = \frac{\theta_{ij}(H)}{\sum_{k=1}^{N} \theta_{ik}(H)},
    \qquad \sum_{j=1}^{N} \tilde{\theta}_{ij}(H) = 1,
    \quad \sum_{i,j=1}^{N} \tilde{\theta}_{ij}(H) = N .
    \label{eq:gfevd_normalized}
\end{equation}
We use $H = 10$ trading days, as in \citet{DieboldYilmaz2012}, and report $H \in \{5, 20\}$ as robustness checks. Since the decomposition converges as $H$ grows, the choice of horizon affects the level of the measures more than the ranking of the instruments.

\subsection{Connectedness Measures}
\label{sec:measures}

All measures are expressed in percent of forecast-error variance. The total spillover index is the share of forecast-error variance in the system that is due to cross-variable shocks,
\begin{equation}
    \TSI(H) = \frac{100}{N} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \tilde{\theta}_{ij}(H) .
    \label{eq:tsi}
\end{equation}
The directional spillover received by $i$ from all other variables and the spillover transmitted by $j$ to all other variables are the off-diagonal row and column sums,
\begin{equation}
    \mathrm{FROM}_i(H) = 100 \sum_{\substack{j=1 \\ j \neq i}}^{N} \tilde{\theta}_{ij}(H),
    \qquad
    \mathrm{TO}_j(H) = 100 \sum_{\substack{i=1 \\ i \neq j}}^{N} \tilde{\theta}_{ij}(H),
    \label{eq:directional}
\end{equation}
and the net directional spillover of $i$ is
\begin{equation}
    \mathrm{NET}_i(H) = \mathrm{TO}_i(H) - \mathrm{FROM}_i(H) .
    \label{eq:net_spillover}
\end{equation}
$\mathrm{FROM}_i$ lies between 0 and 100; $\mathrm{TO}_j$ is unbounded above and can exceed 100 when a variable explains a large share of many others; and $\mathrm{NET}_i$ is positive for a net transmitter and negative for a net receiver, with $\sum_i \mathrm{NET}_i = 0$. We express $\mathrm{NET}_i$ in percentage points. Finally, the net pairwise directional spillover from $i$ to $j$,
\begin{equation}
    \mathrm{NET}_{i \to j}(H) = 100\,\big[\tilde{\theta}_{ji}(H) - \tilde{\theta}_{ij}(H)\big],
    \label{eq:net_pairwise}
\end{equation}
is the signed edge of the connectedness network of \citet{DieboldYilmaz2014}.

\subsection{Inference}
\label{sec:inference}

The connectedness measures are smooth functions of $(\boldsymbol{\Phi}_1, \ldots, \boldsymbol{\Phi}_p, \boldsymbol{\Sigma})$, and we approximate their sampling distribution with a residual bootstrap \citep{Kilian1998,KilianLutkepohl2017}. Let $\hat{\boldsymbol{\varepsilon}}_t$ denote the centred least-squares residuals of \eqref{eq:var}. For each of $B = 500$ replications, residuals are resampled with replacement, an artificial sample is generated recursively from the estimated coefficients and the first $p$ observed values, the VAR is re-estimated and the measures are recomputed. We report 90\% percentile intervals. The bootstrap treats the residuals as independent and identically distributed and therefore does not capture conditional heteroskedasticity in the innovations of log realized volatility; the intervals should be read as indicative of estimation uncertainty rather than as exact. For the calendar-year sub-samples we use $B = 300$.

\subsection{Rolling Windows and Sub-Samples}
\label{sec:rolling}

Time variation is examined in three complementary ways. First, the VAR and the decomposition are re-estimated over rolling windows of $w$ trading days advancing by one day,
\begin{equation}
    \TSI_t(H; w) = \TSI\big(H; \mathbf{y}_{t-w+1}, \ldots, \mathbf{y}_t\big),
    \label{eq:rolling_tsi}
\end{equation}
and similarly for $\mathrm{NET}_{i,t}(H; w)$. A ten-variable VAR(4) with intercepts has 41 coefficients per equation, so that a 200-day window leaves fewer than five observations per coefficient; the resulting over-fitting biases the total index upwards and makes it noisy. We therefore report $w \in \{200, 250, 500\}$ and treat the 500-day window as the more reliable description of slow-moving trends. Second, four regime sub-samples are estimated separately with bootstrap intervals: the pre-COVID period (January 2019--February 2020), the COVID crash (March--June 2020), the 2022 rate-hiking cycle (January--December 2022) and the spot-ETF era (11 January 2024--March 2026). For the 85-day COVID window the lag order is reduced to $p = 1$; the VAR(4) estimates on the same window are reported in the appendix for comparison. Third, calendar-year estimates with $p = 2$ track the integration of the cryptocurrencies through the shares of bitcoin's forecast-error variance explained by the traditional instruments and the net pairwise spillovers between bitcoin and ES.

\subsection{Granger Causality}
\label{sec:granger}

For a pair of log realized volatilities $(y_{i,t}, y_{j,t})$ we test whether $p$ lags of $y_{j}$ help predict $y_{i}$ beyond $p$ own lags \citep{Granger1969}, with the regression
\begin{equation}
    y_{i,t} = \alpha + \sum_{\ell=1}^{p} \beta_\ell\, y_{i,t-\ell} + \sum_{\ell=1}^{p} \gamma_\ell\, y_{j,t-\ell} + u_t,
    \qquad H_0 : \gamma_1 = \cdots = \gamma_p = 0 .
    \label{eq:granger}
\end{equation}
We report the conventional $F$-statistic, $F = [(\mathrm{RSS}_R - \mathrm{RSS}_U)/p] / [\mathrm{RSS}_U/(T - 2p - 1)]$, and a Wald statistic based on a heteroskedasticity- and autocorrelation-consistent covariance matrix with ten lags \citep{NeweyWest1987}, which is more appropriate for volatility series. Granger causality is a statement about incremental one-day-ahead predictability, not about the economic size of the effect; we therefore also report the sum $\sum_\ell \hat{\gamma}_\ell$, which measures the long-run response of $y_i$ to a unit change in $y_j$ in this bivariate regression.

\subsection{Intraday Periodicity and Session-Ordered Predictability}
\label{sec:sessions_method}

The intraday analysis uses the 5-minute returns directly. Exchange time stamps are converted to Coordinated Universal Time (UTC) with the civil daylight-saving rules of the exchange's time zone, and three sessions are defined in UTC: Asian (00:00--08:00), European (08:00--14:30) and U.S. (14:30--21:00). For each instrument the hourly profile of volatility is the mean absolute 5-minute return by hour of the day, and the session realized variance on day $t$ is
\begin{equation}
    \RV^{s}_{i,t} = \sum_{j \in s} r_{i,t,j}^2, \qquad s \in \{\text{Asian}, \text{European}, \text{U.S.}\},
    \label{eq:session_rv}
\end{equation}
retained when the session contains at least twelve bars. The share of daily variance generated in session $s$ is the average over days of $\RV^{s}_{i,t} / \sum_{s'} \RV^{s'}_{i,t}$.

The chronological ordering of sessions within a UTC day makes it possible to test whether volatility realized in an earlier session helps predict volatility realized in a later one. We estimate
\begin{equation}
    \ln \RV^{\mathrm{US}}_{i,t} = \alpha + \sum_{\ell=1}^{4} \beta_\ell \ln \RV^{\mathrm{US}}_{i,t-\ell} + \sum_{\ell=0}^{3} \gamma_\ell \ln \RV^{s}_{j,t-\ell} + u_t,
    \label{eq:session_predictive}
\end{equation}
where $i$ is ES, NQ or the VIX, $j$ is bitcoin and $s$ is the Asian or the European session of the same day, and test $\gamma_0 = \cdots = \gamma_3 = 0$ with a HAC Wald statistic. The reverse direction replaces the dependent variable by bitcoin's Asian-session variance on day $t$ and the predictor by the U.S.-session variance of ES, NQ or the VIX on the previous trading day. A final regression asks whether bitcoin's weekend variance (Saturday plus Sunday) predicts the Monday U.S.-session variance of ES once Friday's ES and bitcoin variances are controlled for. These regressions complement the daily VAR by exploiting an ordering of information that daily data cannot provide.

\subsection{Network Representation}
\label{sec:network_method}

The full-sample decomposition is represented as a directed graph whose nodes are the ten instruments and whose edges are the net pairwise directional spillovers of \eqref{eq:net_pairwise}, drawn from the net transmitter to the net receiver when $|\mathrm{NET}_{i \to j}|$ exceeds the 60th percentile of all pairwise values. Node size is proportional to $\mathrm{TO}_j$ and nodes are coloured by asset class.