'Perform an iterative two-stage procedure to estimate a GLS model.
The model is assumed to have AR(p) errors, AR(p) parameters and
regression coefficients are estimated iteratively.
Parameters
maxiter : integer, optional
the number of iterations
rtol : float, optional
Relative tolerance between estimated coefficients to stop the
estimation.  Stops if
max(abs(last - current) / abs(last)) < rtol'

def iterative_fit(self, maxiter=3, rtol=0.0001, **kwds):

    converged = False
    i = (-1)
    history = {'params': [], 'rho': [self.rho]}
    for i in range((maxiter - 1)):
        if hasattr(self, 'pinv_wexog'):
            del self.pinv_wexog
        self.initialize()
        results = self.fit()
        history['params'].append(results.params)
        if (i == 0):
            last = results.params
        else:
            diff = np.max((np.abs((last - results.params)) / np.abs(last)))
            if (diff < rtol):
                converged = True
                break
            last = results.params
        (self.rho, _) = yule_walker(results.resid, order=self.order, df=None)
        history['rho'].append(self.rho)
    if ((not converged) and (maxiter > 0)):
        if hasattr(self, 'pinv_wexog'):
            del self.pinv_wexog
        self.initialize()
    results = self.fit(history=history, **kwds)
    results.iter = (i + 1)
    if (not converged):
        results.history['params'].append(results.params)
        results.iter += 1
    results.converged = converged
    return results
