Nonlinear optimization with JuMP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nonlinear optimization using JuMP (Julia for Mathematical Programming)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By Tyler Ransom, Duke University Social Science Research Institute"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's examine how to use a new `Julia` packaged called `JuMP` to illustrate some of Julia's powerful resources for optimizing objective functions of all kinds.\n",
    "\n",
    "For an overview of the `JuMP` package, see [here](http://jump.readthedocs.org/en/latest/). Essentially, `JuMP` makes use of `Julia`'s macros to greatly simplify construction and optimization of various problems.\n",
    "\n",
    "Here are some features of `JuMP`:\n",
    "- interfaces seamlessly with many industry-grade solvers (e.g. AMPL)\n",
    "- can be used to solve linear programming, nonlinear programming, and\n",
    "many other types of problems (including constrained optimization)\n",
    "- automatically differentiates the objective function (*not* numerical\n",
    "differentiation), resulting in speed gains\n",
    "- user-friendly model construction: user simply writes the objective\n",
    "function and any constraints in mathematical notation; JuMP then\n",
    "translates this into binaries at compile time\n",
    "\n",
    "In this illustration, we will use JuMP for nonlinear programming problems (i.e. maximum likelihood estimation). We will be using the open-source optimizer `Ipopt` (pronounced eye-PEE-opt)\n",
    "\n",
    "So let's load the packages we'll need and get started!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[34mINFO: Nothing to be done\n",
      "\u001b[0m\u001b[1m\u001b[34mINFO: METADATA is out-of-date — you may not have the latest version of JuMP\n",
      "\u001b[0m\u001b[1m\u001b[34mINFO: Use `Pkg.update()` to get the latest versions of your packages\n",
      "\u001b[0m\u001b[1m\u001b[34mINFO: Nothing to be done\n",
      "\u001b[0m\u001b[1m\u001b[34mINFO: METADATA is out-of-date — you may not have the latest version of Ipopt\n",
      "\u001b[0m\u001b[1m\u001b[34mINFO: Use `Pkg.update()` to get the latest versions of your packages\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "Pkg.add(\"JuMP\")\n",
    "Pkg.add(\"Ipopt\")\n",
    "using JuMP\n",
    "using Ipopt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now let's create some data inside of a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "datagen (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function datagen(N,T)\n",
    "    ## Generate data for a linear model to test optimization\n",
    "    srand(1234)\n",
    "    \n",
    "    N = convert(Int64,N) #inputs to functions such as -ones- need to be integers!\n",
    "    T = convert(Int64,T) #inputs to functions such as -ones- need to be integers!\n",
    "    const n = convert(Int64,N*T) # use -const- as a way to declare a variable to be global (so other functions can access it)\n",
    "    \n",
    "    # generate the Xs\n",
    "    const X = cat(2,ones(N*T,1),5+3*randn(N*T,1),rand(N*T,1),\n",
    "               2.5+2*randn(N*T,1),15+3*randn(N*T,1),\n",
    "               .7-.1*randn(N*T,1),5+3*randn(N*T,1),\n",
    "               rand(N*T,1),2.5+2*randn(N*T,1),\n",
    "               15+3*randn(N*T,1),.7-.1*randn(N*T,1),\n",
    "               5+3*randn(N*T,1),rand(N*T,1),2.5+2*randn(N*T,1),\n",
    "               15+3*randn(N*T,1),.7-.1*randn(N*T,1));\n",
    "    \n",
    "    # generate the betas (X coefficients)\n",
    "    const ßans = [ 2.15; 0.10;  0.50; 0.10; 0.75; 1.2; 0.10;  0.50; 0.10; 0.75; 1.2; 0.10;  0.50; 0.10; 0.75; 1.2 ]\n",
    "    \n",
    "    # generate the std dev of the errors\n",
    "    const σans = 0.3\n",
    "    \n",
    "    # generate the Ys from the Xs, betas, and error draws\n",
    "    draw = 0 + σans*randn(N*T,1)\n",
    "    const y = X*ßans+draw\n",
    "    \n",
    "    # return generated data so that other functions (below) have access\n",
    "    return X,y,ßans,σans,n\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now let's evaluate the function so that the function outputs are in the current scope:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X,y,ßans,σans,n = datagen(1e4,5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## OLS estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can estimate the `bAns` parameter vector using OLS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ßhat = (X'*X)\\(X'*y);\n",
    "σhat = sqrt((y-X*ßhat)'*(y-X*ßhat)/(n-size(X,2)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can compare the estimates with the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17×2 Array{Float64,2}:\n",
       " 2.101  2.15\n",
       " 0.1    0.1 \n",
       " 0.497  0.5 \n",
       " 0.1    0.1 \n",
       " 0.751  0.75\n",
       " 1.215  1.2 \n",
       " 0.1    0.1 \n",
       " 0.505  0.5 \n",
       " 0.1    0.1 \n",
       " 0.75   0.75\n",
       " 1.208  1.2 \n",
       " 0.1    0.1 \n",
       " 0.509  0.5 \n",
       " 0.101  0.1 \n",
       " 0.75   0.75\n",
       " 1.221  1.2 \n",
       " 0.3    0.3 "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[round([ßhat;σhat],3) [ßans;σans]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimator performs well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Estimating the parameters of a normal MLE using `JuMP`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now let's estimate this using maximum likelihood under the (valid) assumption that the errors in our `datagen()` function are independently drawn from a normal distribution with mean 0 and standard deviation to be estimated.\n",
    "\n",
    "Recall that the likelihood function for this scenario is:\n",
    "$$ L_{i} = \\prod_{i=1}^{N} \\frac{1}{\\sqrt{2\\pi\\sigma^{2}}} \\exp \\left( -\\frac{\\left(y_{i}-x_{i}\\beta\\right)^{2}}{2\\sigma^{2}}\\right) $$\n",
    "\n",
    "or, more simply,\n",
    "\n",
    "$$ L_{i} = \\prod_{i=1}^{N} \\frac{1}{\\sigma} \\phi \\left(\\frac{y_{i}-x_{i}\\beta}{\\sigma}\\right) $$\n",
    "\n",
    "where $\\phi\\left(\\cdot\\right)$ is the probability density function of the standard normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In maximum likelihood estimation, we search for the parameters that maximize the `log` likelihood function. Our function above becomes:\n",
    "\n",
    "$$ l_{i} = \\sum_{i=1}^{N} -\\frac{1}{2}\\log(2\\pi) - \\log(\\sigma) -\\frac{\\left(y_{i}-x_{i}\\beta\\right)^{2}}{2\\sigma^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### `JuMP` syntax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`JuMP` requires the following syntax elements:\n",
    "- model name\n",
    "- solver\n",
    "- variable (i.e. a parameter to search over)\n",
    "- objective function\n",
    "- `solve()` call\n",
    "\n",
    "We will construct these line-by-line below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "First, let's decleare the model name and solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myMLE = Model(solver=IpoptSolver(tol=1e-6));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now let's tell `JuMP` which variables to search over. For the example above, this is the vector $\\beta$ and the scalar $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "@variable(myMLE, ß[i=1:size(X,2)]);\n",
    "@variable(myMLE, σ>=0.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now we will write the objective function, which is the log likelihood function from above. We will also tell `JuMP` that we are maximizing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "@NLobjective(myMLE, Max, (n/2)*log(1/(2*π*σ^2))-sum((y[i]-sum(X[i,k]*ß[k] \n",
    "                        for k=1:size(X,2)))^2 for i=1:size(X,1))/(2σ^2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we specify the summation term over $i$ using `sum(... for i=1:size(X,1))`. We also have to specify the linear algebra implied by `X*b` as another `sum`, this time over $k$. This is because `JuMP` doesn't understand matrix multiplication. This last fact adds a little bit of syntax to the objective function, but overall it is still very readable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now that we have all of our elements set up, we can issue a `solve()` call and ask `JuMP` to perform the optimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "******************************************************************************\n",
      "\n",
      "This is Ipopt version 3.12.1, running with linear solver mumps.\n",
      "NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:        0\n",
      "Number of nonzeros in Lagrangian Hessian.............:      153\n",
      "\n",
      "Total number of variables............................:       17\n",
      "                     variables with only lower bounds:        1\n",
      "                variables with lower and upper bounds:        0\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:        0\n",
      "Total number of inequality constraints...............:        0\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        0\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  4.3317673e+11 0.00e+00 8.66e+05  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  2.0586430e+11 0.00e+00 3.34e+05  -1.0 1.21e-01   4.0 1.00e+00 1.00e+00f  1\n",
      "   2  7.2604941e+10 0.00e+00 1.07e+05  -1.0 2.36e-01   3.5 1.00e+00 1.00e+00f  1\n",
      "   3  4.9336093e+09 0.00e+00 1.21e+04  -1.0 7.63e-01   3.0 1.00e+00 5.00e-01f  2\n",
      "   4  1.6112915e+09 0.00e+00 6.71e+03  -1.0 8.34e-02   3.5 1.00e+00 1.00e+00f  1\n",
      "   5  7.4680240e+08 0.00e+00 2.71e+03  -1.0 1.75e-02   3.0 1.00e+00 1.00e+00f  1\n",
      "   6  3.8887145e+08 0.00e+00 1.11e+03  -1.0 1.60e-02   2.5 1.00e+00 1.00e+00f  1\n",
      "   7  2.0253796e+08 0.00e+00 4.55e+02  -1.0 1.93e-02   2.0 1.00e+00 1.00e+00f  1\n",
      "   8  1.0604602e+08 0.00e+00 1.88e+02  -1.0 1.98e-02   1.6 1.00e+00 1.00e+00f  1\n",
      "   9  5.6063504e+07 0.00e+00 7.90e+01  -1.0 1.66e-02   1.1 1.00e+00 1.00e+00f  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  10  2.9618543e+07 0.00e+00 3.35e+01  -1.0 1.22e-02   0.6 1.00e+00 1.00e+00f  1\n",
      "  11  1.5335879e+07 0.00e+00 1.44e+01  -1.0 1.56e-02   0.1 1.00e+00 1.00e+00f  1\n",
      "  12  9.7927449e+06 0.00e+00 3.32e+01  -1.0 2.65e+01    -  1.00e+00 4.06e-02f  3\n",
      "  13  4.0837562e+06 0.00e+00 4.70e+00  -1.0 6.50e-01    -  1.00e+00 1.00e+00f  1\n",
      "  14  1.4963984e+06 0.00e+00 4.15e+00  -1.0 9.52e-02    -  1.00e+00 1.00e+00f  1\n",
      "  15  5.8425937e+05 0.00e+00 7.66e-01  -1.0 5.80e-02    -  1.00e+00 1.00e+00f  1\n",
      "  16  2.4222369e+05 0.00e+00 1.41e-01  -1.7 6.18e-02    -  1.00e+00 1.00e+00f  1\n",
      "  17  1.0772293e+05 0.00e+00 3.12e-02  -2.5 1.37e-01    -  1.00e+00 1.00e+00f  1\n",
      "  18  3.2318192e+04 0.00e+00 2.49e-02  -2.5 6.53e-02    -  1.00e+00 1.00e+00f  1\n",
      "  19  1.9786699e+04 0.00e+00 2.66e-03  -3.8 1.96e-01    -  9.23e-01 1.00e+00f  1\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "  20  1.1765538e+04 0.00e+00 1.32e-03  -3.8 5.44e-02    -  1.00e+00 1.00e+00f  1\n",
      "  21  1.0821210e+04 0.00e+00 1.65e-04  -5.7 2.77e-02    -  9.22e-01 1.00e+00f  1\n",
      "  22  1.0715195e+04 0.00e+00 1.68e-05  -5.7 1.23e-02    -  1.00e+00 1.00e+00f  1\n",
      "  23  1.0714713e+04 0.00e+00 9.00e-08  -8.6 9.13e-04    -  9.98e-01 1.00e+00f  1\n",
      "  24  1.0714713e+04 0.00e+00 7.30e-12  -8.6 8.84e-06    -  1.00e+00 1.00e+00f  1\n",
      "  25  1.0714713e+04 0.00e+00 5.23e-14 -12.9 7.51e-07    -  1.00e+00 1.00e+00f  1\n",
      "\n",
      "Number of Iterations....: 25\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:   1.0714712841850888e-04    1.0714712841850887e+04\n",
      "Dual infeasibility......:   5.2259077749931546e-14    5.2259077749931544e-06\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   1.3171685520953427e-13    1.3171685520953427e-05\n",
      "Overall NLP error.......:   1.3171685520953427e-13    1.3171685520953427e-05\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 37\n",
      "Number of objective gradient evaluations             = 26\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 0\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 0\n",
      "Number of Lagrangian Hessian evaluations             = 25\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      0.332\n",
      "Total CPU secs in NLP function evaluations           =     14.136\n",
      "\n",
      "EXIT: Optimal Solution Found.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       ":Optimal"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "status = solve(myMLE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We see that the solver finished at a solution that is a local maximum. To store the estimated values of the parameters, we can use the `getValue()` functions, applied to each separate parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ßopt = getvalue(ß[:]);\n",
    "σopt = getvalue(σ);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And if we compare them with the solution, we see that we have indeed obtained the correct parameter estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17×2 Array{Float64,2}:\n",
       " 2.101  2.15\n",
       " 0.1    0.1 \n",
       " 0.497  0.5 \n",
       " 0.1    0.1 \n",
       " 0.751  0.75\n",
       " 1.215  1.2 \n",
       " 0.1    0.1 \n",
       " 0.505  0.5 \n",
       " 0.1    0.1 \n",
       " 0.75   0.75\n",
       " 1.208  1.2 \n",
       " 0.1    0.1 \n",
       " 0.509  0.5 \n",
       " 0.101  0.1 \n",
       " 0.75   0.75\n",
       " 1.221  1.2 \n",
       " 0.3    0.3 "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[round([ßopt;σopt],3) [ßans;σans]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also return other helpful values such as the objective function value at the optimum, which is the log likelihood of the estimated model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-10714.712841850887"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getobjectivevalue(myMLE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Additionally, we can obtain the Hessian matrix of the model, which serves as a key component to the variance-covariance matrix that is used for statistical inference. The syntax for this is a bit messy, but I will put it below for reference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "this_par = myMLE.colVal\n",
    "m_eval = JuMP.NLPEvaluator(myMLE);\n",
    "MathProgBase.initialize(m_eval, [:ExprGraph, :Grad, :Hess])\n",
    "hess_struct = MathProgBase.hesslag_structure(m_eval)\n",
    "hess_vec = zeros(length(hess_struct[1]))\n",
    "numconstr = length(m_eval.m.linconstr) + length(m_eval.m.quadconstr) + length(m_eval.m.nlpdata.nlconstr)\n",
    "dimension = length(myMLE.colVal)\n",
    "MathProgBase.eval_hesslag(m_eval, hess_vec, this_par, 1.0, zeros(numconstr))\n",
    "this_hess_ld = sparse(hess_struct[1], hess_struct[2], hess_vec, dimension, dimension)\n",
    "hOpt = this_hess_ld + this_hess_ld' - sparse(diagm(diag(this_hess_ld)));\n",
    "hOpt = -full(hOpt); #since we are maximizing\n",
    "seOpt = sqrt(diag(full(hOpt)\\eye(size(hOpt,1))));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The above code obtained the estimated hessian matrix, took the negative of it (because we are maximizing), and then obtained standard errors of our parameter estimates from the diagonal elements of the inverse of the negative hessian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17×2 Array{Float64,2}:\n",
       " 2.101  0.021\n",
       " 0.1    0.0  \n",
       " 0.497  0.005\n",
       " 0.1    0.001\n",
       " 0.751  0.0  \n",
       " 1.215  0.013\n",
       " 0.1    0.0  \n",
       " 0.505  0.005\n",
       " 0.1    0.001\n",
       " 0.75   0.0  \n",
       " 1.208  0.013\n",
       " 0.1    0.0  \n",
       " 0.509  0.005\n",
       " 0.101  0.001\n",
       " 0.75   0.0  \n",
       " 1.221  0.013\n",
       " 0.3    0.001"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[round([ßopt;σopt],3) round(seOpt,3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can compare these with the OLS standard errors and confirm that they are the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16×2 Array{Float64,2}:\n",
       " 2.101  0.021\n",
       " 0.1    0.0  \n",
       " 0.497  0.005\n",
       " 0.1    0.001\n",
       " 0.751  0.0  \n",
       " 1.215  0.013\n",
       " 0.1    0.0  \n",
       " 0.505  0.005\n",
       " 0.1    0.001\n",
       " 0.75   0.0  \n",
       " 1.208  0.013\n",
       " 0.1    0.0  \n",
       " 0.509  0.005\n",
       " 0.101  0.001\n",
       " 0.75   0.0  \n",
       " 1.221  0.013"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seHat = sqrt((σhat^2).*diag((X'*X)\\eye(size(X,2))));\n",
    "[round(ßhat,3) round(seHat,3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Constrained optimization in JuMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constraints can easily be added to the objective function. There are two ways of doing this:\n",
    "1. Directly to the `@variable` statement, e.g. `@variable(myMLE, σ>=0.0);` restricts `σ` to be weakly positive.\n",
    "2. Additional lines after the set of `@variable` statements, but before the `@NLobjective` statement. The syntax for this is, e.g. `@constraint(myMLE, ß[15] == 0)` which would restrict the 15th element of the `ß` vector to equal 0.\n",
    "\n",
    "`JuMP` accepts up to literally thousands of constraints, and its simple syntax makes coding much easier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Final Word"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Pros of `JuMP`:\n",
    "- Versatile modeling language that takes advantage of `Julia`'s features.\n",
    "- Interfaces to a variety of optimizers, meaning the user can simply call a different optimizer without adjusting any underlying code.\n",
    "- Very fast. Roughly 3+ times faster than `Matlab`'s `fminunc` for a simple large-scale problem.\n",
    "- User can specify maximization or minimization on the fly without needing to re-express the objective function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Cons of `JuMP`:\n",
    "- Because JuMP does constrained optimization if the user specifies constraints on any of the parameters. Thus, the hessian it returns is the hessian of the corresponding *Lagrangian*, not the hessian of the *objective*. This makes statistical inference of the parameter estimates a bit trickier.\n",
    "- No matrix multiplication is allowed, so users must specify matrix operations in scalar form"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}