2. Monte Carlo Methods for Risk Management.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo Methods for Risk Management"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Value at Risk\n",
    "\n",
    "##### Risk Metrics\n",
    "\n",
    "A technique for quantifying the risk in a portfolio is known as a risk metric. This gives us a way of numerically describing the amount of uncertainty in that portfolio.\n",
    "\n",
    "Examples: Volatility/Variance and Correlation which give an idea of the extent to which a portfolio's value could decrease and how susceptible it is to massive swings in value.\n",
    "\n",
    "However, to fully characterize these risks, we would need the distribution of returns to be multivariate normal (very unlikely). Therefore, we use a general metric as described below.\n",
    "\n",
    "##### Value at Risk (VaR)\n",
    "\n",
    "Value at Risk gives an indication of how much you stand to lose on a portfolio with a given probability, over a specific time period. This can be measured in absolute terms (we could lose $100 in portfolio value) or in relative terms (the value of our portfolio will decrease by 5%)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let our current portfolio value be $X(t)$, then we can define the return on the portfolio until some fixed time T as:\n",
    "\n",
    "$$ R_{t,T} = \\frac{X(T-t) - X(t)}{X(t)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of VaR using Alpha-Quantile of a Distribution\n",
    "\n",
    "We can define the $\\alpha$-quantile of a distribution, $x_{\\alpha}(t,T)$ as:\n",
    "\n",
    "$$ P[R_{t,T} < x_{\\alpha}(t,T)] $$\n",
    "\n",
    "which translates into our VaR values being:\n",
    "\n",
    "$$VaR_{\\alpha}(t,T) = \\begin{cases}\n",
    "               -x_{\\alpha}(t,T)\\\\\n",
    "                x_{\\alpha}(t,T)X(t)\\\\\n",
    "            \\end{cases} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can interpret VaR as being a quantile value: what values does the portfolio take on in the worst 5% of cases? We can also interpret this in terms of confidence intervals. We are (1 − α)-percent certain that the value of the portfolio will not fall below 𝑉𝑎𝑅$_{\\alpha}$(𝑡, 𝑇).\n",
    "\n",
    "**N.B.** the definition of VaR is the negative of a quantile such that the actual VaR value is positive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example:\n",
    "\n",
    "For example, suppose we have a portfolio that is worth 100 USD. Say we know that there is a 5% chance that the portfolio will have a value of 20 USD or less. Then the 5% VaR for the portfolio is 80 USD, or 80%. We are 95% confident that the value of our portfolio will not fall by more than 80 USD.\n",
    "\n",
    "This changes according to the significance level chosen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expected Shortfall\n",
    "\n",
    "Other names for this include conditional VaR or expected tail loss. Expected shortfall tells us how much we stand to lose assuming that our portfolio value does actually fall below our VaR. Put more simply, it is how much we stand to lose, on average, if we only consider the worst α-percentile of cases.\n",
    "\n",
    "$$ EL_{\\alpha}(t,T) = - E[R_{t,T} | R_{t,T} < -VaR_{\\alpha}(t,T)]\\times X(t) $$\n",
    "\n",
    "Which is the negative of our expected return, given the expected return is less than the VaR, multiplied by the portfolio value. Note that this means that we are dealing with absolute terms, and not relative terms.\n",
    "\n",
    "**N.B.** the VaR value must also be in terms of the portfolio returns, and not the portfolio value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating VaR\n",
    "\n",
    "### a. Monte Carlo VaR Models\n",
    "\n",
    "Monte Carlo VaR models make use of Monte Carlo methods to project asset values, in order to get some idea of what happens to your portfolio value in the future. In doing so, you will be able to build a distribution for your portfolio value. You can then use this distribution to estimate the α- percentile of the portfolio, and to estimate a VaR.\n",
    "\n",
    "- Determine dynamics for the assets in the portfolio\n",
    "\n",
    "\n",
    "\n",
    "- Use Monte Carlo simulation to simulate the value of these assets n-days into the future\n",
    "\n",
    "\n",
    "\n",
    "- Determine the portfolio value in each scenario\n",
    "\n",
    "\n",
    "\n",
    "- Determine the $\\alpha$-quantile of these projected portfolio values\n",
    "\n",
    "\n",
    "Monte Carlo methods don't fall prey to data history issues once you have calibrated your return model; but a model needs to be implied which immediately introduces model risk into your VaR estimate.\n",
    "\n",
    "**Monte Carlo simulation tries to find a balance between what has happened in the past, and what could potentially happen in the future. It makes use of models which can allow for events which haven't explicitly occurred before, and it also uses historical data to calibrate its model parameters, which allows for it to incorporate what has happened in the past into its estimates.**\n",
    "\n",
    "### b. Historical VaR Models\n",
    "\n",
    "This simulation process makes no parametric assumptions, and instead entirely uses data to project what may happen. \n",
    "\n",
    "The historical method doesn't require any parametric assumptions regarding the returns of your assets; however, it does need a considerable amount of data to be applied. \n",
    "\n",
    "This may force you to use data from further into the past than you would be comfortable (arguably the most recent data would be a better reflection of the future asset performance). In fact, many regulators require that the data used for historical VaR calculations include periods of “stress”. Thus, the firm may be forced to include data from high stress periods if the most recent data doesn't meet regulatory standards.\n",
    "\n",
    "Selecting appropriate data can be difficult.\n",
    "\n",
    "**Historical simulation focuses on assessing what has happened in the past in order to determine the risk we face in the future.**\n",
    "\n",
    "\n",
    "##### Example:\n",
    "\n",
    "Suppose we have three days’ worth of share data.\n",
    "\n",
    "We can work out two days’ worth of returns for that share. If we then sample from these two returns, we can then project the value of our share one day into the future. \n",
    "\n",
    "In doing so, we are assuming that the company will be subject to similar conditions as it was in the past, so the past returns will be reflective of future returns. If we were to repeat this process several times (using significantly more data), we could apply a similar process as with Monte Carlo VaR models in order to estimate VaR.\n",
    "\n",
    "\n",
    "- Source price data for the assets in the portfolio\n",
    "\n",
    "\n",
    "\n",
    "- Calculate the daily returns for all assets in the portfolio\n",
    "\n",
    "\n",
    "\n",
    "- Randomly sample $n$ many log-returns and use these to grow the respective asset values in the portfolio $n$-days into the future\n",
    "\n",
    "\n",
    "\n",
    "- Determine your portfolio value in each scenario\n",
    "\n",
    "\n",
    "- Determine the $\\alpha$-quantile of these projected portfolio values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c. Conditional Tail VaR\n",
    "\n",
    "In order to determine Conditional Tail VaR, given that you have already determined the VaR, you would just need to work out the average value of the portfolio values which are less the VaR value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Value at Risk in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import numpy.matlib\n",
    "from scipy.stats import uniform\n",
    "import matplotlib.pyplot as plt\n",
    "import math as m\n",
    "import random as r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting Up Relevant Share Information\n",
    "\n",
    "We will create column vectors for the initial share values and the volatilities for each share.\n",
    "\n",
    "- Assets Prices: `[[100], [95], [50]]`\n",
    "\n",
    "\n",
    "$$ \\begin{bmatrix}\n",
    "100 \\\\\n",
    "95 \\\\\n",
    "50  \\\\\n",
    "\\end{bmatrix} $$\n",
    "\n",
    "\n",
    "- Volatility: `[[0.15],[0.2],[0.3]]`\n",
    "\n",
    "$$ \\begin{bmatrix}\n",
    "0.15 \\\\\n",
    "0.2 \\\\\n",
    "0.3  \\\\\n",
    "\\end{bmatrix} $$\n",
    "\n",
    "- 3x3 correlation matrix: `[[1, 0.2, 0.4],[0.2, 1, 0.8],[0.4, 0.8, 1]]`\n",
    "\n",
    "\n",
    "$$ \\begin{bmatrix}\n",
    "1 & 0.2 & 0.4 \\\\\n",
    "0.2 & 1 & 0.8 \\\\\n",
    "0.4 & 0.8 & 1 \\\\\n",
    "\\end{bmatrix} $$\n",
    "\n",
    "\n",
    "- Cholesky Decomposition of this matrix as L\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General share info\n",
    "S_0 = np.array([[100],[95],[50]])\n",
    "sigma = np.array([[0.15],[0.2],[0.3]])\n",
    "corr_mat = np.array([[1, 0.2, 0.4],[0.2, 1, 0.8],[0.4, 0.8, 1]])\n",
    "L = np.linalg.cholesky(corr_mat)\n",
    "r = 0.1\n",
    "T = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to decide the **total number of simulations** and the **level of confidence for out VaR estimate**.\n",
    "\n",
    "Once this is done, we can calculate our initial portfolio value (assuming that we hold one unit of each of the three assets) as the **sum of the current asset values**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a. Applying Monte Carlo Estimation of VaR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "t_simulations = 10000\n",
    "alpha = 0.05\n",
    "# current portfolio value\n",
    "portval_current = np.sum(S_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Geometric Brownian Motion Solution\n",
    "\n",
    "$$S_T = S_0 e^{(r- \\frac{\\sigma^{2}}{2})*T+ \\sigma \\sqrt{T}Z}$$\n",
    "\n",
    "where T is the terminal price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def terminal_shareprice(S_0, risk_free_rate, sigma, Z, T):\n",
    "    \"\"\"\n",
    "    Generates the terminal share price given some random normal values, z\n",
    "    \"\"\"\n",
    "    # It returns an array of terminal stock prices.\n",
    "    return S_0*np.exp((risk_free_rate-sigma**2/2)*T+sigma*np.sqrt(T)*Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate future share values and portfolio values (assuming that we keep our holdings constant). The “Z” variable is a correlated multivariate normal sample.\n",
    "\n",
    "Where `norm.rvs(size=[3, t_simulations]`\n",
    "\n",
    "creates a 3xt_simulations matrix so 3x10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating 10000 simulations for future portfolio values\n",
    "Z = np.matmul(L, norm.rvs(size= [3, t_simulations]))\n",
    "portval_future = np.sum(terminal_shareprice(S_0, r, sigma, Z, T), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([195.2392385 , 314.83192009, 233.49453515, ..., 221.4042465 ,\n",
       "       309.34317272, 238.4214139 ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "portval_future"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the Expected Portfolio Returns Given the Simulations\n",
    "\n",
    "Before trying to calculate the VaR, it would make things significantly easier to sort this list. Otherwise, we would need to run some form of search algorithm to find the smallest portfolio returns. Lastly, we determine our VaR as the negative of the α-quantile, where our α is 5%. We use the int function here because, when using matrix indexation, python can only index using integer values. The int function just makes sure that the script doesn't error out and we use the floor function to round down any floats to integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating portfolio returns\n",
    "port_return = (portval_future - portval_current)/portval_current\n",
    "\n",
    "# Sorting the Returns\n",
    "port_return = np.sort(port_return)\n",
    "\n",
    "# Determining VaR\n",
    "mVar_estimate = -port_return[int (np.floor(alpha*t_simulations))-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b. Applying Historical Simulation to Estimate VaR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import uniform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get Preliminary information. Instead of sourcing historical data here we will generate our own.\n",
    "\n",
    "##### 1. Function that produces a Share Price Path\n",
    "It functions in a similar way to the terminal share value function with the key difference being that this keeps track of the share value on a **daily basis** (through the use of the cumsum function)\n",
    "\n",
    "This is a very convenient function which returns the sum of the elements in a matrix at that point in the matrix (so the first element of the resulting matrix will be the same as the original matrix, but the last element in the cumsum matrix will be the sum of the elements up until that point). This allows us to write share prices as being a function of the individual changes up until that point.\n",
    "\n",
    "**Price Path** can be viewed as a synthetic version of the historical data that we would normally use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Historical simulation\n",
    "np.random.seed(0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_simulations = 10000\n",
    "alpha = 0.05\n",
    "\n",
    "def share_path(S_0, risk_free_rate, sigma, Z, dT):\n",
    "    \"\"\"\n",
    "    Generates the terminal share price given some random normal values, Z\n",
    "    \"\"\"\n",
    "    return S_0*np.exp(np.cumsum((risk_free_rate - sigma**2/2)*dT + sigma*np.sqrt(dT)*Z,1))\n",
    "\n",
    "# generate synthetic share data\n",
    "\n",
    "Z = norm.rvs(size = [3,5*365])\n",
    "corr_Z = np.transpose(np.matmul(L,Z))\n",
    "price_path = share_path(S_0, r, sigma, corr_Z, dT)\n",
    "\n",
    "hist_S_0 = pice_path[-1]\n",
    "hist_portval = np.sum(hist_S_0)\n",
    "hist_portret = [None]*t_simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Credit Valuation Adjustment\n",
    "\n",
    "We now introduce a concept known as exposure. Your exposure at any time 𝑡 is what you are at risk of losing, were the counterparty to default. It is the higher of the value of the portfolio and 0, i.e.\n",
    "\n",
    "$$ E(t) = max(V(t), 0)$$\n",
    "\n",
    "if your portfolio has a negative value (meaning you would owe the other party money), and they default, you don't lose out on any profit. If the value of the portfolio was positive and they defaulted, you would lose the value of the portfolio (the amount they owed you)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Different Types of Exposure:\n",
    "\n",
    "- Current exposure, which is the exposure when t is the current time\n",
    "\n",
    "\n",
    "- Expected Exposure, which is what you predict the exposure to be at a future time t\n",
    "\n",
    "\n",
    "- Potential Future Exposure, the confidence level on exposure at some future time (similar to VaR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of Credit Valuation Adjustments\n",
    "\n",
    "CVA is calculated as the difference between the value of a portfolio which we assume is risk-free, and a portfolio where we account for default risk. As such, CVA can be thought of as the market value of the counterparty credit risk. As you might expect, the higher your exposure, the higher your CVA, since you stand to lose more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let:\n",
    "\n",
    "- $\\delta$ be the recovery rate - the fraction of the portfolio that we would receive if the counterparty defaults\n",
    "\n",
    "\n",
    "\n",
    "- $\\tau$ be a stopping time\n",
    "\n",
    "\n",
    "If our portfolio is worth $V(t)$ at time t and the counterparty defaults at time t, we would:\n",
    "\n",
    "$$ \\begin{cases}\n",
    "               \\text{receive } \\delta V(t)\\\\\n",
    "                \\text{lose } (1-\\delta)V(t)\\\\\n",
    "            \\end{cases} \n",
    "$$\n",
    "\n",
    "\n",
    "N.B. if the counterparty never defaults,$ \\tau = \\infty$\n",
    "\n",
    "Therefore, we only lose out if the counterparty defaults before the time we close-out our position with them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$CVA = E^Q[e^{-r \\tau}(1-\\delta)V(\\tau)I_{\\tau \\leq T}]$$\n",
    "\n",
    "Where:\n",
    "\n",
    "- T is the close-out time for our portfolio \n",
    "\n",
    "\n",
    "- $e^{-r \\tau}$: we are discounting our value to the present time.\n",
    "\n",
    "\n",
    "- $ 1- \\delta $: is the proportion of our portfolio that we lose if there is a default.\n",
    "\n",
    "\n",
    "- $V(\\tau)$: is how much our portfolio is worth in total at the time of default\n",
    "\n",
    "\n",
    "- The Indicator function, I, tells us that we only lose out if the counterparty defaults before we have closed-out our position with them.\n",
    "\n",
    "In order to measure CVA better we have to introduce the Merton Model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merton Model\n",
    "\n",
    "In order to account for the risk of default, we introduce the notion of firm value at time 𝑡. We keep all the assumptions which we had in the Black-Scholes model (besides for the no risk of default) and treat our firm value as an asset. This means it follows geometric Brownian motion and has the following SDE under the risk-neutral measure:\n",
    "\n",
    "$$dV_t = rV_tdt + \\sigma V_t dW_t$$\n",
    "$$ V_0 = V $$\n",
    "\n",
    "Where $V_t$ is the value of the firm at time t. We assume that the value of the firm is the maximum amount which they can pay towards settling their debt.\n",
    "\n",
    "#### Assumption\n",
    "\n",
    "We assume that the total debt which the firm has is some fixed amount, 𝐷, which is due at some future time 𝑇, and that default can only happen at time 𝑇. This means that, at time 𝑇, the debt holders receive:\n",
    "\n",
    "$$min(V_T, D) = D - (D-V_T)^+$$\n",
    "\n",
    "Where the positive-only portion is the same payoff as that of a put on the firm value.\n",
    "\n",
    "**The Expected PV of the Debt is:**\n",
    "\n",
    "$$ De^{-rT} - P(V_0, r \\sigma, D, T) $$\n",
    "\n",
    "Where,\n",
    "\n",
    "$P(V_0, r \\sigma, D, T) $ is the Black-Scholes price for a put with initial asset value $V_0$, risk-free rate, $r$, volatility $\\sigma$, strike $D$ and marturity $T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the closed-form soution for the price of a put in the Black-Scholes framework as follows:\n",
    "\n",
    "$$ V_0 \\Phi(-d_1) + (1 - \\Phi(-d_2))De^{-rT} $$\n",
    "\n",
    "with,\n",
    "\n",
    "$$ d_1 = \\frac{ln(\\frac{V_0}{D}) + (r + \\frac{\\sigma^2}{2})T}{\\sigma \\sqrt{T}} $$\n",
    "$$ d_2 = d_1 - \\sigma \\sqrt{T} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Probability of Default:**\n",
    "\n",
    "is given by:\n",
    "\n",
    "$$ Q[V_T < D] = \\Phi(-d_2) $$\n",
    "\n",
    "These two formulas allow us to give values to the total debt of the firm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "\n",
    "Suppose our portfolio has value $X_t$ at time 𝑡, where $X_t$ is some random variable with maturity time 𝑇. Now, implementing our assumption that default only happens at time 𝑇, our CVA given by equation for the CVA becomes:\n",
    "\n",
    "$$CVA = E^Q[e^{-r \\tau}(1-\\delta)V(\\tau)I_{\\tau = T}]$$\n",
    "$$ = E^Q[e^{-r \\tau}(1-\\delta)V(\\tau)I_{V_T < D}] $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Independence of the Value of the Portfolio with the Value of the Counterparty Firm:** \n",
    "\n",
    "If this is the case we can calculate the expectation as:\n",
    "\n",
    "$$ (1-\\delta)E^Q[e^{-r \\tau}X_T]Q[V_T < D] = (1-\\delta)X_0 \\Phi (-d_2)$$\n",
    "\n",
    "which follows from the fact that discounted assets are martingales under the risk-neutral measure.\n",
    "\n",
    "In the case of correlated portfolio values and firm values, however, we will need to resort to Monte Carlo methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. CVA in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assumptions \n",
    "\n",
    "- Option Type: Call Option\n",
    "\n",
    "\n",
    "- Initial Stock Value: 100\n",
    "\n",
    "\n",
    "- Volatility: 30%\n",
    "\n",
    "\n",
    "- Strike: 110\n",
    "\n",
    "\n",
    "- T: 1 Year\n",
    "\n",
    "\n",
    "- Continuously-compounded risk-free rate: 10%\n",
    "\n",
    "\n",
    "\n",
    "- Counterparty Risk is present\n",
    "\n",
    "\n",
    "\n",
    "- Firm/Counterparty Value: 200\n",
    "\n",
    "\n",
    "\n",
    "- Firm/Counterparty Volatility: 25%\n",
    "\n",
    "\n",
    "\n",
    "- Firm/Counterparty Debt: 180 owed in 1 Year\n",
    "\n",
    "\n",
    "\n",
    "- Value Received if Default Occurs: 20% of what is owed (Recovery Rate)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "import random as r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Market Information\n",
    "risk_free = 0.1\n",
    "\n",
    "# Share Specific Information\n",
    "S_0 = 100\n",
    "sigma = 0.3\n",
    "\n",
    "# Call Option Specific Information\n",
    "strike = 110\n",
    "T = 1\n",
    "\n",
    "# Firm Specific Information\n",
    "V_0 = 200\n",
    "sigma_firm = 0.25\n",
    "debt = 180\n",
    "recovery_rate = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is necessary to create the functions which transform our random standard normals into terminal values. Note that we are going to be modeling our firm value in the same way which we have been modeling stock values – hence, we are going to be using the terminal_value function for both our stock and our firm. The call_payoff function finds the payoff for a call option at terminal time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions for later valuations\n",
    "def terminal_value(S_0, risk_free_rate, sigma, Z, T):\n",
    "    \"\"\"\n",
    "    Generates the terminal share price given some random normal values, Z\n",
    "    \"\"\"\n",
    "    return S_0*np.exp((risk_free_rate-sigma**2/2)*T+sigma*np.sqrt(T)*Z)\n",
    "\n",
    "def call_payoff(S_T, K):\n",
    "    \"\"\"\n",
    "    Function for evaluating the call price in Monte Carlo Estimation\n",
    "    \"\"\"\n",
    "    return np.maximum(S_T-K, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set our seed so that we can easily compare our answers. Line 3 creates an array of the correlations we are going to be testing.This code will have numbers from −1 to 1, going up in steps of 0.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each correlation we test, we will have a CVA estimate and a CVA standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "corr_tested = np.linspace(-1,1,21)\n",
    "cva_estimates = [None]*len(corr_tested)\n",
    "cva_std = [None]*len(corr_tested)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,  0. ,\n",
       "        0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_tested"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following loop we are running our Monte Carlo estimation. \n",
    "\n",
    "We are going to be using a sample size of 50 000 – this means that we have 50 000 estimates for our stock price and 50 000 estimates for our firm value. If our correlation is equal to 1 or −1, our standard normal numbers for each asset (the stock and the firm) are generated as follows:\n",
    "\n",
    "\n",
    "- Generate 50 000 standard normal random variables,\n",
    "\n",
    "\n",
    "\n",
    "- Take the vector of these standard normals and multiply them by the correlation\n",
    "\n",
    "\n",
    "-  Combine these two vectors into a matrix with each array making up a row.\n",
    "\n",
    "\n",
    "\n",
    "If the correlation lies between −0.9 and 0.9 inclusive, we generate the correlated standard normal numbers using a Cholesky decomposition.\n",
    "\n",
    "Once we have our matrix of standard normal random numbers, we can generate our terminal stock and firm values. Line 14 creates an array of stock values using the first row of the matrix of correlated standard normals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Amount Lost\n",
    "\n",
    "$$ e^{-rT}(1-\\delta)X_T I_{V_T < D} $$\n",
    "\n",
    "where,\n",
    "\n",
    "- $X_T$ are the call values\n",
    "\n",
    "\n",
    "- $V$ is the firm value\n",
    "\n",
    "\n",
    "- $D$ is the debt of the firm\n",
    "\n",
    "\n",
    "- $\\delta$ is the recovery rate\n",
    "\n",
    "\n",
    "- 𝑟 is the risk-free rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case a: Correlation is +1 or -1\n",
    "\n",
    "**Correlation = Corr_tested:** \n",
    "\n",
    "Scalar\n",
    "\n",
    "**norm_vec_0**\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "0  \\\\\n",
    "1  \\\\\n",
    ".  \\\\\n",
    ". \\\\\n",
    ". \\\\\n",
    "50 000\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "**norm_vec_1:** Scalar x norm_vector_0\n",
    "\n",
    "$$\\text{scalar} \\times \n",
    "\\begin{bmatrix}\n",
    "0  \\\\\n",
    "1  \\\\\n",
    ".  \\\\\n",
    ". \\\\\n",
    ". \\\\\n",
    "50 000\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "**corr_norm_matrix**: matrix of norm_vec_0 and norm_vec_1\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "0_1 & 0_2 \\\\\n",
    "1_1 & 1_2 \\\\\n",
    ".  & . \\\\\n",
    ". & . \\\\\n",
    ". & . \\\\\n",
    "50 000_1 & 50 000_2 \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "### Case b: Correlation is between +1 and -1\n",
    "\n",
    "**corr_matrix**:\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "1 & scalar \\\\\n",
    "scalar & 1 \\\\\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "**norm_matrix:**\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "0_1 & 0_2 \\\\\n",
    "1_1 & 1_2 \\\\\n",
    ".  & . \\\\\n",
    ". & . \\\\\n",
    ". & . \\\\\n",
    "50 000_1 & 50 000_2 \n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(corr_tested)):\n",
    "    \n",
    "    correlation = corr_tested[i]\n",
    "    if (correlation == 1 or correlation == -1):\n",
    "        norm_vec_0 = norm.rvs(size = 50000)\n",
    "        norm_vec_1 = correlation*norm_vec_0\n",
    "        corr_norm_matrix = np.array([norm_vec_0, norm_vec_1])\n",
    "        \n",
    "    else:\n",
    "        corr_matrix = np.array([[1, correlation], [correlation,1]])\n",
    "        norm_matrix = norm.rvs(size = np.array([2,50000]))\n",
    "        corr_norm_matrix = np.matmul(np.linalg.cholesky(corr_matrix), norm_matrix)\n",
    "     \n",
    "    # create an array of stock values using the first row of the matrix of\n",
    "    # correlated standard normals. (creates an array of 50 000 stock values)\n",
    "    term_stock_val = terminal_value(S_0, risk_free, sigma, corr_norm_matrix[0,], T)\n",
    "    \n",
    "    # create an array of call values for the given stock values.\n",
    "    call_val = call_payoff(term_stock_val, strike)\n",
    "    \n",
    "    # create an array of terminal firm values.\n",
    "    term_firm_val = terminal_value(V_0, risk_free, sigma_firm, corr_norm_matrix[1,], T)\n",
    "    \n",
    "    # Using the call values in the formula\n",
    "    amount_lost = np.exp(-risk_free*T)*(1-recovery_rate)*(term_firm_val < debt)*call_val\n",
    "    cva_estimates[i] = np.mean(amount_lost)\n",
    "    cva_std[i] = np.std(amount_lost)/np.sqrt(50000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Probability of Default\n",
    "\n",
    "We assume our CVA assumes correlations equal to zero with the probability of default given by:\n",
    "\n",
    "$$ \\Phi(-d_2) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to calculate probability of default\n",
    "\n",
    "d_1 = (np.log(V_0/debt) + (risk_free + sigma_firm**2/2)*(T))/(sigma_firm*np.sqrt(T))\n",
    "d_2 = d_1 - sigma_firm*np.sqrt(T)\n",
    "\n",
    "default_prob = norm.cdf(-d_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code for analytical solution for Vanilla European Call Options\n",
    "\n",
    "d_1_stock = (np.log(V_0/debt) + (risk_free + sigma**2/2)*(T))/(sigma_firm*np.sqrt(T))\n",
    "d_2_stock = d_1_stock - sigma*np.sqrt(T)\n",
    "\n",
    "analytic_callprice = S_0*norm.cdf(d_1_stock) - strike*np.exp(-risk_free*(T))*norm.cdf(d_2_stock)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With closed-form analytical solution for the call option on the share being: \n",
    "\n",
    "$$ (1-\\delta)X_0 \\Phi(-d_2)$$\n",
    "\n",
    "Where, $X_0$ is the closed-from solution to price the call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncorr_cva = (1-recovery_rate)*default_prob*analytic_callprice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the Monte Carlo CVA estimates for different correlations, alongside three SD error bounds and the zero CVA we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xm81fP2x/HXOqe5pCKllFOkNNAskQxNIhma0BVNuH6UeciUeZ5dZIp7CZUU10wyFGmgQaIURVKEaK7P74+1z213hs6pzh7O3u/n43EenbP3Z++9+nZa+7s/389nLQshICIiqS8j0QGIiEh8KOGLiKQJJXwRkTShhC8ikiaU8EVE0oQSvohImlDCFxFJE0r4IiJpQglfRCRNlEh0ANH23HPPkJWVlegwRESKjenTp68MIVQtzNikSvhZWVlMmzYt0WGIiBQbZvZ9YcdqSkdEJE0o4YuIpAklfBGRNKGELyKSJpTwRUTShBK+iEiaUMIXEUkTxT/hhwA33QQzZyY6EhGRpFb8E/6qVTBiBHTuDN98k+hoRESSVvFP+FWqwDvv+PcdO8LSpYmNR0QkSRX/hA9Qvz68+Sb8/rsn/ZUrEx2RiEjSSY2ED9C8Obz6KixeDMceC6tXJzoiEZGkkhoJf8lU+OhuqFMGRo/2C7jdu8O6dYmOTEQkaSRVtcydsmQqPHMCbN4AmaWg3wR45hno2xf69IExY6BE8f9riojsqpie4ZvZhWY218zmmNkoMytT5C+y+CNP9mGz/7n4Izj9dHjwQRg/HgYOhC1bivxlRUSKm5glfDOrCVwAtAwhNAYygT5F/kJZ7fzM3jL9z6x2fvv//R/ccIOf7V90ka/XFxFJY7Ge6ygBlDWzjUA54Kcif4VarX0aZ/FHnuxrtd5639VXw2+/wX33wR57wDXXFPnLi4gUFzFL+CGEH83sLuAHYC3wdgjh7ZzjzGwwMBigdu3aO/ditVpvm+i3Pjncfbdvzrr2Wqhc2c/8RUTSUCyndCoD3YE6QA2gvJn1zTkuhDAihNAyhNCyatVCtWXc1vr10KsXvPJK3vdnZMATT/iqnfPPh+ee2/HXEBFJAbG8aNsBWBRCWBFC2Ai8DLQt8lfZuBGWLIHeveG99/IeU6IEvPACHHUU9Ovn6/VFRNJMLBP+D0AbMytnZgYcA8wr8lepUAH++1/fbdu9O3z6ad7jypTxVTvNmvkngkmTijwUEZFkFrOEH0L4DBgDzABmR15rRExerEoVePttqF7dd9nOmpX3uN12gzfegKws6NYNZsyISTgiIskopuvwQwjXhRAahBAahxD+EUJYH7MXq14d3n0XypeHTp1gwYK8x+25pxdbq1zZK2x+/XXMQhIRSSapUVohW1aWJ/PNm6FDh/wrZ+6zj4/LyPA3hx9+iGuYIiKJkFoJH+DAA+Gtt3wpZseOsGJF3uMOOMDH/fEHdOnif4qIpLDUS/jglTNfe80rZ3bunH8yb9rUl3N++61fyN20Ka5hiojEU2omfIB27eDll2HOHDj+eFizJu9xRx0FjzziF30vvDC+MYqIxFHqJnzwFTvPPQeTJ8Mpp8CGDXmPGzjQ6+089BA8/HB8YxQRiZPUrxvcsyf8+acn9b59YdQoyMzcev+SqV6H5/we3hN3yBCoV88v5oqIpJDUT/gAAwb4PP7FF0PFivD4415nJ2ct/TtG+Yqdnj1hyhRo2DDRkYuIFJnUntKJdtFFXi3zySfhkku8XHLOWvorZ8CECVC2rG/MUm9cEUkh6XGGn234cG90fs89UKkSnNnZz+yzz/Cz2kGtfX3lzpFHwskn+3r90qUTHbmIyC5Lr4Rv5rXx//jDyyXvvnvetfTbtIGRI+HUU+Hss+Hpp/2xIiLFWHolfPDdtU8+6RdyhwyB3UdCv4tzj+vTx8suDB/um7kuvzzuoYqIFKX0mcOPll0uuUMH6N8fxo3Le9x113nZ5SuvzL/evohIMZGeCR98Xn7cOGjd2qdupkzJPcbMp3NatfLG6DNnxj9OEZEikr4JH7yW/oQJXkyte3f47rvcY8qW9Tr6e+zhK3d+Kvq2vCIi8ZDeCR+galVvoLJpExx3nBddy6l6de+S9fvv/saQX5kGEZEkpoQP3i1r3DhYuBB69Mi7BMPBB/su3enTvU3ili3xj1NEZBco4Wdr395X77z/Ppxzjm/MyqlbN7jjDhgzxi/oiogUI+m3LHN7/vEPP8sfPhz23x+uuir3mIsv9uWaN90EDRr4xVwRkWJACT+n667z9ojDhkHdur4eP5oZ/OtfPqZ/f6hTB9q2TUysIiI7QFM6OZn51M7hh8OZZ8Inn+QeU6oUjB0LtWvDiSd6AxURkSSnhJ+X0qV9o1Xt2r4qJ6+G6Hvs4at7QvBSysuWxT9OEZEdoISfn+iEftxx8NtvuccccAC8/rr3zVVfXBFJckr421Ovnp/pL17slTPzWq7ZqpW3Upw3zz8NrFsX9zBFRApDCb8g7dp5eYVJk2DQoG2Xay6ZCh/dDQdWgmee8TGnnQabNycuXhGRfGiVTmGcdpov17z2Wl+uec01ubtl9ZsA99/vFTj/+U949FGVVBaRpKKEX1hXX+0Xb6+91pdr1v55225Ziz+CCy6Gn3+GW2+FatXghhsSHbWIyP8o4ReWmffC/eEHX3//3AO5u2UB3Hwz/PIL3HijJ/3zzkts3CIiEUr4OyJ7/X3btnD2VTDuMbDvt+2WZebTOStWwPnne3G2Xr0SG7eICLpou+OqVPHlmhkZMOBKaHjm1mSfLbvBymGHQd++8N57CQlVRCSaEv7O2G8/r5G/ZInvtF27NveYsmW91n79+j5mxoz4xykiEkUJf2e1bQvPPuulF3r2hI0bc4+pXBneess3cR17bN47dkVE4kQJf1f06gWPPOJTPGeemXeN/Bo1POlv2QKdO/sqHhGRBFDC31Vnn+3LMJ9/3i/S5lVHv359L8GwfLlKMIhIwijhF4UrroDLLvOyyddck/eY7BIMc+f6nL5KMIhInCnhF5XbbvPSCzffDHffnfeYTp28BMMHH3jjFJVgEJE40jr8omLm8/l//AGXXAKVKsGAAbnHnXaar9EfOtQ3ZT3yiEowiEhcKOEXpcxM+Pe/4c8/YfBg2H13b4qe05AhPp9/661edvmpp6BChfjHKyJpRVM6RS17N+6hh/rZ/Ntv5z3u5pvhzjt9bJs2WrIpIjEX04RvZpXMbIyZfW1m88zs0Fi+XtIoVw5eew0aNoSTToLJk3OPMfOpn7fe8m5ZLVv68k4RkRiJ9Rn+/cCbIYQGwMHAvBi/XvKoVMmTec2a3jFr1qxt78+upV+/Ikyf7s3Qu3Xzomt5recXEdlFMUv4ZlYROAJ4EiCEsCGE8HusXi8pVasG77zj8/OdOm1tdp5dS//9m/3PzF98x+7pp3v55ZNP9usAIiJFKJZn+HWBFcDTZjbTzJ4ws/I5B5nZYDObZmbTVqxYEcNwEmTffT3pb94MHTvCjz967fyctfTLlfNSDfff79NBrVvD118nOnoRSSGxTPglgObAIyGEZsDfwBU5B4UQRoQQWoYQWlatWjWG4SRQgwbw5pu+IqdjR6jYxGvoW+a2tfTN4IILvLrmb7950n/llcTGLiIpI5YJfymwNITwWeTnMfgbQHpq0QJefRUWLYKBw+Dk5+HoYd4aMWd55fbtfV6/QQO/6HvNNdqkJSK7LGYJP4TwM7DEzOpHbjoG+CpWr1cstG8Po0fDF1/AeTdCy3/mTvbZatWCDz/07lo33eQXdFetim+8IpJSYr1K53zgOTObBTQFbonx6yW/44/38gqTJkHv3rB+ff5jy5SBJ57w3bjvvuv1eObMiV+sIpJSYprwQwhfRObnDwohnBhC0Ckq+Iashx/2KZ7DD4fvv89/rBmcc47X3/n7bzjkEHjppbiFKiKpQzttE+Xcc2HcOPjmG2je3C/qbk/btt41q2lT/2Rw2WWwaVN8YhWRlKCEn0gnnugXZ/fZB7p2heuv3/7F2b33hokT/c3izju9X66SvogUkhJ+ou2/P0yZAmecAcOHe+JfuTL/8aVKed3922+HF1/0TltawSMihaCEnwzKlYOnn4YRI3yuvnlzmDp1+4+57DK45RZ47jlfyaOkLyIFUMJPFmbeQGXyZC+zfPjhfiafV8tE8PIMh5eCiwf7Dt1Bg1SDR0S2S/Xwk02LFj6v/49/eIOUyZPhscegfFRViuxaPJs3QKVSMHQg3PeEv1E89hhk6H1cRHJTZkhGVar4ks0bb/Tm6IccAvPnb70/Zy2ek+rDsGG+Zv+88/L/VCAiaU0JP1llZMDVV3uJ5eXLfdPV2LF+X1a7bWvx1DnC3xwuvxwefdTr8Sjpi0gOmtJJdh07+vr7nj29XeJFF3nD9H4T/Ew/q93W8gy33urLNO++26d37r1X/XJF5H+U8IuD7Lo6l1wC99zjK3hefBHaXbztODNfn795M9x3nyf9u+5S0hcRQAm/+ChVCh54wHvlDhzoSzfHj/f5/Whm/qawaZP/WaKEfyJQ0hdJe5rDL25OPdXP8CtU8Omejz/OPcbM3xzOPRfuuMOvBWhOXyTtKeEXR40aebXNGjWgSxffrJWTGTz0kK/Pv+UWL9sgImlNCb+4qlnTE/2++3o5hnffzT0mI8NX7fTvDzfc4F8ikraU8Iuz6tU96der53X286q4mZEBjz8O/frBddf52b6IpCUl/OKualV4/31o2BC6d/cNWzllZMCTT8Lpp/sGrTvuiH+cIpJwSvipYI89vPH5wQfDySfDyy/nHpOZCSNHQp8+vkHrzjvjHqaIJJYSfqqoXBneecd35Pbq5ev0cypRAv79760NVK6/Xqt3RNKI1uGnkt1391IMxx3nbRQ3bvQmKdFKlPCSyuXKef391au1OUskTSjhp5rddoM33oBu3bypyqZN3iQl2k/Tod+BEHr55qy///ZSzKqyKZLSlPBTUfny8NprcNJJcNZZfqY/aJDfF11auW5JOK8fPPyYJ/2nn/ZPACKSknRKl6rKlfPSC127wuDBfgYP25ZW3rIRejeBm2+G//zH5/7Xr09s3CISMzqdS2VlyviKnd69vU7+hg1wSqS08uYN/mdWO2jX2ks1DBniSztfftnfMEQkpSjhp7rSpWH0aL+Ie+GFsPGOvEsrX3CBJ/1Bg+DYY309f8WKiY1dRIqUEn46KFkSRo3yPy+7DDbc5Buwcurf3+f/+/aFDh18526VKvGPV0RiQgk/XWSvwS9Rwqtnrl3rXbJyLsfs3dunc3r2hCOP9LX91aolJGQRKVq6aJtOMjN9Jc6gQX6h9pxzvFlKTt26wX//CwsXQrt2sGRJ/GMVkSKnhJ9uMjPhscfgqqtgxAhfmbNuXe5xxxzjZ/e//OJJf8GC+McqIkVKCT8dmfkZ/v33+4qcY4+FP/7IPa5tWy/M9tdfcMQR8NVX8Y9VRIrMTiV8MytjZj2LOhiJswsu8DILH3/s8/U//5x7TPPm3k8XPOnPmBHXEEWk6BQ64ZtZppkda2bPAt8DvWMXlsTNaaf5rtxvvoHDDvN5+5waNoSPPvJlm0cdBZ98Ev84RWSXFZjwzewIM3sUWAwMBDoBdUIIPWIcm8RL584+dfPHH570Z87MPWa//TzpV68OnTrl3WFLRJLadhO+mS0FbgM+ARqGEE4B1oYQ1sQjOImjQw7xqZ1SpaB9+7z75Naq5dM7++/vFTknTIh7mCKy8wo6wx8L1MSnb7qZWXlABdRTVYMGMHmyJ/bOnfNupFKtGkycCE2berOVUaPiH6eI7JTtJvwQwhAgC7gHOAr4BqhqZr3MrELsw5O422cfn7pp0cI3X40Yse39S6bC3Kdh5G1w+OHeNvGJJxITq4jskAJ32oYQAvA+8L6ZlQS6AKcC/wL2jG14khBVqvga/J494eyzfS3+sGGw9POtpZUzS8GIl2Dobb6Ra/Vqr9UjIklruwnfzC4BXgwhLAEIIWwEXgVeNbOycYhPEqV8eS+vPGAAXHMNLF8OPWpvLa28eQMs/xxeecVX+lx0ka/Xv/pqdc8SSVIFneHXBCab2SJgFDA6hLASIISwNtbBSYKVLOmNz/faC+6+G37oCM1KgrG1tHKpUvDCCzBwIFx7rZ/p3367kr5IEtpuwg8hXGhmFwFHAH2Aa8zsSzz5jwshrC7oBcwsE5gG/BhCOL4IYpZ4ysjwnrfVqnmlzVWtYdgJ0LDj1tLKJUrAU0/5Ov077/Sk//DDapkokmQKO4c/CZhkZv8HdMCXaj4KFKZLxhBgHqDi6sXZpZfCnnv6fP3VW+D1wdven5EBDz7oSf/22316Ry0TRZLKjuy0bQLcADwMbACuKsRj9gGOA7SMIxWcdRaMGwdz5vgGrUWLtr3fDG67TS0TRZJUQRuv6pnZNWb2FfA8sAboFEI4JIRwXyGe/z7gMmDLrocqSaFbN99lu3KlF1f78svcY666yguzjRsHJ5wAa7RPTyQZFHSGPx84EugdQmgSQrg5hPCdmbUzs/2290AzOx74JYQwvYBxg81smplNW7FixQ4FLwly2GG+Vr9ECS+olteu3AsugCef9OWdXbrAn3/GPUwR2VZBCf8N4MIQwuwct6/Fz9635zDgBDNbDLwAHG1m/8k5KIQwIoTQMoTQsmrVqoUMWxKuUSPflVuzpu/KHTs295j+/X0n7pQpXl//11/jH6eI/E9BCX/fEMKsnDeGEKbhO3DzFUK4MoSwTwghC1/h834Ioe/OBipJqFYtr7+TvSv30Udzj+nd26d2Zs/2Gj3LlsU/ThEBCk74ZbZznzZeie/Kffdd6NoVzj0Xrr8eQo5yS8cf7y0TFy/27llz5iQiUpG0V1DC/9zMBuW80cwGANudm48WQvhAa/BTWLlyfhZ/5pkwfHjevXKPOcbfGP76yytzPvdcQkIVSWcFLZIeCowzs9PZmuBbAqWAk2IZmBQzJUv65qu994Zbb4UVK+D556FM1IfENm281n7v3tC3rzdSufdeKF06cXGLpJGCdtouB9qa2VFA48jN/w0hvB/zyKT4MYNbbvFduUOH+sXc8eOhUiW/f8lUWPyRV9p8ZJzv4J02DUaPhn33TWzsImmgUNsgQwgTgYkxjkVSxZAhnvTPOMOXbb75Jmxeum2lzQsmwKGH+mau5s19iqdLl0RHLpLSVOxEYqNPH79Qu2iRb9D6cOy2lTYXf+QNVKZN86WdXbvCddflnvsXkSKjhC+x07Gjb8paswbOfwyWGVjm1kqbAPXqwaef+qeBG27wxL9yZULDFklVSvgSWy1a+MXZSlXg2TVQ/hToN2FrpU3wVT5PP+3dtSZN8imezz5LXMwiKUoJX2KvXj3flVu/AVzxNDz7Tu6pGzOvxPnJJ5CZ6ev1H3oo95p+EdlpSvgSH9Wr+9l7r17eFatjR/jpp9zjWrSA6dOhUyc4/3zvpvXXX/GPVyQFKeFL/FSs6KtxnnrKp2wOOgheey33uCpVYMIEL7P80kvQujXMmxf/eEVSjBK+xJeZL8WcPh322cfLLQ8dmrtufkaGl1l++22/iNuqlbdSFJGdpoQvidGgga/OueACr53fpg3Mn5973DHH+O7cgw+GU0/1+vpz58Y/XpEUoIQviVOmjCf7CRNgyRKfvx85MveF2po1fXnnLbf4dYCDDoIBA2Dp0kRELVJsKeFL4nXr5p2zWrXy6Z7TT8/dMKVkSeh7DPznYjirF/z7377654or4PffExO3SDGjhC/JoWZNr6Z5441+obZZM5g6dev9S6Z6aYYZ90HWRJj4EpxyijdM328/uOceWLcucfGLFANK+JI8MjN9yeakSbBpk7dSvOMO2LLFSzFEl2bYstAbpc+YAS1bwsUXQ/36fuav8gwieVLCl+Rz2GHwxRfQvTtcfrkXVSvX0Esy5CzN0KwZvPWW987dc08v0dC8uRds06YtkW0o4UtyqlzZyyY/9pg3TO96JmRdBUcPy12aAaBDB/j8c6/Bv3o1HHusr/CZNi0h4YskIyV8SV5mMHiwJ+299oJ/DIH7p8Cf5fMen5HhSze//tpX/8ye7ReC+/SBhQvjG7tIElLCl+TXqJFfwL32Wt+I1aSJl1zIa90+QKlSvr5/4UK/JvDqq77u//rrNb8vaU0JX4qHsmW9X+6iRT6vP348NGzofXTzO3uvWNFX/SxY4Gf5w4d7DZ+ff45r6CLJQglfipc99vCeuYsWeUmGF1/0s/dBg+D77/N+zN57wy3nw5W9YcpkaNoU3nsvvnGLJAElfCme9toL7r7bz+7POQeefdY3Yv3zn/Djj9uOzV7DX/ptGFgBKpbzM3112JI0o4QvxVuNGvDggz5t078/PP64b8QaOnTr1E30Gv49t8BDA6FfP++w1bEjLFuW2L+DSJwo4UtqqFULHn0Uvv3WSzM89BDUrQuXXQa7Nd52Df+BHbzD1siRXqa5aVPf5SuS4pTwJbVkZcGTT3r9/FNOgbvugnY9YPmx0Pqibdfw9+vna/f33NMbrlx7raZ4JKUp4UtqqlfPyyzMneuN0R8aCec8Dt+t3XZcw4a+5PPMM31FT4cOmuKRlKWEL6ntwAN9Jc+nn/r6/KOO8mWd0Q1Xypf3LlwjR3ryb9rUSzWIpBglfEkPhxzijVQGDfKCbG3a5G6kkj3FU7UqdO6sKR5JOUr4kj4qVPDaPOPH+9LNFi3ggQe8Gme27Cmes87yKZ5jjsm72bpIMaSEL+nnhBO8zk6HDjBkiBdai07q5cr5hd9nnvEz/qZN4fXXExevSBFRwpf0VK2a19h55BGvxtmkCYwdu+2YM87wwm3VqsFxx0HPnrk3dYkUI0r4kr7MfJfuzJm+Zr9HD5/KiW6veOCBMO4RGNhlaxG2e+7xBi0ixYwSvkj9+jB5slfWfPZZn8L55BO/b8lUeKEH7PMZ/LMCtD7Iu2u1aOGPESlGlPBFwJuk33gjfPih/3zEEf4GsOCDrWUZdt8C150E48bBqlXemWvgQPj114SGLlJYSvgi0bLbK/brBzffDENHwm8ZW8sy1DkCTjwRvvoKLr3UL+zWr+8XeaNX+4gkISV8kZwqVvSNWGPHwtLlMOIv+PUIOO3lrWUZKlTw9fwzZ/pSzoEDoV07mDUrsbGLbIcSvkh+Tj7Zl28eeRQ8OB6OH+jN0aM1bgyTJvku3W++8QbqF13kfXVFkowSvsj21Kjha/BfeQU2bPA1+127et/cbGY+BTR/vp/p33efr+YZPRpCSFzsIjnELOGbWS0zm2hm88xsrpkNidVricSUGXTv7qUY7rrLV/A0buybtn77beu4KlW8RPPkyd6gpVcvf4NYsCBxsYtEieUZ/ibg4hDCgUAb4DwzaxjD1xOJrdKlfUnmt9/6mfxDD3lVzocego0bt45r08Z36N5/vyf/Jk28O5fq8kiCxSzhhxCWhRBmRL5fDcwDasbq9UTiZq+9/Ex+5kxfs3/++XDwwfDWW1vHlCgBF1wA742CFvvBJZfA4YdvOxUkEmdxmcM3syygGfBZPF5PJC4OOsg7ZWXP73fp4iUYspP6kqnw1gDo8CP0qAhff+VvEHfeqbN9SYiYJ3wzqwCMBYaGEP7M4/7BZjbNzKatWLEi1uGIFK2c8/sff+xTOEOGwKy3fNMWW6BxJjw91Of0L7sM2rb1tfwicRTThG9mJfFk/1wI4eW8xoQQRoQQWoYQWlatWjWW4YjETvT8/oABPq/f9y74fBNsyfBNWy2Og5dfhlGjYOFCaNYMbrtNdXkkbmK5SseAJ4F5IYR7YvU6Ikklen6/eUv472r4d2loMNw3bZlBnz7+ieD44+HKK+HQQ3M3YxGJgVie4R8G/AM42sy+iHx1jeHriSSP6Pn9UAZ6ngP9+8PKlX5/tWowZoy3X1y82Dds3XKLzvYlpmK5SufjEIKFEA4KITSNfKmLhKSP6Pn9yy/3puoNGsDTT/uGLDNfq//VV16fZ9gwX9I5e3aiI5cUpZ22IrFWvrzP1c+c6Qm/f39o337rNE7Vqn6mP3o0/PCDl16+8cZt1/aLFAELSbT1u2XLlmHatGk7/Ljhr87lq59yLQASSTq2ZQtHTvkvp7/8MGXX/s2rnU7n5a5nsqFUGQB2++t3znrhHg6b9i6Lah3AI/2G8f0+9RIctcRawxoVua5bo516rJlNDyG0LMxYneGLxFHIyGDiYd248PpRfNy6Mye9+Sx33dCXpnOmALC6QiUeGHgDd599C5V/X8Ett/RnwPN3UumPlQmOXFJBSpzhixRbH3wA557rm7V69PDCazUjG9JnvQPXXQ+vfQYlS8HQob6Gv1KlREYsSUZn+CLFxZFHesOVm26C117zHroPPACLp8CEf0DTeXB+Jeh8BNx6q/feveMOWLs20ZFLMaSEL5JopUv7Cp05c3xN/pAh0Lkn/LBma2vFizr6Rd82bXzFz/77w+OPaxmn7BAlfJFksd9+3mDlhRdg1Tp4YjW8vh7WZkJWO6/D8/rrPg20774weDA0auSre9ReUQpBCV8kmZhB797w7ULo1xOmb4D7VsNdz8HSpT6mfXuvyT9+vDdf79ULWreGd95JbOyS9JTwRZLR7rvD0y/BvK/htNPgX//y+fuzz4bvvvM3hhNOgC+/9EbqK1dCp07QoYPX4hfJgxK+SDI74ABvqJ7ddGXkSL/tjDN8ZU9mpn8/f76v8Jk1y8/2e/RQ7X3JRQlfpDjIyvKz/EWL/KLu2LHQsKFP/3z5pV/4HTLEq3Bef703Y2nUyHvtfvttoqOXJKGEL1Kc1Kjh7RIXL/ZKm2+84RdzTzgBpk6F3XaD667zaZ8BfeDFUV7O4Ywz4JtvEh29JJgSvkhxVLUq3HwzfP893HCDX8Q95BCfx//wQ1i3CGq/B+eXg0PLwOiXfI1/376a6kljSvgixVnlynDNNX7Gf8cdPr3Tvj0c3wfm/w3ltkDH0vDCFXDRRTBunE8FnX46zJuX6OglzpTwRVLBbrvBpZd64n/wQVjxN/znLxixBuZshiadvJfuokU+bvx4n+M/9VS1WkwjSvgiqaRsWfi//4PFS+COYVBqDxjzJxx1Ktx7r99/++2e+C+7DF59FRo39ou/c+YkOno24lSWAAAPHElEQVSJMSV8kVRUqhRcehMsWuZJvU4dn9KpVQuuuMJr7d92m38iuOIK38HbpIlv4lIDlpSlhC+SyjIyvHfuBx/AZ5/5Rd077/Rlnv37w/Ll3lpx8WK46iov7XDQQb6OX4k/5Sjhi6SL1q3hpZd8Xf7ZZ3uXrcaNoWtX37B1002e+K++Gt56Aw4+GE7srHX8KUQJXyTd1K3rF3Z/+MFbKU6fDkcfDa1awdtvQ/9j4fzycFhpeP1tX845cKCPl2JNDVBE0t26dd5g/e67vURD9crQdB00KwEbMmBxE3h1qo89+2yf+qlePbExy/+oAYqIFF6ZMjBokC/PHD8eatWGN9fCvX/B9E2+zv/bb323bnYRt8svh19/TXTksoOU8EXEZWRESjR8Aa88Aa0awQd/Q9uT4J57PPF//TWcfLJf+K1bF4YPhz//THTkUkhK+CKSW/cB8OFsmDsXevaEhx/2Bi033eRTOrNmwTHHeKG2unX9DWDNmkRHLQVQwheR/DVs6CWZFy6E887z7lqNGvnZ/uWXe+39li19E9d++/kbw/r1iY5a8qGELyIFq13b6+1//71X45w0yfvrXnopXHih/1yvnu/y3b8OXHwKzHkv0VFLDlqlIyI77q+/vIn63XfDjz9Cs2Z+xr9hGVx5Kfy4CUoAXTrAuRf6hq8SJRIddUrakVU6SvgisvM2bIDnnvP6PPPnQ809odnfUNVg7maYnwl/rIFq1bxCZ79+vpNXioyWZYpIfJQqBWed5Rd3x46FParCa2th1FrYaHDX9b67t21b3+x18MHesOXee72sg8SVzvBFpOiEAC8+AiOfginz4c+/oGJFX+7ZuTOsWAHPPw/Tpnk/3i5d/Ky/WzffDyA7TFM6IpJ4GzbAu+/CmDHwyiuwahVUqODJ/ZBDvFTDCy/ATz9BpUpeorlfP78YbJbo6IsNJXwRSS4bN8LEib6sc9w436VbvrwXbqtf33fyTpgAa9dCzT2gS0c4dSC0a+fTRpIvJXwRSV6bNvkyzjFj4OWX4ZdfvDFLu1aw6nNYuQGWbIZNeCevTp3guOP8zaFatURHn3SU8EWkeNi8GT76yJP/qGfht9W+lKRGJtRrDOX2gRkzYNkyH9+qldf3P+44XwqaoXUnSvgiUvwsngI3Hwfz1sDizbBsC2zZAqVLezeu3Xf3xJ/dg3fvvf2s//jjoUMHvz6QhpTwRaR4WjIVFn8EWe2gYn348EOf+584Eb780lcBlS3rZRzM4Lvv4O+/fZ7/yCPh0EawT4Cje0Ldton+28SFEr6IpJ5ff/W5/+w3gLlz/fby5X1u/49V8Osqv60k0KwptO/oq37atIEaNRIWeiwp4YtI6lu+3Hv1TpwI77+/tRWjAWUNSpaF1et8Wgi8gfuhh259A2jWLCXW/idNwjezLsD9QCbwRAjhtu2NV8IXkZ32+WtwZ19Yth5+2QKrK8GyX7beX7q0TwOtW+c/lyjhu37btvU3gNplYfMCqHME1GqdmL/DTkiKhG9mmcA3QEdgKfA5cGoI4av8HqOELyK7JPoaQK3W8NtvMHu21++fNcuvA8yevTXpg78JZOfBDKBCBtRvAk1a+B6BunUhKwvq1IEqVZJuU1iyJPxDgetDCJ0jP18JEEK4Nb/HKOGLSMxt3uz1/aPfBCZPgpV/FPzYkiU96e+9N+y7r5eEbtwYqmXC2m8gqy1kHeKfHkqW9K+MjO2/SeR8k9pBO5LwY1mvtCawJOrnpcAhMXw9EZGCZWbCAQf4V48eftuSqTDieFi1HtZkQrMLYWN5WLzY3xyWLPENYr//7tcOli+HL77YsdfNTvpmW78IEDZ7Jr66GvSbENPppFgm/Lze0nJ9nDCzwcBggNq1a8cwHBGRfNRqDYNfK9yZ9tq1nvB/+slLQr/9DHwxGdZv8QxXsTZU2Ns/SeT82hLZW5D9/ZpVsOZ3X1W0eYO/fjFN+EuBWlE/7wP8lHNQCGEEMAJ8SieG8YiI5K9W68Il27JlfU4/K8sv+HZoBM+c4Ak7sxT0e6HwSXvJ1G0fm9VuV/4GBYplwv8cqGdmdYAfgT7AaTF8PRGR+KvV2qdidmYeflceuxNilvBDCJvM7P+At/BlmU+FEObG6vVERBKmsJ8OivqxOyimTSZDCK8Dr8fyNUREpHBUak5EJE0o4YuIpAklfBGRNKGELyKSJpTwRUTSRFKVRzazFcD3O/nwPYGVRRhOUVFcO0Zx7RjFtWNSMa59QwhVCzMwqRL+rjCzaYUtIBRPimvHKK4do7h2TLrHpSkdEZE0oYQvIpImUinhj0h0APlQXDtGce0YxbVj0jqulJnDFxGR7UulM3wREdmOYpXwzaynmc01sy1mlu8VbTPrYmbzzWyBmV0RdXsdM/vMzL41sxfNrFQRxVXFzN6JPO87ZlY5jzFHmdkXUV/rzOzEyH0jzWxR1H1N4xVXZNzmqNeeEHV7Io9XUzObEvn3nmVmvaPuK9Ljld/vS9T9pSN//wWR45EVdd+Vkdvnm1nnXYljJ+K6yMy+ihyf98xs36j78vw3jVNcZ5rZiqjXHxh1X7/Iv/u3ZtYvznHdGxXTN2b2e9R9MTleZvaUmf1iZnPyud/M7IFIzLPMrHnUfUV/rEIIxeYLOBCoD3wAtMxnTCawEKgLlAK+BBpG7nsJ6BP5/lHg3CKK6w7gisj3VwC3FzC+CvAbUC7y80igRwyOV6HiAv7K5/aEHS/gAKBe5PsawDKgUlEfr+39vkSN+SfwaOT7PsCLke8bRsaXBupEniczjnEdFfU7dG52XNv7N41TXGcCD+Xx2CrAd5E/K0e+rxyvuHKMPx8v2R7r43UE0ByYk8/9XYE38A6BbYDPYnmsitUZfghhXghhfgHDWgMLQgjfhRA2AC8A3c3MgKOBMZFxzwAnFlFo3SPPV9jn7QG8EUJYU0Svn58djet/En28QgjfhBC+jXz/E/ALUKjNJTsoz9+X7cQ7Bjgmcny6Ay+EENaHEBYBCyLPF5e4QggTo36HPsW7ysVaYY5XfjoD74QQfgshrALeAbokKK5TgVFF9Nr5CiF8iJ/c5ac78GxwnwKVzGxvYnSsilXCL6S8mqfXBPYAfg8hbMpxe1GoFkJYBhD5c68Cxvch9y/bzZGPdPeaWek4x1XGzKaZ2afZ00wk0fEys9b4WdvCqJuL6njl9/uS55jI8fgDPz6FeWws44o2AD9TzJbXv2k84zol8u8zxsyyW50mxfGKTH3VAd6PujlWx6sg+cUdk2MV0wYoO8PM3gWq53HXsBDC+MI8RR63he3cvstxFfY5Is+zN9AE7wSW7UrgZzypjQAuB26IY1y1Qwg/mVld4H0zmw38mce4RB2vfwP9QghbIjfv9PHK6yXyuC3n3zMmv1MFKPRzm1lfoCXQPurmXP+mIYSFeT0+BnG9CowKIaw3s3PwT0dHF/KxsYwrWx9gTAhhc9RtsTpeBYnr71bSJfwQQoddfIr8mqevxD8ulYicpeXZVH1n4jKz5Wa2dwhhWSRB/bKdp+oFjAshbIx67mWRb9eb2dPAJfGMKzJlQgjhOzP7AGgGjCXBx8vMKgL/Ba6OfNzNfu6dPl55yO/3Ja8xS82sBLA7/jG9MI+NZVyYWQf8TbR9CGF99u35/JsWRQIrMK4Qwq9RPz4O3B712CNzPPaDIoipUHFF6QOcF31DDI9XQfKLOybHKhWndP7XPN18VUkfYELwKyET8flzgH5AYT4xFMaEyPMV5nlzzR1Gkl72vPmJQJ5X9GMRl5lVzp4SMbM9gcOArxJ9vCL/duPw+c3ROe4ryuOV5+/LduLtAbwfOT4TgD7mq3jqAPWAqbsQyw7FZWbNgMeAE0IIv0Tdnue/aRzj2jvqxxOAeZHv3wI6ReKrDHRi20+6MY0rElt9/CLolKjbYnm8CjIBOCOyWqcN8EfkhCY2xyoWV6Zj9QWchL/zrQeWA29Fbq8BvB41rivwDf4OPSzq9rr4f8gFwGigdBHFtQfwHvBt5M8qkdtbAk9EjcsCfgQycjz+fWA2nrj+A1SIV1xA28hrfxn5c0AyHC+gL7AR+CLqq2ksjldevy/4FNEJke/LRP7+CyLHo27UY4dFHjcfOLaIf98LiuvdyP+D7OMzoaB/0zjFdSswN/L6E4EGUY/tHzmOC4Cz4hlX5OfrgdtyPC5mxws/uVsW+V1eil9rOQc4J3K/AQ9HYp5N1OrDWBwr7bQVEUkTqTilIyIieVDCFxFJE0r4IiJpQglfRCRNKOGLiKQJJXwp1sysupm9YGYLzStHvm5mB8TgdT6w7VRojYwZamblon5+3cwqFXUsIjtLCV+KrcjGq3HAByGE/UIIDYGrgGqFeGxmzucys139/zAU+F/CDyF0DSH8vp3xInGlhC/F2VHAxhDCo9k3hBC+AD42szvNbI6ZzbZILX0zO9LMJprZ88BsM8sys3lm9i9gBlDLzDqZ1+GfYWajzaxCzhc1s0cihbbmmtnwyG0X4BsAJ5rZxMhtiyM7N7Nr18+JfA2N3Jb9+o9HnuttMysb0yMmaU0JX4qzxsD0PG4/GWgKHAx0AO6M2u7fGt+F2TDyc328fEMz4G/gaqBDCKE5MA24KI/nHxZCaAkcBLQ3s4NCCA/gNVCOCiEcFT3YzFoAZwGH4DXPB0XKIoCXY3g4hNAI+B04ZUcPgkhhKeFLKjocr9a4OYSwHJgEtIrcNzV47fps34ethdna4E1NPjGzL/D6OfuSWy8zmwHMBBpFHlNQPONCCH+HEP4CXgbaRe5bFPlUAv7mlVXYv6TIjkq6apkiO2AuW4u7RcurtGy2v7fzs+FNJ07N78GRImmXAK1CCKvMbCRea2d7thfP+qjvNwOa0pGY0Rm+FGfvA6XNbFD2DWbWClgF9DazTDOrireZK0wVy0+Bw8xs/8hzlctjxU9F/E3iDzOrBhwbdd9qYLc8nvdD4MTI85XHiwB+VKi/oUgR0hm+FFshhGBmJwH3mTetXgcsxlfLVMCrHwbgshDCz2bWoIDnW2FmZwKjbGsXravxCozZY740s5n4p4vvgE+inmIE8IaZLYuexw8hzIh8Esh+03kihDDTopqhi8SDqmWKiKQJTemIiKQJJXwRkTShhC8ikiaU8EVE0oQSvohImlDCFxFJE0r4IiJpQglfRCRN/D/dynfm3McvOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(corr_tested, [uncorr_cva]*21)\n",
    "plt.plot(corr_tested, cva_estimates, \".\")\n",
    "plt.plot(corr_tested, cva_estimates + 3*np.array(cva_std), \"r\")\n",
    "plt.plot(corr_tested, cva_estimates - 3*np.array(cva_std), \"r\")\n",
    "plt.xlabel(\"Correlation\")\n",
    "plt.ylabel(\"CVA\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}