econmt-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(regression)=\n",
    "# Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn how to run linear regressions with code.\n",
    "\n",
    "If you're running this code (either by copying and pasting it, or by downloading it using the icons at the top of the page), you may need to install the packages it uses first. There's a brief guide to installing packages in the chapter on {ref}`code-preliminaries`.\n",
    "\n",
    "Most of this chapter will rely on [**pyfixest**](https://github.com/s3alfisc/pyfixest) and we are indebted to [Alexander Fischer](https://github.com/s3alfisc) for his efforts in bringing a convenient and flexible way to absorb high dimensional fixed effects to Python via **pyfixest**. Some of the material in this chapter follows [Grant McDermott](https://grantmcdermott.com/)'s excellent notes and the [Library of Statistical Translation](https://lost-stats.github.io/).\n",
    "\n",
    "## Quick Review of Econometrics\n",
    "\n",
    "### Notation\n",
    "\n",
    "Greek letters, like $\\beta$, are the truth and represent parameters. Modified Greek letters are an estimate of the truth, for example $\\hat{\\beta}$. Sometimes Greek letters will stand in for vectors of parameters. Most of the time, upper case Latin characters such as $X$ will represent random variables (which could have more than one dimension). Lower case letters from the Latin alphabet denote realised data, for instance $x$ (which again could be multi-dimensional).  Modified Latin alphabet letters denote computations performed on data, for instance $\\bar{x} = \\frac{1}{n} \\displaystyle\\sum_{i} x_i$ where $n$ is number of samples.\n",
    "\n",
    "Ordinary least squares (OLS) regression can be used to *estimate* the parameters of certain types of model, most typically models of the form\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 \\cdot x_1 + \\beta_2 \\cdot x_2 .\n",
    "$$\n",
    "\n",
    "This generic model says that the value of an outcome variable $y$ is a linear function of one or more input predictor variables $x_i$, where the $x_i$ could be transforms of original data. But the above equation is a platonic ideal, what we call a data generating process (DGP). OLS allows us to recover *estimates* of the parameters of the model , i.e. to find $\\hat{\\beta_i}$ and to enable us to write an estimated model:\n",
    "\n",
    "$$\n",
    "y = \\hat{\\beta_0} + \\hat{\\beta_1} \\cdot x_1 + \\hat{\\beta_2} \\cdot x_2 + \\epsilon .\n",
    "$$\n",
    "\n",
    "This equation can also be expressed in matrix form as\n",
    "\n",
    "$$\n",
    "y = x'\\cdot \\hat{\\beta} + \\epsilon\n",
    "$$\n",
    "\n",
    "where $x' = (1, x_1, \\dots, x_{n})'$ and $\\hat{\\beta} = (\\hat{\\beta_0}, \\hat{\\beta_1}, \\dots, \\hat{\\beta_{n}})$.\n",
    "\n",
    "Given data $y_i$ stacked to make a vector $y$ and $x_{i}$ stacked to make a matrix $X$, this can be solved for the coefficients $\\hat{\\beta}$ according to\n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = \\left(X'X\\right)^{-1} X'y .\n",
    "$$\n",
    "\n",
    "### Gauss-Markov Conditions\n",
    "\n",
    "To be sure that the estimates of these parameters are the *best linear unbiased estimate*, a few conditions need to hold: the Gauss-Markov conditions:\n",
    "\n",
    "1. $y$ is a linear function of the $\\beta_i$\n",
    "2. $y$ and the $x_i$ are randomly sampled from the population.\n",
    "3. There is no perfect multi-collinearity of variables.\n",
    "4. $\\mathbb{E}(\\epsilon | x_1, \\dots, x_n) = 0$ (unconfoundedness)\n",
    "5. $\\text{Var}(\\epsilon | x_1, \\dots, x_n) = \\sigma^2$ (homoskedasticity)\n",
    "\n",
    "(1)-(4) also guarantee that OLS estimates are unbiased and $\\mathbb{E}(\\hat{\\beta}_i) = \\beta_i$.\n",
    "\n",
    "The classic linear model requires a 6th assumption; that $\\epsilon \\thicksim \\mathcal{N}(0, \\sigma^2)$.\n",
    "\n",
    "The interpretation of regression coefficients depends on what their units are to begin with, but you can always work it out by differentiating both sides of the model equation with respect to the $x_i$. For example, for the first model equation above\n",
    "\n",
    "$$\n",
    "\\frac{\\partial y}{\\partial x_i} = \\beta_i\n",
    "$$\n",
    "\n",
    "so we get the interpretation that $\\beta_i$ is the rate of change of y with respect to $x_i$. If $x_i$ and $y$ are in levels, this means that a unit increase in $x_i$ is associated with a $\\beta_i$ units increase in $y$. If the right-hand side of the model is $\\ln x_i$ then we get\n",
    "\n",
    "$$\n",
    "\\frac{\\partial y}{\\partial x_i} = \\beta_i \\frac{1}{x_i} \n",
    "$$\n",
    "\n",
    "with some abuse of notation, we can rewrite this as $\\partial y = \\beta_i \\partial x_i/x_i$, which says that a percent change in $x_i$ is associated with a $\\beta_i$ unit change in $y$. With a logged $y$ variable, it's a percent change in $x_i$ that is associated with a percent change in $y$, or $\\partial y/y = \\beta_i \\partial x_i/x_i$ (note that both sides of this equation are unitless in this case). Finally, another example that is important in practice is that of log differences, eg $y = \\beta_i (\\ln x_i - \\ln x_i')$. Again, we will abuse notation and say that this case may be represented as $\\partial y = \\beta_i (\\partial x_i/x_i - \\partial x_i'/x_i')$, i.e. the difference in two percentages, a *percentage point* change, in $x_i$ is associated with a $\\beta_i$ unit change in $y$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports\n",
    "\n",
    "Let's import some of the packages we'll be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from lets_plot import *\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import pyfixest as pf\n",
    "\n",
    "LetsPlot.setup_html()\n",
    "\n",
    "# Set seed for random numbers\n",
    "seed_for_prng = 78557\n",
    "prng = np.random.default_rng(seed_for_prng)  # prng=probabilistic random number generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression basics\n",
    "\n",
    "There are two ways to run regressions: passing the data directly as objects, and using formulae. For most situations, the formula API is more convenient and so most of what we'll cover uses this.\n",
    "\n",
    "For our workhorse regression package, we'll be using **pyfixest**, which has [impressive speed](https://github.com/s3alfisc/pyfixest/issues/202#issuecomment-1819958417).\n",
    "\n",
    "We'll use the starwars dataset to run a regression of mass on height for star wars characters. This example borrows very heavily from notes by [Grant McDermott](https://grantmcdermott.com/). First, let's bring the dataset in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/data/starwars.csv\",\n",
    "    index_col=0,\n",
    ")\n",
    "# Look at first few rows\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, now let's do a regression using OLS and a formula that says our y-variable is mass and our regressor is height:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_sw = pf.feols(\"mass ~ height\", data=df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, where are the results!? They're stored in the object we created. To peek at them we need to call the summary function (and, for easy reading, I'll print it out too using `print`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_sw.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we're seeing here are really several bits of information glued together. To just grab the coefficient results in a tidy format, use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_sw.tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll have noticed that we got an intercept, even though we didn't specify one in the formula. **pyfixest** adds in an intercept by default because, most of the time, you will want one. To turn it off, add a `-1` at the end of the formula command, eg in this case you would call `feols(\"mass ~ height -1\", data=df)`.\n",
    "\n",
    "The fit we got in the case with the intercept was pretty terrible; a low $R^2$ and both of our confidence intervals are large and contain zero. What's going on? If there's one adage in regression that's always worth paying attention to, it's *always plot your data*. Let's see what's going on here (using **matplotlib** for plotting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(df[\"height\"], df[\"mass\"], s=200, alpha=0.8)\n",
    "ax.annotate(\n",
    "    \"Jabba the Hutt\",\n",
    "    df.iloc[df[\"mass\"].idxmax()][[\"height\", \"mass\"]],\n",
    "    xytext=(0, -50),\n",
    "    textcoords=\"offset points\",\n",
    "    arrowprops=dict(\n",
    "        arrowstyle=\"fancy\",\n",
    "        color=\"k\",\n",
    "        connectionstyle=\"arc3,rad=0.3\",\n",
    "    ),\n",
    ")\n",
    "ax.set_ylim(0, None)\n",
    "ax.set_title(\"Always plot the data\", loc=\"left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh dear, Jabba's been on the paddy frogs again, and he's a bit of different case. When we're estimating statistical relationships, we have all kinds of choices and should be wary about arbitrary decisions of what to include or exclude in case we fool ourselves about the generality of the relationship we are capturing. Let's say we knew that we weren't interested in Hutts though, but only in other species: in that case, it's fair enough to filter out Jabba and run the regression without this obvious outlier. We'll exclude any entry that contains the string 'Jabba' in the `name` column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_outlier_free = pf.feols(\n",
    "    \"mass ~ height\", data=df[~df[\"name\"].str.contains(\"Jabba\")]\n",
    ")\n",
    "print(results_outlier_free.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks a lot more healthy. Not only is the model explaining a *lot* more of the data, but the coefficients are now significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard errors\n",
    "\n",
    "You'll have seen that there's a column for the standard error of the estimates in the regression table and a message saying that the covariance type of these is 'nonrobust'. Let's say that, instead, we want to use Eicker-Huber-White robust standard errors, aka \"HC2\" standard errors. We can specify to use these up front standard errors up front in the fit method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"mass ~ height\", data=df, vcov=\"HC2\").summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, alternatively, we can go back to our existing results and get new standard errors \"on the fly\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_sw.vcov(\"HC2\").summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several different types of standard errors available in **pyfixest**:\n",
    "\n",
    "- \"iid\", \"HC1\", \"HC2\", \"HC3\"\n",
    "- \"hetero\", for heteroskedasticity and autocorrelation consistent standard errors, for which you may want to also use some keyword arguments\n",
    "- clustered standard errors: see below. CRV1 and CRV3 available.\n",
    "\n",
    "You can find information on all of these [here](https://s3alfisc.github.io/pyfixest/tutorial/#standard-errors-and-inference).\n",
    "\n",
    "For now, let's look more closely at those last ones: clustered standard errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clustered standard errors\n",
    "\n",
    "Often, we know something about the structure of likely errors, namely that they occur in groups. In the below example we use one-way clusters to capture this effect in the errors.\n",
    "\n",
    "Note that in the below example, we grab a subset of the data for which a set of variables we're interested in are defined, otherwise the below example would execute with an error because of missing cluster-group values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf = df.dropna(subset=[\"homeworld\", \"mass\", \"height\", \"species\"])\n",
    "pf.feols(\"mass ~ height\", data=xf, vcov={\"CRV1\": \"homeworld\"}).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add two-way clustering of standard errors using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"mass ~ height\", data=xf, vcov={\"CRV1\": \"homeworld + species\"}).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you would generally expect, the addition of clustering has increased the standard errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed effects and categorical variables\n",
    "\n",
    "Fixed effects are a way of allowing the intercept of a regression model to vary freely across individuals or groups. It is, for example, used to control for any individual-specific attributes that do not vary across time in panel data.\n",
    "\n",
    "Let's use the 'mtcars' dataset to demonstrate this. We'll read it in and set the datatypes of some of the columns at the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg = pd.read_csv(\n",
    "    \"https://raw.githubusercontent.com/LOST-STATS/lost-stats.github.io/source/Data/mtcars.csv\",\n",
    "    dtype={\"model\": str, \"mpg\": float, \"hp\": float, \"disp\": float, \"cyl\": \"category\"},\n",
    ")\n",
    "mpg.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have our data in we want to regress mpg (miles per gallon) on hp (horsepower) with fixed effects for cyl (cylinders). Now we *could* just pop in a formula like this `'mpg ~ hp + cyl'` because we took the trouble to declare that `cyl` was of datatype category when reading it in from the csv file. This means that statsmodels will treat it as a category and use it as a fixed effect by default.\n",
    "\n",
    "But when I read that formula I get nervous that `cyl` might not have been processed correctly (ie it could have been read in as a float, which is what it looks like) and it might just be treated as a float (aka a continuous variable) in the regression. Which is not what we want at all. So, to be safe, and make our intentions explicit (even when the data is of type 'category'), it's best to use the syntax `C(cyl)` to ask for a fixed effect.\n",
    "\n",
    "Here's a regression which does that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"mpg ~ hp + C(cyl)\", data=mpg).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see here that two of the three possible values of `cyl`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"cyl\"].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "have been added as fixed effects regressors. The way that `+C(cyl)` has been added makes it so that the coefficients given are relative to the coefficient for the intercept. We can turn the intercept off to get a coefficient per unique `cyl` value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"mpg ~ hp + C(cyl) -1\", data=mpg).tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When there is an intercept, the coefficients of fixed effect variables can be interpreted as being the average of $y$ for that class *compared* to the excluded classes holding all other categories and variables fixed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### High dimensional fixed effects, aka absorbing regression\n",
    "\n",
    "Sometimes, you just have a LOT of fixed effects (and perhaps you don't particularly care about them individually). A common example is having a large number of firms as part of a panel. Fortunately, there are ways to make regression with high dimensional fixed effects be both fast and concise. (In Stata, this is provided by the `reghdfe` package.)\n",
    "\n",
    "Let's say we have a regression of the form\n",
    "\n",
    "$$\n",
    "y_i = x_i\\cdot \\beta + z_i\\cdot \\gamma +\\epsilon_i\n",
    "$$\n",
    "\n",
    "where $y_i$ are observations indexed by $i$, $x_i$ are vectors of exogenous variables we care about the coefficients ($\\beta$), $z_i$ are vectors of fixed effects we don't care too much about the coefficients ($\\gamma$) for, and the $\\epsilon_i$ are errors. Then we can use an *absorbing regression* to solve for the $\\beta$ while ignoring the $\\gamma$.\n",
    "\n",
    "Here's an example using simulated data on workers taken from the docs for a different regression package ([**linearmodels**](https://bashtage.github.io/linearmodels//)). Let's simulate some data first, with two fixed effects (state and firm) alongside the two exogenous variables we're interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import default_rng\n",
    "\n",
    "rng = default_rng()  # Random number generator\n",
    "\n",
    "# Create synthetic input data\n",
    "nobs = 1_000_000  # No. observations\n",
    "state_id = rng.integers(50, size=nobs)  # State identifier\n",
    "firm_id = rng.integers(nobs // 5, size=nobs)  # Firm identifier (mean of 5 workers/firm)\n",
    "x = rng.standard_normal((nobs, 2))  # Exogenous variables\n",
    "sim = pd.DataFrame(\n",
    "    {\n",
    "        \"state_id\": pd.Categorical(state_id),\n",
    "        \"firm_id\": pd.Categorical(firm_id),\n",
    "        \"exog_0\": x[:, 0],\n",
    "        \"exog_1\": x[:, 1],\n",
    "    }\n",
    ")\n",
    "\n",
    "# Create synthetic relationship\n",
    "beta = [1, 3]  # coefficients of interest\n",
    "state_effects = rng.standard_normal(state_id.max() + 1)\n",
    "state_effects = state_effects[state_id]  # Generate state fixed effects\n",
    "firm_effects = rng.standard_normal(firm_id.max() + 1)\n",
    "firm_effects = firm_effects[firm_id]  # Generate firm fixed effects\n",
    "eps = rng.standard_normal(nobs)  # Generate errors\n",
    "# Generate endogeneous outcome variable\n",
    "sim[\"y\"] = (\n",
    "    sim[\"exog_0\"] * beta[0]\n",
    "    + sim[\"exog_1\"] * beta[1]\n",
    "    + firm_effects\n",
    "    + state_effects\n",
    "    + eps\n",
    ")\n",
    "sim.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we pass this to **pyfixest** and with the `state_id` and `firm_id` variables entered after a vertical bar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_hdfe = pf.feols(\"y ~ exog_0 + exog_1 | state_id + firm_id\", data=sim)\n",
    "results_hdfe.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, from our 1,000,000 observations, we have roughly 200,000 fixed effects that have been scooped up and packed away, leaving us with just the coefficients, $\\beta$, on the exogenous variables of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations of regressors\n",
    "\n",
    "This chapter is showcasing *linear* regression. What that means is that the model is linear in the regressors: but it doesn't mean that those regressors can't be some kind of (potentially non-linear) transform of the original features $x_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logs and arcsinh\n",
    "\n",
    "You have two options for adding in logs: do them before, or do them in the formula. Doing them before just makes use of standard dataframe operations to declare a new column:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"lnhp\"] = np.log(mpg[\"hp\"])\n",
    "pf.feols(\"mpg ~ lnhp\", data=mpg).tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can specify the log directly in the formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_ln = pf.feols(\"mpg ~ np.log(hp)\", data=mpg)\n",
    "results_ln.tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the first method will work for `arcsinh(x)` and `log(x+1)`, but you can also pass both of these into the formula directly too. (For more on the pros and cons of arcsinh, see {cite:t}`bellemare2020elasticities`.) Here it is with arcsinh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"mpg ~ np.arcsinh(hp)\", data=mpg).tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interaction terms and powers\n",
    "\n",
    "This chapter is showcasing *linear* regression. What that means is that the model is linear in the regressors: but it doesn't mean that those regressors can't be some kind of non-linear transform of the original features $x_i$. Two of the most common transformations that you might want to use are *interaction terms* and *polynomial terms*. An example of an interaction term would be\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 x_1 \\cdot x_2\n",
    "$$\n",
    "\n",
    "while an example of a polynomial term would be\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 x_1^2\n",
    "$$\n",
    "\n",
    "i.e. the last term enters only after it is multiplied by itself.\n",
    "\n",
    "One note of warning: the interpretation of the effect of a variable is no longer as simple as was set out at the start of this chapter. To work out *what* the new interpretation is, the procedure is the same though: just take the derivative. In the case of the interaction model above, the effect of a unit change in $x_1$ on $y$ is now going to be a function of $x_2$. In the case of the polynomial model above, the effect of a unit change in $x_1$ on $y$ will be $2\\beta_1 \\cdot x_1$. For more on interaction terms, see {cite:t}`balli2013interaction`.\n",
    "\n",
    "Alright, with all of that preamble out of the way, let's see how we actual do some of this! Let's try including a linear and squared term in the regression of `mpg` on `hp` making use of the numpy power function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_poly = pf.feols(\"mpg ~ hp + np.power(hp, 2)\", data=mpg)\n",
    "res_poly.tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's include the original term in hp, a term in disp, and the interaction between them, which is represented by hp:disp in the table. We can do this by multiplying out the terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_inter = pf.feols(\"mpg ~ hp * disp\", data=mpg)\n",
    "res_inter.tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the unusual case that you want *only* the interaction term, you write it as it appears in the table above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_only_inter = pf.feols(\"mpg ~ hp : disp\", data=mpg)\n",
    "res_only_inter.tidy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The formula API explained\n",
    "\n",
    "As you will have seen `~` separates the left- and right-hand sides of the regression. `+` computes a set union, which will also be familiar from the examples above (ie it includes two terms as long as they are distinct). `-` computes a set difference; it adds the set of terms to the left of it while removing any that appear on the right of it. As we've seen, `a*b` is a short-hand for `a + b + a:b`, with the last term representing the interaction. `/` is short hand for `a + a:b`, which is useful if, for example `b` is nested within `a`, so it doesn't make sense to control for `b` on its own. Actually, the `:` character can interact multiple terms so that `(a + b):(d + c)` is the same as `a:c + a:d + b:c + b:d`. `C(a)` tells statsmodels to treat `a` as a categorical variable that will be included as a fixed effect. Finally, as we saw above with powers, you can also pass in vectorised functions, such as `np.log()` and `np.power()`, directly into the formulae.\n",
    "\n",
    "One gotcha with the formula API is ensuring that you have sensible variable names in your dataframe, i.e. ones that do *not* include whitespace or, to take a really pathological example, have the name 'a + b' for one of the columns that you want to regress on. You can dodge this kind of problem by passing in the variable name as, for example, `Q(\"a + b\")` to be clear that the *column name* is anything within the `Q(\"...\")`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple regression models\n",
    "\n",
    "As is so often the case, you're likely to want to run more than one model at once with different specifications. There are two  ways you may wish to do this.\n",
    "\n",
    "### Comparing a set of models you've already created\n",
    "\n",
    "In the code above, we created three different models: `results_ln`, `res_poly`, and `res_inter`. Let's say we now want to compare them in a table: we can! By using **pyfixest**'s `etable()` function. (If you need to do this in **statsmodels**, use `summary_col()` as explained [here](https://aeturrell.com/blog/posts/running-many-regressions-alongside-pandas/) or [**stargazer**](https://github.com/mwburke/stargazer).)\n",
    "\n",
    "Let's see `etable` in action:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.etable([results_ln, res_poly, res_inter], type=\"df\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few different `type=` options, including dataframe, markdown, and latex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.etable([results_ln, res_poly, res_inter], type=\"md\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And of course the latex one (`type=\"tex\"`) can be written to a file using `open('regression.tex', 'w').write(...)` where the ellipsis is the relevant `etable()` command."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stepwise multiple models\n",
    "\n",
    "Sometimes you want to use a formula and add variables in step-wise. There are 4 stepwise functions: `sw`, `sw0`, `csw`, `csw0`. \n",
    "\n",
    "Assume you have the following formula: `y ~ x1 + sw(x2, x3)`. The stepwise function sw will estimate the following two models: `y ~ x1 + x2` and `y ~ x1 + x3`. That is, each element in `sw()` is sequentially, and separately, added to the formula. If you use `sw0` in lieu of `sw`, then the model `y ~ x1` would also be estimated: the `0` in the name implies that the model without any stepwise element will also be estimated. This explains how to use `sw` and `sw0`.\n",
    "\n",
    "The prefix c in the remaining two options, `csw` and `csw0`, means \"cumulative\": each stepwise element is added to the next. That is, `y ~ x1 + csw(x2, x3)` would cause `y ~ x1 + x2` and `y ~ x1 + x2 + x3` to be estimated. The `0` has the same meaning, ie it estimates the most basic version. As an example, `y ~ x1 + csw0(x2, x3)` leads to the following three models being estimated: `y ~ x1`, `y ~ x1 + x2` and `y ~ x1 + x2 + x3`.\n",
    "\n",
    "Let's see some examples (courtesy of the brilliant **fixest** documentation for the statistical language R). We'll be using the iris dataset, which comes with the **plotly** package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.express as px\n",
    "iris = px.data.iris()\n",
    "iris.columns = [\"y\", \"x1\", \"x2\", \"x3\", \"species\", \"id\"]\n",
    "iris.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try a regular (not cumulative) stepwise regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"y ~ sw(x1, x2, x3)\", data=iris).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a cumulative stepwise regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.feols(\"y ~ csw(x1, x2, x3)\", data=iris).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a very verbose way to get regression coefficients out! We can get a summary of the info by using `etable` again, but this time as a *method* rather than a stand-alone function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_iris_csw = pf.feols(\"y ~ csw(x1, x2, x3)\", data=iris)\n",
    "reg_iris_csw.etable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even handier, **pyfixest** lets you visualise these results (via **Lets-Plot**):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_iris_csw.coefplot().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instrumental variables\n",
    "\n",
    "We will again be using **pyfixest** for instrumental variables (IV) regression. This sub-section is indebted to another package, [**linearmodels**](https://bashtage.github.io/linearmodels//doc/index.html), for the following example.\n",
    "\n",
    "\n",
    "Recall that a good instrumental variable $z$ has zero covariance with the error from the regression (which is untestable) and non-zero covariance with the variable of interest (which is).\n",
    "\n",
    "\n",
    "Recall that in IV regression, we have a model of the form\n",
    "\n",
    "$$\n",
    "\\begin{split}y_i    & = x_{1i}\\hat{\\beta}_1 + x_{2i}\\hat{\\beta}_2 + \\epsilon_i \\\\\n",
    "x_{2i} & = z_{1i}\\hat{\\delta} + z_{2i}\\hat{\\gamma} + \\nu_i\\end{split}\n",
    "$$\n",
    "\n",
    "where $x_{1i}$ is a set of $k_1$ exogenous regressors and $x_{2i}$ is a set of $k_2$ endogenous regressors such that $\\text{Cov}(x_{2i}, \\epsilon_i)\\neq 0$. This is a problem for the usual OLS assumptions (the right-hand side should be exogenous). \n",
    "\n",
    "\n",
    "To get around this, in 2-stage least squares IV, we first regress $x_{2i}$ on instruments that explain $x_{2i}$ *but not* $y_i$, and then regress $y_i$ only on the predicted/estimated left-hand side from the first regression, ie on $\\hat{x}_{2i}$.\n",
    "\n",
    "\n",
    "\n",
    "It's always easiest to see an example, so let's estimate what might cause (realised) cigarette demand for the 48 continental US states in 1995."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfiv = pd.read_csv(\n",
    "    \"https://vincentarelbundock.github.io/Rdatasets/csv/AER/CigarettesSW.csv\",\n",
    "    dtype={\"state\": \"category\", \"year\": \"category\"},\n",
    ").assign(\n",
    "    rprice=lambda x: x[\"price\"] / x[\"cpi\"],\n",
    "    rincome=lambda x: x[\"income\"] / x[\"population\"] / x[\"cpi\"],\n",
    ")\n",
    "dfiv.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll specify the model. It's going to be in the form `dep ~ exog + [endog ~ instruments]`, where endog will be regressed on instruments and dep will be regressed on both exog and the predicted values of endog.\n",
    "\n",
    "In this case, the model will be\n",
    "\n",
    "$$\n",
    "\\text{Price}_i = \\hat{\\pi}_0 + \\hat{\\pi}_1 \\text{SalesTax}_i + v_i \n",
    "$$\n",
    "\n",
    "in the first stage regression and\n",
    "\n",
    "$$\n",
    "\\text{Packs}_i = \\hat{\\beta}_0 + \\hat{\\beta}_2\\widehat{\\text{Price}_i} + \\hat{\\beta}_1 \\text{RealIncome}_i + u_i\n",
    "$$\n",
    "\n",
    "in the second stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_iv = pf.feols(\"np.log(packs) ~ np.log(rincome) + 1 | year + state | np.log(rprice) ~ taxs \", data=dfiv, vcov={\"CRV1\": \"year\"})\n",
    "results_iv.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the IV results against (naive) OLS. First, run the OLS equivalent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_cigs_ols = pf.feols(\"np.log(packs) ~ np.log(rprice) + np.log(rincome) | year + state\", data=dfiv, vcov={\"CRV1\": \"year\"})\n",
    "results_cigs_ols.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do a proper comparison of these models using `etable()`. Note the order: the IV results are the second column of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.etable([results_cigs_ols, results_iv], type=\"md\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we take into account the fact that the real price is endogenous to (realised) demand, we find that its coefficient is more negative; i.e. an increase in the real price of cigarettes creates a bigger fall in number of packs bought."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3.8.13 ('codeforecon')",
   "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.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "caf5ac9f613b176c5984ad2a1a4525760eb7d898a3291351da4c152dc719ffa1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}