> # 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_USL", "log_OIL")]
> # 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_USL : num [1:3031] 3.14 3.15 3.14 3.14 3.13 ...
 $ log_OIL : num [1:3031] 5.64 5.64 5.64 5.63 5.62 ...
 - 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_USL      -9.306146 0.5907100
Dickey-Fuller Z(alpha)4  log_OIL      -7.929339 0.6674934
> 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_USL       20.71650    0.01
KPSS Level4  log_OIL       24.82409    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_USL      -3208.247    0.01
Dickey-Fuller Z(alpha)4  log_OIL      -3212.187    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_USL     0.07285208     0.1
KPSS Level4  log_OIL     0.10775410     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 AIC
> 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="trace", ecdet="const", K=p, spec="transitory")
> summary(johansen_test)

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

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

Eigenvalues (lambda):
[1] 1.376445e-01 5.485084e-02 1.460022e-02 6.534633e-03 7.806915e-04 2.873217e-17

Values of teststatistic and critical values of test:

           test 10pct  5pct  1pct
r <= 4 |   2.36  7.52  9.24 12.97
r <= 3 |  22.20 17.85 19.96 24.60
r <= 2 |  66.71 32.00 34.91 41.07
r <= 1 | 237.41 49.65 53.12 60.16
r = 0  | 685.53 71.86 76.07 84.45

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

              log_Spot.l1 log_Fut.l1 log_USO.l1 log_USL.l1 log_OIL.l1  constant
log_Spot.l1  1.0000000000  1.0000000    1.00000  1.0000000   1.000000  1.000000
log_Fut.l1  -1.0017520247 -1.3244152   16.80517  0.7523886  -8.329693 -1.820012
log_USO.l1   0.0023907033 -0.5014035  -74.82706  8.6459957  15.122284 -9.148170
log_USL.l1  -0.0043867524  0.5028687  -35.29678 -6.3180603  -4.415480  5.683194
log_OIL.l1   0.0005538807  0.2823013   76.35024 -4.8355656  -8.218566  3.311713
constant     0.0069055220  0.5034032  -47.90226 -0.2647043  16.978031 10.890075

Weights W:
(This is the loading matrix)

           log_Spot.l1 log_Fut.l1   log_USO.l1    log_USL.l1   log_OIL.l1      constant
log_Spot.d -0.24199295 0.02295218 0.0003860446 -0.0045865854 0.0001805646 -3.849043e-15
log_Fut.d   0.18580099 0.04212123 0.0004176988 -0.0042320772 0.0001764496  1.940801e-15
log_USO.d  -0.15331674 0.06760043 0.0005287744 -0.0024737038 0.0001574821 -2.212215e-15
log_USL.d  -0.07075644 0.02880841 0.0006447661 -0.0007056863 0.0001551454 -1.849213e-15
log_OIL.d  -0.18973739 0.07470091 0.0002518853 -0.0017978899 0.0001975052 -1.958137e-15

> 
> # 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)
> print(vecm_var)

Coefficient matrix of lagged endogenous variables:

A1:
         log_Spot.l1  log_Fut.l1  log_USO.l1 log_USL.l1 log_OIL.l1
log_Spot   1.0314300 -0.11368162 -0.44715044  0.1682500  0.3250480
log_Fut    0.2645245  0.59622236 -0.42963580  0.1877418  0.3369964
log_USO    0.2434639  0.05889863  0.03556935  0.3508134  0.2977061
log_USL    0.2032903 -0.02635639 -0.42589721  1.0087496  0.1832933
log_OIL    0.2104761  0.08089268 -0.91049287  0.3401844  1.2535131


A2:
         log_Spot.l2  log_Fut.l2 log_USO.l2  log_USL.l2 log_OIL.l2
log_Spot -0.25695033  0.17157095  0.3068445  0.04945633 -0.1870758
log_Fut  -0.04085083 -0.04446826  0.3808412 -0.03100964 -0.1976518
log_USO  -0.40256986  0.06270538  0.6244124 -0.07113705 -0.1625828
log_USL  -0.34630385  0.08910242  0.4026364  0.05302569 -0.1230994
log_OIL  -0.37715993  0.04294832  0.7207870 -0.06653196 -0.2329700


A3:
         log_Spot.l3 log_Fut.l3 log_USO.l3 log_USL.l3   log_OIL.l3
log_Spot   0.1458521 -0.1510641  0.3339362 -0.3276682 -0.027737351
log_Fut    0.2034098 -0.1089163  0.1496027 -0.2465002 -0.015756806
log_USO    0.1856833 -0.2871774  0.2521153 -0.1912245  0.005717099
log_USL    0.1226135 -0.2350356  0.1789151 -0.1471299  0.038250937
log_OIL    0.2326656 -0.3121597  0.2278128 -0.1744500 -0.025562538


