14_IV_estimation_spatial_lag.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "89145993",
   "metadata": {},
   "source": [
    "# Instrumental Variables Estimation - Spatial Lag Model\n",
    "\n",
    "### Luc Anselin\n",
    "\n",
    "### (revised 09/19/2024)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a85427f9",
   "metadata": {},
   "source": [
    "## Preliminaries\n",
    "\n",
    "An alternative to maximum likelihood is to tackle the endogeneity of the spatially lagged dependent variable by means of instrumental variables (IV) estimation. This is implemented in the `spreg.GM_Lag` class. As before, the estimation of the Spatial Durbin model is achieved through the inclusion of the `slx_lags` argument. \n",
    "\n",
    "Distinct from what is possible for maximum likelihood, other endogenous variables can be included as well, using the familiar **yend** (for the endogenous variables) and **q** (for the instruments) arguments.\n",
    "\n",
    "The treatment in this notebook will focus on the specific properties of the IV estimation. Generic properties of the spatial lag model, such as the different predicted values, residuals, and the impact measures will not be treated in detail again. Technical aspects pertaining to these issues are covered in the maximum likelihood notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9a2c9eb",
   "metadata": {},
   "source": [
    "### Modules Needed\n",
    "\n",
    "As before, the main module is *spreg* for spatial regression analysis. From this, `OLS` and `GM_Lag` are imported. In addition, the utilities in *libpysal* (to open spatial weights and access the sample data set), *pandas* and *geopandas* are needed. All of these rely on *numpy* as a dependency. \n",
    "\n",
    "The usual *numpy* `set_printoptions` is included as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5ac490b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import os\n",
    "os.environ['USE_PYGEOS'] = '0'\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from libpysal.io import open\n",
    "from libpysal.examples import get_path\n",
    "from libpysal import weights\n",
    "\n",
    "from spreg import OLS, GM_Lag\n",
    "np.set_printoptions(legacy=\"1.25\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee18820",
   "metadata": {},
   "source": [
    "### Functions Used\n",
    "\n",
    "- from pandas/geopandas:\n",
    "  - read_file\n",
    "  \n",
    "- from libpysal:\n",
    "  - io.open\n",
    "  - examples.get_path\n",
    "  - weights.distance.Kernel\n",
    " \n",
    "- from spreg:\n",
    "  - spreg.OLS\n",
    "  - spreg.GM_Lag"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47934d1b-7587-4ddf-be38-83904eede8e8",
   "metadata": {},
   "source": [
    "### Variable definition and data input"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfb04af4-4521-4ac0-8e55-d3b92f54a403",
   "metadata": {},
   "source": [
    "As in the maximum likelihood notebook, the data set and spatial weights are from the **chicagoSDOH** sample data set:\n",
    "\n",
    "- **Chi-SDOH.shp,shx,dbf,prj**: socio-economic indicators of health for 2014 in 791 Chicago tracts\n",
    "- **Chi-SDOH_q.gal**: queen contiguity weights\n",
    "- **Chi-SDOH_k10tri.kwt**: triangular kernel weights based on a variable bandwidth with 10 nearest neighbors from `GeoDa`\n",
    "\n",
    "The same descriptive model is used that relates the rate of uninsured households in a tract(for health insurance, **EP_UNINSUR**) to the lack of high school education (**EP_NOHSDP**), the economic deprivation index (**HIS_ct**), limited command of English (**EP_LIMENG**) and the lack of access to a vehicle (**EP_NOVEH**). This is purely illustrative of a spatial lag specification and does not have a particular theoretical or policy motivation.\n",
    "\n",
    "In an alternative specification, **HIS_ct** is considered to be endogenous, with, as before, **COORD_X** and **COORD_Y** as instruments.\n",
    "\n",
    "The file names and variable names are set in the usual manner. Any customization for different data sets/weights and different variables should be specified in this top cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "14a1e98e",
   "metadata": {},
   "outputs": [],
   "source": [
    "infileshp = get_path(\"Chi-SDOH.shp\")     # input shape file with data\n",
    "infileq = get_path(\"Chi-SDOH_q.gal\")     # queen contiguity weights created with GeoDa\n",
    "infilek = get_path(\"Chi-SDOH_k10tri.kwt\") # triangular kernel weights\n",
    "\n",
    "y_name = 'EP_UNINSUR'\n",
    "x_names = ['EP_NOHSDP','HIS_ct','EP_LIMENG','EP_NOVEH']\n",
    "xe_names = ['EP_NOHSDP','EP_LIMENG','EP_NOVEH']\n",
    "yend_names = ['HIS_ct']\n",
    "q_names = ['COORD_X','COORD_Y']\n",
    "ds_name = 'Chi-SDOH'\n",
    "w_name = 'Chi-SDOH_q'\n",
    "wk_name = 'Chi-SDOH_k10tri'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9514984-a22e-4d5b-a993-f40771a290a0",
   "metadata": {},
   "source": [
    "The `read_file` and `open` functions are used to access the sample data set and contiguity weights. The contiguity weights are row-standardized, the `class` of the kernel weights adjusted, and the data frames for the dependent and explanatory variables are constructed. As before, this functionality is agnostic to the actual data sets and variables used, since it relies on the specification given in the initial block above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3248d23d-2247-40bf-a3a5-3ba21b477aa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs = gpd.read_file(infileshp)\n",
    "wq =  open(infileq).read()    \n",
    "wq.transform = 'r'    # row-transform the weights\n",
    "wk = open(infilek).read()\n",
    "wk.__class__ = weights.distance.Kernel\n",
    "y = dfs[y_name]\n",
    "x = dfs[x_names]\n",
    "yend = dfs[yend_names]\n",
    "xe = dfs[xe_names]\n",
    "q = dfs[q_names]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aea09bf",
   "metadata": {},
   "source": [
    "## OLS and SLX with Spatial Diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb5eb774-814b-467e-ad51-11d21d3cee50",
   "metadata": {},
   "source": [
    "For easy reference, standard OLS and SLX regressions with spatial diagnostics are repeated here to provide a point of reference. The results are identical to those reported in the maximum likelihood notebook\n",
    "\n",
    "Refer to the specific OLS notebook for further details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb2b8ce-6ae8-4a15-a6c7-7ec23700dfc7",
   "metadata": {},
   "source": [
    "### OLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "379724ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "ols1 = OLS(y,x,w=wq,spat_diag=True,moran=True,\n",
    "                 name_w=w_name,name_ds=ds_name)\n",
    "print(ols1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d0e0f74-5877-4497-8e7b-79394b6701fd",
   "metadata": {},
   "source": [
    "The specification achieves an acceptable $R^2$ of about 0.63 and all coefficients are positive and highly significant.\n",
    "\n",
    "The non-spatial diagnostics suggest non-normality as well as a hight degree of heteroskedasticity. There is no problem with multicollinearity.\n",
    "\n",
    "The spatial diagnostics against the SARERR alternatives show very significant LM-Lag and LM-Error, but of the two robust tests, only RLM-Lag is highly significant (RLM-Error only at p < 0.03). Hence, there is a strong indication that a Lag rather than an Error alternative may be appropriate. While the joint LM test is also highly significant, this is likely due to a strong one-sided (Lag) alternative.\n",
    "\n",
    "Interestingly, the diagnostics against a spatial Durbin alternative strongly support the latter as well. Both LM tests and their robust forms are highly significant, and so is the joint test. Moreover, the value for the robust forms of the test is smaller than the original, which is the expected behavior (although not always reflected in empirical practice).\n",
    "\n",
    "In sum, in addition to a spatial Lag model as an alternative, the spatial Durbin specification deserves consideration as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfc5d4e-a1b3-4b07-8120-2def922d8ef6",
   "metadata": {},
   "source": [
    "### SLX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53f9b90d-7f6d-4b3e-9e99-83430957fe3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "slx1 = OLS(y,x,w=wq,slx_lags=1,spat_diag=True,moran=True,\n",
    "                 name_w=w_name,name_ds=ds_name)\n",
    "print(slx1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ee2be15-2067-4caf-b02b-aca6b46118c5",
   "metadata": {},
   "source": [
    "Relative to the classic regression model, the fit improves slightly, but the constant, **EP_NOHSDP** and **HIS_CT** become non-significant at p = 0.01 (they are marginally signifcant at p=0.05). All but one coefficient of the SLX terms are significant (**W_EP_NOVEH** is not). The signs and magnitudes of the SLX coefficients relative to their unlagged counterparts remain a bit confusing. Only for **EP_LIMENG** and **W_EP_LIMENG** are they the same, with the lag coefficient smaller than the unlagged one, in accordance with Tobler's law. The coefficient for **W_HIS_ct** is significant and larger than that of **HIS_ct**, while the latter is not significant at p = 0.01. In other words, the interpretation of these results in terms of distance decay and Tobler's law may be a bit problematic.\n",
    "\n",
    "In terms of diagnostics, there is a slight problem with multicollinearity (often the case in SLX specifications), strong non-normality and evidence of heteroskedasticity. Moran's I is significant, as are both LM-tests, but neither of the robust forms is significant. Based on the relative magnitudes of the test statistics, there is a slight indication of a possible Lag alternative, i.e., a spatial Durbin specification. However, this indication is not as strong as that provided by the LM-SDM test statistics in the classic specification."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8a0c186-dd30-4f64-9a5b-d95bb562ebde",
   "metadata": {},
   "source": [
    "## IV Estimation of the Lag Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cc74112-a2eb-4efa-9116-d1e5197d20ed",
   "metadata": {},
   "source": [
    "### Principle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0385d805-bfb4-4e79-b0ed-0fc59c811926",
   "metadata": {},
   "source": [
    "The lag model $y = \\rho Wy + X\\beta + u$ can also be written as $y = Z\\delta + u$, where $Z = [ X, Wy ]$ and the coefficient vector is rearranged as $\\delta = [\\beta, \\rho ]$. This setup is the same as that considered in the treatment of endogeneity by means of 2SLS estimation.\n",
    "\n",
    "The spatially lagged dependent variable is indeed endogenous. Based on the reduced form, \n",
    "$$y = (I - \\rho W)^{-1} X\\beta + (I - \\rho W)^{-1}u,$$ \n",
    "$Wy$ follows as:\n",
    "$$Wy = W(I - \\rho W)^{-1} X\\beta + W(I - \\rho W)^{-1}u.$$\n",
    "As a result, after some algebra, $E[(Wy)'u] = tr[W(I - \\rho W)'^{-1}].E[uu'] \\neq 0$. This follows because of the presence of non-zero diagonal terms in the inverse matrix (see the discussion of impacts in the maximum likelihood notebook). Consequently, one of the fundamental assumptions of OLS estimation is violated, and Wy must be treated as an endogenous variable.\n",
    "\n",
    "This is carried out by means of 2SLS estimation, whereby the main question becomes the choice of good (optimal) instruments for Wy. The conditional expectation of Wy, given X, provides the basis for this:\n",
    "$$E[Wy | X] = W(I - \\rho W)^{-1} X\\beta = W(I + \\rho W + \\rho^2 W + \\dots) X \\beta.$$\n",
    "As a result, the spatially lagged explanatory variables $WX, W^2X, \\dots$ suggest themselves as instruments. Since the $X$ are assumed to be uncorrelated with the error terms, so are the $WX$. Also, from the conditional expectation, it follows that they are related to the endogenous Wy, which satisfies a second requirement for instruments. Note that WX is *not* applied to the constant term, since this would result in perfect multicollinearity.\n",
    "\n",
    "Estimation then follows as a straightforward application of 2SLS, referred to as *spatial* 2SLS, or S2SLS. The estimator, variance-covariance matrix, and robust standard errors are the same as for the classis 2SLS (see the 2SLS notebook for details).\n",
    "\n",
    "In a straightforward manner, additional endogenous variables with associated instruments can be incorporated as well. One question remains whether the additional instruments should be lagged or not. In `spreg.GM_Lag`, this is handled by the `lag_q` option (the default is `lag_q = True`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bfd7648-bdd1-4ebb-bcf5-4dc74782693d",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ada2a0d-81fc-4b07-8ba9-d7e1bc101b1a",
   "metadata": {},
   "source": [
    "IV estimation of the spatial lag model is carried out by means of `spreg.GM_Lag`. This is a customized implementation of the two stage least squares estimation, with the instruments for the spatially lagged dependent variable computed internally. As mentioned, additional endogenous variables can be specified as well.\n",
    "\n",
    "The default setup requires the dependent variable, `y`, the explanatory variables (without a constant term), `x`, and the spatial weights `w`. The instruments are the spatially lagged explanatory variables, WX. They do not need to be specified separately. The order of spatial weights used as instruments is set by means of `w_lags` (the default is `w_lags = 1`).\n",
    "\n",
    "As is customary, the main results are listed using the `summary` method. For now, no impacts are listed by setting `spat_impacts = None` (note that `spat_impacts = \"simple\"` is the default). The AK-test for remaining residual spatial autocorrelation is included by default (i.e., `spat_diag = True`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8887e6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag1 = GM_Lag(y,x,w=wq,\n",
    "              name_w=w_name,name_ds=ds_name,\n",
    "              spat_impacts = None)\n",
    "print(lag1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee6e4e6",
   "metadata": {},
   "source": [
    "In this example, the estimates are very similar to the ML results (see the ML notebook). This is not always the case, but it is encouraging when it is (in the sense of not having other misspecification issues). The autoregressive coefficient of 0.377 compares to 0.392 in the ML case and is highly significant. As in the ML case, the coefficients of the constant term and of **EP_NOHSDP** become non-significant. The other regression coefficients are of similar magnitudes as in the ML case and thus also generally smaller than the corresponding OLS estimates.\n",
    "\n",
    "The results show the **Instrumented** variable as **W_EP_UNINSUR** and the **Instruments** as the spatial lags of the explanatory variables.\n",
    "\n",
    "As in the ML case, there are two measures of fit, the **Pseudo R-squared** (based on the *naive* residuals) and the **Spatial Pseudo R-squared** (based on the predicted values from the reduced form). The results are essentially the same as for the ML estimates, e.g., respectively 0.684 (vs. 0.685) and 0.647 (vs. 0.647). In constrast to the ML results, there is no Likelihood, AIC or SC.\n",
    "\n",
    "Note that with S2SLS estimation, it is possible for the spatial autoregressive parameter to take on a value larger than one. Unlike maximum likelihood estimation, where the parameter space is constrained in the optimization routine, there is no such constraint in S2SLS. When the parameter estimate is outside the bounds, a warning is given and any properties based on the reduced form predicted values will not be computed. This includes **predy_e** and **e_pred**, as well as the **Spatial Pseudo R-squared** and the spatial impacts (for an example, see below).\n",
    "\n",
    "The AK test shows no evidence of remaining residual spatial autocorrelation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14ac009",
   "metadata": {},
   "source": [
    "### Instrument lag order"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1aea751",
   "metadata": {},
   "source": [
    "The order of the spatial lags used for the instruments is determined by the `w_lags` argument. By default, this is set to 1. The use of higher order lags results in greater precision, but often at the cost of an increase in multicollinearity. In practice, using lags larger than 2 is not recommended.\n",
    "\n",
    "The effect of setting `w_lags = 2` is illustrated next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc994370",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag2 = GM_Lag(y,x,w=wq,\n",
    "              name_w=w_name,name_ds=ds_name,\n",
    "              w_lags = 2,\n",
    "              spat_impacts = None)\n",
    "print(lag2.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb8126b",
   "metadata": {},
   "source": [
    "The additional instruments are listed below the parameter estimates with the prefix **W2**. In this example, the effect of including the second order lags is minimal. The autoregressive coefficient becomes 0.381 (vs. 0.377) and the standard errors are marginally smaller. The fit is unchanged."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f747b04",
   "metadata": {},
   "source": [
    "### Predicted values and residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1306bd65",
   "metadata": {},
   "source": [
    "The treatment of predicted values and residuals is the same as for maximum likelihood estimation. The *naive* results are stored in the regression object attributes **predy** and **u**, and the values that are based on the reduced form prediction are in **predy_e** and **e_pred**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8a6d92f",
   "metadata": {},
   "source": [
    "## Robust Standard Errors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a1bb5e8",
   "metadata": {},
   "source": [
    "As in the standard 2SLS estimation, it is possible to obtain robust standard errors by means of the `robust = \"white\"` and `robust = \"hac\"` options. Given the common presence of heteroskedastic errors in cross-sectional regression, this is highly recommended. The results are given below. As before, the `hac` option requires that a kernel spatial weights object is specified as the argument `gwk`, with, optionally, `name_gwk`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e17afc1",
   "metadata": {},
   "source": [
    "### White standard errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "791f7485",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag3 = GM_Lag(y,x,w=wq,\n",
    "              name_w=w_name,name_ds=ds_name,\n",
    "              robust = 'white',\n",
    "              spat_impacts = None)\n",
    "print(lag3.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61ac8385",
   "metadata": {},
   "source": [
    "The estimates are the same as before, but the standard errors are slightly larger. For example, for the spatial autoregressive coefficient, the standard error becomes 0.07377, compared to 0.06455. The overall impact is minimal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a174f964",
   "metadata": {},
   "source": [
    "### HAC standard errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b422888e",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag4 = GM_Lag(y,x,w=wq,\n",
    "              name_w=w_name,name_ds=ds_name,\n",
    "              robust = 'hac',gwk=wk,name_gwk=wk_name,\n",
    "              spat_impacts = None)\n",
    "print(lag4.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7697c10",
   "metadata": {},
   "source": [
    "In this case, the HAC standard errors are very similar to the heteroskedasticity-robust ones, suggesting the main source of misspecification comes from the latter. This is in accordance with the lack of significance of the AK-test. The standard error for the spatial autoregressive coefficient, 0.07238, is even slightly smaller than the White standard error (0.07337)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9f084bf-276f-44a7-863b-41835831d7ee",
   "metadata": {},
   "source": [
    "## Spatial Multipliers - Impacts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6ed0760-3ba2-4f0e-a06d-a295c13e29b4",
   "metadata": {},
   "source": [
    "Similar to what holds for ML estimation, impact measures are computed for spatial lag models estimated by means of instrumental variables. This is implemented through the `spat_impacts` argument. The default setting is `spat_impacts = \"simple\"` for the Kim et al approach (see the ML notebook for technical details). As before, the other options are `full`, `power`, `all` or ``None`.\n",
    "\n",
    "In the example below, the argument is set as `spat_impacts = [\"simple\",\"full\"]`. Only classic standard errors are considered. Since the multipliers and impacts only depend on the coefficient estimates, the type of standard error is immaterial.\n",
    "\n",
    "Note that the reported impacts are only *average* effects. See the spatial multipliers notebook for a more extensive analysis of the associated spatial pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cb376e5-d8b2-4aaf-aad9-e3c659613419",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag5 = GM_Lag(y,x,w=wq,\n",
    "                     name_w=w_name,name_ds=ds_name,\n",
    "                     spat_impacts=['simple','full'])\n",
    "print(lag5.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b392f8",
   "metadata": {},
   "source": [
    "The interpretation of the direct, indirect and total effects is the same as before. Since the coefficient estimates were very similar to the results obtained for ML, the impacts are similar as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f83c898",
   "metadata": {},
   "source": [
    "## Additional Endogenous Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5ec68eb",
   "metadata": {},
   "source": [
    "Additional endogenous variables are included by means of the **yend** and **q** arguments, in the same way as for classic 2SLS estimation. The only relevant option in this regard is whether the instruments should be lagged as well. The default is to include their spatial lags, through `lag_q = True`. In practice, there is no good reason not to lag them, since they provide additional information. For clarity, the `lag_q` argument is included, even though it is not needed, since it is the default.\n",
    "\n",
    "Note that here `x = xe`, which contains only the exogenous variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b87f919",
   "metadata": {},
   "outputs": [],
   "source": [
    "lag6 = GM_Lag(y,x=xe,w=wq,yend=yend,q=q,\n",
    "                     lag_q = True,\n",
    "                     name_w=w_name,name_ds=ds_name,\n",
    "                     spat_impacts=['simple','full'])\n",
    "print(lag6.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0d983e4",
   "metadata": {},
   "source": [
    "The output listing now includes **HIS_ct** as **Instrumented**, and both **COORD_X** and **COORD_Y**, as well as their spatial lags among the **Instruments**.\n",
    "\n",
    "The impact on the estimates is relatively minor, although the spatial autoregressive coefficient decreases slightly to 0.32995 (compared to 0.37698). The treatment of **HIS_ct** as endogenous results in its coefficient becoming only marginally significant (p=0.02). The model impacts reflect the slightly different coefficient estimates. \n",
    "\n",
    "Note that impacts are only computed for the exogenous variables. In a strict sense, the additional endogenous variables should not be part of the reduced form. Moreover, since they are not determined in a fully simultaneous equation system, there is no practical way to include other exogenous variables from such a system. Therefore, the multiplier effect of additional endogenous variables is not considered."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "708ffa42-763f-4ee2-acf0-2ab079795cfd",
   "metadata": {},
   "source": [
    "## IV Estimation of Spatial Durbin Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a27b8d3-13cb-4e6d-aef2-549f19c6dc6f",
   "metadata": {},
   "source": [
    "IV estimation of the Spatial Durbin model is a special case of `spreg.GM_Lag`, with the additional argument of `slx_lags=1` (or a larger value). Everything else remains the same. \n",
    "\n",
    "In the example, `spat_impacts` is set to `[\"simple\",\"full\"]` (note, the default setting remains `spat_impacts = \"simple\"`).\n",
    "\n",
    "Another default setting is `spat_diag = True`, which yields the results for the both the AK test and the Common Factor Hypothesis test. To avoid these tests, `spat_diag` must be set to `False`. In the illustration, both arguments are listed explicitly for clarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d5bc8e4-3683-4bfd-a904-f944cc69a71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "spdur1 = GM_Lag(y,x,w=wq,slx_lags=1,\n",
    "                    name_w=w_name,name_ds=ds_name,\n",
    "                    spat_impacts = ['simple','full'],\n",
    "                    spat_diag = True)\n",
    "print(spdur1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b650f7-2272-4f1a-8d2b-e6c007f45d3a",
   "metadata": {},
   "source": [
    "Note how the list of **Instruments** now consists of the second order spatial lags, even though `w_lags = 1`. This is because the spatial Durbin model already includes the first order lags as explanatory variables. Using first order lags as instruments would result in perfect multicollinearity. This is detected internally in `GM_Lag`, so that the computation of the lagged instruments is adjusted accordingly.\n",
    "\n",
    "As in the case of ML estimation, the inclusion of the lagged explanatory variables has quite an effect, even though only **W_EP_NOHSDP** turns out to be significant (for ML, **W_HIS_ct** was significant as well). Also, the spatial autoregressive coefficient (0.509) is no longer significant. Three of the four lag coefficients (including the non-significant ones) have a negative sign, the opposite of the original regression coefficients. This raises the suspicion that the proper specification may be a spatial error model. The results of the Common Factory Hypothesis Test bear this out. At 4.083, it is *not* significant, which means that the common factor constraint can *not* be rejected. This contrasts greatly with the result for ML estimation, where the constraint was strongly rejected. However, given the lack of significance of the spatial lag terms, this needs to be interpreted with caution. A better strategy would be to consider further refinements of the model specification by eliminating some lag terms by means of `slx_vars`, as in the SLX model. This is not further pursued here.\n",
    "\n",
    "The impacts are computed as before. Also, robust standard errors can be implemented in the usual fashion."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b180e95a",
   "metadata": {},
   "source": [
    "### Spatial Durbin model with endogenous variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a66754f",
   "metadata": {},
   "source": [
    "Additional endogenous variables can be included in a spatial Durbin specification in the same way as in the standard spatial lag model. With `spat_diag = True`, the default, only the AK test is produced. Since there is no spatial lag for the endogenous variable, a common factor test is not meaningful.\n",
    "\n",
    "The example below uses the same specification as in the spatial lag model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e70ee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "spdur2 = GM_Lag(y,x=xe,w=wq,yend=yend,q=q,\n",
    "                     lag_q = True,\n",
    "                     slx_lags=1,\n",
    "                     name_w=w_name,name_ds=ds_name,\n",
    "                     spat_impacts=['simple','full'],\n",
    "                     spat_diag = True)\n",
    "print(spdur2.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d0d987f",
   "metadata": {},
   "source": [
    "In this example, the estimate for the spatial autoregressive coefficient turns out to be larger than one (1.01654), which is outside the accepted parameter space. This raises a warning and precludes the computation of the spatial pseudo R-squared and the model impacts.\n",
    "\n",
    "When this happens, further refinement of the model and/or estimation is required. There are several options, such as using higher lags for the instruments, using different/more instruments for the additional endogenous variable, dropping some lag terms (using `slx_vars`), or changing the original specification altogether. This is left as an exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b086a22",
   "metadata": {},
   "source": [
    "## Practice\n",
    "\n",
    "IV estimation of the spatial lag model is sensitive to several aspects of the model specification. However, in large(r) samples, it generally yields more robust results than maximum likelihood, especially when using robust standard errors.\n",
    "\n",
    "As practice, different model specifications could be considered, including adding additional explanatory variables, selectively removing some lag terms, and using different spatial weights. Make sure to carefully consider the interpretation of the estimated coefficients and associated direct and indirect effects, as well as the robustness of the standard errors."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}