> # Assuming log_data is already loaded in the environment
> print("Loaded log_data and created Y matrix:")
[1] "Loaded log_data and created Y matrix:"
> Y <- log_data[, c("log_Spot", "log_Fut", "log_USO", "log_USL", "log_OIL")]
> print(head(Y))
# A tibble: 6 × 5
  log_Spot log_Fut log_USO log_USL log_OIL
     <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
1     4.12    4.12    4.64    3.14    5.64
2     4.12    4.12    4.64    3.15    5.64
3     4.12    4.12    4.64    3.14    5.64
4     4.11    4.11    4.63    3.14    5.63
5     4.10    4.10    4.62    3.13    5.62
6     4.10    4.10    4.62    3.13    5.61
> 
> # Step 2: Determine the optimal lag length using VARselect
> print("Determining optimal lag length:")
[1] "Determining optimal lag length:"
> VARselect_result <- VARselect(Y, lag.max=15, type="const")
> print(VARselect_result)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
    15      9      5     15 

$criteria
                   1             2             3             4             5
AIC(n) -4.711556e+01 -4.740194e+01 -4.747100e+01 -4.753801e+01 -4.759086e+01
HQ(n)  -4.709405e+01 -4.736252e+01 -4.741366e+01 -4.746275e+01 -4.749768e+01
SC(n)  -4.705576e+01 -4.729231e+01 -4.731154e+01 -4.732872e+01 -4.733173e+01
FPE(n)  3.451226e-21  2.591776e-21  2.418841e-21  2.262069e-21  2.145624e-21
                   6             7             8             9            10
AIC(n) -4.762770e+01 -4.765247e+01 -4.769426e+01 -4.771974e+01 -4.772879e+01
HQ(n)  -4.751659e+01 -4.752345e+01 -4.754732e+01 -4.755488e+01 -4.754601e+01
SC(n)  -4.731874e+01 -4.729368e+01 -4.728564e+01 -4.726129e+01 -4.722051e+01
FPE(n)  2.068028e-21  2.017422e-21  1.934861e-21  1.886189e-21  1.869205e-21
                  11            12            13            14            15
AIC(n) -4.774500e+01 -4.776265e+01 -4.777702e+01 -4.780304e+01 -4.781737e+01
HQ(n)  -4.754430e+01 -4.754403e+01 -4.754048e+01 -4.754859e+01 -4.754499e+01
SC(n)  -4.718688e+01 -4.715470e+01 -4.711924e+01 -4.709543e+01 -4.705993e+01
FPE(n)  1.839162e-21  1.806993e-21  1.781230e-21  1.735487e-21  1.710822e-21

> 
> p <- VARselect_result$selection["HQ(n)"]
> print(paste("Optimal lag length (based on HQ criterion):", p))
[1] "Optimal lag length (based on HQ criterion): 9"
> 
> K <- p - 1  # Lagged differences in Johansen test
> print(paste("Using K =", K, "for lagged differences in Johansen test"))
[1] "Using K = 8 for lagged differences in Johansen test"
> 
> # Step 3: Johansen cointegration test with adjusted settings
> print("Performing Johansen cointegration test:")
[1] "Performing Johansen cointegration test:"
> johansen_test <- ca.jo(Y, type="trace", ecdet="const", K=K, spec="longrun")
> print(summary(johansen_test))

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

Test type: trace statistic , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 9.194013e-02 6.137328e-02 1.082481e-02 5.234867e-03 7.584923e-04
[6] 9.339895e-17

Values of teststatistic and critical values of test:

           test 10pct  5pct  1pct
r <= 4 |   2.29  7.52  9.24 12.97
r <= 3 |  18.16 17.85 19.96 24.60
r <= 2 |  51.06 32.00 34.91 41.07
r <= 1 | 242.53 49.65 53.12 60.16
r = 0  | 534.08 71.86 76.07 84.45

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

             log_Spot.l8 log_Fut.l8 log_USO.l8 log_USL.l8 log_OIL.l8
