def surrogate_connectivity(measure_names, data, var, nfft=512, repeats=100,
                           n_jobs=1, verbose=0, random_state=None):
    """Calculate surrogate connectivity for a multivariate time series by phase
    randomization [1]_.

    .. note:: Parameter `var` will be modified by the function. Treat as
    undefined after the function returns.

    Parameters
    ----------
    measures : str or list of str
        Name(s) of the connectivity measure(s) to calculate. See
        :class:`Connectivity` for supported measures.
    data : array, shape (trials, channels, samples) or (channels, samples)
        Time series data (2D or 3D for multiple trials)
    var : VARBase-like object
        Instance of a VAR model.
    nfft : int, optional
        Number of frequency bins to calculate. Note that these points cover the
        range between 0 and half the sampling rate.
    repeats : int, optional
        Number of surrogate samples to take.
    n_jobs : int | None, optional
        Number of jobs to run in parallel. If set to None, joblib is not used
        at all. See `joblib.Parallel` for details.
    verbose : int, optional
        Verbosity level passed to joblib.

    Returns
    -------
    result : array, shape (`repeats`, n_channels, n_channels, nfft)
        Values of the connectivity measure for each surrogate. If
        `measure_names` is a list of strings a dictionary is returned, where
        each key is the name of the measure, and the corresponding values are
        arrays of shape (`repeats`, n_channels, n_channels, nfft).

    .. [1] J. Theiler et al. Testing for nonlinearity in time series: the
           method of surrogate data. Physica D, 58: 77-94, 1992.
    """
    par, func = parallel_loop(_calc_surrogate, n_jobs=n_jobs, verbose=verbose)
    output = par(func(randomize_phase(data, random_state=random_state), var,
                      measure_names, nfft) for _ in range(repeats))
    return convert_output_(output, measure_names)