How to make script run more efficiently "Modified Baek and Brock Nonlinear Granger Causality Test in R"

Question

I am trying to implement the test propesed by Hiemstra and Jones (1994) for nonlinear granger causality test between two time series xt and yt.

I am struggling to find any relevant packages or similar forum discussions on the topic.

The paper: "Testing for Linear and Nonlinear Granger Causality in the Stock Price- Volume Relation" Author(s): Craig Hiemstra and Jonathan D. Jones Source: The Journal of Finance , Dec., 1994, Vol. 49, No. 5 (Dec., 1994), pp. 1639-1664

I am struggling with efficiency. The script as is, takes about 1 hour to run, I need to run the script 16 times. Does anyone know of how to make the script run more efficiently?

To clarify I have added screenshots of the equations I want to code. Equations 5:6 with descriptions Equations 7:9 with descriptinosAppendix variance estimator A hats

I tried this code, however I am sure if it is coded correctly wrt what I described:

obs <- 800
set.seed(123)
x <- rnorm(obs)
y <- rnorm(obs)
Lx <- 8
Ly <- Lx

x         <- scale(x, center = TRUE, scale = TRUE)
y         <- scale(y, center = TRUE, scale = TRUE)

# Parameters
T <- length(x)
n <- T + 1 - m - max(Lx, Ly)
e <- sd(x) * e
t <- max(Lx, Ly) + 1
s <- t


# Function for Maximum Norm Distance
dist.func <- function(a, b) {
  max(abs(a-b))
}

# Indicator Function
I.funct <- function(a, b, c) {
  as.numeric(dist.func(a, b) < c)
}

## Calculating Correlation Integral
corr.int <- function(x, y, m, Lx, Ly, e){
  T <- length(x) #nrow()!
  n <- T + 1 - m - max(Lx, Ly)
  # t <- max(Lx, Ly) + 1
  # s <- t
  
  C1 <- 0
  C2 <- 0
  C3 <- 0
  C4 <- 0
  for (t in (max(Lx, Ly) + 1):(T-m+1)) {
    for (s in (max(Lx, Ly) + 1):(T-m+1)) {
      if (t < s){
        C1 <- C1 + (I.funct(x[(t-Lx):(t-Lx+m+Lx-1)], x[(s-Lx):(s-Lx+m+Lx-1)], e) * I.funct(y[(t-Ly):(t-Ly+Ly-1)], y[(s-Ly):(s-Ly+Ly-1)] , e))
        C2 <- C2 + (I.funct(x[(t-Lx):(t-Lx+Lx-1)],   x[(s-Lx):(s-Lx+Lx-1)] , e)  * I.funct(y[(t-Ly):(t-Ly+Ly-1)], y[(s-Ly):(s-Ly+Ly-1)] , e))
        C3 <- C3 + (I.funct(x[(t-Lx):(t-Lx+m+Lx-1)], x[(s-Lx):(s-Lx+m+Lx-1)], e))
        C4 <- C4 + (I.funct(x[(t-Lx):(t-Lx+Lx-1)],   x[(s-Lx):(s-Lx+Lx-1)] , e))
      }
    }
  }
  C1 <- C1 * (2/(n*(n-1)))
  C2 <- C2 * (2/(n*(n-1)))
  C3 <- C3 * (2/(n*(n-1)))
  C4 <- C4 * (2/(n*(n-1)))
  
  test.stat <- sqrt(n)*((C1/C2)-(C3/C4))
  test.stat
  
  d_hat  <- as.matrix(cbind((1/C2), (-C2/(C2^2)), (-1/C4), (C3/(C4^2))))
  d_hatt <- t(d_hat)
  
  return(list(C1=C1, C2=C2, C3=C3, C4=C4, test.stat=test.stat))
  
}

eq.9 <- corr.int(x = x, y = y, m = m, Lx = Lx, Ly = Ly, e = e)

