'igls : bool
Iterate until estimates converge if sigma is None instead of
two-step GLS, which is the default is sigma is None.
tol : float
maxiter : int
Notes
This ia naive implementation that does not exploit the block
diagonal structure. It should work for ill-conditioned `sigma`
but this is untested.'

def fit(self, igls=False, tol=1e-05, maxiter=100):

    if (not np.any(self.sigma)):
        self.sigma = self._compute_sigma(self.endog, self.exog)
    M = self._M
    beta = np.dot(self.pinv_wexog, self.wendog)
    self._update_history(beta)
    self.iterations += 1
    if (not igls):
        sur_fit = SysResults(self, beta, self.normalized_cov_params)
        return sur_fit
    conv = self.history['params']
    while (igls and np.any((np.abs((conv[(-2)] - conv[(-1)])) > tol)) and (self.iterations < maxiter)):
        fittedvalues = (self.sp_exog * beta).reshape(M, (-1))
        resids = (self.endog - fittedvalues)
        self.sigma = self._compute_sigma(resids)
        self.wendog = self.whiten(self.endog)
        self.wexog = self.whiten(self.sp_exog)
        self.pinv_wexog = np.linalg.pinv(self.wexog)
        self.normalized_cov_params = np.dot(self.pinv_wexog, np.transpose(self.pinv_wexog))
        beta = np.dot(self.pinv_wexog, self.wendog)
        self._update_history(beta)
        self.iterations += 1
    sur_fit = SysResults(self, beta, self.normalized_cov_params)
    return sur_fit
