How to use the pysal.model.spreg.regimes function in pysal

# LR test on spatial autoregressive coefficients
        if spat_diag:
            liklambda = 2.0 * (self.surerrllik - self.llik)
            plik2 = stats.chisqprob(liklambda, self.n_eq)
            self.likrlambda = (liklambda,self.n_eq,plik2)
        else:
            self.likrlambda = None
        
        # asymptotic variance for spatial coefficient
        if vm:
            self.vm = surerrvm(self.n,self.n_eq,w,self.lamsur,self.sig)
            vlam = self.vm[:self.n_eq,:self.n_eq]
            self.lamsetp = lam_setp(self.lamsur,vlam)
            # test on constancy of lambdas
            R = REGI.buildR(kr=1,kf=0,nr=self.n_eq)
            w,p = REGI.wald_test(self.lamsur,R,np.zeros((R.shape[0],1)),vlam)
            self.lamtest = (w,R.shape[0],p)
            if spat_diag:  # test on joint significance of lambdas
                Rj = np.identity(self.n_eq)
                wj,pj = REGI.wald_test(self.lamsur,Rj,np.zeros((Rj.shape[0],1)),vlam)
                self.joinlam = (wj,Rj.shape[0],pj)
            else:
                self.joinlam = None
        else:
            self.vm = None
            self.lamsetp = None
            self.lamtest = None
            self.joinlam = None

        # test on constancy of regression coefficients across equations
        if check_k(self.bigK):   # only for equal number of variables

if regime_err_sep == True:
            if set(cols2regi) == set([True]):
                self._error_regimes_multi(y, x, regimes, w, cores,
                                          cols2regi, vm, name_x)
            else:
                raise Exception("All coefficients must vary accross regimes if regime_err_sep = True.")
        else:
            self.x, self.name_x = REGI.Regimes_Frame.__init__(self, x_constant,
                                                              regimes, constant_regi=None, cols2regi=cols2regi, names=name_x)
            ols = BaseOLS(y=y, x=self.x)
            self.k = ols.x.shape[1]
            moments = _momentsGM_Error(w, ols.u)
            lambda1 = optim_moments(moments)
            xs = get_spFilter(w, lambda1, x_constant)
            ys = get_spFilter(w, lambda1, y)
            xs = REGI.Regimes_Frame.__init__(self, xs,
                                             regimes, constant_regi=None, cols2regi=cols2regi)[0]
            ols2 = BaseOLS(y=ys, x=xs)

            # Output
            self.predy = spdot(self.x, ols2.betas)
            self.u = y - self.predy
            self.betas = np.vstack((ols2.betas, np.array([[lambda1]])))
            self.sig2 = ols2.sig2n
            self.e_filtered = self.u - lambda1 * lag_spatial(w, self.u)
            self.vm = self.sig2 * ols2.xtxi
            self.title = "SPATIALLY WEIGHTED LEAST SQUARES - REGIMES"
            self.name_x.append('lambda')
            self.kf += 1
            self.chow = REGI.Chow(self)
            self._cache = {}
            SUMMARY.GM_Error(reg=self, w=w, vm=vm, regimes=True)

