Asian Option valuation with MC and Control Variate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#A Python Class to price Asian Options with Control Variate\n",
    "****"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About this notebook: \n",
    "Notebook prepared by **Jesus Perez Colino** Version 0.2 (July-2015), First Released: 25/11/2012, Alpha\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 shown an example of class in Python, to estimate the price of **Arithmetic Asian option** using **Monte-Carlo Simulation**, *with* and *without* **Control Variate** over the Geometric average option value.\n",
    "\n",
    "\n",
    "- **Python & packages versions** to reproduce the results of this notebook: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================================================================================\n",
      "Python version:     2.7.10 |Anaconda 2.3.0 (x86_64)| (default, May 28 2015, 17:04:42) \n",
      "[GCC 4.2.1 (Apple Inc. build 5577)]\n",
      "Numpy version:      1.9.2\n",
      "IPython version:    3.2.0\n",
      "=====================================================================================\n"
     ]
    }
   ],
   "source": [
    "import IPython\n",
    "import numpy as np\n",
    "from sys import version \n",
    "print '='*85\n",
    "print 'Python version:     ' + version\n",
    "print 'Numpy version:      ' + np.__version__\n",
    "print 'IPython version:    ' + IPython.__version__\n",
    "print '='*85"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.special import erf\n",
    "\n",
    "\n",
    "class AsianOptionMC_MC(object):\n",
    "    \"\"\" Class for Asian options pricing using control variate\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",
    "    >>> myAsianCall = AsianOptionMC_MC('call', 4., 4., 1., 100., .03, 0, .25, 10000)\n",
    "    >>> myAsianCall.value\n",
    "    (0.26141622329842962, 0.25359138249114327, 0.26924106410571597)\n",
    "    >>> myAsianCall.value_with_control_variate\n",
    "    (0.25813771282805958, 0.25771678775128265, 0.25855863790483652)\n",
    "    \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",
    "            self.T = float(T)\n",
    "            self.M = int(M)\n",
    "            self.r = float(r)\n",
    "            self.div = float(div)\n",
    "            self.sigma = float(sigma)\n",
    "            self.simulations = int(simulations)\n",
    "        except ValueError:\n",
    "            print('Error passing Options parameters')\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.T)\n",
    "\n",
    "    @property\n",
    "    def GeometricAsianOption(self):\n",
    "        sigsqT = ((self.sigma ** 2 * self.T * (self.M + 1) * (2 * self.M + 1))\n",
    "                  / (6 * self.M * self.M))\n",
    "        muT = (0.5 * sigsqT + (self.r - 0.5 * self.sigma ** 2)\n",
    "               * self.T * (self.M + 1) / (2 * self.M))\n",
    "        d1 = ((np.log(self.S0 / self.strike) + (muT + 0.5 * sigsqT))\n",
    "              / np.sqrt(sigsqT))\n",
    "        d2 = d1 - np.sqrt(sigsqT)\n",
    "        N1 = 0.5 * (1 + erf(d1 / np.sqrt(2)))\n",
    "        N2 = 0.5 * (1 + erf(d2 / np.sqrt(2)))\n",
    "        geometric_value = self.discount * (self.S0 * np.exp(muT) * N1 - self.strike * N2)\n",
    "        return geometric_value\n",
    "\n",
    "    @property\n",
    "    def price_path(self, seed = 100):\n",
    "        np.random.seed(seed)\n",
    "        price_path = (self.S0 *\n",
    "                      np.cumprod (np.exp ((self.r - 0.5 * self.sigma ** 2) * self.time_unit +\n",
    "                                    self.sigma * np.sqrt(self.time_unit) \n",
    "                                          * np.random.randn(self.simulations, self.M)), 1))\n",
    "        return price_path\n",
    "\n",
    "    @property\n",
    "    def MCPayoff(self): \n",
    "        if self.option_type == 'call':\n",
    "            MCpayoff = self.discount \\\n",
    "                              * np.maximum(np.mean(self.price_path,1)-self.strike, 0)\n",
    "        else:\n",
    "            MCpayoff = self.discount \\\n",
    "                              * np.maximum(self.strike - np.mean(self.price_path,1), 0)\n",
    "        return MCpayoff\n",
    "    \n",
    "    @property \n",
    "    def value(self):\n",
    "        MCvalue = np.mean(self.MCPayoff)\n",
    "        MCValue_std = np.std(self.MCPayoff)\n",
    "        upper_bound = MCvalue + 1.96 * MCValue_std/np.sqrt(self.simulations)\n",
    "        lower_bound = MCvalue - 1.96 * MCValue_std/np.sqrt(self.simulations)\n",
    "        return MCvalue, lower_bound, upper_bound\n",
    "    \n",
    "    @property\n",
    "    def value_with_control_variate(self):\n",
    "        \n",
    "        geometric_average = np.exp( (1/float(self.M)) * np.sum(np.log(self.price_path), 1) )\n",
    "        if self.option_type == 'call':\n",
    "            MCpayoff_geometric = self.discount * np.maximum(geometric_average - self.strike, 0)\n",
    "        else:\n",
    "            MCpayoff_geometric = self.discount * np.maximum(self.strike - geometric_average, 0)\n",
    "        value_with_CV = self.MCPayoff + self.GeometricAsianOption - MCpayoff_geometric        \n",
    "        value_with_control_variate = np.mean(value_with_CV , 0)\n",
    "        value_with_control_variate_std = np.std(value_with_CV, 0)\n",
    "        upper_bound_CV = value_with_control_variate + 1.96 * value_with_control_variate_std/np.sqrt(self.simulations)\n",
    "        lower_bound_CV = value_with_control_variate - 1.96 * value_with_control_variate_std/np.sqrt(self.simulations)        \n",
    "        return value_with_control_variate, lower_bound_CV, upper_bound_CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TestResults(failed=0, attempted=3)"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import doctest\n",
    "doctest.testmod()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myAsianCall = AsianOptionMC_MC('call', 4., 4., 1., 100., .03, 0, .25, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.26141622329842962, 0.25359138249114327, 0.26924106410571597)"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myAsianCall.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.24604398839697852"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myAsianCall.GeometricAsianOption"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.25813771282805958, 0.25771678775128265, 0.25855863790483652)"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myAsianCall.value_with_control_variate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}