ipynb/arnav_sood-lcp_simple.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Optimal Stopping Problems\n",
                "\n",
                "> Presented by [Arnav Sood](email:arnav.sood@ubc.ca) and [Jesse Perla](email:jesseperla@gmail.com)\n",
                "\n",
                "This notebook provides an expansion of lecture notes and Julia code to solve optimal stopping\n",
                "problems numerically as a [Differential Variational Inequality](https://en.wikipedia.org/wiki/Differential_variational_inequality) (DVI).  For more, see\n",
                "\n",
                "* [Huang, J. and J.-S. Pang (1998): \"Option Pricing and Linear Complementarity,\"](Journal of Computational Finance)\n",
                "* [Ben Moll's excellent notes](http://www.princeton.edu/~moll/HACTproject/option_simple.pdf)\n",
                "\n",
                "## A Univariate Example\n",
                "\n",
                "Assume the following\n",
                "\n",
                "* A stochastic process $x$, with [infinitesimal generator](https://en.wikipedia.org/wiki/Infinitesimal_generator_(stochastic_processes)) $\\mathcal{L}_x$\n",
                "  and boundary conditions.  (e.g. the infinitesimal generator for Brownian motion with drift $\\mu$ and variance $\\sigma$ is\n",
                "  $\\mathcal{L}_x \\equiv \\mu \\partial_x + \\frac{\\sigma^2}{2}\\partial_{xx}$.)\n",
                "* Boundary conditions for the $x$ stochastic process (e.g., reflecting barriers at some minimum and maximum.)\n",
                "* Flow payoffs $u(x)$ discounted at rate $\\rho > 0$.\n",
                "* Stopping payoff $S(x)$ which gives the value of stopping at state $x$.\n",
                "\n",
                "The control problem is: choose when to stop the process and gain the payoff $S(x)$.  If there are only two regions, then this becomes the choice of a stopping threshold $\\hat{x}$.\n",
                "\n",
                "#### Classical Formulation\n",
                "\n",
                "The \"classical\" formulation (e.g., [Stokey's \"The Economics of Inaction\"](https://press.princeton.edu/titles/8766.html)) of this problem is as a free-boundary problem.\n",
                "\n",
                "First, in a stationary environment, the Bellman equation within the interior (continuation) region is\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  \\rho v(x) = u(x) + \\mathcal{L}_x v(x) \\tag{1}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "Define a new operator and rewrite,\n",
                "\n",
                "$$\n",
                "\\begin{align}\n",
                "  \\mathcal{L} v(x) &= u(x)\\label{eq:bellman}\\tag{2} \\\\\n",
                "  \\mathcal{L} &\\equiv \\rho - \\mathcal{L}_x\\label{eq:L-op} \\tag{3}\n",
                "\\end{align}\n",
                "$$\n",
                "\n",
                "Along with boundary conditions from the stochastic process (e.g., a reflecting barrier at an $x_{\\min}$ and another at some $x_{\\max}$, or a transversality condition.)\n",
                "\n",
                "In special cases with sufficient monotonicity, there is a single stopping region, defined by a threshold at $\\hat{x}$.\n",
                "\n",
                "With this, the necessary conditions for optimality are the value matching condition\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  v(\\hat{x}) = S(\\hat{x})\\label{eq:vm} \\tag{4}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "and the smooth pasting condition\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  \\partial_x v(\\hat{x}) = 0\\label{eq:sp} \\tag{5}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "These conditions are necessary for optimality for any free-boundary problem.  Typically, you would solve $\\eqref{eq:bellman}, \\eqref{eq:vm},\\eqref{eq:sp},$ to find a value function $v(x)$ and a stopping threshold $\\hat{x}$.\n",
                "\n",
                "While this approach can be analytically tractable, it is inconvenient for numerical solutions if: (1) the equations cannot be solved in closed form; (2) there are multiple stopping regions; or (3) the dimension increases, leading to complicated stopping regions separated by non-linear manifolds.\n",
                "\n",
                "#### HJB Variational Inequality\n",
                "\n",
                "As discussed in [Ben Moll's notes](http://www.princeton.edu/~moll/HACTproject/option_simple.pdf), the [DVI](https://en.wikipedia.org/wiki/Differential_variational_inequality) form of the above problem is to solve\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  0 = \\min\\{\\mathcal{L} v(x) - u(x), v(x) - S(x)\\}\\label{eq:DVI}\\tag{6}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "for $v(x)$.  A few important notes:\n",
                "\n",
                "* The $\\mathcal{L} v(x) - u(x)$ term is just the Bellman equation $\\eqref{eq:bellman}$ if it is binding and $=0$.\n",
                "* The $v(x) - S(x)$ term is the value matching condition in $\\eqref{eq:vm}$ if it is binding and $=0$.\n",
                "* You are **not directly solving for the boundaries** of the stopping regions (i.e., there is no $\\hat{x}$ in the system.)\n",
                "\n",
                "Given a solution for $v(x)$, $\\eqref{eq:DVI}$ determines which regime a particular $x$ is in.  Those boundaries can be used to determine where an $\\hat{x}$ occurs.\n",
                "\n",
                "Furthermore, you can prove that at any of those points, a solution to $\\eqref{eq:DVI}$ implies the smooth-pasting condition $\\eqref{eq:sp}$.\n",
                "\n",
                "## Numerical Solution\n",
                "\n",
                "To solve the equation numerically, we first discretize the state-space, value function, and operators with  upwind finite difference methods.\n",
                "\n",
                "Abusing notation, define\n",
                "* The grid of $M$ values in the domain, $x \\equiv \\{x_1, x_2, \\ldots x_M\\}$.\n",
                "* Vectors for the values of functions at grid points: $v \\equiv \\{v(x_m)\\}_{m=1}^M$ and $S \\equiv \\{S(x_m)\\}_{m=1}^M$.\n",
                "* The matrix $L \\in \\mathbb{R}^{M \\times M}$ as the discretization of the $\\mathcal{L}$ in $\\eqref{eq:L-op}$, such that any boundary conditions on the stochastic process are fulfilled (i.e. properties of the $x$ such as reflecting barriers.)\n",
                "\n",
                "Note: See the [document](https://quantecon.github.io/SimpleDifferentialOperators.jl/dev/generated/discretized-differential-operator-derivation.pdf) in this repository for notes on the discretization of the $\\mathcal{L}$ to form the $L$ operator - subject to various boundary conditions.\n",
                "\n",
                "Given these definitions, the DVI in $\\eqref{eq:DVI}$ becomes the variational inequality\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  0 = \\min\\{L v - u, v - S\\}\\label{eq:VI}\\tag{7}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "To solve $\\eqref{eq:VI}$ for the $v\\in\\mathbb{R}^M$ vector, we can express it as a [linear complementarity problem](https://en.wikipedia.org/wiki/Linear_complementarity_problem) (or LCP).\n",
                "\n",
                "First, do a change of variables\n",
                "\n",
                "$$\n",
                "\\begin{align}\n",
                "  z &\\equiv v - S\\tag{8} \\\\ \n",
                "  q &\\equiv -u + L S\\tag{9} \\\\ \n",
                "  w &\\equiv L z + q\\tag{10}\n",
                "\\end{align}\n",
                "$$\n",
                "\n",
                "Substituting into $\\eqref{eq:VI}$ yields\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  0 = \\min \\{L z + q, z\\}\\tag{11}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "Or, in the typical notation of a complementarity problem\n",
                "\n",
                "$$\n",
                "\\begin{equation}\n",
                "  0 \\leq (L z + q) \\perp z \\geq 0\\tag{12}\n",
                "\\end{equation}\n",
                "$$\n",
                "\n",
                "## Implementation\n",
                "\n",
                "We will solve an instance of this model where\n",
                "\n",
                "$$\n",
                "\\begin{align}\n",
                "  u(x) &\\equiv x^{\\gamma}\\tag{13}\\\\\n",
                "  \\mathcal{L}_x &\\equiv \\mu \\partial_x + \\frac{\\sigma^2}{2}\\partial_{xx}\\tag{14}\\\\\n",
                "  S(x) &\\equiv S_0\\tag{15}\n",
                "\\end{align}\n",
                "$$\n",
                "\n",
                "First, we install and load necessary packages.\n",
                "\n",
                "The key one is the [`PATHSolver.jl`](https://github.com/chkwon/PATHSolver.jl) package for complementarity problems.\n",
                "\n",
                "This is a commercial solver due to Dirkse, Ferris, and Munson, who have generously made it [available](http://pages.cs.wisc.edu/~ferris/path.html) for non-commercial academic use."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "metadata": {},
            "outputs": [],
            "source": [
                "# install if required,\n",
                "# ] add SimpleDifferentialOperators PATHSolver Plots Suppressor\n",
                "\n",
                "using SimpleDifferentialOperators, PATHSolver, Plots, LinearAlgebra, Suppressor"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Next, define the necessary parameters and equations"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {},
            "outputs": [],
            "source": [
                "μ = -0.1\n",
                "σ = 0.1\n",
                "ρ = 0.05\n",
                "γ = 0.5\n",
                "ρ = 0.05\n",
                "S₀ = 20.0\n",
                "u(x) = x^γ  # (13)\n",
                "S(x) = S₀  # (15)\n",
                "@assert μ <= 0 # otherwise we need to swap the L₁₋ operator"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Next, discretize the grid, operators, and change of variables"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "metadata": {},
            "outputs": [],
            "source": [
                "M = 300\n",
                "x = range(0.0, 5.0, length = M)\n",
                "bc = (Reflecting(), Reflecting())  # assume x_t reflects at 0.0 and 5.0\n",
                "L = I * ρ - μ*L₁₋(x, bc) - σ^2 / 2 * L₂(x, bc)  # (2) and (14)\n",
                "\n",
                "q = -u.(x) + L*S.(x)  # (9)\n",
                "w(z) = L*z + q; # (10)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Finally, solve the model by passing to the PATH solver"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/svg+xml": [
                            "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
                            "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
                            "<defs>\n",
                            "  <clipPath id=\"clip9300\">\n",
                            "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
                            "  </clipPath>\n",
                            "</defs>\n",
                            "<defs>\n",
                            "  <clipPath id=\"clip9301\">\n",
                            "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
                            "  </clipPath>\n",
                            "</defs>\n",
                            "<polygon clip-path=\"url(#clip9301)\" points=\"\n",
                            "0,1600 2400,1600 2400,0 0,0 \n",
                            "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
                            "<defs>\n",
                            "  <clipPath id=\"clip9302\">\n",
                            "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
                            "  </clipPath>\n",
                            "</defs>\n",
                            "<polygon clip-path=\"url(#clip9301)\" points=\"\n",
                            "211.005,1440.48 2321.26,1440.48 2321.26,125.984 211.005,125.984 \n",
                            "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
                            "<defs>\n",
                            "  <clipPath id=\"clip9303\">\n",
                            "    <rect x=\"211\" y=\"125\" width=\"2111\" height=\"1315\"/>\n",
                            "  </clipPath>\n",
                            "</defs>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  270.729,1440.48 270.729,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  668.89,1440.48 668.89,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  1067.05,1440.48 1067.05,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  1465.21,1440.48 1465.21,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  1863.37,1440.48 1863.37,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  2261.54,1440.48 2261.54,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  211.005,1320.88 2321.26,1320.88 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  211.005,1073.68 2321.26,1073.68 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  211.005,826.473 2321.26,826.473 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  211.005,579.269 2321.26,579.269 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
                            "  211.005,332.066 2321.26,332.066 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,1440.48 2321.26,1440.48 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,1440.48 211.005,125.984 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  270.729,1440.48 270.729,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  668.89,1440.48 668.89,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  1067.05,1440.48 1067.05,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  1465.21,1440.48 1465.21,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  1863.37,1440.48 1863.37,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  2261.54,1440.48 2261.54,1420.77 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,1320.88 242.659,1320.88 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,1073.68 242.659,1073.68 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,826.473 242.659,826.473 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,579.269 242.659,579.269 \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  211.005,332.066 242.659,332.066 \n",
                            "  \"/>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 270.729, 1494.48)\" x=\"270.729\" y=\"1494.48\">0</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 668.89, 1494.48)\" x=\"668.89\" y=\"1494.48\">1</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1067.05, 1494.48)\" x=\"1067.05\" y=\"1494.48\">2</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1465.21, 1494.48)\" x=\"1465.21\" y=\"1494.48\">3</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1863.37, 1494.48)\" x=\"1863.37\" y=\"1494.48\">4</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2261.54, 1494.48)\" x=\"2261.54\" y=\"1494.48\">5</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 187.005, 1338.38)\" x=\"187.005\" y=\"1338.38\">21</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 187.005, 1091.18)\" x=\"187.005\" y=\"1091.18\">24</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 187.005, 843.973)\" x=\"187.005\" y=\"843.973\">27</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 187.005, 596.769)\" x=\"187.005\" y=\"596.769\">30</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 187.005, 349.566)\" x=\"187.005\" y=\"349.566\">33</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1266.13, 73.2)\" x=\"1266.13\" y=\"73.2\">Value Function for Optimal Stopping</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1266.13, 1590.4)\" x=\"1266.13\" y=\"1590.4\">x</text>\n",
                            "</g>\n",
                            "<g clip-path=\"url(#clip9301)\">\n",
                            "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 57.6, 783.233)\" x=\"57.6\" y=\"783.233\">v(x)</text>\n",
                            "</g>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  270.729,1403.28 277.387,1403.28 284.046,1403.28 290.704,1403.28 297.362,1403.28 304.02,1403.28 310.678,1403.28 317.337,1403.28 323.995,1403.28 330.653,1403.28 \n",
                            "  337.311,1403.28 343.97,1403.28 350.628,1403.28 357.286,1403.28 363.944,1403.28 370.602,1403.28 377.261,1403.28 383.919,1403.28 390.577,1403.28 397.235,1403.28 \n",
                            "  403.894,1403.28 410.552,1403.28 417.21,1403.28 423.868,1403.28 430.526,1403.28 437.185,1403.28 443.843,1403.28 450.501,1403.28 457.159,1403.28 463.817,1403.28 \n",
                            "  470.476,1403.28 477.134,1403.28 483.792,1403.28 490.45,1403.28 497.109,1403.28 503.767,1403.28 510.425,1403.28 517.083,1403.28 523.741,1403.28 530.4,1403.28 \n",
                            "  537.058,1403.28 543.716,1403.28 550.374,1403.28 557.032,1403.28 563.691,1403.28 570.349,1403.28 577.007,1403.28 583.665,1403.28 590.324,1403.28 596.982,1403.28 \n",
                            "  603.64,1403.28 610.298,1403.28 616.956,1403.28 623.615,1403.28 630.273,1403.28 636.931,1403.28 643.589,1403.28 650.247,1403.27 656.906,1403.16 663.564,1402.93 \n",
                            "  670.222,1402.6 676.88,1402.16 683.539,1401.61 690.197,1400.97 696.855,1400.22 703.513,1399.38 710.171,1398.43 716.83,1397.39 723.488,1396.26 730.146,1395.03 \n",
                            "  736.804,1393.72 743.462,1392.31 750.121,1390.81 756.779,1389.22 763.437,1387.55 770.095,1385.8 776.754,1383.96 783.412,1382.04 790.07,1380.04 796.728,1377.96 \n",
                            "  803.386,1375.8 810.045,1373.56 816.703,1371.25 823.361,1368.87 830.019,1366.41 836.678,1363.87 843.336,1361.27 849.994,1358.59 856.652,1355.85 863.31,1353.04 \n",
                            "  869.969,1350.16 876.627,1347.22 883.285,1344.21 889.943,1341.14 896.601,1338 903.26,1334.8 909.918,1331.54 916.576,1328.22 923.234,1324.84 929.893,1321.4 \n",
                            "  936.551,1317.91 943.209,1314.36 949.867,1310.75 956.525,1307.09 963.184,1303.37 969.842,1299.6 976.5,1295.78 983.158,1291.91 989.816,1287.98 996.475,1284.01 \n",
                            "  1003.13,1279.99 1009.79,1275.92 1016.45,1271.8 1023.11,1267.63 1029.77,1263.42 1036.42,1259.16 1043.08,1254.86 1049.74,1250.51 1056.4,1246.12 1063.06,1241.68 \n",
                            "  1069.72,1237.21 1076.37,1232.69 1083.03,1228.13 1089.69,1223.54 1096.35,1218.9 1103.01,1214.22 1109.66,1209.51 1116.32,1204.75 1122.98,1199.97 1129.64,1195.14 \n",
                            "  1136.3,1190.28 1142.96,1185.38 1149.61,1180.45 1156.27,1175.48 1162.93,1170.48 1169.59,1165.45 1176.25,1160.38 1182.9,1155.28 1189.56,1150.15 1196.22,1144.99 \n",
                            "  1202.88,1139.8 1209.54,1134.58 1216.2,1129.33 1222.85,1124.05 1229.51,1118.74 1236.17,1113.41 1242.83,1108.04 1249.49,1102.65 1256.15,1097.23 1262.8,1091.79 \n",
                            "  1269.46,1086.32 1276.12,1080.83 1282.78,1075.31 1289.44,1069.76 1296.09,1064.19 1302.75,1058.6 1309.41,1052.98 1316.07,1047.34 1322.73,1041.68 1329.39,1036 \n",
                            "  1336.04,1030.29 1342.7,1024.57 1349.36,1018.82 1356.02,1013.05 1362.68,1007.27 1369.33,1001.46 1375.99,995.63 1382.65,989.784 1389.31,983.919 1395.97,978.035 \n",
                            "  1402.63,972.133 1409.28,966.214 1415.94,960.277 1422.6,954.323 1429.26,948.352 1435.92,942.365 1442.58,936.362 1449.23,930.343 1455.89,924.308 1462.55,918.258 \n",
                            "  1469.21,912.193 1475.87,906.114 1482.52,900.02 1489.18,893.913 1495.84,887.791 1502.5,881.656 1509.16,875.508 1515.82,869.347 1522.47,863.173 1529.13,856.987 \n",
                            "  1535.79,850.789 1542.45,844.579 1549.11,838.357 1555.76,832.124 1562.42,825.88 1569.08,819.625 1575.74,813.36 1582.4,807.084 1589.06,800.797 1595.71,794.501 \n",
                            "  1602.37,788.196 1609.03,781.881 1615.69,775.556 1622.35,769.223 1629.01,762.881 1635.66,756.53 1642.32,750.171 1648.98,743.804 1655.64,737.429 1662.3,731.046 \n",
                            "  1668.95,724.655 1675.61,718.258 1682.27,711.853 1688.93,705.441 1695.59,699.023 1702.25,692.597 1708.9,686.166 1715.56,679.728 1722.22,673.285 1728.88,666.835 \n",
                            "  1735.54,660.38 1742.19,653.92 1748.85,647.454 1755.51,640.983 1762.17,634.507 1768.83,628.026 1775.49,621.541 1782.14,615.051 1788.8,608.557 1795.46,602.059 \n",
                            "  1802.12,595.557 1808.78,589.051 1815.44,582.541 1822.09,576.028 1828.75,569.511 1835.41,562.992 1842.07,556.469 1848.73,549.943 1855.38,543.414 1862.04,536.883 \n",
                            "  1868.7,530.349 1875.36,523.812 1882.02,517.274 1888.68,510.733 1895.33,504.19 1901.99,497.646 1908.65,491.099 1915.31,484.551 1921.97,478.002 1928.62,471.451 \n",
                            "  1935.28,464.899 1941.94,458.346 1948.6,451.791 1955.26,445.236 1961.92,438.68 1968.57,432.123 1975.23,425.566 1981.89,419.008 1988.55,412.45 1995.21,405.891 \n",
                            "  2001.87,399.333 2008.52,392.774 2015.18,386.215 2021.84,379.657 2028.5,373.099 2035.16,366.541 2041.81,359.984 2048.47,353.427 2055.13,346.871 2061.79,340.316 \n",
                            "  2068.45,333.762 2075.11,327.21 2081.76,320.658 2088.42,314.108 2095.08,307.561 2101.74,301.015 2108.4,294.472 2115.05,287.932 2121.71,281.396 2128.37,274.865 \n",
                            "  2135.03,268.341 2141.69,261.823 2148.35,255.316 2155,248.821 2161.66,242.343 2168.32,235.886 2174.98,229.457 2181.64,223.066 2188.3,216.724 2194.95,210.449 \n",
                            "  2201.61,204.261 2208.27,198.192 2214.93,192.28 2221.59,186.58 2228.24,181.164 2234.9,176.129 2241.56,171.605 2248.22,167.769 2254.88,164.857 2261.54,163.187 \n",
                            "  \n",
                            "  \"/>\n",
                            "<polyline clip-path=\"url(#clip9303)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
                            "  650.247,132890 650.247,-131324 \n",
                            "  \"/>\n",
                            "</svg>\n"
                        ]
                    },
                    "execution_count": 4,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "lb = zeros(M)\n",
                "ub = 300*ones(M) # Need to have upper bounds for the z.\n",
                "code, z, w_val = @suppress solveLCP(w, lb, ub)  # Solves (12)\n",
                "@assert code == :Solved  # otherwise, an error in convergence\n",
                "\n",
                "v = z + S.(x)  # (8)\n",
                "x̂_index = findfirst(z .> 0)  # i.e. find first x index where non-binding\n",
                "if isnothing(x̂_index)\n",
                "  x̂ = Inf # i.e. never non-binding = always stop\n",
                "else\n",
                "  x̂ = x[x̂_index]\n",
                "end\n",
                "\n",
                "# plot the solution and the stopping threshold\n",
                "plot(x, v, title = \"Value Function for Optimal Stopping\",\n",
                "     legend = false, ylabel = \"v(x)\", xlabel = \"x\")\n",
                "vline!([x̂])"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Performance Note** As always in Julia, only test performance of code within functions (i.e. the global variables here are only for simplicity)\n",
                "\n",
                "### References\n",
                "* [Github Source](https://github.com/QuantEcon/SimpleDifferentialOperators.jl/blob/master/docs/examples/LCP_simple.jmd)\n",
                "* [Huang, J. and J.-S. Pang (1998): \"Option Pricing and Linear Complementarity,\"](https://ecommons.cornell.edu/handle/1813/5453)\n",
                "* [Ben Moll's notes](http://www.princeton.edu/~moll/HACTproject/option_simple.pdf)"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Julia 1.1.0",
            "language": "julia",
            "name": "julia-1.1"
        },
        "language_info": {
            "file_extension": ".jl",
            "mimetype": "application/julia",
            "name": "julia",
            "version": "1.1.0"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}