Include the most important results in my paper:
Introduction

This paper examines the dynamics between the United States 12-Month Oil Fund (USL) and futures oil prices (FUT) using time series econometric techniques. Focusing on volatility modeling and causal relationships, we aim to provide insights into how oil price movements impact market conditions and policy decisions. Given oil's central role in economic stability and geopolitical strategies, this analysis is valuable for traders, policymakers, and risk managers seeking to understand market behavior and improve price discovery processes.

We apply advanced econometric tools, including Vector Autoregression (VAR), Forecast Error Variance Decomposition (FEVD), and Generalized Autoregressive Conditional Heteroskedasticity (GARCH), to explore interdependencies, volatility clustering, and the dynamic responses of oil prices to market shocks. These methods enable a comprehensive analysis of the relationships between USL and FUT, contributing to both academic research and practical decision-making in volatile energy markets.

Literature Review

Key Studies on Oil Price Modeling

Oil price modeling has been extensively studied due to its economic importance. Traditional approaches, such as ARIMA and VAR models, forecast oil price movements and their linkages to macroeconomic variables. Structural models analyze long-run equilibria between oil prices and global economic indicators.

VAR Models in Financial and Commodity Markets

VAR models are well-established in capturing interdependencies among financial and commodity market variables. They are particularly useful for analyzing dynamic relationships between different forms of oil prices (e.g., futures and ETFs like USL) and broader macroeconomic variables such as exchange rates and inflation.

Volatility Modeling Using GARCH

GARCH models capture volatility clustering in financial time series, a prominent feature in oil markets. These models highlight the persistence of volatility and the impact of external shocks, making them relevant for analyzing oil price dynamics.

Gaps in Literature

While oil price forecasting and volatility modeling are well-researched, limited focus exists on the interplay between futures oil prices (FUT) and medium-term price proxies like USL, particularly through combined VAR and GARCH modeling. This paper addresses this gap by examining these relationships and their economic implications, providing a novel contribution to the field.

Methods and Data

Data Description and Transformation

We used daily prices for USL (a mix of futures oil contracts ETF) and FUT (futures oil prices) from December 6, 2007, to December 30, 2019, sourced from www.investing.com. After merging by date, the dataset contains 3,148 observations. We computed the natural logarithm of both USL and FUT prices to stabilize variance and then calculated first differences of the log-transformed data to obtain returns, ensuring stationarity.

Visual inspection revealed substantial volatility and significant spikes during periods of economic uncertainty, such as the 2008 financial crisis. The differenced log data, representing returns, appeared stationary and suitable for time series modeling.

Econometric Methods

To analyze the dynamic relationships between USL and FUT prices, we employed several econometric techniques. Stationarity was tested using Dickey-Fuller tests, confirming that the log-transformed data were non-stationary, while the differenced log data were stationary.

A VAR model was estimated to examine interdependencies between USL and FUT returns, with six lags selected based on statistical criteria. The VAR results revealed significant interrelationships between the two variables. FEVD quantified the extent to which each variable's variance was attributable to shocks in the other, highlighting the strong influence of FUT on USL.

Impulse Response Functions (IRF) analyzed the effects of shocks to FUT and USL over time, showing that shocks to FUT had a significant and immediate impact on USL, reinforcing the futures market's role in leading price discovery. Lastly, GARCH(1,1) models were applied to capture persistent volatility and clustering effects in both FUT and USL.

Empirical Results

Data Visualization

The raw data for USL and FUT from 2009 to 2020 show overall trends and volatility of oil prices, highlighting key economic events such as the 2008 financial crisis. The strong co-movement between the series reflects their economic linkage. FUT appears slightly more volatile than USL due to its sensitivity to market expectations.

The logarithmic transformation stabilizes variance and allows interpreting results in terms of percentage changes. The first differences of the logarithmic data, representing daily returns, appear stationary and are suitable for econometric modeling like VAR and GARCH.

Stationarity Analysis

Augmented Dickey-Fuller (ADF) tests indicate that the log-transformed levels of FUT and USL are non-stationary, while the first-differenced log data (returns) are stationary. This confirms that returns are mean-reverting and appropriate for time series modeling.

VAR Model

The optimal lag length for the VAR model was determined to be six, capturing immediate interdependencies between FUT and USL. The VAR(6) model shows significant lagged relationships between the returns of FUT and USL, indicating bi-directional influence.

Hypothesis tests confirm that the coefficients of lagged variables are jointly significant for both FUT and USL, reinforcing the necessity of including lagged terms to adequately model dynamic interrelations.

Residual diagnostics confirm that the model is stable, with roots of the characteristic polynomial inside the unit circle.

Forecast Error Variance Decomposition (FEVD)

FEVD results highlight the dynamics between FUT and USL. When diff_log_FUT is ordered first, its variance is almost entirely explained by its own shocks, while diff_log_USL's variance is largely explained by shocks from diff_log_FUT. Conversely, when diff_log_USL is ordered first, its variance is driven by its own shocks, and a significant portion of diff_log_FUT's variance is explained by shocks to diff_log_USL.

Impulse Response Function (IRF) Results

IRF analysis shows that shocks to FUT result in a strong and immediate positive response in USL, demonstrating the significant influence of futures prices on the ETF. Shocks to USL have a smaller and short-lived impact on FUT, highlighting the futures market's superior ability to incorporate forward-looking information.

GARCH Model

The GARCH(1,1) model captures volatility clustering in the residuals of FUT and USL. Significant spikes in volatility occur during periods of market stress, such as the 2008 financial crisis. High values of the GARCH parameters indicate that volatility is highly persistent.

Standardized residuals are centered around zero, indicating that the GARCH(1,1) model effectively captures conditional heteroskedasticity in the data.

Discussion and Conclusion

Discussion

The results demonstrate a strong interdependency between FUT and USL prices, with futures markets playing a dominant role in price discovery. Shocks to FUT lead to immediate and significant impacts on USL, reflecting the forward-looking nature of futures prices. While the VAR model reveals bidirectional relationships, IRFs show that the influence of USL on FUT is weaker and short-lived.

Stationarity testing confirmed that returns are stationary, aligning with financial theory. GARCH(1,1) models revealed persistent volatility clustering in both FUT and USL returns, highlighting the sensitivity of oil markets to macroeconomic and geopolitical shocks.

Conclusion

This study confirms that FUT significantly drives USL, emphasizing the futures market's role in price discovery. VAR and IRF results reveal a strong dynamic relationship where FUT shocks propagate efficiently to USL, while GARCH models capture persistent volatility clustering. These findings contribute to understanding the interplay between immediate futures prices and ETFs representing medium-term trends, offering insights for market participants and policymakers.

Future research could expand on these results by analyzing higher-frequency data or incorporating macroeconomic variables. Alternative models, such as structural VAR or regime-switching GARCH, may further clarify the complex dynamics of oil markets.

