{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Likelihood estimation\n",
    "\n",
    "This notebook shows how to do a simple maximum likelihood (ml) estimation with estimagic. As an illustrating example, we implement a simple linear regression model. This is the same example model used as in the method of moments notebook.\n",
    "\n",
    "We proceed in 4 steps:\n",
    "\n",
    "\n",
    "1. Create a data generating process\n",
    "2. Set up a likelihood function\n",
    "3. Maximize the likelihood function\n",
    "4. Calculate standard errors, confidence intervals, and p-values\n",
    "\n",
    "The user only needs to do step 1 and 2. The rest is done by `estimate_ml`. \n",
    "\n",
    "To be very clear: Estimagic is not a package to estimate linear models or other models that are implemented in Stata, statsmodels or anywhere else. Its purpose is to estimate parameters with custom likelihood or method of simulated moments functions. We just use an ordered logit model as an example of a very simple likelihood function.\n",
    "\n",
    "\n",
    "### Model:\n",
    "\n",
    "$$ y = \\beta_0 + \\beta_1 x + \\epsilon, \\text{ where } \\epsilon \\sim N(0, \\sigma^2)$$\n",
    "\n",
    "We aim to estimate $\\beta_0, \\beta_1, \\sigma^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from patsy import dmatrices\n",
    "from scipy import stats\n",
    "from scipy.stats import norm\n",
    "\n",
    "import estimagic as em\n",
    "\n",
    "rng = np.random.default_rng(seed=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Create a data generating process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_data(params, n_draws):\n",
    "    x = rng.normal(0, 1, size=n_draws)\n",
    "    e = rng.normal(0, params.loc[\"sd\", \"value\"], size=n_draws)\n",
    "    y = params.loc[\"intercept\", \"value\"] + params.loc[\"slope\", \"value\"] * x + e\n",
    "    return pd.DataFrame({\"y\": y, \"x\": x})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>lower_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>2</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>slope</th>\n",
       "      <td>-1</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>1</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           value   lower_bound\n",
       "intercept      2          -inf\n",
       "slope         -1          -inf\n",
       "sd             1  1.000000e-10"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_params = pd.DataFrame(\n",
    "    data=[[2, -np.inf], [-1, -np.inf], [1, 1e-10]],\n",
    "    columns=[\"value\", \"lower_bound\"],\n",
    "    index=[\"intercept\", \"slope\", \"sd\"],\n",
    ")\n",
    "true_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = simulate_data(true_params, n_draws=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Define the `loglike` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normal_loglike(params, data):\n",
    "    norm_rv = norm(\n",
    "        loc=params.loc[\"intercept\", \"value\"] + params.loc[\"slope\", \"value\"] * data[\"x\"],\n",
    "        scale=params.loc[\"sd\", \"value\"],\n",
    "    )\n",
    "    contributions = norm_rv.logpdf(data[\"y\"])\n",
    "    return {\"contributions\": contributions, \"value\": contributions.sum()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few remarks before we move on:\n",
    "\n",
    "1. There are numerically better ways to calculate the likelihood; we chose this implementation for brevity and readability. \n",
    "2. The loglike function takes params and other arguments. You are completely flexible with respect to the number and names of the other arguments as long as the first argument is params. \n",
    "3. The loglike function returns a dictionary with the entries \"contributions\" and \"value\". The \"contributions\" are the log likelihood evaluations of each individual in the dataset. The \"value\" are their sum. The \"value\" entry could be omitted, the \"contributions\" entry, however, is mandatory. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Estimate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_params = true_params.assign(value=[100, 100, 100])\n",
    "\n",
    "res = em.estimate_ml(\n",
    "    loglike=normal_loglike,\n",
    "    params=start_params,\n",
    "    optimize_options={\"algorithm\": \"scipy_lbfgsb\"},\n",
    "    loglike_kwargs={\"data\": data},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>standard_error</th>\n",
       "      <th>ci_lower</th>\n",
       "      <th>ci_upper</th>\n",
       "      <th>p_value</th>\n",
       "      <th>free</th>\n",
       "      <th>stars</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>1.945</td>\n",
       "      <td>0.104</td>\n",
       "      <td>1.742</td>\n",
       "      <td>2.148</td>\n",
       "      <td>0.0</td>\n",
       "      <td>True</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>slope</th>\n",
       "      <td>-0.945</td>\n",
       "      <td>0.113</td>\n",
       "      <td>-1.167</td>\n",
       "      <td>-0.723</td>\n",
       "      <td>0.0</td>\n",
       "      <td>True</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>0.954</td>\n",
       "      <td>0.079</td>\n",
       "      <td>0.799</td>\n",
       "      <td>1.109</td>\n",
       "      <td>0.0</td>\n",
       "      <td>True</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           value  standard_error  ci_lower  ci_upper  p_value  free stars\n",
       "intercept  1.945           0.104     1.742     2.148      0.0  True   ***\n",
       "slope     -0.945           0.113    -1.167    -0.723      0.0  True   ***\n",
       "sd         0.954           0.079     0.799     1.109      0.0  True   ***"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. What's in the results?\n",
    "\n",
    "`LikelihoodResult` objects provide attributes and methods to calculate standard errors, confidence intervals, and p-values. For all three, several methods are available. You can even calculate cluster robust standard errors. \n",
    "\n",
    "A few examples are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>lower_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>1.944957</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>slope</th>\n",
       "      <td>-0.944914</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>0.954233</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              value   lower_bound\n",
       "intercept  1.944957          -inf\n",
       "slope     -0.944914          -inf\n",
       "sd         0.954233  1.000000e-10"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>intercept</th>\n",
       "      <th>slope</th>\n",
       "      <th>sd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>0.008986</td>\n",
       "      <td>0.000426</td>\n",
       "      <td>-0.001904</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>slope</th>\n",
       "      <td>0.000426</td>\n",
       "      <td>0.007734</td>\n",
       "      <td>0.000303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>-0.001904</td>\n",
       "      <td>0.000303</td>\n",
       "      <td>0.003748</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           intercept     slope        sd\n",
       "intercept   0.008986  0.000426 -0.001904\n",
       "slope       0.000426  0.007734  0.000303\n",
       "sd         -0.001904  0.000303  0.003748"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.cov(method=\"robust\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>lower_bound</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>0.103760</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>slope</th>\n",
       "      <td>0.113343</td>\n",
       "      <td>-inf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sd</th>\n",
       "      <td>0.078961</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              value   lower_bound\n",
       "intercept  0.103760          -inf\n",
       "slope      0.113343          -inf\n",
       "sd         0.078961  1.000000e-10"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.se()"
   ]
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}