{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Instrumental-Variables Regression\n",
    "These examples follow those in **Chapter 6** of _Microeconometrics Using Stata_ by Cameron & Trivedi.\n",
    "\n",
    "The first step is to import the main estimator for linear IV models:\n",
    "\n",
    "* `IV2SLS` - standard two-stage least squares\n",
    "* `IVLIML` - Limited information maximum likelihood and k-class estimators\n",
    "* `IVGMM` - Generalized method of moment estimation\n",
    "* `IVGMMCUE` - Generalized method of moment estimation using continuously updating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels import IV2SLS, IVGMM, IVGMMCUE, IVLIML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing data\n",
    "The data uses comes from the Medical Expenditure Panel Survey (MEPS) and includes data on out-of-pocket drug expenditure (in logs), individual characteristics, whether an individual was insured through an employer or union (a likely endogenous variable), and some candidate instruments including the percentage of income from Social Security Income, the size of the individual\"s firm and whether the firm has multiple locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.datasets import meps\n",
    "\n",
    "data = meps.load()\n",
    "data = data.dropna()\n",
    "print(meps.DESCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the data -- dependent, endogenous and controls -- are summarized.  The controls are grouped into a list to simplify model building."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "controls = [\"totchr\", \"female\", \"age\", \"linc\", \"blhisp\"]\n",
    "print(data[[\"ldrugexp\", \"hi_empunion\"] + controls].describe(percentiles=[]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also worth examining the instruments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instruments = [\"ssiratio\", \"lowincome\", \"multlc\", \"firmsz\"]\n",
    "print(data[instruments].describe(percentiles=[]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally the simple correlation between the endogenous variable and the instruments.  Instruments must be correlated to be relevant (but also must be exogenous, which can\"t be examined using simple correlation).  The correlation of `firmsz` is especially low, which might lead to the weak instruments problem if used exclusively.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[[\"hi_empunion\"] + instruments].corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``add_constant`` from ``statsmodels`` is used to simplify the process of adding a constant column to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.api import OLS, add_constant\n",
    "\n",
    "data[\"const\"] = 1\n",
    "controls = [\"const\"] + controls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2SLS as OLS\n",
    "Before examining the IV estimators, it is worth noting that 2SLS nests the OLS estimator, so that a call to ``IV2SLS`` using ``None`` for the endogenous and instruments will produce OLS estimates of parameters.\n",
    "\n",
    "The OLS estimates indicate that insurance through an employer or union leads to an **increase** in out-of-pocket drug expenditure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivolsmod = IV2SLS(data.ldrugexp, data[[\"hi_empunion\"] + controls], None, None)\n",
    "res_ols = ivolsmod.fit()\n",
    "print(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Just identified 2SLS\n",
    "\n",
    "The just identified two-stage LS estimator uses as many instruments as endogenous variables.  In this example there is one of each, using the SSI ratio as the instrument.  The with the instrument, the effect of insurance through employer or union has a strong negative effect on drug expenditure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IV2SLS(data.ldrugexp, data[controls], data.hi_empunion, data.ssiratio)\n",
    "res_2sls = ivmod.fit()\n",
    "print(res_2sls.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Instruments\n",
    "\n",
    "Using multiple instruments only requires expanding the data array in the instruments input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IV2SLS(\n",
    "    data.ldrugexp, data[controls], data.hi_empunion, data[[\"ssiratio\", \"multlc\"]]\n",
    ")\n",
    "res_2sls_robust = ivmod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative covariance estimators\n",
    "\n",
    "All estimator allow for three types of parameter covariance estimator:\n",
    "\n",
    "* ``\"unadjusted\"`` is the classic homoskedastic estimator\n",
    "* ``\"robust\"``  is robust to heteroskedasticity\n",
    "* ``\"clustered\"`` allows one- or two-way clustering to account for additional sources of dependence between the model scores\n",
    "* ``\"kernel\"`` produces a heteroskedasticity-autocorrelation robust covariance estimator\n",
    "\n",
    "The default is ``\"robust\"``.\n",
    "\n",
    "These are all passed using the keyword input ``cov_type``.  Using clustered requires also passing the clustering variable(s).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IV2SLS(\n",
    "    data.ldrugexp, data[controls], data.hi_empunion, data[[\"ssiratio\", \"multlc\"]]\n",
    ")\n",
    "res_2sls_std = ivmod.fit(cov_type=\"unadjusted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GMM Estimation\n",
    "\n",
    "GMM estimation can be more efficient than 2SLS when there are more than one instrument.  By default, 2-step efficient GMM is used (assuming the weighting matrix is correctly specified).  It is possible to iterate until convergence using the optional keyword input ``iter_limit``, which is naturally 2 by default.  Generally, GMM-CUE would be preferred to using multiple iterations of standard GMM.\n",
    "\n",
    "The default weighting matrix is robust to heteroskedasticity (but not clustering)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IVGMM(\n",
    "    data.ldrugexp, data[controls], data.hi_empunion, data[[\"ssiratio\", \"multlc\"]]\n",
    ")\n",
    "res_gmm = ivmod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing the weighting matrix structure in GMM estimation\n",
    "The weighting matrix in the GMM objective function can be altered when creating the model.  This example uses clustered weight by age.  The covariance estimator should usually match the weighting matrix, and so clustering is also used here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IVGMM(\n",
    "    data.ldrugexp,\n",
    "    data[controls],\n",
    "    data.hi_empunion,\n",
    "    data[[\"ssiratio\", \"multlc\"]],\n",
    "    weight_type=\"clustered\",\n",
    "    clusters=data.age,\n",
    ")\n",
    "res_gmm_clustered = ivmod.fit(cov_type=\"clustered\", clusters=data.age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuously updating GMM \n",
    "The continuously updating GMM estimator simultaneously optimizes the moment conditions and the weighting matrix.  It can be more efficient (in the second order sense) than standard 2-step GMM, although it can also be fragile.  Here the optional input ``display`` is used to produce the output of the non-linear optimizer used to estimate the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IVGMMCUE(\n",
    "    data.ldrugexp, data[controls], data.hi_empunion, data[[\"ssiratio\", \"multlc\"]]\n",
    ")\n",
    "res_gmm_cue = ivmod.fit(cov_type=\"robust\", display=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing results\n",
    "The function ``compare`` can be used to compare the results of multiple models, possibly with different variables, estimators and/or instruments. Usually a dictionary or ``OrderedDict`` is used to hold results since the keys are used as model names.  The advantage of an ``OrderedDict`` is that it will preserve the order of the models in the presentation.\n",
    "\n",
    "With the expectation of the OLS estimate, the parameter estimates are fairly consistent. Standard errors vary slightly although the conclusions reached are not sensitive to the choice of covariance estimator either. T-stats are reported in parentheses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "from linearmodels.iv.results import compare\n",
    "\n",
    "res = OrderedDict()\n",
    "res[\"OLS\"] = res_ols\n",
    "res[\"2SLS\"] = res_2sls\n",
    "res[\"2SLS-Homo\"] = res_2sls_std\n",
    "res[\"2SLS-Hetero\"] = res_2sls_robust\n",
    "res[\"GMM\"] = res_gmm\n",
    "res[\"GMM Cluster(Age)\"] = res_gmm_clustered\n",
    "res[\"GMM-CUE\"] = res_gmm_cue\n",
    "print(compare(res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing endogeneity\n",
    "\n",
    "The Durbin test is a classic of endogeneity which compares OLS estimates with 2SLS and exploits the fact that OLS estimates will be relatively efficient.  Durbin\"s test is not robust to heteroskedasticity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_2sls.durbin()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Wu-Hausman test is a variant of the Durbin test that uses a slightly different form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_2sls.wu_hausman()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test statistic can be directly replicated using the squared t-stat in a 2-stage approach where the first stage regresses the endogenous variable on the controls and instrument and the second stage regresses the dependent variable on the controls, the endogenous regressor and the residuals.  If the regressor was in fact exogenous, the residuals should not be correlated with the dependent variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "step1 = IV2SLS(data.hi_empunion, data[[\"ssiratio\"] + controls], None, None).fit()\n",
    "resids = step1.resids\n",
    "exog = pd.concat([data[[\"hi_empunion\"] + controls], resids], axis=1)\n",
    "step2 = IV2SLS(data.ldrugexp, exog, None, None).fit(cov_type=\"unadjusted\")\n",
    "print(step2.tstats.residual**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wooldridge\"s regression-based test of exogeneity is robust to heteroskedasticity since it inherits the covariance estimator from the model. Here there is little difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_2sls.wooldridge_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wooldridge\"s score test is an alternative to the regression test, although it usually has slightly less power since it is an LM rather than a Wald type test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_2sls.wooldridge_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exogeneity Testing\n",
    "When there is more than one instrument (the model is overidentified), the J test can be used in GMM models to test whether the model is overidentified -- in other words, whether the instruments are actually exogenous (assuming they are relevant).  In the case with 2 instruments there is no evidence that against the null."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_gmm.j_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When all instruments are included the story changes, and some of the additional instrument (`lowincome` or `firmsz`) appear to be endogenous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IVGMM(data.ldrugexp, data[controls], data.hi_empunion, data[instruments])\n",
    "res_gmm_all = ivmod.fit()\n",
    "res_gmm_all.j_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single Instrument Regressions\n",
    "It can be useful to run the just identified regressions to see how the IV estimate varies by instrument. The OLS model is included for comparison. The coefficient when using ``lowincome`` is very similar to the OLS as is the $R^2$ which indicates this variable may be endogenous. The coefficient using ``firmsz`` is also very different, but this is probably due to the low correlation between ``firmsz`` and the endogenous regressor so that this is a weak instrument. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "od = OrderedDict()\n",
    "for col in instruments:\n",
    "    od[col] = IV2SLS(data.ldrugexp, data[controls], data.hi_empunion, data[col]).fit(\n",
    "        cov_type=\"robust\"\n",
    "    )\n",
    "od[\"OLS\"] = res_ols\n",
    "print(compare(od))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Stage Diagnostics\n",
    "First stage diagnostics are available to assess whether the instruments appear to be credible for the endogenous regressor. The Partial F-statistic is the F-statistic for all instruments once controls have been partialed out.  In the case of a single instrument, it is just the squared t-stat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_2sls.first_stage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The F-statistic actually has a $chi^2$ distribution since it is just a Wald test that all of the coefficients are 0.  This breaks the \"rule-of-thumb\" but it can be applied by dividing the F-stat by the number of instruments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IV2SLS(data.ldrugexp, data[controls], data.hi_empunion, data[instruments])\n",
    "res_2sls_all = ivmod.fit()\n",
    "print(res_2sls_all.first_stage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LIML\n",
    "The LIML estimator and related k-class estimators can be used through ``IVLIML``.  LIML can have better finite sample properties if the model is not strongly identified. By default the $\\kappa$ parameter is estimated. In this dataset it is very close to 1 and to the results for LIML are similar to 2SLS (they would be exact if $\\kappa=1$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivmod = IVLIML(\n",
    "    data.ldrugexp, data[controls], data.hi_empunion, data[[\"ssiratio\", \"multlc\"]]\n",
    ")\n",
    "res_liml = ivmod.fit(cov_type=\"robust\")\n",
    "print(compare({\"2SLS\": res_2sls_robust, \"LIML\": res_liml, \"GMM\": res_gmm}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated value of $\\kappa$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_liml.kappa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IV2SLS as OLS\n",
    "\n",
    "As one final check, the \"OLS\" version of ``IV2SLS`` is compared to ``statsmodels`` ``OLS`` command. The parameters are identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "ivolsmod = IV2SLS(data.ldrugexp, data[[\"hi_empunion\"] + controls], None, None)\n",
    "res_ivols = ivolsmod.fit()\n",
    "sm_ols = res_ols.params\n",
    "sm_ols.name = \"sm\"\n",
    "print(pd.concat([res_ivols.params, sm_ols], axis=1))"
   ]
  }
 ],
 "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.0"
  },
  "nbsphinx": {
   "allow_errors": true
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}