> # 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")]
> print("Selected columns for Y:")
[1] "Selected columns for Y:"
> print(Y)
# A tibble: 3,031 × 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
 7     4.12    4.11    4.63    3.13    5.63
 8     4.11    4.11    4.62    3.13    5.63
 9     4.11    4.11    4.62    3.13    5.63
10     4.10    4.10    4.61    3.12    5.61
# ℹ 3,021 more rows
# ℹ Use `print(n = ...)` to see more rows
> print("Dimensions of Y:")
[1] "Dimensions of Y:"
> print(dim(Y))
[1] 3031    5
> 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
Dickey-Fuller Z(alpha)  log_Spot
Dickey-Fuller Z(alpha)1  log_Fut
Dickey-Fuller Z(alpha)2  log_USO
Dickey-Fuller Z(alpha)3  log_USL
Dickey-Fuller Z(alpha)4  log_OIL
                        Test_Statistic
Dickey-Fuller Z(alpha)      -10.237353
Dickey-Fuller Z(alpha)1      -9.894884
Dickey-Fuller Z(alpha)2      -8.401415
Dickey-Fuller Z(alpha)3      -9.306146
Dickey-Fuller Z(alpha)4      -7.929339
                          p_value
Dickey-Fuller Z(alpha)  0.5387773
Dickey-Fuller Z(alpha)1 0.5578765
Dickey-Fuller Z(alpha)2 0.6411661
Dickey-Fuller Z(alpha)3 0.5907100
Dickey-Fuller Z(alpha)4 0.6674934
> print("Dimensions of PP results:")
[1] "Dimensions of PP results:"
> print(dim(pp_results))
[1] 5 3
> 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
> print("Dimensions of KPSS results:")
[1] "Dimensions of KPSS results:"
> print(dim(kpss_results))
[1] 5 3
> 
> # 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
Dickey-Fuller Z(alpha)  log_Spot
Dickey-Fuller Z(alpha)1  log_Fut
Dickey-Fuller Z(alpha)2  log_USO
Dickey-Fuller Z(alpha)3  log_USL
Dickey-Fuller Z(alpha)4  log_OIL
                        Test_Statistic
Dickey-Fuller Z(alpha)       -3041.242
Dickey-Fuller Z(alpha)1      -3164.641
Dickey-Fuller Z(alpha)2      -3217.069
Dickey-Fuller Z(alpha)3      -3208.247
Dickey-Fuller Z(alpha)4      -3212.187
                        p_value
Dickey-Fuller Z(alpha)     0.01
Dickey-Fuller Z(alpha)1    0.01
Dickey-Fuller Z(alpha)2    0.01
Dickey-Fuller Z(alpha)3    0.01
Dickey-Fuller Z(alpha)4    0.01
> print("Dimensions of PP results for differences:")
[1] "Dimensions of PP results for differences:"
> print(dim(pp_results_diff))
[1] 5 3
> 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
> print("Dimensions of KPSS results for differences:")
[1] "Dimensions of KPSS results for differences:"
> print(dim(kpss_results_diff))
[1] 5 3
> 
> # 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 
> print("Dimensions of VARselect result:")
[1] "Dimensions of VARselect result:"
> print(dim(VARselect_result$selection))
NULL
> 
> # 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)"]
> print("Chosen lag order based on SC:")
[1] "Chosen lag order based on SC:"
> print(p)
SC(n) 
    5 
> 
> # 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
[4] 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_Spot.l1  1.0000000000  1.0000000
log_Fut.l1  -1.0017520247 -1.3244152
log_USO.l1   0.0023907033 -0.5014035
log_USL.l1  -0.0043867524  0.5028687
log_OIL.l1   0.0005538807  0.2823013
constant     0.0069055220  0.5034032
            log_USO.l1 log_USL.l1 log_OIL.l1
log_Spot.l1    1.00000  1.0000000   1.000000
log_Fut.l1    16.80517  0.7523886  -8.329693
log_USO.l1   -74.82706  8.6459957  15.122284
log_USL.l1   -35.29678 -6.3180603  -4.415480
log_OIL.l1    76.35024 -4.8355656  -8.218566
constant     -47.90226 -0.2647043  16.978031
             constant
log_Spot.l1  1.000000
log_Fut.l1  -1.820012
log_USO.l1  -9.148170
log_USL.l1   5.683194
log_OIL.l1   3.311713
constant    10.890075

