From 07d019592aa26fcaa53695ea48fc02f00fa414cc Mon Sep 17 00:00:00 2001 From: Pedro Amaral Date: Wed, 3 Jul 2024 10:24:11 -0300 Subject: [PATCH] Minor adjustments to printouts and spatial impacts --- spreg/dgp.py | 4 ++-- spreg/ml_lag.py | 12 ++++++----- spreg/ml_lag_regimes.py | 17 ++++++++------- spreg/output.py | 21 +++++++++++-------- spreg/sputils.py | 42 ++++++++++++++++++++++++++++++++++++-- spreg/twosls_sp.py | 13 +++++++----- spreg/twosls_sp_regimes.py | 14 +++++++------ 7 files changed, 87 insertions(+), 36 deletions(-) diff --git a/spreg/dgp.py b/spreg/dgp.py index 7053e01f..36576d79 100644 --- a/spreg/dgp.py +++ b/spreg/dgp.py @@ -291,7 +291,7 @@ def make_wxg(wx,gamma): def dgp_errproc(u,w,lam=0.5,model='sar',imethod='power_exp'): """ - dgp_sar: generates pure spatial error process + dgp_errproc: generates pure spatial error process Arguments: ---------- @@ -303,7 +303,7 @@ def dgp_errproc(u,w,lam=0.5,model='sar',imethod='power_exp'): Returns: -------- - y : vector of observations following a pure SAR process + y : vector of observations following a spatial AR or MA process Examples -------- diff --git a/spreg/ml_lag.py b/spreg/ml_lag.py index 5d350635..32287e76 100755 --- a/spreg/ml_lag.py +++ b/spreg/ml_lag.py @@ -340,9 +340,11 @@ class ML_Lag(BaseML_Lag): tolerance criterion in mimimize_scalar function and inverse_product spat_diag : boolean If True, then compute Common Factor Hypothesis test when applicable - spat_impacts : boolean - If True, include average direct impact (ADI), average indirect impact (AII), - and average total impact (ATI) in summary results + spat_impacts : string + Include average direct impact (ADI), average indirect impact (AII), + and average total impact (ATI) in summary results. + Options are 'simple', 'full', 'power', or None. + See sputils.spmultiplier for more information. vm : boolean if True, include variance-covariance matrix in summary results @@ -594,7 +596,7 @@ def __init__( slx_lags=0, method="full", epsilon=0.0000001, - spat_impacts=True, + spat_impacts="simple", vm=False, spat_diag=True, name_y=None, @@ -642,7 +644,7 @@ def __init__( if spat_diag and slx_lags==1: diag_out = _spat_diag_out(self, w, 'yend', ml=True) if spat_impacts and slx_lags == 0: - impacts = _summary_impacts(self, _spmultiplier(w, self.rho)) + impacts = _summary_impacts(self, _spmultiplier(w, self.rho, method=spat_impacts), spat_impacts) try: diag_out += impacts except TypeError: diff --git a/spreg/ml_lag_regimes.py b/spreg/ml_lag_regimes.py index 93da512f..0eb990ca 100644 --- a/spreg/ml_lag_regimes.py +++ b/spreg/ml_lag_regimes.py @@ -66,9 +66,11 @@ class ML_Lag_Regimes(BaseML_Lag, REGI.Regimes_Frame): the spatial parameter is fixed accross regimes. spat_diag : boolean If True, then compute Common Factor Hypothesis test when applicable - spat_impacts : boolean - If True, include average direct impact (ADI), average indirect impact (AII), - and average total impact (ATI) in summary results + spat_impacts : string + Include average direct impact (ADI), average indirect impact (AII), + and average total impact (ATI) in summary results. + Options are 'simple', 'full', 'power', or None. + See sputils.spmultiplier for more information. cores : boolean Specifies if multiprocessing is to be used Default: no multiprocessing, cores = False @@ -318,7 +320,7 @@ def __init__( regime_lag_sep=False, cores=False, spat_diag=True, - spat_impacts=True, + spat_impacts="simple", vm=False, name_y=None, name_x=None, @@ -448,8 +450,9 @@ def __init__( diag_out = _spat_diag_out(self, w, 'yend', ml=True) else: self.title = ("MAXIMUM LIKELIHOOD SPATIAL LAG - REGIMES"+ " (METHOD = "+ method+ ")") + if spat_impacts and slx_lags == 0: - impacts = _summary_impacts(self, _spmultiplier(w, self.rho)) + impacts = _summary_impacts(self, _spmultiplier(w, self.rho, method=spat_impacts), spat_impacts, regimes=True) try: diag_out += impacts except TypeError: @@ -605,7 +608,7 @@ def ML_Lag_Regimes_Multi( if spat_diag and slx_lags == 1: results[r].other_mid += _spat_diag_out(results[r], None, 'yend', ml=True) if spat_impacts and slx_lags == 0: - results[r].other_mid += _summary_impacts(results[r], _spmultiplier(results[r].w, results[r].rho)) + results[r].other_mid += _summary_impacts(results[r], _spmultiplier(results[r].w, results[r].rho, method=spat_impacts), spat_impacts) counter += 1 self.multi = results self.chow = REGI.Chow(self) @@ -644,7 +647,7 @@ def _work( model.schwarz = DIAG.schwarz(reg=model) model.slx_lags = slx_lags model.w = w_r - model.rho = model.betas[-1] + model.rho = model.betas[-1] return model diff --git a/spreg/output.py b/spreg/output.py index 389fd12e..5c98f586 100644 --- a/spreg/output.py +++ b/spreg/output.py @@ -5,10 +5,10 @@ import textwrap as TW import numpy as np import pandas as pd -import math from . import diagnostics as diagnostics from . import diagnostics_tsls as diagnostics_tsls from . import diagnostics_sp as diagnostics_sp +from .sputils import _sp_effects __all__ = [] @@ -568,7 +568,7 @@ def _summary_iteration(reg): return txt -def _summary_impacts(reg, spmult): +def _summary_impacts(reg, spmult, spat_impacts, slx_lags=0, regimes=False): """ Spatial direct, indirect and total effects in spatial lag model. Uses multipliers computed by sputils._spmultipliers. @@ -577,21 +577,24 @@ def _summary_impacts(reg, spmult): ---------- reg: spreg regression object spmult: spatial multipliers as a dictionary + spat_impacts: spatial impacts method as string + slx_lags: int, number of spatial lags of X in the model + regimes: boolean, True if regimes model Returns ------- strings with direct, indirect and total effects """ - variables = reg.output[reg.output['var_type'].isin(['x', 'wx', 'yend']) & (reg.output.index != 0)] + variables = reg.output.query("var_type in ['x', 'yend'] and index != 0") + if regimes: + variables = variables[~variables['var_names'].str.endswith('_CONSTANT')] variables_index = variables.index - m1 = spmult['ati'] - btot = m1 * reg.betas[variables_index] - m2 = spmult['adi'] - bdir = m2 * reg.betas[variables_index] - m3 = spmult['aii'] - bind = m3 * reg.betas[variables_index] + + btot, bdir, bind = _sp_effects(reg, variables, spmult, slx_lags) + strSummary = "\nSPATIAL LAG MODEL IMPACTS\n" + strSummary += "Impacts computed using the '" + spat_impacts + "' method.\n" strSummary += " Variable Direct Indirect Total\n" for i in range(len(variables)): strSummary += "%20s %12.4f %12.4f %12.4f\n" % ( diff --git a/spreg/sputils.py b/spreg/sputils.py index 65cae18c..e0e7cece 100644 --- a/spreg/sputils.py +++ b/spreg/sputils.py @@ -307,8 +307,7 @@ def _spmultiplier(w, rho, method="simple", mtol=0.00000001): adii0 = np.sum(np.diag(invirw0)) multipliers["adi"] = adii0 / n elif method == "power": - wf = w.full()[0] - ws3 = SP.csr_array(wf) + ws3 = w.to_sparse(fmt='csr') rhop = rho ww = ws3 pow = 1 @@ -329,6 +328,45 @@ def _spmultiplier(w, rho, method="simple", mtol=0.00000001): multipliers["aii"] = multipliers["ati"] - multipliers["adi"] return (multipliers) +def _sp_effects(reg, variables, spmult, slx_lags=0): + """ + Calculate spatial lag, direct and indirect effects + + Attributes + ---------- + reg : regression object + variables : chunk of self.output with variables to calculate effects + spmult : dictionary with spatial multipliers + slx_lags : number of SLX lags + + Returns + ------- + btot : total effects + bdir : direct effects + bind : indirect effects + """ + variables_index = variables.index + m1 = spmult['ati'] + btot = m1 * reg.betas[variables_index] + m2 = spmult['adi'] + bdir = m2 * reg.betas[variables_index] + + # Assumes all SLX effects are indirect effects. Needs revision by LA. + if slx_lags > 0: + variables_wx = reg.output.query("var_type == 'wx'") + variables_wx_index = variables_wx.index + + chunk_size = len(variables) + for i in range(slx_lags): + start_idx = i * chunk_size + end_idx = start_idx + chunk_size + chunk_indices = variables_wx_index[start_idx:end_idx] + btot += m1 * reg.betas[chunk_indices] + bind = btot - bdir + else: + m3 = spmult['aii'] + bind = m3 * reg.betas[variables_index] + return btot, bdir, bind def _test(): import doctest diff --git a/spreg/twosls_sp.py b/spreg/twosls_sp.py index a58139bc..be57004c 100755 --- a/spreg/twosls_sp.py +++ b/spreg/twosls_sp.py @@ -238,9 +238,11 @@ class GM_Lag(BaseGM_Lag): If True, then use n-k to estimate sigma^2. If False, use n. spat_diag : boolean If True, then compute Anselin-Kelejian test and Common Factor Hypothesis test (if applicable) - spat_impacts : boolean - If True, include average direct impact (ADI), average indirect impact (AII), - and average total impact (ATI) in summary results + spat_impacts : string + Include average direct impact (ADI), average indirect impact (AII), + and average total impact (ATI) in summary results. + Options are 'simple', 'full', 'power', or None. + See sputils.spmultiplier for more information. vm : boolean If True, include variance-covariance matrix in summary results @@ -505,7 +507,7 @@ def __init__( gwk=None, sig2n_k=False, spat_diag=True, - spat_impacts=True, + spat_impacts="simple", vm=False, name_y=None, name_x=None, @@ -585,10 +587,11 @@ def __init__( self.output['regime'], self.output['equation'] = (0, 0) self.other_top = _spat_pseudo_r2(self) diag_out = None + if spat_diag: diag_out = _spat_diag_out(self, w, 'yend') if spat_impacts and slx_lags == 0: - impacts = _summary_impacts(self, _spmultiplier(w, self.rho)) + impacts = _summary_impacts(self, _spmultiplier(w, self.rho, method=spat_impacts), spat_impacts, slx_lags) try: diag_out += impacts except TypeError: diff --git a/spreg/twosls_sp_regimes.py b/spreg/twosls_sp_regimes.py index 185a0878..18a443dd 100755 --- a/spreg/twosls_sp_regimes.py +++ b/spreg/twosls_sp_regimes.py @@ -88,9 +88,11 @@ class GM_Lag_Regimes(TSLS_Regimes, REGI.Regimes_Frame): matrix must have ones along the main diagonal. sig2n_k : boolean If True, then use n-k to estimate sigma^2. If False, use n. - spat_impacts : boolean - If True, include average direct impact (ADI), average indirect impact (AII), - and average total impact (ATI) in summary results + spat_impacts : string + Include average direct impact (ADI), average indirect impact (AII), + and average total impact (ATI) in summary results. + Options are 'simple', 'full', 'power', or None. + See sputils.spmultiplier for more information. spat_diag : boolean If True, then compute Anselin-Kelejian test and Common Factor Hypothesis test (if applicable) vm : boolean @@ -455,7 +457,7 @@ def __init__( gwk=None, sig2n_k=False, spat_diag=True, - spat_impacts=True, + spat_impacts="simple", constant_regi="many", cols2regi="all", regime_lag_sep=False, @@ -643,7 +645,7 @@ def __init__( if spat_diag: diag_out = _spat_diag_out(self, w, 'yend') if spat_impacts and slx_lags == 0: - impacts = _summary_impacts(self, _spmultiplier(w, self.rho)) + impacts = _summary_impacts(self, _spmultiplier(w, self.rho, method=spat_impacts), spat_impacts, regimes=True) try: diag_out += impacts except TypeError: @@ -847,7 +849,7 @@ def GM_Lag_Regimes_Multi( if spat_diag: results[r].other_mid += _spat_diag_out(results[r], results[r].w, 'yend') if spat_impacts and slx_lags == 0: - results[r].other_mid += _summary_impacts(results[r], _spmultiplier(results[r].w, results[r].rho)) + results[r].other_mid += _summary_impacts(results[r], _spmultiplier(results[r].w, results[r].rho, method=spat_impacts), spat_impacts) counter += 1 self.multi = results if robust == "hac":