files/Friday Workshop/2_3_20/option_pricing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/israeldi/quantlab/master/assets/images/Program-Logo.png\" width=\"400px\" align=\"right\">\n",
    "\n",
    "# QuantLab: Stochastics\n",
    "### [(Go to Quant Lab)](https://israeldi.github.io/quantlab/)\n",
    "\n",
    "&copy; Dr. Yves J. Hilpisch | The Python Quants GmbH\n",
    "\n",
    "<img src=\"http://hilpisch.com/images/py4fi_2nd_shadow.png\" width=\"200px\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially import all the modules we will be using for our notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import numpy.random as npr \n",
    "\n",
    "# COMMAND PROMPT: pip install matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Pricing European Options by Monte Carlo\n",
    "\n",
    "The payoff of a European call option on an index at maturity is given by $h(S_T)\\equiv max\\{S_T-K, 0\\}$, where $S_T$ is the index level at maturity date $T$ and $K$ is the strike price. Given a risk-neutral measure for the relevant stochastic process (e.g., geometric Brownian motion), the price of such an option is given by the formula:\n",
    "\n",
    "$$C_{0}=e^{-rT}\\mathbb{E}_{0}^{Q}[h(S_{T})]=e^{-rT}\\intop_{0}^{\\infty}h(s)q(s)ds$$\n",
    "\n",
    "The equation below provides the respective Monte Carlo estimator for the European option, where $\\tilde{S}_{T}^{i}$ is the $T$-th simulated index level at maturity.\n",
    "\n",
    "$$\\tilde{C}_{0}=e^{-rT}\\frac{1}{I}\\sum_{i=1}^{I}h(\\tilde{S}_{T}^{i})$$\n",
    "\n",
    "Consider the following parameterization for the geometric Brownian motion and the valuation function `gbm_mcs_stat()`, taking as a parameter only the strike price. Here, only the index level at maturity is simulated. As a reference, consider the case with a strike price of $K = 105$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bsm_call_value(S0, K, T, r, sigma):\n",
    "    ''' Valuation of European call option in BSM model.\n",
    "    Analytical formula.\n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    S0: float\n",
    "        initial stock/index level\n",
    "    K: float\n",
    "        strike price\n",
    "    T: float\n",
    "        maturity date (in year fractions)\n",
    "    r: float\n",
    "        constant risk-free short rate\n",
    "    sigma: float\n",
    "        volatility factor in diffusion term\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    value: float\n",
    "        present value of the European call option\n",
    "    '''\n",
    "    from math import log, sqrt, exp\n",
    "    from scipy import stats\n",
    "\n",
    "    S0 = float(S0)\n",
    "    d1 = (log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))\n",
    "    d2 = (log(S0 / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))\n",
    "    # stats.norm.cdf --> cumulative distribution function\n",
    "    #                    for normal distribution\n",
    "    value = (S0 * stats.norm.cdf(d1, 0.0, 1.0) -\n",
    "             K * exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))\n",
    "    return value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_sn(M, I):\n",
    "    ''' Function to generate random numbers for simulation.\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    M: int\n",
    "        number of time intervals for discretization\n",
    "    I: int\n",
    "        number of paths to be simulated\n",
    "    '''\n",
    "    sn = npr.standard_normal((M + 1, I))\n",
    "    return sn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "S0 = 100.\n",
    "r = 0.05\n",
    "sigma = 0.25\n",
    "T = 1.0\n",
    "I = 50000\n",
    "K = 105"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "uuid": "693f44be-b3dd-4820-9610-a127f0e9b31b"
   },
   "outputs": [],
   "source": [
    "def gbm_mcs_stat(K):\n",
    "    ''' Valuation of European call option in Black-Scholes-Merton\n",
    "    by Monte Carlo simulation (of index level at maturity)\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    K: float\n",
    "        (positive) strike price of the option\n",
    "    \n",
    "    Returns\n",
    "    =======\n",
    "    C0: float\n",
    "        estimated present value of European call option\n",
    "    '''\n",
    "    sn = gen_sn(1, I)\n",
    "    # simulate index level at maturity\n",
    "    ST = S0 * np.exp((r - 0.5 * sigma ** 2) * T \n",
    "                 + sigma * math.sqrt(T) * sn[1])\n",
    "    # calculate payoff at maturity\n",
    "    hT = np.maximum(ST - K, 0)\n",
    "    # calculate MCS estimator\n",
    "    C0 = math.exp(-r * T) * np.mean(hT)\n",
    "    return C0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "uuid": "f325da52-3e45-4e9e-a4a2-067efb1c3bb7"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.962393419194543"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The Monte Carlo estimator value for the European call option.\n",
    "gbm_mcs_stat(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.00220211715488"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The analytical solution is\n",
    "bsm_call_value(S0, K, T, r, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, consider the dynamic simulation approach and allow for European put options in addition to the call option. The function `gbm_mcs_dyna()` implements the algorithm. The code also compares option price estimates for a call and a put stroke at the same level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The number of time intervals for the discretization.\n",
    "M = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "uuid": "511974d5-5ceb-4b68-bf7f-e01eaa43f7c6"
   },
   "outputs": [],
   "source": [
    "def gbm_mcs_dyna(K, option='call'):\n",
    "    ''' Valuation of European options in Black-Scholes-Merton\n",
    "    by Monte Carlo simulation (of index level paths)\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    K: float\n",
    "        (positive) strike price of the option\n",
    "    option : string\n",
    "        type of the option to be valued ('call', 'put')\n",
    "    \n",
    "    Returns\n",
    "    =======\n",
    "    C0: float\n",
    "        estimated present value of European call option\n",
    "    '''\n",
    "    dt = T / M\n",
    "    \n",
    "    # simulation of index level paths\n",
    "    S = np.zeros((M + 1, I))\n",
    "    S[0] = S0\n",
    "    sn = gen_sn(M, I)\n",
    "    \n",
    "    for t in range(1, M + 1):\n",
    "        S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt \n",
    "                + sigma * math.sqrt(dt) * sn[t])\n",
    "    \n",
    "    # case-based calculation of payoff\n",
    "    if option == 'call':\n",
    "        hT = np.maximum(S[-1] - K, 0)\n",
    "    else:\n",
    "        # Put option\n",
    "        hT = np.maximum(K - S[-1], 0)\n",
    "    # calculation of MCS estimator\n",
    "    C0 = math.exp(-r * T) * np.mean(hT)\n",
    "    return C0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "uuid": "44ae2961-ec7c-4e69-b6ff-17b8093a894b"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.088345691009375"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The Monte Carlo estimator value for the European call option (Dynamic Simulation Approach)\n",
    "gbm_mcs_dyna(K, option='call')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "uuid": "bedb79ae-4f01-41ea-b16a-22ea9781fc0e"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.964122862849353"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The Monte Carlo estimator value for the European put option (Dynamic Simulation Approach)\n",
    "gbm_mcs_dyna(K, option='put')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Euler Discretization\n",
    "\n",
    "Consider the stochastic differential equation $\\mathrm{d} X_{t}=a\\left(X_{t}\\right) \\mathrm{d} t+b\\left(X_{t}\\right) \\mathrm{d} W_{t}$ with initial condition $X_{0}=x_{0},$ where $W_{t}$ stands for the Wiener process, and suppose that we wish to solve this SDE on some internal of time $[0, T]$. Then the Euler-Maruyama approximation to the true solution $X$ is the Markov chain $Y$ defined as follows:\n",
    "\n",
    "- partition the interval $[0,7 \\text { into } N \\text { equal subintervals of width } \\Delta t>0$ :\n",
    "$$\n",
    "0=\\tau_{0}<\\tau_{1}<\\cdots<\\tau_{N}=T \\text { and } \\Delta t=T / N\n",
    "$$\n",
    "- set $Y_{0}=x_{0}$\n",
    "- recursively define $Y_{n}$ for $1 \\leq n \\leq N$ by\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "{Y_{n+1}=Y_{n}+a\\left(Y_{n}\\right) \\Delta t+b\\left(Y_{n}\\right) \\Delta W_{n}} \\\\\n",
    "{\\Delta W_{n}=W_{\\tau_{n+1}}-W_{\\tau_{n}}}\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "S0 = 100.\n",
    "r = 0.05\n",
    "sigma = 0.25\n",
    "T = 1.0\n",
    "I = 50000\n",
    "K = 105"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_sims = 5000  # Display five runs\n",
    "\n",
    "N = 100  # Compute 100 grid points\n",
    "dt = T / N\n",
    "y_init = S0\n",
    "\n",
    "def dW(delta_t):\n",
    "    \"\"\"Sample a random number at each call.\"\"\"\n",
    "    return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))\n",
    "\n",
    "# Note by default all all values are initialized to zero\n",
    "ys = np.zeros((N, num_sims))\n",
    "ys[0,:] = y_init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "for j in range(num_sims):\n",
    "    for i in range(1, N):\n",
    "        y = ys[i-1, j]\n",
    "        ys[i, j] = y + (r * y * dt) + (sigma * y* dW(dt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.012785056220235"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hT = np.maximum(ys[-1,:] - K, 0)\n",
    "C0 = math.exp(-r * T) * np.mean(hT)\n",
    "C0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite Difference Scheme\n",
    "(Next Time!)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}