A Time-Varying Gerber Statistic: Application of a Novel Correlation Metric to Commodity Price Co-Movements ^†

Abstract

This study investigates the daily co-movements in commodity prices over the period 2006–2020 using a novel approach based on a time-varying Gerber correlation. The statistic is computed considering a set of probabilities estimated via non-traditional models that give a time-varying structure to the measure. The results indicate that there are several co-movements across commodities, that these co-movements change over time, and that they are tendentially positive. Conditional auto-regressive multithreshold logit models show higher forecasting accuracy for agricultural returns, while dynamic conditional correlation models are more accurate for energy products and metals. The proposed models are shown to be superior in terms of forecasting power to the benchmark method which is based on estimating the Gerber correlation moving a rolling window.

1. Introduction

The new millennium witnessed a rapid growth in commodity trading and a progressive process of financialization of commodity markets. This process was fueled by the diversification opportunities and the ‘safe haven’ characteristics commodities offer [1,2,3,4,5]. It was further accelerated by the financial deregulation measures adopted in the US in the early 2000s. Indeed, the Commodity Futures Modernization Act [6] increased the presence and importance of financial investors in commodity markets. Since then, several agricultural, metal and energy commodities have experienced large sequences of price swings, high volatility and an unprecedented surge in the dependence across returns on different commodities [7,8]. Commodity price dynamics and co-movements (see [9] for a discussion on the definition of co-movement) have several consequences for the economy. Sharp price fluctuations, for instance, could reduce the ability of consumers to secure supplies and increase the risks of producers in terms of low returns [10]. Price hikes and volatility could discourage investments and reduce economic and political stability [11,12,13]. Movements in commodity prices could affect a country’s external and internal balance (i.e., the current account and employment/unemployment rates), its fiscal and monetary policies as well as the business cycle [14,15]. The tendency of commodity prices to co-move has important welfare implications for both investors and commodity importers/exporters. For example, investors could have problems in managing commodity portfolios given the easy transmission of shocks from one commodity market to another. Likewise, a synchronized increase in commodity prices is likely to place commodity import-dependent countries under considerable inflation pressure.

The widespread and undesirable consequences associated with commodity price fluctuations and their joint wavering underline the critical importance of investigating the co-movements in commodity prices. The present study goes in this direction and enriches the literature by introducing a novel method to assess co-movements across commodities and examine their common price dynamics. Specifically, 14 commodities encompassing all major metals, energy and agricultural commodities are examined on a daily basis for the period from 2006 to 2020. To assess the pervasiveness of co-movements in commodity prices, we propose a set of models that allow the computation of a time-varying Gerber correlation statistic [16]. The latter provides a robust measurement of common price movements and is well suited to quantify the degree of dependency in financial markets. The Gerber correlation statistic was recently applied to commodity spot markets by Zaremba et al. [17] who examined monthly price changes over a period of 170 years. The authors estimated the Gerber statistic by moving a rolling window (see Section 2.1.1) and found that cross-commodity correlations are not so ‘unprecedented’ as usually expected, thereby casting doubt on the link with the phenomenon of financialisation. Differently from Zaremba et al. [17], we explicitly examine commodity futures markets and consider daily instead of monthly data. In addition, we propose a novel approach based on a set of models that enable us to compute a time-varying Gerber correlation in a more robust way. Specifically, the conditional auto-regressive multithreshold logit (CARML), the dynamic conditional correlation (DCC) and the filtered historical simulation (FHS) models are considered. Rather than using traditional methodologies based on time-varying parameters, such as state space models (see, e.g., [18]), we apply dynamic models that permit us to model the conditional variance matrix (driving the dependence between commodity returns) or the logit transformation of the probability appearing in the Gerber statistic. The proposed models turn out to be superior to the benchmark method, i.e., historical simulation. Therefore, these models are better suited to gauge the extent to which the degree of dependency in financial markets changes over time.

The remainder of the study is organized as follows. Section 2 describes the main feature of the Gerber statistic and the models we use to make the Gerber correlation vary over time. Section 3 shows the empirical results for the considered commodities and the robustness checks. Section 4 concludes.

2. Materials and Methods

The Gerber correlation statistic [16] is a measure well-suited to capture the characteristics of financial time series, namely volatility clustering, leptokurtosis and outliers, that could distort the classic Pearson’s coefficient. Furthermore, it is insensitive to extremely large co-movements that distort product-moment-based measures, while also being insensitive to small noise fluctuations. The Gerber statistic, g, is, hence, a robust measure of pairwise movements defined for a sample of T observations as

$$g_{ij}=\frac{{\sum }_{t=1}^T{m}_{ij}\left(t\right)}{{\sum }_{t=1}^T\left|m_{ij}\left(t\right)\right|}$$

(1)

where

$$\begin{array}{c}{m}_{ij}\left(t\right)=I\left(r_{i,t}\ge Q_i\right)I\left(r_{j,t}\ge Q_j\right)+I\left(r_{i,t}\le -Q_i\right)I\left(r_{j,t}\le -Q_j\right)\\ -I\left(r_{i,t}\ge Q_i\right)I\left(r_{j,t}\le -Q_j\right)-I\left(r_{i,t}\le -Q_i\right)I\left(r_{j,t}\ge Q_j\right).\end{array}$$

