Generates a correlated noise vector using a multivariate normal
    random number generator with zero mean and covariance
    
        Sigma_ij = s^2 exp(-|t_i - t_j|/l),
    
    where s is the variance and l is the time scale. 
    The Power spectral density associated to  this covariance is 
    
        S(f) = 2*l*s^2/(4*pi^2*f^2*l^2 +1),
        
    red noise spectrum is defined as proportional to 1/f^2. This 
    covariance is the one expected from a first order markov process 
    (Reference?)
    
    Parameters
    ---------
    t: ndarray
        A time vector for which the red noise vector will be sampled
    variance: positive float
        variance of the resulting red noise vector
    time_scale: positive float
        Parameter of the covariance matrix
        
    Returns
    -------
    red_noise: ndarray
        Vector containing the red noise realizations
        
    See also
    --------
    power_law_noise

def first_order_markov_process(t, variance, time_scale, rseed=None):
    """
    Generates a correlated noise vector using a multivariate normal
    random number generator with zero mean and covariance
    
        Sigma_ij = s^2 exp(-|t_i - t_j|/l),
    
    where s is the variance and l is the time scale. 
    The Power spectral density associated to  this covariance is 
    
        S(f) = 2*l*s^2/(4*pi^2*f^2*l^2 +1),
        
    red noise spectrum is defined as proportional to 1/f^2. This 
    covariance is the one expected from a first order markov process 
    (Reference?)
    
    Parameters
    ---------
    t: ndarray
        A time vector for which the red noise vector will be sampled
    variance: positive float
        variance of the resulting red noise vector
    time_scale: positive float
        Parameter of the covariance matrix
        
    Returns
    -------
    red_noise: ndarray
        Vector containing the red noise realizations
        
    See also
    --------
    power_law_noise
        
    """
    if variance < 0.0:
        raise ValueError("Variance must be positive")
    if time_scale < 0.0:
        raise ValueError("Time scale must be positive")
    np.random.seed(rseed)
    N = len(t)
    mu = np.zeros(shape=(N,))
    if variance == 0.0:
        return mu
    dt = np.repeat(np.reshape(t, (1, -1)), N, axis=0)
    dt = np.absolute(dt - dt.T)  # This is NxN
    S = variance*np.exp(-np.absolute(dt)/time_scale)
    red_noise = np.random.multivariate_normal(mu, S)
    return red_noise