slides_09_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7b8f817c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Optimization\n",
    "\n",
    "### Zhentao Shi\n",
    "\n",
    "<img src=\"graph/chongchicollege.png\" width=\"1000\">\n",
    "\n",
    "<!-- code is tested on SCRP -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c41b618",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(magrittr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "856666bf",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Optimization \n",
    "\n",
    "* Econometrics curriculum does not pay enough attention to numerical optimization\n",
    "* Most estimators solve optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29e2cbc9",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Operational research\n",
    "* Understand the essence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c874e9d0",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### User Cases\n",
    "\n",
    "* Maximum likelihood estimation\n",
    "* Discrete / mixed data type\n",
    "* Machine learning / regularization\n",
    "* Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "998fdf34",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Big data stochastic algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b199d2a7",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "* Generic optimization problem\n",
    "\n",
    "$$\n",
    "\\min_{\\theta \\in \\Theta } f(\\theta) \\,\\, \\mathrm{ s.t. }  g(\\theta) = 0, h(\\theta) \\leq 0,\n",
    "$$\n",
    "\n",
    "* $f(\\cdot)\\in \\mathbb{R}$: criterion function\n",
    "* $g(\\theta) = 0$: a vector of equality constraints\n",
    "* $h(\\theta)\\leq 0$: a vector of inequality constraints.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0081bf3f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* unconstrained\n",
    "* constrained\n",
    "* Lagrangian "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "157b4b97",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convenience vs. Efficiency\n",
    "\n",
    "* Convenience: readability of the mathematical expressions and the code\n",
    "* Efficiency:  computing speed\n",
    "\n",
    "* Put convenience as priority at the trial-and-error stage, \n",
    "* Improve efficiency when necessary at a later stage for full-scale execution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0ae8ac2",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Methods\n",
    "\n",
    "* Many optimization algorithms\n",
    "* Variants of a few fundamental principles.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d36f3d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Newton's Method\n",
    "\n",
    "* Essential idea for optimizing a twice-differentiable objective function\n",
    "* Necessary condition: the first-order condition\n",
    "\n",
    "$$\n",
    "s(\\theta) = \\partial f(\\theta) / \\partial \\theta = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "630674fa",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "f <- function(x) 0.1 * (x - 5)^2 + cos(x) # criterion\n",
    "s <- function(x) 0.2 * (x - 5) - sin(x) # gradient\n",
    "h <- function(x) 0.2 - cos(x) # Hessian\n",
    "\n",
    "# plot\n",
    "par(mfrow = c(3, 1))\n",
    "par(mar = c(2, 4, 1, 2))\n",
    "\n",
    "x_base <- seq(0.1, 10, by = 0.1)\n",
    "plot(y = f(x_base), x = x_base, type = \"l\", lwd = 2, ylab = \"f\")\n",
    "plot(y = s(x_base), x = x_base, type = \"l\", ylab = \"score\")\n",
    "abline(h = 0, lty = 2)\n",
    "plot(y = h(x_base), x = x_base, type = \"l\", ylab = \"Hessian\")\n",
    "abline(h = 0, lty = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a3e511e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Iteartion\n",
    "\n",
    "* Initial trial value $\\theta_0$, \n",
    "* If $s(\\theta_0) \\neq 0$, updated by\n",
    "\n",
    "$$\n",
    "\\theta_{t+1} = \\theta_{t} -  \\left( H(\\theta_t)  \\right)^{-1}  s(\\theta_t)\n",
    "$$\n",
    "\n",
    "for the index of iteration $t=0,1,\\cdots$\n",
    "* $H(\\theta) = \\frac{ \\partial s(\\theta )}{ \\partial \\theta}$ is the Hessian. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "805cb59d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mechanism\n",
    "\n",
    "* Taylor expansion\n",
    "at $\\theta_t$ round  $\\theta_{\\star}$, a root of $s(\\cdot)$. Because $\\theta_{ \\star }$  is a root,\n",
    "\n",
    "$$\n",
    "0 = s(\\theta_{\\star}) = s(\\theta_t) + H(\\theta_t) (\\theta_{\\star} - \\theta_t) + O( (\\theta_{\\star} - \\theta_t)^2 ).\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1daa97b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Update\n",
    "\n",
    "* Ignore the high-order term and rearrange,\n",
    "\n",
    "$$\n",
    "\\theta_{\\star} = \\theta_{t} -  \\left( H(\\theta_t)  \\right)^{-1}  s(\\theta_t)\n",
    "$$ \n",
    "\n",
    "* iteration formula by replacing $\\theta_{\\star}$ with the updated $\\theta_{t+1}$. \n",
    "* Iterate until $|\\theta_{t+1} -\\theta_{t}| < \\epsilon$ (absolute criterion) and/or\n",
    "$|\\theta_{t+1} -\\theta_{t}|/|\\theta_{t}| < \\epsilon$ (relative criterion), \n",
    "* $\\epsilon$ is a small positive number chosen as a tolerance level.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de6ce8c4",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Newton's method\n",
    "Newton <- function(x) {\n",
    "  x - s(x) / h(x)\n",
    "} # update formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2413f934",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "x_init <- 1 # can experiment with various initial values\n",
    "\n",
    "gap <- 1\n",
    "epsilon <- 0.000001 # tolerance\n",
    "while (gap > epsilon) {\n",
    "  x_new <- Newton(x_init) %>% print()\n",
    "  gap <- abs(x_init - x_new)\n",
    "  x_init <- x_new\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5cf0c57",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Features of Newton's Method\n",
    "\n",
    "\n",
    "* It seeks the solution to $s(\\theta) = 0$. \n",
    "* The first-order condition is necessary but not sufficient\n",
    "* Verify the second-order condition\n",
    "* Compare the value of multiple minima for global minimum\n",
    "\n",
    "* It requires gradient $s(\\theta)$ and the Hessian $H(\\theta)$.\n",
    "* It numerically converges at quadratic rate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ab1220f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Quasi-Newton Method\n",
    "\n",
    "* Most well-known quasi-Newton algorithm is [BFGS](http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)\n",
    "* It avoids explicit calculation of Hessian\n",
    "* It starts from an initial (inverse) Hessian\n",
    "* Updates Hessian by an explicit formula via quadratic approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b848391",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Derivative-Free Method\n",
    "\n",
    "* [Nelder-Mead](http://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method)\n",
    "* Simplex method\n",
    "* Search a local minimum \n",
    "  * reflection\n",
    "  * expansion\n",
    "  * contraction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "204ed7d9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "## Package\n",
    "\n",
    "* [R Optimization Task View](http://cran.r-project.org/web/views/Optimization.html) \n",
    "* Package [`optimx`](http://cran.r-project.org/web/packages/optimx/index.html) ([Nash, 2014](https://www.jstatsoft.org/article/view/v060i02))\n",
    "  * Unified interface for various widely-used algorithms. \n",
    "  * Facilitates comparison among optimization algorithms\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dec6362e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "* Use `optimx` to solve pseudo Poisson maximum likelihood estimation (PPML)\n",
    "* Popular estimator for cross-country bilateral trade\n",
    "* Conditional mean model\n",
    "\n",
    "$$\n",
    "E[y_i | x_i] = \\exp( x_i' \\beta),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6abe9de",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Poisson MLE\n",
    "\n",
    "If $Z \\sim Poisson(\\lambda)$, then \n",
    "\n",
    "$$\n",
    "\\Pr(Z = k) = \\frac{\\mathrm{e}^{-\\lambda} \\lambda^k}{k!}, \\mathrm{ for }\\, \\, k=0,1,2,\\ldots,\n",
    "$$\n",
    "\n",
    "and the log-likelihood\n",
    "\n",
    "$$\n",
    "\\log \\Pr(Y = y | x) =  -\\exp(x'\\beta) + y\\cdot x'\\beta - \\log k!\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3925c61",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Log-likelihood function of the sample\n",
    "$\n",
    "\\ell(\\beta) = \\log \\Pr( \\mathbf{y} | \\mathbf{x};\\beta ) =\n",
    "-\\sum_{i=1}^n \\exp(x_i'\\beta) + \\sum_{i=1}^n y_i x_i'\\beta.\n",
    "$\n",
    "\n",
    "* gradient:\n",
    "$s(\\beta) =\\frac{\\partial \\ell(\\beta)}{\\partial \\beta} =\n",
    "-\\sum_{i=1}^n \\exp(x_i'\\beta)x_i + \\sum_{i=1}^n y_i x_i.$\n",
    "\n",
    "* Hessian:\n",
    "$H(\\beta) = \\frac{\\partial^2 \\ell(\\beta)}{\\partial \\beta \\partial \\beta'} =\n",
    "-\\sum_{i=1}^n \\exp(x_i'\\beta)x_i x_i'$\n",
    "\n",
    "is negative definite."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e94f29d8",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* $\\ell(\\beta)$ is strictly concave in $\\beta$.\n",
    "\n",
    "* Default optimization is minimization\n",
    "* Use *negative* log-likelihood\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcecb096",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Poisson likelihood\n",
    "poisson.loglik <- function(b) {\n",
    "  b <- as.matrix(b)\n",
    "  lambda <- exp(X %*% b)\n",
    "  ell <- -sum(-lambda + y * log(lambda))\n",
    "  return(ell)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15921bb",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Write the criterion as a function of the parameter to be optimized \n",
    "* Data can be fed inside or outside of the function.\n",
    "  * If the data is provided as additional arguments, these arguments must be explicit.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6e4501a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# implement both BFGS and Nelder-Mead for comparison.\n",
    "\n",
    "library(AER)\n",
    "\n",
    "## prepare the data\n",
    "data(\"RecreationDemand\")\n",
    "y <- RecreationDemand$trips\n",
    "X <- with(RecreationDemand, cbind(1, income))\n",
    "\n",
    "## estimation\n",
    "b.init <- c(2, 1) # initial value\n",
    "b.hat <- optimx::optimx(b.init, poisson.loglik,\n",
    "  method = c(\"BFGS\", \"Nelder-Mead\"),\n",
    "  control = list(\n",
    "    reltol = 1e-7,\n",
    "    abstol = 1e-7\n",
    "  )\n",
    ")\n",
    "print(b.hat)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ed425ea",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Alternative Formulation\n",
    "\n",
    "* Nonlinear least squares (NLS) is also valid in theory.\n",
    "* NLS minimizes\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^n (y_i - \\exp(x_i \\beta))^2\n",
    "$$\n",
    "\n",
    "* PPML's optimization for the linear index is globally convex.\n",
    "* Numerical optimization of PPML is easier and more robust"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e502f7e0",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Caveats\n",
    "\n",
    "* No algorithm suits all problems. \n",
    "* Simulation is helpful to check the accuracy of optimization\n",
    "* Contour plot helps visualize the function surface/manifold"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2464c1d1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e534db3",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "x.grid <- seq(0, 1.8, 0.02)\n",
    "x.length <- length(x.grid)\n",
    "y.grid <- seq(-.5, .2, 0.01)\n",
    "y.length <- length(y.grid)\n",
    "\n",
    "z.contour <- matrix(0, nrow = x.length, ncol = y.length)\n",
    "\n",
    "for (i in 1:x.length) {\n",
    "  for (j in 1:y.length) {\n",
    "    z.contour[i, j] <- poisson.loglik(c(x.grid[i], y.grid[j]))\n",
    "  }\n",
    "}\n",
    "\n",
    "contour(x.grid, y.grid, z.contour, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7964e0b4",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### NLOPT\n",
    "\n",
    "* Third-party standalone solvers \n",
    "* [`NLopt`](http://ab-initio.mit.edu/wiki/index.php/NLopt_Installation)\n",
    "* [extensive list of algorithms](http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms#SLSQP)\n",
    "* Package [`nloptr`](http://cran.r-project.org/web/packages/nloptr/index.html)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba0ec5c2",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "We first carry out the Nelder-Mead algorithm in NLOPT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99b14a64",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "## optimization with NLoptr\n",
    "\n",
    "opts <- list(\n",
    "  \"algorithm\" = \"NLOPT_LN_NELDERMEAD\",\n",
    "  \"xtol_rel\" = 1.0e-7,\n",
    "  maxeval = 500\n",
    ")\n",
    "\n",
    "res_NM <- nloptr::nloptr(\n",
    "  x0 = b.init,  eval_f = poisson.loglik,   opts = opts\n",
    ")\n",
    "print(res_NM)\n",
    "\n",
    "# \"SLSQP\" is indeed the BFGS algorithm in NLopt,\n",
    "# though \"BFGS\" doesn't appear in the name\n",
    "opts <- list(\"algorithm\" = \"NLOPT_LD_SLSQP\", \"xtol_rel\" = 1.0e-7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03769c2d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To invoke BFGS in NLOPT, we must code up the gradient $s(\\beta)$,\n",
    "as in the function `poisson.log.grad()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8e313fa",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "poisson.loglik.grad <- function(b) {\n",
    "  b <- as.matrix(b)\n",
    "  lambda <- exp(X %*% b)\n",
    "  ell <- -colSums(-as.vector(lambda) * X + y * X)\n",
    "  return(ell)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bed9195",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* Compare the analytical gradient with the numerical gradient\n",
    "* Ensure the code is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bea1cf5b",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# check the numerical gradient and the analytical gradient\n",
    "b <- c(0, .5)\n",
    "numDeriv::grad(poisson.loglik, b)\n",
    "poisson.loglik.grad(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ede1a291",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "BFGS with Gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d268d5fc",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "res_BFGS <- nloptr::nloptr(\n",
    "  x0 = b.init,\n",
    "  eval_f = poisson.loglik,\n",
    "  eval_grad_f = poisson.loglik.grad,\n",
    "  opts = opts\n",
    ")\n",
    "print(res_BFGS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f4891ca",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convex Optimization\n",
    "\n",
    "* Local minimum is a global minimum.\n",
    "* Particularly important in high-dimensional problems\n",
    "* [Boyd and Vandenberghe (2004)](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)\n",
    "* \"Convex optimization is technology; all other optimizations are arts.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9c7cfa8",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "f1 <- function(x) x^2\n",
    "f2 <- function(x) abs(x)\n",
    "f3 <- function(x) (-x - 1) * (x <= -1) + (0.4 * x - .2) * (x >= .5)\n",
    "\n",
    "par(mfrow = c(1, 3))\n",
    "par(mar = c(4, 2, 1, 2))\n",
    "\n",
    "x_base <- seq(-3, 3, by = 0.1)\n",
    "plot(y = f1(x_base), x = x_base, type = \"l\", lwd = 2, xlab = \"differentiable\")\n",
    "plot(y = f2(x_base), x = x_base, type = \"l\", lwd = 2, xlab = \"non-differentiable\")\n",
    "plot(y = f3(x_base), x = x_base, type = \"l\", lwd = 2, xlab = \"multiple minimizers\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf870a01",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "* Linear regression model MLE\n",
    "\n",
    "\n",
    "* Normal MLE. The (negative) log-likelihood \n",
    "\n",
    "$$\n",
    "\\ell (\\beta, \\sigma) = \\log \\sigma + \\frac{1}{2\\sigma^2}\\sum_{i=1}^n (y_i - x_i' \\beta)^2\n",
    "$$\n",
    "\n",
    "is not convex"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6da30f5",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Re-parameterize the criterion function by $\\gamma = 1/\\sigma$ and $\\alpha = \\beta / \\sigma$, then\n",
    "\n",
    "$$\n",
    "\\ell (\\alpha, \\gamma) = -\\log \\gamma + \\frac{1}{2}\n",
    "\\sum_{i=1}^n (\\gamma y_i - x_i' \\alpha)^2\n",
    "$$\n",
    "\n",
    "is convex in $\\alpha, \\gamma$\n",
    "\n",
    "* Many MLE estimators in econometric textbooks are convex. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "620629e0",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* [Gao and Shi (2021)](https://link.springer.com/article/10.1007/s10614-020-09995-z) explore the infrastructure in R for convex optimization with two econometric examples. \n",
    "\n",
    "* `CVXR` [(Fu, 2018)](https://arxiv.org/abs/1711.07582) is a convenient convex modeling language \n",
    "  * Proprietary solvers `CLEPX`, `MOSEK`, and `Gurubi`\n",
    "  * Open-source solvers `ECOS` and `SDPT3`. \n",
    "\n",
    "* `MOSEK` offers free academic license\n",
    "* [`Rmosek`](http://rmosek.r-forge.r-project.org/)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.1"
  },
  "rise": {
   "enable_chalkboard": true,
   "scroll": true,
   "theme": "serif"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}