Write the code for this in R for a 5 series using Y, I am doing price discovery forgeneralized information share (GIS) and Gonzalo and Granger permanent-temporary decomposition (PT-GG) based information share measures for the spot, futures, and three ETFs.

> # Load necessary packages
> library(urca)
> library(vars)
> library(tseries)
> library(MASS)
> library(tidyverse)
> 
> # Assuming Y is already loaded in the environment
> # Select the relevant columns
> Y <- log_data[, c("log_Spot", "log_Fut",  "log_USO", "log_OIL", "log_USL")]
> # View structure of Y
> str(Y)
tibble [3,031 × 5] (S3: tbl_df/tbl/data.frame)
 $ log_Spot: num [1:3031] 4.12 4.12 4.12 4.11 4.1 ...
 $ log_Fut : num [1:3031] 4.12 4.12 4.12 4.11 4.1 ...
 $ log_USO : num [1:3031] 4.64 4.64 4.64 4.63 4.62 ...
 $ log_OIL : num [1:3031] 5.64 5.64 5.64 5.63 5.62 ...
 $ log_USL : num [1:3031] 3.14 3.15 3.14 3.14 3.13 ...
 - attr(*, "na.action")= 'omit' Named int [1:6847] 1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "names")= chr [1:6847] "1" "2" "3" "4" ...
> 
> # Step 1: Unit root tests using Phillips-Perron (PP) test and KPSS test
> # Initialize data frames to store test results
> pp_results <- data.frame(Series=character(), Test_Statistic=numeric(), p_value=numeric(), stringsAsFactors=FALSE)
> kpss_results <- data.frame(Series=character(), Test_Statistic=numeric(), p_value=numeric(), stringsAsFactors=FALSE)
> 
> # Loop over each column in Y
> for(series in names(Y)) {
+   # Phillips-Perron test
+   pp_test <- pp.test(Y[[series]])
+   pp_results <- rbind(pp_results, data.frame(Series=series, Test_Statistic=pp_test$statistic, p_value=pp_test$p.value))
+   
+   # KPSS test
+   kpss_test <- kpss.test(Y[[series]], null="Level")
+   kpss_results <- rbind(kpss_results, data.frame(Series=series, Test_Statistic=kpss_test$statistic, p_value=kpss_test$p.value))
+ }
Warning messages:
1: In kpss.test(Y[[series]], null = "Level") :
  p-value smaller than printed p-value
2: In kpss.test(Y[[series]], null = "Level") :
  p-value smaller than printed p-value
3: In kpss.test(Y[[series]], null = "Level") :
  p-value smaller than printed p-value
4: In kpss.test(Y[[series]], null = "Level") :
  p-value smaller than printed p-value
5: In kpss.test(Y[[series]], null = "Level") :
  p-value smaller than printed p-value
> 
> # Print the results
> print("Phillips-Perron Test Results:")
[1] "Phillips-Perron Test Results:"
> print(pp_results)
                          Series Test_Statistic   p_value
Dickey-Fuller Z(alpha)  log_Spot     -10.237353 0.5387773
Dickey-Fuller Z(alpha)1  log_Fut      -9.894884 0.5578765
Dickey-Fuller Z(alpha)2  log_USO      -8.401415 0.6411661
Dickey-Fuller Z(alpha)3  log_OIL      -7.929339 0.6674934
Dickey-Fuller Z(alpha)4  log_USL      -9.306146 0.5907100
> print("KPSS Test Results:")
[1] "KPSS Test Results:"
> print(kpss_results)
              Series Test_Statistic p_value