Here, $r_{i,t}$ denotes the return (change in log prices) of asset i at time t, $Q_i$ is a threshold and $I\left(A\right)$ denotes the indicator function for the event A. $I\left(A\right)$ equals 1 if A is true and 0 otherwise. Hence, $m_{ij}\left(t\right)$ equals 1 if both returns simultaneously pierce their thresholds in the same direction at time t; it equals $-1$ if they pierce their thresholds in opposite directions at time t; or it is zero if at least one of the returns is in absolute value smaller that its threshold.

For simplicity, from now on, we drop the two indices $(i,j)$. Dividing both the numerator and the denominator of (1) by T, we obtain an alternative representation of the Gerber statistic:

$$g=\frac{p_{UU}+p_{DD}-p_{UD}-p_{DU}}{p_{UU}+p_{DD}+p_{UD}+p_{DU}}$$

(2)

where $p_{UU}$ denotes the empirical probability that both returns lie above their upper threshold, $p_{UD}$ the probability that the first return is above its upper threshold and the second one is below its lower threshold (i.e., $-Q_j$), and so on. Note that Equation (2) does not include the probability that at least one return lies between the lower and the upper threshold. As in [16], we use the subscript N (neutral) for this case. Therefore, all nine possible probabilities are identified by $p_{ab}$ with $(a,b)\in \mathcal{J}$, where the set $\mathcal{J}=\{U,N,D\}\times \{U,N,D\}$.

In a similar way, it is possible to derive a different type of Gerber correlation:

$$g^{\left(2\right)}=\frac{p_{UU}+p_{DD}-p_{UD}-p_{DU}}{\sqrt{p^{\left(A\right)}{p}^{\left(B\right)}}}$$

(3)

where

$$\begin{array}{c}{p}^{\left(A\right)}=p_{UU}+p_{UN}+p_{UD}+p_{DU}+p_{DN}+p_{DD}\\ p^{\left(B\right)}=p_{UU}+p_{NU}+p_{UD}+p_{DU}+p_{ND}+p_{DD}.\end{array}$$

The correlation matrix constructed from (3) is positive semidefinite [16].

2.1. Dynamic Models

With the aim of constructing a Gerber correlation statistic that changes over time and is able to capture daily common movements between return pairs of commodities, we propose first to use a set of models to derive the probabilities appearing in (2) and/or (3) and then employ them for the calculation of the statistic itself. This methodology enables us to measure the time-varying co-movements for all the combinations of commodity pairs.

The following sections briefly describe the models considered for the estimations of the aforementioned probabilities.

2.1.1. Historical Simulation

Historical simulation (HS) is the non-parametric traditional method that is used as a benchmark for the comparisons with other models presented in Section 2.1.2, Section 2.1.3 and Section 2.1.4. Based on a window of T observations, the HS method provides the $T+1$-forecast for the probabilities associated to the nine regions identified by the set $\mathcal{J}$. For instance, the $T+1$-forecasts for the probabilities associated to the regions $UD$ and $NU$ are given by:

$$\begin{array}{c}{\widehat{p}}_{UD,T+1|T}=\frac{1}{T}\sum _{t=1}^{T}I\left(r_{i,t}\ge Q_i\right)I\left(r_{j,t}\le -Q_j\right)\\ {\widehat{p}}_{NU,T+1|T}=\frac{1}{T}\sum _{t=1}^{T}I\left(-Q_i<r_{i,t}<Q_i\right)I\left(r_{j,t}\ge Q_j\right).\end{array}$$

2.1.2. Conditional Auto-Regressive Multithreshold Logit Models

As a first alternative to HS, the CARML model proposed by Taylor [19] is considered. The model allows us to estimate the probabilities needed for the dynamic Gerber correlation statistic. To start with, as in a multinomial logit model, we express the probabilities as

$$\begin{array}{c}{p}_{ab,t}=\frac{{\mathrm{e}}^{x_{ab,t}}}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,t}}},(a,b)\in \mathcal{J}\backslash (N,N);\\ p_{NN,t}=\frac{1}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,t}}},\end{array}$$

where the state $(N,N)$ is assumed to be the reference category.

Next, we consider a linear model for each of the eight logit transforms. For instance,

$$x_{ab,t}={\gamma }_{ab}+{\alpha }_{ab}{f}_{ab}(r_{i,t-1},r_{j,t-1})+{\beta }_{ab}{x}_{ab,t-1},(a,b)\in \mathcal{J}\backslash (N,N).$$

(4)

A possible choice of the functions $f_{ab}(·)$ could be as follows:

$$\begin{array}{c}{f}_{UU}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\ge Q_i\right)I\left(r_{j,t-1}\ge Q_j\right)\\ f_{DD}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\le -Q_i\right)I\left(r_{j,t-1}\le -Q_j\right)\\ f_{UD}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\ge Q_i\right)I\left(r_{j,t-1}\le -Q_j\right)\\ f_{DU}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\le -Q_i\right)I\left(r_{j,t-1}\ge Q_j\right)\\ f_{NU}(r_{i,t-1},r_{j,t-1})=I\left(-Q_i<r_{i,t-1}<Q_i\right)I\left(r_{j,t-1}\ge Q_j\right)\\ f_{UN}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\ge Q_i\right)I\left(-Q_j<r_{j,t-1}<Q_j\right)\\ f_{ND}(r_{i,t-1},r_{j,t-1})=I\left(-Q_i<r_{i,t-1}<Q_i\right)I\left(r_{j,t-1}\le -Q_j\right)\\ f_{DN}(r_{i,t-1},r_{j,t-1})=I\left(r_{i,t-1}\le -Q_i\right)I\left(-Q_j<r_{j,t-1}<Q_j\right)\end{array}$$