reg  : BaseSUR object

       Creates
       -------

       reg.bOLS    : dictionary with regression coefficients for each equation
       reg.olsE    : N x n_eq array with OLS residuals for each equation

    '''
    reg.bOLS = {}
    for r in range(reg.n_eq):
        reg.bOLS[r] = np.dot(la.inv(reg.bigXX[(r,r)]),reg.bigXy[(r,r)])
    reg.olsE = sur_resids(reg.bigy,reg.bigX,reg.bOLS)
    return reg

class SUR(BaseSUR, REGI.Regimes_Frame):
    """ User class for SUR estimation, both two step as well as iterated

        Parameters
        ----------

        bigy       : dictionary with vector for dependent variable by equation
        bigX       : dictionary with matrix of explanatory variables by equation
                     (note, already includes constant term)
        w          : spatial weights object, default = None
        regimes    : list; default = None
                     List of n values with the mapping of each
                     observation to a regime. Assumed to be aligned with 'x'.
        nonspat_diag : boolean; flag for non-spatial diagnostics, default = True
        spat_diag  : boolean; flag for spatial diagnostics, default = False
        iter       : boolean; whether or not to use iterated estimation
                     default = False

self.regime_err_sep = regime_err_sep

        if regime_err_sep == True:
            if set(cols2regi) == set([True]):
                self._endog_error_regimes_multi(y, x, regimes, w, yend, q, cores,
                                                cols2regi, vm, name_x, name_yend, name_q, add_lag)
            else:
                raise Exception("All coefficients must vary accross regimes if regime_err_sep = True.")
        else:
            x_constant = USER.check_constant(x)
            q, name_q = REGI.Regimes_Frame.__init__(self, q,
                                                    regimes, constant_regi=None, cols2regi='all', names=name_q)
            x, name_x = REGI.Regimes_Frame.__init__(self, x_constant,
                                                    regimes, constant_regi=None, cols2regi=cols2regi,
                                                    names=name_x)
            yend2, name_yend = REGI.Regimes_Frame.__init__(self, yend,
                                                           regimes, constant_regi=None,
                                                           cols2regi=cols2regi, yend=True, names=name_yend)

            tsls = BaseTSLS(y=y, x=x, yend=yend2, q=q)
            self.k = tsls.z.shape[1]
            self.x = tsls.x
            self.yend, self.z = tsls.yend, tsls.z
            moments = _momentsGM_Error(w, tsls.u)
            lambda1 = optim_moments(moments)
            xs = get_spFilter(w, lambda1, x_constant)
            xs = REGI.Regimes_Frame.__init__(self, xs,
                                             regimes, constant_regi=None, cols2regi=cols2regi)[0]
            ys = get_spFilter(w, lambda1, y)
            yend_s = get_spFilter(w, lambda1, yend)
            yend_s = REGI.Regimes_Frame.__init__(self, yend_s,
                                                 regimes, constant_regi=None, cols2regi=cols2regi,

name_ye = ['end_' + str(i + 1) + "_" + str(r+1) for i in range(ky)]
                name_bigyend[r] = name_ye

        if name_bigq is None:
            name_bigq = {}
            for r in range(self.n_eq):
                ki = bigq[r].shape[1]
                name_i = ['inst_' + str(i + 1) + "_" + str(r+1) for i in range(ki)]
                name_bigq[r] = name_i

        if regimes is not None:
            self.constant_regi = 'many'
            self.cols2regi = 'all'
            self.regime_err_sep = False
            self.name_regimes = USER.set_name_ds(name_regimes)
            self.regimes_set = REGI._get_regimes_set(regimes)
            self.regimes = regimes
            cols2regi_dic = {}
            self.name_bigX,self.name_x_r,self.name_bigq,self.name_bigyend = {},{},{},{}

            for r in range(self.n_eq):
                self.name_x_r[r] = name_bigX[r] + name_bigyend[r]
                cols2regi_dic[r] = REGI.check_cols2regi(self.constant_regi, self.cols2regi, bigX[r], yend=bigyend[r], add_cons=False)
                USER.check_regimes(self.regimes_set, bigy[0].shape[0], bigX[r].shape[1])
                bigX[r], self.name_bigX[r] = REGI.Regimes_Frame.__init__(self, bigX[r],\
                 regimes, constant_regi=None, cols2regi=cols2regi_dic[r], names=name_bigX[r])
                bigq[r], self.name_bigq[r] = REGI.Regimes_Frame.__init__(self, bigq[r],\
                 regimes, constant_regi=None, cols2regi='all', names=name_bigq[r])
                bigyend[r], self.name_bigyend[r] = REGI.Regimes_Frame.__init__(self, bigyend[r],\
                 regimes, constant_regi=None, cols2regi=cols2regi_dic[r], yend=True, names=name_bigyend[r])
        else:
            self.name_bigX,self.name_bigq,self.name_bigyend = name_bigX,name_bigq,name_bigyend

from numpy import linalg as la
from pysal.lib.weights.spatial_lag import lag_spatial
from .utils import power_expansion, set_endog, iter_msg, sp_att
from .utils import get_A1_hom, get_A2_hom, get_A1_het, optim_moments
from .utils import get_spFilter, get_lags, _moments2eqs
from .utils import spdot, RegressionPropsY, set_warn
from .ols import BaseOLS
from .twosls import BaseTSLS
from .error_sp_hom import BaseGM_Error_Hom, BaseGM_Endog_Error_Hom, moments_hom, get_vc_hom, get_omega_hom, get_omega_hom_ols
from . import regimes as REGI
from . import user_output as USER
from . import summary_output as SUMMARY
from platform import system


class GM_Error_Hom_Regimes(RegressionPropsY, REGI.Regimes_Frame):

    '''
    GMM method for a spatial error model with homoskedasticity, with regimes, 
    results and diagnostics; based on Drukker et al. (2013) :cite:`Drukker2013`, following
    Anselin (2011) :cite:`Anselin2011`.

    Parameters
    ----------
    y            : array
                   nx1 array for dependent variable
    x            : array
                   Two dimensional array with n rows and one column for each
                   independent (exogenous) variable, excluding the constant
    regimes      : list
                   List of n values with the mapping of each
                   observation to a regime. Assumed to be aligned with 'x'.

