*2.1. Data*

We select natural gas price indicators for the North American, European, and Asia–Pacific markets. We use Henry Hub (HH) futures as the North American index. HH is the name of a distribution hub in Louisiana for the natural gas pipeline system in North America, referring to natural gas delivered at that hub. HH futures are listed on the New York Mercantile Exchange. We adopt National Balancing Point (NBP) and Title Transfer Facility (TTF) futures as the European index. NBP futures are futures of natural gas at a virtual trading location in the United Kingdom, being listed on the Intercontinental Exchange (ICE). This market is related to production in the North Sea gas field, which is trading with the European continent by pipelines, and consumption within the United Kingdom. TTF futures are the futures of natural gas at a virtual trading point in the Netherlands, being listed on the ICE. This market reflects the continental European market, which consists of trading by long-haul pipelines and LNG. We utilize the Japan/Korea Marker (JKM) as the Asia-Pacific index, which is associated with the short-term trading market for LNG in the Asia–Pacific region. It is provided by Platts, which assesses benchmark prices in physical energy markets. We obtain the daily data from February 2, 2009 to February 28, 2019 from Bloomberg.

Figure 1 shows the time series for these indexes. First, HH fluctuates independently from the other variables. Second, the two European indexes move together. Finally, JKM is linked to the European market over a specific period. As a result of the Great East Japan Earthquake on 11 March 2011, the amount of power generated by nuclear power plants has decreased significantly in Japan. Thus, we can assume that the Asia-Pacific LNG market was tight, and therefore the JKM was soaring. Moreover, these European indexes might have increased as a result of arbitrage trading with the Asia-Pacific LNG market. HH's downtrend around 2012 might be caused by the expectations of increased shale gas production, while the price spike in 2014 was caused by the North American cold wave.

**Figure 1.** Time series plots of natural gas price. (HH = Henry Hub, NBP = National Balancing Point, TTF = Title Transfer Facility, JKM = Japan/Korea Marker).

Table 1 presents the summary statistics of the return and volatility series of each index. There are 2436 observations in each case. The JKM return series has a distribution biased to the right and the other series have distributions biased to the left, as skewness is negative only for the JKM return series. The kurtosis of all series are extremely large. In other words, all series have fat tail distributions. We can reject the hypothesis that each series is normally distributed by the Jarque-Bera statistics calculated from the skewness and kurtosis of each series.



\*: calculated by the stochastic volatility (SV) model.

As described in Toyoshima and Hamori [26], let *xt* be a return with mean zero and variance *exp*(*ht*). The SV model can be expressed as follows:

$$\mathbf{x}\_{t}|h\_{t}\sim\mathcal{N}(0,\,\exp(h\_{t}))$$

$$h\_{t}|\_{\mathbf{l}\_{t-1},\,\,\mu,\,\,\boldsymbol{\varrho},\,\,\boldsymbol{\sigma}\_{\tau}}\sim\mathcal{N}(\boldsymbol{\mu}+\boldsymbol{\varrho}(h\_{t-1}-\boldsymbol{\mu}),\,\boldsymbol{\sigma}\_{\tau}^{2})$$

$$h\_{0}|\boldsymbol{\mu},\,\,\boldsymbol{\varrho},\,\,\boldsymbol{\sigma}\_{\tau}\sim\mathcal{N}(\boldsymbol{\mu},\,\boldsymbol{\sigma}\_{\tau}^{2}/\left(1-\boldsymbol{\varrho}^{2}\right))$$

where μ, ϕ, and στ are the level of log variance, the persistence of log variance, and the volatility of log variance, respectively.

#### *2.2. Preliminary Analyses*

The condition for the VMA representation of the VAR model is that all variables are stationary. Accordingly, we test for the stationarity status of all series by the augmented Dickey–Fuller (ADF) test. Table 2 presents the results. The ADF test rejects the null hypothesis that all variables have a unit root. Therefore, we can utilize the VMA representation.

We estimate the coefficient of the diagonal Baba, Engle, Kraft, and Kroner (BEKK) GARCH model to calculate the dynamic correlation coefficients. The BEKK model is one of several variations of multivariate GARCH models depending on the formulation of the time-varying variance–covariance matrix. The BEKK model was developed by Engle and Kroner [28], and it guarantees a positive estimated variance. Moreover, the BEKK model can dynamically calculate the correlation coefficient. When there are *k* variables, the variance-covariance matrix *Ht* of the BEKK model is as follows:

$$H\_t = L'L + M'\varepsilon\_{t-1}\varepsilon'\_{t-1}M + N'H\_{t-1}N\tag{1}$$

where *M* and *N* are a *k* × *k* matrix, and *L* is a *k* × *k* symmetric matrix.


#### **Table 2.** Unit root tests.

We calculate the correlation coefficients over the entire period, namely from 2 February 2009 to 28 February 2019, to understand the simultaneous relationship between variables as a phenomenon. Tables 3 and 4 present the correlation coefficients for the return and volatility series, respectively. The correlation coefficient between the NBP and TTF is the largest for both the return and volatility series. Furthermore, the correlation coefficients of the volatility series are larger than those of the return series under each combination of variables. This fact indirectly indicates the high sustainability of volatility.


**Table 3.** Correlation coefficient of return series.


Figure 2 provides the time series of the dynamic correlation coefficients between each pair of return series. Although the NBP and TTF have the highest correlation coefficient in Table 4, even that dynamic correlation coefficient temporarily decreases to 0. Conversely, even if a correlation coefficient over the entire period is low, the dynamic correlation coefficient may temporarily increase.

**Figure 2.** Dynamic correlation between each return series.

Figure 3 provides the time series of the dynamic correlation coefficients between each volatility series. We find that the correlation and inverse correlation phases are clear because volatility accumulates at a high rate. However, only the dynamic correlation coefficient between the NBP and TTF continues to be almost +1.

**Figure 3.** Dynamic correlation between each volatility series.