> # Load necessary libraries
> library(readxl)
> library(dplyr)
> library(tidyr)
> library(ggplot2)
> library(MTS)
> library(data.table) 
> library(urca)
> library(vars)
> library(rugarch)
> library(car) 
> library(MCMCpack)
> library(maxLik)
> library(mvtnorm)
> library(lmtest)
> library(xts)
> setwd("C:/Users/abelma2/Desktop/R_Projects/ACE_527_Price_Discovery")
> Fut <- read_excel("OILNF_US0_NEW.xlsx", sheet = "Fut", col_types = c("date", "numeric"))
> USL <- read_excel("OILNF_US0_NEW.xlsx", sheet = "USL", col_types = c("date", "numeric"))
Warning messages:
1: Expecting numeric in B3150 / R3150C2: got '#N/A N/A' 
2: Expecting numeric in B3151 / R3151C2: got '#N/A N/A' 
3: Expecting numeric in B3152 / R3152C2: got '#N/A N/A' 
> # Merge datasets by Date using full join to include all dates
> data <- Fut %>%
+   full_join(USL, by = "Dates") %>% 
+   na.omit(.)
> range.Date(data$Dates)
[1] "2007-12-06 UTC" "2019-12-30 UTC"
> par(mfrow = c(2,1), mar = c(3, 4, 1, 1))  # Adjust margins for each plot
> names = c("Fut", "USL")
> for (name in names) {
+   plot(data$Dates, data[[name]], type = "l", 
+        ylab = "", 
+        main = paste("Plot of", name))
+ }
> # Merge the datasets by Date using left_join
> log_data <- Fut %>%
+   arrange(Dates) %>% 
+   full_join(USL, by = "Dates") %>%
+   na.omit(.) %>%  # Use na.omit() instead of drop_na()
+   mutate(
+     log_Fut = log(Fut),
+     log_USL = log(USL),
+     Time = row_number()  # Adding the Time column with incremental values
+   ) %>%
+   na.omit(.) %>% 
+   dplyr::select(Time, Dates, log_Fut, log_USL)
> par(mfrow = c(2,1), mar = c(3, 4, 1, 1))  # Adjust margins for each plot
> log_names = c("log_Fut", "log_USL")
> for (name in log_names) {
+   plot(log_data$Dates, log_data[[name]], type = "l", 
+        ylab = "", 
+        main = paste("Plot of", name))
+ }
> # Compute first differences
> diff_log_data <- log_data %>%
+   arrange(Dates) %>%
+   mutate(
+     diff_log_Fut = c(NA, diff(log_Fut)),
+     diff_log_USL = c(NA, diff(log_USL))
+   ) %>%
+   drop_na(.) %>% 
+   dplyr::select(Time, Dates, diff_log_Fut, diff_log_USL)
> # Set up plotting area for multiple plots
> par(mfrow = c(2,1), mar = c(3, 4, 1, 1))  # Adjust margins for each plot
> diff_log_names <- c("diff_log_Fut", "diff_log_USL")
> for (name in diff_log_names) {
+   plot(diff_log_data$Dates, diff_log_data[[name]], type = "l", 
+        ylab = "", 
+        main = paste("Plot of", name), 
+        xlab = "Dates")
+ }
> # Compute first differences
> diff_log_dataframe <- as.data.frame(diff_log_data)
> dev.off()
null device 
          1 
> # Dickey Fuller Test
> dicky_fuller.trend <- lapply(log_data[, c("log_Fut", "log_USL")], 
+                              function(x) ur.df(as.vector(x), type = "trend", selectlags = "AIC"))
> dicky_fuller.constant <- lapply(log_data[, c("log_Fut", "log_USL")], 
+                                 function(x) ur.df(as.vector(x), type = "drift", selectlags = "AIC"))
> dicky_fuller.none <- lapply(log_data[, c("log_Fut", "log_USL")], 
+                             function(x) ur.df(as.vector(x), type = "trend", selectlags = "AIC"))
> lapply(dicky_fuller.trend, summary)
$log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.132672 -0.011079  0.000854  0.011452  0.157751 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.454e-02  6.804e-03   2.137   0.0327 * 
z.lag.1     -3.245e-03  1.484e-03  -2.186   0.0289 * 
tt          -5.481e-07  5.420e-07  -1.011   0.3120   
z.diff.lag  -5.053e-02  1.782e-02  -2.836   0.0046 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02362 on 3142 degrees of freedom
Multiple R-squared:  0.004243,	Adjusted R-squared:  0.003292 
F-statistic: 4.463 on 3 and 3142 DF,  p-value: 0.003915


Value of test-statistic is: -2.1862 1.627 2.3997 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


$log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106327 -0.008434  0.000813  0.009516  0.084499 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  9.926e-03  5.406e-03   1.836   0.0664 .
z.lag.1     -2.593e-03  1.344e-03  -1.929   0.0538 .
tt          -7.918e-07  5.879e-07  -1.347   0.1782  
z.diff.lag  -3.500e-02  1.783e-02  -1.963   0.0497 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01826 on 3142 degrees of freedom
Multiple R-squared:  0.00253,	Adjusted R-squared:  0.001578 
F-statistic: 2.657 on 3 and 3142 DF,  p-value: 0.04682


Value of test-statistic is: -1.9289 1.4759 1.9049 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


> lapply(dicky_fuller.constant, summary)
$log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.131682 -0.011082  0.000761  0.011450  0.159001 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.010373   0.005416   1.915  0.05554 . 
z.lag.1     -0.002468   0.001270  -1.943  0.05206 . 
z.diff.lag  -0.050965   0.017811  -2.862  0.00424 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02362 on 3143 degrees of freedom
Multiple R-squared:  0.003919,	Adjusted R-squared:  0.003285 
F-statistic: 6.183 on 2 and 3143 DF,  p-value: 0.002089


Value of test-statistic is: -1.9434 1.9292 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1  6.43  4.59  3.78


$log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106083 -0.008531  0.000711  0.009510  0.084921 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.0037347  0.0028437   1.313   0.1892  
z.lag.1     -0.0011581  0.0008198  -1.413   0.1579  
z.diff.lag  -0.0358395  0.0178196  -2.011   0.0444 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01826 on 3143 degrees of freedom
Multiple R-squared:  0.001954,	Adjusted R-squared:  0.001319 
F-statistic: 3.077 on 2 and 3143 DF,  p-value: 0.04622


Value of test-statistic is: -1.4127 1.3068 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1  6.43  4.59  3.78


