py4fi_L7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 7: Option Pricing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as npr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Binomial tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomialtree(T,N,S,r,sigma,K, typeEA, typeCP):\n",
    "\n",
    "    if typeCP == \"call\": # Call or put option\n",
    "        typeCP = 1\n",
    "    else:\n",
    "        typeCP = -1\n",
    "\n",
    "    dt = np.float(T)/N # step size\n",
    "    u = np.exp(sigma*np.sqrt(dt)) # up and down steps\n",
    "    d = 1/u\n",
    "    p = (np.exp(r*dt)-d)/(u-d) # risk-neutral probabilities\n",
    "\n",
    "    ST = np.zeros(N+1) # Final stock prices vector\n",
    "    option = np.zeros(N+1) # option price vector\n",
    "\n",
    "    for i in range(N+1):\n",
    "        ST[i] = S*u**(N-i)*d**i # Fill in terminal stock prices\n",
    "        option[i] = max(typeCP*(ST[i]-K),0) # Fill in terminal option prices\n",
    "    \n",
    "    for i in range(N-1,-1,-1):\n",
    "        for j in range(i+1):\n",
    "            option[j] = np.exp(-r*dt)*(p*option[j]+(1-p)*option[j+1])\n",
    "            if typeEA == \"American\":\n",
    "                ST[j] = np.exp(-r*dt)*(p*ST[j]+(1-p)*ST[j+1])\n",
    "                option[j] = max(option[j], max(typeCP*(ST[j]-K),0))\n",
    "\n",
    "    return option[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8932270788587007"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binomialtree(2, 500, 100, 0.05, 0.10, 100, \"Euro\", \"put\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.8754584823327676"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binomialtree(2, 500, 100, 0.05, 0.10, 100, \"American\", \"put\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo simulation (static, plain vanilla call & put)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montecarlo(T,N,S,r,sigma,K,typeCP):\n",
    "\n",
    "    typeCP = 1 if typeCP == \"call\" else -1\n",
    "    z = npr.randn(N)\n",
    "    # simulate index level at maturity\n",
    "    ST = S*np.exp((r-0.5*sigma**2)*T + sigma*np.sqrt(T)*z)\n",
    "    # calculate payoff at maturity\n",
    "    option = np.maximum(typeCP*(ST - K), 0)\n",
    "    option = np.exp(-r*T)*option.mean()\n",
    "    return option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8983150606145902"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "montecarlo(2, 10**7, 100, 0.05, 0.10, 100, \"put\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo simulation (dynamic, i.e., Asian options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montecarlo_dynamic(T,N,steps,S,r,sigma,K,typeCP):\n",
    "\n",
    "    typeCP = 1 if typeCP == \"call\" else -1\n",
    "    dt = np.float(T)/steps\n",
    "    SPath = np.zeros((steps+1, N)) # matrix of N paths\n",
    "    SPath[0] = S # first value is S for all paths\n",
    "    m = np.exp((r-0.5*sigma**2)*dt)\n",
    "    for t in range(1, steps+1):\n",
    "        z = npr.randn(N) # random noise for all paths at time t\n",
    "        # SPath[t]=SPath[t-1]*np.exp((r-0.5*sigma**2)*dt + sigma*np.sqrt(dt)*z)\n",
    "        SPath[t] = SPath[t-1]*m*np.exp(sigma*np.sqrt(dt)*z)\n",
    "    # calculate payoff at maturity\n",
    "    option = np.maximum(typeCP*(SPath.mean(axis=0)-K), 0)\n",
    "    option = np.exp(-r*T)*option.mean()\n",
    "    return option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 839 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.253980264129479"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time montecarlo_dynamic(2, 10**4, 1000, 100, 0.05, 0.10, 100, \"put\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least-Squares Monte Carlo (for American options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LSM_montecarlo(T,N,steps,S,r,sigma,K, typeCP):\n",
    "\n",
    "    typeCP = 1 if typeCP == \"call\" else -1    \n",
    "    dt = np.float(T) / steps\n",
    "    df = np.exp(-r * dt)\n",
    "    \n",
    "    # simulation of index levels\n",
    "    SPath = np.zeros((steps + 1, N))\n",
    "    SPath[0] = S\n",
    "    z = npr.randn(steps, N)\n",
    "    for t in range(1, steps + 1):\n",
    "        SPath[t] = SPath[t-1]*np.exp((r - 0.5*sigma**2)*dt + sigma*np.sqrt(dt)*z[t-1])\n",
    "\n",
    "    # case-based calculation of payoff\n",
    "    h = np.maximum(typeCP*(SPath - K), 0)\n",
    "    # LSM algorithm\n",
    "    V = np.copy(h)\n",
    "    for t in range(steps-1, 0, -1):\n",
    "        reg = np.polyfit(SPath[t], V[t+1]*df, 7)\n",
    "        C = np.polyval(reg, SPath[t])\n",
    "        V[t] = np.where(C > h[t], V[t+1]*df, h[t])\n",
    "    # MCS estimator\n",
    "    option = df * 1 / N * np.sum(V[1])\n",
    "    return option"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Old exam applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(2015)** A *lookback* call option pays the difference between the spot price at expiration and the minimum spot price over the term of the contract, if this difference is larger than zero. At maturity, the value of the option is:\n",
    "\n",
    "\\begin{equation}\n",
    "\\textit{Lookback}_T=\\max \\left(0, S_T-\\min_{0\\leq t \\leq T} S_t\\right).\n",
    "\\end{equation}\n",
    "\n",
    "Write a Python function to value the *lookback* call option using a **binomial tree**. \n",
    "\n",
    "1. The function takes as inputs: the stock price at $t=0$, i.e., $S_0$, the stock volatility $\\sigma$, the number of steps in the tree $N$, the risk-free rate $r$ (yearly, continuously compounded) and the maturity $T$ (in years). \n",
    "2. Use the Cox-Ross-Rubinstein model to compute the up and down steps, i.e., $u=\\exp\\left(\\sigma \\sqrt{\\frac{T}{N}}\\right)$ and $ud=1$. \n",
    "3. The output of the function is the price of the \\textit{lookback} option. \n",
    "4. Numerical example: S0=50, $\\sigma$=0.25, N=10, r=0.02, T=1.\n",
    "5. Can you value this option via Monte Carlo simulation of the final price $S_T$? Discuss. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lookback(T,N,S0,r,sigma):\n",
    "\n",
    "    ### Preliminary steps\n",
    "    dt=np.float(T)/N # step size\n",
    "    u=np.exp(sigma*np.sqrt(dt)) # up and down steps\n",
    "    d=1/u\n",
    "    p=(np.exp(r*dt)-d)/(u-d) # risk-neutral probability\n",
    "    \n",
    "    # Path generator (2**N paths, N steps)\n",
    "    Paths=np.zeros([N+1,2**N]) # initialize paths\n",
    "    Paths[0,:]=S0 # first element on path is S0\n",
    "\n",
    "    for j in range(N,0,-1): # fill in all (2^N) possible paths in a N x 2^N matrix\n",
    "        temp=2**(N-j)*([1 for i in range(0,2**(j-1))],[0 for i in range(0,2**(j-1))])\n",
    "        print(temp, \"\\n\")\n",
    "        Paths[N-j+1,:]=[val for s in temp for val in s]\n",
    "        print(Paths, \"\\n\")\n",
    "    \n",
    "    print(Paths[1:])\n",
    "    # risk-neutral probabilities\n",
    "    RN_Probs=p**Paths[1:].sum(axis=0)*(1-p)**(N-Paths[1:].sum(axis=0))\n",
    "\n",
    "    # fill the paths matrix with up and down steps\n",
    "    Paths[1:]=Paths[1:]*(u-d)+d\n",
    "\n",
    "    for i in range(1,np.size(Paths,axis=0)):\n",
    "        Paths[i,:]=Paths[i-1,:]*Paths[i,:]\n",
    "\n",
    "    ### Option price\n",
    "    option_path_values=Paths[-1,:]-Paths.min(axis=0) # option value for each path\n",
    "\n",
    "    option=np.dot(option_path_values,RN_Probs) # scalar product of value and path probability\n",
    "\n",
    "    return option_path_values, option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "([1, 1], [0, 0]) \n",
      "\n",
      "[[100. 100. 100. 100.]\n",
      " [  1.   1.   0.   0.]\n",
      " [  0.   0.   0.   0.]] \n",
      "\n",
      "([1], [0], [1], [0]) \n",
      "\n",
      "[[100. 100. 100. 100.]\n",
      " [  1.   1.   0.   0.]\n",
      " [  1.   0.   1.   0.]] \n",
      "\n",
      "[[1. 1. 0. 0.]\n",
      " [1. 0. 1. 0.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([4.24119019e+01, 1.42108547e-14, 1.62033114e+01, 0.00000000e+00]),\n",
       " 13.406035970903066)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lookback(1,2,100,0.01,0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(2016)** A *power* call option pays a (positive) power of the difference between the spot price at expiration and a strike price, if this difference is larger than zero. At maturity, the value of the option is:\n",
    "\\begin{equation}\n",
    "\\textit{Power}_{T,k}=\\max \\left(0, \\left[S_T-K\\right]^{k}\\right).\n",
    "\\end{equation}\n",
    "\n",
    "1. Write a Python function to value the power call option using a binomial tree. \n",
    "2. The function takes as inputs: the stock price at $t=0$, i.e., $S_0$, the stock volatility $\\sigma$, the number of steps in the tree $N$, the strike price $K$, the power coefficient $k$, the risk-free rate $r$ (yearly, continuously compounded) and the maturity $T$ (in years). \n",
    "3. Use the Cox-Ross-Rubinstein model to compute the up and down steps, i.e., $u=\\exp\\left(\\sigma \\sqrt{\\frac{T}{N}}\\right)$ and $ud=1$. \n",
    "4. The output of the function is the price of the lookback option. \n",
    "5. Numerical example: S0=50, sigma=0.25, N=10, r=0.02, T=1, K=50, k=2\n",
    "6. How would you evaluate a plain vanilla call option using the function you just created? \n",
    "7. Can you evaluate the power option using Monte Carlo simulations? Explain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power(T,N,S0,K,k,r,sigma):\n",
    "    \n",
    "    ### Preliminary steps\n",
    "    dt=np.float(T)/N # step size\n",
    "    u=np.exp(sigma*np.sqrt(dt)) # up and down steps\n",
    "    d=1/u\n",
    "    p=(np.exp(r*dt)-d)/(u-d) # risk-neutral probability\n",
    "\n",
    "    #### Path generator\n",
    "    Paths=np.zeros([N+1,2**N]) # initialize paths\n",
    "    Paths[0,:]=S0 # first element on path is S0\n",
    "\n",
    "    for j in range(N,0,-1): # fill in all (2^N) possible paths in a N x 2^N matrix\n",
    "        temp=2**(N-j)*([1 for i in range(0,2**(j-1))],[0 for i in range(0,2**(j-1))])\n",
    "        Paths[N-j+1,:]=[val for s in temp for val in s]\n",
    "\n",
    "    # risk-neutral probabilities\n",
    "    RN_Probs=p**Paths[1:].sum(axis=0)*(1-p)**(N-Paths[1:].sum(axis=0))\n",
    "\n",
    "    # fill the paths matrix with up and down steps\n",
    "    Paths[1:]=Paths[1:]*(u-d)+d\n",
    "\n",
    "    for i in range(1,np.size(Paths,axis=0)):\n",
    "        Paths[i,:]=Paths[i-1,:]*Paths[i,:]\n",
    "\n",
    "    ### Option price\n",
    "    option_path_values=np.maximum(0,(Paths[-1,:]-K)**k) # option value for each path\n",
    "\n",
    "    option=np.exp(-r*T)*np.dot(option_path_values,RN_Probs) # scalar product of value and path probability\n",
    "\n",
    "    return option_path_values, Paths, option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power2(T,N,S,K,k,r,sigma):\n",
    "\n",
    "    dt=np.float(T)/N # step size\n",
    "    u=np.exp(sigma*np.sqrt(dt)) # up and down steps\n",
    "    d=1/u\n",
    "    p=(np.exp(r*dt)-d)/(u-d) # risk-neutral probabilities\n",
    "\n",
    "    ST=np.zeros(N+1) # Final stock prices vector\n",
    "    option=np.zeros(N+1) # option price vector\n",
    "\n",
    "    for i in range(0,N+1):\n",
    "        ST[i]=S*u**(N-i)*d**i # Fill in terminal stock prices\n",
    "        option[i]=max((ST[i]-K)**k,0) # Fill in terminal option prices\n",
    "    \n",
    "    for i in range(N-1,-1,-1):\n",
    "        for j in range(0,i+1):\n",
    "            option[j]=np.exp(-r*dt)*(p*option[j]+(1-p)*option[j+1])\n",
    "\n",
    "    return option[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(784.3374468302424, 784.337446830242)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power(1,5,100,150,3,0.10,0.25)[2], power2(1,5,100,150,3,0.10,0.25) # two ways to solve the same problem"
   ]
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}