Simulate vector autoregressive (VAR) model.

        This function generates data from the VAR model.

        Parameters
        ----------
        l : int or [int, int]
            Number of samples to generate. Can be a tuple or list, where l[0]
            is the number of samples and l[1] is the number of trials.
        noisefunc : func, optional
            This function is used to create the generating noise process. If
            set to None, Gaussian white noise with zero mean and unit variance
            is used.

        Returns
        -------
        data : array, shape (n_trials, n_samples, n_channels)
            Generated data.

def simulate(self, l, noisefunc=None, random_state=None):
        """Simulate vector autoregressive (VAR) model.

        This function generates data from the VAR model.

        Parameters
        ----------
        l : int or [int, int]
            Number of samples to generate. Can be a tuple or list, where l[0]
            is the number of samples and l[1] is the number of trials.
        noisefunc : func, optional
            This function is used to create the generating noise process. If
            set to None, Gaussian white noise with zero mean and unit variance
            is used.

        Returns
        -------
        data : array, shape (n_trials, n_samples, n_channels)
            Generated data.
        """
        m, n = np.shape(self.coef)
        p = n // m

        try:
            l, t = l
        except TypeError:
            t = 1

        if noisefunc is None:
            rng = check_random_state(random_state)
            noisefunc = lambda: rng.normal(size=(1, m))

        n = l + 10 * p

        y = np.zeros((n, m, t))
        res = np.zeros((n, m, t))

        for s in range(t):
            for i in range(p):
                e = noisefunc()
                res[i, :, s] = e
                y[i, :, s] = e
            for i in range(p, n):
                e = noisefunc()
                res[i, :, s] = e
                y[i, :, s] = e
                for k in range(1, p + 1):
                    y[i, :, s] += self.coef[:, (k - 1)::p].dot(y[i - k, :, s])

        self.residuals = res[10 * p:, :, :].T
        self.rescov = sp.cov(cat_trials(self.residuals).T, rowvar=False)

        return y[10 * p:, :, :].transpose([2, 1, 0])