KPSS Level  log_Spot       10.48372    0.01
KPSS Level1  log_Fut       10.61626    0.01
KPSS Level2  log_USO       24.56629    0.01
KPSS Level3  log_OIL       24.82409    0.01
KPSS Level4  log_USL       20.71650    0.01
> 
> # Test first differences
> # Initialize data frames to store test results for first differences
> pp_results_diff <- data.frame(Series=character(), Test_Statistic=numeric(), p_value=numeric(), stringsAsFactors=FALSE)
> kpss_results_diff <- data.frame(Series=character(), Test_Statistic=numeric(), p_value=numeric(), stringsAsFactors=FALSE)
> 
> # Loop over each column in Y
> for(series in names(Y)) {
+   diff_series <- diff(Y[[series]])
+   
+   # Phillips-Perron test
+   pp_test <- pp.test(diff_series)
+   pp_results_diff <- rbind(pp_results_diff, data.frame(Series=series, Test_Statistic=pp_test$statistic, p_value=pp_test$p.value))
+   
+   # KPSS test
+   kpss_test <- kpss.test(diff_series, null="Level")
+   kpss_results_diff <- rbind(kpss_results_diff, data.frame(Series=series, Test_Statistic=kpss_test$statistic, p_value=kpss_test$p.value))
+ }
Warning messages:
1: In pp.test(diff_series) : p-value smaller than printed p-value
2: In kpss.test(diff_series, null = "Level") :
  p-value greater than printed p-value
3: In pp.test(diff_series) : p-value smaller than printed p-value
4: In kpss.test(diff_series, null = "Level") :
  p-value greater than printed p-value
5: In pp.test(diff_series) : p-value smaller than printed p-value
6: In kpss.test(diff_series, null = "Level") :
  p-value greater than printed p-value
7: In pp.test(diff_series) : p-value smaller than printed p-value
8: In kpss.test(diff_series, null = "Level") :
  p-value greater than printed p-value
9: In pp.test(diff_series) : p-value smaller than printed p-value
10: In kpss.test(diff_series, null = "Level") :
  p-value greater than printed p-value
> 
> # Print the results for first differences
> print("Phillips-Perron Test Results for First Differences:")
[1] "Phillips-Perron Test Results for First Differences:"
> print(pp_results_diff)
                          Series Test_Statistic p_value
Dickey-Fuller Z(alpha)  log_Spot      -3041.242    0.01
Dickey-Fuller Z(alpha)1  log_Fut      -3164.641    0.01
Dickey-Fuller Z(alpha)2  log_USO      -3217.069    0.01
Dickey-Fuller Z(alpha)3  log_OIL      -3212.187    0.01
Dickey-Fuller Z(alpha)4  log_USL      -3208.247    0.01
> print("KPSS Test Results for First Differences:")
[1] "KPSS Test Results for First Differences:"
> print(kpss_results_diff)
              Series Test_Statistic p_value
KPSS Level  log_Spot     0.05516311     0.1
KPSS Level1  log_Fut     0.05715371     0.1
KPSS Level2  log_USO     0.10059645     0.1
KPSS Level3  log_OIL     0.10775410     0.1
KPSS Level4  log_USL     0.07285208     0.1
> 
> # Step 2: Determine the optimal lag length using VARselect
> VARselect_result <- VARselect(Y, lag.max=10, type="const")
> print("Optimal lag length selection:")
[1] "Optimal lag length selection:"
> print(VARselect_result$selection)
AIC(n)  HQ(n)  SC(n) FPE(n) 
    10      9      5     10 
> 
> # Choose the number of lags based on AIC or other criteria
> # We'll use the lag order suggested by SC
> p <- VARselect_result$selection["SC(n)"]
> 
> # Step 3: Johansen cointegration test
> # Since the Johansen test in urca requires the number of lags in levels,
> # we use K = p, where p is the number of lags in the VAR
> johansen_test <- ca.jo(Y, type="eigen", ecdet="none", K= p - 1, spec="transitory")
> summary(johansen_test)

###################### 
# Johansen-Procedure # 
###################### 

Test type: maximal eigenvalue statistic (lambda max) , with linear trend 

Eigenvalues (lambda):
[1] 1.608026e-01 3.401796e-02 1.535435e-02 5.966332e-03
[5] 1.513677e-05

Values of teststatistic and critical values of test:

           test 10pct  5pct  1pct