The model parameters can be estimated by maximizing the likelihood derived from a categorical distribution:

$$L\left(\theta \right)=\prod _{t=1}^T\prod _{(a,b)\in \mathcal{J}}{p}_{ab,t}^{I_{ab,t}}$$

where $I_{ab,t}$ is 1 if returns at time t are in the region identified by the indices $(a,b)$ and 0 otherwise. A model based on (4) comprises 24 (=$3\times 8$) parameters. To reduce the number of parameters to three, we consider the following dynamics:

$$x_{ab,t}=\gamma +\alpha f_{ab}(r_{i,t-1},r_{j,t-1})+\beta x_{ab,t-1},(a,b)\in \mathcal{J}\backslash (N,N).$$

(5)

Once model (4) has been estimated using the return data ${\left\{(r_{i,t},r_{j,t})\right\}}_{t=1,\dots ,T}$, we obtain the time-varying probabilities as

$$\begin{array}{cc}\hfill {\widehat{p}}_{ab,t|T}=& \frac{{\mathrm{e}}^{x_{ab,t|T}}}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,t|T}}},(a,b)\in \mathcal{J}\backslash (N,N);\hfill \\ \hfill {\widehat{p}}_{NN,t|T}=& \frac{1}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,t|T}}},\hfill \end{array}$$

(6)

with

$$x_{ab,t|T}={\gamma }_{ab}+{\alpha }_{ab}{f}_{ab}(r_{i,t},r_{j,t})+{\beta }_{ab}{x}_{ab,t-1|T},(a,b)\in \mathcal{J}\backslash (N,N)\mathrm{and}t=1,\dots ,T$$

(7)

or its reduced form

$$x_{ab,t|T}=\gamma +\alpha f_{ab}(r_{i,t},r_{j,t})+\beta x_{ab,t-1|T},(a,b)\in \mathcal{J}\backslash (N,N)\mathrm{and}t=1,\dots ,T$$

(8)

Furthermore, we forecast $p_{ab,T+1}$ as

$$\begin{array}{c}{\widehat{p}}_{ab,T+1|T}=\frac{{\mathrm{e}}^{x_{ab,T+1|T}}}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,T+1|T}}},(a,b)\in \mathcal{J}\backslash (N,N);\\ {\widehat{p}}_{NN,T+1|T}=\frac{1}{1+{\sum }_{(h,k)\in \mathcal{J}\backslash (N,N)}{\mathrm{e}}^{x_{hk,T+1|T}}},\end{array}$$

with

$$x_{ab,T+1|T}={\gamma }_{ab}+{\alpha }_{ab}{f}_{ab}(r_{i,T},r_{j,T})+{\beta }_{ab}{x}_{ab,T},(a,b)\in \mathcal{J}\backslash (N,N)$$

or its reduced form

$$x_{ab,T+1|T}=\gamma +\alpha f_{ab}(r_{i,T},r_{j,T})+\beta x_{ab,T},(a,b)\in \mathcal{J}\backslash (N,N)$$

2.1.3. Dynamic Conditional Correlation Models

In addition to CARML, the DCC model of Engle [20] is considered. In DCC models, the conditional variance matrix of the returns of k assets is decomposed as follows:

$${\mathit{H}}_t={\mathit{D}}_t{\mathit{R}}_t{\mathit{D}}_t$$

where ${\mathit{D}}_t=diag(\sqrt{h_{11,t}},\dots ,\sqrt{h_{kk,t}})$ and ${\mathit{R}}_t$ is the positive definite conditional correlation matrix. In this way, if ${\rho }_{ij,t}$ is the element of position $i,j$ of the correlation matrix ${\mathit{R}}_t$, then the corresponding element of ${\mathit{H}}_t$ is found to be ${\rho }_{ij,t}\sqrt{h_{ii,t}{h}_{jj,t}}$.

Let ${\widehat{u}}_{i,t}$ be the time t residual from the mean equation of asset i. We denote by ${\mathbf{\eta }}_t={({\eta }_{1,t},\dots ,{\eta }_{k,t})}^{\prime }$ the marginally standardized innovation vector:

$${\eta }_{i,t}=\frac{{\widehat{u}}_{i,t}}{\sqrt{h_{ii,t}}}i=1,\dots ,k.$$

Then, ${\mathit{R}}_t$ is the covariance matrix of ${\mathbf{\eta }}_t$.

Engle [20] proposed modeling the correlation matrix as

$$\begin{array}{c}{\mathit{Q}}_t=(1-a-b)\overline{\mathit{Q}}+a{\mathbf{\eta }}_{t-1}{\mathbf{\eta }}_{t-1}^{\prime }+b{\mathit{Q}}_{t-1}\\ {\mathit{R}}_t={\mathit{J}}_t{\mathit{Q}}_t{\mathit{J}}_t,\end{array}$$

where $\overline{\mathit{Q}}$ is the unconditional covariance matrix of ${\mathbf{\eta }}_t$, a and b are non-negative real numbers satisfying $0<a+b<1$ and ${\mathit{J}}_t=diag(q_{11,t}^{-1/2},\dots ,q_{kk,t}^{-1/2})$, where $q_{ii,t}$ denotes the element of position $(i,i)$ of ${\mathit{Q}}_t$.

