def bivar_central_moment(x, y, pow=1, ddof=1):
    """Computes the bivariate central moment.

    Default parameters compute sample covariance. Two-pass algorithm for stability.

    Parameters
    ----------
    x, y : 1D array-likes
        x and y values to be compared
    pow : int, optional
        Power to raise each sum to, by default 1
    ddof : int, optional
        Degrees of freedom correction, by default 1

    Returns
    -------
    float64
        Product central moment of x, y

    Notes
    -----
    This ddof correction is only valid for the first product central moment (covariance).
    Simply doing ddof=2 does not provide an unbiased estimator for higher order moments.
    """
    n = len(x)

    mean_x = mean_y = 0
    for i in range(n):
        mean_x += x[i]
        mean_y += y[i]
    mean_x /= n
    mean_y /= n

    sum_of_prods = 0
    for x_, y_ in zip(x, y):
        sum_of_prods += ((x_ - mean_x) ** pow) * ((y_ - mean_y) ** pow)

    moment = sum_of_prods / (n - ddof)

    return moment