r <= 4 |   0.05  6.50  8.18 11.65
r <= 3 |  18.11 12.91 14.90 19.19
r <= 2 |  46.84 18.90 21.07 25.75
r <= 1 | 104.76 24.78 27.14 32.14
r = 0  | 530.66 30.84 33.32 38.78

Eigenvectors, normalised to first column:
(These are the cointegration relations)

              log_Spot.l1 log_Fut.l1 log_USO.l1 log_OIL.l1
log_Spot.l1  1.0000000000  1.0000000    1.00000  1.0000000
log_Fut.l1  -1.0043201907 -1.5527447  -14.91075 -0.2909075
log_USO.l1  -0.0028367998 -1.1727863   57.02780  4.0308423
log_OIL.l1   0.0032555363  0.7267839  -57.74833 -2.0965395
log_USL.l1   0.0007628293  0.9123068   26.28275 -3.1961322
            log_USL.l1
log_Spot.l1   1.000000
log_Fut.l1    1.543037
log_USO.l1   33.235074
log_OIL.l1  -10.103925
log_USL.l1  -23.529305

Weights W:
(This is the loading matrix)

           log_Spot.l1 log_Fut.l1    log_USO.l1
log_Spot.d -0.28517146 0.02907439 -3.985293e-04
log_Fut.d   0.15789943 0.03691442 -3.836727e-04
log_USO.d  -0.16255981 0.04662993 -5.139967e-04
log_OIL.d  -0.22665865 0.04940787 -8.929494e-05
log_USL.d  -0.08940389 0.02500726 -6.971265e-04
              log_OIL.l1   log_USL.l1
log_Spot.d -5.397166e-03 1.398140e-05
log_Fut.d  -4.596148e-03 1.349959e-05
log_USO.d  -1.476627e-03 1.089308e-05
log_OIL.d   5.685786e-05 1.303423e-05
log_USL.d   1.197393e-03 1.081638e-05

> 
> # Determine the number of cointegrating vectors (r)
> # From the summary, identify the rank r where the test statistic
> # is less than the critical value. Suppose r = 4
> r <- 4
> 
> # Step 4: Estimate VECM and transform to VAR representation
> vecm_var <- vec2var(johansen_test, r=r)
> 
> # Step 5: Extract beta and alpha matrices
> # Extract the first r columns (cointegrating vectors) and exclude the constant term
> beta <- johansen_test@V[1:ncol(Y), 1:r]
> alpha <- johansen_test@W[, 1:r]
> 
> # Step 6: Compute orthogonal complements beta_perp and alpha_perp using SVD
> # Since we have (n - r) = 1 common stochastic trend, beta_perp and alpha_perp are vectors
> t_beta_svd <- svd(t(beta))
> beta_perp <- t_beta_svd$v[, ncol(t_beta_svd$v)]  # Last column corresponds to smallest singular value
> 
> t_alpha_svd <- svd(t(alpha))
> alpha_perp <- t_alpha_svd$v[, ncol(t_alpha_svd$v)]  # Similarly for alpha
> 
> # Normalize beta_perp and alpha_perp
> beta_perp <- beta_perp / beta_perp[1]
> alpha_perp <- alpha_perp / alpha_perp[1]
> 
> # Step 7: Compute long-run impact matrix Psi(1)
> # Since beta_perp and alpha_perp are vectors, Psi(1) is computed as:
> Psi_1 <- beta_perp %*% t(alpha_perp) / as.numeric(t(alpha_perp) %*% beta_perp)
> 
> # Step 8: Compute residual covariance matrix Omega
> residuals_var <- residuals(vecm_var)
> Omega <- cov(residuals_var)
> 
> # Step 9: Compute residual correlation matrix and its eigenvalues and eigenvectors
> corr_matrix <- cov2cor(Omega)
> eigen_result <- eigen(corr_matrix)
> Lambda <- diag(eigen_result$values)
> G <- eigen_result$vectors
> 
> # Step 10: Compute Cholesky decomposition of Omega
> V <- t(chol(Omega))
> 
> # Step 11: Compute F^M matrix
> Lambda_inv_sqrt <- diag(1 / sqrt(eigen_result$values))
> temp_matrix <- G %*% Lambda_inv_sqrt %*% t(G)
> V_inv <- solve(V)
> F_M_inv <- temp_matrix %*% V_inv
> F_M <- solve(F_M_inv)
> 
> # Step 12: Compute psi_1^G
> psi_1_r <- matrix(Psi_1[1, ], ncol = 1)  # Ensure it's a column vector (5x1)
> psi_1_G <- t(psi_1_r) %*% F_M  # Now dimensions should align
> 
> # Step 13: Compute numerator and denominator for GIS calculation
> numerator <- as.vector(psi_1_G^2)
> denominator <- as.numeric(t(psi_1_r) %*% Omega %*% psi_1_r)
> 
> 
> # Step 14: Compute GIS measures
> GIS_measures <- numerator / denominator
> 
> # Normalize GIS measures to sum to 1
> GIS_measures <- GIS_measures / sum(GIS_measures)
> 
> # Assign names to GIS measures
> names(GIS_measures) <- colnames(Y)
> print("Generalized Information Share (GIS) Measures:")
[1] "Generalized Information Share (GIS) Measures:"
> print(GIS_measures)
  log_Spot    log_Fut    log_USO    log_OIL    log_USL 
