Least Square Monte Carlo Implementation in a Python Class.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****\n",
    "# Least-Square Monte Carlo for American Options in a Python Class\n",
    "****\n",
    "## About this notebook...\n",
    "\n",
    "\n",
    "- Documention prepared by **Jesus Perez Colino**.Version 0.4, Released 01/04/2013, Beta\n",
    "\n",
    "\n",
    "- This work is licensed under a [Creative Commons Attribution-ShareAlike 3.0 Unported License](http://creativecommons.org/licenses/by-sa/3.0/deed.en_US). This work is offered for free, with the hope that it will be useful.\n",
    "\n",
    "\n",
    "- **Summary**: This notebook was made as a documentation to support a course that I gave in Germany 2013, about Scientific Python and Quantitative Finance. This is a Python implementation of the well-know **Least-Square Monte Carlo (LSMC) method** for the valuation of **American Options** as it appears in the **Longstaff-Schwartz (2001)** paper. \n",
    "\n",
    "\n",
    "- **Reference**: [Longstaff-Schwartz (2001): \"Valuing American Options by Simulation: A Simple Least-Squares Approach.\" Review of Financial Studies, Vol. 14, 113-147](https://www.google.de/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0CCEQFjAAahUKEwiXtNSZm4rHAhXHOhQKHTjBD3k&url=https%3A%2F%2Fpeople.math.ethz.ch%2F~hjfurrer%2Fteaching%2FLongstaffSchwartzAmericanOptionsLeastSquareMonteCarlo.pdf&ei=7PO9VZeOBcf1ULiCv8gH&usg=AFQjCNFQr1r_Cf_pxylg_amU3TFOZVDc8w&sig2=ixZnX_wWQ48G66BMuQTPZA&bvm=bv.99261572,d.d24)\n",
    "\n",
    "\n",
    "- **Python & packages versions** to reproduce the results of this notebook: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------- Least-Squares MC for American Options: Conditions for Replication ---------\n",
      "Python version:     2.7.10 |Anaconda 2.3.0 (x86_64)| (default, Sep 15 2015, 14:29:08) \n",
      "[GCC 4.2.1 (Apple Inc. build 5577)]\n",
      "Numpy version:      1.9.2\n",
      "IPython version:    4.0.0\n",
      "-------------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "import IPython\n",
    "import numpy as np\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "from sys import version \n",
    "print ' Least-Squares MC for American Options: Conditions for Replication '.center(85,\"-\")\n",
    "print 'Python version:     ' + version \n",
    "print 'Numpy version:      ' + np.__version__\n",
    "print 'IPython version:    ' + IPython.__version__\n",
    "print '-'*85"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****\n",
    "# A Python Class for American Options pricing using LSMC \n",
    "****"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic idea of the **Longstaff-Schwartz** algorithm, described in detail in **Longstaff and Schwartz (2001)** (and similar approaches like those reported in **Clement and Protter (2002)**), is to use least-squares regression on a finite set of functions as a proxy for condtional expectation estimates. \n",
    "\n",
    "First of all, the time axis has to be discretized—i.e., if the American option is alive within the time horizon $[0, T ]$, early exercise is only allowed at discrete times $0 < t_1 < t_2 < ... < t_J = T$. The American option is thus approximated by a Bermudan option. For a particular exercise date $t_k$, early exercise is performed if the payoff from immediate exercise exceeds the continuation value—i.e., the value of the (remaining) option if it is not exercised at $t_k$. This continuation value can be expressed as conditional expectation of the option payoff with respect to the risk-neutral pricing measure $Q$. \n",
    "\n",
    "- The initial step of the actual algorithm is to determine the cashflow vector $C_{t_J}$ at the last timestep $t_J$ . These cashflows are easy to get because the continuation values are then zero, or in other words, it is directly the payoff of a vanilla call-option in the terminal value of each simulation $i$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{lll}\n",
    "C_{i,t_J}& =& \\max \\left(S_{i,t_J} - K,0\\right)\\\\\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "- Second, we consider the spot prices at time-step $t_{J−1}$, and estimate the **exercise value**, selecting those for which has positive payoff, or \n",
    "\n",
    "\\begin{equation} \\max \\left(S_{i,t_{J-1}}-K,0\\right) > 0 \\end{equation} \n",
    "\n",
    "\n",
    "- In order to obtain the mentioned **continuation values**, **Longstaff-Schwartz** regress the discounted future cash-flows realized from continuing onto a finite set of basis functions of our values for the spot price. The regression is done by using the values from all of the paths. The set of the basis functions for the regression, in this notebook, is a polynomial regression (of 5 degrees) but it could also be *Hermite, Legendre, Chebyshev, Gegenbauer, or Jacobi polynomials*, for example. If $S$ is the spot price, $a_j$ are coefficients and $B_j$ is the set of basis functions, then the **continuation value** for a path $i$ with values $S_{i,t_n}$ at time $t_{n}$ is\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{lll}\n",
    "Cont_{i,t_{n-1}}& =& \\sum_{j=0}^{\\infty} a_j\\left(t_{n}\\right)B_j(S)\\\\\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\n",
    "- Once we have the **continuation values** and **exercise values** , we will perform **early exercise condition of the american option** whenever\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{lll}\n",
    "C_{i,t_{n-1}}& >& Cont_{i,t_{n-1}}\\\\\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "- And finally, we then **step backward** through time, until we reach the *first time-step*. At each time-step, early exercise is performed as described previously. Note that whenever a cashflow at timestep $t_k$ is generated by early exercise in path $i$, all cashflows that occur in this path later than $t_k$ (this is, at most, one) have to be removed.\n",
    "\n",
    "\n",
    "- Once the whole backward process reach the initial point, we can build a cashflow or **value matrix** from the cashflow vectors $C_{i,t_{n}}$ by concatenating the cashflow vectors $C_{i,t_{n}}$, $n = 1,\\ldots,J$, and the option value is given by the arithmetic average of the row sums."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "class AmericanOptionsLSMC(object):\n",
    "    \"\"\" Class for American options pricing using Longstaff-Schwartz (2001):\n",
    "    \"Valuing American Options by Simulation: A Simple Least-Squares Approach.\"\n",
    "    Review of Financial Studies, Vol. 14, 113-147.\n",
    "    S0 : float : initial stock/index level\n",
    "    strike : float : strike price\n",
    "    T : float : time to maturity (in year fractions)\n",
    "    M : int : grid or granularity for time (in number of total points)\n",
    "    r : float : constant risk-free short rate\n",
    "    div :    float : dividend yield\n",
    "    sigma :  float : volatility factor in diffusion term \n",
    "    \n",
    "    Unitest(doctest): \n",
    "    >>> AmericanPUT = AmericanOptionsLSMC('put', 36., 40., 1., 50, 0.06, 0.06, 0.2, 10000  )\n",
    "    >>> AmericanPUT.price\n",
    "    4.4731177017712209\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, option_type, S0, strike, T, M, r, div, sigma, simulations):\n",
    "        try:\n",
    "            self.option_type = option_type\n",
    "            assert isinstance(option_type, str)\n",
    "            self.S0 = float(S0)\n",
    "            self.strike = float(strike)\n",
    "            assert T > 0\n",
    "            self.T = float(T)\n",
    "            assert M > 0\n",
    "            self.M = int(M)\n",
    "            assert r >= 0\n",
    "            self.r = float(r)\n",
    "            assert div >= 0\n",
    "            self.div = float(div)\n",
    "            assert sigma > 0\n",
    "            self.sigma = float(sigma)\n",
    "            assert simulations > 0\n",
    "            self.simulations = int(simulations)\n",
    "        except ValueError:\n",
    "            print('Error passing Options parameters')\n",
    "\n",
    "\n",
    "        if option_type != 'call' and option_type != 'put':\n",
    "            raise ValueError(\"Error: option type not valid. Enter 'call' or 'put'\")\n",
    "        if S0 < 0 or strike < 0 or T <= 0 or r < 0 or div < 0 or sigma < 0:\n",
    "            raise ValueError('Error: Negative inputs not allowed')\n",
    "\n",
    "        self.time_unit = self.T / float(self.M)\n",
    "        self.discount = np.exp(-self.r * self.time_unit)\n",
    "\n",
    "    @property\n",
    "    def MCprice_matrix(self, seed = 123):\n",
    "        \"\"\" Returns MC price matrix rows: time columns: price-path simulation \"\"\"\n",
    "        np.random.seed(seed)\n",
    "        MCprice_matrix = np.zeros((self.M + 1, self.simulations), dtype=np.float64)\n",
    "        MCprice_matrix[0,:] = self.S0\n",
    "        for t in xrange(1, self.M + 1):\n",
    "            brownian = np.random.standard_normal( self.simulations / 2)\n",
    "            brownian = np.concatenate((brownian, -brownian))\n",
    "            MCprice_matrix[t, :] = (MCprice_matrix[t - 1, :]\n",
    "                                  * np.exp((self.r - self.sigma ** 2 / 2.) * self.time_unit\n",
    "                                  + self.sigma * brownian * np.sqrt(self.time_unit)))\n",
    "        return MCprice_matrix\n",
    "\n",
    "    @property\n",
    "    def MCpayoff(self):\n",
    "        \"\"\"Returns the inner-value of American Option\"\"\"\n",
    "        if self.option_type == 'call':\n",
    "            payoff = np.maximum(self.MCprice_matrix - self.strike,\n",
    "                           np.zeros((self.M + 1, self.simulations),dtype=np.float64))\n",
    "        else:\n",
    "            payoff = np.maximum(self.strike - self.MCprice_matrix,\n",
    "                            np.zeros((self.M + 1, self.simulations),\n",
    "                            dtype=np.float64))\n",
    "        return payoff\n",
    "\n",
    "    @property\n",
    "    def value_vector(self):\n",
    "        value_matrix = np.zeros_like(self.MCpayoff)\n",
    "        value_matrix[-1, :] = self.MCpayoff[-1, :]\n",
    "        for t in range(self.M - 1, 0 , -1):\n",
    "            regression = np.polyfit(self.MCprice_matrix[t, :], value_matrix[t + 1, :] * self.discount, 5)\n",
    "            continuation_value = np.polyval(regression, self.MCprice_matrix[t, :])\n",
    "            value_matrix[t, :] = np.where(self.MCpayoff[t, :] > continuation_value,\n",
    "                                          self.MCpayoff[t, :],\n",
    "                                          value_matrix[t + 1, :] * self.discount)\n",
    "\n",
    "        return value_matrix[1,:] * self.discount\n",
    "\n",
    "\n",
    "    @property\n",
    "    def price(self): return np.sum(self.value_vector) / float(self.simulations)\n",
    "    \n",
    "    @property\n",
    "    def delta(self):\n",
    "        diff = self.S0 * 0.01\n",
    "        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0 + diff, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma, self.simulations)\n",
    "        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0 - diff, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma, self.simulations)\n",
    "        return (myCall_1.price - myCall_2.price) / float(2. * diff)\n",
    "    \n",
    "    @property\n",
    "    def gamma(self):\n",
    "        diff = self.S0 * 0.01\n",
    "        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0 + diff, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma, self.simulations)\n",
    "        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0 - diff, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma, self.simulations)\n",
    "        return (myCall_1.delta - myCall_2.delta) / float(2. * diff)\n",
    "    \n",
    "    @property\n",
    "    def vega(self):\n",
    "        diff = self.sigma * 0.01\n",
    "        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma + diff, \n",
    "                                       self.simulations)\n",
    "        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0,\n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma - diff, \n",
    "                                       self.simulations)\n",
    "        return (myCall_1.price - myCall_2.price) / float(2. * diff)    \n",
    "    \n",
    "    @property\n",
    "    def rho(self):        \n",
    "        diff = self.r * 0.01\n",
    "        if (self.r - diff) < 0:        \n",
    "            myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r + diff, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "            myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "            return (myCall_1.price - myCall_2.price) / float(diff)\n",
    "        else:\n",
    "            myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r + diff, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "            myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T, self.M, \n",
    "                                       self.r - diff, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "            return (myCall_1.price - myCall_2.price) / float(2. * diff)\n",
    "    \n",
    "    @property\n",
    "    def theta(self): \n",
    "        diff = 1 / 252.\n",
    "        myCall_1 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T + diff, self.M, \n",
    "                                       self.r, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "        myCall_2 = AmericanOptionsLSMC(self.option_type, self.S0, \n",
    "                                       self.strike, self.T - diff, self.M, \n",
    "                                       self.r, self.div, self.sigma, \n",
    "                                       self.simulations)\n",
    "        return (myCall_2.price - myCall_1.price) / float(2. * diff)\n",
    "    \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TestResults(failed=0, attempted=2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import doctest\n",
    "doctest.testmod()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Price:  4.47311770177\n",
      "Delta:  -0.711225132473\n",
      "Gamma:  0.126152332031\n",
      "Vega:   12.1968358245\n",
      "Rho:    -10.0335229852\n",
      "Theta:  -1.82717282672\n"
     ]
    }
   ],
   "source": [
    "AmericanPUT = AmericanOptionsLSMC('put', 36., 40., 1., 50, 0.06, 0.06, 0.2, 10000  )\n",
    "print 'Price: ', AmericanPUT.price\n",
    "print 'Delta: ', AmericanPUT.delta\n",
    "print 'Gamma: ', AmericanPUT.gamma\n",
    "print 'Vega:  ', AmericanPUT.vega\n",
    "print 'Rho:   ', AmericanPUT.rho\n",
    "print 'Theta: ', AmericanPUT.theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def prices():\n",
    "    for S0 in (36., 38., 40., 42., 44.):  # initial stock price values\n",
    "        for vol in (0.2, 0.4):  # volatility values\n",
    "            for T in (1.0, 2.0):  # times-to-maturity\n",
    "                AmericanPUT = AmericanOptionsLSMC('put', S0, 40., T, 50, 0.06, 0.06, vol, 1500)\n",
    "                print \"Initial price: %4.1f, Sigma: %4.2f, Expire: %2.1f --> Option Value %8.3f\" % (S0, vol, T, AmericanPUT.price)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial price: 36.0, Sigma: 0.20, Expire: 1.0 --> Option Value    4.439\n",
      "Initial price: 36.0, Sigma: 0.20, Expire: 2.0 --> Option Value    4.779\n",
      "Initial price: 36.0, Sigma: 0.40, Expire: 1.0 --> Option Value    7.135\n",
      "Initial price: 36.0, Sigma: 0.40, Expire: 2.0 --> Option Value    8.459\n",
      "Initial price: 38.0, Sigma: 0.20, Expire: 1.0 --> Option Value    3.225\n",
      "Initial price: 38.0, Sigma: 0.20, Expire: 2.0 --> Option Value    3.726\n",
      "Initial price: 38.0, Sigma: 0.40, Expire: 1.0 --> Option Value    6.134\n",
      "Initial price: 38.0, Sigma: 0.40, Expire: 2.0 --> Option Value    7.666\n",
      "Initial price: 40.0, Sigma: 0.20, Expire: 1.0 --> Option Value    2.296\n",
      "Initial price: 40.0, Sigma: 0.20, Expire: 2.0 --> Option Value    2.808\n",
      "Initial price: 40.0, Sigma: 0.40, Expire: 1.0 --> Option Value    5.201\n",
      "Initial price: 40.0, Sigma: 0.40, Expire: 2.0 --> Option Value    6.815\n",
      "Initial price: 42.0, Sigma: 0.20, Expire: 1.0 --> Option Value    1.589\n",
      "Initial price: 42.0, Sigma: 0.20, Expire: 2.0 --> Option Value    2.145\n",
      "Initial price: 42.0, Sigma: 0.40, Expire: 1.0 --> Option Value    4.484\n",
      "Initial price: 42.0, Sigma: 0.40, Expire: 2.0 --> Option Value    6.123\n",
      "Initial price: 44.0, Sigma: 0.20, Expire: 1.0 --> Option Value    1.088\n",
      "Initial price: 44.0, Sigma: 0.20, Expire: 2.0 --> Option Value    1.646\n",
      "Initial price: 44.0, Sigma: 0.40, Expire: 1.0 --> Option Value    3.838\n",
      "Initial price: 44.0, Sigma: 0.40, Expire: 2.0 --> Option Value    5.438\n",
      "Duration in Seconds 12.214\n"
     ]
    }
   ],
   "source": [
    "from time import time\n",
    "t0 = time()\n",
    "optionValues = prices()  # calculate all values\n",
    "t1 = time(); d1 = t1 - t0\n",
    "print \"Duration in Seconds %6.3f\" % d1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}