diff --git a/spreg/dgp.py b/spreg/dgp.py
index 7053e01..36576d7 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 5d35063..32287e7 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 93da512..0eb990c 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 389fd12..5c98f58 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 65cae18..e0e7cec 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 a58139b..be57004 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 185a087..18a443d 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":