R implementation of the Matlab Mann-Kendall test function

matlab_mann_kendall.R

matlab_mann_kendall <- function(V,
                                alpha = 0.05) {
  
  alpha <- alpha/2
  n <- length(V)
  S <- 0
  
  for (i in 1:(n-1)) {
    for(j in (i+1):n) {
      S <- S + sign(V[j] - V[i])
    }
  }
  
  VarS <- (n * (n - 1) * (2 * n + 5)) / 18
  
  StdS <- sqrt(VarS)
  
  # Note: ties are not considered 
  if (S >= 0) {
    Z <- (S - 1) / StdS * (sum(S != 0))
  } else {
    Z <- (S + 1) / StdS
  }
  
  p_value <- 2 * (1 - pnorm(abs(Z), 0, 1)) # Two-tailed test 
  pz <- qnorm(1 - alpha, 0, 1) 
  H <- as.integer(abs(Z) > pz)
  
  return(list(H = H,
              VarS = VarS,
              p_value = p_value))
}

Matlab version:
https://www.mathworks.com/matlabcentral/fileexchange/25531-mann-kendall-test

As noted by Nick Brown on Twitter:
https://twitter.com/sTeamTraen/status/1118218894098030593

This version of the test gives a different result than the Kendall::MannKendall() function.

I haven't been able to pin down the difference, but found that this version and Kendall::MannKendall() have different values for VarS (varS in MannKendall() results). If the MannKendall() varS value is injected into this function in place of the VarS calculation, the Kendall::MannKendall() p-value is returned, so I suspect that this or the S counting step provide the key difference.