% For submission to the Journal of Futures Markets (JFM)

\section{Econometric Methodology} \label{sec:methodolog}

We analyze volatility transmission between commodity ETFs and their iNAVs with two complementary frameworks. Heterogeneous Autoregressive (HAR) models \citep{corsi2009simple} capture the long memory of realized volatility across daily, weekly, and monthly horizons and, in cross-market form, yield interpretable spillover coefficients. Bayesian Vector Autoregression (BVAR) models \citep{koop2011forecasting, carriero2015forecasting} treat ETF and iNAV volatility as jointly endogenous, using Minnesota-prior shrinkage to control parameter proliferation and providing impulse responses and variance decompositions with proper uncertainty bands. We favor realized-variance models over GARCH because realized variance is a nearly model-free volatility estimate that forecasts better \citep{andersen2001distribution}, and HAR over fractionally integrated (ARFIMA) specifications because it captures long memory parsimoniously through a horizon cascade \citep{corsi2009simple}.

\subsection{Theoretical Framework for Volatility Transmission}

Under frictionless arbitrage, ETF prices equal their NAV and volatility transmission is instantaneous and bidirectional \citep{ackert2000arbitrage, petajisto2017inefficiencies}. Real-world frictions---transaction costs, inventory and funding constraints---can instead create asymmetric transmission that varies across horizons and market conditions. We study volatility transmission rather than price discovery because it reveals how risk propagates and whether arbitrage stabilizes or amplifies fluctuations, with direct implications for hedging effectiveness and systemic risk. Examining multiple horizons follows the heterogeneous market hypothesis \citep{muller1997volatilities}, under which day traders, institutions, and longer-horizon participants generate volatility that persists over their characteristic time scales. Finally, the contrast between physically-backed ETFs (GLD, SLV), whose arbitrage requires physical delivery, and futures-based ETFs (USO, UNG), which face roll costs and basis risk \citep{todorov2021etf}, motivates our cross-sectional comparison.

\subsection{Heterogeneous Autoregressive (HAR) Models}

The HAR model \citep{corsi2009simple} is the dominant framework for realized-volatility dynamics \citep{andersen2007roughing}. The baseline specification is

\begin{equation} \label{eq:har_baseline}
\log(RV_{i,t}) = \beta_0 + \beta_1 \log(RV_{i,t-1}) + \beta_2 \log(\overline{RV}_{i,t-5:t-1}) + \beta_3 \log(\overline{RV}_{i,t-22:t-1}) + \varepsilon_{i,t}
\end{equation}

where $\overline{RV}_{i,t-h:t-1} = \frac{1}{h}\sum_{j=1}^{h} RV_{i,t-j}$ averages realized variance over the previous $h$ days, capturing daily, weekly ($h=5$), and monthly ($h=22$) persistence. The log transform keeps fitted volatility positive and stabilizes the residual variance.

To measure transmission, we add cross-market terms. For iNAV volatility,

\begin{align} \label{eq:har_x_nav}
\log(RV_{t,\text{NAV}}) &= \beta_0 + \beta_1 \log(RV_{t-1,\text{NAV}}) + \beta_2 \log(\overline{RV}_{t-5:t-1,\text{NAV}}) + \beta_3 \log(\overline{RV}_{t-22:t-1,\text{NAV}}) \\
&\quad + \alpha_1 \log(RV_{t-1,\text{ETF}}) + \varepsilon_{t,\text{NAV}} \nonumber
\end{align}

and symmetrically for ETF volatility,

\begin{align} \label{eq:har_x_etf}
\log(RV_{t,\text{ETF}}) &= \gamma_0 + \gamma_1 \log(RV_{t-1,\text{ETF}}) + \gamma_2 \log(\overline{RV}_{t-5:t-1,\text{ETF}}) + \gamma_3 \log(\overline{RV}_{t-22:t-1,\text{ETF}}) \\
&\quad + \delta_1 \log(RV_{t-1,\text{NAV}}) + \varepsilon_{t,\text{ETF}} \nonumber
\end{align}

Each equation includes the own daily, weekly, and monthly terms and a single cross-market term at the daily lag, following the parsimony of the HAR cascade. The coefficient $\alpha_1$ measures the daily ETF-to-iNAV spillover and $\delta_1$ the daily iNAV-to-ETF spillover; the own weekly ($\beta_2,\gamma_2$) and monthly ($\beta_3,\gamma_3$) terms capture longer-horizon volatility persistence.

To separate the transmission of continuous and jump volatility, we estimate a HAR-CJ-X specification using the components $C$ and $J$ defined in Section~\ref{sec:data}:

