{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Training \"inspecteur modèles\" BNP Paribas - Monte-Carlo methods - December 12-17/2019" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TP 3 - Continuous time processes and their approximation schemes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1. Euler scheme\n", "\n", "We consider the process $(S_t)_{t \\ge 0}$ satisfying the equation\n", "\n", "$$\n", "dS_t = r \\, S_tdt + \\sigma \\, S_t dW_t,\n", "\\qquad S_0=s_0>0,\n", "\\qquad r,\\sigma>0,\n", "$$\n", "\n", "where $(W_t)_{t \\ge 0}$ is a real valued Brownian motion.\n", "Recall that the stochastic differential equation (SDE) above admits an explicit solution:\n", "\n", "$$\n", "S_t = s_0 \\exp \\left( \\bigl(r-\\frac{\\sigma^2}{2} \\bigr)t + \\sigma W_t \\right), \\qquad t \\ge 0,\n", "$$\n", "\n", "called geometric Brownian motion.\n", "\n", "The explicit knowledge of the solution allows to\n", "\n", "+ _Exactly_ simulate both the solution of the equation and its Euler scheme.\n", "\n", "+ Compare their trajectories, and estimate the error of the Euler scheme.\n", "\n", "We fix a time horizon $T$ and a number $n$ of discretization steps, yielding the grid\n", "\n", "$$ t_k = k \\ \\frac T n $$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 1 (a) (fundamental) :\n", "\n", "Write the evolution of the true process \n", "\n", "$$ \\bigl( S_{t_k} \\bigr)_{0\\leq k\\leq n} $$\n", "\n", "on the time grid $t_k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 1 (b) (fundamental as well) :\n", "\n", "Write the evolution of the Euler scheme\n", "\n", "$$ \\bigl(S^n_{t_k} \\bigr)_{0\\leq k\\leq n} $$\n", "\n", "on the same time grid $t_k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Implementation: \n", "\n", "Questions 1 (a) and (b) provide recursive simulation schemes for both the true process $S$ and its Euler scheme, which can be simulated using the same Brownian increments.\n", "\n", "In the cells below, implement the simulation schemes given above, and compare their trajectories of $S$ and the Euler scheme $S^n$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = 1. \n", "sigma = 0.3\n", "r = 0.02\n", "S_0 = 100" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#########################################\n", "## Discretisation dates t_1, ..., t_n\n", "n = 52\n", "dates = np.linspace(0, T, n+1)\n", "\n", "# Time step\n", "Delta_T = T/n\n", "\n", "#################################################\n", "# Simulation of the trajectories\n", "# S: the true process\n", "# Se: the Euler scheme\n", "#################################################\n", "# We initialize vectors with n+1 components\n", "S = S_0 * np.ones(n+1)\n", "Se = S_0 * np.ones(n+1)\n", "\n", "for i in range(n):\n", " ## Simulation of a standard Gaussian distribution N(0,1)\n", " G = np.random.randn()\n", " \n", " ##############################################################\n", " ## TO DO: generate the Brownien increment on [t_i, t_{i+1}]\n", " ## using the random variable G\n", " ##############################################################\n", " accr_brownien = np.sqrt(Delta_T) * G\n", " \n", " ############################################################\n", " ## TO DO: complete the simulation\n", " # + of the trajectory of the process S\n", " # + of the trajectory of the Euler scheme Se\n", " ############################################################\n", " \n", " ############################################################\n", " # Process S\n", " S[i+1] = S[i]*np.exp( (r-sigma*sigma/2)*Delta_T + sigma*accr_brownien )\n", " \n", " ############################################################\n", " # Euler scheme\n", " Se[i+1] = Se[i] * (1 + r*Delta_T + sigma*accr_brownien)\n", "\n", "###########################\n", "## Plot the trajectories\n", "###########################\n", "import matplotlib.pyplot as plt\n", "\n", "###########################################################\n", "plt.plot(dates, S, color=\"b\", label=\"Process S\")\n", "plt.plot(dates, Se, color=\"r\", label=r\"Euler scheme $S^n$\",\n", " marker = \".\", linewidth=0.,)\n", "###########################################################\n", "\n", "plt.xlabel(\"discretization dates t_i\", fontsize=12)\n", "\n", "plt.legend(loc=\"best\", fontsize=12)\n", "\n", "######################################\n", "# Comparison of final values: \n", "# remove the if 0\n", "if 1:\n", " print(\"Final value S_T of the true process: %1.2f\" %S[-1])\n", " print(\"Final value S^n_T of the Euler scheme: %1.2f\" %Se[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Question 2:\n", "\n", "Our goal is now to observe the error at final time $T$:\n", "\n", "$$\n", "error(n) = \\mathbb{E} \\Bigl[ \\bigl|S_T - S^n_T \\bigr| \\Bigr]\n", "$$\n", "\n", "called the _strong_ error of the Euler scheme.\n", "\n", "We can estimate $error(n)$ with a standard Monte-Carlo: replace the expectation with the empirical mean over a large number $N$ of simulated trajectories\n", "\n", "$$\n", "\\frac 1 N \\sum_{i=1}^N \\bigl|S^i_T - S^{n,i}_T \\bigr|\n", "$$\n", "\n", "of the process and the Euler scheme." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def errorEuler(T, sigma, r, S_0, n, N):\n", " \"\"\"\n", " Euler scheme with n steps for geometric Brownian motion.\n", " \n", " This function: estimation of the strong error E[|S_T - S^n_T|]\n", " from a sample of N iid simulations of the Euler scheme.\n", " \n", " Output of the function: the couple (strong error, 95% confidence interval)\n", " \"\"\"\n", " S = S_0*np.ones(N)\n", " Se = S_0*np.ones(N)\n", " strong_error = np.zeros(N)\n", " \n", " Delta_T = T / n\n", " \n", " #################################################### \n", " # We can now recursively compute the evolution of the\n", " # whole iid sample of N simultaneously: \n", " # recall S and Se are vectors of size N\n", " \n", " for k in range(n):\n", " g = np.random.normal(loc=0, scale=1, size=N)\n", " #######################################################\n", " # To do: implement the evolution of the Euler scheme Se\n", " # and the process S\n", " #\n", " # REMARK: we do not need to stock the whole trajectories.\n", " # We just require the terminal values S_T and S^n_T\n", " # in order to evaluate the final error\n", " #######################################################\n", " brownian_increment = np.sqrt(Delta_T) * g\n", " \n", " S = S * np.exp( (r-sigma*sigma/2)*Delta_T + sigma*brownian_increment )\n", " \n", " Se = Se * (1 + r*Delta_T + sigma*brownian_increment)\n", " \n", " #######################################################\n", " # The loop for ends here.\n", " # Now: evaluate the empirical approximation of the error \n", " #######################################################\n", " sample = np.abs(S - Se) \n", " \n", " strong_error = np.mean( sample )\n", " \n", " #############################################\n", " # To Do: evaluate the 95% confidence interval\n", " ############################################\n", " empirical_variance = np.mean(sample**2) - strong_error**2\n", " \n", " # ANOTHER POSSIBILITY\n", " empirical_variance = np.var(sample)\n", " \n", " radius_conf_int = 1.96 * np.sqrt(empirical_variance / N)\n", " \n", " return strong_error, radius_conf_int" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cell below\n", "\n", "+ Estimates the strong error for different values of $n$ (number of time steps): $n = 1, 2, 4, ..., 2^A$ (here for $A = 6$).\n", "\n", "\n", "+ Plots the graph of the strong error +/- the confidence interval, as a function of the time step $\\frac T n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = int(1.e4) # size of the iid sample\n", "\n", "n = 1 # initial number of time steps for the Euler scheme\n", "\n", "## The number of time steps will evolve as n = 1, 2, 4, ..., 2^A\n", "A = 6\n", "\n", "number_steps = np.zeros(A)\n", "\n", "strong_errors = np.zeros(A) \n", "radii_conf_intervals = np.zeros(A) \n", "\n", "for i in range(A):\n", " number_steps[i] = n\n", " \n", " #############################################\n", " # UPDATE here\n", " #############################################\n", " strong_errors[i], radii_conf_intervals[i] = errorEuler(T, sigma, r, S_0, n, N)\n", "\n", " n = 2*n # multiplication of the number of time steps by 2\n", "\n", "#############################################\n", "# Plot the error +/- the confidence interval\n", "# as functions of the time step T / n\n", "#############################################\n", "\n", "\n", "plt.plot(T/ number_steps, strong_errors, color=\"r\", label=\"Strong error\")\n", "plt.plot( T/ number_steps , strong_errors - radii_conf_intervals, color=\"b\", label=\"95% conf interval\")\n", "plt.plot( T/ number_steps , strong_errors + radii_conf_intervals, color=\"b\")\n", "\n", "plt.xlabel('Discretisation step')\n", "plt.legend(loc=\"best\", fontsize = 12)\n", "plt.show()\n", "\n", "\n", "print(\"Number of steps n:\"); print(number_steps)\n", "print(\"Errors of the Euler scheme with step n:\"); print(strong_errors)\n", "print(\"95% conf interval:\"); print (radii_conf_intervals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercice 2. Multi-level Monte-Carlo method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now consider the application of the Multi-level Monte-Carlo method to the evaluation of\n", "\n", "$$\\mathbb{E}[f(S_T)]$$\n", "\n", "where $f(s)=(K-s)^+$ and $K$ is fixed. $S$ is the process considered in the previous exercise.\n", "\n", "Reminders: the Multi-level method is based on a combination of Euler schemes with decreasing step\n", "\n", "$$h_l = \\frac T {n_l} = \\frac T{2^l}.$$\n", "\n", "The corresponding time grids are given by $t_k^l = k \\,h_l = k \\frac T{2^l}$, for $k = 0, \\dots, 2^l$.\n", "\n", "Precisely, the multi-level estimator of $\\mathbb{E}[f(S_T)]$ is given by\n", "\n", "$$\n", "\\widehat{f(S_T)}_{\\mathrm{ML}}\n", "= \\frac1{M_0} \\sum_{m=1}^{M_0} (K - S_T^{h_0,0,m})^+\n", "+ \\sum_{l=1}^L \\frac1{M_l} \\sum_{m=1}^{M_l}\n", "\\left( (K - S_T^{h_l,l,m})^+ - (K - S_T^{h_{l-1},l,m})^+ \\right)\n", "$$ \n", "\n", "where\n", "- for each $h$ and $l$, the $(S^{h,l,m})_{1 \\le m \\le M_l}$ are iid samples of the Euler scheme with step $h$;\n", "\n", "\n", "- all the simulations are independent between different levels $l$;\n", "\n", "\n", "- the schemes $S^{h_l,l,m}$ and $S^{h_{l-1},l,m}$ are constructed from the same Brownian increments.\n", "\n", "\n", "We recall that, using the the properties of the Euler scheme (weak error of order $=1$, strong error of order $=1/2$), one can show the following upper bound for the quadratic error of the multi-level estimator:\n", "\n", "$$\n", "\\mathbb{E} \\Bigl[\n", "\\bigl(\n", "\\widehat{f(S_T)}_{\\mathrm{ML}} - \\mathbb{E}[(K - S_T)^+]\n", "\\bigr)^2 \\Bigr]\n", "\\le\n", "C \\biggl( h_L^2 + \\sum_{l=0}^L \\frac{h_l}{M_l} \\biggr)\n", "\\qquad \\forall M_0, \\dots, M_L,\n", "$$\n", "\n", "where $C$ is a constant independent from $L$ and the $M_l$.\t\t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Find the smallest $L$ so to have \n", "\n", "$$h_L^2 \\le \\varepsilon^2$$\n", "\n", "where $\\varepsilon$ is a fixed desired precision.\n", "\n", "The solution $(M_0, \\dots, M_L)$ of the optimization problem\n", "\n", "$$\n", "\\min_{M_0, \\dots, M_L} \\sum_{l=0}^L \\frac{M_l}{h_l} \\qquad \\mbox{ such that } \\qquad\n", "\\sum_{l=0}^L \\frac{h_l}{M_l} = \\varepsilon^2\n", "$$\n", "\n", "yields\n", "\n", "$$\n", "M_0 \\sim_c \\frac{|\\log \\varepsilon|}{\\varepsilon^2},\n", "\\qquad\n", "M_l = M_0 \\, 2^{-l}.\n", "$$\n", "\n", "What is the interest of the solution of this problem, for the application of the Multi-level method?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. The code in the cells below generates $P$ iid samples of the multi-level estimator\n", "\n", "Complete the code:\n", "\n", "+ implement the evolution between the times $t^{l-1}_k = T \\frac{k}{2^{l-1}}$ and $t^{l-1}_{k+1} = T \\frac{k+1}{2^{l-1}}$ \n", "\n", "+ of a matrix *Se\\_fin* (resp. *Se\\_gr*) of size $M_l\\times P$, containing the Euler scheme of finer discretization step $\\frac T{2^{l}}$ (resp. of coarser discretization step $\\frac T{2^{l-1}}$)\n", "\n", "+ using the same Brownian increments.\n", "\n", "\n", "+ Finally, add to the variable *estim\\_multi* the sum\n", "\n", "$$\\sum_{m = 1}^{M_l} \\bigl( (K - S_T^{h_l,l,m})^+ - (K - S_T^{h_{l-1},l,m})^+ \\bigr)$$\n", "\n", "divided by $M_l$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "T = 1.\n", "sig = 0.2\n", "r = 0.05\n", "S0 = 100.\n", "K = 110" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def estimMultiniveaux(eps, P):\n", " \"\"\"\n", " This function: \n", " \n", " outputs an iid sample of P simulations of the Multi-level estimator\n", " \n", " having a desired target precision of order eps\n", " \"\"\" \n", " h_0 = T # initial value of the discretization step: just one step at level l=0 !\n", " \n", " ####################################################\n", " # To do: fix the total number L of levels \n", " # and the number of samples M_0 wished at level 0\n", " ####################################################\n", " \n", " # Number of levels\n", " L = int( np.log( T/eps ) / np.log(2) )\n", " \n", " # Sample size at level l=0\n", " M_0 = int( np.abs(np.log(eps)) / eps**2 )\n", " \n", " # We start with the simulation at level l=0\n", " g = np.random.randn(M_0, P)\n", " \n", " Se_coarse = S0*(1. + r*h_0 + sig * np.sqrt(h_0)*g)\n", "\n", " estim_multi = np.mean((K-Se_coarse)*(K>Se_coarse), axis=0)\n", " \n", " ##########################################\n", " # This is the loop on the different levels \n", " for l in range(1, L+1):\n", " M_l = int(M_0/2.**l)\n", " \n", " #############################################\n", " # Parameters of Euler scheme with finer step\n", " #############################################\n", " h_fine = h_0/2.**l\n", " \n", " Se_fine = S0*np.ones((M_l, P)) # Euler scheme: fine step h_l\n", " Se_coarse = S0*np.ones((M_l, P)) # Euler scheme: coarse step h_{l-1} \n", " \n", " ######################################################\n", " # This is the loop over the coarse discretization grid t_k \n", " # \n", " # Two Euler schemes evolve: the fine and the coarse \n", " #\n", " # Here, they are both arrays of size M_l x P\n", " #\n", " for k in range(2**(l-1)):\n", " \n", " g1 = np.random.randn(M_l, P)\n", " g2 = np.random.randn(M_l, P)\n", " ####################################################\n", " # Implement the evolution of the two Euler schemes\n", " ###################################################\n", " brown_increment_1 = g1 * np.sqrt( h_fine )\n", " brown_increment_2 = g2 * np.sqrt( h_fine )\n", " \n", " # Two steps for the finer Euler scheme\n", " Se_fine = Se_fine * (1 + r*h_fine + sig*brown_increment_1)\n", " \n", " Se_fine = Se_fine * (1 + r*h_fine + sig*brown_increment_2)\n", " \n", " # Only one step for the coarse Euler scheme\n", " brown_increment = brown_increment_1 + brown_increment_2\n", " \n", " Se_coarse = Se_coarse * (1 + r*2*h_fine + sig*brown_increment)\n", " \n", " ###################################################################\n", " # Add to the variable estim_multi the contribution at level l\n", " # \n", " # (See formulas above: we want the mean over the columns\n", " # of the M_l x P arrays)\n", " ###################################################################\n", " term1 = np.mean((K - Se_fine)*(K > Se_fine), axis=0)\n", " term2 = np.mean((K - Se_coarse)*(K > Se_coarse), axis=0)\n", " \n", " estim_multi = estim_multi + term1 - term2\n", " \n", " return estim_multi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. The code in the cell below represents the histogram of $P=1000$ simulations of the Multi-level estimator with a fixed desired precision threshold $\\varepsilon$ divided by $2$ from an histogram to the other.\n", "\n", "The variable *RMSE* contains the empirical estimation of the root mean square error as a function of $\\varepsilon$, where\n", "the true value of the expectation $\\mathbb{E}[(K - S_T)^+]$ is (in this special case!) computed exactly using the Black-Scholes formula\n", "\n", "$$\n", "\\mathbb{E}[(K-S_T)^+] = K \\, {\\cal N}(-d+\\sigma\\sqrt{T}) - s_0e^{rT} {\\cal N}(-d)\n", "\\ \\mbox{ where } \\ d = \\frac{\\ln(s_0/K)+ (r+\\frac{\\sigma^2}2) T}{\\sigma\\sqrt{T}}\n", "$$\n", "\n", "and with ${\\cal N}(x)=\\int_{-\\infty}^xe^{-\\frac{y^2}{2}}\\frac{dy}{\\sqrt{2\\pi}}$.\n", "\n", "\n", "Does the error evolve as expected in terms of $\\varepsilon$ ?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Root mean square error (eps=0.20) : 2.13\n", "Root mean square error (eps=0.10) : 0.92\n", "Root mean square error (eps=0.05) : 0.41\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAH/tJREFUeJzt3Xl0VeXZ/vHvbQhTAiSYAAIBIgooIFQoWKwVbYsoFXAoMooWYYna1xZbRa2ILaWIA9RVhxeHgohT1VegTCrCz2EVBKwog5RBAiGWKTgAJgS4f38kHDOckAOc5CQ712etszh77yf73Dk5uXjy7OExd0dERILltFgXICIi0adwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgFUI1YvnJKS4q1atYrVy0s1s2rVKgC6dOkS40pETs2qVav2uHtqWe0sVrcf6Nq1q69cuTImry3Vj5kBoNttSFVnZqvcvWtZ7TQsIyISQAp3EZEAUriLiASQwl1EJIBidraMSHVw9OhRMjMzOXDgQKxLkSoiPj6eRo0aUb9+/VPaj8JdpBzt2bMHM6Nt27acdpr+UJbjc3e+++47duzYAXBKAa9Pm0g5+uqrr2jcuLGCXSJiZtStW5dmzZqxa9euU9qXPnEi5ejIkSPEx8fHugypYurUqUNeXt4p7UPhLlLOjl1AJRKpaHxmFO5SrXSc0THWJUiAvP/++7Rt2zbWZYSlcBepplq1asU777wT6zJiZunSpTRv3vyEvsbM2LRpU2j5oosuYsOGDdEuLSoU7iIS1uHDh2NdgpyCMsPdzJ4zs11mtqaU7WZmj5nZJjP71MzOj36ZIhJNw4YNY9u2bVx55ZUkJiYyefJktm7dipnx7LPP0qJFCy699NKwvdvCPf6jR48yadIkWrduzemnn86AAQPIzs4u9XX/+c9/0rlzZ5KSkujRoweffvopAJs3b6Zhw4Z8/PHHAGRlZZGSksLSpUsB6NmzJ3fffTfdunWjQYMG9OvXr8jrLFu2jB49epCUlESnTp1CXweQnZ3NjTfeSNOmTUlOTqZ///4cOHCAyy+/nKysLBITE0lMTCQrK4uPPvqIH/3oRyQlJXHGGWdw2223cejQIQB+8pOfANCpUycSExN55ZVXSrw/69evp2fPniQlJdG+fXvmzJkT2nbDDTdw66230qdPH+rVq0f37t3ZvHnzif7oIufux30APwHOB9aUsv0KYAFgwAXA8rL26e506dLFRSoK4IB3mN6hQl933bp1Ffp6J6Jly5b+9ttvh5a/+OILB3zYsGG+f/9+P3jwoC9ZssSbNWtW6tdNmTLFu3fv7tu3b/ecnBwfNWqUDxw4MOzrrVq1ylNTU33ZsmV++PBhnz59urds2dJzcnLc3X3atGnerl07P3DggPfq1cvvuOOO0NdefPHF3rRpU//ss898//79fvXVV/uQIUPc3T0zM9MbNmzo8+bN8yNHjvhbb73lDRs29F27drm7+xVXXOEDBgzw7OxsP3TokC9dutTdPez3tnLlSv/Xv/7leXl5/sUXX3i7du18ypQpoe2Ab9y4MbRceB+HDh3y1q1b+5///GfPzc31xYsXe2Jion/++efu7j58+HBPTk725cuXe15eng8ePNivu+66Un8+pX12gJUeQcaW2SB/X7Q6Trj/LzCo0PIG4Iyy9qlwl4pUWcL9WB3l/YhEaeG+efPm0Lqywr1du3b+zjvvhLZlZWV5jRo1PC8vr8Tr3Xzzzf6HP/yhyLo2bdqEwtbd/corr/QOHTp4x44dQ6Hvnh/ud911V2h57dq1Hh8f74cPH/ZJkyb50KFDi+y3V69ePn36dM/KynIz8+zs7BL1hPveipsyZYr3798/tHy8cH/vvfe8cePGfuTIkdD2gQMH+v333+/u+eE+YsSI0LZ58+Z527ZtS33tUw33aIy5NwO2F1rOLFhXgpmNMrOVZrZy9+7dUXhpEYm2tLS0iNtmZGRw1VVXkZSURFJSEueccw5xcXHs3LkzbNtHHnkk1DYpKYnt27eTlZUVajNy5EjWrFnDr3/9a2rVqlVqXS1btiQvL489e/aQkZHBP/7xjyL7/eCDD/jyyy/Zvn07DRs2JDk5OaLv5z//+Q+/+MUvaNKkCfXr1+eee+5hz549EX1tVlYWaWlpRS5Ya9myZehqU4AmTZqEntetW5f9+/dHtO+TEY1wD3dCZtgZEdx9mrt3dfeuqallTiQiEjiR9Lii8YhEaedSF16fkJDAwYMHQ8tHjhyhcMcsLS2NBQsW8NVXX4UeOTk5NGtWsn+XlpbGvffeW6TtwYMHGTRoEAD79+/nN7/5DSNGjGD8+PElxu63b/++D7lt2zbi4+NJSUkhLS2NYcOGFdnvgQMHGDt2LGlpaWRnZ/PVV19F9P2PHj2adu3asXHjRr755hsmTpwY8fvZtGlTtm/fztGjR4vUGe69qAjRCPdMoPB/9c2BrFLaikgl0bhxY7Zs2XLcNm3atCEnJ4d58+aRl5fHhAkTyM3NDW2/+eabuffee8nIyABg9+7dzJ49O+y+Ro4cyVNPPcXy5ctxdw4cOMC8efP49ttvAbj99tvp0qULzzzzDH369OHmm28u8vUvvPAC69at4+DBg4wbN45rr72WuLg4hg4dyty5c1m0aBFHjhwhJyeHpUuXkpmZyRlnnMHll1/OLbfcwr59+8jLy+O9994Lff979+7l66+/Dr3Gt99+S/369UlMTOTzzz/nySefjPg96969OwkJCUyePJm8vDyWLl3K3LlzGThw4HHf4/ISjXCfA1xfcNbMBcDX7v5lFPYrIuXo7rvvZsKECSQlJfHwww+HbdOgQQOeeOIJbrrpJpo1a0ZCQkKRs0Nuv/12+vbtS69evahXrx4XXHABy5cvD7uvrl278vTTT3PbbbeRnJzMWWedxfTp0wGYPXs2Cxcu5KmnngLg0Ucf5eOPP2bWrFmhrx82bBg33HADTZo0IScnh8ceewzI/4tg9uzZTJw4kdTUVNLS0njooYdCPeiZM2cSHx9Pu3btaNSoEVOnTgWgXbt2DBo0iDPPPJOkpCSysrJ4+OGHefHFF6lXrx4jR47kuuuuK/I9jB8/nuHDh5OUlMSrr75aZFvNmjWZM2cOCxYsICUlhVtuuYXnn3+edu3aRfojiaoy51A1s5eAnkAKsBO4H4gHcPenLP9vm78BvYGDwI3uXubkqJpDVSrSsT/BO0zvwGfDP6uw112/fj3nnHNOhb1eUPXs2ZOhQ4dy0003xbqUClPaZyfSOVTLvOWvuw8qY7sDt5a1HxERqTi6QlVEJIA0WYeIVHqFrziVyKjnLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgA6Tx3kQrWauy8ct3/1kl9ynX/FSE3N5fRo0fz2muvUbduXe68807GjBkTtu2MGTN47LHH2LhxI/Xr12fw4MFMnDiRGjXy4y07O5sRI0bw1ltvkZKSwl/+8hcGDx5ckd9OTKjnLoHUcUbHcmkrFWP8+PFs3LiRjIwMlixZwuTJk1m4cGHYtgcPHmTq1Kns2bOH5cuXs3jx4iI3Qrv11lupWbMmO3fuZNasWYwePZq1a9dW1LcSMwp3kWoqKyuLa665htTUVNLT00N3WYT8cL322mu57rrrqFevHueffz6rV68ObX/wwQdp1qwZ9erVo23btixevDiqtT3//PPcd999JCcnc8455zBy5MjQHSSLGz16NBdddBE1a9akWbNmDBkyhA8//BCAAwcO8Prrr/OnP/2JxMREfvzjH9O3b19mzpwZ1XorI4W7SDV09OhRrrzySjp16sSOHTtYvHgxU6dOZdGiRaE2s2fP5pe//CXZ2dkMHjyY/v37k5eXx4YNG/jb3/7GihUr+Pbbb1m0aBGtWrUK+zqTJk0qMkNS8Uc4+/btIysri06dOoXWderUKeLe9nvvvUf79u2B/JmV4uLiaNOmzUntqypTuItUQytWrGD37t2MGzeOmjVrcuaZZzJy5EhefvnlUJsuXbpw7bXXEh8fz5gxY8jJyWHZsmXExcWRm5vLunXryMvLo1WrVrRu3Trs64wdO7bIDEnFH+Ecm3quQYMGoXUNGjQITepxPH//+99ZuXIlv/vd70L7KryfE9lXVadwF6mGMjIyyMrKKtKLnjhxYpG5TwvPWXraaafRvHlzsrKyOOuss5g6dSrjx4+nUaNGDBw4sMg8qKcqMTERgG+++Sa07ptvvqFevXrH/bo333yTsWPHhibLOLavwvuJdF9BoHAXqYbS0tJIT08v0ov+9ttvmT9/fqhN4TlLjx49SmZmJk2bNgVg8ODBfPDBB2RkZGBm3HXXXWFfZ+LEiSQmJpb6CCc5OZkzzjijyBj/6tWrQ0Mt4SxcuJCRI0cyd+5cOnb8/gB5mzZtOHz4MBs3box4X0GhcBephrp160b9+vV58MEH+e677zhy5Ahr1qxhxYoVoTarVq3ijTfe4PDhw0ydOpVatWpxwQUXsGHDBt59911yc3OpXbs2derUIS4uLuzr3HPPPezfv7/UR2muv/56JkyYwL59+/j88895+umnueGGG8K2fffddxkyZAivv/463bp1K7ItISGBq6++mnHjxnHgwAE+/PBDZs+ezbBhw078TatqKmo29uKPLl26uEh56TC9Q5FlwIES68O1jaZ169aV275P1Y4dO3zgwIHeuHFjT0pK8u7du/vbb7/t7u7333+/X3PNNT5gwABPTEz0zp07+6pVq9zdffXq1f7DH/7QExMTPTk52fv06eM7duyIam05OTl+4403er169bxRo0b+yCOPhLZlZGR4QkKCZ2RkuLt7z549PS4uzhMSEkKP3r17h9rv3bvX+/Xr53Xr1vW0tDSfNWtWVGstL6V9doCVHkHG6iImkWqqadOmvPTSS6Vur127Ni+88EKJ9eeddx4fffRReZZGrVq1eO6553juuedKbGvRokWRXv+SJUuOu6+GDRvy5ptvRr3Gyk7DMiIiAaRwFxEJIA3LiEgJ48ePj3UJcorUcxcRCSCFu4hIACncRUQCSOEugVf8lr4dZ3TUbX4l8BTuIiIBpHAXkUonNzeXX/3qV9SvX58mTZrw6KOPltp2zZo1XHbZZaSkpGBmJbZnZ2dz1VVXkZCQQMuWLXnxxRfLs/RKQ6dCilS08Q3KbnNK+/+6fPdfAQrPxPTf//6XSy65hHPPPZfevXuXaBsfH8+AAQO45ZZb6N+/f4nthWdi+uSTT+jTpw+dOnUK/M3D1HMXqaaCMhNT27ZtGTFiRNiw1kxMIgGhg6WRCfJMTIVpJqYymFlvM9tgZpvMbGyY7S3MbImZ/dvMPjWzK6JfqohES1BnYgq3L83EVAoziwMeBy4HzgUGmdm5xZr9AXjV3X8ADASeiHahIhI9QZyJqbR9aSam0nUDNrn7Fnc/BLwM9CvWxoH6Bc8bANH7SYtI1AVtJqbSaCam42sGbC+0nFmwrrDxwFAzywTmA78OtyMzG2VmK81s5e7du0+iXBGJhiDNxOTu5OTkcOjQIQBycnLIzc0FqvdMTJGEe8kTR/N76oUNAqa7e3PgCmCmmZXYt7tPc/eu7t41NTX1xKsVkaiIi4tj7ty5fPLJJ6Snp5OSksJNN93E119/fxplv379eOWVV0hOTmbmzJm88cYbxMfHk5uby9ixY0lJSaFJkybs2rWLiRMnRrW+Bx54gNatW9OyZUsuvvhifv/734dOg9y2bRuJiYls27YNyB9iqlOnTqg3XqdOHdq2bRva1xNPPMF3331Ho0aNGDRoEE8++WS16LlHcp57JpBWaLk5JYddRgC9Adz9X2ZWG0gBdkWjSJFAqSTnoQdlJqZWrVqRP/tceJqJqXQrgLPNLN3MapJ/wHROsTbbgJ8CmNk5QG1A4y4iIjFSZs/d3Q+b2W3AIiAOeM7d15rZH8mfqHUOcAfwtJn9lvwhmxv8eP+VikSZzm0XKSqi2w+4+3zyD5QWXjeu0PN1wIXRLU1EYkUzMVV9ureMBJZ681Kd6fYDIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEA0qmQIhWsvE/R/Gz4Z+W6/4qQm5vL6NGjee2116hbty533nknY8aMKbX9lClTQjdBu+aaa3jyySepVasWkH97gp07d4ZubtajRw/eeuutCvk+Ykk9dxGpdArPobpkyRImT57MwoULw7ZdtGgRkyZNYvHixWzdupUtW7Zw//33F2kzd+7c0J0oq0Owg8JdpNoKyhyqM2bMCM2hmpyczH333Vdq2+pE4S5SDQVpDtW1a9eWaLtz50727t0bWjdkyBBSU1Pp1atXkf+kgkzhLlINBWkO1eLzpB57fqz9rFmz2Lp1KxkZGVxyySVcdtllpb52kCjcRaqhIM2hWnye1GPPj7W/8MILqVOnDnXr1uXuu+8mKSmJ999/P2r1VlYKd5FqKEhzqLZv375E28aNG3P66aeHbW9mx53cIygU7iLVUJDmUL3++ut59tlnWbduHfv27WPChAmhttu2bePDDz/k0KFD5OTk8NBDD7Fnzx4uvDD4dyjXee4iFawynId+bA7VO+64g/T0dHJzc2nbti0TJkwItTk2h+rw4cM566yzSsyhun79euLj4+nRowfTpk2Lan0PPPAAo0ePpmXLltSpU4e77rqryByq5557LuvWraNFixb07t2bO++8k0suuSR0nvsDDzwA5I+7jx49ms2bN1O7dm06d+7MggULSu3VB4nCXaSaCsocqgBjxowJe5FT+/bt+fTTT8utzspMwzIiIgGkcBcRCSANy4hICZpDtepTz11EJIAU7iLlrDqcUy3RdfTo0VPeh8JdpBzVrl2bvXv3KuAlIu7OoUOH2LFjBwkJCae0L425S5VW3vdGP1XNmzcnMzOT3bt3x7oUqSJq1KhBgwYNSElJObX9RKkeEQkjPj6e9PT0WJch1ZCGZUREAkjhLiISQAp3EZEA0pi7VFuV/WCsyKlQz11EJIAU7iIiARRRuJtZbzPbYGabzGxsKW0GmNk6M1trZi9Gt0wRETkRZY65m1kc8DjwcyATWGFmc9x9XaE2ZwN3Axe6+z4za1ReBYuISNki6bl3Aza5+xZ3PwS8DPQr1mYk8Li77wNw913RLVNERE5EJOHeDNheaDmzYF1hbYA2ZvahmS0zs97RKlBERE5cJKdCWph1xe+CVAM4G+gJNAfeN7MO7v5VkR2ZjQJGQf5UWSIiUj4i6blnAmmFlpsDWWHazHb3PHf/AthAftgX4e7T3L2ru3dNTU092ZpFRKQMkYT7CuBsM0s3s5rAQGBOsTZvApcAmFkK+cM0W6JZqIiIRK7McHf3w8BtwCJgPfCqu681sz+aWd+CZouAvWa2DlgC/N7d95ZX0SIicnwR3X7A3ecD84utG1fouQNjCh4iIhJjureMVEm6L4zI8en2AyIiAaRwFxEJIIW7iEgAKdxFRAJI4S4SRscZHXXQVqo0hbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkCarEOkEN1PRoJCPXcRkQBSuIuIBJDCXUQkgBTuIiIBpAOqUqXogKdIZNRzFxEJIIW7iEgAKdxFRAJI4S4iEkA6oCoCML5B/r/pLWJbh0iUqOcuIhJACncRkQBSuIuIBFBE4W5mvc1sg5ltMrOxx2l3rZm5mXWNXokiInKiyjygamZxwOPAz4FMYIWZzXH3dcXa1QP+B1heHoWKlKeOOpAqARNJz70bsMndt7j7IeBloF+Ydn8CJgM5UaxPREROQiTh3gzYXmg5s2BdiJn9AEhz939GsTYRETlJkYS7hVnnoY1mpwFTgDvK3JHZKDNbaWYrd+/eHXmVIiJyQiIJ90wgrdBycyCr0HI9oAOw1My2AhcAc8IdVHX3ae7e1d27pqamnnzVIiJyXJGE+wrgbDNLN7OawEBgzrGN7v61u6e4eyt3bwUsA/q6+8pyqVhERMpUZri7+2HgNmARsB541d3XmtkfzaxveRcoIiInLqJ7y7j7fGB+sXXjSmnb89TLEhGRU6ErVEWOo+OMjpr9SaokhbuISAAp3EVEAkjhLiISQJqsQyq/YxNpgCbTEImQeu4iIgGkcBcRCSCFu1QZui2vSOQU7iIiAaRwFxEJIIW7iEgAKdxFRAJI57mLRKLwufbjv45dHSIRUs9dRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkh3hZTKpfDdF0XkpKnnLiISQAp3EZEAUriLiASQwl1EJIAiOqBqZr2BvwJxwDPuPqnY9jHATcBhYDfwK3fPiHKtIpVDuIO+mnpPKpkye+5mFgc8DlwOnAsMMrNzizX7N9DV3c8DXgMmR7tQERGJXCTDMt2ATe6+xd0PAS8D/Qo3cPcl7n6wYHEZ0Dy6ZYqIyImIZFimGbC90HIm0P047UcAC06lKKkmdE67SLmJJNwtzDoP29BsKNAVuLiU7aOAUQAtWrSIsESRyq/V2HkAbJ3UJ8aViOSLZFgmE0grtNwcyCreyMx+BtwL9HX33HA7cvdp7t7V3bumpqaeTL0iIhKBSHruK4CzzSwd2AEMBAYXbmBmPwD+F+jt7ruiXqVIjHVMb8FnX2wrdfvW2gW/EuMLVujsGYmxMnvu7n4YuA1YBKwHXnX3tWb2RzPrW9DsISAR+IeZfWJmc8qtYhERKVNE57m7+3xgfrF14wo9/1mU6xIRkVOgK1RFItQxvQUd03UigFQNuuWvSHnQVawSY+q5i4gEkMJdRCSAFO4iIgGkMXc5aceuyixL6KrNk7zdgA5iipw49dxFRAJIPXcpd6H7rtSOcSEi1YjCXcKKdMhFRConDcuIiASQwl1EJIA0LCMSAyd8ppHICVK4i1QQHceQiqRwF6nE1MOXk6UxdxGRAFLPvZrR0IBI9aCeu4hIACncpdLRpBgip07DMhIVoQmiC7TKeTFGlVRPkQy36aBr9aJwl3JRPOxFpGIp3KXSqqxDM8Xr+uyLbTGq5MTotMrqReEeEDoLRkQKU7iLVJCTGarSsQs5WQr3Sk49chE5GQp3ESlCY/PBoHCXE1ZeZ8JU1gOoZTlWd0UdWA33/mv4RorTRUwiIgGknruUqTx66hXd262qKvP1Ahq+qdzUcxcRCSD13GNEZ8FIdaEefmyo5y4SZVX1wLAES0Q9dzPrDfwViAOecfdJxbbXAp4HugB7gevcfWt0S6061CsXkVgrM9zNLA54HPg5kAmsMLM57r6uULMRwD53P8vMBgIPAteVR8EQuz/zFNrRF6RebpC+l1jQ8E10RdJz7wZscvctAGb2MtAPKBzu/YDxBc9fA/5mZubuHsVaJcoq85kYEn2R/Lyrwvnyur1xZCIJ92bA9kLLmUD30tq4+2Ez+xo4HdgTjSJPVnXuaSu4qxfdT7+oyv67XxH/+UQS7hZmXfEeeSRtMLNRwKiCxf1mtiGC14+FFGL8H9NJKFJzuB9I5bIGYvA+r7lhTTR2U2bdsX//f1F8RYpF9F6X+LpYqoq/hxDJ5+PBU9p/y0gaRRLumUBaoeXmQFYpbTLNrAbQAMguviN3nwZMi6SwWDKzle7eNdZ1nAjVXHGqYt2queJUlrojORVyBXC2maWbWU1gIDCnWJs5wPCC59cC72q8XUQkdsrsuReMod8GLCL/VMjn3H2tmf0RWOnuc4BngZlmton8HvvA8ixaRESOL6Lz3N19PjC/2LpxhZ7nAL+MbmkxVemHjsJQzRWnKtatmitOpajbNHoiIhI8uv2AiEgAKdwLMbMkM3vNzD43s/Vm9qNY11QWM/utma01szVm9pKZ1Y51TeGY2XNmtsvM1hRa19DM3jazjQX/JseyxuJKqfmhgs/Hp2b2f2aWFMsawwlXd6FtvzMzN7OUWNRWmtJqNrNfm9mGgs/45FjVV5pSPiOdzWyZmX1iZivNrFssalO4F/VXYKG7twM6AetjXM9xmVkz4H+Aru7egfwD3pX1YPZ0oHexdWOBxe5+NrC4YLkymU7Jmt8GOrj7ecB/gLsruqgITKdk3ZhZGvm3EamMN9GfTrGazewS8q9+P8/d2wMPx6Cuskyn5Hs9GXjA3TsD4wqWK5zCvYCZ1Qd+Qv6ZP7j7IXf/KrZVRaQGUKfg+oK6lLwGoVJw9/coee1DP2BGwfMZQP8KLaoM4Wp297fc/XDB4jLyr/uoVEp5rwGmAHcS5gLDWCul5tHAJHfPLWizq8ILK0MpdTtQv+B5A2L0O6lw/96ZwG7g72b2bzN7xswSYl3U8bj7DvJ7M9uAL4Gv3f2t2FZ1Qhq7+5cABf82inE9J+pXwIJYFxEJM+sL7HD31bGu5QS0AS4ys+Vm9v/M7IexLihCvwEeMrPt5P9+xuSvO4X792oA5wNPuvsPgANUvmGCIgrGqPsB6UBTIMHMhsa2qurBzO4FDgOzYl1LWcysLnAv+UMEVUkNIBm4APg98KqZxf7ODmUbDfzW3dOA31IwGlDRFO7fywQy3X15wfJr5Id9ZfYz4At33+3uecAbQI8Y13QidprZGQAF/1a6P7vDMbPh5N+EZUgVuRK7NfkdgNVmtpX8oaSPzaxJTKsqWybwhuf7CDhK/n1bKrvh5P8uAvyD/DvrVjiFewF3/y+w3czaFqz6KUVva1wZbQMuMLO6BT2an1LJDwIXU/i2FcOB2TGsJSIFE9fcBfR194OxricS7v6Zuzdy91bu3or80Dy/4DNfmb0JXApgZm2AmlSNG4llARcXPL8U2BiTKtxdj4IH0BlYCXxK/gcrOdY1RVDzA8Dn5N9mcSZQK9Y1lVLnS+QfF8gjP1xGkH9b6MXkf/gXAw1jXWcENW8i//bWnxQ8nop1nZHUXWz7ViAl1nVG8F7XBF4o+Gx/DFwa6zojrPvHwCpgNbAc6BKL2nSFqohIAGlYRkQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiATQ/weX4dbmU8lIYgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "eps = 0.2 # target precision\n", "P = 1000 # wished number of samples of the Multi-level estimator\n", "\n", "import scipy.stats as sps\n", "\n", "###################################################\n", "# Black Scholes explicit formula for E[(K-S_T)^+]\n", "###################################################\n", "d = (np.log(S0/K) + r*T) / (sig*np.sqrt(T)) + sig*np.sqrt(T)/2.\n", "d2 = d - sig*np.sqrt(T)\n", "expectation = K * sps.norm.cdf(-d2) - S0*np.exp(r*T) * sps.norm.cdf(-d)\n", "\n", "plt.clf()\n", "\n", "for i in range(3):\n", " estim_multi = estimMultiniveaux(eps, P)\n", " \n", " plt.hist(estim_multi, density=\"True\", bins=int(np.sqrt(P)), label='eps = %1.2f' %(eps))\n", " \n", " #########################################################\n", " # To do: estimation of the RMSE from the sample estim_multi\n", " RMSE = np.sqrt( np.mean( (estim_multi - expectation)**2 ) )\n", " #########################################################\n", " print(\"Root mean square error (eps=%1.2f) : %1.2f\" %(eps, RMSE))\n", " \n", " eps = eps/2\n", "\n", "plt.axvline(expectation, linewidth=2.0, color='k',label=\"true expectation\")\n", "plt.legend(loc=\"best\", fontsize=12)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }