Lecture9.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "lecture = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#Lecture 9: Linear Programming and ALM: Cash-Flow Matching\n",
    "\n",
    "## Topics\n",
    "* I. Linear Programming Theory\n",
    "\n",
    "    * Formulation\n",
    "    * Duality\n",
    "    * Solution Methods\n",
    "        * Simplex Method\n",
    "        * Primal-Dual Interior-Point Method "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* II. Asset Liability Management: Cash-Flow Matching\n",
    "\n",
    "    * Basic setup\n",
    "    * Example\n",
    "    * Nonlinear scenario\n",
    "    * Stochastic liabilities\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# I. Linear Programming Theory\n",
    "\n",
    "\n",
    "## Formulation\n",
    "\n",
    "\\begin{aligned}\n",
    "\\renewcommand{bx}{\\boldsymbol x}\n",
    "\\renewcommand{bs}{\\boldsymbol}\n",
    "\\min_{\\bx \\in \\mathbb{R}^n } \\;\\; &z = \\bs{c}^T \\bx\n",
    "\\\\\n",
    "s.t. \\;\\;\\bs{A}\\bx  &\\geq \\bs{b};\\;\n",
    "x_j \\geq 0, \\;\\forall j\n",
    "\\end{aligned}\n",
    "\n",
    "\n",
    "* Simpler than any other optimization problem\n",
    "\n",
    "* With plenty of application (even in Finance)\n",
    "\n",
    "* Serve as a good starting point\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Duality\n",
    "\n",
    "Primal (P)\n",
    "\n",
    "\\begin{aligned}\n",
    "\\min_{\\bx \\in \\mathbb{R}^n } \\;\\; &z = \\bs{c}^T \\bx\n",
    "\\\\\n",
    "s.t. \\;\\;\\bs{A}\\bx  \\geq \\bs{b};\\;\n",
    "& x_j \\geq 0, \\;\\forall j = \\overline{1, n}\n",
    "\\end{aligned}\n",
    "\n",
    "\n",
    "\n",
    "Dual (D)\n",
    "\n",
    "\\begin{aligned}\n",
    "\\max_{\\bs{y} \\in \\mathbb{R}^m } \\;\\; & z = \\bs{b}^T \\bs{y}\n",
    "\\\\\n",
    "s.t. \\;\\;\\bs{A}^T\\bs{y}  \\leq \\bs{c}; \\;\\;\n",
    "& y_k \\geq 0, \\;\\; \\forall k = \\overline{1, m}\n",
    "\\end{aligned}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Derivation\n",
    "\n",
    "Lagrangian:\n",
    "\n",
    "$$\n",
    "\\renewcommand{ml}{\\mathbb{\\mathcal L}}\n",
    "\\renewcommand{bA}{\\boldsymbol A}\n",
    "\\renewcommand{by}{\\boldsymbol y}\n",
    "\\renewcommand{bb}{\\boldsymbol b}\n",
    "\\renewcommand{bc}{\\boldsymbol c}\n",
    "\\renewcommand{bw}{\\boldsymbol w}\n",
    "\\begin{array}\n",
    "\\\\\n",
    "\\ml(\\bx, \\by, \\bw) & = & \\bc^T \\bx + \\by^T (\\bb-\\bA\\bx) - \\bw^T\\bx\n",
    "\\\\\n",
    "& = & \\by^T\\bb + (\\bc^T  - \\by^T \\bA - \\bw^T)\\bx\n",
    "\\\\\n",
    "& & \\mbox{with} \\;\\; \\by \\geq \\bf{0}, \\;\\;\\; \\bw \\geq \\bf{0}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Lagrangian dual function\n",
    "\n",
    "$$\n",
    "\\renewcommand{mD}{\\mathbb{\\mathcal D}}\n",
    "\\renewcommand{df}{\\hat{f}}\n",
    "\\df(\\by, \\bw) = \\inf_{\\bx\\in \\mD}\\; \\ml(\\bx, \\by, \\bw) =\\by^T\\bb + \\inf_{\\bx\\in\\mD}  (\\bc^T  - \\by^T \\bA - \\bw^T)\\bx\n",
    "$$\n",
    "\n",
    "Obviously the last term is linear in $\\bx$, so the only scenario where the infimum can be reached is when the slop is zero:\n",
    "\n",
    "$$\n",
    "\\bc  - \\bA^T \\by - \\bw = 0\n",
    "$$\n",
    "\n",
    "write-off the slack variable $\\bw$ will generate the dual problem (D).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Duality Theorem:\n",
    "\n",
    "   * If $\\bx$ is feasible for (P) and $\\bs{y}$ is feasible for (D), then $\\bs{c}^T \\bx \\geq \\bs{b}^T \\bs{y}$, \n",
    "    and only at optimality the equality holds $\\bs{c}^T \\bx^* = \\bs{b}^T \\bs{y}^*$ .\n",
    "\n",
    "\n",
    "* The quantity: $\\bs{c}^T \\bx - \\bs{b}^T \\bs{y}$  is called **duality gap**. This value is used commonly as the stopping criteria for Interior LP method.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Solution Methods\n",
    "\n",
    "* **Simplex** (Dantzig 1947)\n",
    "\n",
    "* **Ellipsoid** (Kachian 1979, the first algorithm known to be in **polynomial time**)\n",
    "\n",
    "* **Interior Point** (Karmakar 1984, the first  **practical** polynomial time algorithm)\n",
    "    * Projection method (Karmakar 1984)\n",
    "    * Affine Method (Dikin 1967)\n",
    "    * Log-Barrier Method (many ...)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "* **Interior Point** method has been extended to NLP problems, has been the focus of research for optimization in the last few decades.\n",
    "\n",
    "* Although asymptotically superior, there is no clear winner between Simplex and Interior Point for LP problems: very much depends on the problem\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Simplex Methods\n",
    "\n",
    "* Optimal points must be at the intersections of constraints:\n",
    "\n",
    "<center><img src = \"img/LP_GeomView.png\" height=600 width=600></center>\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Algorithm (Simplex Method)\n",
    "\n",
    "* Input\n",
    "\\begin{aligned}\n",
    "\\\\\n",
    "\\min_{\\bx \\in \\mathbb{R}^n } \\;\\;   z & = \\bs{c}^T \\bx\n",
    "\\\\\n",
    "s.t. \\;\\; \\bx \\in \\mathcal{P} & = \\{\\bx: \\bs{A}\\bx  \\leq \\bs{b}, \\bx \\geq 0 \\}\n",
    "\\end{aligned}\n",
    "\n",
    "\n",
    "* Look at the polytope $\\mathcal{P}$ as a graph $\\mathcal{G} = (V,E)$ with $V$ as vertices and $E$ as edges.\n",
    "\n",
    "\n",
    "* **Algorithm**\n",
    "    1. Find any vertex $\\bs{v}$ in  $\\mathcal{P}$;\n",
    "    2. While $\\exists $ a neighbor $\\bs{u}$ of $\\bs{v}$ in $\\mathcal{G}$ with $\\bs{c}^T \\bs{u} < \\bs{c}^T \\bs{v}$, set $\\bs{v} = \\bs{u}$;\n",
    "    3. Output $\\bs{v}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Primal-dual Interior-point Method\n",
    "\n",
    "* For convenience, reformulate the LP problme into an equivalent form\n",
    "\n",
    "Primal (P)\n",
    "\n",
    "\\begin{aligned}\n",
    "\\min_{\\bx \\in \\mathbb{R}^n } \\;\\; &  \\bs{c}^T \\bx\n",
    "\\\\\n",
    "s.t. \\;\\;  \\bs{A}\\bx  &= \\bs{b}, \\;\\; \\bx \\geq \\bf{0}\n",
    "\\end{aligned}\n",
    "\n",
    "\n",
    "Dual (D)\n",
    "\\begin{aligned}\n",
    "\\max_{\\bs{y} \\in \\mathbb{R}^m } \\;\\; & \\bs{b}^T \\bs{y}\n",
    "\\\\\n",
    "s.t. \\;\\; \\bs{A}^T\\bs{y} + \\bs{s} & = \\bs{c}, \\;\\; \\bs{s} \\geq \\bf{0}\n",
    "\\end{aligned}\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* And the optimality (KKT) conditions are\n",
    "\n",
    "\\begin{aligned} \n",
    " \\bs{A}^T\\bs{y}  + \\bs{s} &= \\bs{c}, \n",
    "\\\\\n",
    " \\bs{A}\\bx  &= \\bs{b}, \n",
    "\\\\\n",
    " \\bx\\circ \\bs{s} &= \\bf{0}, \n",
    "\\\\\n",
    "\\bx \\geq \\bf{0}, \\;\\; \\bs{s} &\\geq \\bf{0}\n",
    "\\end{aligned}\n",
    "\n",
    "The third complementary slackness condition relates to the duality gap.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "Primal-dual interior-point methods find primal-dual solutions $(\\bx^*, \\by^*, \\bs{s}^*)$ by applying variants of Newton's method to the first three KKT conditions above, and modifying search directions and step lengths, so that the last KKT condition is strictly satisfied at each iteration.\n",
    "\n",
    "Rewrite the KKT conditions as:\n",
    "$$\n",
    "\\begin{matrix}\n",
    "\\bs{F}(\\bx,\\by,\\bs{s}) = \\left[\n",
    "\\begin{matrix}\n",
    "\\bs{A}^T\\bs{y}  + \\bs{s} - \\bs{c} \\\\\n",
    "\\bs{A}\\bx  - \\bs{b} \\\\\n",
    "\\bs{X} \\bs{S} \\bs{e}\n",
    "\\end{matrix}\n",
    "\\right] = \\bs{0}\n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "where $\\bs{X}, \\bs{S}$ are diagonal matrices, and $\\bs{e}$ is the vector with all elements equal to 1.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "And the search direction is generated by solving,\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    "\\bs{F}'(\\bx,\\by,\\bs{s})  \\left[\n",
    "\\begin{matrix}\n",
    "\\delta\\bx \\\\\n",
    "\\delta\\by \\\\\n",
    "\\delta\\bs{s}\n",
    "\\end{matrix}\n",
    "\\right] = \\bs{F}(\\bx,\\by,\\bs{s})\n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "Looking for iterates that are strictly feasible, this becomes\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "0 & \\bs{A}^T & \\bs{I} \\\\\n",
    "\\bs{A}^T & 0 & 0  \\\\\n",
    "\\bs{S} & 0 & \\bs{X} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\delta\\bx \\\\\n",
    "\\delta\\by \\\\\n",
    "\\delta\\bs{s}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "= \\left[\n",
    "\\begin{matrix}\n",
    "\\bs{0} \\\\\n",
    "\\bs{0} \\\\\n",
    "-\\bs{X} \\bs{S} \\bs{e}\n",
    "\\end{matrix}\n",
    "\\right] \n",
    "\\end{matrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "In Interior-point method, the critical concept is the central path: instead of the third KKT condition, the path of iterative points satisfies,\n",
    "\n",
    "$$\n",
    "x_i s_i = \\tau\\gamma, \\;\\;\\; i=1,2,\\cdots, n\n",
    "$$\n",
    "\n",
    "where $\\tau\\in  [0,1] $ is the centering parameter and $\\gamma = \\frac{\\bs{s}^T\\bx}{n}$ is the duality gap parameter.\n",
    "\n",
    "With this, the equation for search direction becomes\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "0 & \\bs{A}^T & \\bs{I} \\\\\n",
    "\\bs{A}^T & 0 & 0  \\\\\n",
    "\\bs{S} & 0 & \\bs{X} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\delta\\bx \\\\\n",
    "\\delta\\by \\\\\n",
    "\\delta\\bs{s}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "= \\left[\n",
    "\\begin{matrix}\n",
    "\\bs{0} \\\\\n",
    "\\bs{0} \\\\\n",
    "-\\bs{X} \\bs{S} \\bs{e} + \\tau\\gamma \\bs{e}\n",
    "\\end{matrix}\n",
    "\\right] \\hspace{2in} (*)\n",
    "\\end{matrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### **Algorithm** (Primal-dual Interior-point method)\n",
    "\n",
    "1. Determine $(\\bx^0, \\by^0, \\bs{s}^0)$ strictly feasible\n",
    "2. Set k = 0\n",
    "3. Do\n",
    "    1. Set $\\tau^k\\in [0,1]$ and $\\gamma^k = \\frac{\\bs{s}^T\\bx}{n}$\n",
    "    2. Solve (*) for $(\\delta\\bx^k, \\delta\\by^k, \\delta\\bs{s}^k)$\n",
    "    3. Set $(\\bx^{k+1}, \\by^{k+1}, \\bs{s}^{k+1}) = (\\bx^k, \\by^k, \\bs{s}^k) + \\alpha^k(\\delta\\bx^k, \\delta\\by^k, \\delta\\bs{s}^k)$\n",
    "       choosing $\\alpha^k$ so that $(\\bx^{k+1}, \\bs{s}^{k+1}) \\geq 0 $\n",
    "    4. k = k + 1\n",
    "4. While  $\\gamma^k > \\epsilon$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "* Software packages for LP\n",
    "    * Most optimization packages will have the LP solver implemented, and most likely include both the Simplex and Interior method.\n",
    "    * Scipy just released the first LP solver (in Simplex): Scipy version 0.15.1 (with Numpy version: 1.9.2)\n",
    "    * Here is the code applying the Scipy LP optimizer to the example in the geometric view slide:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: -38.000000  \n",
      "         Iterations: 2\n",
      "The solution from Scipy Simplex method: [ 4. 10.]\n",
      "The Scipy version:\n",
      "1.1.0\n",
      "The Numpy version:\n",
      "1.14.3\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "\n",
    "c = np.array([-2, -3])\n",
    "A = [[1, -2], [2, 1], [0, 1]]\n",
    "b = np.array([4,18,10])\n",
    "res = linprog(c, A, b,\n",
    "                options={ 'disp': True})\n",
    "\n",
    "print(\"The solution from Scipy Simplex method:\", res.x)\n",
    "import scipy\n",
    "print(\"The Scipy version:\")\n",
    "print(scipy.__version__)\n",
    "print(\"The Numpy version:\")\n",
    "print(np.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## II. Asset Liability Management: Cash-FLow Matching\n",
    "\n",
    "* Asset Liability Management (ALM) is the practice of managing a business so that decisions and actions taken with respect to assets and liabilities are coordinated. It's the sound management of the finances to match the asset cash flows with its future liability cash flow needs and capital requirements. \n",
    "\n",
    "* ALM practice now considers a wide range of risks: liquidity, interest rate, equity, currency, legal, etc.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We will be considering the interest rate risk part of the ALM.\n",
    "\n",
    "* ALM usually takes the form of \"immunization\": matching the present value, duration and convexity of the liabilities and assets in the portfolio. This is typically done at each time bucket of the planning horizon.\n",
    "\n",
    "* Here, we will focus on the form of \"dedication\": choose an investment asset portfolio such that, for each period of the planning horizon, the investment cash flow will be sufficient to meet the projected liability cash flows. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basic Setup for \"Dedication\"\n",
    "\n",
    "Given:\n",
    "\n",
    "*  $t_k, \\;\\; k = 1,\\cdots, K $ are the time periods\n",
    "*  $L_k, \\;\\; k = 1,\\cdots, K $ are the future liability payment at time $t_k$\n",
    "\n",
    "*  $P_j, \\;\\; j = 1,\\cdots, N $ are the prices of $N$ investable assets\n",
    "*  $C_{j, k}, \\;\\; j = 1,\\cdots, N, \\;\\; k= 1, \\cdots, K $ are the cash flows of asset $j$ at time  period $t_k$\n",
    "\n",
    "\n",
    "Assuming:\n",
    "\n",
    "*  $x_j, \\;\\; j = 1,\\cdots, N $ are the amount to be invested in asset $j$\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "\n",
    "The problem is to find the asset portfolio such that the total cost\n",
    "\n",
    "$$\n",
    "\\sum_j x_j P_j \n",
    "$$\n",
    "\n",
    "is minimized, while the projected liability cash flows are covered by the cash flows generated by the asset portfolio\n",
    "\n",
    "$$\n",
    "\\sum_j x_j C_{j,k} \\geq L_k \\;\\;  k = 1,\\cdots,K\n",
    "$$\n",
    "\n",
    "and $x_j \\geq 0$.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "\n",
    "\n",
    "* A pension fund has estimated that for the next 8 years it has the following schedule of benefits to be paid (in $millions ):\n",
    "\n",
    "\n",
    "|  Years | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |\n",
    "| :-----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: |\n",
    "|  Benefits($millions) | 12 | 18 | 20 | 20 | 16 | 15 | 18 | 20 |\n",
    "\n",
    "* The fund's investment advisors recommended the following set of bonds to invest:\n",
    "\n",
    "|  Bonds | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |\n",
    "| :-----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: | :----: |\n",
    "|  Price | 102 | 99 | 101 | 98 | 98 | 104 | 100 | 101 | 102 | 94 |\n",
    "|  Coupon(%) | 5 | 6.5 | 5 | 3.5 | 4 | 9 | 6 | 8 | 9 | 7 |\n",
    "|  Maturity(Years) | 1 | 2 | 2 | 3 | 4 | 5 | 5 | 6 | 7 | 8 |\n",
    "\n",
    "* The bonds all have face values $\\$100$, the coupons are annualized percentage.\n",
    "\n",
    "* Can you help find the **least cost** portfolio of bonds so that the pension fund can meet its future benefit payments?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### LP formulation set up:\n",
    "\n",
    "\n",
    "* $L_k, k = 1,\\cdots, 8 $ are the future benefits (liabilities) payment at year $k$\n",
    "* $P_i, i = 1,\\cdots, 10 $ are the prices of bond $i$\n",
    "* $C_i, i = 1,\\cdots, 10 $ are the coupons of bond $i$\n",
    "* $M_i, i = 1,\\cdots, 10 $ are the maturity of bond $i$\n",
    "\n",
    "\n",
    "* We will assume the variables are:\n",
    "    * $x_i, i = 1,\\cdots, 10 $ are the amount of bond $i$ to be purchased\n",
    "    * $s_k, k = 1,\\cdots, 8 $ are the free cash at the end of year $k$ \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "The LP problem can be formulated as:\n",
    "\n",
    "$$\n",
    "\\min_{\\bx} \\left( \\sum_i x_i P_i \\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mbox{s.t.} & \\sum_{i=1,\\cdots,10} ({\\bf 1}_{M_i\\geq k}\\cdot C_i(t) + {\\bf 1}_{M_i=k}\\cdot 100)\\; x_i - s_k + s_{k-1}= L_k, \\;\\;  k = 1,2,\\cdots, 8\n",
    "\\\\\n",
    "     & x_i \\geq 0 \\;\\; i = 1,2,\\cdots, 10; \\;\\;\\; s_k \\geq 0 \\;\\;  k = 1,2,\\cdots, 8\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where we assume $s_0 = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 101.615025  \n",
      "         Iterations: 9\n",
      "The solution from Scipy Simplex method:\n",
      " [ 0.04860098  0.2588899   0.          0.          0.14785887  0.11377323\n",
      "  0.          0.11401282  0.15313384  0.18691589  0.         14.78588715\n",
      "  0.          0.          0.          0.          0.          0.        ]\n"
     ]
    }
   ],
   "source": [
    "p = np.array([102, 99, 101, 98, 98, 104, 100, 101, 102, 94, 0, 0, 0, 0, 0, 0, 0, 0])\n",
    "b = np.array([12, 18, 20, 20, 16, 15, 18, 20])\n",
    "A = [[100+5, 6.5, 5, 3.5, 4, 9, 6, 8, 9, 7, -1.0, 0, 0, 0, 0, 0, 0, 0], \n",
    "     [0, 100+6.5, 100+5, 3.5, 4, 9, 6, 8, 9, 7, 1.0, -1.0, 0, 0, 0, 0, 0, 0],  \n",
    "     [0, 0, 0, 100+3.5, 4, 9, 6, 8, 9, 7, 0, 1.0, -1.0, 0, 0, 0, 0, 0],  \n",
    "     [0, 0, 0, 0, 100+4, 9, 6, 8, 9, 7, 0, 0, 1.0, -1.0, 0, 0, 0, 0],  \n",
    "     [0, 0, 0, 0, 0, 100+9, 100+6, 8, 9, 7, 0, 0, 0, 1.0, -1.0, 0, 0, 0],  \n",
    "     [0, 0, 0, 0, 0, 0, 0, 100+8, 9, 7, 0, 0, 0, 0, 1.0, -1.0, 0, 0],  \n",
    "     [0, 0, 0, 0, 0, 0, 0, 0, 100+9, 7, 0, 0, 0, 0, 0, 1.0, -1.0, 0],  \n",
    "     [0, 0, 0, 0, 0, 0, 0, 0, 0, 100+7, 0, 0, 0, 0, 0, 0, 1.0, -1.0] ]\n",
    "\n",
    "res = linprog(p, A_eq=A, b_eq=b,\n",
    "                options={ 'disp': True})\n",
    "\n",
    "print(\"The solution from Scipy Simplex method:\\n\", res.x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Extensions of the Asset liability cash flow matching problem\n",
    "\n",
    "\n",
    "The constraint condition can be turned into equality based by assuming the surplus free cash flow at each time period and consider them to be variables\n",
    "\n",
    "*  $s_k, \\;\\; k = 1,\\cdots, K $ are the amount of excess cash after paying the liability\n",
    "\n",
    "With that, the constraints are now\n",
    "\n",
    "$$\n",
    "\\sum_j x_j C_{j,k} - s_k = L_k, \\;\\;  k = 1,\\cdots,K\n",
    "$$\n",
    "and $x_j \\geq 0, \\;\\; s_k \\geq 0$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* With these excess cash flows, it is then natural to assume that they will be re-invested, and the return is simply $r$, therefore the constraints becomes\n",
    "\n",
    "$$\n",
    "\\sum_j x_j C_{j,k} - s_k + (1+r) \\cdot s_{k-1} = L_k, \\;\\;  k = 1,\\cdots,K\n",
    "$$\n",
    "where we set $s_0 = 0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* The setup can be made even more realistic  by allowing the reinvestment return to be time period dependent: $r \\rightarrow r_k$\n",
    "\n",
    "* In fact, if borrowing is allowed, the constraints $ s_k \\geq 0$ for all time periods are unnecessary, except we probably want to keep the excess cash to be positive at the very last time period.\n",
    "\n",
    "* Which introduces another dimension: the borrowing rate and lending (re-investment) rate are typically different in real life. \n",
    "\n",
    "* In practice, these assumptions (or approximations) can certainly make a big difference in terms of the complexity and realistism of your model.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Undeterministic future liabilities and asset returns\n",
    "\n",
    "* Cash flow matching so far require the cash flows from the liabilities can be forcasted deterministically\n",
    "\n",
    "* This is hardly the case: liabilities of insurance companies are not fixed\n",
    "    * Policyholders can withdraw cash\n",
    "    * Hurricane frequency and severity are highly unpredictable\n",
    "    * Pension payments are hardly deterministic\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}