monte_carlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b2dd0f2d",
   "metadata": {},
   "source": [
    "\n",
    "<a id='monte-carlo'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45ff922f",
   "metadata": {},
   "source": [
    "# Monte Carlo and Option Pricing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c44907ad",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Simple probability calculations can be done either\n",
    "\n",
    "- with pencil and paper, or  \n",
    "- by looking up facts about well known probability distributions, or  \n",
    "- in our heads.  \n",
    "\n",
    "\n",
    "For example, we can easily work out\n",
    "\n",
    "- the probability of three heads in five flips of a fair coin  \n",
    "- the expected value of a random variable that equals $ -10 $ with probability\n",
    "  $ 1/2 $ and $ 100 $ with probability $ 1/2 $.  \n",
    "\n",
    "\n",
    "But some probability calculations are very complex.\n",
    "\n",
    "Complex calculations concerning probabilities and expectations occur in many\n",
    "economic and financial problems.\n",
    "\n",
    "Perhaps the most important tool for handling complicated probability\n",
    "calculations is [Monte Carlo methods](https://en.wikipedia.org/wiki/Monte_Carlo_method).\n",
    "\n",
    "In this lecture we introduce Monte Carlo methods for computing expectations,\n",
    "with some applications in finance.\n",
    "\n",
    "We will use the following imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370df99b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import randn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e0248c3",
   "metadata": {},
   "source": [
    "## An introduction to Monte Carlo\n",
    "\n",
    "In this section we describe how Monte Carlo can be used to compute\n",
    "expectations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42260d50",
   "metadata": {},
   "source": [
    "### Share price with known distribution\n",
    "\n",
    "Suppose that we are considering buying a share in some company.\n",
    "\n",
    "Our plan is either to\n",
    "\n",
    "1. buy the share now, hold it for one year and then sell it, or  \n",
    "1. do something else with our money.  \n",
    "\n",
    "\n",
    "We start by thinking of the share price in one year as a random variable $ S $.\n",
    "\n",
    "Before deciding whether or not to buy the share, we need to know some features\n",
    "of the distribution of $ S $.\n",
    "\n",
    "For example, suppose the mean of $ S $ is high relative to the price of buying\n",
    "the share.\n",
    "\n",
    "This suggests we have a good chance of selling at a relatively high price.\n",
    "\n",
    "Suppose, however, that the variance of $ S $ is also high.\n",
    "\n",
    "This suggests that buying the share is risky, so perhaps we should refrain.\n",
    "\n",
    "Either way, this discussion shows the importance of understanding the\n",
    "distribution of $ S $.\n",
    "\n",
    "Suppose that, after analyzing the data, we guess that $ S $ is well\n",
    "represented by a lognormal distribution with parameters $ \\mu, \\sigma $ .\n",
    "\n",
    "- $ S $ has the same distribution as $ \\exp(\\mu + \\sigma Z) $ where $ Z $ is standard normal.  \n",
    "- We write this statement as $ S \\sim LN(\\mu, \\sigma) $.  \n",
    "\n",
    "\n",
    "Any good reference on statistics (such as\n",
    "[Wikipedia](https://en.wikipedia.org/wiki/Log-normal_distribution)) will tell\n",
    "us that the mean and variance are\n",
    "\n",
    "$$\n",
    "\\mathbb E S\n",
    "        = \\exp \\left(\\mu + \\frac{\\sigma^2}{2} \\right)\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\mathop{\\mathrm{Var}} S\n",
    "    = [\\exp(\\sigma^2) - 1] \\exp(2\\mu + \\sigma^2)\n",
    "$$\n",
    "\n",
    "So far we have no need for a computer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5128a52d",
   "metadata": {},
   "source": [
    "### Share price with unknown distribution\n",
    "\n",
    "But now suppose that we study the distribution of $ S $ more carefully.\n",
    "\n",
    "We decide that the share price depends on three variables, $ X_1 $, $ X_2 $, and\n",
    "$ X_3 $ (e.g., sales, inflation, and interest rates).\n",
    "\n",
    "In particular, our study suggests that\n",
    "\n",
    "$$\n",
    "S = (X_1 + X_2 + X_3)^p\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ p $ is a positive number, which is known to us  (i.e., has been estimated),  \n",
    "- $ X_i \\sim LN(\\mu_i, \\sigma_i) $ for $ i=1,2,3 $,  \n",
    "- the values $ \\mu_i, \\sigma_i $ are also known, and  \n",
    "- the random variables $ X_1 $, $ X_2 $ and $ X_3 $ are independent.  \n",
    "\n",
    "\n",
    "How should we compute the mean of $ S $?\n",
    "\n",
    "To do this with pencil and paper is hard (unless, say, $ p=1 $).\n",
    "\n",
    "But fortunately there’s an easy way to do this, at least approximately.\n",
    "\n",
    "This is the Monte Carlo method, which runs as follows:\n",
    "\n",
    "1. Generate $ n $ independent draws of $ X_1 $, $ X_2 $ and $ X_3 $ on a computer,  \n",
    "1. use these draws to generate $ n $ independent draws of $ S $, and  \n",
    "1. take the average value of these draws of $ S $.  \n",
    "\n",
    "\n",
    "This average will be close to the true mean when $ n $ is large.\n",
    "\n",
    "This is due to the law of large numbers, which we discussed in [LLN and CLT](https://intro.quantecon.org/lln_clt.html).\n",
    "\n",
    "We use the following values for $ p $ and each $ \\mu_i $ and $ \\sigma_i $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f2626a8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 1_000_000\n",
    "p = 0.5\n",
    "μ_1, μ_2, μ_3 = 0.2, 0.8, 0.4\n",
    "σ_1, σ_2, σ_3 = 0.1, 0.05, 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b535de03",
   "metadata": {},
   "source": [
    "#### A routine using loops in python\n",
    "\n",
    "Here’s a routine using native Python loops to calculate the desired mean\n",
    "\n",
    "$$\n",
    "\\frac{1}{n} \\sum_{i=1}^n S_i\n",
    "    \\approx \\mathbb E S\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05891f1e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "S = 0.0\n",
    "for i in range(n):\n",
    "    X_1 = np.exp(μ_1 + σ_1 * randn())\n",
    "    X_2 = np.exp(μ_2 + σ_2 * randn())\n",
    "    X_3 = np.exp(μ_3 + σ_3 * randn())\n",
    "    S += (X_1 + X_2 + X_3)**p\n",
    "S / n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97368403",
   "metadata": {},
   "source": [
    "We can also construct a function that contains these operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32649205",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_mean(n=1_000_000):\n",
    "    S = 0.0\n",
    "    for i in range(n):\n",
    "        X_1 = np.exp(μ_1 + σ_1 * randn())\n",
    "        X_2 = np.exp(μ_2 + σ_2 * randn())\n",
    "        X_3 = np.exp(μ_3 + σ_3 * randn())\n",
    "        S += (X_1 + X_2 + X_3)**p\n",
    "    return (S / n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cced911",
   "metadata": {},
   "source": [
    "Now let’s call it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51d256e4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "compute_mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22fcebd1",
   "metadata": {},
   "source": [
    "### A vectorized routine\n",
    "\n",
    "If we want a more accurate estimate we should increase $ n $.\n",
    "\n",
    "But the code above runs quite slowly.\n",
    "\n",
    "To make it faster, let’s implement a vectorized routine using NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a58d5032",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_mean_vectorized(n=1_000_000):\n",
    "    X_1 = np.exp(μ_1 + σ_1 * randn(n))\n",
    "    X_2 = np.exp(μ_2 + σ_2 * randn(n))\n",
    "    X_3 = np.exp(μ_3 + σ_3 * randn(n))\n",
    "    S = (X_1 + X_2 + X_3)**p\n",
    "    return S.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce672e0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "compute_mean_vectorized()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91712ef8",
   "metadata": {},
   "source": [
    "Notice that this routine is much faster.\n",
    "\n",
    "We can increase $ n $ to get more accuracy and still have reasonable speed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15f24842",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "compute_mean_vectorized(n=10_000_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43eb49a4",
   "metadata": {},
   "source": [
    "## Pricing a European call option under risk neutrality\n",
    "\n",
    "Next we are going to price a European call option under risk neutrality.\n",
    "\n",
    "Let’s first discuss risk neutrality and then consider European options."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36f43d7a",
   "metadata": {},
   "source": [
    "### Risk-neutral pricing\n",
    "\n",
    "When we use risk-neutral pricing, we determine the price of a given asset\n",
    "according to its expected payoff:\n",
    "\n",
    "$$\n",
    "\\text{cost } = \\text{ expected benefit}\n",
    "$$\n",
    "\n",
    "For example, suppose someone promises to pay you\n",
    "\n",
    "- 1,000,000 dollars if “heads” is the outcome of a fair coin flip  \n",
    "- 0 dollars if “tails” is the outcome  \n",
    "\n",
    "\n",
    "Let’s denote the payoff as $ G $, so that\n",
    "\n",
    "$$\n",
    "\\mathbb P\\left\\{G = 10^6 \\right\\} = \\mathbb P\\{G = 0\\} = \\frac{1}{2}\n",
    "$$\n",
    "\n",
    "Suppose in addition that you can sell this promise to anyone who wants it.\n",
    "\n",
    "- First they pay you $ P $, the price at which you sell it  \n",
    "- Then they get $ G $, which could be either 1,000,000 or 0.  \n",
    "\n",
    "\n",
    "What’s a fair price for this asset (this promise)?\n",
    "\n",
    "The definition of “fair” is ambiguous, but we can say that the\n",
    "**risk-neutral price** is 500,000 dollars.\n",
    "\n",
    "This is because the risk-neutral price is just the expected payoff of the\n",
    "asset, which is\n",
    "\n",
    "$$\n",
    "\\mathbb E G = \\frac{1}{2} \\times 10^6 + \\frac{1}{2} \\times 0 = 5 \\times 10^5\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc366806",
   "metadata": {},
   "source": [
    "### A comment on risk\n",
    "\n",
    "As suggested by the name, the risk-neutral price ignores risk.\n",
    "\n",
    "To understand this, consider whether you would pay 500,000 dollars for such a\n",
    "promise.\n",
    "\n",
    "Would you prefer to receive 500,000 for sure or 1,000,000 dollars with\n",
    "50% probability and nothing with 50% probability?\n",
    "\n",
    "At least some readers will strictly prefer the first option — although some\n",
    "might prefer the second.\n",
    "\n",
    "Thinking about this makes us realize that 500,000 is not necessarily the\n",
    "“right” price — or the price that we would see if there was a market for\n",
    "these promises.\n",
    "\n",
    "Nonetheless, the risk-neutral price is an important benchmark, which economists\n",
    "and financial market participants try to calculate every day."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9435daa",
   "metadata": {},
   "source": [
    "### Discounting\n",
    "\n",
    "Another thing we ignored in the previous discussion was time.\n",
    "\n",
    "In general, receiving $ x $ dollars now is preferable to receiving $ x $ dollars\n",
    "in $ n $ periods (e.g., 10 years).\n",
    "\n",
    "After all, if we receive $ x $ dollars now, we could put it in the bank at\n",
    "interest rate $ r > 0 $ and receive $ (1 + r)^n x $ in $ n $ periods.\n",
    "\n",
    "Hence future payments need to be discounted when we consider their present\n",
    "value.\n",
    "\n",
    "We will implement discounting by\n",
    "\n",
    "- multiplying a payment in one period by $ \\beta < 1 $  \n",
    "- multiplying a payment in $ n $ periods by $ \\beta^n $, etc.  \n",
    "\n",
    "\n",
    "The same adjustment needs to be applied to our risk-neutral price for the\n",
    "promise described above.\n",
    "\n",
    "Thus, if $ G $ is realized in $ n $ periods, then the risk-neutral price is\n",
    "\n",
    "$$\n",
    "P = \\beta^n \\mathbb E G\n",
    "      = \\beta^n 5 \\times 10^5\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdfa48e8",
   "metadata": {},
   "source": [
    "### European call options\n",
    "\n",
    "Now let’s price a European call option.\n",
    "\n",
    "The option is described by three things:\n",
    "\n",
    "1. $ n $, the **expiry date**,  \n",
    "1. $ K $, the **strike price**, and  \n",
    "1. $ S_n $, the price of the **underlying** asset at date $ n $.  \n",
    "\n",
    "\n",
    "For example, suppose that the underlying is one share in Amazon.\n",
    "\n",
    "The owner of this option has the right to buy one share in Amazon at price $ K $ after $ n $ days.\n",
    "\n",
    "If $ S_n > K $, then the owner will exercise the option, buy at $ K $, sell at\n",
    "$ S_n $, and make profit $ S_n - K $.\n",
    "\n",
    "If $ S_n \\leq K $, then the owner will not exercise the option and the payoff is zero.\n",
    "\n",
    "Thus, the payoff is $ \\max\\{ S_n - K, 0 \\} $.\n",
    "\n",
    "Under the assumption of risk neutrality,  the price of the option is\n",
    "the expected discounted payoff:\n",
    "\n",
    "$$\n",
    "P = \\beta^n \\mathbb E \\max\\{ S_n - K, 0 \\}\n",
    "$$\n",
    "\n",
    "Now all we need to do is specify the distribution of $ S_n $, so the expectation\n",
    "can be calculated.\n",
    "\n",
    "Suppose we know that $ S_n \\sim LN(\\mu, \\sigma) $ and $ \\mu $ and $ \\sigma $ are known.\n",
    "\n",
    "If $ S_n^1, \\ldots, S_n^M $ are independent draws from this lognormal distribution then, by the law of large numbers,\n",
    "\n",
    "$$\n",
    "\\mathbb E \\max\\{ S_n - K, 0 \\}\n",
    "    \\approx\n",
    "    \\frac{1}{M} \\sum_{m=1}^M \\max \\{S_n^m - K, 0 \\}\n",
    "$$\n",
    "\n",
    "We suppose that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36625e4e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "μ = 1.0\n",
    "σ = 0.1\n",
    "K = 1\n",
    "n = 10\n",
    "β = 0.95"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225eeed0",
   "metadata": {},
   "source": [
    "We set the simulation size to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "761a4459",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "M = 10_000_000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e42ac1d",
   "metadata": {},
   "source": [
    "Here is our code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82a5595e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "S = np.exp(μ + σ * np.random.randn(M))\n",
    "return_draws = np.maximum(S - K, 0)\n",
    "P = β**n * np.mean(return_draws)\n",
    "print(f\"The Monte Carlo option price is approximately {P:3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f4162fc",
   "metadata": {},
   "source": [
    "## Pricing via a dynamic model\n",
    "\n",
    "In this exercise we investigate a more realistic model for the share price $ S_n $.\n",
    "\n",
    "This comes from specifying the underlying dynamics of the share price.\n",
    "\n",
    "First we specify the dynamics.\n",
    "\n",
    "Then we’ll compute the price of the option using Monte Carlo."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31c2c90c",
   "metadata": {},
   "source": [
    "### Simple dynamics\n",
    "\n",
    "One simple model for $ \\{S_t\\} $ is\n",
    "\n",
    "$$\n",
    "\\ln \\frac{S_{t+1}}{S_t} = \\mu + \\sigma \\xi_{t+1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ S_0 $ is lognormally distributed and  \n",
    "- $ \\{ \\xi_t \\} $ is IID and standard normal.  \n",
    "\n",
    "\n",
    "Under the stated assumptions, $ S_n $ is lognormally distributed.\n",
    "\n",
    "To see why, observe that, with $ s_t := \\ln S_t $, the price dynamics become\n",
    "\n",
    "\n",
    "<a id='equation-s-mc-dyms'></a>\n",
    "$$\n",
    "s_{t+1} = s_t + \\mu + \\sigma \\xi_{t+1} \\tag{20.1}\n",
    "$$\n",
    "\n",
    "Since $ s_0 $ is normal and $ \\xi_1 $ is normal and IID, we see that $ s_1 $ is\n",
    "normally distributed.\n",
    "\n",
    "Continuing in this way shows that $ s_n $ is normally distributed.\n",
    "\n",
    "Hence $ S_n = \\exp(s_n) $ is lognormal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10a6573b",
   "metadata": {},
   "source": [
    "### Problems with simple dynamics\n",
    "\n",
    "The simple dynamic model we studied above is convenient, since we can work out\n",
    "the distribution of $ S_n $.\n",
    "\n",
    "However, its predictions are counterfactual because, in the real world,\n",
    "volatility (measured by $ \\sigma $) is not stationary.\n",
    "\n",
    "Instead it rather changes over time, sometimes high (like during the GFC) and sometimes low.\n",
    "\n",
    "In terms of our model above, this means that $ \\sigma $ should not be constant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c83abb",
   "metadata": {},
   "source": [
    "### More realistic dynamics\n",
    "\n",
    "This leads us to study the improved version:\n",
    "\n",
    "$$\n",
    "\\ln \\frac{S_{t+1}}{S_t} = \\mu + \\sigma_t \\xi_{t+1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\sigma_t = \\exp(h_t),\n",
    "    \\quad\n",
    "        h_{t+1} = \\rho h_t + \\nu \\eta_{t+1}\n",
    "$$\n",
    "\n",
    "Here $ \\{\\eta_t\\} $ is also IID and standard normal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bcc92f0",
   "metadata": {},
   "source": [
    "### Default parameters\n",
    "\n",
    "For the dynamic model, we adopt the following parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c4e2ef4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "default_μ  = 0.0001\n",
    "default_ρ  = 0.1\n",
    "default_ν  = 0.001\n",
    "default_S0 = 10\n",
    "default_h0 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a9e932",
   "metadata": {},
   "source": [
    "(Here `default_S0` is $ S_0 $ and `default_h0` is $ h_0 $.)\n",
    "\n",
    "For the option we use the following defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "535c39f6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "default_K = 100\n",
    "default_n = 10\n",
    "default_β = 0.95"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11702310",
   "metadata": {},
   "source": [
    "### Visualizations\n",
    "\n",
    "With $ s_t := \\ln S_t $, the price dynamics become\n",
    "\n",
    "$$\n",
    "s_{t+1} = s_t + \\mu + \\exp(h_t) \\xi_{t+1}\n",
    "$$\n",
    "\n",
    "Here is a function to simulate a path using this equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85f1be3d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def simulate_asset_price_path(μ=default_μ, S0=default_S0, h0=default_h0, n=default_n, ρ=default_ρ, ν=default_ν):\n",
    "    s = np.empty(n+1)\n",
    "    s[0] = np.log(S0)\n",
    "\n",
    "    h = h0\n",
    "    for t in range(n):\n",
    "        s[t+1] = s[t] + μ + np.exp(h) * randn()\n",
    "        h = ρ * h + ν * randn()\n",
    "\n",
    "    return np.exp(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6cbac91",
   "metadata": {},
   "source": [
    "Here we plot the paths and the log of the paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c7b5287",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1)\n",
    "\n",
    "titles = 'log paths', 'paths'\n",
    "transforms = np.log, lambda x: x\n",
    "for ax, transform, title in zip(axes, transforms, titles):\n",
    "    for i in range(50):\n",
    "        path = simulate_asset_price_path()\n",
    "        ax.plot(transform(path))\n",
    "    ax.set_title(title)\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82b99275",
   "metadata": {},
   "source": [
    "### Computing the price\n",
    "\n",
    "Now that our model is more complicated, we cannot easily determine the\n",
    "distribution of $ S_n $.\n",
    "\n",
    "So to compute the price $ P $ of the option, we use Monte Carlo.\n",
    "\n",
    "We average over realizations $ S_n^1, \\ldots, S_n^M $ of $ S_n $ and appealing to\n",
    "the law of large numbers:\n",
    "\n",
    "$$\n",
    "\\mathbb E \\max\\{ S_n - K, 0 \\}\n",
    "    \\approx\n",
    "    \\frac{1}{M} \\sum_{m=1}^M \\max \\{S_n^m - K, 0 \\}\n",
    "$$\n",
    "\n",
    "Here’s a version using Python loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bee8207f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_call_price(β=default_β,\n",
    "                       μ=default_μ,\n",
    "                       S0=default_S0,\n",
    "                       h0=default_h0,\n",
    "                       K=default_K,\n",
    "                       n=default_n,\n",
    "                       ρ=default_ρ,\n",
    "                       ν=default_ν,\n",
    "                       M=10_000):\n",
    "    current_sum = 0.0\n",
    "    # For each sample path\n",
    "    for m in range(M):\n",
    "        s = np.log(S0)\n",
    "        h = h0\n",
    "        # Simulate forward in time\n",
    "        for t in range(n):\n",
    "            s = s + μ + np.exp(h) * randn()\n",
    "            h = ρ * h + ν * randn()\n",
    "        # And add the value max{S_n - K, 0} to current_sum\n",
    "        current_sum += np.maximum(np.exp(s) - K, 0)\n",
    "\n",
    "    return β**n * current_sum / M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82c63436",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "compute_call_price()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "060435e0",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd1731e",
   "metadata": {},
   "source": [
    "## Exercise 20.1\n",
    "\n",
    "We would like to increase $ M $ in the code above to make the calculation more\n",
    "accurate.\n",
    "\n",
    "But this is problematic because Python loops are slow.\n",
    "\n",
    "Your task is to write a faster version of this code using NumPy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "145fea78",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.1](https://intro.quantecon.org/#monte_carlo_ex1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "273bbf96",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_call_price_vector(β=default_β,\n",
    "                       μ=default_μ,\n",
    "                       S0=default_S0,\n",
    "                       h0=default_h0,\n",
    "                       K=default_K,\n",
    "                       n=default_n,\n",
    "                       ρ=default_ρ,\n",
    "                       ν=default_ν,\n",
    "                       M=10_000):\n",
    "\n",
    "    s = np.full(M, np.log(S0))\n",
    "    h = np.full(M, h0)\n",
    "    for t in range(n):\n",
    "        Z = np.random.randn(2, M)\n",
    "        s = s + μ + np.exp(h) * Z[0, :]\n",
    "        h = ρ * h + ν * Z[1, :]\n",
    "    expectation = np.mean(np.maximum(np.exp(s) - K, 0))\n",
    "\n",
    "    return β**n * expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51bdcfe4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "compute_call_price_vector()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cf51b31",
   "metadata": {},
   "source": [
    "Notice that this version is faster than the one using a Python loop.\n",
    "\n",
    "Now let’s try with larger $ M $ to get a more accurate calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da12baec",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "compute_call_price(M=10_000_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee8e959b",
   "metadata": {},
   "source": [
    "## Exercise 20.2\n",
    "\n",
    "Consider that a European call option may be written on an underlying with spot price of \\$100 and a knockout barrier of \\$120.\n",
    "\n",
    "This option behaves in every way like a vanilla European call, except if the spot price ever moves above \\$120, the option “knocks out” and the contract is null and void.\n",
    "\n",
    "Note that the option does not reactivate if the spot price falls below \\$120 again.\n",
    "\n",
    "Use the dynamics defined in [(20.1)](#equation-s-mc-dyms) to price the European call option."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fab56f0",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 20.2](https://intro.quantecon.org/#monte_carlo_ex2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2183202a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "default_μ  = 0.0001\n",
    "default_ρ  = 0.1\n",
    "default_ν  = 0.001\n",
    "default_S0 = 10\n",
    "default_h0 = 0\n",
    "default_K = 100\n",
    "default_n = 10\n",
    "default_β = 0.95\n",
    "default_bp = 120"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dadd4d04",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_call_price_with_barrier(β=default_β,\n",
    "                                    μ=default_μ,\n",
    "                                    S0=default_S0,\n",
    "                                    h0=default_h0,\n",
    "                                    K=default_K,\n",
    "                                    n=default_n,\n",
    "                                    ρ=default_ρ,\n",
    "                                    ν=default_ν,\n",
    "                                    bp=default_bp,\n",
    "                                    M=50_000):\n",
    "    current_sum = 0.0\n",
    "    # For each sample path\n",
    "    for m in range(M):\n",
    "        s = np.log(S0)\n",
    "        h = h0\n",
    "        payoff = 0\n",
    "        option_is_null = False\n",
    "        # Simulate forward in time\n",
    "        for t in range(n):\n",
    "            s = s + μ + np.exp(h) * randn()\n",
    "            h = ρ * h + ν * randn()\n",
    "            if np.exp(s) > bp:\n",
    "                payoff = 0\n",
    "                option_is_null = True\n",
    "                break\n",
    "\n",
    "        if not option_is_null:\n",
    "            payoff = np.maximum(np.exp(s) - K, 0)\n",
    "        # And add the payoff to current_sum\n",
    "        current_sum += payoff\n",
    "\n",
    "    return β**n * current_sum / M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea98f5c3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time compute_call_price_with_barrier()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d143d75b",
   "metadata": {},
   "source": [
    "Let’s look at the vectorized version which is faster than using Python loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c473ac",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_call_price_with_barrier_vector(β=default_β,\n",
    "                                           μ=default_μ,\n",
    "                                           S0=default_S0,\n",
    "                                           h0=default_h0,\n",
    "                                           K=default_K,\n",
    "                                           n=default_n,\n",
    "                                           ρ=default_ρ,\n",
    "                                           ν=default_ν,\n",
    "                                           bp=default_bp,\n",
    "                                           M=50_000):\n",
    "    s = np.full(M, np.log(S0))\n",
    "    h = np.full(M, h0)\n",
    "    option_is_null = np.full(M, False)\n",
    "    for t in range(n):\n",
    "        Z = np.random.randn(2, M)\n",
    "        s = s + μ + np.exp(h) * Z[0, :]\n",
    "        h = ρ * h + ν * Z[1, :]\n",
    "        # Mark all the options null where S_n > barrier price\n",
    "        option_is_null = np.where(np.exp(s) > bp, True, option_is_null)\n",
    "\n",
    "    # mark payoff as 0 in the indices where options are null\n",
    "    payoff = np.where(option_is_null, 0, np.maximum(np.exp(s) - K, 0))\n",
    "    expectation = np.mean(payoff)\n",
    "    return β**n * expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a023a3e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%time compute_call_price_with_barrier_vector()"
   ]
  }
 ],
 "metadata": {
  "date": 1730689204.0845785,
  "filename": "monte_carlo.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Monte Carlo and Option Pricing"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}