Weights W:
(This is the loading matrix)

           log_Spot.l1 log_Fut.l1
log_Spot.d -0.24199295 0.02295218
log_Fut.d   0.18580099 0.04212123
log_USO.d  -0.15331674 0.06760043
log_USL.d  -0.07075644 0.02880841
log_OIL.d  -0.18973739 0.07470091
             log_USO.l1    log_USL.l1
log_Spot.d 0.0003860446 -0.0045865854
log_Fut.d  0.0004176988 -0.0042320772
log_USO.d  0.0005287744 -0.0024737038
log_USL.d  0.0006447661 -0.0007056863
log_OIL.d  0.0002518853 -0.0017978899
             log_OIL.l1      constant
log_Spot.d 0.0001805646 -3.849043e-15
log_Fut.d  0.0001764496  1.940801e-15
log_USO.d  0.0001574821 -2.212215e-15
log_USL.d  0.0001551454 -1.849213e-15
log_OIL.d  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
> print("Number of cointegrating vectors (r):")
[1] "Number of cointegrating vectors (r):"
> print(r)
[1] 4
> 
> # Step 4: Estimate VECM and transform to VAR representation
> vecm_var <- vec2var(johansen_test, r=r)
> print("Dimensions of the VECM VAR representation:")
[1] "Dimensions of the VECM VAR representation:"
> print(dim(vecm_var))
NULL
> 
> # 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]
> print("Extracted beta matrix:")
[1] "Extracted beta matrix:"
> print(beta)
              log_Spot.l1 log_Fut.l1
log_Spot.l1  1.0000000000  1.0000000
log_Fut.l1  -1.0017520247 -1.3244152
log_USO.l1   0.0023907033 -0.5014035
log_USL.l1  -0.0043867524  0.5028687
log_OIL.l1   0.0005538807  0.2823013
            log_USO.l1 log_USL.l1
log_Spot.l1    1.00000  1.0000000
log_Fut.l1    16.80517  0.7523886
log_USO.l1   -74.82706  8.6459957
log_USL.l1   -35.29678 -6.3180603
log_OIL.l1    76.35024 -4.8355656
> print("Dimensions of beta matrix:")
[1] "Dimensions of beta matrix:"
> print(dim(beta))
[1] 5 4
> print("Extracted alpha matrix:")
[1] "Extracted alpha matrix:"
> print(alpha)
           log_Spot.l1 log_Fut.l1
log_Spot.d -0.24199295 0.02295218
log_Fut.d   0.18580099 0.04212123
log_USO.d  -0.15331674 0.06760043
log_USL.d  -0.07075644 0.02880841
log_OIL.d  -0.18973739 0.07470091
             log_USO.l1    log_USL.l1
log_Spot.d 0.0003860446 -0.0045865854
log_Fut.d  0.0004176988 -0.0042320772
log_USO.d  0.0005287744 -0.0024737038
log_USL.d  0.0006447661 -0.0007056863
log_OIL.d  0.0002518853 -0.0017978899
> print("Dimensions of alpha matrix:")
[1] "Dimensions of alpha matrix:"
> print(dim(alpha))
[1] 5 4
> 
> # 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
> print("Computed beta_perp:")
[1] "Computed beta_perp:"
> print(beta_perp)
[1] -0.6856252 -0.6210203  0.2081991
[4] -0.1225053  0.2930828
> print("Dimensions of beta_perp:")
[1] "Dimensions of beta_perp:"
> print(dim(beta_perp))
NULL
> 
> t_alpha_svd <- svd(t(alpha))
> alpha_perp <- t_alpha_svd$v[, ncol(t_alpha_svd$v)]  # Similarly for alpha
> print("Computed alpha_perp:")
[1] "Computed alpha_perp:"
> print(alpha_perp)
[1]  0.01412683  0.02298025 -0.10644853
[4] -0.89827793  0.42548538
> print("Dimensions of alpha_perp:")
[1] "Dimensions of alpha_perp:"
> print(dim(alpha_perp))
NULL
> 
> # Normalize beta_perp and alpha_perp
> beta_perp <- beta_perp / beta_perp[1]
> alpha_perp <- alpha_perp / alpha_perp[1]
> print("Normalized beta_perp:")
[1] "Normalized beta_perp:"
> print(beta_perp)
[1]  1.0000000  0.9057723 -0.3036631
[4]  0.1786767 -0.4274679
> print("Normalized alpha_perp:")
[1] "Normalized alpha_perp:"
> print(alpha_perp)
[1]   1.000000   1.626709  -7.535201
[4] -63.586644  30.118949
> 
> # 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)
> print("Computed long-run impact matrix Psi(1):")
[1] "Computed long-run impact matrix Psi(1):"
> print(Psi_1)
             [,1]        [,2]        [,3]