> lapply(dicky_fuller.none, summary)
$log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.132672 -0.011079  0.000854  0.011452  0.157751 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.454e-02  6.804e-03   2.137   0.0327 * 
z.lag.1     -3.245e-03  1.484e-03  -2.186   0.0289 * 
tt          -5.481e-07  5.420e-07  -1.011   0.3120   
z.diff.lag  -5.053e-02  1.782e-02  -2.836   0.0046 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02362 on 3142 degrees of freedom
Multiple R-squared:  0.004243,	Adjusted R-squared:  0.003292 
F-statistic: 4.463 on 3 and 3142 DF,  p-value: 0.003915


Value of test-statistic is: -2.1862 1.627 2.3997 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


$log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106327 -0.008434  0.000813  0.009516  0.084499 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  9.926e-03  5.406e-03   1.836   0.0664 .
z.lag.1     -2.593e-03  1.344e-03  -1.929   0.0538 .
tt          -7.918e-07  5.879e-07  -1.347   0.1782  
z.diff.lag  -3.500e-02  1.783e-02  -1.963   0.0497 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01826 on 3142 degrees of freedom
Multiple R-squared:  0.00253,	Adjusted R-squared:  0.001578 
F-statistic: 2.657 on 3 and 3142 DF,  p-value: 0.04682


Value of test-statistic is: -1.9289 1.4759 1.9049 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


> # Difference Data Dickey Fuller Test
> diff_df.trend <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")], 
+                         function(x) ur.df(as.vector(x), type = "trend", selectlags = "AIC"))
> diff_df.constant <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")],
+                            function(x) ur.df(as.vector(x), type = "drift", selectlags = "AIC"))
> diff_df.none <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")], 
+                        function(x) ur.df(as.vector(x), type = "none", selectlags = "AIC"))
> lapply(diff_df.trend , summary)
$diff_log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.128257 -0.011036  0.000608  0.011333  0.155890 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.274e-04  8.429e-04  -0.270  0.78737    
z.lag.1     -1.102e+00  2.585e-02 -42.618  < 2e-16 ***
tt           6.542e-08  4.639e-07   0.141  0.88786    
z.diff.lag   4.711e-02  1.782e-02   2.644  0.00824 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02362 on 3141 degrees of freedom
Multiple R-squared:  0.5272,	Adjusted R-squared:  0.5267 
F-statistic:  1167 on 3 and 3141 DF,  p-value: < 2.2e-16


Value of test-statistic is: -42.6185 605.4453 908.1679 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


$diff_log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106586 -0.008681  0.000631  0.009393  0.085398 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.206e-04  6.523e-04  -0.645    0.519    
z.lag.1     -1.030e+00  2.569e-02 -40.106   <2e-16 ***
tt           1.058e-07  3.589e-07   0.295    0.768    
z.diff.lag  -6.015e-03  1.784e-02  -0.337    0.736    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01827 on 3141 degrees of freedom
Multiple R-squared:  0.5182,	Adjusted R-squared:  0.5177 
F-statistic:  1126 on 3 and 3141 DF,  p-value: < 2.2e-16


Value of test-statistic is: -40.1064 536.1738 804.2607 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.96 -3.41 -3.12
phi2  6.09  4.68  4.03
phi3  8.27  6.25  5.34


> lapply(diff_df.constant, summary)
$diff_log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.128341 -0.011063  0.000597  0.011342  0.155806 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0001244  0.0004211  -0.295  0.76766    
z.lag.1     -1.1018014  0.0258488 -42.625  < 2e-16 ***
z.diff.lag   0.0471053  0.0178177   2.644  0.00824 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02361 on 3142 degrees of freedom
Multiple R-squared:  0.5272,	Adjusted R-squared:  0.5269 
F-statistic:  1752 on 2 and 3142 DF,  p-value: < 2.2e-16


Value of test-statistic is: -42.6249 908.4414 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1  6.43  4.59  3.78


$diff_log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106745 -0.008659  0.000601  0.009359  0.085443 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.0002541  0.0003259  -0.780    0.436    
z.lag.1     -1.0300903  0.0256809 -40.111   <2e-16 ***
z.diff.lag  -0.0060460  0.0178361  -0.339    0.735    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01827 on 3142 degrees of freedom
Multiple R-squared:  0.5182,	Adjusted R-squared:  0.5179 
F-statistic:  1690 on 2 and 3142 DF,  p-value: < 2.2e-16


Value of test-statistic is: -40.1111 804.4511 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1  6.43  4.59  3.78


> lapply(diff_df.none, summary)
$diff_log_Fut

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.128467 -0.011187  0.000472  0.011218  0.155686 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -1.10175    0.02584 -42.630  < 2e-16 ***
z.diff.lag  0.04708    0.01781   2.643  0.00827 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02361 on 3143 degrees of freedom
Multiple R-squared:  0.5272,	Adjusted R-squared:  0.5269 
F-statistic:  1752 on 2 and 3143 DF,  p-value: < 2.2e-16


Value of test-statistic is: -42.6301 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62


$diff_log_USL

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.106998 -0.008917  0.000341  0.009105  0.085184 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -1.029696   0.025674  -40.11   <2e-16 ***
z.diff.lag -0.006246   0.017833   -0.35    0.726    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01827 on 3143 degrees of freedom
Multiple R-squared:  0.5181,	Adjusted R-squared:  0.5178 
F-statistic:  1690 on 2 and 3143 DF,  p-value: < 2.2e-16


Value of test-statistic is: -40.106 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62