0.29548849 0.18456008 0.15505181 0.08007934 0.28482028 
> 
> # Step 15: Compute PT-GG information shares
> # mu is proportional to alpha_perp
> mu <- alpha_perp
> 
> # Check for negative elements in mu
> if(any(mu < 0)) {
+   # Apply equation (15) from the document
+   mu_star <- mu + abs(min(mu))
+   PTGG <- mu_star / sum(mu_star)
+   print("PT-GG Information Shares (Adjusted mu):")
+   PTGG_adjusted <- PTGG
+   names(PTGG_adjusted) <- colnames(Y)
+   print(PTGG_adjusted)
+ } else {
+   # Apply equation (13) from the document
+   PTGG <- mu / sum(mu)
+   print("PT-GG Information Shares:")
+   PTGG_normal <- PTGG
+   names(PTGG_normal) <- colnames(Y)
+   print(PTGG_normal)
+ }
[1] "PT-GG Information Shares (Adjusted mu):"
 log_Spot   log_Fut   log_USO   log_OIL   log_USL 
0.2160799 0.1696278 0.2091605 0.0000000 0.4051318 
> 
> # Alternatively, use absolute values (equation 14)
> PTGG_abs <- abs(mu) / sum(abs(mu))
> names(PTGG_abs) <- colnames(Y)
> print("PT-GG Information Shares (Absolute mu):")
[1] "PT-GG Information Shares (Absolute mu):"
> print(PTGG_abs)
   log_Spot     log_Fut     log_USO     log_OIL 
0.087591745 0.007723486 0.073393823 0.355783118 
    log_USL 
0.475507828 
> 
> # Check for negative elements in mu
> if(any(mu < 0)) {
+   # Apply equation (15) from the document
+   mu_star <- mu + abs(min(mu))
+   PTGG <- mu_star / sum(mu_star)
+   print("PT-GG Information Shares (Adjusted mu):")
+   PTGG_adjusted <- PTGG
+   names(PTGG_adjusted) <- colnames(Y)
+   print(PTGG_adjusted)
+ } else {
+   # Apply equation (13) from the document
+   PTGG <- mu / sum(mu)
+   print("PT-GG Information Shares:")
+   PTGG_normal <- PTGG
+   names(PTGG_normal) <- colnames(Y)
+   print(PTGG_normal)
+ }
[1] "PT-GG Information Shares (Adjusted mu):"
 log_Spot   log_Fut   log_USO   log_OIL   log_USL 
0.2160799 0.1696278 0.2091605 0.0000000 0.4051318 
> 
> # Alternatively, use absolute values (equation 14)
> PTGG_abs <- abs(mu) / sum(abs(mu))
> names(PTGG_abs) <- colnames(Y)
> print("PT-GG Information Shares (Absolute mu):")
[1] "PT-GG Information Shares (Absolute mu):"
> print(PTGG_abs)
   log_Spot     log_Fut     log_USO     log_OIL 