\begin{align} \label{eq:har_cj_x_nav}
\log(RV_{t,\text{NAV}}) &= \beta_0 + \beta_1 \log(C_{t-1,\text{NAV}}) + \beta_2 \log(\overline{C}_{t-5:t-1,\text{NAV}}) + \beta_3 \log(\overline{C}_{t-22:t-1,\text{NAV}}) \\
&\quad + \beta_4 \log(1+J_{t-1,\text{NAV}}) + \beta_5 \log(1+\overline{J}_{t-5:t-1,\text{NAV}}) + \beta_6 \log(1+\overline{J}_{t-22:t-1,\text{NAV}}) \nonumber \\
&\quad + \alpha_1 \log(C_{t-1,\text{ETF}}) + \alpha_2 \log(1+J_{t-1,\text{ETF}}) + \varepsilon_{t,\text{NAV}} \nonumber
\end{align}

with an analogous ETF equation. The own continuous and jump components enter at all three horizons through $\beta_1$--$\beta_3$ and $\beta_4$--$\beta_6$; the daily cross-market continuous and jump spillovers are $\alpha_1$ and $\alpha_2$. The transform $\log(1+J_{t})$ keeps jump terms defined when $J_{t}=0$.

\subsection{Bayesian Vector Autoregression (BVAR)}

The HAR cascade imposes a fixed horizon structure. As a flexible complement, we estimate a VAR treating ETF and iNAV volatility as jointly endogenous:

\begin{equation} \label{eq:var_baseline}
\mathbf{y}_t = \mathbf{c} + \sum_{k=1}^{p} \mathbf{A}_k \mathbf{y}_{t-k} + \mathbf{u}_t
\end{equation}

where $\mathbf{y}_t = [\log(RV_{t,\text{ETF}}), \log(RV_{t,\text{NAV}})]'$ and $\mathbf{u}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{\Sigma})$. Information criteria (BIC) select $p=2$. Because unrestricted VARs over-parameterize, we apply the Minnesota prior \citep{litterman1986forecasting}:

\begin{align} \label{eq:minnesota_prior}
\beta_{ij}^{(k)} &\sim \mathcal{N}(0, \lambda_1^2 \cdot k^{-\lambda_3}) \quad \text{for } i \neq j \\
\beta_{ii}^{(1)} &\sim \mathcal{N}(1, \lambda_1^2) \nonumber \\
\beta_{ii}^{(k)} &\sim \mathcal{N}(0, \lambda_1^2 \cdot k^{-\lambda_3}) \quad \text{for } k > 1 \nonumber
\end{align}

where $\beta_{ij}^{(k)}$ is the coefficient on variable $j$ at lag $k$ in equation $i$; $\lambda_1$ controls overall tightness, $\lambda_2$ cross-variable shrinkage, and $\lambda_3$ lag decay. We set $\lambda_1=0.2$, $\lambda_2=0.5$, $\lambda_3=2$ and vary them in sensitivity analysis. The prior encodes that own lags matter more than others, recent lags more than distant ones, and coefficients are not extreme---assumptions well suited to persistent, mean-reverting volatility.

We estimate by Gibbs sampling, alternating between the coefficients (multivariate normal given $\mathbf{\Sigma}$) and the covariance (inverse-Wishart given the coefficients):