We estimate bivariate ($k=2$) DCC models (i.e., one for each pair of commodity returns) based on the bivariate normal distribution, ${(r_{i,t},r_{j,t})}^{\prime }\sim N(\mathbf{\mu },{\mathit{H}}_t)$. The vector of means $\mathbf{\mu }$ is estimated using the sample means and the individual volatilities are assumed to follow a GARCH(1,1) model. Once the DCC model for the pair $(r_{i,t},r_{j,t})$ has been estimated, it is possible to calculate the probabilities appearing in the different types of Gerber statistics. For instance, the forecast for time $T+1$ for $p_{UD}$ is obtained as:

$${\widehat{p}}_{UD,T+1|T}={\int }_{Q_i}^{+\infty }{\int }_{-\infty }^{-Q_j}f(\mathit{x},\mathbf{\mu },{\widehat{\mathit{H}}}_{T+1|T})\mathrm{d}\mathit{x}$$

(9)

where $f(\mathit{x},\mathbf{\mu },{\widehat{\mathit{H}}}_{T+1|T})$ denotes the bivariate normal density with vector of means $\mathbf{\mu }$ and variance matrix ${\widehat{\mathit{H}}}_{T+1|T}$. Similarly,

$${\widehat{p}}_{NU,T+1|T}={\int }_{-Q_i}^{Q_i}{\int }_{Q_j}^{\infty }f(\mathit{x},\mathbf{\mu },{\widehat{\mathit{H}}}_{T+1|T})\mathrm{d}\mathit{x}.$$

(10)

Formulae (9) and (10) make use of ${\widehat{\mathit{H}}}_{T+1|T}$, the predicted variance matrix for time $T+1$ based on the information up to time T. It is also possible to use, for $t=1,\dots ,T$, formulae such as the two above with ${\widehat{\mathit{H}}}_{t|T}$ instead of ${\widehat{\mathit{H}}}_{T+1|T}$ to obtain in-sample estimates of the probabilities appearing in Equation (2) or Equation (3).

2.1.4. Filtered Historical Simulation

Filtered historical simulation is a semi-parametric method that combines HS and a parametric model (the normal-DCC model in our case). The first step consists in estimating the same DCC model of the previous section. Let

$$\left(\begin{array}{c}{z}_{i,t}\\ z_{j,t}\end{array}\right)={\mathit{H}}_t^{-1/2}\left(\left(\begin{array}{c}{r}_{i,t}\\ r_{j,t}\end{array}\right)-\mathbf{\mu }\right),t=1,\dots ,T$$

be the residuals from the model. Here, ${\mathit{H}}_t^{-1/2}$ denotes the inverse of the square-root matrix of ${\mathit{H}}_t$, obtained via the Cholesky decomposition.

Next, we obtain the thresholds, $K_i$ and $K_j$, from the residuals based on the predicted covariance matrix for time $T+1$:

$$\left(\begin{array}{c}{K}_i\\ K_j\end{array}\right)={\widehat{\mathit{H}}}_{T+1}^{-1/2}\left(\left(\begin{array}{c}{Q}_i\\ Q_j\end{array}\right)-\mathbf{\mu }\right).$$

Finally, we can use the same formulae used in the context of HS, but based on the residuals $\mathit{z}$ and the new calculated thresholds $\mathit{K}$:

$$\begin{array}{c}{\widehat{p}}_{UD,T+1|T}=\frac{1}{T}\sum _{t=1}^{T}I\left(z_{i,t}\ge K_i\right)I\left(z_{j,t}\le -K_j\right)\\ {\widehat{p}}_{NU,T+1|T}=\frac{1}{T}\sum _{t=1}^{T}I\left(-K_i<z_{i,t}<K_i\right)I\left(z_{j,t}\ge K_j\right).\end{array}$$

The remaining probabilities are predicted in a similar way.

2.1.5. The Time-Varying Gerber Correlation

Once one of the model of the previous section has been estimated, the $T+1$-Gerber statistic is obtained as

$${\widehat{g}}_{T+1|T}=\frac{{\widehat{p}}_{UU,T+1|T}+{\widehat{p}}_{DD,T+1|T}-{\widehat{p}}_{UD,T+1|T}-{\widehat{p}}_{DU,T+1|T}}{{\widehat{p}}_{UU,T+1|T}+{\widehat{p}}_{DD,T+1|T}+{\widehat{p}}_{UD,T+1|T}+{\widehat{p}}_{DU,T+1|T}}$$

(11)

or

$${\widehat{g}}_{T+1|T}^{\left(2\right)}=\frac{{\widehat{p}}_{UU,T+1|T}+{\widehat{p}}_{DD,T+1|T}-{\widehat{p}}_{UD,T+1|T}-{\widehat{p}}_{DU,T+1|T}}{\sqrt{{\widehat{p}}_{T+1|T}^{\left(A\right)}{\widehat{p}}_{T+1|T}^{\left(B\right)}}}.$$

(12)

The above formulae are the ones that we use in Section 3.3. In Section 3.2, we instead use (for the CARML and the DCC models) formulae such as the two above, but with ${\widehat{p}}_{ab,T+1|T}$ replaced by ${\widehat{p}}_{ab,t|T}$, for $t=1,\dots ,T$.