> optimal = VARorder(diff_log_dataframe[, c("diff_log_Fut","diff_log_USL")], maxp = 100)
selected order: aic =  80 
selected order: bic =  6 
selected order: hq =  22 
Summary table:  
         p      AIC      BIC       HQ     M(p) p-value
  [1,]   0 -16.7500 -16.7500 -16.7500   0.0000  0.0000
  [2,]   1 -16.7787 -16.7710 -16.7760  95.0670  0.0000
  [3,]   2 -16.7969 -16.7815 -16.7914  62.9809  0.0000
  [4,]   3 -16.7999 -16.7768 -16.7916  16.9747  0.0020
  [5,]   4 -16.8080 -16.7773 -16.7970  32.4076  0.0000
  [6,]   5 -16.8192 -16.7807 -16.8054  41.5388  0.0000
  [7,]   6 -16.8327 -16.7865 -16.8161  48.6636  0.0000
  [8,]   7 -16.8357 -16.7818 -16.8164  16.8449  0.0021
  [9,]   8 -16.8368 -16.7752 -16.8147  10.9595  0.0270
 [10,]   9 -16.8466 -16.7773 -16.8217  37.3495  0.0000
 [11,]  10 -16.8511 -16.7742 -16.8235  21.4362  0.0003
 [12,]  11 -16.8509 -16.7663 -16.8206   7.1599  0.1277
 [13,]  12 -16.8487 -16.7563 -16.8155   0.8028  0.9381
 [14,]  13 -16.8592 -16.7592 -16.8233  39.4787  0.0000
 [15,]  14 -16.8593 -16.7516 -16.8207   8.0655  0.0892
 [16,]  15 -16.8677 -16.7523 -16.8263  32.8795  0.0000
 [17,]  16 -16.8670 -16.7439 -16.8228   5.6517  0.2267
 [18,]  17 -16.8652 -16.7344 -16.8182   2.1150  0.7146
 [19,]  18 -16.8663 -16.7278 -16.8166  11.1239  0.0252
 [20,]  19 -16.8648 -16.7186 -16.8124   3.1554  0.5322
 [21,]  20 -16.8652 -16.7113 -16.8100   8.7858  0.0667
 [22,]  21 -16.8646 -16.7030 -16.8066   5.8255  0.2126
 [23,]  22 -16.8970 -16.7277 -16.8363 104.9180  0.0000
 [24,]  23 -16.8961 -16.7191 -16.8326   4.8184  0.3064
 [25,]  24 -16.8957 -16.7110 -16.8294   6.4107  0.1705
 [26,]  25 -16.8936 -16.7012 -16.8246   1.2921  0.8627
 [27,]  26 -16.8956 -16.6956 -16.8238  13.7252  0.0082
 [28,]  27 -16.8957 -16.6880 -16.8212   7.8870  0.0958
 [29,]  28 -16.8935 -16.6780 -16.8161   0.8049  0.9378
 [30,]  29 -16.8932 -16.6700 -16.8131   6.6982  0.1527
 [31,]  30 -16.8926 -16.6618 -16.8098   6.0703  0.1940
 [32,]  31 -16.8910 -16.6524 -16.8054   2.5795  0.6305
 [33,]  32 -16.8899 -16.6437 -16.8015   4.3895  0.3558
 [34,]  33 -16.9022 -16.6483 -16.8111  44.3111  0.0000
 [35,]  34 -16.9017 -16.6400 -16.8078   5.9068  0.2062
 [36,]  35 -16.9018 -16.6325 -16.8052   8.0750  0.0889
 [37,]  36 -16.8996 -16.6226 -16.8002   0.8952  0.9253
 [38,]  37 -16.9036 -16.6189 -16.8014  19.4636  0.0006
 [39,]  38 -16.9125 -16.6201 -16.8076  33.9413  0.0000
 [40,]  39 -16.9144 -16.6143 -16.8067  13.2798  0.0100
 [41,]  40 -16.9173 -16.6095 -16.8068  15.9442  0.0031
 [42,]  41 -16.9168 -16.6013 -16.8036   6.1810  0.1860
 [43,]  42 -16.9146 -16.5914 -16.7986   0.9979  0.9101
 [44,]  43 -16.9135 -16.5826 -16.7948   4.3855  0.3563
 [45,]  44 -16.9128 -16.5742 -16.7913   5.2466  0.2629
 [46,]  45 -16.9112 -16.5650 -16.7870   2.9815  0.5609
 [47,]  46 -16.9122 -16.5582 -16.7852  10.2469  0.0365
 [48,]  47 -16.9138 -16.5522 -16.7841  12.4503  0.0143
 [49,]  48 -16.9186 -16.5492 -16.7861  21.5136  0.0003
 [50,]  49 -16.9263 -16.5492 -16.7910  30.0773  0.0000
 [51,]  50 -16.9318 -16.5471 -16.7938  23.8803  0.0001
 [52,]  51 -16.9340 -16.5416 -16.7932  13.9777  0.0074
 [53,]  52 -16.9354 -16.5352 -16.7918  11.4445  0.0220
 [54,]  53 -16.9356 -16.5278 -16.7893   8.2329  0.0834
 [55,]  54 -16.9383 -16.5227 -16.7892  15.1948  0.0043
 [56,]  55 -16.9398 -16.5165 -16.7879  11.9075  0.0181
 [57,]  56 -16.9422 -16.5112 -16.7876  14.4935  0.0059
 [58,]  57 -16.9401 -16.5015 -16.7827   1.4072  0.8429
 [59,]  58 -16.9407 -16.4943 -16.7805   9.0484  0.0599
 [60,]  59 -16.9397 -16.4856 -16.7767   4.5029  0.3422
 [61,]  60 -16.9378 -16.4761 -16.7722   2.0657  0.7237
 [62,]  61 -16.9369 -16.4675 -16.7684   4.6784  0.3219
 [63,]  62 -16.9374 -16.4603 -16.7662   9.0409  0.0601
 [64,]  63 -16.9373 -16.4525 -16.7634   7.1520  0.1281
 [65,]  64 -16.9429 -16.4504 -16.7661  23.4895  0.0001
 [66,]  65 -16.9510 -16.4508 -16.7715  31.0555  0.0000
 [67,]  66 -16.9517 -16.4438 -16.7694   9.4799  0.0502
 [68,]  67 -16.9516 -16.4360 -16.7666   7.0597  0.1328
 [69,]  68 -16.9506 -16.4273 -16.7628   4.5732  0.3340
 [70,]  69 -16.9515 -16.4205 -16.7609   9.9364  0.0415
 [71,]  70 -16.9502 -16.4116 -16.7570   3.8405  0.4280
 [72,]  71 -16.9509 -16.4046 -16.7549   9.3914  0.0520
 [73,]  72 -16.9497 -16.3956 -16.7509   3.7202  0.4452
 [74,]  73 -16.9524 -16.3907 -16.7509  15.3714  0.0040
 [75,]  74 -16.9526 -16.3831 -16.7482   7.7822  0.0999
 [76,]  75 -16.9514 -16.3743 -16.7443   3.9613  0.4113
 [77,]  76 -16.9546 -16.3697 -16.7447  16.4973  0.0024
 [78,]  77 -16.9550 -16.3624 -16.7423   8.4799  0.0755
 [79,]  78 -16.9538 -16.3536 -16.7385   4.0897  0.3940
 [80,]  79 -16.9557 -16.3478 -16.7376  12.7633  0.0125
 [81,]  80 -16.9567 -16.3411 -16.7358  10.1130  0.0386
 [82,]  81 -16.9556 -16.3323 -16.7320   4.3413  0.3618
 [83,]  82 -16.9544 -16.3234 -16.7280   3.6966  0.4486
 [84,]  83 -16.9532 -16.3145 -16.7240   3.9598  0.4115
 [85,]  84 -16.9519 -16.3055 -16.7199   3.5175  0.4752
 [86,]  85 -16.9530 -16.2989 -16.7183  10.5424  0.0322
 [87,]  86 -16.9508 -16.2891 -16.7134   1.0846  0.8967
 [88,]  87 -16.9492 -16.2797 -16.7090   2.6556  0.6170
 [89,]  88 -16.9488 -16.2716 -16.7058   5.9786  0.2008
 [90,]  89 -16.9486 -16.2637 -16.7028   6.7379  0.1504
 [91,]  90 -16.9474 -16.2549 -16.6989   4.0114  0.4045
 [92,]  91 -16.9457 -16.2454 -16.6944   2.1719  0.7042
 [93,]  92 -16.9444 -16.2364 -16.6903   3.5563  0.4694
 [94,]  93 -16.9428 -16.2271 -16.6860   2.8128  0.5896
 [95,]  94 -16.9404 -16.2170 -16.6808   0.3291  0.9879
 [96,]  95 -16.9395 -16.2085 -16.6772   4.8324  0.3049
 [97,]  96 -16.9398 -16.2011 -16.6748   8.1648  0.0857
 [98,]  97 -16.9385 -16.1920 -16.6706   3.3971  0.4937
 [99,]  98 -16.9379 -16.1838 -16.6673   5.6507  0.2268