0.087591745 0.007723486 0.073393823 0.355783118 
    log_USL 
0.475507828 

Write the code for this in R:

\documentclass{article}
\usepackage{amsmath}

\begin{document}

Therefore, the series have the following vector error-correction (VEC) representation (Engle and Granger, 1987):

\begin{equation}
\Delta Y_t = \Pi Y_{t-1} + \sum_{i=1}^{k} A_i \Delta Y_{t-i} + \varepsilon_t, \quad \Pi = \alpha \beta^T
\tag{1}
\end{equation}

where $\beta$ and $\alpha$ are $n \times (n - 1)$ matrices of rank $n - 1$. The columns of $\beta$ consist of the $n - 1$ cointegrating vectors, and each column of $\alpha$ consists of the adjustment coefficients. The matrix $\Pi$ is decomposed in such a way that $\beta^T Y_t$ represents the vector of $n - 1$ stationary series. Let $\Omega$ denote the $n \times n$ covariance matrix of the innovation vector, i.e., $E[\varepsilon_t \varepsilon_t^T] = \Omega$. Following Stock and Watson (1988), equation (1) can be transformed into the following two equivalent vector moving average (VMA) representations (Hasbrouck, 1995):

\begin{equation}
\Delta Y_t = \Psi(L) \varepsilon_t
\tag{2}
\end{equation}

\begin{equation}
Y_t = Y_0 + \Psi(1) \sum_{i=1}^{t} \varepsilon_i + \Psi^*(L) \varepsilon_t
\tag{3}
\end{equation}

Then, the Engle-Granger representation theorem (Engle and Granger, 1987) implies the following (De Jong, 2002 and Lehmann, 2002):

\begin{equation}
\beta^T \Psi(1) = 0 \quad \text{and} \quad \Psi(1) \alpha = 0
\tag{4}
\end{equation}

Based on the above representations, $\Psi(1) \varepsilon_t$ represents the long-run impact of innovations on the unit-root series (Hasbrouck, 1995). Different information share measures considered by Hasbrouck (1995), Lien and Shrestha (2009), and Lien and Shrestha (2014) are based on this term.

\textbf{GENERALIZED INFORMATION SHARE (GIS) MEASURES}

Based on the above framework, we first discuss the GIS measure. Note that, in this study, we have five non-stationary series with four cointegrating vectors. Therefore, the cointegrating vector represented by matrix $\beta$ can be expressed as follows:

[
\beta^T =
\begin{bmatrix}
1 & -\gamma_1 & 0 & 0 & 0 \
1 & 0 & -\gamma_2 & 0 & 0 \
1 & 0 & 0 & -\gamma_3 & 0 \
1 & 0 & 0 & 0 & -\gamma_4
\end{bmatrix}
\tag{5}
]

Let $\psi_j^r$ be the $j$th row of $\Psi(1)$. Then, (n-1) cointegrating relations imply the following:

[
\psi_1^r = \gamma_{j-1} \psi_j^r, \quad j = 2, \ldots, 5 \tag{6}
]

In other words, equation (4) implies that the first row of $\Psi(1)$ is equal to the $\gamma_{j-1}$ times the $j$th row of $\Psi(1)$. Therefore, the long-run impact of innovations on the $i$th series is respectively given by

[
\psi_i^r \varepsilon_t = \psi_j^r \gamma_{i-1} \varepsilon_t, \quad i = 1, \ldots, 5 \tag{7}
]

where $\gamma_0 = 1$. When the innovations are \textit{independent} (i.e., $\Omega$ is diagonal), the variances of long-run impact on the $i$th series is given by:

[
\psi_i^r \Omega \psi_i^{rT} = \gamma_{i-1}^2 \sum_{j=1}^5 \psi_{1j}^2 \Omega_{jj} \tag{8}
]

where $\psi_{ij}$ is the $j$th element of $\psi_i^r$ and $\Omega_{ij}$ is the $(i,j)$th element of $\Omega$.