\begin{equation} \label{eq:coefficient_posterior}
\boldsymbol{\beta} | \mathbf{\Sigma}, \mathbf{Y} \sim \mathcal{N}(\hat{\boldsymbol{\beta}}, \mathbf{\Sigma} \otimes (\mathbf{X}'\mathbf{X} + \mathbf{V}_0^{-1})^{-1}), \qquad
\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X} + \mathbf{V}_0^{-1})^{-1}(\mathbf{X}'\text{vec}(\mathbf{Y}) + \mathbf{V}_0^{-1}\boldsymbol{\beta}_0)
\end{equation}

\begin{equation} \label{eq:covariance_posterior}
\mathbf{\Sigma} | \boldsymbol{\beta}, \mathbf{Y} \sim \text{IW}(\mathbf{S} + \mathbf{S}_0, T + \nu_0)
\end{equation}

where $\mathbf{Y}$ and $\mathbf{X}$ are the stacked data, $\mathbf{V}_0$, $\mathbf{S}_0$, and $\nu_0$ are prior parameters, $\mathbf{S}$ the residual sum of squares, and $T$ the sample size. We run 50{,}000 iterations, discard 10{,}000 as burn-in, and thin every tenth draw; convergence is assessed with trace plots and the Geweke test.

From the posterior we compute orthogonalized impulse responses and forecast error variance decompositions:

\begin{equation} \label{eq:impulse_response}
\text{IRF}(h) = \mathbf{C}_h \mathbf{P}, \qquad
\text{FEVD}_{i,j}(h) = \frac{\sum_{k=0}^{h-1} [\mathbf{C}_k \mathbf{P}]_{i,j}^2}{\sum_{k=0}^{h-1} [\mathbf{C}_k \mathbf{\Sigma} \mathbf{C}_k']_{i,i}}
\end{equation}

where $\mathbf{C}_h$ is the $h$-step moving-average matrix and $\mathbf{P}$ the Cholesky factor of $\mathbf{\Sigma}$. We order iNAV before ETF, so iNAV innovations may affect ETF volatility contemporaneously but not the reverse, reflecting that underlying price movements lead ETF adjustments through arbitrage. The confidence bands follow \citet{sims1999error}. A high $\text{FEVD}_{\text{ETF},\text{NAV}}(h)$ indicates that iNAV innovations explain ETF volatility, and conversely for $\text{FEVD}_{\text{NAV},\text{ETF}}(h)$.

\subsection{Hypotheses and Testing}

We test four hypotheses, mapping each to coefficients in the equations above.

\textbf{H1 (Direction and cross-commodity heterogeneity).} Theory predicts unidirectional iNAV-to-ETF transmission, but a more liquid ETF can reverse the flow \citep{glosten2021etf}. We expect physically-backed precious metals to show unidirectional iNAV-to-ETF transmission and futures-based energy ETFs more balanced bidirectional transmission. In equations~\eqref{eq:har_x_nav}--\eqref{eq:har_x_etf}, ETF-to-iNAV transmission tests $H_0:\alpha_1=0$ against $H_a:\alpha_1\neq 0$, and iNAV-to-ETF transmission tests $H_0:\delta_1=0$ against $H_a:\delta_1\neq 0$. Transmission is unidirectional when only $\delta_1$ is significant and bidirectional when both $\alpha_1$ and $\delta_1$ are significant. Heterogeneity means these outcomes differ across commodities.

\textbf{H2 (Frequency dependence).} If arbitrage occurs within minutes, daily aggregation would understate transmission. Using the daily spillover coefficients $\delta_1$ and $\alpha_1$ estimated at the 1-, 5-, and 30-minute frequencies, we test $H_0: \delta_1^{(1\text{m})}=\delta_1^{(5\text{m})}=\delta_1^{(30\text{m})}$ against $H_a:$ the coefficients differ across frequencies.

\textbf{H3 (Jumps versus continuous component).} Tail risk from discrete shocks is harder to hedge than smooth volatility, which is why we decompose volatility into continuous and jump parts. In the HAR-CJ-X model~\eqref{eq:har_cj_x_nav}, we compare the daily cross-market jump and continuous spillovers, testing $H_0:\alpha_2=\alpha_1$ (jump and continuous transmission equal) against $H_a:\alpha_2>\alpha_1$ (jump transmission dominates), and analogously for the ETF equation.

\textbf{H4 (Stability over time).} Over a data sample spanning several crises, patterns of transmission may evolve. We re-estimate the BVAR over sub-periods (2010--2014, 2015--2019, 2020--2023) and apply the structural-break test of \citet{bai2003computation}, testing $H_0:$ transmission coefficients are constant across regimes against $H_a:$ they change.

\subsection{Estimation, Model Selection, and Robustness}

We select lag lengths and compare nested models with the AIC and BIC (the DIC for Bayesian models, which accounts for shrinkage), and evaluate out-of-sample forecasts with rolling windows (a 1{,}000-day estimation window re-estimated every 250 days) at the 1-, 5-, and 22-day horizons using RMSE and MAE. Residual diagnostics include the Ljung--Box (serial correlation), Breusch--Pagan (heteroskedasticity), and Jarque--Bera (normality) tests, plus posterior predictive checks for the BVAR. We also gauge economic magnitude through the effect of one-standard-deviation shocks on forecast volatility.

We assess robustness along three dimensions, consolidating the checks noted in Section~\ref{sec:data}. First, \emph{specifications}: alternative HAR lag structures, BVAR lags from one to four, level and square-root (rather than log) transforms, and time-varying-parameter versions. Second, \emph{sample stability}: the sub-period and regime estimation underlying H4, structural-break tests \citep{bai2003computation}---which flag the 2014 oil collapse, the 2016 Brexit referendum, and the 2020 pandemic---and bootstrap inference. Our transmission results hold across regimes, indicating they capture general mechanisms rather than period-specific effects. Third, \emph{data construction}: additional 15- and 60-minute frequencies, alternative estimators (truncated realized variance, realized kernels, range-based measures), and alternative iNAV constructions (currency conversion, cash treatment), including validation against the GLD IOPV. Our main findings are robust throughout.
