docs: def correlation_matvec(P, obs1, obs2=None, times=[1]):
      r"""Time-correlation for equilibrium experiment - via matrix vector products.
  
      Parameters
      ----------
      P : (M, M) ndarray
          Transition matrix
      obs1 : (M,) ndarray
          Observable, represented as vector on state space
      obs2 : (M,) ndarray (optional)
          Second observable, for cross-correlations
      times : list of int (optional)
          List of times (in tau) at which to compute correlation
  
      Returns
      -------
      correlations : ndarray
          Correlation values at given times
  
      """
      if obs2 is None:
          obs2 = obs1
  
      """Compute stationary vector"""
      mu = statdist(P)
      obs1mu = mu * obs1
  
      times = np.asarray(times)
      """Sort in increasing order"""
      ind = np.argsort(times)
      times = times[ind]
  
      if times[0] < 0:
          raise ValueError("Times can not be negative")
      dt = times[1:] - times[0:-1]
  
      nt = len(times)
  
      correlations = np.zeros(nt)
  
      """Propagate obs2 to initial time"""
      obs2_t = 1.0 * obs2
      obs2_t = propagate(P, obs2_t, times[0])
      correlations[0] = np.dot(obs1mu, obs2_t)
      for i in range(nt - 1):
          obs2_t = propagate(P, obs2_t, dt[i])
          correlations[i + 1] = np.dot(obs1mu, obs2_t)
  
      """Cast back to original order of time points"""
      correlations = correlations[ind]
  
      return correlations