'Perform an iterative two-step procedure to estimate a WLS model.
The model is assumed to have heteroscedastic errors.
The variance is estimated by OLS regression of the link transformed
squared residuals on Z, i.e.::
link(sigma_i) = x_i*gamma.
Parameters
maxiter : integer, optional
the number of iterations
Notes
maxiter=1: returns the estimated based on given weights
maxiter=2: performs a second estimation with the updated weights,
this is 2-step estimation
maxiter>2: iteratively estimate and update the weights
TODO: possible extension stop iteration if change in parameter
estimates is smaller than x_tol
Repeated calls to fit_iterative, will do one redundant pinv_wexog
calculation. Calling fit_iterative(maxiter) ones does not do any
redundant recalculations (whitening or calculating pinv_wexog).'

def iterative_fit(self, maxiter=3):

    import collections
    self.history = collections.defaultdict(list)
    res_resid = None
    for i in range(maxiter):
        if hasattr(self, 'pinv_wexog'):
            del self.pinv_wexog
        results = self.fit()
        self.history['self_params'].append(results.params)
        if (not (i == (maxiter - 1))):
            self.results_old = results
            res_resid = OLS(self.link((results.resid ** 2)), self.exog_var).fit()
            self.history['ols_params'].append(res_resid.params)
            self.weights = (1.0 / self.linkinv(res_resid.fittedvalues))
            self.weights /= self.weights.max()
            self.weights[(self.weights < 1e-14)] = 1e-14
            self.initialize()
    results._results.results_residual_regression = res_resid
    return results