[100,]  99 -16.9371 -16.1753 -16.6638   4.9952  0.2878
[101,] 100 -16.9394 -16.1699 -16.6633  13.7237  0.0082
> detach("package:MTS", unload = TRUE)
> var.m1<- VAR(diff_log_dataframe[, c("diff_log_Fut","diff_log_USL")], p = 6, type = "const")
> # display the results for each equation
> summary(var.m1)

VAR Estimation Results:
========================= 
Endogenous variables: diff_log_Fut, diff_log_USL 
Deterministic variables: const 
Sample size: 3141 
Log Likelihood: 17509.433 
Roots of the characteristic polynomial:
0.7642 0.7642 0.7006 0.7006 0.591 0.591 0.5634 0.5634 0.5558 0.5558 0.5512 0.5302
Call:
VAR(y = diff_log_dataframe[, c("diff_log_Fut", "diff_log_USL")], 
    p = 6, type = "const")


Estimation results for equation diff_log_Fut: 
============================================= 
diff_log_Fut = diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + diff_log_USL.l6 + const 

                  Estimate Std. Error t value Pr(>|t|)    
diff_log_Fut.l1 -0.1659467  0.0339424  -4.889 1.06e-06 ***
diff_log_USL.l1  0.1673456  0.0436141   3.837 0.000127 ***
diff_log_Fut.l2 -0.0857581  0.0344167  -2.492 0.012763 *  
diff_log_USL.l2  0.0802019  0.0442190   1.814 0.069813 .  
diff_log_Fut.l3  0.1346629  0.0346410   3.887 0.000103 ***
diff_log_USL.l3 -0.1591189  0.0443706  -3.586 0.000341 ***
diff_log_Fut.l4  0.1272179  0.0347195   3.664 0.000252 ***
diff_log_USL.l4 -0.1414052  0.0444419  -3.182 0.001478 ** 
diff_log_Fut.l5 -0.0857248  0.0344829  -2.486 0.012970 *  
diff_log_USL.l5  0.0898457  0.0442255   2.032 0.042285 *  
diff_log_Fut.l6 -0.0585723  0.0339752  -1.724 0.084812 .  
diff_log_USL.l6  0.0800061  0.0436148   1.834 0.066693 .  
const           -0.0001114  0.0004182  -0.266 0.789902    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.02341 on 3128 degrees of freedom
Multiple R-Squared: 0.02475,	Adjusted R-squared: 0.02101 
F-statistic: 6.616 on 12 and 3128 DF,  p-value: 7.879e-12 


Estimation results for equation diff_log_USL: 
============================================= 
diff_log_USL = diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + diff_log_USL.l6 + const 

                  Estimate Std. Error t value Pr(>|t|)    
diff_log_Fut.l1  0.0425544  0.0264010   1.612 0.107096    
diff_log_USL.l1 -0.0881794  0.0339240  -2.599 0.009385 ** 
diff_log_Fut.l2  0.0364277  0.0267700   1.361 0.173686    
diff_log_USL.l2 -0.0396952  0.0343944  -1.154 0.248541    
diff_log_Fut.l3  0.1121559  0.0269445   4.162 3.23e-05 ***
diff_log_USL.l3 -0.1239459  0.0345124  -3.591 0.000334 ***
diff_log_Fut.l4  0.0857349  0.0270055   3.175 0.001514 ** 
diff_log_USL.l4 -0.0743118  0.0345678  -2.150 0.031652 *  
diff_log_Fut.l5  0.0409723  0.0268215   1.528 0.126715    
diff_log_USL.l5 -0.0427415  0.0343995  -1.243 0.214144    
diff_log_Fut.l6  0.0447080  0.0264266   1.692 0.090787 .  
diff_log_USL.l6 -0.0120778  0.0339245  -0.356 0.721849    
const           -0.0003127  0.0003253  -0.961 0.336587    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.01821 on 3128 degrees of freedom
Multiple R-Squared: 0.01117,	Adjusted R-squared: 0.007379 
F-statistic: 2.945 on 12 and 3128 DF,  p-value: 0.0004354 



Covariance matrix of residuals:
             diff_log_Fut diff_log_USL
diff_log_Fut    0.0005479    0.0003631
diff_log_USL    0.0003631    0.0003315

Correlation matrix of residuals:
             diff_log_Fut diff_log_USL
diff_log_Fut       1.0000       0.8521
diff_log_USL       0.8521       1.0000


> vars <- c("diff_log_Fut", "diff_log_USL")
> names <- c("diff_log_Fut.l", "diff_log_USL.l")
> # Loop through each variable and run the linear hypothesis test
> for (n in names) {
+   for (var in vars) {
+     print(paste("Testing for:", var))
+     print(paste("Hypothesis matrix used:", n))
+     
+     # Construct the hypothesis matrix
+     hypothesis_matrix <- paste0(n, 1:6)  # Assuming 7 lags for the VAR model
+     
+     # Run the linear hypothesis test
+     result <- linearHypothesis(
+       model = var.m1$varresult[[var]], 
+       hypothesis.matrix = hypothesis_matrix, 
+       rhs = rep(0, 6)
+     )
+     
+     # Print the result of the hypothesis test
+     print(result)  # Assuming result[6] contains the test statistic or p-value
+   }
+ } 
[1] "Testing for: diff_log_Fut"
[1] "Hypothesis matrix used: diff_log_Fut.l"

Linear hypothesis test:
diff_log_Fut.l1 = 0
diff_log_Fut.l2 = 0
diff_log_Fut.l3 = 0
diff_log_Fut.l4 = 0
diff_log_Fut.l5 = 0
diff_log_Fut.l6 = 0

Model 1: restricted model
Model 2: y ~ -1 + (diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + 
    diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + 
    diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + 
    diff_log_USL.l6 + const)

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   3134 1.7564                                  
2   3128 1.7138  6  0.042675 12.982 1.496e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Testing for: diff_log_USL"
[1] "Hypothesis matrix used: diff_log_Fut.l"

