16_GMM_estimation_spatial_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "89145993",
   "metadata": {},
   "source": [
    "# GMM Estimation - Spatial Error Model\n",
    "\n",
    "### Luc Anselin\n",
    "\n",
    "### (revised 09/26/2024)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a85427f9",
   "metadata": {},
   "source": [
    "## Preliminaries\n",
    "\n",
    "This module covers the estimation of spatial error models by means of the general method of moments (GMM). As for the IV estimation of the spatial lag model, this allows for both exogenous and endogenous explanatory variables. The estimates for the spatial error parameter are obtained as the solution of a set of moment conditions. As of version 1.4 of *spreg*, all spatial error estimation is implemented through the `spreg.GMM_Error` command. This is essentially a wrapper around the original implementations of the various `GM_Error`, `GM_Endog_Error` and `GM_Combo` functions (that still work exactly as before) with a more simplified interface.\n",
    "Beyond the classic spatial error specification, more complex models can be estimated by including `slx_lags` (for additional WX) or `add_wy` (for inclusion of a spatially lagged dependent variable). This yields a range of higher order models, such as SLX Error, SARSAR (spatial lag with spatial autoregressive errors), and the generalized nested specification (GNS), i.e., a Spatial Durbin model with spatial autoregressive errors. The latter two models are treated in a separate notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9a2c9eb",
   "metadata": {},
   "source": [
    "### Modules Needed\n",
    "\n",
    "The main module continues to be *spreg* for spatial regression analysis. From this, `GMM_Error` is 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. "
   ]
  },
  {
   "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.weights import lag_spatial\n",
    "\n",
    "from spreg import GMM_Error"
   ]
  },
  {
   "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",
    "\n",
    "- from spreg:\n",
    "  - spreg.GMM_Error"
   ]
  },
  {
   "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": [
    "The data set and spatial weights are again 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",
    "\n",
    "To illustrate the methods, the same descriptive model is used as in the ML notebook. It 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 error 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",
    "\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'"
   ]
  },
  {
   "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 weights are row-standardized 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",
    "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": "a8a0c186-dd30-4f64-9a5b-d95bb562ebde",
   "metadata": {},
   "source": [
    "## GMM Estimation of the Error Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cc74112-a2eb-4efa-9116-d1e5197d20ed",
   "metadata": {},
   "source": [
    "### Principle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbed835b",
   "metadata": {},
   "source": [
    "As in the ML case, estimation of the spatial error model by means of GMM is based on spatially weighted least squares, where a consistent estimate for $\\lambda$ is used to create the spatially filtered dependent and explanatory variables, $(I - \\hat{\\lambda} W)y$ and $(I - \\hat{\\lambda} W)X$. If the only interest is in obtaining consistent estimates for $\\beta$, then any consistent estimate for $\\lambda$ will do. The only property that matters is consistency, since, unlike what holds in the spatial lag case, the precision of the $\\lambda$ estimate does not affect the precision of the $\\beta$ estimates.\n",
    "\n",
    "In two classic papers by Kelejian and Prucha (1998,1999), it was shown how a consistent estimate can be obtained as the solution of a set of moment equations formulated in terms of the regression residuals and their spatial lags. However, this did not provide an asymptotic variance for the spatial parameter, and thus did not allow for inference. In a series of later papers (Kelejian and Prucha 2010, Arraiz et al. 2010, 2013, Drukker et al 2013), the approach was extended to also allow for heteroskedasticity of unknown form and to include an asymptotic variance matrix.\n",
    "\n",
    "The technical details are quite complex and are outlined in Chapter 9 of Anselin and Rey (2014). The upshot is that there are three main estimation methods: the original Kelejian-Prucha generalized momens (GM) approach, the generalized method of moments (GMM) approach with heteroskedastic errors; and the GMM approach with homoskedastic errors. In practice, the GMM-heteroskedastic approach is greatly preferred. It is also the default used by `GMM_Error`.\n",
    "\n",
    "The GM and GMM approaches can be extended to models that include endogenous variables by means of spatially weighted two stage least squares (SW2SLS). This estimator uses the same expression as standard 2SLS, but replaces the $Z$ matrix of exogenous and endogenous regressors by its spatially filtered counterpart, $Z - \\hat{\\lambda} WZ$. The instrument matrix $Q$ is unaffected.\n",
    "\n",
    "In both cases, inference in based on two asymptotic variance matrices: one for the regression coefficients and one for the spatial error coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bfd7648-bdd1-4ebb-bcf5-4dc74782693d",
   "metadata": {},
   "source": [
    "### Implementation methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eec12f7",
   "metadata": {},
   "source": [
    "GMM estimation of the spatial error model is implemented in `spreg.GMM_Error`. This requires the standard regression arguments (i.e., at a minimum, `y`, `x` and `w`, as well as `yend` and `q` for the endogenous case). Three methods are implemented, specified by means of the `estimator` argument.\n",
    "\n",
    "The default is `estimator = \"het\"` for the heteroskedastic-robust GMM method. Other options are `estimator = \"hom\"`, for GMM with homoskedastic errors, and `estimator = \"kp98\"`, for the legacy Kelejian-Prucha GM approach.\n",
    "\n",
    "In addition, there are some more technical options for the GMM methods. In practice, these are rarely needed. \n",
    "\n",
    "For GMM-heteroskedastic, the extra arguments include `step1c`, `max_iter` and `epsilon`. The default is `step1c = False`, `max_iter = 1` and `epsilon = 0.00001`. When set to `True`, `step1c` carries out an additional estimation step for the autoregressive coefficient after the solution of the initial set of moment equations, as in Arraiz et. al (2010). The default follows the later paper by Arraiz et al. (2013) and skips this step. It is possible to iterate the procedure by using the new/updated residuals in additional rounds of estimation by setting `max_iter` to a larger value than `1`. Typically, this is not needed. When there are additional iterations, `epsilon` is used as a convergence criterion to stop the procedure.\n",
    "\n",
    "For GMM-homoskedastic, the option `A1` determines the exact manner in which the first set of moment equations is constructed. This pertains to very technical details regarding the trace of a matrix. The default is `A1 = \"hom_sc\"`, which applies a scaling factor suggested in Drukker et al. (2013). Other options are `A1 = \"hom\"` for no scaling, and `A1 = \"het\"` for a slightly different matrix expression (details are given on p. 217 of Anselin and Rey 2014). In practice, the options seldom make much of a difference and the default is fine. There are no `step1c` or `max_iter` options for this case.\n",
    "\n",
    "The legacy `kp98` method has no special options."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71307ef9-cc53-46a6-866e-5bd521e4503e",
   "metadata": {},
   "source": [
    "#### Estimator `het`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "799e1559",
   "metadata": {},
   "source": [
    "GMM-heteroskedastic is the default, so the `estimator` argument does not need to be specified. The first illustration is for all default settings with only exogenous regressors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1975d6b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "err1 = GMM_Error(y,x,w=wq,\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(err1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "debc6278",
   "metadata": {},
   "source": [
    "The estimate of $\\lambda$, 0.4677, is almost identical to that obtained by means of ML (0.451). The other characteristics of the model are very similar as well, with a pseudo $R^2$ of 0.634 (compared to 0.633). The interpretation of the other features of the regression object are the same as for the ML estimation, and will not be repeated here.\n",
    "\n",
    "Again, there is only one type of predicted value, contained in **predy**. There are two types of residuals, the classic residual, $u = y - X \\hat{\\beta}$, stored in **u**, and the spatially filtered residuals, $u - \\lambda W u$, stored in **e_filtered**. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "418df371",
   "metadata": {},
   "source": [
    "To illustrate the (minor) effect of the additional arguments, `step1c` is next set to `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0a22c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "err2 = GMM_Error(y,x,w=wq,\n",
    "                     step1c=True,\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(err2.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "343e01de",
   "metadata": {},
   "source": [
    "The estimate for $\\lambda$ is only marginally different, at 0.4722, with a slightly smaller standard error, but the fit is unaffected (in fact, in terms of pseudo $R^2$, it is slightly worse than before, 0.6337 vs. 0.6340).\n",
    "\n",
    "With `max_iter = 10`, the differences are again marginal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30acd3e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "err3 = GMM_Error(y,x,w=wq,\n",
    "                     max_iter=10,\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(err3.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d3b5229-f169-4d15-aafe-416f325183a5",
   "metadata": {},
   "source": [
    "#### Estimator `hom`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ba06a23-111a-430e-bd60-ce54ee60a660",
   "metadata": {},
   "source": [
    "In practice, ignoring heteroskedasticiy is typically not a good idea for cross-sectional regressions. The `hom` estimator is included here for the sake of completeness, but should usually not be considered. Only the default settings are illustrated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ec4ab7-5eed-4c14-b046-2f9efff13642",
   "metadata": {},
   "outputs": [],
   "source": [
    "err4 = GMM_Error(y,x,w=wq,\n",
    "                     estimator=\"hom\",\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(err4.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff1ea981-669d-4ba0-a200-6f41c17717f0",
   "metadata": {},
   "source": [
    "The spatial coefficient is again slightly different, but the effect on the $\\beta$ coefficient is marginal. The main differences are in the results for the standard errors. As mentioned, ignoring heteroskedasticity may yield an overly optimistic impression of precision."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0439c7c-771a-446f-8677-51570d7f8dea",
   "metadata": {},
   "source": [
    "#### Estimator `kp98`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "846b5ac1-e4b0-4b99-8a56-7663403935d1",
   "metadata": {},
   "source": [
    "The final estimator is the legacy `kp98` method, which is only included for completeness. Since it does not provide inference for the spatial parameter, it is otherwise not recommended for use in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e485647-951e-4b62-80be-b1bb505b5027",
   "metadata": {},
   "outputs": [],
   "source": [
    "err5 = GMM_Error(y,x,w=wq,\n",
    "                     estimator=\"kp98\",\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(err5.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34682841-1194-4f32-a14d-1aea4ce6b9e3",
   "metadata": {},
   "source": [
    "Note how the output now does not list standard errors, z-statistics and p-value for the $\\lambda$ estimate. The regression estimates are essentially the same."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "708ffa42-763f-4ee2-acf0-2ab079795cfd",
   "metadata": {},
   "source": [
    "## GMM Estimation of the SLX-Error Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a27b8d3-13cb-4e6d-aef2-549f19c6dc6f",
   "metadata": {},
   "source": [
    "As for ML, GMM estimation of the SLX-Error model is a special case of `spreg.GMM_Error`, with the additional argument of `slx_lags=1` (or a larger value). Everything else remains the same. More specifically, the three estimators of `het`, `hom` and `kp98` are again available. Only the default `het` is considered here. The results show only minor differences for the other methods. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d5bc8e4-3683-4bfd-a904-f944cc69a71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "slxerr1 = GMM_Error(y,x,w=wq,slx_lags=1,\n",
    "                    name_w=w_name,name_ds=ds_name)\n",
    "print(slxerr1.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b650f7-2272-4f1a-8d2b-e6c007f45d3a",
   "metadata": {},
   "source": [
    "In general, as before, there is little support for including all the WX terms after controlling for spatial error autocorrelation. Only the coefficient of **W_HIS_ct** is significant, but, as before, the value of the estimate is larger than that for the unlagged regressor.\n",
    "\n",
    "The spatial autoregressive coefficient, 0.422, is highly significant.\n",
    "\n",
    "As usual, further refinements of the model specification can be carried out by eliminating some lag terms by means of `slx_vars`, as in the standard SLX model. This is not further pursued here."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba39aea",
   "metadata": {},
   "source": [
    "## Exogenous and Endogenous Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e9b6a82",
   "metadata": {},
   "source": [
    "Additional endogenous variables and associated instruments are added in the standard way by including `yend` and `q` as arguments. Note that in the example, the regression matrix is now set to `xe`, for only the exogenous variables. Everything else is the same as before. Only the default `estimator = \"het\"` is illustrated, with its default settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7a9c39b",
   "metadata": {},
   "outputs": [],
   "source": [
    "enderr = GMM_Error(y,xe,w=wq,yend=yend,q=q,\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(enderr.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b05d3b36",
   "metadata": {},
   "source": [
    "The estimate for $\\lambda$ is again in the same general ballpark, but **HIS_ct** is no longer significant. As usual, the endogenous variable is listed as **Instrumented** as well as the **Instruments**.\n",
    "\n",
    "Finally, the SLX model can be estimated with additional endogenous variables, by including `yend` and `q` with `slx_lags=1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45855356",
   "metadata": {},
   "outputs": [],
   "source": [
    "slxenderr = GMM_Error(y,xe,w=wq,yend=yend,q=q,\n",
    "                      slx_lags=1,\n",
    "                     name_w=w_name,name_ds=ds_name)\n",
    "print(slxenderr.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6973349",
   "metadata": {},
   "source": [
    "In this example, none of the WX terms end up being significant. In addition, of the original regressors, only **EP_LIMENG** remains as significant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b086a22",
   "metadata": {},
   "source": [
    "## Practice\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, in the same way as for the other models."
   ]
  }
 ],
 "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
}