docs: def _sigmainf(N, h, m, dW, Km0, Pm0):
      """Asymptotic covariance matrix \Sigma_\infty  Wiktorsson2001 eqn (4.5)"""
      M = m*(m-1)//2
      Im = broadcast_to(np.eye(m), (N, m, m))
      IM = broadcast_to(np.eye(M), (N, M, M))
      Ims0 = np.eye(m**2)
      factor1 = broadcast_to((2.0/h)*np.dot(Km0, Ims0 - Pm0), (N, M, m**2))
      factor2 = _kp2(Im, _dot(dW, _t(dW)))
      factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
      return 2*IM + _dot(_dot(factor1, factor2), factor3)