Linear hypothesis test:
diff_log_Fut.l1 = 0
diff_log_Fut.l2 = 0
diff_log_Fut.l3 = 0
diff_log_Fut.l4 = 0
diff_log_Fut.l5 = 0
diff_log_Fut.l6 = 0

Model 1: restricted model
Model 2: y ~ -1 + (diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + 
    diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + 
    diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + 
    diff_log_USL.l6 + const)

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   3134 1.0454                                  
2   3128 1.0368  6 0.0086045 4.3265 0.0002356 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Testing for: diff_log_Fut"
[1] "Hypothesis matrix used: diff_log_USL.l"

Linear hypothesis test:
diff_log_USL.l1 = 0
diff_log_USL.l2 = 0
diff_log_USL.l3 = 0
diff_log_USL.l4 = 0
diff_log_USL.l5 = 0
diff_log_USL.l6 = 0

Model 1: restricted model
Model 2: y ~ -1 + (diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + 
    diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + 
    diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + 
    diff_log_USL.l6 + const)

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   3134 1.7422                                  
2   3128 1.7138  6  0.028437 8.6506 2.339e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Testing for: diff_log_USL"
[1] "Hypothesis matrix used: diff_log_USL.l"

Linear hypothesis test:
diff_log_USL.l1 = 0
diff_log_USL.l2 = 0
diff_log_USL.l3 = 0
diff_log_USL.l4 = 0
diff_log_USL.l5 = 0
diff_log_USL.l6 = 0

Model 1: restricted model
Model 2: y ~ -1 + (diff_log_Fut.l1 + diff_log_USL.l1 + diff_log_Fut.l2 + 
    diff_log_USL.l2 + diff_log_Fut.l3 + diff_log_USL.l3 + diff_log_Fut.l4 + 
    diff_log_USL.l4 + diff_log_Fut.l5 + diff_log_USL.l5 + diff_log_Fut.l6 + 
    diff_log_USL.l6 + const)

  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1   3134 1.0438                                