[1,] -0.051348548 -0.08352914  0.38692162
[2,] -0.046510093 -0.07565839  0.35046289
[3,]  0.015592659  0.02536472 -0.11749381
[4,] -0.009174791 -0.01492471  0.06913389
[5,]  0.021949856  0.03570603 -0.16539657
           [,4]       [,5]
[1,]  3.2650818 -1.5465643
[2,]  2.9574207 -1.4008351
[3,] -0.9914848  0.4696345
[4,]  0.5833942 -0.2763351
[5,] -1.3957177  0.6611066
> print("Dimensions of Psi(1):")
[1] "Dimensions of Psi(1):"
> print(dim(Psi_1))
[1] 5 5
> 
> # 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_Spot       0.0005644515
resids of log_Fut        0.0005321534
resids of log_USO        0.0004577378
resids of log_USL        0.0003636763
resids of log_OIL        0.0004834179
                   resids of log_Fut
resids of log_Spot      0.0005321534
resids of log_Fut       0.0005445740
resids of log_USO       0.0004583258
resids of log_USL       0.0003635811
resids of log_OIL       0.0004826371
                   resids of log_USO
resids of log_Spot      0.0004577378
resids of log_Fut       0.0004583258
resids of log_USO       0.0004424266
resids of log_USL       0.0003554061
resids of log_OIL       0.0004709637
                   resids of log_USL
resids of log_Spot      0.0003636763
resids of log_Fut       0.0003635811
resids of log_USO       0.0003554061
resids of log_USL       0.0003310452
resids of log_OIL       0.0003825727
                   resids of log_OIL
resids of log_Spot      0.0004834179
resids of log_Fut       0.0004826371
resids of log_USO       0.0004709637
resids of log_USL       0.0003825727
resids of log_OIL       0.0005346706
> print("Dimensions of Omega:")
[1] "Dimensions of Omega:"
> print(dim(Omega))
[1] 5 5
> 
> # 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_Spot          1.0000000
resids of log_Fut           0.9598318
resids of log_USO           0.9159733
resids of log_USL           0.8413140
resids of log_OIL           0.8799668
                   resids of log_Fut
resids of log_Spot         0.9598318
resids of log_Fut          1.0000000
resids of log_USO          0.9337385
resids of log_USL          0.8563065
resids of log_OIL          0.8944358
                   resids of log_USO
resids of log_Spot         0.9159733
resids of log_Fut          0.9337385
resids of log_USO          1.0000000
resids of log_USL          0.9286681
resids of log_OIL          0.9683306
                   resids of log_USL
resids of log_Spot         0.8413140
resids of log_Fut          0.8563065
resids of log_USO          0.9286681
resids of log_USL          1.0000000
resids of log_OIL          0.9093419
                   resids of log_OIL
resids of log_Spot         0.8799668
resids of log_Fut          0.8944358
resids of log_USO          0.9683306
resids of log_USL          0.9093419
resids of log_OIL          1.0000000
> print("Dimensions of correlation matrix:")
[1] "Dimensions of correlation matrix:"
> print(dim(corr_matrix))
[1] 5 5
> print("Eigenvalues of correlation matrix:")
[1] "Eigenvalues of correlation matrix:"
> print(eigen_result$values)
[1] 4.63622032 0.21050792 0.09032045
[4] 0.04009674 0.02285457
> print("Dimensions of eigenvalues:")
[1] "Dimensions of eigenvalues:"
> print(length(eigen_result$values))
[1] 5
> print("Eigenvectors of correlation matrix:")
[1] "Eigenvectors of correlation matrix:"
> print(G)
           [,1]       [,2]        [,3]
