docs: def _calc_q_statistic(x, h, nt):
      """Calculate Portmanteau statistics up to a lag of h.
      """
      t, m, n = x.shape
  
      # covariance matrix of x
      c0 = acm(x, 0)
  
      # LU factorization of covariance matrix
      c0f = sp.linalg.lu_factor(c0, overwrite_a=False, check_finite=True)
  
      q = np.zeros((3, h + 1))
      for l in range(1, h + 1):
          cl = acm(x, l)
  
          # calculate tr(cl' * c0^-1 * cl * c0^-1)
          a = sp.linalg.lu_solve(c0f, cl)
          b = sp.linalg.lu_solve(c0f, cl.T)
          tmp = a.dot(b).trace()
  
          # Box-Pierce
          q[0, l] = tmp
  
          # Ljung-Box
          q[1, l] = tmp / (nt - l)
  
          # Li-McLeod
          q[2, l] = tmp
  
      q *= nt
      q[1, :] *= (nt + 2)
  
      q = np.cumsum(q, axis=1)
  
      for l in range(1, h+1):
          q[2, l] = q[0, l] + m * m * l * (l + 1) / (2 * nt)
  
      return q