Julia/notebooks/PortfolioOpt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization Example: A Portfolio Problem\n",
    "**Author:    David Childers**\n",
    "\n",
    "**Date:      3-29-2023**\n",
    "\n",
    "In this document I present several approaches to numerical optimization through the context of a portfolio choice problem with CARA utility. The approach will use in particular the Julia toolbox Optim.jl, along with a few related tools, for which this is intended to serve as an introduction. The problem is taken from Judd, *Computational Methods for Economics Ch. 4* and the code is based on my previous Matlab code for the same problem. Another useful introduction to these tools, to which this example is related, is available from QuantEcon at <https://julia.quantecon.org/more_julia/optimization_solver_packages.html>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Activating\u001b[22m\u001b[39m project at `~/Library/CloudStorage/Dropbox/Computational Methods/Code/Julia`\n"
     ]
    }
   ],
   "source": [
    "# #Set up the Julia environment to use standard packages using Project and Manifest in document's directory\n",
    "import Pkg; \n",
    "Pkg.activate(\".\")\n",
    "Pkg.instantiate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1mStatus\u001b[22m\u001b[39m `~/Library/CloudStorage/Dropbox/Computational Methods/Code/Julia/Project.toml`\n",
      " \u001b[90m [28f2ccd6] \u001b[39mApproxFun v0.13.17\n",
      " \u001b[90m [a134a8b2] \u001b[39mBlackBoxOptim v0.6.2\n",
      " \u001b[90m [e2554f3b] \u001b[39mClp v1.0.2\n",
      " \u001b[90m [0c46a032] \u001b[39mDifferentialEquations v7.7.0\n",
      " \u001b[90m [31c24e10] \u001b[39mDistributions v0.25.87\n",
      " \u001b[90m [442a2c76] \u001b[39mFastGaussQuadrature v0.5.1\n",
      " \u001b[90m [f6369f11] \u001b[39mForwardDiff v0.10.35\n",
      " \u001b[90m [7073ff75] \u001b[39mIJulia v1.24.0\n",
      " \u001b[90m [a98d9a8b] \u001b[39mInterpolations v0.14.7\n",
      " \u001b[90m [b6b21f68] \u001b[39mIpopt v1.2.1\n",
      " \u001b[90m [4076af6c] \u001b[39mJuMP v1.10.0\n",
      " \u001b[90m [0fc2ff8b] \u001b[39mLeastSquaresOptim v0.8.4\n",
      " \u001b[90m [d41bc354] \u001b[39mNLSolversBase v7.8.3\n",
      " \u001b[90m [2774e3e8] \u001b[39mNLsolve v4.5.1\n",
      " \u001b[90m [429524aa] \u001b[39mOptim v1.7.5\n",
      " \u001b[90m [d96e819e] \u001b[39mParameters v0.12.3\n",
      " \u001b[90m [91a5bcdd] \u001b[39mPlots v1.38.9\n",
      " \u001b[90m [8a4e6c94] \u001b[39mQuasiMonteCarlo v0.2.19\n",
      " \u001b[90m [f2b01f46] \u001b[39mRoots v2.0.10\n",
      " \u001b[90m [e88e6eb3] \u001b[39mZygote v0.6.60\n",
      " \u001b[90m [37e2e46d] \u001b[39mLinearAlgebra\n",
      " \u001b[90m [de0858da] \u001b[39mPrintf\n",
      " \u001b[90m [9a3f8284] \u001b[39mRandom\n",
      " \u001b[90m [10745b16] \u001b[39mStatistics\n"
     ]
    }
   ],
   "source": [
    "#Check status of packages\n",
    "Pkg.status()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Uncomment and execute following line if any of the following packages are missing\n",
    "#Pkg.add(LinearAlgebra, Statistics, ForwardDiff, Zygote, Optim, JuMP, Ipopt, BlackBoxOptim, Roots, NLsolve, LeastSquaresOptim)\n",
    "#Pkg.status()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load packages to be used for Optimization\n",
    "using LinearAlgebra, Statistics #Basic math\n",
    "#Specialized differentiation, optimization, and equation solving libraries\n",
    "using ForwardDiff, Zygote, Optim, JuMP, Ipopt, BlackBoxOptim, Roots, NLsolve, LeastSquaresOptim\n",
    "using Optim: converged, maximum, maximizer, minimizer, iterations #some functions from the Optim library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we first declare the function we want to optimize. The function represents a one-period portfolio optimization problem over a vector of $n$ assets, which each have a current price `assetprice` and pay out a certain defined payoff in each of $S$ total states, stored as $n\\times S$ array `assetreturns`. The investor is endowed with a vector of `endowments` of each asset, and receives payoff given by the utility of initial period consumption plus expected utility of returns, where expectation is taken with respect to a probability distribution over states with probability mass vector `stateprobs` and utility of returns is CARA with coefficient of absolute risk aversion `cara`. The choice the investor faces is to find the optimal vector of portfolio allocations to each asset, `pweights` by maximizing this utility(respectively, minimizing $-1$ times this utility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "portfolio_return (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Our minus expected utility function, that the investor would like to minimize\n",
    "\n",
    "function portfolio_return(pweights; stateprobs, endowment, assetprice, cara = 2)\n",
    "    stateutil = -exp.(-cara*(pweights'*assetreturns))\n",
    "    return mutility = exp(cara*assetprice'*(pweights-endowment))-sum(stateutil*stateprobs)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try some example parameter values to test the function out\n",
    "\n",
    "cara = 2 #This was set as a default, so you don't have to pass it to the function, but you can if you want\n",
    "stateprobs = [0.1, 0.2, 0.7] #Probability of each state\n",
    "\n",
    "endowment=[0.1, 0.1]    #Endowment of each asset\n",
    "assetprice=[1, 1]                #Price of each asset\n",
    "assetreturns=[2 1 0.1;         #Returns of each asset in each state\n",
    "             0 1 2] \n",
    "\n",
    "guessweight=[0.1,0.1]; #Initial guess of possible portfolio weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6610287876812315"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Evaluate directly\n",
    "stateutil = -exp.(-cara*(guessweight'*assetreturns))\n",
    "guessutility = exp(cara*assetprice'*(guessweight-endowment))-sum(stateutil*stateprobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6610287876812315"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mutilityguess = portfolio_return(guessweight, stateprobs=stateprobs, endowment=endowment, assetprice=assetprice, cara = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     Nelder-Mead\n",
       "\n",
       " * Convergence measures\n",
       "    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    28\n",
       "    f(x) calls:    59\n"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Declare an anonymous function, which takes only the portfolio choice as input\n",
    "portreturn = x->portfolio_return(x,stateprobs=stateprobs, endowment=endowment, assetprice=assetprice, cara = 2)\n",
    "# Use default optimization, with nothing more than an initial guess\n",
    "result = optimize(portreturn, guessweight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With no gradient, `optimize` runs the Nelder-Mead \"Amoeba\" method, which takes quite a few function calls to reach an optimum.  Unless you know your function is not differentiable, or derivatives are hard to compute, it may be strongly preferable to use a first or even second order method.  To try this, we need a Jacobian, and, for second order, a Hessian. In this problem, these have analytical forms, so I will provide them. I will then compare to results produced by automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "portGrad! (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Declare Gradient function: use an in-place, or \"non-allocating\" function\n",
    "# ! means it takes a preallocated vector and changes it, rather than creating a new vector\n",
    "#This saves memory, and can also improve speed due to reducing time in allocating new vectors\n",
    "# Note here we are defining this function without the additional arguments: they are taken from the memory environment\n",
    "# We would have to define again, or use an anonymous function, to change the constant choices\n",
    "function portGrad!(J,pweights)\n",
    "    stategrad = repeat(cara*exp.(-cara*(pweights'*assetreturns)),length(assetprice)).*assetreturns\n",
    "    return J .= cara*exp.(cara*assetprice'*(pweights-endowment))*assetprice-stategrad*stateprobs\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       "  1.3717574083973805\n",
       " -0.10785911389641445"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Evaluate at a particular parameter choice\n",
    "GradG = copy(guessweight)\n",
    "portGrad!(GradG,guessweight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     L-BFGS\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 4.62e-09 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 1.80e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n",
       "    |g(x)|                 = 1.37e-13 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    6\n",
       "    f(x) calls:    17\n",
       "    ∇f(x) calls:   17\n"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Call optimizer with gradient function\n",
    "result1 = optimize(portreturn,portGrad!,guessweight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With one derivative given, by default `Optim.jl` calls a first order method, here L-BFGS, a version of BFGS which does not store the entire Broyden approximation of the Jacobian. This is probably the most commonly used first order quasi-Newton method, with good tradeoff between number of iterations and iteration complexity. We see it achieved the same minimum in many fewer iterations and function calls than Nelder Mead, though it did require gradient calls.\n",
    "\n",
    "Analytic gradients can sometimes be a pain to code up, so it would be useful to also apply automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       "  1.3717574083973805\n",
       " -0.10785911389641445"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compare to result from ForwardDiff, which computes Gradient by forward mode automatic differentiation\n",
    "\n",
    "GG = ForwardDiff.gradient(portreturn,guessweight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pgrad! (generic function with 1 method)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To pass to optimizer, call as an anonymous function again, and also make non-allocating\n",
    "# ! means it takes a preallocated vector and changes it, rather than creating a new vector\n",
    "function pgrad!(G,x)\n",
    "    return G.=ForwardDiff.gradient(portreturn,x)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       "  1.3717574083973805\n",
       " -0.10785911389641445"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Evaluate Gradient\n",
    "pgrad!(GradG,guessweight) #Should be same as analytical formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     L-BFGS\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 4.62e-09 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 1.80e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n",
       "    |g(x)|                 = 1.37e-13 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    6\n",
       "    f(x) calls:    17\n",
       "    ∇f(x) calls:   17\n"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Call optimizer with Forward Diff gradient function\n",
    "result2 = optimize(portreturn,pgrad!,guessweight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That produced exactly the same result, at more or less similar speed, with no coding. You can call `optimize` without supplying a gradient function, but still choose a first order method like L-BFGS, and it will use a finite difference procedure to compute gradients. If the function is amenable to automatic differentiation, you can also call ForwardDiff as an option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     L-BFGS\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 4.62e-09 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 1.80e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 2.22e-16 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 1.57e-16 ≰ 0.0e+00\n",
       "    |g(x)|                 = 0.00e+00 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    6\n",
       "    f(x) calls:    18\n",
       "    ∇f(x) calls:   18\n"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Call optimizer with finite difference gradient function, by choosing option LBFGS()\n",
    "result3 = optimize(portreturn,guessweight,LBFGS())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     L-BFGS\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 4.62e-09 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 1.80e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n",
       "    |g(x)|                 = 1.37e-13 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    6\n",
       "    f(x) calls:    17\n",
       "    ∇f(x) calls:   17\n"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result4 = optimize(portreturn,guessweight,LBFGS(),autodiff=:forward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Autodiff saved us one function call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     Conjugate Gradient\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 1.96e-09 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 7.60e-09 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 2.22e-16 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 1.57e-16 ≰ 0.0e+00\n",
       "    |g(x)|                 = 3.33e-09 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    7\n",
       "    f(x) calls:    16\n",
       "    ∇f(x) calls:   9\n"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We can use first order methods other than L-BFGS. Try, e.g., ConjugateGradient()\n",
    "result5 = optimize(portreturn,guessweight,ConjugateGradient(),autodiff=:forward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different first order methods may have differing degrees of speed and reliability on different problems (e.g., conjugate gradient works well with sparse or low rank problems), but require similar user inputs. Second order methods, like Newton's method, may gain speed at the cost of requiring a Hessian. This can be user provided or automatic as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 4.0  0.0\n",
       " 0.0  0.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M=assetreturns[:,1]*assetreturns[:,1]' #A matrix the same size as the Hessian just as a placeholder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "portHess! (generic function with 1 method)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Calculate Analytical Hessian formula\n",
    "function portHess!(H,pweights)\n",
    "    #Calculate Hessian of expected utility by adding over states\n",
    "    nextperHess=zeros(size(H))\n",
    "    for i in 1:length(stateprobs)\n",
    "        nextperHess += stateprobs[i]*cara^2*exp(-cara*(pweights'*assetreturns[:,i]))*(assetreturns[:,i]*assetreturns[:,i]')\n",
    "    end\n",
    "    #Calculate Hessian of current period utility\n",
    "    hess1= cara^2*exp(cara*assetprice'*(pweights-endowment))*(assetprice*assetprice')\n",
    "    #Hessian of total utility of choice\n",
    "    return H .= nextperHess + hess1\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 5.62717   4.9042\n",
       " 4.9042   11.8952"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "portHess!(M,guessweight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     Newton's Method\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 2.57e-08 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 9.97e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 2.66e-15 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 1.89e-15 ≰ 0.0e+00\n",
       "    |g(x)|                 = 5.11e-15 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    4\n",
       "    f(x) calls:    12\n",
       "    ∇f(x) calls:   12\n",
       "    ∇²f(x) calls:  4\n"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Use this in optimize to apply a second order method\n",
    "result6 = optimize(portreturn,portGrad!,portHess!,guessweight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton's method converged in only 4 iterations, with 12 function and 12 gradent calls, but also required a Hessian. It was also notably slower than a quasi-Newton method, even with user-provided Hessian, which is a common occurence with Newton's method, since it requires Hessian evaluations and linear system solves each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 5.62717   4.9042\n",
       " 4.9042   11.8952"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#As with gradients, we can get Hessians from automatic differentiation\n",
    "ForwardDiff.hessian(portreturn,guessweight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " * Status: success\n",
       "\n",
       " * Candidate solution\n",
       "    Final objective value:     1.410094e+00\n",
       "\n",
       " * Found with\n",
       "    Algorithm:     Newton's Method\n",
       "\n",
       " * Convergence measures\n",
       "    |x - x'|               = 2.57e-08 ≰ 0.0e+00\n",
       "    |x - x'|/|x'|          = 9.97e-08 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|         = 2.66e-15 ≰ 0.0e+00\n",
       "    |f(x) - f(x')|/|f(x')| = 1.89e-15 ≰ 0.0e+00\n",
       "    |g(x)|                 = 5.11e-15 ≤ 1.0e-08\n",
       "\n",
       " * Work counters\n",
       "    Seconds run:   0  (vs limit Inf)\n",
       "    Iterations:    4\n",
       "    f(x) calls:    12\n",
       "    ∇f(x) calls:   12\n",
       "    ∇²f(x) calls:  4\n"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#We can use these in our optimizer by turning into an anonymous, in-place function\n",
    "function phess!(H,x)\n",
    "    return H.=ForwardDiff.hessian(portreturn,x)\n",
    "end\n",
    "\n",
    "#Call Newton's method\n",
    "result7 = optimize(portreturn,pgrad!,phess!,guessweight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structured Optimization and JuMP\n",
    "For generic optimization problems, algorithms which make few assumptions beyond the presence of local information can be reasonable. However, many problems have global structure which permits specialized optimization methods, which can perform better than generic methods, often overwhelmingly so. Problems with linear, quadratic, conic, semidefinite matrix, etc objectives, and linear, quadratic, integer, etc, constraints, have specialized methods available, often implemented in specialized software, including high quality commercial libraries like *Knitro* and *Gurobi*. For these programs, it is helpful to have a specialized modeling language which can express a problem in a suitable format for such a solver. **AMPL** is a commercial program which provides such a language, with interface to many solvers, available to try at <https://ampl.com/try-ampl/try-ampl-online/>. In Julia, the **JuMP** library has similar function, providing a modeling language and interface with external solvers.\n",
    "\n",
    "To see how it looks like, I will translate a linear programming example from the AMPL library: for more details on JuMP, see [the documentation](https://www.juliaopt.org/JuMP.jl/stable/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1:2"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set problem parameters\n",
    "\n",
    "avail = 40 #Total hours available\n",
    "#For each product, set tons produced per hour, profit per ton, and limit on tons sold in a week\n",
    "param = Dict(\"bands\"=>( rate = 200, profit = 25, market = 6000 ), \n",
    "             \"coils\"=>( rate = 140, profit = 30, market = 4000 ) )\n",
    "PROD = collect(keys(param)) #Products\n",
    "pindex = 1:length(PROD) #Set of products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{VariableRef}:\n",
       " make[1]\n",
       " make[2]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Build model in JuMP, with interior point optimizer Ipopt\n",
    "steel=Model(Ipopt.Optimizer) # with_optimizer deprecated\n",
    "\n",
    "#Declare variable to optimize, make, which is nonnegative \n",
    "@variable(steel,make[pindex]>=0,container=Array) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape},1,...} with index sets:\n",
       "    Dimension 1, 1:2\n",
       "And data, a 2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:\n",
       " marketconstraint[1] : make[1] ≤ 4000.0\n",
       " marketconstraint[2] : make[2] ≤ 6000.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Add market constraints\n",
    "@constraint(steel,marketconstraint[p=pindex],make[p]<=param[PROD[p]][:market])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "timeconstraint : $ 0.007142857142857143 make_{1} + 0.005 make_{2} \\leq 40.0 $"
      ],
      "text/plain": [
       "timeconstraint : 0.007142857142857143 make[1] + 0.005 make[2] ≤ 40.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Add time constraint\n",
    "@constraint(steel,timeconstraint,sum(make[p]/param[PROD[p]][:rate] for p in pindex)<=avail)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ 30 make_{1} + 25 make_{2} $$"
      ],
      "text/plain": [
       "30 make[1] + 25 make[2]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Set objective function: maximize profit\n",
    "@objective(steel,Max,sum(make[p]*param[PROD[p]][:profit] for p in pindex))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "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 https://github.com/coin-or/Ipopt\n",
      "******************************************************************************\n",
      "\n",
      "This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.\n",
      "\n",
      "Number of nonzeros in equality constraint Jacobian...:        0\n",
      "Number of nonzeros in inequality constraint Jacobian.:        4\n",
      "Number of nonzeros in Lagrangian Hessian.............:        0\n",
      "\n",
      "Total number of variables............................:        2\n",
      "                     variables with only lower bounds:        2\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...............:        3\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:        0\n",
      "        inequality constraints with only upper bounds:        3\n",
      "\n",
      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n",
      "   0  5.4999945e-01 0.00e+00 1.50e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0\n",
      "   1  2.1299934e+01 0.00e+00 3.00e+01  -1.0 4.00e-01    -  3.19e-02 1.00e+00f  1\n",
      "   2  3.8193949e+04 0.00e+00 2.97e+01  -1.0 1.23e+03    -  3.33e-04 1.00e+00f  1\n",
      "   3  1.2685907e+05 0.00e+00 2.87e+01  -1.0 8.55e+04    -  6.72e-03 3.21e-02f  1\n",
      "   4  1.7674002e+05 0.00e+00 2.83e+01  -1.0 1.40e+05    -  2.17e-03 1.42e-02f  1\n",
      "   5  1.7723897e+05 0.00e+00 1.55e+02  -1.0 4.24e+04    -  5.89e-02 4.68e-04f  1\n",
      "   6  1.7725419e+05 0.00e+00 1.15e+03  -1.0 8.60e+01    -  2.38e-01 3.20e-02f  1\n",
      "   7  1.9185238e+05 0.00e+00 2.08e+03  -1.0 5.04e+03    -  5.45e-03 7.24e-01f  1\n",
      "   8  1.9199843e+05 0.00e+00 5.40e+02  -1.0 1.16e+02    -  1.00e+00 3.13e-01f  1\n",
      "   9  1.9199980e+05 0.00e+00 1.00e-06  -1.0 3.44e-01    -  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  1.9200000e+05 0.00e+00 1.50e-09  -3.8 2.50e-02    -  1.00e+00 1.00e+00f  1\n",
      "  11  1.9200000e+05 0.00e+00 2.67e-14  -8.6 3.76e-05    -  1.00e+00 1.00e+00f  1\n",
      "\n",
      "Number of Iterations....: 11\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -1.9200000191999497e+05    1.9200000191999497e+05\n",
      "Dual infeasibility......:   2.6695555158452723e-14    2.6695555158452723e-14\n",
      "Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00\n",
      "Complementarity.........:   2.5082050994613992e-09    2.5082050994613992e-09\n",
      "Overall NLP error.......:   2.9831173875610469e-10    2.5082050994613992e-09\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 12\n",
      "Number of objective gradient evaluations             = 12\n",
      "Number of equality constraint evaluations            = 0\n",
      "Number of inequality constraint evaluations          = 12\n",
      "Number of equality constraint Jacobian evaluations   = 0\n",
      "Number of inequality constraint Jacobian evaluations = 1\n",
      "Number of Lagrangian Hessian evaluations             = 1\n",
      "Total seconds in IPOPT                               = 0.273\n",
      "\n",
      "EXIT: Optimal Solution Found.\n"
     ]
    }
   ],
   "source": [
    "#optimize using declared solver\n",
    "JuMP.optimize!(steel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "192000.00191999497"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Optimized objective: total profit\n",
    "objective_value(steel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Use 1400.0000140003542 units of coils in optimal production plan.\n",
      "Use 6000.000059999374 units of bands in optimal production plan.\n"
     ]
    }
   ],
   "source": [
    "#Optimal values of decisions\n",
    "for p in pindex\n",
    "    println(\"Use $(value(make[p])) units of $(PROD[p]) in optimal production plan.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#If Clp optimizer is missing, uncomment the following to install it.\n",
    "#Pkg.add(Clp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "using Clp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A JuMP Model\n",
       "Feasibility problem with:\n",
       "Variables: 0\n",
       "Model mode: AUTOMATIC\n",
       "CachingOptimizer state: EMPTY_OPTIMIZER\n",
       "Solver name: Clp"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Set optimizer as Clp, with Primal Simplex algorithm as solver: see https://github.com/JuliaOpt/Clp.jl\n",
    "steel2 = Model(optimizer_with_attributes(Clp.Optimizer, \"Algorithm\" => 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A JuMP Model\n",
       "Maximization problem with:\n",
       "Variables: 2\n",
       "Objective function type: AffExpr\n",
       "`AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 3 constraints\n",
       "`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 2 constraints\n",
       "Model mode: AUTOMATIC\n",
       "CachingOptimizer state: EMPTY_OPTIMIZER\n",
       "Solver name: Clp\n",
       "Names registered in the model: make, marketconstraint, timeconstraint"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Declare model again\n",
    "@variable(steel2,make[pindex]>=0,container=Array) \n",
    "@constraint(steel2,marketconstraint[p=pindex],make[p]<=param[PROD[p]][:market])\n",
    "@constraint(steel2,timeconstraint,sum(make[p]/param[PROD[p]][:rate] for p in pindex)<=avail)\n",
    "@objective(steel2,Max,sum(make[p]*param[PROD[p]][:profit] for p in pindex))\n",
    "steel2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "JuMP.optimize!(steel2) = nothing\n",
      "Coin0506I Presolve 1 (-2) rows, 2 (0) columns and 2 (-2) elements\n",
      "Clp0006I 0  Obj -0 Dual inf 65.714284 (2)\n",
      "Clp0006I 1  Obj 192000\n",
      "Clp0000I Optimal - objective value 192000\n",
      "Coin0511I After Postsolve, objective 192000, infeasibilities - dual 0 (0), primal 0 (0)\n",
      "Clp0032I Optimal objective 192000 - 1 iterations time 0.002, Presolve 0.00\n"
     ]
    }
   ],
   "source": [
    "#optimize using declared solver\n",
    "@show JuMP.optimize!(steel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "192000.0"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective_value(steel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Use 1400.000000000001 units of coils in optimal production plan.\n",
      "Use 5999.999999999999 units of bands in optimal production plan.\n"
     ]
    }
   ],
   "source": [
    "#Optimal values of decisions\n",
    "for p in pindex\n",
    "    println(\"Use $(value(make[p])) units of $(PROD[p]) in optimal production plan.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Are these integers, to machine precision?\n",
    "value(make[1])≈1400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value(make[2])≈6000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Switching from an interior point solver, as in Ipopt, to a simplex method solver resulted in exact solutions, up to machine precision. The reason is that the simplex method moves along constraints, while an interior point solver, as the name suggests, gradually relaxes a penalty that keeps it always strictly inside the constraints. A simplex method is fast and exact \"most of the time\" in a formal sense, but can behave poorly on worst case problems. Interior point methods for linear programs instead give good approximation guarantees in all cases."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}