{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Training \"inspecteur modèles\" BNP Paribas - Monte-Carlo methods - January 7, 2020" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Notebook 4 - Examples of Nested Monte-Carlo and empirical regression methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 1. Nested Monte-Carlo for expected loss\n", "\n", "\n", "#### $\\blacktriangleright$ Setting\n", "We consider a Black-Scholes model at discrete dates $t_i$:\n", "\n", "$$\n", "S_{t_i} = S_0 \\, e^{\\sigma \\, W_{t_i} - \\frac 12 \\sigma^2 t_i},\n", "\\qquad t_0 < t_1 < \\dots < t_n,\n", "$$\n", "\n", "which we use to generate option prices (note we consider zero interest rate $r=0$). \n", "\n", "In particular, the price at time $t_1$ of a Put option with maturity $t_2$ and strike price $K$ is given by \n", "\n", "$$\n", "P(t_i, S_{t_i}) = \\mathbb E \\bigl[ (K - S_{t_2})^+ \\big| S_{t_1} \\bigr]\n", "$$\n", "\n", "where $x^+$ denotes the positive part of $x$; the explicit Black-Scholes formula is of course available for the function $P(t,S)$\n", "\n", "$$\n", "P(t,S) = K \\mathcal N \\biggl(\\frac{\\ln \\frac K {S}}{\\sigma \\sqrt{t_2 - t}} + \\frac 12 \\sigma \\sqrt{t_2 - t} \\biggr)\n", "- S \\, \\mathcal N \\biggl(\\frac{\\ln \\frac K {S}}{\\sigma \\sqrt{t_2 - t}} - \\frac 12 \\sigma \\sqrt{t_2 - t} \\biggl),\n", "\\qquad t \\le t_2, \\ S > 0\n", "$$\n", "\n", "where $\\mathcal N$ is the standard gaussian cdf." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We wish to apply Monte-Carlo simulation to the evaluation of the _nested_ expectation\n", "\n", "$$\n", "I = \\mathbb E \\Bigl[\\Bigl(P(t_1, S_{t_1}) - P(t_0, S_0) \\Bigr)^+ \\Bigr]\n", "$$\n", "\n", "which represents the _expected loss_ (or exposure) at time $t_1$ of a short position on the put option of maturity $t_2$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as sps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1. Benchmark value (semi-explicit in this case)\n", "\n", "In the framework of the model above, it is possible to evaluate the nested expectation $I$ using numerical evaluation of a 1D integral with a deterministic quadrature formula, since\n", "\n", "$$\n", "I = \\int_{\\mathbb R} \\Bigl(P \\bigl(t_1, S_0 \\, e^{\\sigma \\, \\sqrt{t_1} \\, y \\, - \\, \\frac 12 \\sigma^2 t_1} \\bigr) - P(t_0, S_0) \\Bigr)^+\n", "\\frac{e^{-\\, y^2/2}}{\\sqrt{2 \\pi}} \\, dy\n", "$$\n", "\n", "and the integrand function is explicit.\n", "\n", "We can numericall evaluate the integral on the right hand side using the function `scipy.integrate.quad` (we suggest to check out the input and output formats of this function on the online documentation).\n", "\n", "We will use this value as a benchmark value for $I$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Option parameters\n", "t_2 = 2.\n", "K = 100.\n", "S_0 = 100.\n", "\n", "# Model parameter\n", "sigma = 0.3\n", "\n", "# Intermediate date\n", "t_1 = t_2/2\n", "\n", "def putPrice(t, S, T, K, sigma):\n", " \"\"\"\n", " Black-Scholes price at time t for put option with maturity T.\n", " \"\"\"\n", " deltaT = T - t\n", " sigmaSqrtDeltaT = sigma * np.sqrt(deltaT)\n", " \n", " d = np.log(K/S) / sigmaSqrtDeltaT + sigmaSqrtDeltaT/2\n", " \n", " return K * sps.norm.cdf(d) - S * sps.norm.cdf(d - sigmaSqrtDeltaT)\n", "\n", "##############################################################\n", "# Evaluation of the benchmark value with scipy.integrate.quad \n", "##############################################################\n", " \n", "##############################################################\n", "# TO DO: evaluate the initial price P(0, S_0)\n", "price_at_zero = 0. # modify here\n", "\n", "##########################################\n", "# To Do: complete the integrand function\n", "##########################################\n", "def to_integrate(y):\n", " \n", " S_t_1 = S_0 * np.exp(sigma*np.sqrt(t_1) * y - 0.5*sigma*sigma*t_1)\n", " \n", " loss_t_1 = putPrice(t=t_1, S=S_t_1, T=t_2, K=K, sigma=sigma) - price_at_zero\n", "\n", " density = np.exp(-0.5*y*y) / np.sqrt(2*np.pi)\n", "\n", " return 0 # update here\n", " \n", "###################################################\n", "# To Do:\n", "# + Import the function quad from scipy.integrate\n", "# + Evalue the value of the integral I\n", "\n", "benchmark_value = 0. # modify here\n", "\n", "print(\"\\n Benchmark value for the expected loss at time t_1: %1.3f\" %benchmark_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2. Nested Monte-Carlo estimator\n", "\n", "In the definition of $I$, we can recognize a nested expectation of the form\n", "\n", "$$\n", "\\mathbb E \\bigl[ g \\bigl( E[f(X,Y)|X] \\bigr) \\bigr]\n", "$$\n", "\n", "where $f,g$ are two given functions and $X$ and $Y$ are two independent random variables, since\n", "\n", "$$\n", "I = \\mathbb E \\biggl[ \n", "\\biggl( \\mathbb E \\Bigl[ \\Bigl(K - S_{t_1} e^{\\sigma \\, (W_{t_2} - W_{t_1}) - \\frac 12 \\sigma^2 (t_2 - t_1)} \\Bigr)^+ \\Big| S_{t_1} \\Bigr] - P(t_0, S_0) \\biggr)^+\n", "\\biggr]\n", "$$\n", "\n", "and the Brownian increment $W_{t_2} - W_{t_1}$ is independent of $S_{t_1}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider\n", " \n", "+ $M$ i.i.d. samples $S_m$ of $S_{t_1}$\n", "\n", " \n", "+ $M \\times N$ i.i.d. samples of the Brownian increment , denoted $\\Delta W_{m,n}$\n", "\n", "The idea behind the __Nested Monte-Carlo estimator__ is the following: for every value $S_m$ sampled at $t_1$, we use the vector of $N$ samples $(\\Delta W_{m,n})_{n = 1, \\dots, N}$ of the Brownian increment to approximate the _inner_ conditional expectation in the expression above.\n", "Then, we have to take the outer mean with respect to the simulated values of $S_{t_1}$.\n", "\n", "The resulting Nested Monte-Carlo estimator is\n", "\n", "$$\n", "I_{M,N} = \\frac 1 M \\sum_{m=1}^M\n", "\\biggl(\n", "\\frac 1 N \\sum_{n=1}^N\n", "\\Bigl(K - S_m \\, e^{\\sigma \\, \\Delta W_{m,n} \\, - \\, \\frac 12 \\, \\sigma^2 (t_2 - t_1)} \\Bigr)^+\n", "- P(t_0, S_0)\n", "\\biggr)^+.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Implementation: \n", "\n", "In the cell below, complete the function `NestedEstimator` so that it outputs one single sample of the estimator $I_{M,N}$.\n", "\n", "Then, generate a sample of $k$ i.i.d. values of the estimator $I_{M,N}$, choosing\n", "\n", "+ first $M=N$ (and for example $M=1000$)\n", "\n", "+ then a smaller value of $N$, of the order of $\\sqrt M$.\n", "\n", "In both case, give an estimate of the RMSE \n", "\n", "$$\n", "\\sqrt{ \\mathbb E \\Bigl[ \\bigl(I - I_{M,N} \\bigr)^2 \\Bigr]}\n", "$$\n", "\n", "and plot the histogram of the $k$ simulated values of the estimator $I_{M,N}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def NestedEstimator(t_1, t_2, K, S_0, sigma, M, N):\n", " \"\"\"\n", " This function outputs one sample of the \n", " Nested MC estimator I_{M,N}\n", " \"\"\"\n", " ###############################################\n", " # To Do: generate the M i.i.d. samples of S_t_1\n", " G = np.random.normal(loc=0, scale=1, size=M)\n", " \n", " S_1 = 0 # update here\n", " \n", " ######################################\n", " # To Do: simulate M*N i.i.d. samples \n", " # of the Brownian increment Delta W\n", " ######################################\n", " Delta_W = 0 # update here\n", " \n", " ##############################################\n", " # We evaluate the Nested Monte-Carlo estimator\n", " ##############################################\n", " \n", " # This variable will contain the outer sum over {m = 1, ..., M}\n", " outerSum = 0.\n", " \n", " for m in range(M):\n", " # We work with one-dimensional arrays of size N\n", " exponential_factor = np.exp(sigma * Delta_W[m] - 0.5*sigma*sigma*(t_2 - t_1))\n", " \n", " S_2 = S_1[m] * exponential_factor\n", " \n", " ##################################################\n", " ## To Do: implement the inner Monte-Carlo mean 1/N sum_{n = 1, ..., N}\n", " innerMean = 0 # update here\n", " \n", " ##################################################\n", " ## To do: update the outer sum_{m = 1, ..., M}\n", " outerSum += 0 # update here\n", " \n", " Nested_estimator = outerSum / M\n", " \n", " return Nested_estimator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\blacktriangleright$ Simulation of the estimator" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "################################################\n", "# Parameters M et N of the Nested MC estimator\n", "M = int(1.e3)\n", "N = M\n", "\n", "################################################\n", "# We generate a sample of k i.i.d. values tirages\n", "# of the estimator \n", "k = 100\n", "\n", "sample_nested_estim = np.zeros(k)\n", "\n", "from time import time\n", "time0 = time()\n", "\n", "##########################################\n", "## To Do: generate the samples of I_{M,N}\n", "for i in range(k):\n", " sample_nested_estim[i] = 0 # update here\n", "\n", "time1 = time()\n", "\n", "###########################################\n", "# To Do: estimate the RMSE from the sample\n", "# of k values of I_{M,N} generated above\n", "\n", "MSE = 0 # update here\n", "\n", "RMSE = np.sqrt(MSE)\n", "\n", "#######################################\n", "# Some printing of the results\n", "print(\"**Nested estimator \\n\")\n", "print(\"M = %1.0f outer samples and N = %1.0f inner samples \\n\" %(M,N) )\n", "print(\"k = %1.0f samples of the Nested estimator \\n\" %k )\n", "print(\"Estimateur I_{M,N} (1 sample) = %1.3f \\n\" %sample_nested_estim[0] )\n", "print(\"RMSE = %1.3f \\n\" %RMSE )\n", "print(\"Time = %1.2f \\n\" %(time1-time0) )\n", "\n", "#####################################\n", "# To Do: histogram of the estimator\n", "\n", "plt.hist( ??? , density=True, histtype='step', bins=int(np.sqrt(k)), label=\"Nested MC estimator\")\n", "\n", "plt.axvline(benchmark_value, linewidth=1.5, color='k', label=\"True value\")\n", "\n", "plt.legend(loc=9, fontsize=15, bbox_to_anchor=(1.4, 1.0), ncol=1)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Remark (pricing measure and historical measure)\n", "\n", "The estimator above can also be applied to the evaluation of the (proper) expected loss\n", "\n", "$$\n", "\\mathbb E^{\\mathbb P^{\\mathrm{hist}}} \\Bigl[\\Bigl(P(t_1, S_{t_1}) - P(t_0, S_0) \\Bigr)^+ \\Bigr]\n", "$$\n", "\n", "where $\\mathbb P^{\\mathrm{hist}}$ is the physical measure of the market, under which we assume that the underlying spot price follows the dynamics\n", "\n", "$$\n", "S_{t_i} = S_0 \\, e^{\\sigma \\, W_{t_i} \\, + \\, (\\mu - \\frac 12 \\sigma^2 ) t_i},\n", "\\hspace{30mm} (1)\n", "$$\n", "\n", "where $\\mu$ is a drift parameter.\n", "\n", "Note that we still have a nested expectation problem, the difference being that now the outer expectation $\\mathbb E^{\\mathbb P^{\\mathrm{hist}}}$ and the inner expectation $P(t_1, S_{t_1}) = \\mathbb E \\bigl[ (K - S_{t_2})^+ \\big| S_{t_1} \\bigr]$ are computed according to different probability measures.\n", "\n", "$\\blacktriangleright$ __Question__: How would you modify the estimator and the code above so to take the historical model $(1)$ into account? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2. Dynamic programming with empirical regression: the example of the Put Bermuda" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Setting\n", "\n", "In the framework of the discrete-time Black-Scholes model of the previous exercise\n", "\n", "$$\n", "S_i = S_0 \\, e^{\\sigma \\, W_i - \\frac 12 \\sigma^2 t_i},\n", "\\qquad t_0 < t_1 < \\dots < t_n\n", "$$\n", "\n", "we wish to evaluate the price of a __Bermudan put option__ of maturity $T=t_n$, which \n", "\n", "+ can be exerced by its owner at any date $t_i$ previous to $T$\n", "\n", "\n", "+ if exercised at $t_i$, generates the gain $(K - S_i)^+$.\n", "\n", "This pricing problem for such an option is an example within the more general family of _optimal stopping problems_ (the owner looks for an optimal exercise time)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ The recursive dynamic programming equation\n", "\n", "The theory of optimal stopping allow to represent the price $Y_i$ of the Bermudan option at any date $t_i$ in terms of the corresponding _dynamical programming equation_: starting from the final value $Y_n = (K - S_n)^+$, we have to evaluate recursively\n", "\n", "$$\n", "\\left\\{\n", "\\begin{aligned}\n", "&Y_n = (K - S_n)^+\n", "\\\\\n", "&Y_i = \\max \\left( (K - S_i)^+, \\mathbb E \\bigl[Y_{i+1} \\big| S_i \\bigr] \\right)\n", "\\qquad i = n-1, \\dots, 0.\n", "\\end{aligned}\n", "\\right.\n", "$$\n", "\n", "The equation above\n", "\n", "+ Requires to evaluate the conditional expectation $\\mathbb E[Y_{i+1} | S_i]$ ($Y_{i+1}$ being known from the evaluation already performed at date $t_{i+1}$).\n", "\n", "\n", "+ is arguably still simpler to implement than the representation of the option price as\n", "$$\n", "Y_0 = \\sup_{\\tau \\in \\mathcal{T}_n } {\\mathbb E}\\bigl[ (K - S_\\tau)^+ \\bigr],\n", "$$\n", "where $\\mathcal{T}_n$ is the set of all possible (random) stopping times with values in $\\{0,t_1,\\dots, t_n\\}$.\n", "\n", "\n", "Note that the dynamic programming for optimal stopping is a specific example of more general recursive equations of the form\n", "\n", "$$\n", "Y_i = f\\left( t_i, S_i, \\mathbb E \\bigl[Y_{i+1} \\big| S_i \\bigr] \\right)\n", "$$\n", "\n", "where $f$ is a given function, that appear when simulating _backward_ stochastic differential equations (for example in the context of Credit valuation adjuments, as in this afternoon lecture)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ The pricing problem in a simple two-period framework\n", "\n", "When $n=2$, the Bermudan put price is given by\n", "\n", "$$\n", "\\begin{aligned}\n", "Y_0 &= \\max \\left\\{ (K-S_0)^+, \\mathbb E \\left[ Y_1 \\right] \\right\\}\n", "\\\\\n", "&= \\max \\left\\{ (K-S_0)^+, \\mathbb E \\left[ \\max \\Bigl( (K - S_1)^+, \\mathbb E \\left[ (K - S_2)^+ |S_1 \\right] \\Bigr) \\right] \\right\\}\n", "\\end{aligned}\n", "$$\n", "\n", "The problem is analogous to one considered in the previous exercise: we are in front of a nested expectation of the form $\\mathbb E \\bigl[ g \\bigl(X, \\mathbb E[f(X,Y)|X] \\bigr) \\bigr]$ (now for a different outer function $g$). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Benchmark price\n", "\n", "Exactly as in the previous exercise, in the present two-period setting we have a semi-explicit formula for the Bermudan option price, which we can use a reference value for other estimators.\n", "\n", "Here again, we exploit the explicit knowledge of the function (the European put price)\n", "\n", "$$\n", "\\mathbb E \\left[ (K - S_2)^+ |S_1 \\right] = P(t_1, S_1)\n", "$$\n", "\n", "which allows to define the estimator\n", "\n", "$$\n", "Y_0^{M} = \\max \\Bigl\\{\n", "(K-x_0)^+, \\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, P(t_1, S_1^m) \\right) \\Bigr\\}\n", "$$\n", "\n", "based on a sample $(S_{1}^m)_{m = 1, \\dots, M}$ of i.i.d. values of $S_1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the cells below:\n", "\n", "+ For future use of empirical regression and for visualisation purposes, we rewrite the put price as a function of the Brownian motion $W_1$ at time $t_1$ rather than the spot price $S_1$ (the one being a function of the other): that is, we define\n", "$$\n", "v_1(W_1) := P(t_1, S_1) = P \\Bigl( t_1, S_0 \\, e^{\\sigma \\, W_1 \\, - \\, \\frac 12 \\sigma^2 t_1} \\Bigr)\n", "$$\n", "so that the estimator above reads\n", "$$\n", "Y_0^{M} = \\max \\Bigl\\{\n", "(K-x_0)^+, \\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, v_1(W_1^m) \\right) \\Bigr\\}\n", "$$\n", "$S_{1}^m$ and $W_1^m$ being related by $S_{1}^m = S_0 \\, e^{\\sigma \\, W_1^m \\, - \\, \\frac 12 \\sigma^2 t_1}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Problem parameters\n", "K = 1.2\n", "S_0 = 1.\n", "sigma = 0.2\n", " \n", "t_2 = 2.\n", "t_1 = t_2 / 2.\n", " \n", "M = int(5e3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def putBlackScholes(S_0, t_1, t_2, K, W_1, sigma):\n", " \"\"\"\n", " Black Scholes price at t_1 of a put option with maturity t_2,\n", " as a function of the underlying Brownian motion\n", " at time t_1.\n", " \n", " W_1: numpy array\n", " \"\"\"\n", " S_1 = S_0 * np.exp(sigma*W_1 - 0.5*sigma*sigma*t_1)\n", " \n", " price = putPrice(t=t_1, S=S_1, T=t_2, K=K, sigma=sigma)\n", " \n", " return price\n", "\n", "###################################################\n", "# The function below is useful to obtain the values\n", "# of S at time t once the values of W_t are given.\n", "#\n", "# GBM stands for \"Geometric Brownian Motion\"\n", "###################################################\n", "def GBM(S_0, t, sigma, W):\n", " return 0 # update here\n", " \n", "time0 = time()\n", " \n", "#################################################\n", "# To Do: generate a sample of M values of W_1\n", "# and the corresponding sample of values of S_1\n", "W_1 = np.zeros(M) # update here\n", "S_1 = np.zeros(M) # update here\n", "\n", "################################################\n", "\n", "time0_1 = time()\n", " \n", "timeSimulations = time0_1 - time0\n", " \n", "###################################################\n", "### 1. Prix Benchmark\n", "###################################################\n", " \n", "time1 = time()\n", " \n", "v_1 = putBlackScholes(S_0, t_1, t_2, K, W_1, sigma)\n", " \n", "################################################\n", "## TO DO: evalute the benchmark price Y_0^M\n", "\n", "benchmark_price = 0.0 # update here\n", "\n", "var_TCL = 0.0 # update here\n", "################################################\n", " \n", "radiusIC = 1.96 * np.sqrt(var_TCL/M)\n", " \n", "time2 = time()\n", " \n", "print(\"Benchmark price= %1.4f +/- %1.4f \\n\" %(benchmark_price, radiusIC) )\n", "print(\"Erreur relative (TCL) = %1.3f \\n\" %(radiusIC / benchmark_price) )\n", "print(\"Time: %1.4f \\n\" %(time2 - time1 + timeSimulations) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Price by empirical regression\n", "\n", "$\\blacktriangleright$ In this question, we proceed as if we did not know explicitly the conditional expectation $\\mathbb E \\left[ (K - S_2)^+ \\big|S_1 \\right]$ (which is the realistic case in practice, for more general stochastic models and more complex portfolios).\n", "\n", "$\\blacktriangleright$ We therefore have to _approximate_ the conditional expectation, which we can do applying __least square linear regression__ to a simulated i.i.d. sample of the couple $(W_1, (K - S_2)^+)$.\n", "\n", "\n", "$\\blacktriangleright$ We will therefore need to simulate a sample $(W_1^m, S_2^m)_{1 \\le m \\le M}$ of i.i.d. values of the couple $(W_1, S_2)$.\n", "\n", "One can note that, in this part, we work with $M$ samples in total - as opposed to the $M \\times N$ samples of the Nested Monte-Carlo procedure." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Choice of the regression space\n", "\n", "+ We are going to build an empirical approximation $\\tilde v_1$ of the function\n", "$$\n", "x \\mapsto v_1(x) = \n", "P \\left(T_1, e^{\\sigma \\, x \\, - \\, \\frac12 \\sigma^2 T_1 } \\right)\n", "$$\n", "on the state space of the Brownian motion $W_1$. \n", "\n", "\n", "+ We choose as an approximation space the space generated by the indicator functions of disjoint intervals $I_k$\n", "$$\n", "\\varphi_k(x) = 1_{I_k}(x),\n", "\\quad\n", "\\quad k = 1, \\dots, n,\n", "$$\n", "where\n", "$$\n", "I_k = \\left[-a + (k-1) \\delta, -a + k \\delta \\right[, \\qquad \\delta = \\frac{2a}n\n", "$$\n", "so that the union of the $I_k$'s is the (large) interval $[-a,a]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The regression problem for $v(x) = \\mathbb E \\left[ (K - S_2)^+ \\big|W_1 = x \\right]$\n", "\n", "+ As stated above, we approximate the function $v_1$ using the elements of the vector space\n", "$$\n", "\\Phi = \\text{Vect}(\\varphi_1,\\dots, \\varphi_n)\n", "= \\left\\{\\sum_{k=1}^n \\alpha_k 1_{I_k}(\\cdot): \\alpha_1, \\dots, \\alpha_n \\in \\mathbb{R} \\right\\}.\n", "$$\n", "that is, using piece-wise constant functions of the form $\\sum_{k=1}^n \\alpha_k 1_{I_k}(\\cdot)$.\n", "\n", "\n", "\n", "+ The __least square regression coefficients__ for the function $v(x) = \\mathbb E \\left[ (K - S_2)^+ |W_1 = x \\right]$ are defined by the minimization problem\n", "$$ \\label{e:regrEmp}\n", "\\begin{aligned}\n", "\\alpha^* = \\underset{\\alpha\\in \\mathbb{R}^n}{\\rm arg\\min}\n", "\\ \\frac 1 M \\sum_{m = 1}^M\n", "\\left(\n", "\\bigl(K - S_2^m \\bigr)^+ - \\sum_{k=1}^n \\alpha_k 1_{I_k}(W_1^m)\n", "\\right)^2.\n", "\\end{aligned}\n", "$$\n", "Once the $\\alpha^*_k$ are known, the regression function (that is, the empirical approximation of $v_1(\\cdot)$) is given by\n", "$$\n", "\\tilde v_1(\\cdot) = \\sum_{k=1}^{n} \\alpha_k^* \\ 1_{I_k}(\\cdot)\n", "$$\n", "\n", "\n", "+ Finally, once the regression function $\\tilde v_1$ has been evaluated, the new estimator of the Bermudan put is\n", "$$\n", "\\tilde{Y}_0^{M} = \n", "\\max \\Bigl\\{\n", "(K - x_0)^+, \\frac 1M \\sum_{m=1}^M \\max \\left( (K - S_{1}^m)^+, \\tilde{v}_1(W_1^m) \\right)\n", "\\Bigr\\}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\blacktriangleright$ __Question__: do we have an explicit expression for the regression coefficients $\\alpha^*_k$? \n", "\n", "In the cell below, evaluate\n", "\n", "+ the coefficients $\\alpha^*_k$ in the Python function `EmpiricalRegression`\n", "+ the Bermudan price estimator $\\tilde Y_0^{M}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "##################################################\n", "### 2. Bermudan put price by empirical regression \n", "##################################################\n", "\n", "#################################################\n", "# First of all: \n", "# + Generate a sample of M values of W_2\n", "# and the corresponding values of S_1\n", "# + Taking into account that the values of W_1\n", "# have already been simulated above\n", "#################################################\n", "W_2 = np.zeros(M) # update here\n", "S_2 = np.zeros(M) # update here\n", "################################################\n", "\n", "###############################\n", "## We generate the cells I_k\n", "###############################\n", "a = 3. * np.sqrt(t_1)\n", " \n", "#################################################\n", "# n = the dimension of the approximation space\n", "# \n", "# It can be shown that n = M*{1/3} is a pertinent\n", "# choice for a problem in dimension d = 1, as here\n", "#################################################\n", "n = int(M**(1./3))\n", " \n", "partition_points = np.linspace(-a, a, n+1)\n", " \n", "####################################################\n", "## TO DO: evaluate\n", "## + The regression coefficients alpha\n", "## + The vector v_1_tilde(W1) of size M \n", "## inside the array v_1_tilde\n", "####################################################\n", "def EmpiricalRegression(W_1, W_2, partition_points, t_2, K, S_0, sigma):\n", " step = partition_points[1] - partition_points[0]\n", " \n", " alpha = np.zeros(partition_points.size - 1)\n", " \n", " v_1_tilde = np.zeros(M)\n", " \n", " for k in range(n):\n", " \n", " leftPoint = partition_points[k]\n", " \n", " insideCell = np.logical_and( leftPoint <= W_1, W_1 < leftPoint+step ) \n", " \n", " GBM_t_2 = GBM(S_0, t_2, sigma, W_2[insideCell])\n", " \n", " #############################################\n", " ## TO DO: evaluate the coefficients alpha[k]\n", " #############################################\n", " \n", " \n", " #############################################\n", " ## TO DO: evaluate the empirical approximation\n", " ## v_1_tilde(W1) inside the array v_1_tilde\n", " #############################################\n", " \n", " return alpha, v_1_tilde\n", " \n", "time3 = time()\n", "\n", "alpha, v_1_tilde = EmpiricalRegression(W_1, W_2, partition_points, t_2, K, S_0, sigma) \n", "\n", "################################################\n", "# TO DO: evaluate the new estimator Y_tilde_0^M \n", "# of the Bermudan put price\n", "\n", "price_EmpRegr = 0.0 #update here\n", "\n", "################################################\n", " \n", "time4 = time()\n", " \n", "print(\"Price by empirical regression = %1.4f\" %price_EmpRegr)\n", "print(\"Time: %1.4f \\n\" %(time4 - time3 + timeSimulations))\n", " \n", "#######################################################\n", "## We can also plot the true function v_1 and its\n", "## empirical approximation v_1_tilde for comparison\n", "#######################################################\n", "x = np.linspace(-a, a, 100)\n", " \n", "v_1_exact = putBlackScholes(S_0, t_1, t_2, K, x, sigma)\n", " \n", "plt.plot(x, v_1_exact, color=\"b\", label=\"$v_1$\")\n", " \n", "plt.step(intervals, np.append(alpha, alpha[-1]), where=\"post\", color=\"r\", label=r\"$\\tilde{v}_1$\")\n", " \n", "plt.xlabel(\"values of $W_1$\", fontsize=17)\n", "plt.legend(loc=\"best\", fontsize=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Nested Monte-Carlo price\n", "\n", "In the spirit of Exercise 1, we can also propose a Nested Monte-Carlo estimator that does not require the knowledge of the inner conditional expectation $\\mathbb E \\left[ (K - S_2)^+ \\big|W_1 = x \\right]$.\n", "\n", "This third estimator is\n", "\n", "$$\n", "\\hat{Y}_0^{M} = \\max \\Bigl\\{ (K-x_0)^+,\n", "\\frac 1M \\sum_{m = 1}^M \\max\\left( (K - S_{1}^m)^+, \\hat{v}_{1}^m \\right)\n", "\\Bigr\\}\n", "$$\n", "\n", "where now, for every $m$, the estimate $\\hat{v}_{1}^m$ of $v_1(W_m^1)$ is obtained from a sample of $N$ i.i.d. values $(\\Delta W^{m, n})_{1 \\le n \\le N}$ of the Brownian increment between $t_1$ and $t_2$:\n", "\n", "$$\n", "\\hat{v}_{1}^m = \\frac1{N} \\sum_{n = 1}^N \\Bigl(K - S_0 \\exp\n", "\\Bigl(-\\frac12 \\sigma^2 T_2 \\ + \\sigma (W_1^m + \\Delta W^{m, n} \\Bigr)\n", "\\Bigr)^+.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\blacktriangleright$ Simulate the Nested estimator $\\hat{Y}_0^{M}$ in the cell below, and compare with the previous methods." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "##################################\n", "## 3. Nested Monte-Carlo price\n", "##################################\n", " \n", "###################################################\n", "# We now have to generate N samples of S_2 for every\n", "# point in the sample of M values of W_1.\n", "#\n", "# We will therefore need to repeat N gaussian iid samples\n", "# of the Browian increment W_2 - W_1\n", "# for every value of W_1 previously simulated\n", "###################################################\n", "N = int(np.sqrt(M))\n", "\n", "# This variable will contain the sum over m in {1,...,M}\n", "sum_1 = 0.\n", " \n", "time7 = time()\n", " \n", "for m in range(M):\n", " w1 = W_1[m]\n", " s1 = GBM(S_0=S_0, t=t_1, sigma=sigma, W=w1)\n", " \n", " ###################################################\n", " ## To Do: \n", " ## - Simulation of N Brownian increments\n", " ## - Generation of W_2, therefore of S_2, conditionally to W1 = w1\n", " ## - update the sum over m with the current contribution\n", " ###################################################\n", " W2 = np.zeros(N) # update here\n", " \n", " sum_1 += 0 # update here\n", " ###################################################\n", "\n", "Y_0_hat = np.maximum( np.maximum(K-S_0, 0.), sum_1 / M )\n", " \n", "time8 = time()\n", " \n", "print(\"Nested MC price = %1.4f\" %Y_0_hat)\n", "print(\"Time: %1.4f\" %(time8 - time7), \"\\n\")" ] } ], "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }