econmt-generalised-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(generalised-models)=\n",
    "# Generalised regression models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn how to run more generalised regressions with code. This will include logit, probit, generalised linear models, random effects, quantile regressions, and rolling + recursive regressions.\n",
    "\n",
    "Most of this chapter will rely on the amazing [**statsmodels**](https://www.statsmodels.org/stable/index.html) package, including for some of the examples. 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/)."
   ]
  },
  {
   "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",
    "\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\", 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logit, probit, and generalised linear models\n",
    "\n",
    "For this section, we'll be using the [**statsmodels**](https://www.statsmodels.org/) package. **statsmodels** has been around a long time (and is more mature than **pyfixest**), but it doesn't absorb fixed effects, and it doesn't have an out-of-the-box instrumental variables regression function. The two packages are similar in that you can use formulae in both.\n",
    "\n",
    "### Logit\n",
    "\n",
    "A logistical regression, aka a logit, is a statistical method for a best-fit line between a regressors $X$ and an outcome varibale $y$ that takes on values in $(0, 1)$.\n",
    "\n",
    "The function that we're assuming links the regressors and the outcome has a few different names but the most common is the sigmoid function or the logistic function. The data generating process is assumed to be\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mathbb{P}(Y=1\\mid X) = \\frac{1}{1 + e^{-X'\\beta}}}\n",
    "$$\n",
    "\n",
    "we can also write this as $\\ln\\left(\\frac{p}{p-1}\\right) = \\beta_0 + \\sum_i \\beta_i x_i$ to get a 'log-odds' relationship. The coefficients from a logit model do not have the same interpration as in an OLS estimation, and you can see this from the fact that $\\partial y/\\partial x_i \\neq \\beta_i$ for logit. Of course, you can work out what the partial derivative is for yourself but most packages offer a convenient way to quickly recover the marginal effects. \n",
    "\n",
    "Logit models are available in **scikit-learn** and **statsmodels** but bear in mind that the **scikit-learn** logit model is, ermm, extremely courageous in that regularisation is applied by default. If you don't know what that means, don't worry, but it's probably best to stick with **statsmodels** as we will do in this example.\n",
    "\n",
    "We will predict a target `GRADE`, representing whether a grade improved or not, based on some regressors including participation in a programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data from Spector and Mazzeo (1980)\n",
    "df = sm.datasets.spector.load_pandas().data\n",
    "# Look at info on data\n",
    "print(sm.datasets.spector.NOTE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_logit = smf.logit(\"GRADE ~ GPA + TUCE + PSI\", data=df).fit()\n",
    "print(res_logit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, did participation (`PSI`) help increase a grade? Yes. But we need to check the marginal effect to say exactly how much. We'll use `get_margeff()` to do this, we'd like the derivative of y with respect to x effect, and we'll take it at the mean of each regressor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marg_effect = res_logit.get_margeff(at=\"mean\", method=\"dydx\")\n",
    "marg_effect.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So participation gives almost half a grade increase."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probit\n",
    "\n",
    "Probit is very similar to logit: it's a statistical method for a best-fit line between regressors $X$ and an outcome varibale $y$ that takes on values in $(0, 1)$. And, just like with logit, the function that we're assuming links the regressors and the outcome has a few different names! \n",
    "\n",
    "The data generating process is assumed to be\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mathbb{P}(Y=1\\mid X)=\\Phi (X^{T}\\beta )}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\Phi (x)={\\frac {1}{\\sqrt {2\\pi }}}\\int _{-\\infty }^{x}e^{-{\\frac {y^{2}}{2}}}dy.}\n",
    "$$\n",
    "\n",
    "is the cumulative standard normal (aka Gaussian) distribution. The coefficients from a probit model do not have the same interpration as in an OLS estimation, and you can see this from the fact that $\\partial y/\\partial x_i \\neq \\beta_i$ for probit. And, just as with logit, although you can derive the marginal effects, most packages offer a convenient way to quickly recover them. Of particular note is the [**pymarginaleffects**](https://github.com/vincentarelbundock/pymarginaleffects) package by Vincent Arel-Bundock.\n",
    "\n",
    "We can re-use our previous example of predicting a target `GRADE`, representing whether a grade improved or not, based on some regressors including participation (PSI) in a programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_probit = smf.probit(\"GRADE ~ GPA + TUCE + PSI\", data=df).fit()\n",
    "print(res_probit.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_marg_effect = res_probit.get_margeff(at=\"mean\", method=\"dydx\")\n",
    "p_marg_effect.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's no coincidence that we find very similar results here because the two functions we're using don't actually look all that different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "support = np.linspace(-6, 6, 1000)\n",
    "ax.plot(support, st.logistic.cdf(support), ls=\"--\", label=\"Logistic\")\n",
    "ax.plot(support, st.norm.cdf(support), label=\"Probit\")\n",
    "ax.legend()\n",
    "ax.set_ylim(0, None)\n",
    "ax.set_ylim(0, None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What difference there is, is that logistic regression puts more weight into the tails of the distribution. Arguably, logit is easier to interpret too. With logistic regression, a one unit change in $x_i$ is associated with a $\\beta_i$ change in the log odds of a 1 outcome or, alternatively, an $e^{\\beta_i}$-fold change in the odds, all else being equal. With a probit, this is a change of $\\beta_i z$ for $z$ a normalised variable that you'd have to convert into a predicted probability using the normal CDF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generalised linear models\n",
    "\n",
    "Logit and probit (and OLS for that matter) as special cases of a class of models such that $g$ is a 'link' function connects a function of regressors to the output, and $\\mu$ is the mean of a conditional response distribution at a given point in the space of regressors. When $g(\\mu) = X'\\beta$, we just get regular OLS. When it's logit, we have\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mu= \\mathbb{E}(Y\\mid X=x) =g^{-1}(X'\\beta)= \\frac{1}{1 + e^{-X'\\beta}}.}\n",
    "$$\n",
    "\n",
    "But as well as the ones we've seen, there are many possible link functions one can use via the catch-all `glm` function. These come in different 'families' of distributions, with the default for the binomial family being logit. So, running `smf.glm('GRADE ~ GPA + TUCE + PSI', data=df, family=sm.families.Binomial()).fit()` will produce exactly the same as we got both using the `logit` function. For more on the families of distributions and possible link functions, see the [relevant part](https://www.statsmodels.org/stable/glm.html#) of the **statsmodels** documentation.\n",
    "\n",
    "If you need a library dedicated to GLMs that has all the bells and whistles you can dream of, you might want to check out [glum](https://glum.readthedocs.io/en/latest/). At the time of writing, it is [faster](https://glum.readthedocs.io/en/latest/benchmarks.html) than either GLMnet or H2O (two other popular GLM libraries)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Mixed Effects Regressions (LMER)\n",
    "\n",
    "Linear mixed models are an extension of simple linear models to allow both fixed and random effects. It's worth us having a quick recap of random effects before we dig into it.\n",
    "\n",
    "### Fixed vs Random Effects\n",
    "\n",
    "As this is going to come up in this chapter, it's worth having a quick review of the difference between fixed and random effects in regression-based models. There's a bit of disagreement about what people really mean when they say \"fixed effects\" or \"random effects\", so let's define them clearly here.\n",
    "\n",
    "> Random effects are estimated with partial pooling, while fixed effects are not.\n",
    "\n",
    "What does this mean? Well, technically, if an effect is assumed to be a realised value of a random variable, it is called a random effect. We could say that, in frequentist land (where parameters have true values), a single fixed effect model might look like:\n",
    "\n",
    "$$\n",
    "y = \\gamma\\cdot d\n",
    "$$\n",
    "\n",
    "where $d \\in \\{0, 1\\}$, which could be whether an individual belongs to, say, group 0 or group 1 and $\\gamma$ is a coefficient that is to be estimated. $\\gamma$ has a single \"true\" value, and we can admire the value of $\\hat{\\gamma}$ in a regression table once we have done our estimation.\n",
    "\n",
    "In contrast, even in a frequentist model, random effects are allowed to vary: they are drawn from a distribution (or, more accurately, their deviations from the true value are drawn from a distribution). This might be written as\n",
    "\n",
    "$$\n",
    "y_{i}= w_{i}\n",
    "$$\n",
    "\n",
    "and it could be that, instead of a group-specific intercept as in the fixed effects model above, the intercept can take a *range of values* depending on the individual. Although it's not always true, another way of thinking about this is that fixed effects are constant across individuals while random effects can vary.\n",
    "\n",
    "So what's this business about partial pooling? Imagine we have a situation where one group is under-sampled relative to the other groups. With few data points, we may not be able to precisely estimate fixed effects. However, partial pooling means that, if you do have few data points within a group, the group effect estimate can be partially based on the more abundant data from other groups. Random effects can provide a good compromise between estimating an effect by completely pooling all groups, which masks group-level variation, and estimating an effect for all groups completely separately, which could give imprecise estimates for low-sample groups. For example, when estimating means using a random effects approach, subgroup means are allowed to deviate a bit from the mean of the larger group, but not by an arbitrary amount—those deviations follow a distribution, typically a Gaussian (aka Normal) distribution. That's where the \"random\" in random effects comes from: the deviations of subgroups from a \"parent group\" follow the distribution of a random variable.\n",
    "\n",
    "The rule of thumb when thinking about how many different categories should exist for a variable to get the random, rather than fixed, effects treatment is that there should be more than 5 or 6 levels at a minimum in order to apply random effects—fewer than this and it's hard to estimate the variance around the central estimate. Because you could potentially precisely estimated a fixed effect if you had lots of samples in those 5 or 6 levels, you're more likely to reach for random effects when you have fewer samples per class.\n",
    "\n",
    "Random effects are the extension of this partial pooling technique to a wide range of situations: including multiple predictors, mixed continuous and categorical variables, and more.\n",
    "\n",
    "As a concrete example, taken from a [forum post](https://stats.stackexchange.com/questions/4700/what-is-the-difference-between-fixed-effect-random-effect-and-mixed-effect-mode) by Paul of MassMutual, suppose you want to estimate average US household income by ZIP code. You have a large dataset containing observations of household incomes and ZIP codes. Some ZIP codes are well represented in the dataset, but others have only a couple of households.\n",
    "\n",
    "For your initial model you would most likely take the mean income in each ZIP. This will work well when you have lots of data for all ZIPs, but the estimates for your poorly sampled ZIPs will suffer from high variance. You can mitigate this by using a shrinkage estimator (aka partial pooling), which will push extreme values towards the mean income across all ZIP codes.\n",
    "\n",
    "But how much shrinkage/pooling should you do for a particular ZIP? Intuitively, it should depend on the following:\n",
    "\n",
    "1. How many observations you have in that ZIP\n",
    "2. How many observations you have overall\n",
    "3. The individual-level mean and variance of household income across all ZIP codes\n",
    "4. The group-level variance in mean household income across all ZIP codes\n",
    "\n",
    "Now, if you model ZIP code as a random effect, the mean income estimate in all ZIP codes will be subjected to a statistically well-founded shrinkage, taking into account all of these factors.\n",
    "\n",
    "Random effects models automatically handle 4, the variability estimation, for all random effects in the model. This is useful, as it would be hard to do manually. Having accounted for (1)–(4), a random effects model is able to determine the appropriate shrinkage for low-sample groups. It can also handle much more complicated models with many different predictors.\n",
    "\n",
    "**Mixed effects models** are models that combine random and fixed effects.\n",
    "\n",
    "### Example with Random Effects\n",
    "\n",
    "We'll be using the **statsmodels** package to demonstrate this, along with the dietox dataset, which has data on the growth curves of pigs. The data are taken from a study of the \"Influence of Dietary Rapeseed Oli, Vitamin E, and Copper on Performance and Antioxidant and Oxidative Status of Pigs\" {cite:ps}`lauridsen1999influence`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.get_rdataset('dietox', 'geepack').data\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to regress the weight of the pigs on the time and a couple of random effects.\n",
    "\n",
    "We fit two random effects for each pig: a random intercept, and a random slope (with respect to time). This means that each pig may have a different baseline weight, as well as grow at a different rate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = smf.mixedlm(\"Weight ~ Time\", data, groups=data[\"Pig\"], re_formula=\"~Time\")\n",
    "mdf = md.fit(method=[\"lbfgs\"])\n",
    "mdf.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is saying that we have 72 pigs with 861 observations between them. The \"Group Var\" here is the random effect intercept for pigs, and it has a variance of 19.5. Meanwhile the slope random effect has a variance of 0.416. In addition to these, there are the two fixed effects. It's easier to see where all of these estimates come out by plotting them. We'll be honest, plotting the estimates hasn't been made easy---so getting the data out and configuring it isn't as easy as it could be.\n",
    "\n",
    "First we pop the random effects results, and the errors, into dataframes in the tidy format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "value_col_titles = [\"Group\", \"Time\"]\n",
    "group_time_ests = pd.DataFrame.from_dict(mdf.random_effects).T\n",
    "group_time_ests = pd.melt(group_time_ests, value_vars=value_col_titles)\n",
    "# add on the parameters\n",
    "for value_col, add_on in zip(value_col_titles, [mdf.fe_params[\"Intercept\"], mdf.fe_params[\"Time\"]]):\n",
    "    group_time_ests.loc[group_time_ests[\"variable\"] == value_col, \"value\"] = group_time_ests.loc[group_time_ests[\"variable\"] == value_col, \"value\"] + add_on\n",
    "\n",
    "group_time_ests.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pig_errors = mdf.conf_int().loc[[\"Intercept\", \"Time\"]]\n",
    "pig_errors.columns = [\"min\", \"max\"]\n",
    "pig_errors[\"mean\"] = pig_errors.mean(axis=1)\n",
    "pig_errors = pig_errors.reset_index().rename(columns={\"index\": \"variable\"}).replace({\"Intercept\": \"Group\"})\n",
    "pig_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "  ggplot() +\n",
    "  geom_errorbar(aes(y='variable', xmin='min', xmax='max'), data=pig_errors, color='black', size=2, height=0.1, alpha=1) +\n",
    "  geom_point(aes(y='variable', x=\"mean\"), data=pig_errors, color=\"black\", size=5) +\n",
    "  geom_jitter(aes(y=\"variable\", x=\"value\"), data=group_time_ests, height=0.2, size=2, alpha=0.6, color=\"blue\") +\n",
    "  ggtitle(\"Random effects: pig model\") +\n",
    "  xlim(5, 30) +\n",
    "  labs(\n",
    "      x=\"Estimate\", y=\"Variable\",\n",
    "  )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Violations of the classical linear model (CLM)\n",
    "\n",
    "### Heteroskedasticity\n",
    "\n",
    "If an estimated model is homoskedastic then its random variables have equal (finite) variance. This is also known as homogeneity of variance. Another way of putting it is that, for all *observations* $i$ in an estimated model $y_i = X_i\\hat{\\beta} + \\epsilon_i$ then\n",
    "\n",
    "$$\n",
    "\\mathbb{E}(\\epsilon_i \\epsilon_i) = \\sigma^2\n",
    "$$\n",
    "\n",
    "When this relationship does not hold, an estimated model is said to be heteroskedastic.\n",
    "\n",
    "To test for heteroskedasticity, you can use **statsmodels**' versions of the [Breusch-Pagan](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.het_breuschpagan.html#statsmodels.stats.diagnostic.het_breuschpagan) or [White](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.het_white.html#statsmodels.stats.diagnostic.het_white) tests with the null hypothesis that the estimated model is homoskedastic. If the null hypothesis is rejected, then standard errors, t-statistics, and F-statistics are invalidated. In this case, you will need HAC (heteroskedasticity and auto-correlation consistent) standard errors, t- and F-statistics.\n",
    "\n",
    "To obtain HAC standard errors from existing regression results in a variable `results`, you can use (for 1 lag):\n",
    "\n",
    "```python\n",
    "results.get_robustcov_results('HAC', maxlags=1).summary()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantile regression\n",
    "\n",
    "Quantile regression estimates the conditional quantiles of a response variable. In some cases, it can be more robust to outliers and, in the case of the $q=0.5$ quantile it is equivalent LAD (Least Absolute Deviation) regression. Let's look at an example of quantile regression in action, lifted direct from the **statsmodels** [documentation](https://www.statsmodels.org/dev/examples/notebooks/generated/quantile_regression.html) and based on a Journal of Economic Perspectives paper by Koenker and Hallock {cite:ps}`koenker2001quantile`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = sm.datasets.engel.load_pandas().data\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have here are two sets of related data. Let's perform several quantile regressions from 0.1 to 0.9 in steps of 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = smf.quantreg(\"foodexp ~ income\", df)\n",
    "quantiles = np.arange(0.1, 1.0, 0.1)\n",
    "q_results = [mod.fit(q=x) for x in quantiles]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $q=0.5$ entry will be at the `4` index; let's take a look at it: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(q_results[4].summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the results for all of the regressions *and* let's add in OLS for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_res = smf.ols(\"foodexp ~ income\", df).fit()\n",
    "\n",
    "get_y = lambda a, b: a + b * x\n",
    "x = np.arange(df.income.min(), df.income.max(), 50)\n",
    "# Just to make the plot clearer\n",
    "x_max = 3000\n",
    "x = x[x < x_max]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "df.plot.scatter(\n",
    "    ax=ax, x=\"income\", y=\"foodexp\", alpha=0.7, s=10, zorder=2, edgecolor=None\n",
    ")\n",
    "for i, res in enumerate(q_results):\n",
    "    y = get_y(res.params[\"Intercept\"], res.params[\"income\"])\n",
    "    ax.plot(x, y, color=\"grey\", lw=0.5, zorder=0, linestyle=(0, (5, 10)))\n",
    "    ax.annotate(f\"$q={quantiles[i]:1.1f}$\", xy=(x.max(), y.max()))\n",
    "y = get_y(ols_res.params[\"Intercept\"], ols_res.params[\"income\"])\n",
    "ax.plot(x, y, color=\"red\", label=\"OLS\", zorder=0)\n",
    "ax.legend()\n",
    "ax.set_xlim(0, x_max)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chart shows very clearly how quantile regression differs from OLS. The line fitted by OLS is trying to be all things to all points whereas the line fitted by quantile regression is focused only on its quantile. You can also see how points far from the median (not all shown) may be having a large influence on the OLS line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rolling and recursive regressions\n",
    "\n",
    "Rolling ordinary least squares applies OLS (ordinary least squares) across a fixed window of observations and then rolls (moves or slides) that window across the data set. They key parameter is `window`, which determines the number of observations used in each OLS regression. Recursive regression is equivalent to rolling regression but with a window that expands over time.\n",
    "\n",
    "Let's first create some synthetic data to perform estimation on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.regression.rolling import RollingOLS\n",
    "import statsmodels.api as sm\n",
    "from sklearn.datasets import make_regression\n",
    "\n",
    "X, y = make_regression(n_samples=200, n_features=2, random_state=0, noise=4.0, bias=0)\n",
    "df = pd.DataFrame(X).rename(columns={0: \"feature0\", 1: \"feature1\"})\n",
    "df[\"target\"] = y\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's fit the model using a formula and a `window` of 25 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roll_reg = RollingOLS.from_formula(\n",
    "    \"target ~ feature0 + feature1 -1\", window=25, data=df\n",
    ")\n",
    "model = roll_reg.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that -1 in the formala suppresses the intercept. We can see the parameters using `model.params`. Here are the params for time steps between 20 and 30:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.params[20:30]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there aren't parameters for entries between 0 and 23 because our window is 25 steps wide. We can easily look at how any of the coefficients are changing over time. Here's an example for 'feature0'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = model.plot_recursive_coefficient(variables=[\"feature0\"])\n",
    "plt.xlabel(\"Time step\")\n",
    "plt.ylabel(\"Coefficient value\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A rolling regression with an *expanding* rather than *moving* window is effectively a recursive least squares model. We can do this instead using the `RecursiveLS` function from **statsmodels**. Let's fit this to the whole dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_rls = sm.RecursiveLS.from_formula(\"target ~ feature0 + feature1 -1\", df)\n",
    "model_rls = reg_rls.fit()\n",
    "print(model_rls.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now we can look back at how the values of the coefficients changed over time too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = model_rls.plot_recursive_coefficient(\n",
    "    range(reg_rls.k_exog), legend_loc=\"upper right\"\n",
    ")\n",
    "ax_list = fig.axes\n",
    "for ax in ax_list:\n",
    "    ax.set_xlim(0, None)\n",
    "ax_list[-1].set_xlabel(\"Time step\")\n",
    "ax_list[0].set_title(\"Coefficient value\");"
   ]
  }
 ],
 "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
}