# test on constancy of coefficients across equations
        if check_k(self.bigK):   # only for equal number of variables
            self.surchow = sur_chow(self.n_eq,self.bigK,self.b3SLS,self.varb)
        else:
            self.surchow = None

        #Listing of the results
        self.title = "THREE STAGE LEAST SQUARES (3SLS)"
        if regimes is not None:
            self.title += " - REGIMES"
            self.chow_regimes = {}
            varb_counter = 0
            for r in range(self.n_eq):
                counter_end = varb_counter+self.b3SLS[r].shape[0]
                self.chow_regimes[r] = REGI._chow_run(len(cols2regi_dic[r]),0,0,len(self.regimes_set),self.b3SLS[r],self.varb[varb_counter:counter_end,varb_counter:counter_end])
                varb_counter = counter_end     
            regimes=True   

        SUMMARY.SUR(reg=self, tsls=True, ml=False, nonspat_diag=nonspat_diag, regimes=regimes)

(counter * self.kr):((counter + 1) * self.kr)] = results[r].vm
            self.betas[
                (counter * self.kr):((counter + 1) * self.kr), ] = results[r].betas
            self.u[regi_ids[r], ] = results[r].u
            self.predy[regi_ids[r], ] = results[r].predy
            self.e_filtered[regi_ids[r], ] = results[r].e_filtered
            self.name_y += results[r].name_y
            self.name_x += results[r].name_x
            counter += 1
        self.chow = REGI.Chow(self)
        self.multi = results
        SUMMARY.GM_Error_Hom_multi(
            reg=self, multireg=self.multi, vm=vm, regimes=True)


class GM_Endog_Error_Hom_Regimes(RegressionPropsY, REGI.Regimes_Frame):

    '''
    GMM method for a spatial error model with homoskedasticity, regimes and
    endogenous variables.
    Based on Drukker et al. (2013) :cite:`Drukker2013`, following Anselin (2011)
    :cite:`Anselin2011`.

    Parameters
    ----------
    y            : array
                   nx1 array for dependent variable
    x            : array
                   Two dimensional array with n rows and one column for each
                   independent (exogenous) variable, excluding the constant
    yend         : array
                   Two dimensional array with n rows and one column for each

if not cores:
                results[r] = results_p[r]
            else:
                results[r] = results_p[r].get()

            self.vm[(counter * self.kr):((counter + 1) * self.kr),
                    (counter * self.kr):((counter + 1) * self.kr)] = results[r].vm
            self.betas[
                (counter * self.kr):((counter + 1) * self.kr), ] = results[r].betas
            self.u[regi_ids[r], ] = results[r].u
            self.predy[regi_ids[r], ] = results[r].predy
            self.e_filtered[regi_ids[r], ] = results[r].e_filtered
            self.name_y += results[r].name_y
            self.name_x += results[r].name_x
            counter += 1
        self.chow = REGI.Chow(self)
        self.multi = results
        SUMMARY.ML_Error_multi(
            reg=self, multireg=self.multi, vm=vm, spat_diag=spat_diag, regimes=True, w=w)

wu2 = (wu * wu).sum(0)
    wwu = w * wu
    uwwu = (u * wwu).sum(0)
    wwu2 = (wwu * wwu).sum(0)
    wwuwu = (wwu * wu).sum(0)
    trWtW = w.multiply(w).sum()
    moments = {}
    for r in range(u.shape[1]):        
        g = np.array([[u2[r], wu2[r], uwu[r]]]).T / n
        G = np.array(
            [[2 * uwu[r], -wu2[r], n], [2 * wwuwu[r], -wwu2[r], trWtW],
             [uwwu[r] + wu2[r], -wwuwu[r], 0.]]) / n
        moments[r] = [G, g]
    return moments

class SURerrorGM(BaseSURerrorGM, REGI.Regimes_Frame):
    """User class for SUR Error estimation by Maximum Likelihood
    
    Parameters
    ----------
    bigy         : dictionary with vectors of dependent variable, one for
                   each equation
    bigX         : dictionary with matrices of explanatory variables,
                   one for each equation
    w            : spatial weights object
    regimes      : list
                   List of n values with the mapping of each
                   observation to a regime. Assumed to be aligned with 'x'.
    nonspat_diag : boolean; flag for non-spatial diagnostics, default = False
    spat_diag    : boolean; flag for spatial diagnostics, default = False (to be implemented)
    vm           : boolean; flag for asymptotic variance for lambda and Sigma,
                   default = False (to be implemented)