Let $S_{j,i}$ denote the contribution of the innovation of series $j$ to the total variance of the \textit{long-run impact} of innovation on series $i$. Then, we have the following:

[
S_{j,i} = \frac{\psi_{1j}^2 \Omega_{jj}}{\psi_1^r \Omega \psi_1^{rT}} \tag{9}
]

Note that $S_{j,i}$ is independent of $i$. Therefore, the contribution of the innovation of series $j$ to the total variance of the \textit{long-run impact} of innovation on any series will be the same. It is important to note that $S_{j,i}$ given by equation (9) is valid only when the innovations are independent. However, in general, the innovations are not independent. When the innovations are not independent, we diagonalize the correlation matrix and end up with the GIS measure proposed by Lien and Shrestha (2014). Let $\Lambda$ be a diagonal matrix containing the eigenvalues of the innovation correlation matrix on the diagonal, where the corresponding eigenvectors are given by the columns of matrix $G$. Then, we can calculate the information share of the $j$th series as follows:

[
S_j^G = \frac{(\psi_1^G)^2}{\psi_1^r \Omega \psi_1^{rT}} \tag{10}
]

where $\psi_1^G = \psi_1^r F^M, F^M = [G\Lambda^{-1/2} G^T V^{-1}]^{-1}$ and $\psi_j^G$ is the $j$th element of $\psi_1^G$. The information share measure given by equation (10) is referred to as the generalized information share (GIS) measure. It can be shown that the GIS measure is independent of ordering. Therefore, the GIS method leads to a unique information share, unlike the upper and lower bound for Hasbrouck IS measure.

\vspace{1cm}

\textbf{GONZALO-GRANGER PERMANENT TEMPORARY DECOMPOSITION (PT/GG) INFORMATION SHARE}

Here, we briefly describe the PT/GG method. Gonzalo and Granger (1995) propose a way of decomposing the vector of non-stationary series $Y_t$ into permanent component $f_t$ (non-stationary series) and transitory (stationary) component $\hat{Y}_t$, where the identification of these components is achieved by assuming that (i) the permanent component is a linear function of the original series and that (ii) the transitory component does not Granger cause the permanent component in the long-run.

The permanent component $f_t$ (under the linearity condition) can be written as:

[
f_t = \mu^T Y_t \tag{11}
]

where $\mu$ is an $n \times 1$ (or, 5 $\times$ 1 in this study) permanent component coefficient vector which can be shown to be orthogonal to the adjustment coefficient matrix $\alpha$, i.e., $\mu = \alpha_\perp$.

In our study, since we have five unit-root series, the permanent component, $f_t$ can be represented by:

[
f_t = \mu_1 Y_1 t + \mu_2 Y_2 t + \cdots + \mu_5 Y_5 t \tag{12}
]

If all the elements of $\mu$ are non-negative, then the PT/GG based information share of the $i$th series is given by:

[
PTGG_i = \frac{\mu_i}{\sum_{j=1}^5 \mu_j} \tag{13}
]

In the case of non-negative $\mu$, its elements correspond to the contribution of individual series to the permanent component. Therefore, the above definition of price discovery makes sense. However, there is no guarantee that all the elements of $\mu$ will be non-negative. It is important to note that the negative elements of $\mu$ may not necessarily lead to problems. For example, if some of the unit-root processes are negatively related to the common-stochastic trend, we expect the corresponding elements of $\mu$ to be negative. In such cases, where the negative elements of $\mu$ are valid, we can use the following definition of information share:

[
PTGG_i^+ = \frac{|\mu_i|}{\sum_{j=1}^5 |\mu_j|} \tag{14}
]

However, if some of the elements of $\mu$ are unexpectedly negative, one way to compute the PT/GG based information share is to replace the negative elements with zero and use equation (13). However, this will lead to all negative elements having zero information share regardless of their absolute value. Alternatively, we can add a constant number to each element of $\mu$ to obtain $\mu^*$ so that the series with the most negative element in $\mu$ will have zero information share. This method will result in the following information share measure:

