def _construct_var_eqns(data, p, delta=None):
        """Construct VAR equation system (optionally with RLS constraint).
        """
        t, m, l = np.shape(data)
        n = (l - p) * t  # number of linear relations
        rows = n if delta is None else n + m * p

        # Construct matrix x (predictor variables)
        x = np.zeros((rows, m * p))
        for i in range(m):
            for k in range(1, p + 1):
                x[:n, i * p + k - 1] = np.reshape(data[:, i, p - k:-k].T, n)
        if delta is not None:
            np.fill_diagonal(x[n:, :], delta)

        # Construct vectors yi (response variables for each channel i)
        y = np.zeros((rows, m))
        for i in range(m):
            y[:n, i] = np.reshape(data[:, i, p:].T, n)

        return x, y