2   3128 1.0368  6 0.0069364 3.4877 0.001936 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> serial.test(var.m1, lags.pt=20, type="PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object var.m1
Chi-squared = 213.34, df = 56, p-value < 2.2e-16

> roots(var.m1)
 [1] 0.7641953 0.7641953 0.7006276 0.7006276 0.5910127 0.5910127 0.5633950
 [8] 0.5633950 0.5558323 0.5558323 0.5512034 0.5301539
> N = var.m1$obs
> nreg = var.m1$K * (var.m1$p * var.m1$K + 1)
> res = residuals(var.m1)
> det.cov.m1 = det(crossprod(res) / N) # "/N" is applied in RATS
> aic.var.m1 = N * log(det.cov.m1) + 2 * nreg
> sbc.var.m1 = N * log(det.cov.m1) + log(N) * nreg
> print(paste("AIC.var.m1 =", sprintf("%.3f", aic.var.m1), 
+             "SBC.var.m1 =", sprintf("%.3f", sbc.var.m1)))
[1] "AIC.var.m1 = -52794.410 SBC.var.m1 = -52637.050"
> fevd.m1 <- fevd(var.m1, n.ahead = 30)
> fevd.m1   
$diff_log_Fut
      diff_log_Fut diff_log_USL
 [1,]    1.0000000  0.000000000
 [2,]    0.9953946  0.004605381
 [3,]    0.9951656  0.004834400
 [4,]    0.9892554  0.010744612
 [5,]    0.9877543  0.012245737
 [6,]    0.9842010  0.015798969
 [7,]    0.9840035  0.015996545
 [8,]    0.9834838  0.016516178
 [9,]    0.9833767  0.016623260
[10,]    0.9831760  0.016823971
[11,]    0.9830962  0.016903799
[12,]    0.9830460  0.016954029
[13,]    0.9830117  0.016988338
[14,]    0.9830105  0.016989522
[15,]    0.9829907  0.017009308
[16,]    0.9829907  0.017009319
[17,]    0.9829841  0.017015909
[18,]    0.9829838  0.017016249
[19,]    0.9829825  0.017017488
[20,]    0.9829818  0.017018159
[21,]    0.9829816  0.017018424
[22,]    0.9829812  0.017018782
[23,]    0.9829812  0.017018807
[24,]    0.9829811  0.017018937
[25,]    0.9829811  0.017018938
[26,]    0.9829810  0.017018986
[27,]    0.9829810  0.017018991
[28,]    0.9829810  0.017019003
[29,]    0.9829810  0.017019008
[30,]    0.9829810  0.017019009

$diff_log_USL
      diff_log_Fut diff_log_USL
 [1,]    0.7260972    0.2739028
 [2,]    0.7246687    0.2753313
 [3,]    0.7245852    0.2754148
 [4,]    0.7224995    0.2775005
 [5,]    0.7225915    0.2774085
 [6,]    0.7225376    0.2774624
 [7,]    0.7229777    0.2770223
 [8,]    0.7229960    0.2770040
 [9,]    0.7229848    0.2770152
[10,]    0.7229881    0.2770119
[11,]    0.7229883    0.2770117
[12,]    0.7229848    0.2770152
[13,]    0.7229843    0.2770157
[14,]    0.7229836    0.2770164
[15,]    0.7229833    0.2770167
[16,]    0.7229834    0.2770166
[17,]    0.7229833    0.2770167
[18,]    0.7229833    0.2770167
[19,]    0.7229832    0.2770168
[20,]    0.7229832    0.2770168
[21,]    0.7229832    0.2770168
[22,]    0.7229832    0.2770168
[23,]    0.7229832    0.2770168
[24,]    0.7229832    0.2770168
[25,]    0.7229832    0.2770168
[26,]    0.7229832    0.2770168
[27,]    0.7229832    0.2770168
[28,]    0.7229832    0.2770168
[29,]    0.7229832    0.2770168
[30,]    0.7229832    0.2770168

> plot(fevd.m1)
> #Swap order of VAR
> var.m2<- VAR(diff_log_dataframe[, c("diff_log_USL","diff_log_Fut")], p = 6, type = "const")
> fevd.m2 <- fevd(var.m2, n.ahead = 30)
> fevd.m2
$diff_log_USL
      diff_log_USL diff_log_Fut
 [1,]    1.0000000  0.000000000
 [2,]    0.9991822  0.000817754
 [3,]    0.9988863  0.001113672
 [4,]    0.9943982  0.005601758
 [5,]    0.9929427  0.007057270
 [6,]    0.9928333  0.007166672
 [7,]    0.9924442  0.007555826
 [8,]    0.9923989  0.007601128
 [9,]    0.9923306  0.007669447
[10,]    0.9922940  0.007705995
[11,]    0.9922937  0.007706331
[12,]    0.9922852  0.007714785
[13,]    0.9922850  0.007715037
[14,]    0.9922825  0.007717509
[15,]    0.9922812  0.007718849
[16,]    0.9922811  0.007718865
[17,]    0.9922810  0.007718974
[18,]    0.9922808  0.007719154
[19,]    0.9922808  0.007719184
[20,]    0.9922807  0.007719272
[21,]    0.9922807  0.007719273
[22,]    0.9922807  0.007719273
[23,]    0.9922807  0.007719284
[24,]    0.9922807  0.007719284
[25,]    0.9922807  0.007719289
[26,]    0.9922807  0.007719289
[27,]    0.9922807  0.007719289
[28,]    0.9922807  0.007719289
[29,]    0.9922807  0.007719290
[30,]    0.9922807  0.007719290

$diff_log_Fut
      diff_log_USL diff_log_Fut
 [1,]    0.7260972    0.2739028
 [2,]    0.7206964    0.2793036
 [3,]    0.7202415    0.2797585
 [4,]    0.7149681    0.2850319
 [5,]    0.7132798    0.2867202
 [6,]    0.7097518    0.2902482
 [7,]    0.7096735    0.2903265
 [8,]    0.7091628    0.2908372
 [9,]    0.7090623    0.2909377
[10,]    0.7088264    0.2911736
[11,]    0.7087587    0.2912413
[12,]    0.7086871    0.2913129
[13,]    0.7086611    0.2913389
[14,]    0.7086587    0.2913413
[15,]    0.7086381    0.2913619
[16,]    0.7086382    0.2913618
[17,]    0.7086306    0.2913694
[18,]    0.7086305    0.2913695
[19,]    0.7086290    0.2913710
[20,]    0.7086284    0.2913716
[21,]    0.7086281    0.2913719
[22,]    0.7086277    0.2913723
[23,]    0.7086277    0.2913723
[24,]    0.7086276    0.2913724
[25,]    0.7086276    0.2913724
[26,]    0.7086275    0.2913725
[27,]    0.7086275    0.2913725
[28,]    0.7086275    0.2913725
[29,]    0.7086275    0.2913725
[30,]    0.7086275    0.2913725

> irf.m1 <- irf(var.m1, impulse = NULL, response = NULL, n.ahead = 30, 
+               ortho = TRUE, boot = TRUE, runs = 100)
> irf.m1$irf$diff_log_Fut
       diff_log_Fut  diff_log_USL
 [1,]  2.340682e-02  1.551383e-02
 [2,] -1.288114e-03 -3.719370e-04
 [3,] -6.115692e-04  2.148154e-04
 [4,]  9.015601e-04  6.252130e-04
 [5,]  6.944432e-04  7.079794e-04
 [6,] -8.644448e-04  9.307967e-05
 [7,]  1.721780e-05  7.296452e-04
 [8,]  3.254104e-04 -1.651580e-04
 [9,] -1.356885e-04 -1.212078e-04
[10,] -2.745433e-04 -1.152630e-04
[11,]  9.573615e-05 -1.432873e-05
[12,]  1.839405e-04 -3.036338e-05
[13,] -4.908459e-05  7.778953e-06
[14,] -4.375014e-05  1.883627e-05
[15,]  7.048200e-05  1.712201e-05
[16,]  1.364260e-05  4.358071e-06
[17,] -4.768913e-05  3.343718e-06
[18,] -2.105890e-06  4.878584e-06
[19,]  2.215310e-05 -1.515468e-06
[20,] -9.626805e-06 -4.138786e-06
[21,] -1.225239e-05 -1.038523e-06
[22,]  8.557822e-06 -2.640744e-08
[23,]  5.617265e-06 -1.060679e-06
[24,] -5.491403e-06 -2.985269e-07
[25,] -9.261189e-07  8.203589e-07
[26,]  3.945376e-06  3.591763e-07
[27,] -3.782125e-07 -1.320297e-07
[28,] -2.260713e-06  1.292497e-07
[29,]  6.890068e-07  1.931904e-07
[30,]  9.942990e-07 -1.045998e-07
[31,] -7.666123e-07 -1.282351e-07
> irf.m1$irf$diff_log_USL
       diff_log_Fut  diff_log_USL
 [1,]  0.000000e+00  9.528396e-03
 [2,]  1.594535e-03 -8.402086e-04
 [3,]  3.589824e-04 -2.362876e-04
 [4,] -1.819392e-03 -1.053456e-03
 [5,] -9.230522e-04 -3.871693e-04
 [6,]  1.423613e-03 -1.681628e-04
 [7,]  3.355733e-04  1.816704e-05
 [8,] -5.460310e-04  4.421748e-05
 [9,]  2.478630e-04  1.038960e-04
[10,]  3.404401e-04  5.978565e-05
[11,] -2.139160e-04  3.743966e-06
[12,] -1.711198e-04  4.406772e-05
[13,]  1.401688e-04  1.571447e-05
[14,]  2.663875e-05 -2.234601e-05
[15,] -1.067375e-04 -1.446620e-05
[16,]  3.106532e-06 -8.124814e-09
[17,]  6.168690e-05 -5.077740e-06
[18,] -1.393201e-05 -6.158776e-06
[19,] -2.677381e-05  2.795606e-06
[20,]  1.962517e-05  3.871846e-06
[21,]  1.240847e-05 -1.953779e-07
[22,] -1.434304e-05 -9.619717e-08
[23,] -3.868848e-06  1.587509e-06
[24,]  8.650868e-06  2.224901e-07
[25,] -7.918335e-07 -9.659804e-07
[26,] -5.259957e-06 -1.480540e-07
[27,]  1.729897e-06  2.652866e-07
[28,]  2.619388e-06 -2.458957e-07
[29,] -1.619795e-06 -2.292594e-07
[30,] -9.792148e-07  1.705630e-07
[31,]  1.292688e-06  1.133225e-07

> plot.data <- rbindlist(
+   lapply(
+     irf.m1[c("irf", "Lower", "Upper")],
+     function(dt) {
+       temp <- rbindlist(
+         lapply(dt, function(x) as.data.table(x)[, Step := 0:(.N-1)]), idcol = "Shock"
+       )
+       temp <- melt(temp, id.vars = c("Shock", "Step"), variable.name = "Response")
+     }
+   ),
+   idcol = "data.type"
+ )
> # Reorder 'Shock' and 'Response' factors in the desired order
> plot.data[, Shock := factor(Shock, levels = c("diff_log_Fut", "diff_log_USL"))]
> plot.data[, Response := factor(Response, levels = c("diff_log_Fut", "diff_log_USL"))]
> # Create the plot with reordered factors
> ggplot(data = plot.data,
+        aes(x = Step, y = value, linetype = data.type)) + 
+   geom_line() + 
+   facet_grid(rows = vars(Response), cols = vars(Shock), scales = "free_y") +
+   geom_hline(yintercept = 0, colour = "red") + 
+   theme_bw() + 
+   theme(legend.position = "none") +
+   ylab("") + 
+   ggtitle("Responses (in rows) to Shocks (in columns)")
> res <- residuals(var.m1)
> str(res)
 num [1:3141, 1:2] -0.010675 -0.000145 0.01327 -0.002603 0.021582 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:3141] "1" "2" "3" "4" ...
  ..$ : chr [1:2] "diff_log_Fut" "diff_log_USL"