[
PTGG_i^* = \frac{\mu_i^}{\sum_{j=1}^5 \mu_j^}, \quad \mu_i^* = \mu_i + (-\min {\mu_i}) \tag{15}
]

\section*{Empirical Results}

The data starts from 6 December 2007 to 31 December 2019, with a total number of 3,012 observations. We chose the starting date to be 6 December 2007 because this is the earliest date for which the data for all five prices are available. The data set includes the daily WTI crude oil spot and near month crude oil futures prices. The data set also includes prices of three crude oil future-based ETFs. We focus on oil ETFs that replicate the WTI spot Cushing through the nearby month futures contracts. The three oil ETFs that replicate the WTI spot Cushing through the nearby month futures contracts are:

\begin{itemize}
\item United States Oil Fund (USO),
\item iPath S&P GSCI Crude Oil Total Return (OIL), and
\item United States 12 Month Oil Fund (USL).
\end{itemize}

Throughout the empirical analyses, we use the logarithm of the prices. Table 1 presents detailed information on the three oil-based ETFs, including their ticker codes, inception dates, asset under management, and trading volumes. The USO accounts for around 77% of the assets under management, while OIL and USL account for 19% and 4%, respectively.

\begin{table}[ht]
\centering
\caption{Information of ETFs}
\begin{tabular}{|l|l|l|l|l|}
\hline
\textbf{Ticker Code} & \textbf{ETFs} & \textbf{Inception Date} & \textbf{Asset Under Management} & \textbf{Trading Volume} \
& & & \textbf{(in USD million)} & \ \hline
U:USO & United States Oil Fund & 10-04-2006 & 2,510 & 16,012,332 \
U:OIL & iPath S&P GSCI Crude Oil Total Return ETN & 16-08-2006 & 628.87 & 7,109,239 \
U:USL & United States 12 Month Oil Fund Limited Partnership & 06-12-2007 & 104.84 & 26,720 \ \hline
\end{tabular}
\end{table}

The price discovery is measured using GIS and PT/GG information share measures. To compute these two IS measures, we need to validate that each of the log-price series used in this study is non-stationary with a single unit-root. We use the Phillips-Perron (PP) unit-root test, where the unit-root is the null hypothesis. Table 2 presents the unit-root test results for the crude oil spot, futures, and the three ETFs.

\begin{table}[ht]
\centering
\caption{Unit Root Test for Crude Oil Spot, Futures, and ETFs}
\begin{tabular}{|l|l|l|l|l|}
\hline
\textbf{Series} & \textbf{Phillip-Perron Test} & \textbf{KPSS Test} \ \hline
& Log of Price & First Difference & Log of Price & First Difference \ \hline
Spot & -2.150 & -58.935*** & 3.592*** & 0.0481 \
Futures & -2.195 & -58.052*** & 3.611*** & 0.0487 \
USO & -1.661 & -56.787*** & 3.583*** & 0.0490 \
OIL & -1.493 & -56.979*** & 3.579*** & 0.0485 \
USL & -1.525 & -56.004*** & 3.617*** & 0.0495 \ \hline
\end{tabular}
\end{table}

***, **, * indicate the test statistic to be significant at the 1%, 5%, and 10% significance levels, respectively.

All the PP test statistics for log-price series are insignificant, even at the 10% level. Also, all the PP test statistics for the first differenced log-price series are significant even at the 1% level. These results indicate that all the series are integrated of order one, i.e., each of the series consists of a single unit-root. To check the robustness of our test results, we also apply the KPSS test, where the null hypothesis is stationarity. The KPSS test statistics for log prices are highly significant even at the 1% level, which rejects the null hypothesis of stationarity.

However, none of the KPSS statistics is significant for the first differenced series even at the 10% level. Therefore, both PP and KPSS tests reveal a consistent result where each of the series consists of a single unit-root.