A4:
         log_Spot.l4 log_Fut.l4 log_USO.l4  log_USL.l4 log_OIL.l4
log_Spot -0.18863718  0.1710641 0.19081976  0.01677350 -0.1791175
log_Fut  -0.25618937  0.1942828 0.20075854  0.09335346 -0.1992789
log_USO  -0.18182951  0.2347955 0.16972329 -0.03358236 -0.1772673
log_USL  -0.04536394  0.1022300 0.09913327  0.04443748 -0.1707412
log_OIL  -0.31741773  0.3324992 0.18859322 -0.04652779 -0.1442159


A5:
         log_Spot.l5   log_Fut.l5 log_USO.l5  log_USL.l5 log_OIL.l5
log_Spot  0.04506413  0.137166040 -0.4650790  0.12114404  0.1268813
log_Fut   0.05321376  0.124802235 -0.3900878  0.02877604  0.1400407
log_USO   0.06759086  0.001857149 -0.1770362 -0.02323782  0.1077595
log_USL   0.02375506  0.113090122 -0.3237487  0.03741476  0.1330303
log_OIL   0.13485349 -0.050165498 -0.2990014 -0.01180909  0.1981438


Coefficient matrix of deterministic regressor(s).

             constant
log_Spot -0.007395204
log_Fut   0.003598548
log_USO   0.008296858
log_USL  -0.016685317
log_OIL   0.024704477
> # Step 5: Extract beta and alpha matrices
> beta <- johansen_test@V[1:r+1, 1:r]
> alpha <- johansen_test@W[, 1:r]
> print("Extracted beta matrix:")
[1] "Extracted beta matrix:"
> print(beta)
             log_Spot.l1 log_Fut.l1 log_USO.l1 log_USL.l1
log_Fut.l1 -1.0017520247 -1.3244152   16.80517  0.7523886
log_USO.l1  0.0023907033 -0.5014035  -74.82706  8.6459957
log_USL.l1 -0.0043867524  0.5028687  -35.29678 -6.3180603
log_OIL.l1  0.0005538807  0.2823013   76.35024 -4.8355656
> print("Extracted alpha matrix:")
[1] "Extracted alpha matrix:"
> print(alpha)
           log_Spot.l1 log_Fut.l1   log_USO.l1    log_USL.l1
log_Spot.d -0.24199295 0.02295218 0.0003860446 -0.0045865854
log_Fut.d   0.18580099 0.04212123 0.0004176988 -0.0042320772
log_USO.d  -0.15331674 0.06760043 0.0005287744 -0.0024737038
log_USL.d  -0.07075644 0.02880841 0.0006447661 -0.0007056863
log_OIL.d  -0.18973739 0.07470091 0.0002518853 -0.0017978899
> 
> # Step 6: Compute orthogonal complements beta_perp and alpha_perp
> # Use Null space computation from MASS package
> beta_perp <- Null(t(beta))
> alpha_perp <- Null(t(alpha))
> print("Computed orthogonal complement beta_perp:")
[1] "Computed orthogonal complement beta_perp:"
> print(beta_perp)
    
[1,]
[2,]
[3,]
[4,]
> print("Computed orthogonal complement alpha_perp:")
[1] "Computed orthogonal complement alpha_perp:"
> print(alpha_perp)
    
[1,]
[2,]
[3,]
[4,]
> 
> # Ensure beta_perp and alpha_perp are matrices
> beta_perp <- as.matrix(beta_perp)
> alpha_perp <- as.matrix(alpha_perp)
> print("Converted beta_perp to matrix:")
[1] "Converted beta_perp to matrix:"
> print(beta_perp)
    
[1,]
[2,]
[3,]
[4,]
> print("Converted alpha_perp to matrix:")
[1] "Converted alpha_perp to matrix:"
> print(alpha_perp)
    
[1,]
[2,]
[3,]
[4,]
> 
> # Step 7: Compute long-run impact matrix Psi(1)
> Psi_1 <- beta_perp %*% t(alpha_perp)
> print("Computed long-run impact matrix Psi(1):")
[1] "Computed long-run impact matrix Psi(1):"
> print(Psi_1)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0
> 
> # Step 8: Compute residual covariance matrix Omega
> residuals_var <- residuals(vecm_var)
> Omega <- cov(residuals_var)
> print("Computed residual covariance matrix Omega:")
[1] "Computed residual covariance matrix Omega:"
> print(Omega)
                   resids of log_Spot resids of log_Fut resids of log_USO resids of log_USL resids of log_OIL