[1,] -0.4434738  0.5646226 -0.18415788
[2,] -0.4480840  0.4715324 -0.09924278
[3,] -0.4580114 -0.1500695  0.26353846
[4,] -0.4373571 -0.5864481 -0.66839993
[5,] -0.4488832 -0.3040004  0.66334437
            [,4]        [,5]
[1,]  0.66946900 -0.04928965
[2,] -0.70169019  0.27323925
[3,] -0.16465350 -0.81923222
[4,]  0.03065071  0.13076471
[5,]  0.17717761  0.48442751
> print("Dimensions of eigenvectors:")
[1] "Dimensions of eigenvectors:"
> print(dim(G))
[1] 5 5
> 
> # 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_Spot         0.02375819
resids of log_Fut          0.02239874
resids of log_USO          0.01926653
resids of log_USL          0.01530741
resids of log_OIL          0.02034742
                   resids of log_Fut
resids of log_Spot       0.000000000
resids of log_Fut        0.006547554
resids of log_USO        0.004090063
resids of log_USL        0.003163691
resids of log_OIL        0.004105427
                   resids of log_USO
resids of log_Spot       0.000000000
resids of log_Fut        0.000000000
resids of log_USO        0.007382338
resids of log_USL        0.006440466
resids of log_OIL        0.008418477
                   resids of log_USL
resids of log_Spot      0.0000000000
resids of log_Fut       0.0000000000
resids of log_USO       0.0000000000
resids of log_USL       0.0067260564
resids of log_OIL       0.0005797034
                   resids of log_OIL
resids of log_Spot       0.000000000
resids of log_Fut        0.000000000
resids of log_USO        0.000000000
resids of log_USL        0.000000000
resids of log_OIL        0.005708903
> print("Dimensions of matrix V:")
[1] "Dimensions of matrix V:"
> print(dim(V))
[1] 5 5
> 
> # 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
[1,]          394.28409        -219.65122
[2,]         -552.85929         628.98670
[3,]          -16.70373        -616.82561
[4,]           50.63079          77.57887
[5,]           67.67344         184.29224
     resids of log_USO resids of log_USL
[1,]         -4.123141         -20.28489
[2,]        -86.952894         -21.32204
[3,]       1177.184073        -123.78615
[4,]       -354.546552         371.47613
[5,]       -845.963463        -134.04522
     resids of log_OIL
[1,]         -44.44236
[2,]         -32.12105
[3,]        -349.31488
[4,]         -96.23931
[5,]         607.50859
> print("Dimensions of F^M inverse matrix:")
[1] "Dimensions of F^M inverse matrix:"
> print(dim(F_M_inv))
[1] 5 5
> print("Computed F^M matrix:")
[1] "Computed F^M matrix:"
> print(F_M)
                         [,1]       [,2]
resids of log_Spot 0.01591893 0.01091479
resids of log_Fut  0.01801606 0.01452761
resids of log_USO  0.01750445 0.01438370
resids of log_USL  0.01613659 0.01372343
resids of log_OIL  0.02069707 0.01743461
                          [,3]        [,4]
resids of log_Spot 0.008741038 0.007265960
resids of log_Fut  0.010800059 0.008920024
resids of log_USO  0.013041080 0.010149610
resids of log_USL  0.013367719 0.013021034
resids of log_OIL  0.016859437 0.013491171
                          [,5]
resids of log_Spot 0.007918760
resids of log_Fut  0.009709158
resids of log_USO  0.011147485
resids of log_USL  0.011655206
resids of log_OIL  0.015913318
> print("Dimensions of F^M matrix:")
[1] "Dimensions of F^M matrix:"
> print(dim(F_M))
[1] 5 5
> 
> # 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  # This should now work without dimension errors
> print("Computed psi_1^G:")
[1] "Computed psi_1^G:"
> print(psi_1_G)
           [,1]       [,2]       [,3]
[1,] 0.02512852 0.02163579 0.02126741
          [,4]       [,5]
[1,] 0.0244587 0.01653982
> print("Dimensions of psi_1^G:")
[1] "Dimensions of psi_1^G:"
> print(dim(psi_1_G))
[1] 1 5
> 
> # Step 13: Compute numerator and denominator for GIS calculation
> numerator <- psi_1_G^2
> denominator <- as.numeric((psi_1_r) %*% (Omega) %*% t(psi_1_r))
Error in (psi_1_r) %*% (Omega) %*% t(psi_1_r) : non-conformable arguments