3. Uncertain Mortality Model and American Options.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KLnHYuBAjwA6"
   },
   "source": [
    "# Uncertain Mortality Model vs American Options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "INCWgbq8jwA7"
   },
   "source": [
    "#### Re-insurance Deals\n",
    "\n",
    "In this notebook, we study the *Uncertain Mortality Model* for the pricing of reinsurance deals. For the sake of simplicity, we only consider one type of default: the risk of death. (lapse is ignored)\n",
    "\n",
    "We consider the following reinsurance product:\n",
    "\n",
    "- At maturity, if the insurance subscriber is alive, the issuer delivers a put on the underlying $X$\n",
    "\n",
    "$$\n",
    "u^\\textrm{mat}( x) = (K_\\textrm{mat} - x)_+.\n",
    "$$\n",
    "\n",
    "- At the time of death, if it is before the maturity, the issuer delivers an exit payoff, typically another put on the underlying $X$.\n",
    "\n",
    "$$\n",
    "u^D(t, x) = (K_D - x)_+.\n",
    "$$\n",
    "\n",
    "- The subscriber pays a constant fee $α\\Delta t$ at every time step until death or the maturity or the product\n",
    "\n",
    "- The insurance sells a large number of these contracts to subscribers. We assume that the times of death of the subscribers $\\tau^D$ are independent, and identically distributed, and also independent of the underlying's stock price.\n",
    "\n",
    "- We assume that the underlying's risk neutral price dynamics is the Black-Scholes model with zero interest rate/repo/dividend yield\n",
    "\n",
    "$$\n",
    "dX_t = \\sigma X_t d W_t.\n",
    "$$\n",
    "\n",
    "#### The Insurers' Approach and Risk-Neutral Pricing\n",
    "\n",
    "This contract shows two types of risk: the times of death of the subscribers and the changes in the price of the underlying. \n",
    "\n",
    "In this case, the issuer can apply the insurer's approach to the risk of death times, i.e., the law of large numbers. The more people buy the contract, the less risk.\n",
    "\n",
    "Choosing a risk-neutral measure under which the death times $\\tau^D$ have the same distribution as under the historical probability measure is equivalent to applying the arbitrage-pricing approach to the financial risk insurer's rule on the risk of death. The price of the contract is then\n",
    "\n",
    "$$\n",
    "u(t, x) = \\mathbb{E}^\\mathbb{Q}_{t, x} \\left[u^\\textrm{mat}(X_T) \\mathbb{1}_{\\tau^D \\geq T} -\\int_t^T\\alpha\\mathbb{1}_{\\tau^D>s}ds+ u^D(\\tau^D, X_{\\tau^D}) \\mathbb{1}_{\\tau^D < T} \\right].\n",
    "$$\n",
    "\n",
    "#### Deterministic Death Rate\n",
    "\n",
    "If the death intensity is a deterministic function $\\lambda_t^D$ (i.e. $\\tau^D$ has an exponential distribution with time-dependent intensity $\\lambda_t^D$), then we have seen that $u$ is the solution to the linear PDE\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "\\partial_t u + \\frac{1}{2} \\sigma^2 x^2 \\partial_x^2 u -\\alpha + \\lambda_t^D \\cdot (u^D - u) = 0,\\\\\n",
    "u(T, \\cdot) \\equiv u^\\textrm{mat}.\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "#### Uncertain Mortality Model\n",
    "\n",
    "If the death rate is uncertain, we assume that it is adapted (i.e., it does not look into the future) and belongs to a moving corridor\n",
    "$\\left[\\underline{\\lambda}_t, \\overline{\\lambda}_t\\right]$. The most conservative way to price the contract is to compute the (financially) worst death rate process $\\lambda_t^D$ as being chosen so as to maximize the value of the contract. The resulting HJB equation is\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "\\partial_t u + \\frac{1}{2} \\sigma^2 x^2 \\partial^2_x u -\\alpha + \\Lambda^D(t, u^D - u) \\cdot (u^D - u) = 0,\\\\\n",
    "u(T, \\cdot) \\equiv u^\\textrm{mat},\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "where the function $\\Lambda^D$ is defined by\n",
    "$$\n",
    "\\Lambda^D (t, y) = \\left\\{\\begin{array}{l}\n",
    "    \\overline{\\lambda}^D_t \\quad \\textrm{if} \\ y \\geq 0, \\\\\n",
    "    \\underline{\\lambda}^D_t \\quad \\textrm{otherwise.}\n",
    "\\end{array}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2iTDNyZZjwA8"
   },
   "source": [
    "## Link with $1$-BSDE and Numerical Schemes\n",
    "\n",
    "From the Pardoux-Peng theorem, we know that the solution $u(0, x)$ can be represented as the solution $Y_0^x$ to the $1$-BSDE\n",
    "\n",
    "$$\n",
    "dY_t = -f(t, X_t, Y_t, Z_t) \\, dt + Z_t \\, dW_t,\n",
    "$$\n",
    "\n",
    "with the terminal condition $Y_T = u^\\textrm{mat} (X_T)$, where $X_0 = x$ and\n",
    "\n",
    "$$\n",
    "f(t, x, y, z) = -\\alpha + \\Lambda^D(t, u^D(t, x) - y) \\cdot (u^D(t, x) - y).\n",
    "$$\n",
    "\n",
    "\n",
    "## BSDE Discretization\n",
    "\n",
    "### Explicit Euler Schemes\n",
    "\n",
    "\\begin{align*}\n",
    "Y_{t_n} =& \\ u^{\\text{mat}}\\left(X_{t_n}\\right)\\\\\n",
    "Y_{t_{i - 1}} =& \\ \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] -\\alpha\\Delta t_i + \\Lambda^D \\left(t_{i - 1}, u^D(t_{i - 1} , X_{t_{i - 1}}) - \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ]\\right) \\cdot \\left( u^D(t_{i - 1}, X_{t_{i - 1}}) - \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] \\right) \\Delta t_i\\\\\n",
    "=& -\\alpha\\Delta t_i + \\left( \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] \\left( 1 - \\overline{\\lambda}^D \\Delta t_i \\right) +\n",
    "u^D(t_{i - 1}, X_{t_{i - 1}}) \\overline{\\lambda}^D \\Delta t_i \\right) \\mathbb{1}_{u^D\\left(t_{i - 1}, X_{t_{i - 1}}\\right) \\geq \\mathbb{E}_{i-1}[Y_{t_i}]}\\\\\n",
    "& + \\left( \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] \\left( 1 - \\underline{\\lambda}^D \\Delta t_i \\right) +\n",
    "u^D(t_{i - 1}, X_{t_{i - 1}}) \\underline{\\lambda}^D \\Delta t_i \\right) \\mathbb{1}_{u^D\\left(t_{i - 1}, X_{t_{i - 1}}\\right) < \\mathbb{E}_{i-1}[Y_{t_i}]}.\n",
    "\\end{align*}\n",
    "\n",
    "### Implicit Euler Scheme\n",
    "\n",
    "\\begin{align*}\n",
    "Y_{t_n} =& \\ u^{\\text{mat}}\\left(X_{t_n}\\right)\\\\\n",
    "Y_{t_{i - 1}} =& \\ \\mathbb{E}^\\mathbb{Q}_{i - 1} [Y_{t_i}] -\\alpha\\Delta t_i+ \\Lambda^D \\left(t_{i - 1}, u^D(t_{i - 1}, X_{t_{i - 1}}) - Y_{t_{i - 1}}\\right) \\cdot \\left( u^D(t_{i - 1}, X_{t_{i - 1}}) - Y_{t_{i - 1}} \\right) \\Delta t_i.\n",
    "\\end{align*}\n",
    "\n",
    "The implicit scheme involves $Y_{t_{i - 1}}$ on both sides of the equation, which generally requires a root-finding routine to find $Y_{t_{i - 1}}$. However in this specific case, it can be shown that \n",
    "\n",
    "\\begin{equation}u^D(t_{i-1}, X_{t_{i-1}}) \\geq Y_{t_{i-1}}\\quad\\text{if and only if}\\quad u^D(t_{i-1}, X_{t_{i-1}}) \\geq \\mathbb{E}^{\\mathbb{Q}}_{i-1}\\left[Y_{t_i}\\right]-\\alpha\\Delta t_i.\\end{equation}\n",
    "\n",
    "Thus\n",
    "\n",
    "\\begin{split}\n",
    "Y_{t_{i - 1}} =& \\frac{1}{1 + \\overline{\\lambda}^D \\Delta t_i} \\left( \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] - \\alpha\\Delta t_i +\n",
    "u^D(t_{i - 1}, X_{t_{i - 1}}) \\overline{\\lambda}^D \\Delta t_i \\right) \\mathbb{1}_{u^D\\left(t_{i - 1}, X_{t_{i - 1}}\\right) \\geq \\mathbb{E}_{i-1}[Y_{t_i}]-\\alpha\\Delta t_i}\\\\\n",
    "& + \\frac{1}{1 + \\underline{\\lambda}^D \\Delta t_i} \\left( \\mathbb{E}^\\mathbb{Q}_{i - 1} [ Y_{t_i} ] -\\alpha\\Delta t_i +\n",
    "u^D(t_{i - 1}, X_{t_{i - 1}}) \\underline{\\lambda}^D \\Delta t_i \\right) \\mathbb{1}_{u^D(t_{i - 1}, X_{t_{i - 1}}) < \\mathbb{E}_{i-1}[Y_{t_i}]-\\alpha\\Delta t_i}.\n",
    "\\end{split}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "id": "6hHfZT2fjwA-"
   },
   "outputs": [],
   "source": [
    "## Imports \n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import warnings\n",
    "from scipy.stats import norm\n",
    "from collections import deque\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "id": "GUZFEpYvjwA_"
   },
   "outputs": [],
   "source": [
    "def blackscholes_mc(S=100, vol=0.2, r=0, q=0, ts=np.linspace(0, 1, 13), npaths=10):\n",
    "    \"\"\"Generate Monte-Carlo paths in Black-Scholes model.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    S: scalar\n",
    "        The spot price of the underlying security.\n",
    "    vol: scalar\n",
    "        The implied Black-Scholes volatility.\n",
    "    r: scalar\n",
    "        The annualized risk-free interest rate, continuously compounded.\n",
    "    q: scalar\n",
    "        The annualized continuous dividend yield.\n",
    "    ts: array_like\n",
    "        The time steps of the simualtion\n",
    "    npaths: int\n",
    "        the number of paths to simulate\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    paths: ndarray\n",
    "        The Monte-Carlo paths.\n",
    "    \"\"\"\n",
    "    nsteps = len(ts) - 1\n",
    "    ts = np.asfarray(ts)[:, np.newaxis]\n",
    "    W = np.cumsum(np.vstack((np.zeros((1, npaths), dtype=np.float),\n",
    "                             np.random.randn(nsteps, npaths) * np.sqrt(np.diff(ts, axis=0)))),\n",
    "                  axis=0)\n",
    "    paths = np.exp(-0.5*vol**2*ts + vol*W)*S*np.exp((r-q)*ts)\n",
    "    return paths\n",
    "\n",
    "def blackscholes_price(K, T, S, vol, r=0, q=0, callput='call'):\n",
    "    \"\"\"Compute the call/put option price in the Black-Scholes model\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    K: scalar or array_like\n",
    "        The strike of the option.\n",
    "    T: scalar or array_like\n",
    "        The maturity of the option, expressed in years (e.g. 0.25 for 3-month and 2 for 2 years)\n",
    "    S: scalar or array_like\n",
    "        The current price of the underlying asset.\n",
    "    vol: scalar or array_like\n",
    "        The implied Black-Scholes volatility.\n",
    "    r: scalar or array_like\n",
    "        The annualized risk-free interest rate, continuously compounded.\n",
    "    q: scalar or array_like\n",
    "        The annualized continuous dividend yield.\n",
    "    callput: str\n",
    "        Must be either 'call' or 'put'.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    price: scalar or array_like\n",
    "        The price of the option.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    >>> blackscholes_price(95, 0.25, 100, 0.2, r=0.05, callput='put')\n",
    "    1.5342604771222823\n",
    "    \"\"\"\n",
    "    F = S*np.exp((r-q)*T)\n",
    "    v = np.sqrt(vol**2*T)\n",
    "    d1 = np.log(F/K)/v + 0.5*v\n",
    "    d2 = d1 - v\n",
    "    try:\n",
    "        opttype = {'call':1, 'put':-1}[callput.lower()]\n",
    "    except:\n",
    "        raise ValueError('The value of callput must be either \"call\" or \"put\".')\n",
    "    price = opttype*(F*norm.cdf(opttype*d1)-K*norm.cdf(opttype*d2))*np.exp(-r*T)\n",
    "    return price\n",
    "\n",
    "\n",
    "def pwlinear_basis(xknots):\n",
    "    \"\"\"Basis that represent a piecewise linear function with given knots\"\"\"\n",
    "    fs = [lambda x: np.ones_like(x, dtype=np.float), lambda x: x-xknots[0]]\n",
    "    fs.extend([lambda x, a=xknots[i]: np.maximum(x-a, 0) for i in range(len(xknots))])\n",
    "    return fs\n",
    "\n",
    "\n",
    "def pwlinear_fit(xdata, ydata, xknots):\n",
    "    \"\"\"Fit a piecewise linear function with xknots to xdata and ydata\"\"\"\n",
    "    fs = pwlinear_basis(xknots)\n",
    "    A = np.vstack([f(xdata) for f in fs]).T\n",
    "    ps = np.linalg.lstsq(A, ydata, rcond=None)[0]\n",
    "    return ps, fs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LeHhqW7njwBA"
   },
   "source": [
    "To be more specific about the payouts of the product, we make the following assumptions:\n",
    "\n",
    "- We assume that all transactions occur at the end of the month, i.e. the insurance subsriber is expected to pay a fee of $\\alpha/12$ at the end of each month; if the subscriber dies in the middle of the month, the death payment will be made at the end of that month. \n",
    "\n",
    "- If the death occurs at the maturity, we assume that the death comes first and the death payment $u^D$ will be made."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LxQwJZoGjwBA"
   },
   "source": [
    "<b style=\"color:darkorange\"> Example 1 : </b> (Deterministic Mortality Rate Model)\n",
    "\n",
    "We assume that the underlying's risk neutral price dynamics is the Black-Scholes model with zero interest rate/repo/dividend yield\n",
    "\n",
    "$$dX_t = \\sigma X_t d W_t,\\quad X_0=100,\\quad\\sigma=0.3$$\n",
    "\n",
    "We shall price the reinsurance deal using two Monte Carlo methods explained below :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "id": "Mq_KdRrBjwBB"
   },
   "outputs": [],
   "source": [
    "# note that lambda is a reserved keyword in python, so we use lambd instead\n",
    "def exponential_samples(n, lambd):\n",
    "    \"\"\"Generate n samples from exponential distribution with rate lambd\"\"\"\n",
    "    return -np.log(np.random.rand(n))/lambd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qI3mfrUcjwBC"
   },
   "source": [
    "<b>(a).</b> (Direct simulation of death times) Monte Carlo simulation by directly simulating death time $\\tau^D$ for each path to estimate the quantity\n",
    "\n",
    "$$u(t, x) = \\mathbb{E}^\\mathbb{Q}_{t, x} \\left[u^\\textrm{mat}(X_T) \\mathbb{1}_{\\tau^D \\geq T} - \\int_t^T\\alpha \\mathbb{1}_{\\tau^D> s}ds\n",
    "+ u^D(\\tau^D, X_{\\tau^D}) \\mathbb{1}_{\\tau^D < T} \\right],$$\n",
    "\n",
    "when $T=10$, $K_{\\text{mat}}=90$, $K_D=100$, $\\alpha=3$, and $\\tau^D$ has an exponential distribution of constant intensity $\\lambda^D=0.025$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "id": "3qFBADxajwBC"
   },
   "outputs": [],
   "source": [
    "S = 100\n",
    "vol = 0.3\n",
    "T = 10\n",
    "ts = np.linspace(0, T, int(np.round(T*12))+1)\n",
    "alpha = 3\n",
    "lambda_d = 0.025\n",
    "K_mat = 90.0\n",
    "K_D = 100.0\n",
    "u_mat = lambda x : np.maximum(K_mat - x, 0)\n",
    "u_D   = lambda x : np.maximum(K_D - x, 0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "1Z28hQGMnGsU"
   },
   "outputs": [],
   "source": [
    "npaths = 100000\n",
    "np.random.seed(2954)\n",
    "paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "tau_d = exponential_samples(npaths,lambda_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "skvegWtwm1yp",
    "outputId": "d79bbe49-2c54-4aeb-c0c8-7d46e62d0e5f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.1197249917261867"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind,_ind = np.where(tau_d < T)[0],np.where(tau_d >= T)[0]\n",
    "payoff_D = np.array([u_D(paths[np.argmax(ts>tau_d[i])][i]) - alpha*ts[np.argmax(ts>tau_d[i])] for i in ind])\n",
    "payoff_mat = u_mat(paths[-1][_ind])-alpha*T\n",
    "fairPrice = (np.sum(payoff_D) + np.sum(payoff_mat))/npaths\n",
    "fairPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-U03w54djwBE"
   },
   "source": [
    "<b>(b).</b> (Averaing over death events) Using the Feynman-Kac formula, a stochastic representation of the linear PDE \n",
    "\n",
    "$$\\left\\{\\begin{array}{l}\\partial_t u + \\frac{1}{2}\\sigma^2 x^2 \\partial_x^2 u - \\alpha + \\lambda^D \\cdot (u^D - u) = 0\\\\\n",
    "        u(T, x) \\equiv u^\\textrm{mat}(x)\\end{array}\\right.$$\n",
    "        \n",
    "that does not involve simulating death times $\\tau^D$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9NVujMTNjwBE",
    "outputId": "33ae69ed-bed6-438a-f199-dfcc5a4f0251"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.1077793704462353\n"
     ]
    }
   ],
   "source": [
    "Cp = np.zeros_like(paths)\n",
    "Cp[-1] = u_mat(paths[-1])\n",
    "dt = 1/12\n",
    "for t in range(len(paths)-2,-1,-1) :\n",
    "    Cp[t] = np.exp(-lambda_d*dt)*Cp[t+1] - alpha*dt + (1-np.exp(-lambda_d*dt))*u_D(paths[t+1])\n",
    "    \n",
    "fairPrice = np.mean(Cp[0])\n",
    "print(fairPrice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GzT1BV9HVEol"
   },
   "source": [
    "$\\text{Explanation of how we arrive at the stochastic representation :}$\n",
    "\n",
    "Consider $u \\in C^{1,2}\\left([0, T) \\times \\mathbb{R}^{d}\\right) \\cap C^{0}\\left([0, T] \\times \\mathbb{R}^{d}\\right)$ with quadratic growth in $x$ uniformly in $t$ and solution to the parabolic PDE:\n",
    "$$\n",
    "\\partial_{t} u(t, x)+\\mathcal{L} u(t, x)-r(t, x) u(t, x)+f(t, x)=0\n",
    "$$\n",
    "with terminal condition $u(T, x)=g(x)$. Then $u$ admits the stochastic representation\n",
    "\n",
    "$$u(t, x) = \\mathbb{E}^\\mathbb{Q}_{t, x} \\left[u^\\textrm{mat}(X_T) \\mathbb{1}_{\\tau^D \\geq T} - \\int_t^T\\alpha \\mathbb{1}_{\\tau^D> s}ds\n",
    "+ u^D(\\tau^D, X_{\\tau^D}) \\mathbb{1}_{\\tau^D < T} \\right],$$\n",
    "\n",
    "$$ \\mathbb{E}^\\mathbb{Q}_{t, x} \\left[u^\\textrm{mat}(X_T) \\mathbb{1}_{\\tau^D \\geq T}\\right] = e^{-\\int_{t}^{T} λ^D _s ds }\\mathbb{E}^\\mathbb{Q}_{t, x}\\left[u^{mat}(X_T)\\right] \n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\mathbb{E}^\\mathbb{Q}_{t, x}\\left[- ∫_{t}^{T} α1_{τ^D \\geq s}ds \\right] = -α∫_{t}^{T} e^{-\\int_{t}^{s} λ_{v}^{D}dv}\n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\mathbb{E}^\\mathbb{Q}_{t, x}\\left[u^D(τ^D,X_{τ^D})1_{τ^D < T} \\right] = \\mathbb{E}^\\mathbb{Q}_{t, x}\\left[ \\int_{t}^{T} u^D(s,X_s)λ_s^D e^{-\\int_t^sλ_v^D}\\right]\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\textrm{where} \\ \\mathbb{E}_{t, x}\\left[.\\right] = \\mathbb{E}\\left[.|X_t = x,τ>t\\right]   \n",
    "$$\n",
    "\n",
    "we also have\n",
    "\n",
    "$$\n",
    "u(t, x)=\\mathbb{E}^{\\mathbb{Q}}_{t, x}\\left[\\int_{t}^{T} e^{-\\int_{t}^{s} r\\left(u, X_{u}\\right) d u} f\\left(s, X_{s}\\right) d s+e^{-\\int_{t}^{T} r\\left(s, X_{s}\\right) d s} g\\left(X_{T}\\right) \\mid X_{t}=x\\right]\n",
    "$$\n",
    "With f as an intermediate payoff and the Ito generator as $\\mathcal{L}=\\frac{1}{2} \\sigma^2 x^2 \\partial_{x}^2$ here, in the Feynman-Kac formula, we put $r(t,x) = \\lambda^D$ and $f(t,x)= \\lambda^D u^D(t,x) - \\alpha$, $g=u_{mat}$ \n",
    "\n",
    "$$\n",
    "u(t, x)=\\mathbb{E}^{\\mathbb{Q}}\\left[\\int_{t}^{T} e^{-\\lambda^D \\times (s-t)} \\left(\\lambda^D u^D(t,x)- \\alpha \\right) d s+e^{-\\lambda^D \\times (T-t)} u_{mat}(X_{T}) \\mid X_{t}=x\\right]\n",
    "$$\n",
    "\n",
    "we get\n",
    "<br>\n",
    "\n",
    "  \\begin{array}{l}\\partial_t u + \\frac{1}{2}\\sigma^2 x^2 \\partial_x^2 u - \\alpha + \\lambda^D \\cdot (u^D - u) = 0\\\\\n",
    "        u(T, x) \\equiv u^\\textrm{mat}(x)\\end{array}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Aar15jweW2TB"
   },
   "source": [
    "----------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lbaKBOCYjwBE"
   },
   "source": [
    "<b style=\"color:darkorange\">Example 2 : </b> (Uncertain Mortality Rate Model - BSDE)\n",
    "\n",
    "We assume that $\\lambda^D_t \\in \\left[\\underline{\\lambda}^D, \\overline{\\lambda}^D \\right]$, where $\\underline{\\lambda}^D=0.005$ and $\\overline{\\lambda}^D=0.04$ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ucD8VAC5jwBF"
   },
   "source": [
    "<b>(a)</b> We use the implicit Euler scheme for BSDE. Using piecewise-linear fit with 10 knots to estimate the conditional expectation $\\mathbb{E}^{\\mathbb{Q}}_{i-1}\\left[Y_{t_i}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Hnjlmx5wjwBF",
    "outputId": "1d38ab0f-595c-4ef9-ed5d-7caa4b265086"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.315064524680559\n"
     ]
    }
   ],
   "source": [
    "lambda_max = 0.04\n",
    "lambda_min = 0.005\n",
    "LAMBDA_d = lambda x : lambda_max if x>=0 else lambda_min\n",
    "LAMBDA_d = np.vectorize(LAMBDA_d)\n",
    "q = deque()\n",
    "\n",
    "Yt = np.zeros_like(paths)\n",
    "Yt[-1] = u_mat(paths[-1])\n",
    "dt = 1/12\n",
    "\n",
    "for t in range(len(paths)-2,0,-1) :\n",
    "    xknots = np.linspace(np.percentile(paths[t], 1), np.percentile(paths[t], 99), 10)\n",
    "    ps, fs = pwlinear_fit(paths[t], Yt[t+1], xknots)\n",
    "    q.append((ps,fs))\n",
    "    EYt    = sum([f(paths[t])*p for (f, p) in zip(fs, ps)])\n",
    "    y = u_D(paths[t])-(EYt-alpha*dt)\n",
    "    Yt[t] = (EYt - alpha*dt + LAMBDA_d(y)*u_D(paths[t])*dt)/(1+LAMBDA_d(y)*dt)\n",
    "        \n",
    "    \n",
    "fairPrice = np.mean(Yt[1]-alpha*dt)\n",
    "print(fairPrice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oedvogy4jwBF"
   },
   "source": [
    "<b>(b)</b> Running an independent new simuation using the estimates of $\\mathbb{E}^{\\mathbb{Q}}_{i-1}\\left[Y_{t_i}\\right]$ obtained in (a) to get a lower bound estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "TnuK8C8MjwBG",
    "outputId": "cfba57a3-4ea4-4fba-8845-aa82bb4cf338"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.3073301734464806\n"
     ]
    }
   ],
   "source": [
    "### Independent Simulation \n",
    "ind_paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "Yt = np.zeros_like(ind_paths)\n",
    "Yt[-1] = u_mat(ind_paths[-1])\n",
    "for t in range(len(ind_paths)-2,0,-1) :\n",
    "    ps,fs = q.popleft()\n",
    "    EYt    = sum([f(ind_paths[t])*p for (f, p) in zip(fs, ps)])\n",
    "    y = u_D(ind_paths[t])-(EYt-alpha*dt)\n",
    "    Yt[t] = EYt - alpha*dt + LAMBDA_d(y)*y*dt\n",
    "    \n",
    "fairPrice = np.mean(Yt[1]-alpha*dt)\n",
    "print(fairPrice)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gmhQbqeNjwBG"
   },
   "source": [
    "<b>(c).</b> Comparing the price in the model with deterministic mortality rate $\\underline{\\lambda}^D$ and with deterministic mortality rate $\\overline{\\lambda}^D$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "11KzfUJYjwBG"
   },
   "source": [
    "with deterministic mortality rate $\\underline{\\lambda}^D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "5RiZCbdmjwBG",
    "outputId": "a821cae0-87d8-4fca-a1ab-7becce23dd43"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.40906682038746267\n"
     ]
    }
   ],
   "source": [
    "lambda_d = 0.005\n",
    "Cp = np.zeros_like(paths)\n",
    "Cp[-1] = u_mat(paths[-1])\n",
    "dt = 1/12\n",
    "for t in range(len(paths)-2,-1,-1) :\n",
    "    Cp[t] = np.exp(-lambda_d*dt)*Cp[t+1] - alpha*dt + (1-np.exp(-lambda_d*dt))*u_D(paths[t+1])\n",
    "    \n",
    "fairPrice = np.mean(Cp[0])\n",
    "print(fairPrice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l3bf8RQwjwBH"
   },
   "source": [
    "with deterministic mortality rate $\\overline{\\lambda}^D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9-lIe3QBjwBH",
    "outputId": "f51ce3b9-047f-4818-c020-9ea7f3995a21"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.179281680794135\n"
     ]
    }
   ],
   "source": [
    "lambda_d = 0.04\n",
    "Cp = np.zeros_like(paths)\n",
    "Cp[-1] = u_mat(paths[-1])\n",
    "dt = 1/12\n",
    "for t in range(len(paths)-2,-1,-1) :\n",
    "    Cp[t] = np.exp(-lambda_d*dt)*Cp[t+1] - alpha*dt + (1-np.exp(-lambda_d*dt))*u_D(paths[t+1])\n",
    "    \n",
    "fairPrice = np.mean(Cp[0])\n",
    "print(fairPrice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IJCby6m7jwBH"
   },
   "source": [
    "-------------------------------------------------------------------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c_Ykg4DBjwBH"
   },
   "source": [
    "<b style=\"color:darkorange\">Example 3 :</b> (Uncertain Mortality Rate Model - Longstaff and Schwartz) In addition to the BSDE simulation, the reinsurance deals can be priced by adapting the Longstaff-Schwartz algorithm\n",
    "The same lower and upper bounds of $\\lambda_t^D$ are used. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_aNvmCi7jwBH",
    "outputId": "a6e2ab98-1755-4f16-bb85-a7325ddc312d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.139202249155837"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Longstaff-Schwartz algorithm for Uncertain Mortality Rate Model\n",
    "S = 100\n",
    "vol = 0.3\n",
    "T = 10\n",
    "ts = np.linspace(0, T, int(np.round(T*12))+1)\n",
    "alpha = 3\n",
    "lambda_min = 0.005\n",
    "lambda_max = 0.04\n",
    "\n",
    "K_mat = 90.0\n",
    "K_D = 100.0\n",
    "\n",
    "def u_mat(x):\n",
    "    return np.maximum(K_mat - x, 0)\n",
    "\n",
    "def u_D(x):\n",
    "    return np.maximum(K_D - x, 0)\n",
    "\n",
    "\n",
    "# first Monte Carlo run to estimate the value function by regressions\n",
    "npaths = 5000\n",
    "betas = np.zeros((len(ts), 2), dtype=np.float)\n",
    "V = u_mat(paths[-1])\n",
    "lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = 0.5*(K_mat+K_D)\n",
    "for i in range(len(ts)-2, 0, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "    A = np.vstack((np.ones_like(Z), Z)).T\n",
    "    betas[i] = np.linalg.lstsq(A, V, rcond=None)[0]\n",
    "    lambd = np.where(u_D(paths[i]) >= betas[i, 0]+betas[i, 1]*Z, lambda_max, lambda_min)\n",
    "\n",
    "# independent simulation to obtain a lower bound\n",
    "npaths = 100000\n",
    "ind_paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "V = u_mat(ind_paths[-1])\n",
    "lambd = np.where(u_D(ind_paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = 0.5*(K_mat+K_D)\n",
    "for i in range(len(ts)-2, -1, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(ind_paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    Z = blackscholes_price(K, ts[-1]-ts[i], ind_paths[i], vol, callput='put')\n",
    "    lambd = np.where(u_D(ind_paths[i]) >= betas[i, 0]+betas[i, 1]*Z, lambda_max, lambda_min)\n",
    "np.mean(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "79SRI67RjwBI"
   },
   "source": [
    "<b>(a).</b> In the above code, we use European put price (the strike is taken to be the average of $K_{\\text{mat}}$ and $K_D$) and constant as basis functions to estimate the sum of future payoffs. We can modify the code by using the piecewise-linear fit to estimate the sum of future payoffs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "J5csGHMFjwBI",
    "outputId": "da4979c4-d9b7-4af1-c7e6-64883d8c90d5"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.203124235190343"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npaths = 5000\n",
    "betas = np.zeros((len(ts), 2), dtype=np.float)\n",
    "V = u_mat(paths[-1])\n",
    "lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = 0.5*(K_mat+K_D)\n",
    "q = deque()\n",
    "\n",
    "for i in range(len(ts)-2, 0, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "    zknots = np.linspace(np.percentile(Z, 1), np.percentile(Z, 99), 20)\n",
    "    ps, fs = pwlinear_fit(Z, V, zknots)\n",
    "    q.append((ps,fs))\n",
    "    EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "#     A = np.vstack((np.ones_like(Z), Z)).T\n",
    "#     betas[i] = np.linalg.lstsq(A, V, rcond=None)[0]\n",
    "    lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "\n",
    "# independent simulation to obtain a lower bound\n",
    "npaths = 100000\n",
    "ind_paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "V = u_mat(ind_paths[-1])\n",
    "lambd = np.where(u_D(ind_paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = 0.5*(K_mat+K_D)\n",
    "for i in range(len(ts)-2, -1, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(ind_paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    if i!=0 :\n",
    "        Z = blackscholes_price(K, ts[-1]-ts[i], ind_paths[i], vol, callput='put')\n",
    "        ps,fs = q.popleft()\n",
    "        EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "        lambd = np.where(u_D(ind_paths[i]) >= EV, lambda_max, lambda_min)\n",
    "np.mean(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nRje0JEfjwBI"
   },
   "source": [
    "<b>(b).</b> We compare the BSDE scheme with the Longstaff-Schwartz-like scheme in terms of variance. For each method we run the two-step procedure (regression and indedenpent pricing) 20 times to estimate variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "id": "CkBd9vxLjwBI"
   },
   "outputs": [],
   "source": [
    "bsde_values = []\n",
    "ls_values = []\n",
    "\n",
    "for iteration in range(20) :\n",
    "    npaths = 5000\n",
    "    paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "    betas = np.zeros((len(ts), 2), dtype=np.float)\n",
    "    V = u_mat(paths[-1])\n",
    "    lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "    lambda_max = 0.04\n",
    "    lambda_min = 0.005\n",
    "    LAMBDA_d = lambda x : lambda_max if x>=0 else lambda_min\n",
    "    LAMBDA_d = np.vectorize(LAMBDA_d)  \n",
    "    K = 0.5*(K_mat+K_D)\n",
    "    Yt = np.zeros_like(paths)\n",
    "    Yt[-1] = u_mat(paths[-1])\n",
    "    ls_q = deque()\n",
    "    bsde_q = deque()\n",
    "\n",
    "    for i in range(len(ts)-2, 0, -1):\n",
    "        dt = ts[i+1]-ts[i]\n",
    "        p = 1-np.exp(-lambd*dt)\n",
    "        V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "        Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "        zknots = np.linspace(np.percentile(Z, 1), np.percentile(Z, 99), 20)\n",
    "        ps, fs = pwlinear_fit(Z, V, zknots)\n",
    "        ls_q.append((ps,fs))\n",
    "        EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "    #     A = np.vstack((np.ones_like(Z), Z)).T\n",
    "    #     betas[i] = np.linalg.lstsq(A, V, rcond=None)[0]\n",
    "        lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "\n",
    "    for t in range(len(paths)-2,0,-1) :\n",
    "        xknots = np.linspace(np.percentile(paths[t], 1), np.percentile(paths[t], 99), 10)\n",
    "        ps, fs = pwlinear_fit(paths[t], Yt[t+1], xknots)\n",
    "        bsde_q.append((ps,fs))\n",
    "        EYt    = sum([f(paths[t])*p for (f, p) in zip(fs, ps)])\n",
    "        y = u_D(paths[t])-(EYt-alpha*dt)\n",
    "        Yt[t] = (EYt - alpha*dt + LAMBDA_d(y)*u_D(paths[t])*dt)/(1+LAMBDA_d(y)*dt)\n",
    "        \n",
    "\n",
    "\n",
    "    # independent simulation to obtain a lower bound\n",
    "    npaths = 50000\n",
    "    paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "    V = u_mat(paths[-1])\n",
    "    lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "    K = 0.5*(K_mat+K_D)\n",
    "    for i in range(len(ts)-2, -1, -1):  ## iteration for longstaff-schwartz algo\n",
    "        dt = ts[i+1]-ts[i]\n",
    "        p = 1-np.exp(-lambd*dt)\n",
    "        V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "        if i!=0 :\n",
    "            Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "            ps,fs = ls_q.popleft()\n",
    "            EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "            lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "    ls_values.append(np.mean(V))\n",
    "\n",
    "    Yt = np.zeros_like(paths)\n",
    "    Yt[-1] = u_mat(paths[-1])\n",
    "    for t in range(len(paths)-2,0,-1) :  ## iteration for bsde algo\n",
    "        ps,fs = bsde_q.popleft()\n",
    "        EYt    = sum([f(paths[t])*p for (f, p) in zip(fs, ps)])\n",
    "        y = u_D(paths[t])-(EYt-alpha*dt)\n",
    "        Yt[t] = EYt - alpha*dt + LAMBDA_d(y)*y*dt\n",
    "        \n",
    "    bsde_values.append(np.mean(Yt[1]-alpha*dt))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "NG5qCnteh1Uw",
    "outputId": "6d3c9e5a-2db3-4532-9630-e5b1223d3649"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Longstff-Schwartz Algorithm -> Mean :  3.2294429939252227 Variance :  0.011764633785905572\n",
      "BSDE Algorithm -> Mean :  3.2761140369045014 Variance :  0.07019215170355601\n"
     ]
    }
   ],
   "source": [
    "## Calculating Mean and Variance for both methods  \n",
    "print(\"Longstff-Schwartz Algorithm -> Mean : \",np.mean(ls_values),\"Variance : \",np.std(ls_values)**2)\n",
    "print(\"BSDE Algorithm -> Mean : \",np.mean(bsde_values),\"Variance : \",np.std(bsde_values)**2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HPNAshPtibD5"
   },
   "source": [
    "We observe that Longstaff-Schwartz Algorithm has lower variance as compared to BSDE Algorithm.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rnsP5XGajwBJ"
   },
   "source": [
    "<b>(c).</b> The value of the fee $\\alpha$ that makes the deal costless at inception : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_-mO-LiRjwBJ",
    "outputId": "129c42e8-192d-4590-d47f-432b49f255dd"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.3710205078125 0.0\n"
     ]
    }
   ],
   "source": [
    "## We first observe a positive price for alpha=3.3 and negative value for alpha=3.4\n",
    "## Now we perform a binary search between 3.3 and 3.4 till we get fair price of 0 (approx. to 4 decimal places)\n",
    "## We also chose two different seeds for independent simulation otherwise our data will keep changing and we will never be able to reach 0 value since the \n",
    "## fair value will also keep changing\n",
    "\n",
    "value = 1\n",
    "alpha_max = 3.4\n",
    "alpha_min = 3.3\n",
    "\n",
    "while value!=0 :\n",
    "    alpha = (alpha_max+alpha_min)/2\n",
    "    npaths = 5000\n",
    "    np.random.seed(2954)\n",
    "    paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "    betas = np.zeros((len(ts), 2), dtype=np.float)\n",
    "    V = u_mat(paths[-1])\n",
    "    lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "    K = 0.5*(K_mat+K_D)\n",
    "    q = deque()\n",
    "\n",
    "    for i in range(len(ts)-2, 0, -1):\n",
    "        dt = ts[i+1]-ts[i]\n",
    "        p = 1-np.exp(-lambd*dt)\n",
    "        V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "        Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "        zknots = np.linspace(np.percentile(Z, 1), np.percentile(Z, 99), 20)\n",
    "        ps, fs = pwlinear_fit(Z, V, zknots)\n",
    "        q.append((ps,fs))\n",
    "        EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "        lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "\n",
    "\n",
    "\n",
    "    # independent simulation to obtain a lower bound\n",
    "    npaths = 100000\n",
    "    np.random.seed(10)\n",
    "    paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "    V = u_mat(paths[-1])\n",
    "    lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "    K = 0.5*(K_mat+K_D)\n",
    "    for i in range(len(ts)-2, -1, -1):\n",
    "        dt = ts[i+1]-ts[i]\n",
    "        p = 1-np.exp(-lambd*dt)\n",
    "        V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "        if i!=0 :\n",
    "            Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "            ps,fs = q.popleft()\n",
    "            EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "            lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "    \n",
    "    value = np.round(np.mean(V),4)\n",
    "    \n",
    "    if value>0 :\n",
    "        alpha_min = alpha\n",
    "    if value<0:\n",
    "        alpha_max = alpha\n",
    "    \n",
    "print(alpha,value)    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7LStJKRAjwBJ"
   },
   "source": [
    "A value of around 3.37 will give us 0 fair price of the contract"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "b3ShU-pnjwBJ"
   },
   "source": [
    "<b style=\"color:darkorange\">Example 4 :</b> Setting $K_D=0$, and all the other parameters remain the smae, in particualr $\\alpha=3$. Then the default event is equivalent to lapse (with no mortality payoff). It has been observed that insurance subscribers do not lapse optimally, which explains the use of the uncertain lapse model with minimum and maximum lapse rates $\\underline{\\lambda}$, $\\overline{\\lambda}$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aOrp0bnrk6ex"
   },
   "source": [
    " If insurance subscribers were to lapse optimally, i.e., exercise optimally their option to lapse, how would you price the contract? The insurance subscriber can always lapse. A rational policy holder should lapse as soon as the payoff from exercising is greater than continuing the policy. So this is similar to American option and we can price this as an american option. Here, lambda max accounts for the fact that policy holders may not lapse when they should and lambda min accounts for the fact that they may lapse when they should not. We can use lambda max = infinty and lambda min = 0 to price this as an american option."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_kNCoJW5l27-"
   },
   "source": [
    "In the uncertain lapse model we can use lambda max = infinity and lambda min = 0 to price this with optimal exercise. Here the payoff from lapse is 0 and at every time one also has to pay the fees. The values of lambda denote that the policy holder can lapse at any time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "OGKBL0emjwBK",
    "outputId": "2769befc-7345-4dcf-bbbe-0902d631bfff"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.7868522177418984"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npaths = 5000\n",
    "np.random.seed(2954)\n",
    "u_D = lambda x : 0\n",
    "lambda_max = np.inf\n",
    "lambda_min = 0\n",
    "paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "betas = np.zeros((len(ts), 2), dtype=np.float)\n",
    "V = u_mat(paths[-1])\n",
    "lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = K_mat\n",
    "q = deque()\n",
    "\n",
    "for i in range(len(ts)-2, 0, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "    zknots = np.linspace(np.percentile(Z, 1), np.percentile(Z, 99), 20)\n",
    "    ps, fs = pwlinear_fit(Z, V, zknots)\n",
    "    q.append((ps,fs))\n",
    "    EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "    lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "\n",
    "    \n",
    "    \n",
    "# independent simulation to obtain a lower bound\n",
    "npaths = 100000\n",
    "paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "V = u_mat(paths[-1])\n",
    "lambd = np.where(u_D(paths[-1]) >= V, lambda_max, lambda_min)\n",
    "K = K_mat\n",
    "for i in range(len(ts)-2, -1, -1):\n",
    "    dt = ts[i+1]-ts[i]\n",
    "    p = 1-np.exp(-lambd*dt)\n",
    "    V = p*u_D(paths[i+1]) + (1-p)*V - alpha*dt\n",
    "    if i!=0 :\n",
    "        Z = blackscholes_price(K, ts[-1]-ts[i], paths[i], vol, callput='put')\n",
    "        ps,fs = q.popleft()\n",
    "        EV    = sum([f(Z)*p for (f, p) in zip(fs, ps)])\n",
    "        lambd = np.where(u_D(paths[i]) >= EV, lambda_max, lambda_min)\n",
    "np.mean(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m4PDCvEiRXrI"
   },
   "source": [
    "Implementing Longstaff-Schwartz algorithm to price the American option with payoff $u^D=0$ (peforming a second indepedent run to get a clean lower bound price). We compare the price from uncertain lapse model with the American price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "6OpJh36kwcC1",
    "outputId": "f5260101-157a-4b74-e7ee-172f48ffc953"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.707074125017274"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(2954)\n",
    "paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "payoff = np.maximum(K_mat-paths[-1], 0)-alpha*T\n",
    "dt = 1/12\n",
    "q = deque()\n",
    "\n",
    "for i in range(len(ts)-2, 0, -1):\n",
    "    p = np.polyfit(paths[i], payoff, deg=2)\n",
    "    q.append(p)\n",
    "    contval = np.polyval(p, paths[i])\n",
    "    exerval = np.zeros(len(paths[i]))-alpha*ts[i]\n",
    "    # identify the paths where we should exercise\n",
    "    ind = exerval > contval\n",
    "    payoff[ind] = exerval[ind]\n",
    "    \n",
    "np.mean(payoff)\n",
    "    \n",
    "## Independent Monte Carlo Sumulation\n",
    "paths = blackscholes_mc(S=S, vol=vol, r=0, q=0, ts=ts, npaths=npaths)\n",
    "payoff = np.maximum(K_mat-paths[-1], 0)-alpha*T\n",
    "for i in range(len(ts)-2, 0, -1):\n",
    "    p = q.popleft()\n",
    "    contval = np.polyval(p, paths[i])\n",
    "    exerval = np.zeros(len(paths[i]))-alpha*ts[i]\n",
    "    # identify the paths where we should exercise\n",
    "    ind = exerval > contval\n",
    "    payoff[ind] = exerval[ind]\n",
    "\n",
    "np.mean(payoff)\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "ShashankNathani_sn2954_HW3_NonLinearOptionPricing.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  },
  "name": "Uncertain Mortality.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}