Since it is important to provide some measure of uncertainty regarding the estimated statistics (11) and (12), here we give some indication on how to derive the confidence intervals for this measure. One possibility is to use the delta method in conjunction with the covariance matrix of the estimated parameters. The issue with this approach is that the derivatives needed in the delta method are not straightforward to calculate. An alternative approach consists in computing joint confidence sets for the model parameters as done for instance in [21] for stochastic volatility models and in [22] for GARCH models with heavy-tailed innovations. This method constructs simultaneous confidence sets for the parameters of the model by numerically “inverting” some test statistic such as the likelihood ratio test. This is achieved by collecting the parameter values from a grid of admissible candidate values that are not rejected by the test. From the joint confidence set, it is then possible to derive confidence intervals for general functions of the vector of parameters and the observed data using the projection technique of Abdelkhalek and Dufour [23]. This second approach is, however, computationally intensive, given the large number of parameters in the model. For instance, a strategy to reduce the number of parameters in the DCC model could be to focus only on the parameters a and b appearing in the dynamics of the ${\mathit{Q}}_t$ matrix and treat the parameters of the univariate models as nuisance parameters.

3. Empirical Analysis

3.1. Data Description

To measure the degree of dependence across commodities, we compute the daily log-returns for the futures prices of the commodities listed in

Table 1. Commodity futures (Bloomberg Tickers).

Selected Commodities
Ticker	Description
C1 Comdty	Generic 1st Corn No. 2 Yellow futures, US$
S1 Comdty	Generic 1st Soybean No. 2 Yellow futures, US$
W1 Comdty	Generic 1st Wheat futures, US$
CL1 Comdty	WTI crude oil
NG1 Comdty	Natural Gas
GC1 Comdty	Gold
SI1 Comdty	Silver
Remaining Agricultural Commodities
Ticker	Description
KC1 Comdty	Generic 1st Coffee futures contract
SB1 Comdty	Generic 1st Sugar No. 11 (raw) futures
RR1 Comdty	Generic 1st Rice futures
CC1 Comdty	Generic 1st Cocoa
Remaining Energy Commodities
Ticker	Description
CO1 Comdty	Brent Oil
HO1 Comdty	Heating oil
Remaining Metals
Ticker	Description
HG1 Comdty	Copper

. The 14 commodities include major metals, energy and agricultural products. The first panel of Table 1 reports the commodities for which detailed results are discussed in Section 3.2 and Section 3.3. Commodity log-returns are obtained as the log differences of futures prices and are based on first generic futures contracts series extracted from Bloomberg. Specifically, we consider, at each date, the price of the contract closest to maturity. When a given contract approaches the expiration date, Bloomberg calculates a weighted average of the prices of two consecutive contracts, hence making a smooth transition between the contracts. Data cover the period from 11 August 2006 to 4 March 2020 for a total of 3500 observations for each commodity.

The descriptive statistics for the considered log-returns of the 14 commodities are reported in

Table 2. Descriptive statistics for the log-returns of commodities futures.

	Selected Commodities
	C1	S1	W1	CL1	NG1	GC1	SI1
Min.	−0.081	−0.073	−0.098	−0.119	−0.186	−0.098	−0.195
1st Qu	−0.008	−0.007	−0.010	−0.010	−0.015	−0.004	−0.007
Median	0.000	0.000	0.000	0.000	−0.001	0.001	0.001
Mean	0.000	0.000	0.000	0.000	0.000	0.000	0.000
3rd Qu	0.009	0.007	0.010	0.010	0.014	0.006	0.009
Max.	0.086	0.065	0.086	0.135	0.166	0.086	0.124
Std. Dev.	0.016	0.013	0.018	0.020	0.026	0.011	0.019
Skewness	−0.043	−0.252	0.034	0.030	0.141	−0.475	−0.916
Kurtosis	5.411	5.818	5.033	7.013	5.862	10.335	11.132
JB stat.	1214.252	1710.147	863.537	3360.445	1726.199	11,407.754	14,491.883
JB pval.	0	0	0	0	0	0	0
Q1	652.816	667.872	577.179	765.575	796.927	777.274	846.418
Q5	961.155	1241.034	842.120	1436.187	1068.713	935.768	951.451
Q1 pval.	0	0	0	0	0	0	0
Q5 pval.	0	0	0	0	0	0	0
LM1	5073.130	5761.032	4919.449	6467.057	4966.216	8230.659	7457.130
LM5	1567.114	1661.710	1531.263	1845.515	1577.667	2625.688	2431.796
LM1 pval.	0	0	0	0	0	0	0
LM5 pval.	0	0	0	0	0	0	0
	Remaining Agricultural Commodities				Remaining Energy Commodities		Remaining Metals
	KC1	SB1	RR1	CC1	CO1	HO1	HG1