We then proceed to examine the number of cointegrating vectors. As discussed earlier, we are dealing with five series. Therefore, to compute GIS and PT/GG information share measures, we need to establish that there are four cointegrating vectors. We apply the Johansen (1991) cointegration test to find the number of cointegrating vectors. Table 3 reports the $\lambda_{\text{max}}$ and Trace statistics. The $\lambda_{\text{max}}$ statistics for the number of cointegrating vectors (r) less than or equal to 3 is significant even at the 1% level. However, the $\lambda_{\text{max}}$ statistics for the number of cointegrating vectors (r) less than or equal to 4 is insignificant even at the 10% level. Similarly, the Trace statistics for the number of cointegrating vectors (r) less than or equal to 3 is significant at the 1% level. Therefore, we conclude that the number of cointegrating vectors is equal to 4.

\begin{table}[ht]
\centering
\caption{Cointegration Test for Crude Oil Spot, Futures, and ETFs}
\begin{tabular}{|l|c|c|c|c|}
\hline
& \textbf{$\lambda_{\text{max}}$} & \textbf{Critical Values} & \textbf{Trace} & \textbf{Critical Values} \ \hline
r $\leq$ 4 & 3.54 & 7.52 (10%), 9.24 (5%), 12.97 (1%) & 3.54 & 7.52 (10%), 9.24 (5%), 12.97 (1%) \
r $\leq$ 3 & 25.46*** & 13.75 (10%), 15.67 (5%), 20.20 (1%) & 29.47*** & 17.85 (10%), 19.96 (5%), 24.60 (1%) \
r $\leq$ 2 & 45.69*** & 19.77 (10%), 22.00 (5%), 26.81 (1%) & 74.97*** & 19.77 (10%), 22.00 (5%), 26.81 (1%) \
r $\leq$ 1 & 176.39*** & 25.56 (10%), 28.14 (5%), 32.24 (1%) & 251.36*** & 25.56 (10%), 28.14 (5%), 32.24 (1%) \
r = 0 & 551.07*** & 31.66 (10%), 34.40 (5%), 39.79 (1%) & 802.43*** & 31.66 (10%), 34.40 (5%), 39.79 (1%) \ \hline
\end{tabular}
\end{table}

The estimated cointegrating matrix $\beta^T$ is given by:

[
\beta^T = \begin{bmatrix}
1 & -0.9999 & 0 & 0 & 0 \
0 & 1 & -0.4034 & 0 & 0 \
0 & 0 & 1 & -0.3467 & 0 \
0 & 0 & 0 & 1 & -0.6690 \
\end{bmatrix}
]

with $\gamma_1 = 0.9999$, $\gamma_2 = 0.4034$, $\gamma_3 = 0.3467$, and $\gamma_4 = 0.6690$. Even though $\gamma_1$ is approximately equal to 1, other $\gamma$s are significantly different from 1. Therefore, Hasbrouck information shares cannot be computed in this case. The GIS and PT/GG methods are the right models to use in such a case.

The GIS measures for the five series are presented in Table 4. Only approximately 41 percent of the price discovery takes place in the futures market. The remaining 59 percent of the price discovery takes place in the other four markets that include the spot and the three ETF markets. Approximately 25 percent of the price discovery occurs in the USO market. As to the remaining two ETFs, around 11 and 6 percent of the price discovery take place in the USL and OIL markets, respectively. In total, approximately 42 percent of the price discovery takes place in the three ETFs markets.

\begin{table}[ht]
\centering
\caption{Information Share}
\begin{tabular}{|l|l|l|l|}
\hline
& \textbf{Generalized Information Share (GIS)} & $\mu = \alpha_1$ & \textbf{PTGG*} \ \hline
Spot & 0.1732 & -0.0439 & 0.1276 \
Futures & 0.4045 & 0.6450 & 0.3686 \
USO & 0.2453 & 0.6438 & 0.3682 \
OIL & 0.0634 & -0.0487 & 0.0000 \
USL & 0.1136 & -0.0211 & 0.1356 \ \hline
\end{tabular}
\end{table}

In sum, we find that ETFs play a significant price discovery role, and additional ETFs play an additional price discovery role. Therefore, additional ETFs are not redundant.

\end{document}