Barrier_Options_with_Control_Variates.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Barrier Options with Control Variates"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some notes:\n",
    "\n",
    "1. Barrier Options can be priced analytically with the traditional Black Scholes partial differential equation. We assume that the barrier option is being continously monitored, and we adjust the boundary conditions of the PDE to account for the barrier. Doing so willl give us analytical formulas that are very long and complex.\n",
    "\n",
    "2. The reason why we can price barrier options analytically, but we can't price other path dependent options analytically, such as Asian options, is because the path-dependent nature is not the same. The path dependency of the barrier option is if the value of the stochastic process (the stock price trajectory) hits a constant value, not a time dependent value, or a stochastic value. Meanwhile, Asian options' payoffs are determined by the average value of the stock price on that path. This is entirely random and cannot be modeled as an input that results in an adjustment to boundary conditions or the like.\n",
    "\n",
    "However, it is possible to price exotics in general, including Asian options, byh extending the process to a higher dimensionality, and then applying the Feyman Kac formula in that higher dimensional space to obtain a Partial Differential Equation. Then you would solve that partial differential equation.\n",
    "\n",
    "3. If we think about it, an up and out call, as well as a down and out put are both types of barrier options that we really don't want to be long. If we are in the money for a call, we want to be in the money by as much as possible. If doing so causes us to breach a predefined barrier that casues our option to become worthless, that's not good. The same argument can be made for the put.\n",
    "\n",
    "I'd be curious to see if there are a sellers discount on these barrier options, or a buyer's premium on the down and in put, up and in call.\n",
    "\n",
    "4. Some information regarding hedging the up and out call. at $S_t = B$, the option price has a discontinuity from being worth $ B - K$ to $0$. Thus, for $t$ near $T$ and $S_t$ just below $B$, the delta and gamma are both large negative values. Therefore, the option seller will need to start shorting shares to hedge the up and out call, and then after expiration the trader would need to buy back the shares. Because the delta is large and negative, the trader would need to take a large short position. Because the gamma is large and negative, the trader would need to make large adjustments to the position. Thus, a large amount of trading would be required.\n",
    "\n",
    "Therefore, to prevent this, a common industry practice would to be price and hedge the up and out call as if the barrier were at a slightly higher level than B. That way, the regime where delta and gamma are large and negative occurs in the region above the contractual barrier B, and therefore this involved trading is no longer necessary."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References:\n",
    "\n",
    "https://bookdown.org/maxime_debellefroid/MyBook/barrier-options.html\n",
    "\n",
    "Shreve, Stochastic calculus for finance, \n",
    "\n",
    "Hull, Options, Futures, and other Derivatives."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The payoff of an up and out, continuously monitored barrier call is:\n",
    "\n",
    "$ V(t, S_t) = \\mathbb{E}^\\mathbb{Q} \\left[ e^{-r(T-t)}\\mathbb{I}_{\\tau(b)> T}\\max(S_T - K, 0) \\right | \\mathcal{F}_t] $\n",
    "\n",
    "where $\\tau(b) = \\inf\\{t: S(t) \\geq b\\}$.\n",
    "\n",
    "\"Inf\" in this case refers to infimum, that is the greatest element that is lower or equal to the lowest element of the superset that the quantity in question belongs to. In simple English, $\\tau(b)$ refers to the first time the process $S(t)$ is greater than or equal to the barrier. Since this is an up and out call, the indicator function in the payoff means that the first time the process exceeds the barrier occurs after the maturity. That is to say, at and prior to maturity, the barrier is not breached, which we want in order for this option payoff to be worth something."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part 0: The Monte Carlo estimator is 0.6390 with 95 percent CI is [0.5985, 0.6796]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as sp\n",
    "\n",
    "\n",
    "# Set seed\n",
    "np.random.seed(seed=234513)\n",
    "\n",
    "\n",
    "# GBM parameters\n",
    "r = 0.02\n",
    "sigma = 0.2\n",
    "BM_drift = 0 # = r - sigma^2/2;\n",
    "S0 = 50\n",
    "\n",
    "# Generate 10000 sample paths\n",
    "reps = 10000 # of replications (sample paths)\n",
    "m = 2**10 # of time points t_i, not including t=0\n",
    "#Sample from a N(0, 1) with dimension/shape (reps, m)\n",
    "Z  = sp.norm.rvs(loc = 0, scale = 1, size = (reps, m))\n",
    "#Add a reps x 1 column to the beginning of the matrix for the t = 0 point\n",
    "Z = np.concatenate((np.zeros(shape = (reps, 1)), Z), axis = 1)\n",
    "\n",
    "coef = np.sqrt(1/m) # sqrt(t_i - t_{i-1})\n",
    "W = np.cumsum(sigma * coef * Z , axis = 1) # cumulative sum across rows (over time)\n",
    "S = S0*np.exp(W)\n",
    "\n",
    "# Other parameters\n",
    "b = 75 # barrier\n",
    "K = 60 # strike\n",
    "T = 1 # maturity\n",
    "\n",
    "# Calculate discounted payoffs\n",
    "\n",
    "#[maxS, ~] = max(S, [], 2) # maximum for each sample path \n",
    "#disc_payoffs = np.exp(-r*T)*(maxS < b).*max(S(:,m+1) - K, zeros(reps, 1));\n",
    "\n",
    "#compute maximum value attained across every simulation\n",
    "maxS = np.max(S, axis = 1)\n",
    "disc_payoffs = np.exp(-r * T) * (maxS < b) * np.maximum(S[:,-1] - K, 0)\n",
    "\n",
    "# Monte Carlo estimator\n",
    "price_MC = np.mean(disc_payoffs)\n",
    "sd_MC = np.std(disc_payoffs)\n",
    "z_alpha2 = 1.96\n",
    "\n",
    "CI_MC = [price_MC - z_alpha2*sd_MC/np.sqrt(reps), price_MC + z_alpha2*sd_MC/np.sqrt(reps)]\n",
    "print(f'Part 0: The Monte Carlo estimator is {price_MC:.4f} with 95 percent CI is [{CI_MC[0]:.4f}, {CI_MC[1]:.4f}]')\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Control Variates\n",
    "\n",
    "Whenever estimating a random variable on its own, the variance is determined by aspects such as sample size and sampling methodology. When estimating a function of a random variable, the function itself introduces more variance. However, if the behavior of the function is known, then those quantities can be \"filtered\" out to reduce the variance of the statistic that requires estimation.\n",
    "\n",
    "In other words, if the estimate contains some known quantities, then any error associated with that known quantity can be removed to the reduce the error of an estimate of the unknown quantity.\n",
    "\n",
    "Formal definition:\n",
    "\n",
    "Let $(X_1, Y_1), ... , (X_n, Y_n)$ be i.i.d. outputs from $n$ replications of a simulation. Suppose the goal is to estimate $\\mathbb{E}[Y_1]$. The standard Monte Carlo estimator is $\\bar{Y} = \\frac{Y_1 + Y_2 + ... + Y_{n-1} + Y_{n}}{n}$ and it is unbiased and consistent.\n",
    "\n",
    "As a quick review, unbiased means the expected value of the estimator ($\\bar{Y}$) is equal to the population parameter, i.e. has no offset. Consistent means as the sample size $n$ approaches infinity (if we take the limit), the estimator converges to the population parameter.\n",
    "\n",
    "Suppose that $\\mathbb{E}[X_1]$ is known (and probably easier, less difficult to compute). Then we can define for any fixed constant $\\lambda$:\n",
    "\n",
    "$Y_i(\\lambda) = Y_i - \\lambda(X_i - \\mathbb{E}[X_1])$\n",
    "\n",
    "We can then define the control variate estimator as\n",
    "\n",
    "$\\bar{Y}(\\lambda) = \\frac{1}{n} \\sum\\limits_{i=1}^{n} Y_i(\\lambda)$\n",
    "\n",
    "It can be shown that this estimator is consistent and unbiased. Using the properties of variance, the variance can be expanded as:\n",
    "\n",
    "$Var(\\bar{Y}(\\lambda)) = \\frac{1}{n} \\left[ Var(Y1) - 2\\lambda Cov(X_1, Y_1) + \\lambda^2 Var(X_1) \\right]$\n",
    "\n",
    "This function is a convex quadratic for $\\lambda$. Thus we can take the derivative with respect to $\\lambda$ and set it equal to zero to find the global minimum. Doing so yields\n",
    "\n",
    "$\\lambda^{*} = \\frac{Cov(X_1, Y_1)}{Var(X_1)}$\n",
    "\n",
    "This value of lambda minimizes the variance.\n",
    "\n",
    "More computation will show that the mimimum variance $\\sigma^2(\\lambda^{*}) = (1 - \\rho_{XY}^2)Var(Y_1)$\n",
    "\n",
    "One last step is that the correlation coefficient itself is typically unknown, so $\\lambda^{*}$ needs to be estimated. This can be done across the Monte Carlo replicates, and this technique is referred to as regression sampling. The estimator is:\n",
    "\n",
    "$\\lambda_n = \\frac{\\sum_{i=1}^{n} (X_i - \\bar{X}) (Y_i - \\bar{Y})}{\\sum_{i=1}^{n}(X_i - \\bar{X})^2}$\n",
    "\n",
    "which if you look closely is simply the ratio of the sample covariance to the sample variance (the $\\frac{1}{n-1}$ terms have been ommitted from both the numerator and denominator)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the question remains, what do we choose as a control variate? As mentioned, it has to be a quantity such that its expectation can be computed. Also, it needs to be correlated to some degree with the statistic we are trying to estimate (the barrier option payoff).\n",
    "\n",
    "We will examine three candidates\n",
    "\n",
    "1. The price of the underlying at the time of expiration $S(T)$\n",
    "2. The discounted payoff of a European call, with the same exact strike, maturity, and dynamics of the underlying, i.e. $e^{-rT}[S(T) - K]^+$\n",
    "3. An indicator for whether or not the barrier option will be positive: $\\mathbb{I}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\} $\n",
    "\n",
    "To further see the results of variance reduction, we will maintain the random seed so we are in effect using common random numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Control Variate S(T) estimator is 0.6410 with 95 percent CI is [0.6057, 0.6763]\n",
      "The estimated variance reduction for CV S(T) is 24.23%.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "#Estimator a: S(T)\n",
    "Xa = S[:,-1]\n",
    "EXa = np.exp(r*T)*S0\n",
    "lambdana = np.mean((Xa - np.mean(Xa)) * (disc_payoffs - price_MC)) / np.mean((Xa - np.mean(Xa))**2)\n",
    "Ylambdaa = disc_payoffs - lambdana*(Xa - EXa)\n",
    "price_CVa = np.mean(Ylambdaa)\n",
    "sd_CVa = np.std(Ylambdaa)\n",
    "\n",
    "CI_CVa = [price_CVa - z_alpha2*sd_CVa/np.sqrt(reps), price_CVa + z_alpha2*sd_CVa/np.sqrt(reps)]\n",
    "var_red_a = 1 - sd_CVa**2/sd_MC**2\n",
    "\n",
    "print(f'The Control Variate S(T) estimator is {price_CVa:.4f} with 95 percent CI is [{CI_CVa[0]:.4f}, {CI_CVa[1]:.4f}]')\n",
    "print(f'The estimated variance reduction for CV S(T) is {var_red_a:.2%}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Control Variate e^(-rT)*[S(T) - K]^+ estimator is 0.6337 with 95 percent CI is [0.5980, 0.6695]\n",
      "The estimated variance reduction for CV S(T) is 22.31%.\n"
     ]
    }
   ],
   "source": [
    "# (b) Control variate of e^(-rT)*[S(T) - K]^+\n",
    "import BMS\n",
    "\n",
    "Xb = np.exp(-r*T)*np.maximum(S[:,-1] - K, 0)\n",
    "# Calculate Black-Scholes price\n",
    "#temp = (np.log(S0/K) + (r - sigma**2 / 2) * T) / (sigma * np.sqrt(T))\n",
    "#EXb = S0 * normcdf(temp + sigma*sqrt(T)) - K * np.exp(-r * T) * normcdf(temp)\n",
    "\n",
    "EXb = BMS.BMS_price(type_option = 'call', S=S0, K=K, r=r, q=0, sigma=sigma, T=T, t=0)\n",
    "lambdanb = np.mean((Xb - np.mean(Xb)) *(disc_payoffs - price_MC))/np.mean((Xb - np.mean(Xb)) **2);\n",
    "Ylambdab = disc_payoffs - lambdanb*(Xb - EXb)\n",
    "price_CVb = np.mean(Ylambdab)\n",
    "sd_CVb = np.std(Ylambdab)\n",
    "CI_CVb = [price_CVb - z_alpha2*sd_CVb/np.sqrt(reps), price_CVb + z_alpha2*sd_CVb/np.sqrt(reps)]\n",
    "var_red_b = 1 - sd_CVb**2/sd_MC**2\n",
    "\n",
    "print(f'The Control Variate e^(-rT)*[S(T) - K]^+ estimator is {price_CVb:.4f} with 95 percent CI is [{CI_CVb[0]:.4f}, {CI_CVb[1]:.4f}]')\n",
    "print(f'The estimated variance reduction for CV S(T) is {var_red_b:.2%}.')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Control Variate 3, we need to do compute the expected value of our control variate. That is, we need to compute:\n",
    "\n",
    "\n",
    "$\\mathbb{E} \\left[ \\mathbb{I}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\} \\right]$\n",
    "\n",
    "An expectation of an indicator variable is simply the probability. Thus we need to find:\n",
    "\n",
    "$\\mathbb{P}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\}$\n",
    "\n",
    "We will solve this using the facts about hitting times for a Brownian Motion.\n",
    "\n",
    "For a BM $W$, it is true that:\n",
    "\n",
    "$[1] : \\mathbb{P}\\{ \\max_{0\\leq t \\leq T} W(t) \\geq a = 2 \\mathbb{P}(W(T) \\geq a) $\n",
    "\n",
    "and\n",
    "$ [2] : \\mathbb{P}\\{ \\max_{0\\leq t \\leq T} W(t) \\geq a = 2, W(T) \\leq c \\} = \\mathbb{P}(W(T) \\geq 2a- c) $\n",
    "\n",
    "for arbitrary values $0 < c < a $\n",
    "\n",
    "Since the events are disjoint, we can say:\n",
    "\n",
    "$\\mathbb{P}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\} = 1 - \\mathbb{P}(\\max_{0\\leq t \\leq T} S(t) \\geq b ) - \\mathbb{P}(S(T) \\leq K) + \\mathbb{P}(S(T) \\leq K, \\max_{0\\leq t \\leq T} S(t) \\geq b )$\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use the definition of GBM:\n",
    "\n",
    "$S(t) = S(0)\\exp((r-\\sigma^2/2)t - \\sigma W(t)) = 50\\exp(0.2*W(t))$\n",
    "\n",
    "Now we can approach this term by term. The first term:\n",
    "\n",
    "$\\mathbb{P}(\\max_{0\\leq t \\leq T} S(t) \\geq b ) = \\mathbb{P}(\\max_{0\\leq t \\leq T} 50\\exp(0.2*W(t)) \\geq 75 )$\n",
    "\n",
    "$  = \\mathbb{P}(\\max_{0\\leq t \\leq T} *W(t) \\geq 5 \\ln (75/50) ) $\n",
    "\n",
    "where $ln$ denotes the natural logarithm (log base e)\n",
    "\n",
    "$  = 2\\mathbb{P}(\\max_{0\\leq t \\leq T} *W(1) \\geq 5 \\ln (75/50) ) $ since $T = 1$\n",
    "\n",
    "and finally we can then say:\n",
    "\n",
    "$  = 2 (1 - \\Phi(5\\ln(75/50))) = 0.0426 $\n",
    "\n",
    "Since the normal CDF is a known function.\n",
    "\n",
    "\n",
    "\n",
    "For term 2 we have:\n",
    "\n",
    "$ \\mathbb{P}(S(T) \\leq K) = \\mathbb{P}(50 \\exp (0.2W(1)) \\leq 60) = \\mathbb{P}(W(1) \\leq 5 \\ln(60/50)) = \\Phi[5 \\ln(60/50)] = 0.8190.$\n",
    "\n",
    "Finally, for term 3, we have:\n",
    "\n",
    "$ \\mathbb{P}(S(T) \\leq K, \\max_{0\\leq t \\leq T} S(t) \\geq b ) = \\mathbb{P} \\left( \\max_{0\\leq t \\leq T} 50\\exp(0.2W(t)) \\geq 75, 50 \\exp(0.2W(T)) \\leq 60 \\right) $\n",
    "\n",
    "$ = \\mathbb{P} \\left( \\max_{0\\leq t \\leq T} W(t) \\geq 5 \\ln(75/50), W(T) \\leq 5 \\ln(60/50) \\right)$\n",
    "\n",
    "$ = \\mathbb{P}( W(1) \\geq 10 \\ln (75/50) - 5 \\ln (60/50))$\n",
    "$ = 1 - \\Phi(10 \\ln (75/50) - 5 \\ln (60/50))$ = 0.0008\n",
    "\n",
    "Note that the joint probability is the lowest, which intuitively makes sense. It is highly unlikely for the joint event to occur as it is a joint event.\n",
    "\n",
    "\n",
    "All together, we have, when combining the values of the terms:\n",
    "\n",
    "$\\mathbb{P}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\} = 1 - \\mathbb{P}(\\max_{0\\leq t \\leq T} S(t) \\geq b ) - \\mathbb{P}(S(T) \\leq K) + \\mathbb{P}(S(T) \\leq K, \\max_{0\\leq t \\leq T} S(t) \\geq b )$\n",
    "\n",
    "$\\mathbb{P}\\{ \\max_{0\\leq t \\leq T} S(t) < b , S(T) > K\\} = 1 - 0.0426 - 0.8190 + 0.00008 = 0.1392 $\n",
    "\n",
    "This value is used below in the code.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Control Variate estimator for the joint probability is 0.6517 with 95 percent CI is [0.6261, 0.6772]\n",
      "The estimated variance reduction for CV S(T) is 60.34%.\n"
     ]
    }
   ],
   "source": [
    "# (c) Control variate of 1{max S(t) < b and S(T) > K}\n",
    "Xc = (maxS < b) *(S[:,-1] > K)\n",
    "EXc = 0.1392# This computation is done above.\n",
    "lambdanc = np.mean((Xc - np.mean(Xc)) *(disc_payoffs - price_MC))/np.mean((Xc - np.mean(Xc))**2)\n",
    "Ylambdac = disc_payoffs - lambdanc*(Xc - EXc)\n",
    "price_CVc = np.mean(Ylambdac)\n",
    "sd_CVc = np.std(Ylambdac)\n",
    "CI_CVc = [price_CVc - z_alpha2*sd_CVc/np.sqrt(reps), price_CVc + z_alpha2*sd_CVc/np.sqrt(reps)]\n",
    "var_red_c = 1 - sd_CVc**2/sd_MC**2\n",
    "\n",
    "print(f'The Control Variate estimator for the joint probability is {price_CVc:.4f} with 95 percent CI is [{CI_CVc[0]:.4f}, {CI_CVc[1]:.4f}]')\n",
    "print(f'The estimated variance reduction for CV S(T) is {var_red_c:.2%}.')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A couple of points. First, we achieved the largest reduction in variance using the last control variate. This makes sense because it is the most strongly correlated with the discounted payoff the barrier option. How is this the case? Whenever CVc is zero, the discounted payoff is zero. Whenever CVc is positive, the discounted payoff is positive. While the CV b will equal the discounted payoff more option, we need to look at scatter plots to see that their relationship is highly nonlinear, and therefore the correlation lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvoElEQVR4nO3deZhU5Zn38e/dTUM3iCKC7JvEwRghqB0lUScxLpiogFlEo1EzRjJx8gaXOAFjFBOMRB3FSSYzE0ZHMjGKidohMo5bXAYVDQgBFxwFRGlQWQQVW2jgfv84p6C6qH1ffp/rqqurzvKcu053113Pcp5j7o6IiEgm6kodgIiIVB4lDxERyZiSh4iIZEzJQ0REMqbkISIiGVPyEBGRjCl5SM0zs+lmtsHM3i51LJkysy+Y2ZoiHWuwmX1oZvXFOJ6UNyUPSYuZHWtmz5jZFjPbZGZPm9lncizzAjObH7PsDjObnlu0GcUwCLgcONTd+ybYZl8zm2lmb4Yfnq+Hr3vl4fhuZp/ItZwEZS83s7+Ls3yymS3MtDx3f9Pd93H3nWkcu2hJTUpDyUNSMrN9gQeAXwA9gQHAtcC2UsYVj5l1ynCXIcBGd383QXmdgceATwGnAPsCnwM2AkflEGpasng/0WYD58VZ/s1wXbHikGrk7nrokfQBNAObU2xzEfAK8AHwMnBEuHwKsCJq+Rnh8k8CHwM7gQ+BzcAkoB3YHi77U7htf+BeYD2wCvh+1HGnAX8Afgu8D3w7Tmz7Ab8J918NXEXwxelEoA3YFR7vjjj7fht4B9gnyXv/JPBE+B5eAsZFrbsD+BdgXngOngOGh+ueAhzYGh5/IvAFYA3wQ+Bt4L+ALsBMYG34mAl0Ccv4ArAmQVwDgR3AkJhYtwO9gFOBxeF5ewuYFrXd0DC2C4E3w1gjyzqF23wr6ne+EvhOuLxbzHn9MPwd1kX9PWwE7gF6hvs0hr/DjeF5/AvQp9R/+3ok+Z8vdQB6lP+D4Nv2RoJvq18C9o9Z/3WgFfgMYMAnIh9Y4brIB8fE8IOyX7juAmB+TFl3ANOjXtcBi4Crgc7AQeEH1dhw/TSChDMh3LYpTvy/Af4IdA8/AP8PuDBcl/DDN1x/NzA7yfoG4HXgyjC+L4YfpiOi3s8mglpKJ+BO4O6o/R34RNTrLxB84P+cIGk0AT8BFgAHAr2BZ4Cfphn/I8BVUa+vB1qi9h0ZnrdRBElyQrhuaBjbbwiSQRN7J49TgeHh7/zzwEfs+dKwV1zAJeH7GBi+t38H7grXfQf4E9AVqAeOBPYt9d++HokfJQ9Aj8p4EHxjvYPgW/EOYC7hN0PgIWBymuUsAcaHzy8gdfI4GngzZpupwH+Gz6cBTyU5Xj1B89qhUcu+AzwRPk/nw3dGkvXHEdQQ6qKW3UX4LT58P/8Rte7LwPKo1/GSx3agMWrZCuDLUa/HAm+kGf+5wKvh8zqCWsQZCbadCdwSPo8kioOi1ndIHnH2b4n8HSRIHq8AJ0S97keQ+DsBf0eQFEeV+m9dj/Qe6vOQtLj7K+5+gbsPBA4jqE3MDFcPIviA24uZnWdmS8xss5ltDvfNpKN5CNA/sn9YxpVAn6ht3kqyfy+CGsHqqGWrCfpt0rGR4EMukf7AW+6+K0n50aO4PgL2SXHM9e7+ccwxYuPvn6KMiPuAfmY2huADvStBExpmdrSZPW5m681sC/D37P27SXhuzexLZrYgHECxmSAxJvvdDgHuj/o9vkLQbNmHoHnuIeBuM1trZjeYWUOa71FKQMlDMubuywm+UR8WLnqLoPmiAzMbAswCvgcc4O49gBcJmjkg+Ba7V/Exr98CVrl7j6hHd3f/cpJ9om0g+HY7JGrZYIJmtnQ8Cow1s24J1q8FBplZ9P9SJuXHE/t+1rJ3/GvTKsj9I4I+ofMIOsrvdvft4erfEdQgB7n7fsC/sed3kygWAMysC0E/1E0ENdAewH+T/Hf7FvClmN9lo7u3unu7u1/r7ocSDEg4jfid/VImlDwkJTM7xMwuN7OB4etBwNkE7dcA/wH8wMyOtMAnwsTRjeBDZH2437fYk3AgaGMfGI5oil52UNTr54H3zeyHZtZkZvVmdli6w4Q9GFZ6D3CdmXUP47qMoHM2Hf9F8KF3b3ge6szsADO70sy+TNABvhX4RzNrMLMvAKcT9JWkI/b9xnMXcJWZ9Q6HB1+dQfwQ9FVNBL5Kx1FW3YFN7v6xmR0FfCODMjsT9FusB3aY2ZeAk6PWvwMcYGb7RS37N4LfwxCA8P2MD58fb2Yjw2tI3idI+CmHBEvpKHlIOj4g6Ht4zsy2EiSNFwmuj8Ddfw9cR/BN9gOCtu+e7v4y8E/AswQfJiOBp6PK/TPB6KS3zWxDuOw24NCwaaMl/PA/HRhNMNJqA0Gyiv5QSuX/EXzArwTmh3Hens6O7r6NYFTWcoL+j/cJElov4LnwW/w4goEEG4BfAeeFtbN0TANmh+/3zATbTAcWAkuBZcAL4bJ0PQVsAVrd/S9Ryy8GfmJmHxAkpHvSLdDdPwC+H+7zHkHimRu1fjlB0lsZvrf+wK3hNg+Hx1xA8HcF0JeghvQ+QXPWk2SWIKXIzF03gxIRkcyo5iEiIhlT8hARkYwpeYiISMaUPEREJGNVMdlZr169fOjQoaUOQ0SkoixatGiDu/fOZt+qSB5Dhw5l4cKMZ5gWEalpZrY69VbxqdlKREQypuQhIiIZU/IQEZGMVUWfh4hIIbW3t7NmzRo+/vjj1BuXocbGRgYOHEhDQ/4mKlbyEBFJYc2aNXTv3p2hQ4diFjvxcHlzdzZu3MiaNWsYNmxY3spV8hApYy2LW7nxoVdZu7mN/j2auGLsCCYcnu6tSCRfPv7444pMHABmxgEHHMD69evzWq6Sh0iZalncyqVzluy+MUbr5jYunbMEQAmkBCoxcUQUInZ1mIuUqcvvWbLXHZUc+Mc//LUU4Yh0oOQhUoZaFreyM8HdErYnWiFV66233mLYsGFs2rQJgPfee49hw4axenXW1/jlTMlDpAxNm/tSqUOQMjJo0CC++93vMmXKFACmTJnCpEmTGDJkSIo9C0d9HiJlaHNbe6lDkBwUYqDDpZdeypFHHsnMmTOZP38+v/jFL/jwww8ZP3487733Hu3t7UyfPp3x48fn6V0kp+QhUmaOvu6RpOs711dux20taFncytT7ltHWHtyCvXVzG1PvWwbkNtChoaGBG2+8kVNOOYWHH36Yzp07U1dXx/3338++++7Lhg0bGDNmDOPGjStK576arUTKSMviVt75YHvSbW742qeLFI1k48aHXt2dOCLa2ndy40Ov5lz2gw8+SL9+/XjxxReB4BqOK6+8klGjRnHiiSfS2trKO++8k/Nx0qGah0gZ+eG9Szu8frDzFRxirbtfv+oDOeRw9YeUs7Wb2zJanq4lS5bwyCOPsGDBAo499ljOOussHnroIdavX8+iRYtoaGhg6NChRbsKXjUPkTJxVcsytu3Ytfv1050v5hBrxYzdjxF1a+CXR5cwSkmlf4+mjJanw9357ne/y8yZMxk8eDBXXHEFP/jBD9iyZQsHHnggDQ0NPP7440UdfaXkIVImfrvgzd3Px9XNp79tJrbp2gA2LC9qXJKZK8aOoKmhvsOypoZ6rhg7IusyZ82axeDBgznppJMAuPjii1m+fDmjR49m4cKFNDc3c+edd3LIIYfkFHsm1GwlUgbOmfVsh9c/a7htr8QhlSHSKZ7P0VaTJk1i0qRJu1/X19ezaNEiAJ599tlEuxWUkodIGXh6xabdz6/tdDvd2FbCaCRXEw4fUPVTyCh5iJTYVS3Ldj8fVzef8+of3V3rmNetK7fu34O3O9XTd8dOJr+3mVM/UmKR0lPyECmx6L6OGxpmdUgcU3sfgIcL1jV0YmrvA2DoGZxaikBFoqjDXKSEWhbvGYZ7bafb6cKeK8uv7NVzd+KIcDN+/Oa8osUnkoiSh0gJXRZOsQ7wzZjmql0JeszbXVOXSOkpeYiUSMviViJXdTzd+WKiU8X1PfdHw62knJU0eZjZ7Wb2rpm9GLWsp5k9YmavhT/3L2WMIoVySVStI/aaji31+l4neySakn327NmcdtppJYmp1H+hdwCnxCybAjzm7gcDj4WvRarWC50v7PB6es8eSbdvsIYCRiPlqBynZC9p8nD3p4BNMYvHA7PD57OBCcWMSaQYIsNzx9XNZ39r61DrmLNv96RNVj899qeFDk9ytfQeuOUwmNYj+Ln0npyLvPTSS1mwYMHuKdkvv/xyAN5//33OOOMMDj30UP7+7/+eXbt2pSgpP8pxqG4fd18H4O7rzOzAeBuZ2SRgEsDgwYOLGJ5I7iLDc2c2/KpDnpjXrWvS/Zrqmzj1IA3ULWtL74E/fR/aw4kQt7wVvAYYdWbWxcabkh3g+eef5+WXX2bIkCGccsop3HfffXzta1/L9V2kVOpmq6y5+6/dvdndm3v37l3qcETSNuqa/wGCGXNj6xdX9+qZtNZxzeeuKWBkkheP/WRP4ohobwuW5yh2SnaAo446ioMOOoj6+nrOPvts5s+fn/Nx0lGOyeMdM+sHEP58t8TxiOTV+9uCez1EZsyNmN6zB9tTjLBSraMCbFmT2fI0RU/Jfsstt7Bu3TqAvW78VIwbQUF5Jo+5wPnh8/OBP5YwFpG8ikyAeG2n2/dal6qvY/i+wwsWl+TRfgMzW56GRFOyQ9BstWrVKnbt2sWcOXM49thjsz5OJko9VPcu4FlghJmtMbMLgRnASWb2GnBS+FqkKkQmQIyevwrgoj69Uu7bckZLgaKSvDrhamiIuXdHQ1OwPEuJpmR/8skn+exnP8uUKVM47LDDGDZsGGeccUYu0aetpB3m7n52glUnFDUQkSKITEXyf52/sde6BU1NSWsdM47Td6iKEekUf+wnQVPVfgODxJFDZ3myKdlLpRxHW4lUpUvmLGFc3XwaLPOLx9XXUWFGnZlTsqgE5djnIVJ1In0dt8QMzQWY0L9PCSISyY2Sh0gRPL1iE+Pq5sf9h1vRuXPSqsiYvmMKF5hIlpQ8RArs6OseAfa+IBDS6yifNXZWIcISyYmSh0iBvfPB9r1mzY1QR7lUKiUPkQI66eYngL1nzYXUEyCCOsqlfCl5iBTQa+9ujTs0F1JfFHhgY9xp3aQGJZqSffXq1SWLSclDpEDOmfUsv2m4LquhuQCPTXws/0FJRdKU7CI15OkVmziu7qW4iSOdJiupXPNWzuPkP5zMqNmjOPkPJzNvZe73nU80JfsNN9zAyJEj+fSnP707uRSDLhIUKYCrWpbxYOcrEq5P1WTVvVP3QoQlRTBv5TymPTONj3d+DMC6reuY9sw0ILc+rHhTsj/44IO0tLTw3HPP0bVr193NWsWgmodIAfx2wZt7zZobkc7w3GfOeaYAUUkx3PrCrbsTR8THOz/m1hduzbns2CnZH330Ub71rW/RtWtwH5iePXvmfIx0KXmI5NlJNz+RtNaRaniuLgqsbG9vfTuj5emKNyW7uxdtCvZYSh4iefbau1sT1jrS6evQRYGVrW+3vhktT0eiKdlPPvlkbr/9dj766CMANVuJVKqWxa2sTDA0F1L3dTTVNyVcJ5Vh8hGTaaxv7LCssb6RyUdMzrrMRFOyNzU1MW7cOJqbmxk9ejQ33XRTTrFnQh3mInn059//kvEN2Q3NBd1mthpEOsVvfeFW3t76Nn279WXyEZNz6ixPNiX75z//+aKOsopQ8hDJk3NmPct/xZm/KmJet64py9AV5dXh1INOrfrfpZqtRPLk3986Pe78VRFTevXUFeVSNZQ8RPKgZXEr3aw9YW64qE+vlG1ZuqK8vLl7qUPIWiFiV/IQyYPT7j806fpUw3O7WJd8hyR51NjYyMaNGysygbg7GzdupLGxMfXGGVCfh0iult5DfZL5q9IZnrvwvIX5jUnyauDAgaxZs4b169eXOpSsNDY2MnDgwLyWqeQhkiO/96KkLVKaiqTyNTQ0MGzYsFKHUVbUbCWSi18enXMRmopEKpGSh0gOfP3ypLWOEwb2y/6iD5EypuQhkq2bDkm62h3e7ZS8ZVjDc6VSKXmIZMk/WJe0UvGZIQM0PFeqlpKHSDam7Zd0tTtsq0v+7zVxxMR8RiRSVEoeIllwT16pmH5Aj5S1jqvGXJXfoESKSMlDJFNp1Dru2VfDb6W6KXmIZOL6wSlrHd/u2wuSznIFM46bkd+4RIqsbC8SNLM3gA+AncAOd28ubUQiwLYtSROHOzzf1DVV7qj6GVel+pVt8ggd7+4bSh2ECADT9kta63CHXU7KxKGOcqkGarYSyUCq6/1GDfpUyjLUUS7VoJyThwMPm9kiM5uUcmuRQgprHYm4wzavo67zB0mLsVTVEpEKUc7NVse4+1ozOxB4xMyWu/tTkZVhQpkEMHjw4FLFKLXggcuA1LWOI4cMp45tSbdZev7SfEUlUlJlW/Nw97Xhz3eB+4GjYtb/2t2b3b25d+/epQhRasXC25Kudoet3kBdffLEIVJNyjJ5mFk3M+seeQ6cDLxY2qikJv20b1qbjW68JGVHuYbnSjUp12arPsD9FrQTdAJ+5+7/U9qQpCbtbEu62h2W+wAa+89JWZSG50o1Kcvk4e4rgU+XOg6pcWlcSe4OX9p+I91tStJtdZtZqTZl2WwlUikO2v47uo9InjhAt5mV6qPkIRJPOrPmevjvk+K/SPfskGqk5CESK81byx6y/bd0G35dyu10zw6pRkoeIrE2LE+6OrrWkeqiQJFqpeQhEi2DWkf3/n9MuZ2G50q1UvIQiZZGreN/d4XzV+33bMriNDxXqpWSh0hEmkNzz2v/UVrFafZcqWZKHiIZOGj77wDo/snUw3M1e65UMyUPEUi71pGu7p10G1qpbmV5hXmxtCxu5caHXmXt5jb692jiirEjmHD4gFKHJcWWZid5pNbRf+R1fLAj+bbPnPNMrlGJlLWarXm0LG5l6n3LaN3chgOtm9u4ZM4SPvnjB2lZ3Frq8KSY0ugkb4+qdXywI/nw3Lra/beSGlKzNY8bH3qVtvadey1va9/FJXOWcMmcJQDUWXBr0QGqmVSnafuntdnfhLWO0UfPYsX7ybf92XE/yzUqkbJXs8lj7ebks6VG7Aq/cUZqJpGkEnHM8J7cedFn8xydFM+upGvdYXL7xbtfr3h/RcoSNTxXakHNJo/+PZpoTTOBJPP0ik0MnTJv92szOOfowUyfMDLnsqXA0uwkn7vrWABmThzNj1PcCFAXBUqtqNnkccXYEVw6ZwkZDKBJizv8dsGb/HbBm7uXNdYby6/7cp6PJMUQ6SQH+PHSU1Jur1qH1IqaTR4TDh/AwtWbOnzIF8rHO71j7QRYNUMfMiWV4dDcgw/sxtspihzTd0zucYlUiJpNHgDTJ4ykeUhPrv3TS7z3UXvRjuvQIZlA0CSizvgiuX5wWptF1zq29rkSUgzPnTV2Vi5RiVSUlMnDzLoClwOD3f0iMzsYGOHuDxQ8uiKYcPiADh/aLYtbmTb3JTa3FS+ZAHE7499Q7aQwtm1Jutod1nqP3a+PGd6TpSmG59ZTn4/IRCqGeYrLZs1sDrAIOM/dDzOzJuBZdx9dhPjS0tzc7AsXFuZObS2LW/nhvUvZtiP5qJxiUDLJgxTNVRAkj2Hb9tQ6eh16Ddt8W9J9lp2/LOfQRIrNzBa5e3M2+6bTbDXc3Sea2dkA7t5mZpbNwSpRbM0E4JxZz/L0ik1FjyW2qevcMRrVlW+xQ3P7dO/MRykSh0gtSid5bA9rGw5gZsOBmv5vir2uo1TJJHZUF6h2klSGQ3MBrj57O1P+N3mxw/cdno/oRCpKOs1WJwFXAYcCDwPHABe4+xMFjy5NhWy2ykapkkk8SiZR0kge0c1VkN7suWqykkpV0GYrd3/EzF4AxhCMMp3s7huyOVitiK2ZtCxu3aszvFhim7pqNplkMWvuGzNOZeTs5MnjwMYDc41MpCKlU/P423jL3f2pgkSUhXKreaQj9kO9lKo+oVw/OK0RVtG1jj7dO/PRwMtSFq1ah1SyQneYXxH1vBE4imD01RezOaAEYj+wS5lMqr52kkbi2BXzHeq5H53EyNnJizVqZtyIyF5S1jz22sFsEHCDu59dmJAyl23NY97KeVz/3PVs2b7nw6WpvgmAtp3BvFeG4XmfxCQNTimOGldFj62L/vt2Z9nqNUzv2YPf79u9w5SIHpMIKvo9V6CJIybuvvPivJXzuPaZa3f/D0arszp2+S76devH5CMmA3DrC7fy9ta36dutL5OPmBx3iph5K+el3C6dbbJViLLzUWYuNY9skocBS929bMaIZpM85q2cx1Xzr2KHp7hsWKpHpGPDTNmhDE0cMZHDDzycqf87Na0vbA11Dbh7h//hxvpGpn1uWocP0Xkr5zHtmWl8vPPjhNuls022ClF2vsosaPIws1+w50twHTAaeMPdz83mgIWQTfI4+Q8ns27rugJFJGUrkjyk7NRZHX269sn5/7Jft348/LWHd79O9L8evV0622SrEGXnq8xC93lEfyrvAO5y96ezOVg5eXtrqmnuRKSYdvmuvPxfxpaRqMzo5elsk6948lF2IeNNVzpDdVN0G1amvg37sq49eUeqiBRPvmoefbv13et1vDKjt0tnm1ziyXfZhYw3XQlvtmxmy8xsaZzHMjNLcUuc3JnZKWb2qpm9bmapr9TK0OT3NtNpV+nnq5IiivR5ZNjPJ8Xx9b/5OpOPmJz2KLaGugY6Wcfvv431jbs70iMmHzGZxvrGpNuls022ClF2IeNNV7Kax2lFiyKGmdUD/wKcBKwB/mJmc9395Xwd49T1a6BbE9f33J8t9XtyaJMHw5za6oI/YKN8Rj1JjsLRVj/t2YM/dBhtZTjqCiml6NFWQF5HW0VeJ9sunW2yVYiyCxlvujIebVUMZvZZYJq7jw1fTwVw9+vjbZ/VUN1bDoMtb+UYqVSayJ97dKLY7p3o/NV/ZdTi6SlH+eiiQKkmuXSYJ2y2iip8jJn9xcw+NLPtZrbTzN7P5mAZGABEf7KvCZdFxzXJzBaa2cL169dnfoQTroa6hpyClMoSL3EAdLYdcN9FuKsZUyRdKZMH8EvgbOA1oAn4NvCLQgYFcRs9O3wldPdfu3uzuzf37t078yOMOhMm/AqaenZc3tAteCQNpbZFdx2U+pFprImcMLBfyrImjpiYwVkSqW5p3YbW3V83s3p33wn8p5k9U+C41gCDol4PBNbm/SijzgwekpHYdFpO83RBx+lVhsXEtqrLN+Lu826nTik7PaLb5EVqXTrJ4yMz6wwsMbMbgHVAtxT75OovwMFmNgxoBc4C4v/XS8nFmwurnObqirbJ9+EA+zCzAt2pofufiaQlnWarb4bbfQ/YSlAj+Gohg3L3HeHxHgJeAe5x95cKeUzJrzdmnNrhUS6u3XEe27zj/cZHDhmYcr+lq96E2eMKFZZIxUlnepIzgP92L997cVbilOy1rpQ1k3F187mm0284oC6ogYwcOihxk5U75s7S1Wv2LNunH/xgeREiFSmsQs9t9Z8E068/BdwNPBTWDMqGkkd1KHZCGVc3nycOnosnmyjRnWVvJBjSPezzcP7cwgUoUmAFn1XXzBqALwETgWOBR9z929kcsBCUPKpTMZLJPodMSd5Pnix5RDRfCKfdnNe4RIqh0BMj4u7tZvYgwXDZJmA8wZBdkYIp9A2zug2/LvkG7sxYvzF1QQtvg7/eDafP1Og9qRnpNFudQjDa6XjgCWAO8HA5NV2p5lG7ckkoeal1xKrvAuN/qSQiFaHQNY8LCPo6vlPOneZSmzKtnaQ9V5kD7pnf/mPnNvzei9h233donPZeBjuKVJaynNsqU6p5SLqualnGbxe8SbcRU5LeUNAdPlw+g3F185nZ8CuMzJJI5N/qf3d9ivPaf8S5YwYzfULZ3HxTBChCn4dItYh8gLdsTp44Ih/+c3cdy9xtx/Jg5ys4hFYgvSQS2ea4updY0fkbXPr8xQxdcOzu9Z3qjJu+/mkmHD4gQQki5U01D6k5I2cnrwFEah3xXNvpdr5Z/2hWNZG13oNjtv8q4TaqnUixFWSorpn9AJjj7mU/b7mSh2QiVfIA2PrKDJLNsXttp9s5r/7RrJqyUiWRiGOG9+TOiz6b/gFEMlSoZqsBwDNmtgq4C/i9u2/I5iAi5WL07NEpt4l3z46Tbn6C197duvv1NTv+DoDz6h8FMmvK6s9mVnX5Br/ZeeLucuJ5esWmvQYAlNNUL1LbkjZbWTAb3N8SDNUdD/yVIJHc7+4fFCXCNKjmIelKp9aR7g2fzpn1LE+v2LS7FgK5dapnY+bE0eo3kawV/Arz8CD1wInADGCEu3fN5oCFoOQh6Zi+YDpzXp2TdJsDGw/ksYmPZVz2VS3LuOyFk9jf2jK+nW2kg/6g7b/L+LixGuuN5dd9OedypDYUY3qSkQS1j4nARuAud5+ZzQELQclD0pHPWkdCD1wGC28L7omewW75qIXE06d7Z5770Ul5K0+qS6E6zA8muIPgWcBOggsF73L3ldkGWihKHpKOoiSPiDCJZCqftZBEDj6wG49c9oWClS+Vo1DJYyVB/8bd7p6n/6jCUPKQVIqaOKJN2y/jXaL/JfNdE4lHyaR2FWq01VigT2ziMLPjgLXuviKbA4rUlGlbgp83HQIfrktrl+g+k+PqXmJll29wSfvFzN11bOKdcvDau1s7jOpSMpF0JKt5PABc6e5LY5Y3A9e4++lFiC8tqnlIMs2/aWZbimnZClLriLX0Hrjvoqx2LUZzViJKJtWrUDWPobGJA8DdF5rZ0GwOJlIKqRJH0Yw6M3hkUAuJiMzDtarLN9jp8IkiJpHYmklTQx3Xf2WUhgjXuGTJozHJuqZ8ByJSCPNWpp6yffi+w4sQSZTILWx/eTRsyOx2tmbQyeCNxm9Al/04p9ccnl6xqQBBJtbWvotL5izhkjlLgpiAczS1Ss1J1mx1F/Bnd58Vs/xC4GR3n1iE+NKiZitJpGQd5enKoSlrt0i/SmjUNf/D+9t25lZmjjS1SmUo1GirPsD9wHZgUbi4GegMnOHub2dzwEJQ8pBEyj55RGQ5tHePOkhw/5DYqVVKQU1d5amgFwma2fHAYeHLl9z9z9kcqJCUPCSeikkc0XJNIl32g6lvptwscl+TUtm/awPXnP4pJZMSK8r0JOVMyUPiqcjkEfHTvrCzLfv99+m3p28lDZF5ukpBfSalo+Sh5CExRs8ezU6St/uXbeKIyEd/yFdmZXU/9VL0m+h+JsWn5KHkITEqutYRK9emrAxrIfEUo2ZSb8aK679My+JWpt63lLb2PXdUUZ9JYSh5KHlIjFTJo556lpy/pDjB5EsW14d0kGUtJJWWxa1c8fsltCe7e1aaZk4czWVzliS8EVedwS6HAT2auGLsCCWTHCl5KHlIjFTJo2JqHfFcPxi2bUm9XSK9DoHvPZe/eKLEqzWkq96Mvvs10ro5vb6epoZ6vnrkAB5fvp61m9vor4SSMSUPJQ+JUdXJA0raH5KJlsWtTJv7Epvb2lNue+6Ywdy54E0y+UQy6LB9U0M9139lpBJImnJJHnX5DkakHHTv1D3huopPHBB86E/bEtQisnXfRUFTWAFNOHwAS645mTdmnMobM07l3DGD97rPibGns7x/j8wmr4hNNG3tO7nxoVdzCVnSpJqHVK3P3fk5Ptix527J3Tt155lznilhRAWUa1NW84Vw2s35iydLLYtbk/Z5pMOAVTH3em9Z3MqND72q5q0YVdVsZWbTgIuA9eGiK939v5Pto+QhEqqCJJJuv0lsk1XEgB5NPD3lizHlLaOtfc/QYzVvBaoxeXzo7jelu4+Sh0iUIl2lXiyRWkPr5jbqzdjpzoAeTRx/SG/uXdSaMikcM+PPcTvh4yWZWqudFGpKdhGpRKfdHDyy7VTftiW4A2IBR2VlYsLhAxJ+iDcP6ZnyA39tgtFb0ctjayetm9uYet+y3ceXvZVrzeMC4H1gIXC5u+8145uZTQImAQwePPjI1atXFzFKkQqSxdTvHZRBU1Yu0ql5pFs7qTYVN9rKzB41sxfjPMYD/woMB0YD64B/ileGu//a3Zvdvbl3797FC16k0nzvuWBkVvOF2e2/8Da4tld+YyqiK8aOoKmhvsOypoZ6rhg7YvfrdGon0lFJkoe7n+juh8V5/NHd33H3ne6+C5gFHFWKGEWqzmk3Zz+819uDpqwHLst/XAU24fABXP+VkQzo0YQR1CZi+0USDRHOdOhwLSnHZqt+7r4ufH4pcLS7n5VsH3WYi2Qhl471MukPyZdCj8gq1874imu2SuEGM1tmZkuB44FLSx2QSFWK1ETqs/h2vWF5UBNZek/+4yqBdGon2YokptbNbTh7OuNbFrfmXHYplV3NIxuqeYjkaPY4WPVkdvsO+zycPze/8VSRcu6Mr7aah4gU2/lzs6+FrHoyGNElcVVrZ7ySh4js8eO3sxuVFWnGuuWwqmnKypdq7YxX8hCRjnIZlbXlreDCRNVEdktnqHA2Wha3csyMPzNsyjyOmfHnovehKHmISHy5XB8SqYn8fFjN10QK0RlfDp3w6jAXkfTkcqV6XWeY8C8Fv39IrchXJ7w6zEWk8CI1kWyas3ZtD5qzKvAiw3JUDp3wSh4ikpnvPRc0ZVl96m1jLbwNrutf801ZuSqHTnglDxHJ3Gk3wzWbsquJtG9VLSRHheqEz4SSh4jk5nvPBRcKZmrhbfDT3qqFZKGQV8SnSx3mIpI/2c6XpavUS0Id5iJSHk67Gb4yCzp3y2y/VU+qGavCKHmISH6NOhOuXJv5NSKL7ihYSJJ/Sh4iUjiRmkhDGjURD6dDX3pPMM3JtB6a7qSMKXmISGGNOhN+tDZ1p7rVB4niT98PpjnB90x3MntcUUKV9Cl5iEhxnD83eTPWkRfAYz+B9jgXuq16UgmkzCh5iEjxRCZdjL7I0OqD16fdDFvWJN531ZOaK6uMaKiuiJSPWw4Lm6xSaOoJX/q55srKkYbqikh1OOHq9LZr2xT0jagWUjJKHiJSPkadmf7V6u1tQR+JlISSh4iUl/Pnpp9AtrylYb0louQhIuXn/LnB9SFNPVNsaB2H9aopq2iUPESkPI06E364KkkSMSBmwI+asopGyUNEylt0EtlvEGDhzwQjRZMN95W86VTqAERE0jLqzI5DcxMN691vYPFiqmGqeYhIZTrhamiIuXNeQ1P6w30lJ0oeIlKZRp0Jp/9zx6as0/85/oWDmmwx79RsJSKVK7YpK57IZIuRObMio7Ii+0tWVPMQkeoWb7JFjcrKmZKHiFS3RKOvNCorJ0oeIlLdEo2+0qisnJQkeZjZ183sJTPbZWbNMeummtnrZvaqmY0tRXwiUkXSHZWlTvWMlKrD/EXgK8C/Ry80s0OBs4BPAf2BR83sb9wj96cUEclQpFP8sZ8ETVX7DQwSR3RnuTrVM1aS5OHurwCYWeyq8cDd7r4NWGVmrwNHAc8WN0IRqSqpRmUl61RX8oir3Po8BgDRl4yuCZftxcwmmdlCM1u4fv36ogQnIlVKneoZK1jyMLNHzezFOI/xyXaLsyzuBDbu/mt3b3b35t69e+cnaBGpTepUz1jBmq3c/cQsdlsDDIp6PRBYm5+IREQSOOHqjn0eoKlOUii3Zqu5wFlm1sXMhgEHA8+XOCYRqXaZTHUiQIk6zM3sDOAXQG9gnpktcfex7v6Smd0DvAzsAP5BI61EpCjSmepEdivVaKv7gfsTrLsOuK64EYmISCbKrdlKREQqgJKHiIhkTMlDREQypuQhIiIZU/IQESm0Kpx0UXcSFBEppCqddFE1DxGRQqrSOxkqeYiIFFKVTrqo5CEiUkhVOumikoeISCGleyfDCqPkISJSSFU66aJGW4mIFFoVTrqomoeIiGRMyUNERDKm5CEiIhlT8hARkYwpeYiIVKISz5el0VYiIpWmDObLUs1DRKTSlMF8WUoeIiKVpgzmy1LyEBGpNGUwX5aSh4hIpSmD+bKUPEREKk0ZzJel0VYiIpWoxPNlqeYhIiIZU/IQEZGMKXmIiEjGlDxERCRjSh4iIpIxc/dSx5AzM1sPrI5a1AvYUKJwyonOQ0DnIaDzENB5CPQCurl772x2rorkEcvMFrp7c6njKDWdh4DOQ0DnIaDzEMj1PKjZSkREMqbkISIiGavW5PHrUgdQJnQeAjoPAZ2HgM5DIKfzUJV9HiIiUljVWvMQEZECUvIQEZGMVXzyMLNBZva4mb1iZi+Z2eRweU8ze8TMXgt/7l/qWAvNzOrNbLGZPRC+rrlzAGBmPczsD2a2PPy7+GytnQszuzT8f3jRzO4ys8ZaOQdmdruZvWtmL0YtS/jezWyqmb1uZq+a2djSRJ1fCc7BjeH/xFIzu9/MekSty/gcVHzyAHYAl7v7J4ExwD+Y2aHAFOAxdz8YeCx8Xe0mA69Eva7FcwBwK/A/7n4I8GmCc1Iz58LMBgDfB5rd/TCgHjiL2jkHdwCnxCyL+97Dz4qzgE+F+/zKzOqLF2rB3MHe5+AR4DB3HwX8HzAVsj8HFZ883H2du78QPv+A4INiADAemB1uNhuYUJIAi8TMBgKnAv8RtbimzgGAme0L/C1wG4C7b3f3zdTeuegENJlZJ6ArsJYaOQfu/hSwKWZxovc+Hrjb3be5+yrgdeCoYsRZSPHOgbs/7O47wpcLgMg9a7M6BxWfPKKZ2VDgcOA5oI+7r4MgwQAHljC0YpgJ/COwK2pZrZ0DgIOA9cB/hk14/2Fm3aihc+HurcBNwJvAOmCLuz9MDZ2DOBK99wHAW1HbrQmXVbu/Ax4Mn2d1DqomeZjZPsC9wCXu/n6p4ykmMzsNeNfdF5U6ljLQCTgC+Fd3PxzYSvU2z8QVtuePB4YB/YFuZnZuaaMqWxZnWVVfv2BmPyJo7r8zsijOZinPQVUkDzNrIEgcd7r7feHid8ysX7i+H/BuqeIrgmOAcWb2BnA38EUz+y21dQ4i1gBr3P258PUfCJJJLZ2LE4FV7r7e3duB+4DPUVvnIFai974GGBS13UCCJr6qZGbnA6cB5/iei/yyOgcVnzzMzAjat19x95ujVs0Fzg+fnw/8sdixFYu7T3X3ge4+lKDj68/ufi41dA4i3P1t4C0zGxEuOgF4mdo6F28CY8ysa/j/cQJBX2AtnYNYid77XOAsM+tiZsOAg4HnSxBfwZnZKcAPgXHu/lHUquzOgbtX9AM4lqCKtRRYEj6+DBxAMKritfBnz1LHWqTz8QXggfB5rZ6D0cDC8G+iBdi/1s4FcC2wHHgR+C+gS62cA+Augr6edoJv1Rcme+/Aj4AVwKvAl0odfwHPwesEfRuRz8l/y+UcaHoSERHJWMU3W4mISPEpeYiISMaUPEREJGNKHiIikjElDxERyZiSh9QcM/tROOPsUjNbYmZHh8svMbOuWZY5zcx+kGT9F8zs2Zhlncxs98VrCfZ5IJt4RApNyUNqipl9luAK2yM8mF30RPbM63MJwSSChfAUMDCcfy3iROBFD+dcEqkkSh5Sa/oBG9x9G4C7b3D3tWb2fYJ5oB43s8cBzOxsM1sW3hPj55ECzOwUM3vBzP5qZo/FHsDMLjKzB82sKbLM3XcBvwcmRm16FnCXmR1lZs+EEzk+E3V1fHSZHWo2YUxDw+fnmtnzYS3q36tkSnEpc0oeUmseBgaZ2f+Z2a/M7PMA7v7PBPP5HO/ux5tZf+DnwBcJrlj/jJlNMLPewCzgq+7+aeDr0YWb2feA04EJ7t4Wc+y7CBIGZtaFYCaEewmuBP9bDyZyvBr4Wbpvxsw+SZCQjnH30cBO4Jx09xfJVqdSByBSTO7+oZkdCRwHHA/MMbMp7n5HzKafAZ5w9/UAZnYnwX1CdgJPeXDfA9w9+p4J3ySYCmKCBxMSxh77L2a2T1iz+CSwwN3fM7NBwGwzO5hgqp2GDN7SCcCRwF+CaaxoorYmPJQSUfKQmuPuO4EngCfMbBnBRHl3xGwWb5rqyPJEc/q8SFBLGQisSrDN3QS1j08S1EQAfgo87u5nhE1RT8TZbwcdWwoao+KZ7e5TExxPpCDUbCU1xcxGhN/wI0YDq8PnHwDdw+fPAZ83s15hH8LZwJPAs+HyYWF5PaPKWgx8B5gbNnvFcxdwLkFz2Nxw2X5Aa/j8ggT7vUEwtTxmdgTBvTogmOTva2Z2YCQeMxuSoAyRvFHykFqzD0ET0ctmthQ4FJgWrvs18KCZPR6OgJoKPA78FXjB3f8YNmNNAu4zs78Cc6ILd/f5wA+AeWbWK/bg7v4y8BHBtPlbw8U3ANeb2dME9xuP516gp5ktAb5LcA/qSHlXAQ+H7+cRgkEBIgWlWXVFRCRjqnmIiEjGlDxERCRjSh4iIpIxJQ8REcmYkoeIiGRMyUNERDKm5CEiIhn7/5WN+2KPwPMkAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(Xa, Ylambdaa)\n",
    "plt.scatter(Xa, Ylambdab)\n",
    "plt.scatter(Xa, Ylambdac)\n",
    "plt.title('Scatter of Control Variates')\n",
    "plt.legend(['Xa', 'Xb', 'Xc'])\n",
    "plt.xlabel(\"Stock Value\")\n",
    "plt.ylabel(\"CV value\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how for Xb is more negative in the region of $S > 75$ which is the barrier limit, yet as we saw above its variance reduction factor was less. Xc meanwhile remains mostly flat, as it mimics the behavior of reporting zero when the barrier option payoff is zero."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}