Min.	−0.111	−0.124	−0.062	−0.096	−0.103	−0.098	−0.116
1st Qu	−0.010	−0.009	−0.007	−0.009	−0.009	−0.009	−0.007
Median	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Mean	0.000	0.000	0.000	0.000	0.000	0.000	0.000
3rd Qu	0.010	0.009	0.007	0.009	0.009	0.009	0.008
Max.	0.110	0.087	0.054	0.082	0.133	0.103	0.117
Std. Dev.	0.017	0.018	0.012	0.016	0.019	0.017	0.016
Skewness	0.047	-0.313	−0.003	−0.210	0.047	0.064	−0.043
Kurtosis	5.104	6.087	4.298	5.472	7.267	6.115	7.703
JB stat.	925.756	2071.051	351.917	1312.606	3800.539	2028.498	4615.341
JB pval.	0	0	0	0	0	0	0
Q1	519.798	679.805	566.418	393.100	790.971	646.284	930.245
Q5	578.844	779.116	691.245	462.829	1514.100	1230.272	1894.352
Q1 pval.	0	0	0	0	0	0	0
Q5 pval.	0	0	0	0	0	0	0
LM1	5203.071	5511.624	3811.809	6051.277	6428.943	5983.694	6486.196
LM5	1683.965	1809.333	1238.506	1962.049	1828.758	1717.633	1737.276
LM1 pval.	0	0	0	0	0	0	0
LM5 pval.	0	0	0	0	0	0	0

Note: ‘JB stat.’ and ‘JB pval.’ denote the Jarque–Bera test statistic and p-value, respectively. ‘Q1’ and ‘Q5’ are the Ljung–Box test statistics (with 1 and 5 lags) calculated on squared returns and ‘Q1 pval.’, ‘Q5 pval.’ are their p-values. ‘LM1’ and ‘LM5’ are the test statistics (with 1 and 5 lags) for Engle’s ARCH test and ‘LM1 pval.’ and ‘LM5 pval.’ are their p-values.

. For the full period, the mean daily returns are all positive but extremely small. Oil and natural gas are the most volatile return series, with a daily standard deviation of 2.0% and 2.6%, respectively. In addition, commodity returns are negatively skewed and leptokurtic. Hence, the null of normality is strongly rejected by the Jarque–Bera tests. Furthermore, we assess if returns are heteroskedastic by applying the ARCH test ([24]) and the Ljung–Box test to squared returns. For all commodities, the null of homoskedasticity is strongly rejected by both tests.

To calculate the Gerber statistic for the pair of returns $(r_i,r_j)$, we consider two alternative pairs of thresholds: (1) $Q_i={\sigma }_i/2$ and $Q_j={\sigma }_j/2$ where ${\sigma }_i$ and ${\sigma }_j$ are the (unconditional) return volatilities; and (2) $Q_i=q_i(90\%)$ and $Q_j=q_j(90\%)$ where $q_i(90\%)$ and $q_j(90\%)$ are the (unconditional) $90\%$ quantiles of the returns.

3.2. In-Sample Analysis

In this section, we estimate the CARML model in its reduced form, Equation (5), and the DCC model for six selected pairs of commodity returns using the whole sample of 3500 observations. Based on the estimated parameters, it is then possible to obtain a time-varying version of probabilities appearing in Equation (2) or Equation (3). For instance, in the case of the CARML model, we use Equation (8) with the estimated parameters and then Equation (6) to recover the time-varying probabilities. In turn, these are used to obtain the time-varying Gerber statistic. We plot the resulting time-varying Gerber correlations of Equation (3) in Figure 1 and Figure 2.

The shaded areas in the graphs—i.e., the periods from December 2007 to June 2009 and from February 2020 to the end of the period of investigation—describe the recession phases identified by the National Bureau of Economic Research. We find evidence of a significant degree of joint movements in our commodities and document that co-movements vary over time. Figure 1 visually shows the changing patterns of co-movements on a daily basis and the significant amount of correlation pervasiveness in commodity returns, with some degree of heterogeneity. Sharp co-movements are detected for the pairs silver–gold (SI1,GC1), wheat–corn (W1,C1) and crude oil–gold (CL1,GC1).

The high co-movements between silver and gold is due to their characteristics of hedging commodities [25,26]. The study by Jaffe [27] is one of the first analyses to empirically investigate the linkages between these two precious metals from 1971 to 1987. The pair was found to have a high correlation of 0.744. Erb and Harvey [28], using data between December 1982 and May 2004, showed that silver has very low correlation with other major commodities, except gold (0.66). Our findings, based on the robust measure we consider, corroborate previous studies and suggest that, since the two precious metals display high correlation, the inclusion of both commodities in the same portfolio would be redundant from an investment perspective. An interesting result is for the pair gold and copper (not reported graphically). There is a very low positive correlation between the two commodities during the global downturn of 2007–2008, while their correlation significantly rises during normal times. This finding reflects the fact that the yellow precious metal is produced, in the great majority, for investment purposes, while the red metal is destined almost entirely for industrial usage. Hence, copper prices and returns tend to increase when economies are strong and growing, simply because there is a greater global demand [29]. Conversely, gold prices and returns would rise when economies are weak, due to a risk aversion channel.

CARML and DCC models deliver similar results for the pairs crude oil–natural gas (CL1,NG1) and for crude oil–corn (CL1,C1). The joint return dynamics between crude oil and corn increased during the worldwide financial crisis. The rise in co-movements can be explained by the fact that an increase in oil prices generates contemporaneous upsurges in the price of other commodities, such as foodstuffs, via both the cost effect on the energy intensive agriculture sector and the substitution effect due to the increasing biofuel production which utilizes corn and soybeans [30]. As an energy-intensive sector, agriculture is traditionally linked to the energy industry through its input channels. While fuel and electricity are used directly in agricultural production, fertilizers and pesticides represent the two most prominent indirect energy inputs. Through these energy input channels, higher energy prices increase the cost of producing and transporting agricultural commodities.