log_Spot.l8  1.000000000  1.0000000   1.000000  1.0000000   1.000000
log_Fut.l8  -1.003644123 -1.5256307  -3.331999  1.4487665  -5.903652
log_USO.l8  -0.006998315 -0.9992708   8.252360 13.6547508   8.825436
log_USL.l8  -0.002521254  0.7607663   4.716197 -8.9283064  -1.004589
log_OIL.l8   0.008398995  0.6448770  -8.722794 -8.0800039  -5.309850
constant     0.008065544  0.7922468   5.954291  0.1266867  10.436638
             constant
log_Spot.l8  1.000000
log_Fut.l8  -1.197880
log_USO.l8  -2.604882
log_USL.l8   1.305870
log_OIL.l8   1.247534
constant     2.071152

Weights W:
(This is the loading matrix)

            log_Spot.l8  log_Fut.l8   log_USO.l8    log_USL.l8   log_OIL.l8
log_Spot.d -0.061910397 0.004041937 -0.003984373 -0.0027565091 0.0003042509
log_Fut.d   0.338164478 0.016601908 -0.003969795 -0.0024759995 0.0002939173
log_USO.d  -0.022184653 0.037273171 -0.004320462 -0.0016669625 0.0002652621
log_USL.d   0.001740984 0.018450168 -0.005337932 -0.0004077368 0.0002501097
log_OIL.d  -0.060171749 0.039992249 -0.002479446 -0.0012813480 0.0003394962
               constant
log_Spot.d 6.558352e-15
log_Fut.d  2.942492e-14
log_USO.d  1.311205e-14
log_USL.d  1.067847e-14
log_OIL.d  8.429320e-15

> 
> r <- 4  # Correctly set as number of cointegrating vectors
> print(paste("Number of cointegrating vectors set to:", r))
[1] "Number of cointegrating vectors set to: 4"
> 
> # Define the restriction matrix R and vector r_vec
> n_vars <- ncol(Y)
> print(paste("Number of variables:", n_vars))
[1] "Number of variables: 5"
> 
> n_beta_coeffs <- n_vars * r
> print(paste("Total number of beta coefficients:", n_beta_coeffs))
[1] "Total number of beta coefficients: 20"
> 
> R <- matrix(0, nrow=3, ncol=n_beta_coeffs)
> r_vec <- rep(0, nrow(R))  # Replicates zeros to match the number of rows in R
> print("Restriction matrix R initialized:")
[1] "Restriction matrix R initialized:"
> print(R)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,]    0    0    0    0    0    0    0    0    0     0     0     0     0
[2,]    0    0    0    0    0    0    0    0    0     0     0     0     0
[3,]    0    0    0    0    0    0    0    0    0     0     0     0     0
     [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,]     0     0     0     0     0     0     0
[2,]     0     0     0     0     0     0     0
[3,]     0     0     0     0     0     0     0
> 
> # Imposing the restrictions based on the paper
> print("Imposing restrictions based on the paper:")
[1] "Imposing restrictions based on the paper:"
> # Restriction 1
> R[1, (1 - 1) * n_vars + 3] <- 1  # beta[3,1]
> print("Restriction 1 applied on beta[3,1]")
[1] "Restriction 1 applied on beta[3,1]"
> 
> # Restriction 2
> R[2, (1 - 1) * n_vars + 4] <- 1  # beta[4,1]
> print("Restriction 2 applied on beta[4,1]")
[1] "Restriction 2 applied on beta[4,1]"
> 
> # Restriction 3
> R[3, (2 - 1) * n_vars + 5] <- 1  # beta[5,2]
> print("Restriction 3 applied on beta[5,2]")
[1] "Restriction 3 applied on beta[5,2]"
> 
> print("Complete Restriction Matrix R:")
[1] "Complete Restriction Matrix R:"
> print(R)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,]    0    0    1    0    0    0    0    0    0     0     0     0     0
[2,]    0    0    0    1    0    0    0    0    0     0     0     0     0
[3,]    0    0    0    0    0    0    0    0    0     1     0     0     0
     [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,]     0     0     0     0     0     0     0
[2,]     0     0     0     0     0     0     0
[3,]     0     0     0     0     0     0     0
> 
> # Test the restrictions
> print("Testing the restrictions with blrtest:")
[1] "Testing the restrictions with blrtest:"
> blr_test <- blrtest(johansen_test, H=R, r=r)
Error in blrtest(johansen_test, H = R, r = r) : 
Row number of 'H' is unequal to VAR order.