> # Plot all residuals on the same graph
> matplot(res, type = "l", col = 1:5, lty = 1, xlab = "Index", ylab = "Residuals", main = "Residuals of VAR Model")
> # Add a legend
> legend("topright", legend = colnames(res), col = 1:5, lty = 1)
> garch11 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
+                       mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
+                       distribution.model = "norm")
> garch.fit.FUT <- ugarchfit(res[, 1], spec = garch11)
> show(garch.fit.FUT)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,1)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000000    0.000313   0.0000  1.00000
omega   0.000001    0.000000   1.5595  0.11889
alpha1  0.069975    0.007074   9.8924  0.00000
beta1   0.897534    0.010874  82.5392  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000000    0.000548  0.00000  1.00000
omega   0.000001    0.000004  0.16143  0.87175
alpha1  0.069975    0.051558  1.35721  0.17472
beta1   0.897534    0.085996 10.43690  0.00000

LogLikelihood : 7580.285 

Information Criteria
------------------------------------
                    
Akaike       -4.8241
Bayes        -4.8164
Shibata      -4.8241
Hannan-Quinn -4.8214

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.9667  0.3255
Lag[2*(p+q)+(p+q)-1][2]    1.8165  0.2951
Lag[4*(p+q)+(p+q)-1][5]    2.7332  0.4579
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.1065  0.7441
Lag[2*(p+q)+(p+q)-1][5]    1.5464  0.7285
Lag[4*(p+q)+(p+q)-1][9]    4.5521  0.4980
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.2674 0.500 2.000  0.6050
ARCH Lag[5]    2.0053 1.440 1.667  0.4701
ARCH Lag[7]    4.5410 2.315 1.543  0.2752

Nyblom stability test
------------------------------------
Joint Statistic:  56.4381
Individual Statistics:              
mu      0.1526
omega  12.3644
alpha1 31.1555
beta1  29.4955

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.07 1.24 1.6
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value    prob sig
Sign Bias           0.9969 0.31889    
Negative Sign Bias  1.5970 0.11037    
Positive Sign Bias  2.0701 0.03853  **
Joint Effect        9.3448 0.02504  **


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     237.2    1.180e-39
2    30     321.2    5.091e-51
3    40     358.9    2.342e-53
4    50     394.7    1.504e-55


Elapsed time : 0.1627362 

> # Align the data by removing the first few rows
> aligned_diff_log_dataframe <- diff_log_dataframe[-c(1:(nrow(diff_log_dataframe) - length((coredata(sigma(garch.fit.FUT)))^2))), ]
> # Add h_t column
> aligned_diff_log_dataframe$f_t <- (coredata(sigma(garch.fit.FUT)))^2
> par(mfrow = c(2,1))
> plot(aligned_diff_log_dataframe$Dates, aligned_diff_log_dataframe$f_t, type = "l",
+      main = "GARCH(1,1) conditional variance of FUT", ylab = "")
> garch.fit.USL <- ugarchfit(res[, 2], spec = garch11)
> show(garch.fit.USL)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,1)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000000    0.000258   0.0000  1.00000
omega   0.000000    0.000000   1.5258  0.12705
alpha1  0.065192    0.004421  14.7444  0.00000
beta1   0.904753    0.006224 145.3615  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000000    0.000424  0.00000  1.00000
omega   0.000000    0.000002  0.16177  0.87149
alpha1  0.065192    0.042883  1.52023  0.12845
beta1   0.904753    0.072769 12.43322  0.00000

LogLikelihood : 8212.853 

Information Criteria
------------------------------------
                    
Akaike       -5.2269
Bayes        -5.2192
Shibata      -5.2269
Hannan-Quinn -5.2241

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      2.129  0.1445
Lag[2*(p+q)+(p+q)-1][2]     3.033  0.1375
Lag[4*(p+q)+(p+q)-1][5]     3.627  0.3044
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4011  0.5265
Lag[2*(p+q)+(p+q)-1][5]    2.1664  0.5796
Lag[4*(p+q)+(p+q)-1][9]    4.8798  0.4473
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]     0.632 0.500 2.000  0.4266
ARCH Lag[5]     3.355 1.440 1.667  0.2423
ARCH Lag[7]     5.329 2.315 1.543  0.1933

Nyblom stability test
------------------------------------
Joint Statistic:  105.1132
Individual Statistics:               
mu      0.03632
omega   6.77411
alpha1 32.11587
beta1  32.58310

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.07 1.24 1.6
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value     prob sig
Sign Bias            1.797 0.072425   *
Negative Sign Bias   2.228 0.025928  **
Positive Sign Bias   1.641 0.100821    
Joint Effect        11.935 0.007609 ***


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     242.0    1.259e-40
2    30     301.8    3.648e-47
3    40     334.3    1.397e-48
4    50     382.4    3.493e-53


Elapsed time : 0.06641102 

> aligned_diff_log_dataframe <- diff_log_dataframe[-c(1:(nrow(diff_log_dataframe) - length((coredata(sigma(garch.fit.FUT)))^2))), ]
> aligned_diff_log_dataframe$U_t <- (coredata(sigma(garch.fit.USL)))^2
> #par(mfrow = c(1,1))
> plot(aligned_diff_log_dataframe$Dates, aligned_diff_log_dataframe$U_t, type = "l",
+      main = "GARCH(1,1) conditional variance of USL", ylab = "")
> # standardized residuals -----
> stdres_Fut <- coredata(residuals(garch.fit.FUT, standardize = TRUE))
> # coredata --> extracts the core data contained in a more complex object and replace it
> plot(stdres_Fut)
> # standardized residuals -----
> stdres_USL <- coredata(residuals(garch.fit.USL, standardize = TRUE))
> # coredata --> extracts the core data contained in a more complex object and replace it
> plot(stdres_USL)
> Box.test(stdres_Fut, lag = 24, type = "Ljung-Box",fitdf =length(coef(garch.fit.FUT)))

	Box-Ljung test

data:  stdres_Fut
X-squared = 14.18, df = 20, p-value = 0.8213

> Box.test(stdres_USL, lag = 24, type = "Ljung-Box",fitdf =length(coef(garch.fit.USL)))

	Box-Ljung test

data:  stdres_USL
X-squared = 15.493, df = 20, p-value = 0.7475