The case of crude oil–gold (CL1,GC1) is interesting. Specifically, during critical economic and financial phases, gold is confirmed to be a ‘safe haven’ given that the common movements between the two commodities fall significantly during the recession phases (both during the period of the global financial crash and the period of sovereign crisis).

The pair soybean–corn (S1,C1) exhibits relatively high co-movements throughout the sample period. A possible reason for the strong common movements between the two commodities could be the demand for biofuels, since corn and soybean are the main crops that are used in the production of biofuels (biodiesel and ethanol) and are good substitutes.

The result further suggests that spikes in co-movements were more marked during the global financial crisis than in other periods. It also indicates that the joint return dynamics among energy, agricultural and metal commodities is always positive (excluding very few negative relations between oil and natural gas). These positive co-movements point to a reduced possibility of diversification across commodities in the short-run.

Finally, DCC models provide more volatile Gerber estimates than CARML models owing to the fact that they do not directly produce the probabilities needed for the Gerber statistic. DCC models indeed focus on the dynamics of the conditional correlation matrix and the probabilities are subsequently obtained via the integrals reported in Section 2.1.3. For CARML models, instead, the (logit transforms of) the required probabilities are directly modeled via the recursions described by Equation (4) or (5).

3.3. Out-of-Sample Analysis

To assess the out-of-sample forecasting performances of the considered models, we estimate them recursively. We use a moving window of length 1000 and end up with $N=2500$ forecasts for the probabilities appearing in the Gerber correlations. We consider again the reduced form of CARML models, Equation (5). We compare the forecasting accuracy of different models using the Brier Scores defined as:

$$\mathrm{Brier}\mathrm{score}\left(\mathcal{M}\right)=\frac{1}{N}\sum _{t=T+1}^{T+N}\sum _{\begin{array}{c}(h,k)\in \\ \mathcal{J}\backslash (N,N)\end{array}}{\left({\widehat{p}}_{hk,t|t-1}\left(\mathcal{M}\right)-I_{hk,t}\right)}^2.$$

(13)

where ${\widehat{p}}_{hk,t|t-1}\left(\mathcal{M}\right)$ represents the forecast obtained using model $\mathcal{M}$ (conditional on the information up to $t-1$) and $I_{hk,t}$ is the time t indicator associated to the region $(h,k)$. For instance, when the investigated commodity returns are $r_i$ and $r_j$, the indicator variables for the regions $UD$ and $NU$ are

$$\begin{array}{c}{I}_{UD,t}=I\left(r_{i,t}\ge Q_i\right)I\left(r_{j,t}\le -Q_j\right)\\ I_{NU,t}=I\left(-Q_i<r_{i,t}<Q_i\right)I\left(r_{j,t}\ge Q_j\right).\end{array}$$

The Brier score is a tool commonly used to evaluate probability forecasts (see, e.g., [19,31]). Since the expected value of the score is lowest for the true probability, the Brier score represents a proper scoring rule [32]. We compare the estimated models with the reference model (HS) using the Brier skill score defined as

$$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathcal{M}\right)=\left(1-\frac{\mathrm{Brier}\mathrm{score}(\mathcal{M})}{\mathrm{Brier}\mathrm{score}\left(\mathrm{HS}\right)}\right)\times 100.$$

(14)

Models that present a Brier skill score larger than zero have a superior forecasting accuracy compared to the reference model (historical simulation). Brier skill scores can also be used to compare the three alternative models considered in the present study as models with larger Brier skill scores are more accurate.

The Brier skill scores for the selected pairs of commodities are reported in

Table 3. Brier skill scores for the selected pairs of commodities.

	${\mathit{Q}}_{\mathit{i}}={\mathit{\sigma }}_{\mathit{i}}/2$ and ${\mathit{Q}}_{\mathit{j}}={\mathit{\sigma }}_{\mathit{j}}/2$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
(GC1,CL1)	3.245	1.229	1.152
(CL1,NG1)	3.534	1.805	0.870
(GC1,SI1)	0.461	−0.512	−0.270
(CL1,C1)	3.689	1.979	1.120
(C1,W1)	2.325	0.744	−0.415
(C1,S1)	0.220	0.294	−3.201
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{90}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{90}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
(GC1,CL1)	1.553	2.827	1.149
(CL1,NG1)	0.489	2.885	−0.034
(GC1,SI1)	1.818	0.064	−0.125
(CL1,C1)	3.418	4.242	−0.373
(C1,W1)	2.833	1.425	−0.079
(C1,S1)	0.526	0.873	−2.625

Note: The thresholds correspond to half the unconditional volatilities or to the 90% quantiles of the returns. The model with the best score is in bold.

. The results indicate that CARML models perform better in term of forecasting power when the thresholds are half the volatility. DCC models are more suitable when forecasting extreme commodity returns for energy and metals. Moreover, CARML models have higher forecasting power for agricultural commodities in the case in which the thresholds are set equal to the 90% quantiles. In all cases, CARML and DCC are more accurate than filtered historical simulation models.

Table 4. Brier skill scores for all the pairs of commodities.

	${\mathit{Q}}_{\mathit{i}}={\mathit{\sigma }}_{\mathit{i}}/2$ and ${\mathit{Q}}_{\mathit{j}}={\mathit{\sigma }}_{\mathit{j}}/2$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)>0$