C1 <- eq.9$C1
C2 <- eq.9$C2
C3 <- eq.9$C3
C4 <- eq.9$C4

test.stat <- eq.9$test.stat

####
## Variance Estimator
# The bread
d_hat  <- as.matrix(cbind((1/C2), (-C2/(C2^2)), (-1/C4), (C3/(C4^2))))
d_hatt <- t(d_hat)
Kn <- as.integer(n^(1/4))


A.hat <- function(case,t,n) {
  A1 <- 0
  A2 <- 0
  A3 <- 0
  A4 <- 0
  
  for (s in (max(Lx, Ly) + 1):(T-m+1)) {
    if (t != s){
      A1 <- A1 + (I.funct(x[(t-Lx):(t-Lx+m+Lx-1)], x[(s-Lx):(s-Lx+m+Lx-1)], e) * I.funct(y[(t-Ly):(t-Ly+Ly-1)], y[(s-Ly):(s-Ly+Ly-1)] , e))
      A2 <- A2 + (I.funct(x[(t-Lx):(t-Lx+Lx-1)],   x[(s-Lx):(s-Lx+Lx-1)] , e)  * I.funct(y[(t-Ly):(t-Ly+Ly-1)], y[(s-Ly):(s-Ly+Ly-1)] , e))
      A3 <- A3 + (I.funct(x[(t-Lx):(t-Lx+m+Lx-1)], x[(s-Lx):(s-Lx+m+Lx-1)], e))
      A4 <- A4 + (I.funct(x[(t-Lx):(t-Lx+Lx-1)],   x[(s-Lx):(s-Lx+Lx-1)] , e))
    }
    
  }
  
  A1 <- A1 * (1/((n-1))) - C1
  A2 <- A2 * (1/((n-1))) - C2
  A3 <- A3 * (1/((n-1))) - C3
  A4 <- A4 * (1/((n-1))) - C4
  
  if (case == 1)  {
    return(A1)
  }
  
  if (case == 2)  {
    return(A2)
  }
  
  if (case == 3)  {
    return(A3)
  }
  
  if (case == 4)  {
    return(A4)
  }
  
}


omega <- function(k,n) {
  2*(1-((k-1)/Kn))*(k != 1) + (k == 1)*1
}

Kn <- as.integer(n^(1/4))

k <- seq(1:Kn)


sigmaHat.compute <- function(i,j,n) {
  
  term1 <- omega(k,n)*(1/(2*n-k+1))
  
  term2 <- NULL
  for(h in (1:length(k)))  {
    print(h)
    res <- 0
    for(t in (max(Lx, Ly) + k[h]):(T-m+1)) {
      res <- res + A.hat(i,t,n)*A.hat(j,t-k[h]+1,n) + A.hat(i,t-k[h]+1,n)*A.hat(j,t,n) #insert correct indices here
    }
    term2[h] <- res
  }
  
  return(4*sum(term1*term2)) #sigma hat i j n
}


testComp1 <- sigmaHat.compute(1,1,n)
testComp1


sigma_hat <- matrix(c( sigmaHat.compute(1,1,n), sigmaHat.compute(1,2,n), sigmaHat.compute(1,3,n), sigmaHat.compute(1,4,n), 
                       sigmaHat.compute(2,1,n), sigmaHat.compute(2,2,n), sigmaHat.compute(2,3,n), sigmaHat.compute(2,4,n), 
                       sigmaHat.compute(3,1,n), sigmaHat.compute(3,2,n), sigmaHat.compute(3,3,n), sigmaHat.compute(3,4,n), 
                       sigmaHat.compute(4,1,n), sigmaHat.compute(4,2,n), sigmaHat.compute(4,3,n), sigmaHat.compute(4,4,n) ), ncol=4, nrow=4, byrow = T)


sandwich <- d_hat %*% sigma_hat %*% d_hatt


test.statistic <- test.stat/sandwich
test.statistic