Tutorials/PythonRootMin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root Finders and Minimizers in Python\n",
    "### by [Richard W. Evans](https://sites.google.com/site/rickecon/), July 2019\n",
    "\n",
    "The code in this Jupyter notebook was written using Python 3.7. It also uses data file `Econ381totpts.txt`. For the code to run properly, you will need every code cell in this notebook that references the data file to have the correct directory path listed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. General characterization of maximization/minimization problem\n",
    "Root finders and minimizers are the two key tools for solving numerical optimization problems. We will define the general formulation of an optimization problem as a minimization problem in the following way,\n",
    "\n",
    "$$ \\min_{x}\\: f\\left(x,z|\\theta\\right) \\quad\\text{s.t}\\quad g(x,z|\\theta)\\geq 0 \\quad\\text{and}\\quad h(x) = 0 $$\n",
    "\n",
    "where $f$ is a system of potentially nonlinear equations that are a function of the vectors of (potentially dynamic)  variables $x$ and $z$ and parameter vector $\\theta$, subject to the vector of inequality constraints $g$ and the vector of equality constraints $h$.\n",
    "\n",
    "A computational algorithm that searches for the value of $x$ that minimizes the problem above is called a **minimizer**. Sometimes the solution to the minimization problem above can be written as a system of equations in $x$.\n",
    "\n",
    "$$ \\hat{x} = x: \\quad \\phi(x|z,\\theta) = 0 $$\n",
    "\n",
    "The maximization problem can be reduced to this system of characterizing equations when the inequality constraints can be shown to never bind (interior solution). A computational algorithm that searches for the value $\\hat{x}$ that sets the value of each equation in the system $\\phi(x|z)$ to 0 is called a **root finder**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Root Finders\n",
    "Suppose the solution to a system of $N$ nonlinear equations $\\phi(x|z,\\theta)=0$ is the $N\\times 1$ vector $\\hat{x}$. \n",
    "\n",
    "\n",
    "### 2.1. Univariate root finders\n",
    "Suppose that $x$ is a scalar such that the system of equations is one equation $\\phi(x|\\theta)=0$ and one unknown $x$. Let's assume the simple cubic polynomial functional form.\n",
    "\n",
    "$$\\phi(x|\\theta)= x^3 + x + \\theta $$\n",
    "\n",
    "The roots of this function are all the values of $x$ such that $\\phi(x|\\theta)=0$, or rather $x^3 + x + \\theta$. We first define a Python function `cubic_pol()` that takes as inputs a value for $x$ and a value for $\\theta$ and returns the value of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi_pol(xvals, theta):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function returns the value phi(x,theta) given x and theta,\n",
    "    where phi(x,theta) = (x ** 3) + x + theta\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals = scalar or (N,) vector, value or values for x\n",
    "    args  = length 1 tuple, (theta)\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    theta = scalar, constant in the phi function\n",
    "    phi   = scalar or (N,) vector, value of phi(xvals, theta)\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: phi\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    phi = (xvals ** 3) + xvals + theta\n",
    "    \n",
    "    return phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can explore this function $\\phi(x|\\theta)$ graphically by plotting it for a given value of $\\theta$ and for many values of $x$. The following plot is for $N=100$ values of $x$ between -10 and 10 and for $\\theta=10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "# This next command is specifically for Jupyter Notebook\n",
    "%matplotlib notebook\n",
    "\n",
    "theta = 10\n",
    "xmin = -10\n",
    "xmax = 10\n",
    "N = 100\n",
    "xvals = np.linspace(xmin, xmax, N)\n",
    "phi_vals = phi_pol(xvals, theta)\n",
    "\n",
    "# Plot the resulting phi values\n",
    "fig, ax = plt.subplots()\n",
    "plt.plot(xvals, phi_vals)\n",
    "minorLocator = MultipleLocator(1)\n",
    "ax.xaxis.set_minor_locator(minorLocator)\n",
    "plt.grid(b=True, which='major', color='0.65', linestyle='-')\n",
    "plt.title('Phi function for theta=10 and x=[-10,10]', fontsize=20)\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel(r'$\\phi(x|\\theta=10)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the plot above, it looks like the function $\\phi(x|\\theta=10)$ has one root somewhere in the neighborhood of $x\\approx -2$. Plugging `x=2` into the function reveals that it is identically the root. So we know what the answer should be. Let's see how close the root finding algorithms get."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python's `Scipy` library has a sub-library called `scipy.optimize`. The `scipy.optimize` sub-library has all of Python's most standard univariate and multivariate root finder algorithms as well as Python's most standard minimizer algorithms. Scrolling a little more than halfway down the main [`scipy.optimize` reference page](https://docs.scipy.org/doc/scipy/reference/optimize.html), we find the links to the documentation for the 5 univariate root finder algorithms in `scipy.optimize`. These are:\n",
    "\n",
    "* [`brentq()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq)\n",
    "* [`brenth()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brenth.html#scipy.optimize.brenth)\n",
    "* [`ridder()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.ridder.html#scipy.optimize.ridder)\n",
    "* [`bisect()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html#scipy.optimize.bisect)\n",
    "* [`newton()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html#scipy.optimize.newton)\n",
    "\n",
    "The first four univariate root finders require a known bracketing interval for the root. That is, the user must know that $\\hat{x}\\in[a, b]$. For the `newton()` method, the bracketing interval is not required. We will solve for this root with the `bisect()` method and with the `newton()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize as opt\n",
    "\n",
    "a =-5\n",
    "b = 0\n",
    "\n",
    "(xhat, result) = opt.bisect(phi_pol, a, b, args=(theta,), full_output=True)\n",
    "\n",
    "print('xhat: ', xhat)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we saw in the Figure of the $\\phi(x|\\theta=10)$ function, the root for this function is very close to $\\hat{x}\\approx 2$.\n",
    "\n",
    "Note that the first argument of the `scipy.optimize.bisect()` function is the function `phi_pol()`. This function's first input must be the scalar being chosen ($x$ in this case), and it must return a scalar representing the function evaluated at a particular $x$. It is the return scalar of the function `phi_pol()` that the `bisect()` algorithm is trying to set to zero. The next two arguments $a$ and $b$ are the bounds within which the `bisect()` algorithm is searching for the root. One requirement of the `bisect()` algorithm is that the target function must return values with different signs when evaluated at $a$ and $b$. The `args=()` argument passes a tuple of extra arguments to the `phi_pol()` function. In this case, we want to pass the `theta` parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One caveat with the bisection method is that it will only find one root in the bracketing interval $[a,b]$. If the function has two roots in the interval, the `bisect()` method will only find one of them. This is true of all five univariate root finder methods as well as all minimizers and multivariate root finders in `scipy.optimize`. They will only return one root of minimum. Whether that root is unique or whether the minimum is a global min rather than a local min depends on how well the researcher understands the function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `newton()` method uses a standard Newton-Raphson hill climbing (gradient descent) method to find the root of the function. Rather than needing to know the bracketing interval where the root exists, you must supply the function with an initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_init = 10.0\n",
    "\n",
    "xhat = opt.newton(phi_pol, x_init, args=(theta,))\n",
    "print(xhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One thing to note here is that the precision of the tolerance matters. That is, I can get a slightly different answer if I start the search from `x_init=-10`. In fact, when I start below the actual value I get the analytical solution $\\hat{x}=-2.0$. However, the difference between the two answers is extremely small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Multivariate root finders\n",
    "For our example of multivariate root finders, let's use the two characterizing equations for optimal lifetime savings from the 3-period-lived agent overlapping generations model with exogenous labor supply. In this case, the solution $x$ will be a $(2\\times 1)$ vector of savings for middle age $b_2$ and old age $b_3$ given interest rates ($r_2, r_3$) and wages ($w_1, w_2, w_3$) over the lifetime and utility function parameters ($\\beta, \\sigma$).\n",
    "\n",
    "$$ \\phi(b_2, b_3|r_2, r_3, w_1, w_2, w_3, \\beta, \\sigma) = 0 $$\n",
    "\n",
    "$$ \\Rightarrow \\qquad\\:\\:\\:\\, (w_1 - b_2)^{-\\sigma} - \\beta(1 + r_2)([1 + r_2]b_2 + w_2 - b_3)^{-\\sigma} = 0 $$\n",
    "$$ \\Rightarrow ([1 + r_2]b_2 + w_2 - b_3)^{-\\sigma} - \\beta(1 + r_3)([1 + r_3]b_3 + w_3)^{-\\sigma} = 0 $$\n",
    "\n",
    "This is two equations and two unknowns ($b_2, b_3$). Everything other than the two savings levels is known in the two equations. Further, Inada conditions bound the solution for consumption strictly away from $c_t\\leq 0$, which means that the two equation system should be solvable using an unconstrained root finder for ($b_2, b_3$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first write some functions that compute the Euler errors shown above given the values for lifetime savings ($b_2, b_3$), interest rates ($r_2, r_3$), wages ($w_1, w_2, w_3$), and utility function parameters ($\\beta, \\sigma$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ct(bt, btp1, rt, wt):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function returns the value of current period consumption given\n",
    "    b_t, b_{t+1}, r_t, and w_t.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    bt   = scalar, current period savings\n",
    "    btp1 = scalar, savings for next period\n",
    "    rt   = scalar > 0, current period interest rate\n",
    "    wt   = scalar > 0, current period wage\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    ct = scalar, current period consumption\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: ct\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    ct = (1 + rt) * bt + wt - btp1\n",
    "    \n",
    "    return ct\n",
    "\n",
    "\n",
    "def MU_stitch(ct, sigma):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    Generate marginal utility of consumption with CRRA consumption\n",
    "    utility and stitched function at lower bound such that the new\n",
    "    hybrid function is defined over all consumption on the real\n",
    "    line but the function has similar properties to the Inada condition.\n",
    "\n",
    "    u'(c) = c ** (-sigma) if c >= epsilon\n",
    "          = g'(c) = 2 * b2 * c + b1 if c < epsilon\n",
    "\n",
    "        such that g'(epsilon) = u'(epsilon)\n",
    "        and g''(epsilon) = u''(epsilon)\n",
    "\n",
    "        u(c) = (c ** (1 - sigma) - 1) / (1 - sigma)\n",
    "        g(c) = b2 * (c ** 2) + b1 * c + b0\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    ct    = scalar, current consumption\n",
    "    sigma = scalar >= 1, coefficient of relative risk aversion for CRRA\n",
    "            utility function: (c**(1-sigma) - 1) / (1 - sigma)\n",
    "\n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "\n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    epsilon  = scalar > 0, positive value close to zero\n",
    "    c_cnstr  = boolean, =True if ct < epsilon\n",
    "    b1       = scalar, intercept value in linear marginal utility\n",
    "    b2       = scalar, slope coefficient in linear marginal utility\n",
    "    MU_c     = scalar, marginal utility of consumption\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "\n",
    "    RETURNS: MU_c\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    epsilon = 0.0001\n",
    "    c_cnstr = ct < epsilon\n",
    "    if c_cnstr:\n",
    "        b2 = (-sigma * (epsilon ** (-sigma - 1))) / 2\n",
    "        b1 = (epsilon ** (-sigma)) - 2 * b2 * epsilon\n",
    "        MU_c = 2 * b2 * ct + b1\n",
    "    else:\n",
    "        MU_c = ct ** (-sigma)\n",
    "    \n",
    "    return MU_c\n",
    "\n",
    "\n",
    "def get_EulErrs(bvec, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    bvec = (2,) vector, values for lifetime savings (b2, b3)\n",
    "    args = length 7 tuple, (r2, r3, w1, w2, w3, beta, sigma)\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        get_ct()\n",
    "        MU_stitch()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    r2    = scalar > 0, interest rate in period 2\n",
    "    r3    = scalar > 0, interest rate in period 3\n",
    "    w1    = scalar > 0, wage in period 1\n",
    "    w2    = scalar > 0, wage in period 2\n",
    "    w3    = scalar > 0, wage in period 3\n",
    "    beta  = scalar in (0, 1), discount factor\n",
    "    sigma = scalar >= 1, coefficient of relative risk aversion\n",
    "    c1    = scalar, consumption in period 1\n",
    "    c2    = scalar, consumption in period 2\n",
    "    c3    = scalar, consumption in period 3\n",
    "    MU_c1 = scalar > 0, marginal utility of consumption in period 1\n",
    "    MU_c2 = scalar > 0, marginal utility of consumption in period 2\n",
    "    MU_c3 = scalar > 0, marginal utility of consumption in period 3\n",
    "    err1  = scalar, Euler error for savings decision b2\n",
    "    err2  = scalar, Euler error for savings decision b3\n",
    "    err_vec = (2,) vector, Euler errors from two Euler equations\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: err_vec\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    b2, b3 = bvec\n",
    "    r2, r3, w1, w2, w3, beta, sigma = args\n",
    "    c1 = get_ct(0.0, b2, 0.0, w1)\n",
    "    c2 = get_ct(b2, b3, r2, w2)\n",
    "    c3 = get_ct(b3, 0.0, r3, w3)\n",
    "    MU_c1 = MU_stitch(c1, sigma)\n",
    "    MU_c2 = MU_stitch(c2, sigma)\n",
    "    MU_c3 = MU_stitch(c3, sigma)\n",
    "    err1 = MU_c1 - beta * (1 + r2) * MU_c2\n",
    "    err2 = MU_c2 - beta * (1 + r3) * MU_c3\n",
    "    err_vec = np.array([err1, err2])\n",
    "    \n",
    "    return err_vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now run a multivariate root finder to try and find the solution to this problem. The function [`scipy.optimize.root()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html#scipy.optimize.root) contains all of the multivariate root finders in `SciPy`. Within `scipy.optimize.root()`, you can choose among 10 different methods of multivariate root finders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declare parameter values and exogenous variable values\n",
    "beta = 0.96 ** (80 / 3)\n",
    "sigma = 2.1\n",
    "r2 = 0.6\n",
    "r3 = 0.7\n",
    "w1 = 0.5\n",
    "w2 = 0.5\n",
    "w3 = 0.5\n",
    "\n",
    "# Make initial guess for solution of savings values. Note that these\n",
    "# two guesses must be feasible and not violate c_t > 0 for all t\n",
    "b2_init = 0.05\n",
    "b3_init = 0.10\n",
    "b_init = np.array([b2_init, b3_init])\n",
    "b_args = (r2, r3, w1, w2, w3, beta, sigma)\n",
    "b_result = opt.root(get_EulErrs, b_init, args=(b_args))\n",
    "print(b_result)\n",
    "print('Roots: ', b_result.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output from the multivariate minimizer above is stored as a Python dictionary named `b_result`. The elements of the dictionary include:\n",
    "\n",
    "* `fun`: the value of the error functions at the solution\n",
    "* `message`: a message about whether the solution converged\n",
    "* `nfev`: number of function evaluations until convergence\n",
    "* `success`: boolean that you want to equal `True`\n",
    "* `x`: the vector of roots of the function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Minimizers\n",
    "In the `scipy.optimize` library is the [`scipy.optimize.minimize()`](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.minimize.html) function, which contains all of `SciPy`'s main unconstrained and constrained minimizer functions. A minimizer is an algorithm that tries to find the $x$ that minimizes the following problem,\n",
    "\n",
    "$$ \\min_{x}\\: f\\left(x,z|\\theta\\right) \\quad\\text{s.t}\\quad g(x,z|\\theta)\\geq 0 \\quad\\text{and}\\quad h(x) = 0 $$\n",
    "\n",
    "where $f$ is a system of potentially nonlinear equations that are a function of the vectors of (potentially dynamic)  variables $x$ and $z$ and parameter vector $\\theta$, subject to the vector of inequality constraints $g$ and the vector of equality constraints $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Unconstrained minimizers\n",
    "An unconstrained minimizer is an algorithm in which no inequality constraints $g$ and no equality constraints $h$ are present.\n",
    "\n",
    "$$ \\min_{x}\\: f\\left(x,z|\\theta\\right) $$\n",
    "\n",
    "Each of the 12 methods available in the [`scipy.optimize.minimize()`](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.minimize.html) function can be used as an unconstrained minimizer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To present an example of using a minimizer, we will perform a standard maximum likelihood estimation. Import some data from the total points earned by all the students in two sections of my intermediate macroeconomics class for undergraduates at my previous University in a certain year (two semesters). The maximum possible points in both semesters was 450."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts = np.loadtxt('DataFiles/Econ381totpts.txt')\n",
    "# pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a histogram of the test score data to see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count, bins, ignored = plt.hist(pts, 30, normed=True, edgecolor='black')\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the fact that this distribution is truncated in the upper tail at `pnts=450` and would be better modeled by a truncated distribution like a truncated normal, we will fit a normal distribution $N(\\mu,\\sigma)$ to this data for simplicity of exposition.\n",
    "\n",
    "The estimation process is a basic minimization problem as stated elsewhere in this notebook. The maximum likelihood estimation problem is to choose the vector of parameters $(\\mu,\\sigma)$ that maximizes the likelihood of seeing the data in `Econ381totpts.txt`, given the assumption that the data came from a normal distribution $N(\\mu,\\sigma)$. Let the probability density function of the normal distribution be represented by $f(x|\\mu,\\sigma)$. This `pdf` gives the probability density of an observation of data $x$ given values for the distribution parameters ($\\mu,\\sigma$). Then the following is the statement of the maximum likelihood estimation as a maximization problem,\n",
    "\n",
    "$$ \\max_{\\mu,\\sigma}\\: \\ln\\bigl(\\mathcal{L}\\bigr) = \\sum_{n=1}^N\\ln\\Bigl(f(pnts_n|\\mu,\\sigma)\\Bigr) $$\n",
    "\n",
    "where $\\ln(\\mathcal{L})$ is the log likelihood function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve this computationally, it must be restated as a minimization problem. This is as simple as minimizing the negative of the log likelihood function.\n",
    "\n",
    "$$ \\min_{\\mu,\\sigma}\\: -\\ln\\bigl(\\mathcal{L}\\bigr) = -\\sum_{n=1}^N\\ln\\Bigl(f(pnts_n|\\mu,\\sigma)\\Bigr) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as sts\n",
    "\n",
    "def log_lik_norm(xvals, mu, sigma):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    Compute the log likelihood function for data xvals given normal\n",
    "    distribution parameters mu and sigma.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, values of the normally distributed random\n",
    "             variable\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        sts.norm.pdf()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    pdf_vals    = (N,) vector, normal PDF values for mu and sigma\n",
    "                  corresponding to xvals data\n",
    "    ln_pdf_vals = (N,) vector, natural logarithm of normal PDF values\n",
    "                  for mu and sigma corresponding to xvals data\n",
    "    log_lik_val = scalar, value of the log likelihood function\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: log_lik_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    pdf_vals = sts.norm.pdf(xvals, mu, sigma)\n",
    "    ln_pdf_vals = np.log(pdf_vals)\n",
    "    log_lik_val = ln_pdf_vals.sum()\n",
    "    \n",
    "    return log_lik_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.1. The criterion function\n",
    "Now that we have written a function `log_lik_norm()` that gives us the value of the log likelihood function $\\ln(\\mathcal{L})$ for a given set of data $\\{pnts_n\\}_{n=1}^N$ and for given parameter values $\\mu$ and $\\sigma$, we must define a criterion function that will be the object of our minimization algorithm.\n",
    "\n",
    "The first step in writing a criterion function is to define a function that takes two inputs and returns a scalar value.\n",
    "1. The first input is either a scalar or a vector of values (the object `params` in the function `crit()` below). This object is the value or values being chosen to minimize the criterion function.\n",
    "2. The second object is Python's variable length input objects `*args`, which is a tuple of variable length positional arguments. As you will see in the `minimize()` function, all the arguments must be passed into the criterion function in one tuple.\n",
    "3. Lastly, you must make sure that the scalar criterion value that the function returns is the value of the problem stated as a minimization problem and not a maximization problem. In this case of maximum likelihood estimation, you want the negative of the log likelihood function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crit(params, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the negative of the log likelihood function\n",
    "    given parameters and data. This is the minimization problem version\n",
    "    of the maximum likelihood optimization problem\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    params = (2,) vector, ([mu, sigma])\n",
    "    args   = length 1 tuple, (xvals,)\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        log_lik_norm()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mu              = scalar, mean of the normally distributed random\n",
    "                      variable\n",
    "    sigma           = scalar > 0, standard deviation of the normally\n",
    "                      distributed random variable\n",
    "    xvals           = (N,) vector, values of the normally distributed\n",
    "                      random variable\n",
    "    log_lik_val     = scalar, value of the log likelihood function\n",
    "    neg_log_lik_val = scalar, negative of log_lik_val\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: neg_log_lik_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mu, sigma = params\n",
    "    xvals = args\n",
    "    log_lik_val = log_lik_norm(xvals, mu, sigma)\n",
    "    neg_log_lik_val = -log_lik_val\n",
    "    \n",
    "    return neg_log_lik_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.2. The minimize() function\n",
    "The `minimize()` function is shorthand for `scipy.optimize.minimize()`. This function returns a dictionary of objects including the solution to the optimization problem and whether the problem actually solved. The `minimize` function has three mandatory arguments, plus a lot of options. You can experiment with the options on the [`minimize()` documentation page](https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize).\n",
    "1. The first argument of the minimize function is the criterion function (`crit()` in this example) from which the `minimize()` function will test values of the parameters in searching for the minimum value.\n",
    "2. The second argument is an initial guess for the values of the parameters that minimize the criterion function `crit()`.\n",
    "3. The third argument is the tuple of all the objects needed to solve the criterion function in `crit()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_init = 400\n",
    "sig_init = 70\n",
    "params_init = np.array([mu_init, sig_init])\n",
    "mle_args = (pts,)\n",
    "norm_results_uncstr = opt.minimize(crit, params_init, args=(mle_args))\n",
    "mu_MLE, sig_MLE = norm_results_uncstr.x\n",
    "print('mu_MLE=', mu_MLE, ' sig_MLE=', sig_MLE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_results_uncstr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how this estimated distribution with $\\hat{\\mu}$ and $\\hat{\\sigma}$ fits the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "count, bins, ignored = plt.hist(pts, 30, normed=True, edgecolor='black')\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\n",
    "\n",
    "# Plot the MLE estimated distribution\n",
    "dist_pts = np.linspace(0, 550, 500)\n",
    "plt.plot(dist_pts, sts.norm.pdf(dist_pts, mu_MLE, sig_MLE),\n",
    "         linewidth=2, color='k', label='$\\mu$=342,$\\sigma$=88')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As was mentioned before, this distribution does not have an ideal fit because we assumed a normal distribution (non-truncated), whereas the data clearly had an upper bound at 450. However, this section demonstrates how to successfully execute a minimizer using the `scipy.optimize.minimize()` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Constrained minimizers\n",
    "The `minimize()` function has many methods that can be used to find the parameter values that minimize some criterion function. These methods are called using the `method='MethodName'` optional input argument to the minimize function. Three of those methods allow for constrained minimization by providing upper and lower bounds for the parameters being chosen. These three methods are `'L-BFGS-B'`, `'TNC'`, and `'SLSQP'`.\n",
    "\n",
    "Suppose you were trying to estimate $\\mu$ and $\\sigma$ of a normal distribution as is the case in our test scores example. The value of $\\mu$ need not be constrained. However, the value of $\\sigma$ must be strictly positive. You could include these bounds in a constrained minimization by using the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_results_cstr = opt.minimize(crit, params_init, args=(mle_args), method='L-BFGS-B',\n",
    "                                 bounds=((None, None), (1e-10, None)))\n",
    "mu_MLE_cstr, sig_MLE_cstr = norm_results_cstr.x\n",
    "print('mu_MLE=', mu_MLE_cstr, ' sig_MLE=', sig_MLE_cstr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_results_cstr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you must set the lower bound of $\\sigma$ equal to some small positive number close to zero. You cannot set it to zero itself because the bounds are inclusive. That is, the minimizer might try a value of $\\sigma=0$ is the lower bound includes zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Final Notes\n",
    "The little advertised truth is that no root finder is guaranteed to find the right roots or any roots, and no minimizer is guaranteed to find the global minimum for any set of characterizing equations. More accurately, it is a rare and special case to find characterizing functions for which standard root finders and minimizers solve robustly. Press, et al (2007, pp. 442-443) highlight that with nonlinear multidimensional root finding problems,\n",
    "\n",
    "> \"...you can never be sure that the root is there at all until you have found it.... It cannot be overemphasized, however, how crucially success depends on having a good first guess for the solution,....\"\n",
    "\n",
    "Press, et al (2007, p. 473) also state,\n",
    "\n",
    "> \"We make an extreme, but wholly defensible statement: There are no good, general methods for solving systems of more than one nonlinear equation. Furthermore, it is not hard to see why (very likely) there never will be any good, general methods.\"\n",
    "\n",
    "Because initial values are so important to root finders and minimization problems, Judd (1998, p. 172) suggests running a minimization problem with a loose stopping rule on the vector of squared errors. Then one can use that solution as the initial guess for the root finder.\n",
    "\n",
    "Because root finding and minimization is so difficult, we take great care to report the errors and optimization output in our code. One must show not only that the minimizer or root finder stopped iterating, but also that the errors in the characterizing equations are arbitrarily small. One can also perform robustness tests using varying initial values to provide evidence that the solution is unique in the neighborhood of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "* Judd, Kenneth L., *Numerical Methods in Economics*, MIT Press (1998).\n",
    "* Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery, *Numerical Recipes: The Art of Scientific Computing*, 3rd edition, Cambridge University Press (2007)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}