NUM	89	90	84
FREQ	0.978	0.989	0.923
	${\mathit{Q}}_{\mathit{i}}={\mathit{\sigma }}_{\mathit{i}}/\mathbf{2}$and${\mathit{Q}}_{\mathit{j}}={\mathit{\sigma }}_{\mathit{j}}/\mathbf{2}$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
NUM	88	91	73
FREQ	0.967	1.000	0.802
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{90}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{90}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)>0$
NUM	89	91	51
FREQ	0.978	1.000	0.560
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{90}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{90}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
NUM	89	91	51
FREQ	0.978	1.000	0.560

Note: The thresholds correspond to half the unconditional volatilities or to the 90% quantiles of the returns.

instead refers to all 91 pairs of commodities we can form. In the table, we provide the number of pairs (and frequency, i.e., the number divided by 91) for which a model has a positive Brier skill score. In the table, we also give the number (and frequency) of pairs for which a given method presents a Brier skill score larger than an alternative model. The results confirm that both the CARML and DCC models outperform FHS. CARML models perform better than DCC models when the threshold is half the unconditional volatility, whereas the opposite is true when the threshold is set equal to the 90% quantile. From the analysis of all the possible pairs of commodities, we notice that CARML models tend to be more accurate when both commodities are of the agricultural type, while DCC models are more accurate when both commodities are energy products or metals.

As a robustness check, we consider two additional thresholds for the Gerber correlation.

Table 5. Brier skill scores for the selected pairs of commodities.

	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(85\%)$ and ${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(85\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
(GC1,CL1)	3.280	2.568	2.125
(CL1,NG1)	3.489	2.885	0.897
(GC1,SI1)	1.270	0.064	0.970
(CL1,C1)	3.418	4.242	1.373
(C1,W1)	2.833	1.425	1.079
(C1,S1)	1.896	2.843	0.297
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{95}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{95}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
(GC1,CL1)	1.297	2.272	0.791
(CL1,NG1)	0.386	2.269	−0.123
(GC1,SI1)	1.192	0.164	−0.070
(CL1,C1)	2.942	3.722	−0.373
(C1,W1)	2.138	2.675	−0.079
(C1,S1)	0.397	0.973	0.059

Note: The thresholds correspond to the 85% quantiles or to the 95% quantiles of the returns. The model with the best score is in bold.

reports the Brier skill scores for the six selected pairs of commodities when the thresholds are chosen as $Q_i=q_i(85\%)$, $Q_j=q_j(85\%)$, $Q_i=q_i(95\%)$ and $Q_j=q_j(95\%)$.

Table 6. Brier skill scores for all pairs of commodities.

	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(85\%)$ and ${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(85\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)>0$
NUM	90	90	87
FREQ	0.989	0.989	0.956
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{85}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{85}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
NUM	87	91	89
FREQ	0.956	1.000	0.978
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{95}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{95}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>0$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)>0$
NUM	86	91	63
FREQ	0.945	1.000	0.692
	${\mathit{Q}}_{\mathit{i}}={\mathit{q}}_{\mathit{i}}(\mathbf{95}\%)$and${\mathit{Q}}_{\mathit{j}}={\mathit{q}}_{\mathit{j}}(\mathbf{95}\%)$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{CARML}\right)>$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)>$
	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{DCC}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$	$\mathrm{Brier}\mathrm{skill}\mathrm{score}\left(\mathrm{FHS}\right)$
NUM	54	91	91
FREQ	0.593	1.000	1.000

Note: The thresholds correspond to the 85% quantiles or to the 95% quantiles of the returns.

summarizes the results for all the pairs of commodities for these thresholds.

Table 5 and Table 6 confirm the superiority of all three models with respect to HS. Comparing the two tables, we see that, when we increase the quantile from 85% to 95%, the number of pairs for which the CARML models outperforms DCC models seems to decrease.

4. Conclusions

The present study examined the co-movement of price returns for a broad spectrum of commodities for the years 2006–2020. Common movements across agricultural, energy and metal prices returns are a concern for consumers, producers, traders, investors and governments since any shock occurring to a specific commodity can easily propagate to the other commodities causing disruptive joint events. The issue of co-movement was addressed using a time-varying Gerber correlation and a set of models comprising historical simulation, filter historical simulation, DCC and CARML models. The findings reveal that there is evidence of co-movements across commodity returns and these co-movements change over time. The correlation appears more noticeable for the pairs wheat–corn, silver–gold, gold–oil and corn–soybeans. Commodity co-movements tend to rise during extreme market distress. On the basis of the Brier scores, CARML and DCC models show stronger forecasting power than the benchmark method (HS) irrespective of the choice of the thresholds. However, CARML models show higher forecasting accuracy for agricultural commodities, while DCC models are more accurate for energy products and metals. The reason can be linked to the ability of GARCH family models to better capture heteroskedasticity in the series, a feature that is more pronounced for energy and metals than agricultural products. Our results are important for financial market participants and portfolio managers in particular for building an optimal portfolio and forecasting joint futures commodity returns. For instance, the inclusion of gold and silver in the same portfolio would be redundant from an investment perspective given their high Gerber correlation. This research could be extended to stock markets or used to analyze the co-movement among spot and futures markets.