resids of log_Spot       0.0005644515      0.0005321534      0.0004577378      0.0003636763      0.0004834179
resids of log_Fut        0.0005321534      0.0005445740      0.0004583258      0.0003635811      0.0004826371
resids of log_USO        0.0004577378      0.0004583258      0.0004424266      0.0003554061      0.0004709637
resids of log_USL        0.0003636763      0.0003635811      0.0003554061      0.0003310452      0.0003825727
resids of log_OIL        0.0004834179      0.0004826371      0.0004709637      0.0003825727      0.0005346706
> 
> # 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
> print("Computed residual correlation matrix:")
[1] "Computed residual correlation matrix:"
> print(corr_matrix)
                   resids of log_Spot resids of log_Fut resids of log_USO resids of log_USL resids of log_OIL
resids of log_Spot          1.0000000         0.9598318         0.9159733         0.8413140         0.8799668
resids of log_Fut           0.9598318         1.0000000         0.9337385         0.8563065         0.8944358
resids of log_USO           0.9159733         0.9337385         1.0000000         0.9286681         0.9683306
resids of log_USL           0.8413140         0.8563065         0.9286681         1.0000000         0.9093419
resids of log_OIL           0.8799668         0.8944358         0.9683306         0.9093419         1.0000000
> print("Eigenvalues of correlation matrix:")
[1] "Eigenvalues of correlation matrix:"
> print(eigen_result$values)
[1] 4.63622032 0.21050792 0.09032045 0.04009674 0.02285457
> print("Eigenvectors of correlation matrix:")
[1] "Eigenvectors of correlation matrix:"
> print(G)
           [,1]       [,2]        [,3]        [,4]        [,5]
[1,] -0.4434738  0.5646226 -0.18415788  0.66946900 -0.04928965
[2,] -0.4480840  0.4715324 -0.09924278 -0.70169019  0.27323925
[3,] -0.4580114 -0.1500695  0.26353846 -0.16465350 -0.81923222
[4,] -0.4373571 -0.5864481 -0.66839993  0.03065071  0.13076471
[5,] -0.4488832 -0.3040004  0.66334437  0.17717761  0.48442751
> 
> # Step 10: Compute Cholesky decomposition of Omega
> V <- t(chol(Omega))
> print("Computed Cholesky decomposition of Omega (matrix V):")
[1] "Computed Cholesky decomposition of Omega (matrix V):"
> print(V)
                   resids of log_Spot resids of log_Fut resids of log_USO resids of log_USL resids of log_OIL
resids of log_Spot         0.02375819       0.000000000       0.000000000      0.0000000000       0.000000000
resids of log_Fut          0.02239874       0.006547554       0.000000000      0.0000000000       0.000000000
resids of log_USO          0.01926653       0.004090063       0.007382338      0.0000000000       0.000000000
resids of log_USL          0.01530741       0.003163691       0.006440466      0.0067260564       0.000000000
resids of log_OIL          0.02034742       0.004105427       0.008418477      0.0005797034       0.005708903
> 
> # 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)
> print("Computed F^M inverse matrix:")
[1] "Computed F^M inverse matrix:"
> print(F_M_inv)
     resids of log_Spot resids of log_Fut resids of log_USO resids of log_USL resids of log_OIL
[1,]          394.28409        -219.65122         -4.123141         -20.28489         -44.44236
[2,]         -552.85929         628.98670        -86.952894         -21.32204         -32.12105
[3,]          -16.70373        -616.82561       1177.184073        -123.78615        -349.31488
[4,]           50.63079          77.57887       -354.546552         371.47613         -96.23931
[5,]           67.67344         184.29224       -845.963463        -134.04522         607.50859
> print("Computed F^M matrix:")
[1] "Computed F^M matrix:"
> print(F_M)
                         [,1]       [,2]        [,3]        [,4]        [,5]
resids of log_Spot 0.01591893 0.01091479 0.008741038 0.007265960 0.007918760
resids of log_Fut  0.01801606 0.01452761 0.010800059 0.008920024 0.009709158
resids of log_USO  0.01750445 0.01438370 0.013041080 0.010149610 0.011147485
resids of log_USL  0.01613659 0.01372343 0.013367719 0.013021034 0.011655206
resids of log_OIL  0.02069707 0.01743461 0.016859437 0.013491171 0.015913318
> 
> # Step 12: Compute psi_1^G
> psi_1_r <- Psi_1[1, ]  # First row of Psi(1)
> psi_1_G <- psi_1_r %*% F_M
Error in psi_1_r %*% F_M : non-conformable arguments