MC-Heston-European.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimal Hedged Monte Carlo\n",
    "\n",
    "Author: Jerry Xia\n",
    "\n",
    "Date: 2018/06/19"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction\n",
    "\n",
    "This is a Python Notebook about variance reduction Monte Carlo simulations. In this script, I implemented the following variance reduction methods as well as their antithetic variates' version:\n",
    "\n",
    "* regular Monte Carlo\n",
    "* Monte Carlo with delta-based control variates\n",
    "* optimal hedged Monte Carlo\n",
    "\n",
    "Due to the significance and robustness, I mainly focus on the optimal hedged Monte Carlo (OHMC) in option pricing. We invoke this method to price European options and make comparison with other methods.\n",
    "\n",
    "### 1.1 Facts\n",
    "* The option price is not simply the average value of the discounted future pay-off over the objective (or historical) probability distribution\n",
    "* The requirement of absence of arbitrage opportunities is equivalent to the existence of \"risk-neutral measure\", such that the price is indeed its average discounted future pay-off.\n",
    "* Risk in option trading cannot be eliminated\n",
    "\n",
    "### 1.2 Objective\n",
    "* It would be satisfactory to have an option theory where the objective stochastic process of the underlying is used to calculate the option price, the hedge strategy and the *residual risk*.\n",
    "\n",
    "### 1.3 Advantages\n",
    "* It is a versatile methods to price complicated path-dependent options.\n",
    "* Considerable variance reduction scheme for Monte Carlo\n",
    "* It provide not only a numerical estimate of the option price, but also of the optimal hedge strategy and of the residual risk.\n",
    "* This method does not rely on the notion of risk-neutral measure, and can be used to any model of the true dynamics of the underlying"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Underlying dynamics\n",
    "\n",
    "### Black-Scholes Model\n",
    "$$dS = \\mu S dt + \\sigma S dW_t$$\n",
    "$$log S_{t+1} = log S_t +(\\mu - \\frac{\\sigma^2}{2})\\Delta t + \\sigma \\sqrt{\\Delta t} \\epsilon$$\n",
    "where\n",
    "$$\\epsilon \\sim N(0,1)$$\n",
    "In risk neutral measure, $\\mu = r - q$. \n",
    "### Heston Model\n",
    "The basic Heston model assumes that $S_t$, the price of the asset, is determined by a stochastic process:\n",
    "$$\n",
    "dS_t = \\mu S_t dt + \\sqrt{v_t} S_t d W_t^S\\\\\n",
    "dv_t = \\kappa (\\theta - v_t) dt + \\xi \\sqrt{v_t} d W_t^v\n",
    "$$\n",
    "where \n",
    "$$E[dW_t^S,dW_t^v]=\\rho dt$$\n",
    "In risk neutral measure, $\\mu = r - q$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Methodology\n",
    "\n",
    "### 3.1 Notation\n",
    "Option price always requires to work backward. That is because the option price is known exactly at the maturity. As with other schemes, we determine the option price step by step from the maturity $t=K\\tau=T$ to the present time $t=0$. The unit of time being $\\tau$, for example, one day. We simulate $N$ trajectories. In trajectory i, the price of the underlying asset at time $k\\tau$ is denoted as $S_k^{(i)}$. The price of the derivative at time $k\\tau$ is denoted as $C_k$, and the hedge function is $H_k$. We define an optimal hedged portfolio as\n",
    "$$W_k^{(i)} = C_k(S_k^{(i)}) + H_k(S_k^{(i)})S_k^{(i)}$$\n",
    "The one-step change of our portfolio is\n",
    "$$\\Delta W_k^{(i)}= df(k,k+1) C_{k+1}(S_{k+1}^{(i)}) - C_k(S_k^{(i)}) + H_k(S_{k}^{(i)}) (df2(k,k+1) S_{k+1}^{(i)} - S_{k}^{(i)})$$\n",
    "Where $df(k,k+1)$ is the discounted factor from time $k\\tau$ to $(k+1) \\tau$, $df2(k,k+1)$ is the discounted factor considering dividend $e^{-(r-q)(t_{k+1}-t_k)}$\n",
    "\n",
    "### 3.2 Objective\n",
    "The optimal hedged algorithm can be interpreted as the following optimal problem\n",
    "\n",
    "\\begin{align}\n",
    "\\mbox{minimize}\\quad & \\quad Var[\\Delta W_k]\\\\\n",
    "\\mbox{subject to}\\quad & \\quad E[\\Delta W_k]=0\n",
    "\\end{align}\n",
    "\n",
    "It means we should try to minimize the realized volatility of hedged portfolio while maintaining the expected value of portfolio unchanged.\n",
    "\n",
    "### 3.3 Basis Functions\n",
    "The original optimization is very difficult to solve. Thus we assume a set of basis function and solved it in such subspace. We use $N_C$and $N_H$ to denote the number of basis functions for price and hedge.\n",
    "\n",
    "\\begin{align}\n",
    "C_k(\\cdot) &= \\sum_{i=0}^{N_C} a_{k,i} A_i(\\cdot)\\\\\n",
    "H_k(\\cdot) &= \\sum_{i=0}^{N_H} b_{k,i} B_i(\\cdot)\n",
    "\\end{align}\n",
    "\n",
    "The basis functions $A_i$ and $B_i$ are priori determined and need not to be identical. The coefficients $a_i$ and $b_i$ can be calibrated by solving the optimal problem.\n",
    "\n",
    "### 3.4 Numerical Solution\n",
    "\n",
    "\\begin{align}\n",
    "\\mbox{minimize}\\quad & \\quad \\frac{1}{N} \\sum_{i=1}^N \\Delta W_k^{(i)2}\\\\\n",
    "\\mbox{subject to}\\quad & \\quad \\frac{1}{N} \\sum_{i=1}^N \\Delta W_k^{(i)}=0\n",
    "\\end{align}\n",
    "\n",
    "Denote the discounted forward underlying price change at time $k\\tau$ as\n",
    "\n",
    "$$\\Delta S_k = df2(k,k+1) S_{k+1} - S_k$$\n",
    "\n",
    "Define\n",
    "\n",
    "\\begin{align}\n",
    "Q_k &= \\begin{bmatrix}\n",
    "    -A_{k,1}(S_k^{(1)}) & \\cdots & -A_{k,N_C}(S_k^{(1)}) & B_{k,1}(S_k^{(1)})\\Delta S_k^{(1)}& \\cdots  & B_{k,N_H}(S_k^{(1)})\\Delta S_k^{(1)} \\\\\n",
    "    -A_{k,1}(S_k^{(2)}) & \\cdots & -A_{k,N_C}(S_k^{(2)}) & B_{k,1}(S_k^{(2)})\\Delta S_k^{(2)}& \\cdots  & B_{k,N_H}(S_k^{(1)})\\Delta S_k^{(2)} \\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "    -A_{k,1}(S_k^{(N)}) & \\cdots & -A_{k,N_C}(S_k^{(N)}) & B_{k,1}(S_k^{(N)})\\Delta S_k^{(N)}& \\cdots  & B_{k,N_H}(S_k^{(N)})\\Delta S_k^{(N)}\n",
    "    \\end{bmatrix}\\\\\\\\\n",
    "c_k &= (a_{k,1}, \\cdots a_{k,N_C}, b_{k,1}, \\cdots, b_{k,N_H})^T\\\\\\\\\n",
    "v_{k} &= df(k,k+1) C_{k+1}(S_{k+1}^{})\n",
    "\\end{align}\n",
    "\n",
    "As for $v_k$, note that we know the exact value at maturity, which means there is no need to approximate price in terms of basis functions, that is\n",
    "\n",
    "\\begin{align}\n",
    "v_k = \\begin{cases}\n",
    "df(N-1,N)\\ payoff(S_N),\\quad & k=N-1\\\\\n",
    "df(k,k+1)\\ \\sum_{i=1}^{N_C} a_{k+1,i} A_i(S_{k+1}), \\quad & 0<k<N-1\\\\\n",
    "df(0,1)\\ C_1(S_1), \\quad & k=0\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "Then, the optimization problem can be expressed as\n",
    "\n",
    "\\begin{align}\n",
    "\\arg\\min_{c_k}\\quad & \\quad (v_{k} + Q_k c_k)^T (v_{k} + Q_k c_k)\\\\\n",
    "\\mbox{subject to}\\quad & \\quad 1_{[N\\times1]}^T (v_{k}  + Q_k c_k)=0\n",
    "\\end{align}\n",
    "\n",
    "In step k, since we already know the information ($v_{k}$) in step k+1. By canceling the constant term, the optimal problem can be simplified as the following \n",
    "\n",
    "\\begin{align}\n",
    "\\arg\\min_{c_k}\\quad & \\quad 2 v_{k}^T Q_k c_k + c_k^T Q_k^T Q_k c_k\\\\\n",
    "\\mbox{subject to}\\quad & \\quad 1_{[N\\times1]}^T v_{k}  + 1_{[N\\times1]}^T Q_k c_k=0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Convex Optimization Problem\n",
    "\n",
    "Let us first review the standard form of linear constrained quadratic programming problem:\n",
    "\n",
    "\\begin{align}\n",
    "    \\min_{x} \\quad & \\frac{1}{2} x^T P x + q^T x\\\\\n",
    "    \\mbox{subject to} \\quad &G x \\preceq h\\\\\n",
    "    &A x = b\n",
    "\\end{align}\n",
    "\n",
    "Note that $x^T$ means the transpose of vector x, and $G x \\preceq h$denotes the inequality is taken element-wise over the vectors $G x$ and $h$. The objective function is convex if and only if the matrix $P$ is positive-semidefinite(Hermitian matrix all of whose eigenvalues are nonnegative), which is the realm we concern with.\n",
    "\n",
    "Recall that the constrained optimization problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\arg\\min_{c_k}\\quad & \\quad  v_{k}^T Q_k c_k + \\frac{1}{2}c_k^T Q_k^T Q_k c_k\\\\\n",
    "\\mbox{subject to}\\quad & \\quad 1_{[N\\times1]}^T v_{k}  + 1_{[N\\times1]}^T Q_k c_k=0\n",
    "\\end{align}\n",
    "\n",
    "Correspondingly, we make the connection by letting\n",
    "\n",
    "\\begin{align}\n",
    "    x &= c_k\\\\\n",
    "    P &= Q_k^T Q_k\\\\\n",
    "    q &= Q_k^T v_k\\\\\n",
    "    A &= 1_{[N\\times1]}^T Q_k\\\\\n",
    "    b &= -1_{[N\\times1]}^T v_{k}\n",
    "\\end{align}\n",
    "\n",
    "The hard work is almost over right now. As you would always find, formulating the problem is usually the hard step. Invoking a solver is straightforward.\n",
    "\n",
    "Note that when $k=0$, the degree of freedom of the quadratic problem decreases to 2. Because here the only concerns are price and hedge at time zero (we don't need to project them into a high dimension space). Let $x=[C_0, H_0]^T$\n",
    "\n",
    "\\begin{align}\n",
    "    Q_0 &= \\begin{bmatrix}\n",
    "    -1 & \\Delta S_0^{(1)}\\\\\n",
    "    \\vdots & \\vdots\\\\\n",
    "    -1 & \\Delta S_0^{(N)}\n",
    "    \\end{bmatrix}\\\\\n",
    "    P &= Q_0^T Q_0\\\\\n",
    "    q &= Q_0^T v_0\\\\\n",
    "    A &= 1_{[N \\times 1]}^T Q_0\\\\\n",
    "    b &= -1_{[N \\times 1]}^T v_0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Variance reduction and other methods\n",
    "The rate of convergence of the Monte Carlo simulation is $O\\left(\\max \\left( \\Delta t, \\frac{1}{N_x} \\right)\\right)$. The variance reduction techniques are used to reduce the constant factor corresponding to the Monte Carlo approximation $O \\left(\\frac{1}{N_x}\\right)$. Some of the most used variance reduction techniques are:\n",
    "\n",
    "* Control Variates\n",
    "* Antithetic Variates\n",
    "* Moment Matching\n",
    "\n",
    "In this part we selected antithetic variates and delta-based control variates methods as a supplement to optimal hedged monte carlo simulation.\n",
    "\n",
    "### 4.1 Antithetic variates\n",
    "The main idea of this technique is to look at the asset equation that you aretrying to simulate:\n",
    "$$d S_t^{(1)} = r S_t^{(1)} dt + \\sigma S_t^{(1)} d W_t$$\n",
    "and recognize that sinceztis a standard Brownian motion so will be−ztandthey will have the same exact distribution.  This means that the equation:\n",
    "$$d S_t^{(2)} = r S_t^{(2)} dt - \\sigma S_t^{(2)} d W_t$$\n",
    "will also generate paths of the same asset.\n",
    "The variance depends on the sign of the covariance of $payoff(S_t^{(1)})$ and $payoff(S_t^{(2)})$. It can increase the eventual variance or decrease it, both case do arise. One sufficient condition to insure variance reduction is the monotony of the payoff function. Then, when using both in the calculation of the final Monte Carlo value the variance of the estimate will be reduced.\n",
    "\n",
    "### 4.2 Delta-based control variates\n",
    "Delta hedging can be summarized succinctly in the following way:  Suppose that at time $t= 0$, we receive $C_0$ the price of an option that pays $C_T$ at time T.  The price of this option at any time $t$ is a function $C(t,S)$. Then, if we hold at any moment in time $\\frac{\\partial C}{\\partial S}(t,S) = \\frac{\\partial C_t}{\\partial S}$ units of stock, then we will be able to replicate the payout of this option $C_T$ at time T. This is in theory since of course we cannot trade continuously. So in practice we perform a partial hedge where we only rebalance at some discrete moments in time say $t_1,t_2,\\cdots,t_N$. The replicating strategy can be expressed as follow:\n",
    "$$W(t_i,S_i) = C(t_0,S_0) e^{r(t_i - t_0)} + \\sum_{j=0}^{i} \\Delta(t_j,S_j) ( S_{j+1} e^{-r(t_{j+1} - t_j )} - S_{j})e^{r(t_i - t_j)} = C(t_i,S_i)$$\n",
    "which is similar to the strategy in the optimal hedged Monte Carlo simulation where the only difference is that in OHMC, we use option and delta hedging to replicate the cash flow and here we do the opposite operation. But when implementing the delta-based control variates, we should move the hedging term to the right hand side which make it identical to the OHMC strategy. Note that here we are assumed to know the delta hedging function. It explains a lot why OHMC can reduce the variance.\n",
    "\n",
    "### 4.3 Optimal hedged Monte Carlo simulation\n",
    "**In conclusion, OHMC is just a control variates method with an optimization on top and it is more practical because we do not have an analytical formula for the hedge sensitivity (i.e. delta, gamma, etc.)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "from cvxopt import matrix, solvers\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "risk_free_rate = 0.0319\n",
    "dividend = 0\n",
    "time_to_maturity = 1\n",
    "volatility = 0.2\n",
    "strike = 100\n",
    "stock_price = 100\n",
    "V0 = 0.010201\n",
    "kappa = 6.21\n",
    "theta = 0.019\n",
    "xi = 0.61\n",
    "rho = -0.7\n",
    "\n",
    "n_trials = 1000\n",
    "n_steps = 200\n",
    "func_list = [lambda x: x**0, lambda x: x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from MonteCarlo import MonteCarlo   \n",
    "      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mc = MonteCarlo(S0=stock_price,K=strike,T=time_to_maturity,r=risk_free_rate,q=dividend,sigma=volatility,\n",
    "                kappa=kappa,theta=theta,xi=xi,rho=rho,V0=V0,underlying_process=\"Heston model\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "price_matrix = mc.simulate(n_trials=n_trials,n_steps=n_steps,boundaryScheme=\"Higham and Mao\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.4286345152762516"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mc.MCPricer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.0048861512919327"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mc.BSDeltaHedgedPricer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.233493440650192"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mc.OHMCPricer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'OHMC': 7.233493440650192, 'regular MC': 7.4286345152762516, 'Black-Scholes': 9.5095295701347027}\n"
     ]
    }
   ],
   "source": [
    "prices = mc.pricing(func_list=func_list)\n",
    "print(prices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'OHMC optimal hedge': 0.5077634149472071, 'Black-Scholes delta hedge': 0.60237525898994781, 'OHMC residual risk': array([ -2.91038305e-14])}\n"
     ]
    }
   ],
   "source": [
    "hedges = mc.hedging()\n",
    "print(hedges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Methods Comparison\n",
    "\n",
    "In this part, we discuss the performance of OHMC and regular MC. The underlying stochastic process is Heston model. In the following, we compare the absolute error with respect to the Black-Scholes price corresponding to trajactories number and runtime. The result is displayed by scatter plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def performance_comparison(n_trials,n_steps,underlying_process=\"geometric brownian motion\"):\n",
    "    if(underlying_process==\"geometric brownian motion\"):\n",
    "        mc = MonteCarlo(stock_price,strike,time_to_maturity,risk_free_rate,dividend,volatility)\n",
    "        theoretical_value = mc.BlackScholesPricer()\n",
    "    elif(underlying_process==\"Heston model\"):\n",
    "        mc = MonteCarlo(S0=stock_price,K=strike,T=time_to_maturity,r=risk_free_rate,q=dividend,sigma=volatility,\n",
    "                kappa=kappa,theta=theta,xi=xi,rho=rho,V0=V0,underlying_process=\"Heston model\")\n",
    "        theoretical_value = 6.8061\n",
    "    \n",
    "    # black scholes\n",
    "    bs_start = timeit.default_timer()\n",
    "    bs_price = mc.BlackScholesPricer('c')\n",
    "    bs_sde_estimate = 'NA'\n",
    "    bs_rmse = 'NA'\n",
    "    bs_end = timeit.default_timer()\n",
    "    \n",
    "    # regular MC\n",
    "    rmc_start = timeit.default_timer()\n",
    "    mc.simulate(n_trials,n_steps)\n",
    "    rmc_price = mc.MCPricer('c')\n",
    "    rmc_sde_estimate = mc.standard_error()\n",
    "    rmc_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    rmc_end = timeit.default_timer()\n",
    "    \n",
    "    # regular MC with antithetic variates\n",
    "    rmc_anti_start = timeit.default_timer()\n",
    "    mc.simulate(int(n_trials/2),n_steps,antitheticVariates=True)\n",
    "    rmc_anti_price = mc.MCPricer('c')\n",
    "    rmc_anti_sde_estimate = mc.standard_error()\n",
    "    rmc_anti_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    rmc_anti_end = timeit.default_timer()\n",
    "    \n",
    "    # Delta-based MC\n",
    "    dbmc_start = timeit.default_timer()\n",
    "    mc.simulate(n_trials,n_steps)\n",
    "    dbmc_price = mc.BSDeltaHedgedPricer('c')\n",
    "    dbmc_sde_estimate = mc.standard_error()\n",
    "    dbmc_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    dbmc_end = timeit.default_timer()\n",
    "    \n",
    "    # Delta-based MC with antithetic variates\n",
    "    dbmc_anti_start = timeit.default_timer()\n",
    "    mc.simulate(int(n_trials/2),n_steps)\n",
    "    dbmc_anti_price = mc.BSDeltaHedgedPricer('c')\n",
    "    dbmc_anti_sde_estimate = mc.standard_error()\n",
    "    dbmc_anti_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    dbmc_anti_end = timeit.default_timer()\n",
    "    \n",
    "    # OHMC\n",
    "    ohmc_start = timeit.default_timer()\n",
    "    mc.simulate(n_trials,n_steps)\n",
    "    ohmc_price = mc.OHMCPricer('c')\n",
    "    ohmc_sde_estimate = \"NA\"\n",
    "    ohmc_rmse = 'NA'\n",
    "    ohmc_end = timeit.default_timer()\n",
    "    \n",
    "    # OHMC with antithetic varietes\n",
    "    ohmc_anti_start = timeit.default_timer()\n",
    "    mc.simulate(int(n_trials/2),n_steps,antitheticVariates=True)\n",
    "    ohmc_anti_price = mc.OHMCPricer('c')\n",
    "    ohmc_anti_sde_estimate = \"NA\"\n",
    "    ohmc_anti_rmse = 'NA'\n",
    "    ohmc_anti_end = timeit.default_timer()\n",
    "    \n",
    "    bs_runtime = bs_end - bs_start\n",
    "    rmc_runtime = rmc_end - rmc_start\n",
    "    dbmc_runtime = dbmc_end - dbmc_start\n",
    "    ohmc_runtime = ohmc_end-ohmc_start\n",
    "    rmc_anti_runtime = rmc_anti_end - rmc_anti_start\n",
    "    dbmc_anti_runtime = dbmc_anti_end - dbmc_anti_start\n",
    "    ohmc_anti_runtime = ohmc_anti_end-ohmc_anti_start\n",
    "    \n",
    "    bs_err = np.abs(bs_price - theoretical_value)\n",
    "    rmc_err = np.abs(rmc_price - theoretical_value)\n",
    "    dbmc_err = np.abs(dbmc_price - theoretical_value)\n",
    "    ohmc_err = np.abs(ohmc_price - theoretical_value)\n",
    "    rmc_anti_err = np.abs(rmc_anti_price - theoretical_value)\n",
    "    dbmc_anti_err = np.abs(dbmc_anti_price - theoretical_value)\n",
    "    ohmc_anti_err = np.abs(ohmc_anti_price - theoretical_value)\n",
    "    \n",
    "    \n",
    "    result = {\"method\":[\"Black Scholes\",\"MC\",\"antithetic MC\",\"DBMC\",\"antithetic DBMC\",\"OHMC\",\"antithetic OHMC\"],\n",
    "              \"runtime\":[bs_runtime,rmc_runtime,rmc_anti_runtime,dbmc_runtime,dbmc_anti_runtime,\n",
    "                         ohmc_runtime,ohmc_anti_runtime],\n",
    "              \"err\":[bs_err,rmc_err,rmc_anti_err,dbmc_err,dbmc_anti_err,ohmc_err,ohmc_anti_err], \n",
    "              \"standard err estimate\":[bs_sde_estimate,rmc_sde_estimate,rmc_anti_sde_estimate,dbmc_sde_estimate,\n",
    "                              dbmc_anti_sde_estimate,ohmc_sde_estimate,ohmc_anti_sde_estimate],\n",
    "              \"root mean square error\":[bs_rmse, rmc_rmse, rmc_anti_rmse, dbmc_rmse, dbmc_anti_rmse, ohmc_rmse, ohmc_anti_rmse],\n",
    "              \"n_trials\":[n_trials]*7,\n",
    "              \"n_steps\":[n_steps]*7}\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_trials: 100; n_steps: 50\n",
      "n_trials: 100; n_steps: 100\n",
      "n_trials: 100; n_steps: 150\n",
      "n_trials: 600; n_steps: 50\n",
      "n_trials: 600; n_steps: 100\n",
      "n_trials: 600; n_steps: 150\n"
     ]
    }
   ],
   "source": [
    "n_trials_list = np.arange(100,1000,500)\n",
    "n_steps_list = np.arange(50,200,50)\n",
    "performance_df = pd.DataFrame()\n",
    "for n_trials in n_trials_list:\n",
    "    for n_steps in n_steps_list:\n",
    "        print(\"n_trials: {}; n_steps: {}\".format(n_trials,n_steps))\n",
    "        new_df = pd.DataFrame(performance_comparison(n_trials,n_steps,underlying_process=\"Heston model\"))\n",
    "        performance_df = pd.concat((performance_df,new_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>err</th>\n",
       "      <th>root mean square error</th>\n",
       "      <th>runtime</th>\n",
       "      <th>standard err estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>method</th>\n",
       "      <th>n_steps</th>\n",
       "      <th>n_trials</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">Black Scholes</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000480</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000324</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000346</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000333</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000431</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>2.703430</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.000645</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">DBMC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.548824</td>\n",
       "      <td>1.60789</td>\n",
       "      <td>0.607866</td>\n",
       "      <td>0.151893</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.422655</td>\n",
       "      <td>1.53969</td>\n",
       "      <td>2.692323</td>\n",
       "      <td>0.0604935</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.551444</td>\n",
       "      <td>1.38492</td>\n",
       "      <td>1.319747</td>\n",
       "      <td>0.12768</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.384824</td>\n",
       "      <td>1.4216</td>\n",
       "      <td>5.705039</td>\n",
       "      <td>0.0559165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.024386</td>\n",
       "      <td>1.24228</td>\n",
       "      <td>1.670167</td>\n",
       "      <td>0.12483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.204089</td>\n",
       "      <td>1.44522</td>\n",
       "      <td>12.626791</td>\n",
       "      <td>0.0584583</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">MC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>1.089855</td>\n",
       "      <td>8.18064</td>\n",
       "      <td>0.023478</td>\n",
       "      <td>0.814856</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.438077</td>\n",
       "      <td>8.80041</td>\n",
       "      <td>0.007731</td>\n",
       "      <td>0.359129</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.109861</td>\n",
       "      <td>7.34121</td>\n",
       "      <td>0.004958</td>\n",
       "      <td>0.737736</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.665718</td>\n",
       "      <td>8.43635</td>\n",
       "      <td>0.014241</td>\n",
       "      <td>0.343625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.048388</td>\n",
       "      <td>8.21575</td>\n",
       "      <td>0.007147</td>\n",
       "      <td>0.8257</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.422223</td>\n",
       "      <td>8.1912</td>\n",
       "      <td>0.047758</td>\n",
       "      <td>0.334238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">OHMC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.902040</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.015391</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.639134</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.021909</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.032222</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.040364</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.475584</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.078652</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.492817</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.044626</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.032711</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.086919</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">antithetic DBMC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.402666</td>\n",
       "      <td>1.78885</td>\n",
       "      <td>0.273746</td>\n",
       "      <td>0.248992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.410585</td>\n",
       "      <td>1.51329</td>\n",
       "      <td>1.336698</td>\n",
       "      <td>0.0842328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.635056</td>\n",
       "      <td>1.30191</td>\n",
       "      <td>0.649386</td>\n",
       "      <td>0.16236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.249550</td>\n",
       "      <td>1.40571</td>\n",
       "      <td>4.248047</td>\n",
       "      <td>0.0800031</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.012359</td>\n",
       "      <td>1.18782</td>\n",
       "      <td>0.727217</td>\n",
       "      <td>0.169679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.187533</td>\n",
       "      <td>1.37166</td>\n",
       "      <td>5.073001</td>\n",
       "      <td>0.0785805</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">antithetic MC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.562670</td>\n",
       "      <td>11.4092</td>\n",
       "      <td>0.002589</td>\n",
       "      <td>1.14527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.308745</td>\n",
       "      <td>11.6998</td>\n",
       "      <td>0.005625</td>\n",
       "      <td>0.477876</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.391531</td>\n",
       "      <td>10.0882</td>\n",
       "      <td>0.004393</td>\n",
       "      <td>1.01314</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.032050</td>\n",
       "      <td>9.8583</td>\n",
       "      <td>0.010187</td>\n",
       "      <td>0.402797</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>1.626620</td>\n",
       "      <td>18.4021</td>\n",
       "      <td>0.006634</td>\n",
       "      <td>1.84224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.664989</td>\n",
       "      <td>8.67813</td>\n",
       "      <td>0.031840</td>\n",
       "      <td>0.353536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">antithetic OHMC</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>1.194824</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.015838</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.389528</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.022475</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.736264</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.041542</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.339304</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.102591</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.518802</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.044841</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.300874</td>\n",
       "      <td>NA</td>\n",
       "      <td>0.078421</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       err root mean square error    runtime  \\\n",
       "method          n_steps n_trials                                               \n",
       "Black Scholes   50      100       2.703430                     NA   0.000480   \n",
       "                        600       2.703430                     NA   0.000324   \n",
       "                100     100       2.703430                     NA   0.000346   \n",
       "                        600       2.703430                     NA   0.000333   \n",
       "                150     100       2.703430                     NA   0.000431   \n",
       "                        600       2.703430                     NA   0.000645   \n",
       "DBMC            50      100       0.548824                1.60789   0.607866   \n",
       "                        600       0.422655                1.53969   2.692323   \n",
       "                100     100       0.551444                1.38492   1.319747   \n",
       "                        600       0.384824                 1.4216   5.705039   \n",
       "                150     100       0.024386                1.24228   1.670167   \n",
       "                        600       0.204089                1.44522  12.626791   \n",
       "MC              50      100       1.089855                8.18064   0.023478   \n",
       "                        600       0.438077                8.80041   0.007731   \n",
       "                100     100       0.109861                7.34121   0.004958   \n",
       "                        600       0.665718                8.43635   0.014241   \n",
       "                150     100       0.048388                8.21575   0.007147   \n",
       "                        600       0.422223                 8.1912   0.047758   \n",
       "OHMC            50      100       0.902040                     NA   0.015391   \n",
       "                        600       0.639134                     NA   0.021909   \n",
       "                100     100       0.032222                     NA   0.040364   \n",
       "                        600       0.475584                     NA   0.078652   \n",
       "                150     100       0.492817                     NA   0.044626   \n",
       "                        600       0.032711                     NA   0.086919   \n",
       "antithetic DBMC 50      100       0.402666                1.78885   0.273746   \n",
       "                        600       0.410585                1.51329   1.336698   \n",
       "                100     100       0.635056                1.30191   0.649386   \n",
       "                        600       0.249550                1.40571   4.248047   \n",
       "                150     100       0.012359                1.18782   0.727217   \n",
       "                        600       0.187533                1.37166   5.073001   \n",
       "antithetic MC   50      100       0.562670                11.4092   0.002589   \n",
       "                        600       0.308745                11.6998   0.005625   \n",
       "                100     100       0.391531                10.0882   0.004393   \n",
       "                        600       0.032050                 9.8583   0.010187   \n",
       "                150     100       1.626620                18.4021   0.006634   \n",
       "                        600       0.664989                8.67813   0.031840   \n",
       "antithetic OHMC 50      100       1.194824                     NA   0.015838   \n",
       "                        600       0.389528                     NA   0.022475   \n",
       "                100     100       0.736264                     NA   0.041542   \n",
       "                        600       0.339304                     NA   0.102591   \n",
       "                150     100       0.518802                     NA   0.044841   \n",
       "                        600       0.300874                     NA   0.078421   \n",
       "\n",
       "                                 standard err estimate  \n",
       "method          n_steps n_trials                        \n",
       "Black Scholes   50      100                         NA  \n",
       "                        600                         NA  \n",
       "                100     100                         NA  \n",
       "                        600                         NA  \n",
       "                150     100                         NA  \n",
       "                        600                         NA  \n",
       "DBMC            50      100                   0.151893  \n",
       "                        600                  0.0604935  \n",
       "                100     100                    0.12768  \n",
       "                        600                  0.0559165  \n",
       "                150     100                    0.12483  \n",
       "                        600                  0.0584583  \n",
       "MC              50      100                   0.814856  \n",
       "                        600                   0.359129  \n",
       "                100     100                   0.737736  \n",
       "                        600                   0.343625  \n",
       "                150     100                     0.8257  \n",
       "                        600                   0.334238  \n",
       "OHMC            50      100                         NA  \n",
       "                        600                         NA  \n",
       "                100     100                         NA  \n",
       "                        600                         NA  \n",
       "                150     100                         NA  \n",
       "                        600                         NA  \n",
       "antithetic DBMC 50      100                   0.248992  \n",
       "                        600                  0.0842328  \n",
       "                100     100                    0.16236  \n",
       "                        600                  0.0800031  \n",
       "                150     100                   0.169679  \n",
       "                        600                  0.0785805  \n",
       "antithetic MC   50      100                    1.14527  \n",
       "                        600                   0.477876  \n",
       "                100     100                    1.01314  \n",
       "                        600                   0.402797  \n",
       "                150     100                    1.84224  \n",
       "                        600                   0.353536  \n",
       "antithetic OHMC 50      100                         NA  \n",
       "                        600                         NA  \n",
       "                100     100                         NA  \n",
       "                        600                         NA  \n",
       "                150     100                         NA  \n",
       "                        600                         NA  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted_efficiency_df = performance_df.set_index(['method','n_steps','n_trials']).sort_index()\n",
    "sorted_efficiency_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Variance Reduction Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In thie section, we test more case and visualize the empirical error. Note that we do not concern Delta-based MC simulation since we cannot convert it to a matrix form, as a result, the method involes for-loop which is not efficient (see the form above). Moreover, the hypothesis which we know the delta hedge function is not realistic in practice. The result shows as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def performance_comparison2(n_trials,n_steps,underlying_process=\"geometric brownian motion\"):\n",
    "    if(underlying_process==\"geometric brownian motion\"):\n",
    "        mc = MonteCarlo(stock_price,strike,time_to_maturity,risk_free_rate,dividend,volatility)\n",
    "        theoretical_value = mc.BlackScholesPricer()\n",
    "    elif(underlying_process==\"Heston model\"):\n",
    "        mc = MonteCarlo(S0=stock_price,K=strike,T=time_to_maturity,r=risk_free_rate,q=dividend,sigma=volatility,\n",
    "                kappa=kappa,theta=theta,xi=xi,rho=rho,V0=V0,underlying_process=\"Heston model\")\n",
    "        theoretical_value = 6.8061\n",
    "    \n",
    "    # black scholes\n",
    "    bs_start = timeit.default_timer()\n",
    "    bs_price = mc.BlackScholesPricer('c')\n",
    "    bs_sde_estimate = 'NA'\n",
    "    bs_rmse = 'NA'\n",
    "    bs_end = timeit.default_timer()\n",
    "    \n",
    "    # regular MC\n",
    "    rmc_start = timeit.default_timer()\n",
    "    mc.simulate(n_trials,n_steps)\n",
    "    rmc_price = mc.MCPricer('c')\n",
    "    rmc_sde_estimate = mc.standard_error()\n",
    "    rmc_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    rmc_end = timeit.default_timer()\n",
    "    \n",
    "    # regular MC with antithetic variates\n",
    "    rmc_anti_start = timeit.default_timer()\n",
    "    mc.simulate(int(n_trials/2),n_steps,antitheticVariates=True)\n",
    "    rmc_anti_price = mc.MCPricer('c')\n",
    "    rmc_anti_sde_estimate = mc.standard_error()\n",
    "    rmc_anti_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "    rmc_anti_end = timeit.default_timer()\n",
    "    \n",
    "    # OHMC\n",
    "    ohmc_start = timeit.default_timer()\n",
    "    mc.simulate(n_trials,n_steps)\n",
    "    ohmc_price = mc.OHMCPricer('c')\n",
    "    ohmc_sde_estimate = \"NA\"\n",
    "    ohmc_rmse = 'NA'\n",
    "    ohmc_end = timeit.default_timer()\n",
    "    \n",
    "    # OHMC with antithetic varietes\n",
    "    ohmc_anti_start = timeit.default_timer()\n",
    "    mc.simulate(int(n_trials/2),n_steps,antitheticVariates=True)\n",
    "    ohmc_anti_price = mc.OHMCPricer('c')\n",
    "    ohmc_anti_sde_estimate = \"NA\"\n",
    "    ohmc_anti_rmse = 'NA'\n",
    "    ohmc_anti_end = timeit.default_timer()\n",
    "    \n",
    "    bs_runtime = bs_end - bs_start\n",
    "    rmc_runtime = rmc_end - rmc_start\n",
    "    \n",
    "    ohmc_runtime = ohmc_end-ohmc_start\n",
    "    rmc_anti_runtime = rmc_anti_end - rmc_anti_start\n",
    "    \n",
    "    ohmc_anti_runtime = ohmc_anti_end-ohmc_anti_start\n",
    "    \n",
    "    bs_err = np.abs(bs_price - theoretical_value)\n",
    "    rmc_err = np.abs(rmc_price - theoretical_value)\n",
    "    \n",
    "    ohmc_err = np.abs(ohmc_price - theoretical_value)\n",
    "    rmc_anti_err = np.abs(rmc_anti_price - theoretical_value)\n",
    "    \n",
    "    ohmc_anti_err = np.abs(ohmc_anti_price - theoretical_value)\n",
    "    \n",
    "    \n",
    "    result = {\"method\":[\"Black Scholes\",\"MC\",\"antithetic MC\",\"OHMC\",\"antithetic OHMC\"],\n",
    "              \"runtime\":[bs_runtime,rmc_runtime,rmc_anti_runtime,\n",
    "                         ohmc_runtime,ohmc_anti_runtime],\n",
    "              \"err\":[bs_err,rmc_err,rmc_anti_err,ohmc_err,ohmc_anti_err], \n",
    "              \"standard err estimate\":[bs_sde_estimate,rmc_sde_estimate,rmc_anti_sde_estimate,\n",
    "                                       ohmc_sde_estimate,ohmc_anti_sde_estimate],\n",
    "              \"root mean square error\":[bs_rmse, rmc_rmse, rmc_anti_rmse, ohmc_rmse, ohmc_anti_rmse],\n",
    "              \"n_trials\":[n_trials]*5,\n",
    "              \"n_steps\":[n_steps]*5}\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_trials: 100; n_steps: 150\n",
      "n_trials: 600; n_steps: 150\n",
      "n_trials: 1100; n_steps: 150\n",
      "n_trials: 1600; n_steps: 150\n",
      "n_trials: 2100; n_steps: 150\n",
      "n_trials: 2600; n_steps: 150\n",
      "n_trials: 3100; n_steps: 150\n",
      "n_trials: 3600; n_steps: 150\n",
      "n_trials: 4100; n_steps: 150\n",
      "n_trials: 4600; n_steps: 150\n",
      "n_trials: 5100; n_steps: 150\n",
      "n_trials: 5600; n_steps: 150\n",
      "n_trials: 6100; n_steps: 150\n",
      "n_trials: 6600; n_steps: 150\n",
      "n_trials: 7100; n_steps: 150\n",
      "n_trials: 7600; n_steps: 150\n",
      "n_trials: 8100; n_steps: 150\n",
      "n_trials: 8600; n_steps: 150\n",
      "n_trials: 9100; n_steps: 150\n",
      "n_trials: 9600; n_steps: 150\n"
     ]
    }
   ],
   "source": [
    "n_trials_list = np.arange(100,10000,500)\n",
    "# n_steps_list = np.arange(50,200,50)\n",
    "performance_df2 = pd.DataFrame()\n",
    "for n_trials in n_trials_list:\n",
    "#     for n_steps in n_steps_list:\n",
    "    print(\"n_trials: {}; n_steps: {}\".format(n_trials,n_steps))\n",
    "    new_df = pd.DataFrame(performance_comparison2(n_trials,n_steps,underlying_process=\"Heston model\"))\n",
    "    performance_df2 = pd.concat((performance_df2,new_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAELCAYAAAAhuwopAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd8VFXe/99nkkkmk0YSUiCVEjohkIA0FaSILgoIsrju\n+rjProhlVWzLiuta1va4lsXyY3XXtay7WCiC4lJEkCJIrwGSQBqQEDLpM2kz5/fHZIaUmdSZhJjz\nfr14kbn33HNPQrife75VSClRKBQKhaI1aDp7AQqFQqHoeijxUCgUCkWrUeKhUCgUilajxEOhUCgU\nrUaJh0KhUChajRIPhUKhULQaJR4KhUKhaDVKPBQKhULRapR4KBQKhaLVeHb2AlpLz549ZVxcXGcv\nQ6FQKLoU+/fvvySlDHXVfF1OPOLi4ti3b19nL0OhUCi6FEKITFfOp8xWCoVCoWg1SjwUCoVC0Wrc\nKh5CiBlCiFNCiDQhxBIH5wOFEOuEEIeFEMeFEL9253oUCoVC4Rrc5vMQQngAbwPTgBxgrxBirZTy\nRJ1h9wEnpJQ3CSFCgVNCiE+klFXuWpdC0ZWorq4mJyeHioqKzl6Kooug0+mIiopCq9W69T7udJiP\nAdKklGcAhBArgFlAXfGQgL8QQgB+gAGoceOaFIouRU5ODv7+/sTFxWH9b6JQOEdKSUFBATk5OfTp\n08et93Kn2SoSyK7zOaf2WF3eAgYD54GjwINSSosb16RQdCkqKioICQlRwqFoEUIIQkJCOmSn2tkO\n8+uBQ0BvIBF4SwgR0HCQEGKhEGKfEGJffn5+R69RoehUlHAoWkNH/b64UzzOAdF1PkfVHqvLr4FV\n0koacBYY1HAiKeW7UspkKWVyaGjbclwyLpXzwc6zlFZUt+l6hUKhUFzGneKxF4gXQvQRQngBC4C1\nDcZkAVMAhBDhwEDgjDsWk3KhhKfXnSDLYHTH9AqFog34+fm16/pJkyYRExODlNJ+bPbs2fXmPX36\nNDfeeCPx8fGMGjWK+fPnk5eX1677KtwoHlLKGuB+YAOQAnwmpTwuhFgkhFhUO+w5YLwQ4ijwLfB7\nKeUld6wnPFAHwMWSSndMr1D85JFSYrF0nkvS2f179OjBzp07ASgqKuLChQv2cxUVFfzsZz/jnnvu\nITU1lQMHDnDvvfeizN/tx60+DynleinlACllPynl87XHlkspl9d+fV5KOV1KOVxKOUxK+S93rSU8\nwCoeuSUq5FGhaCkZGRkMHDiQO+64g2HDhpGdnc3GjRsZN24co0aN4tZbb6WsrAyA9evXM2jQIJKS\nknjggQeYOXMmAE8//TR/+ctf7HMOGzaMjIyMevcpKytjypQpjBo1iuHDh/Pll186vX9DFixYwIoV\nKwBYtWoVt9xyi/3cv//9b8aNG8dNN91kPzZp0iSGDRvmmh9QN6bL1bZqK2H+3ggBucVKPBRdk2fW\nHefE+RKXzjmkdwB/umlok2NSU1P58MMPGTt2LJcuXeLPf/4zmzdvxtfXl5dffpnXXnuNxx9/nLvv\nvpvvv/+ePn36cNttt7VqHTqdjtWrVxMQEMClS5cYO3YsN998c6P7O2LKlCncddddmM1mVqxYwbvv\nvstzzz0HwLFjx0hKSmrVWhQto9uIh9ZDQ4ivNxdLlXgoFK0hNjbW/uDevXs3J06cYMKECQBUVVUx\nbtw4Tp48Sd++fe25Bbfddhvvvvtui+8hpeSJJ57g+++/R6PRcO7cObtfou79HeHh4cHEiRNZsWIF\nJpMJVXW7Y+g24gEQHuCtdh6KLktzOwR34evra/9aSsm0adP4z3/+U2/MoUOHnF7v6elZz1fhKAfh\nk08+IT8/n/3796PVaomLi7OPq3t/ZyxYsIA5c+bw9NNP1zs+dOhQtm3b1uz1itbT2XkeHUpEgI48\n5TBXKNrM2LFj2blzJ2lpaQCUl5dz+vRpBg4cyJkzZ+y+jE8//dR+TVxcHAcOHADgwIEDnD17ttG8\nxcXFhIWFodVq+e6778jMbF318Kuvvpo//OEPjcxlv/jFL9i1axdff/21/dj333/PsWPHWjW/ojHd\nSjzCAnTkKYe5QtFmQkND+eCDD7jttttISEiwm6x8fHx45513mDFjBklJSfj7+xMYGAjA3LlzMRgM\nDB06lLfeeosBAwY0mvf2229n3759DB8+nI8++ohBgxqlezWJEIJHH32Unj171jvu4+PDV199xZtv\nvkl8fDxDhgzhnXfeoa35YorLiLrx0V2B5ORk2dZmUH/dnMrrm09z+s834OXZrXRT0UVJSUlh8ODB\nnb2MFlFWVoafnx9SSu677z7i4+NZvHhxZy+rW+Lo90YIsV9Kmeyqe3SrJ2h4gDeAcporFG7gvffe\nIzExkaFDh1JcXMzdd9/d2UtSuJHu5TCvTRTMK6kkKkjfyatRKH5aLF68WO00uhHda+fhbxMPtfNQ\nKBSK9tCtxCOiduehwnUVCoWifXQr8QjSa/Hy0JCnfB4KhULRLrqVeAghCAvwJk/tPBQKhaJddCvx\nAJUoqFC4mzVr1nDixOVu00899RSbN28G4I033sBovNwWobUl2bdu3cquXbvsn5cvX85HH33UzhW7\nhu5WXr7biUe4ShRUKNxKQ/F49tlnmTp1KtBYPFpLQ/FYtGgRd9xxR6vnUeXl248SD4VC0SSzZ88m\nKSmJoUOH1it26Ofnx9KlSxkxYgRjx44lLy+PXbt2sXbtWh577DESExNJT0/nzjvv5IsvvmDZsmWc\nP3+eyZMnM3nyZPs8DecAyM/PZ+7cuYwePZrRo0ezc+dOMjIyWL58Oa+//jqJiYls3769Xrn3tLQ0\npk6dyogRIxg1ahTp6en1vg9VXt61dKs8D7AmCpZXmSmtqMZfp+3s5SgULeebJZB71LVzRgyHG15q\ncsj7779PcHAwJpOJ0aNHM3fuXEJCQigvL2fs2LE8//zzPP7447z33ns8+eST3HzzzcycOZN58+bV\nm+eBBx7gtdde47vvvrOXEXE2x4MPPsjixYuZOHEiWVlZXH/99aSkpLBo0SL8/Px49NFHAfj222/t\n899+++0sWbKEOXPmUFFR4fDNXpWXdx1uFQ8hxAzgr4AH8Hcp5UsNzj8G3F5nLYOBUCmlwV1riqiT\nKKjEQ6FonmXLlrF69WoAsrOzSU1NJSQkBC8vL/sbeVJSEps2bWr13M7m2Lx5cz3TV0lJiX1X4IjS\n0lLOnTvHnDlzAOsD3BGqvLzrcJt4CCE8gLeBaUAOsFcIsVZKaf+NkFK+ArxSO/4mYLE7hQMgrE6i\nYP+w9jm4FIoOpZkdgjvYunUrmzdv5ocffkCv1zNp0iR7qXStVosQArA+9Gpqalo9v7M5LBYLu3fv\ndioCbUWVl3cd7vR5jAHSpJRnpJRVwApgVhPjbwP+08R5l3B556H8HgpFcxQXFxMUFIRer+fkyZPs\n3r272Wv8/f0pLS1t9bm6TJ8+nTfffNP+2fZAd3a9v78/UVFRrFmzBoDKyspmHfOqvHz7cKd4RAJ1\nPUI5tccaIYTQAzOAlW5cD3C5OKLqZa5QNM+MGTOoqalh8ODBLFmypEmTi40FCxbwyiuvMHLkyEZO\n64ULFzJjxox6DnNHLFu2jH379pGQkMCQIUNYvnw5ADfddBOrV6+2O8zr8vHHH7Ns2TISEhIYP348\nubm5Td5DlZdvH24ryS6EmAfMkFL+tvbzr4CrpJT3Oxj7c+CXUsqbGp6rPb8QWAgQExOT1Folb8jw\npzdwy8hInpnVOVEKCkVL6Uol2X9KdPXy8l29JPs5ILrO56jaY45YQBMmKynlu1LKZCllsitUViUK\nKhSKplDl5ZvHndFWe4F4IUQfrKKxAPhFw0FCiEDgWuCXblxLPcIDdMpspVAonKLKyzeP23YeUsoa\n4H5gA5ACfCalPC6EWCSEWFRn6Bxgo5Sy3F1raUh4gI6LSjwUCoWizbg1z0NKuR5Y3+DY8gafPwA+\ncOc6GhIe4M3F0kosFolGIzry1gqFQvGToNuVJwFruG6NRVJQXtXZS1EoFIouSbcUjzDVUVChUCja\nRbcUj+6YKHiuyMSst3eSVdD2iqaK7ktOTg6zZs0iPj6efv368eCDD1JVVcXWrVvt5UVs2AohQtcr\nM65oOd1SPDojUdBskWxPzcddeTXNsfesgcPZRfxrT/tyZBTdDyklt9xyC7NnzyY1NZXTp09TVlbG\n0qVLW3R9Vyozrmg53VI8Qv28EYIO7Sj43cmL/OofP3I4p7jD7lmXLIN1x7HqQA7V5s7rY6DoemzZ\nsgWdTsevf/1rwFqD6vXXX+f9999vUW+OrlRmXNFyul1JdgBPDw09/bw7NFHwzCVrRdD0i2UkRvfo\nsPvayKw1V10qq+K7kxeZPjSiw9egaB8v//gyJw0nXTrnoOBB/H7M75scc/z48UalwAMCAoiJiSEt\nLY3t27eTmJhoP5eVlVXPlNWVyowrWk633HmANcu8I81Wtjf/zIIOS2dpcP9yRsX0INTfm8/25XTK\nGhQ/Ta6++moOHTpk/2PrXWGjK5UZV7ScbrnzAGuiYE5hxzmPsw0mADINneOwziwwcs2AUEb3Cebv\n289ysbTCHnWm6Bo0t0NwF0OGDLE7wG2UlJSQlZVF//792bhxY7NzdJUy44qW0213HuEB3h0abZVt\n33l0vHiYqsxcLK0kNljPrUnRmC2S1QeclRlTKOozZcoUjEYjH330EQBms5lHHnmEO++8E71e36I5\nukqZcUXL6bbiERGgo9BYTWWN2e33slgkOYXWnUdWJ+w8smt3WDEhevqH+ZEUG8Rn+7I7LfJL0bUQ\nQrB69Wo+//xz4uPjGTBgADqdjhdeeKFVc3SFMuOKltOtzVYAF0sqiQ5u2dtTW8krraDKbCE62Ids\ng4mSimoCOrAFrm23E1P7fc5PjuL3K49yIKuIpNigDluHousSHR3NunXrGh2fNGkSkyZNqnfsgw8+\nsH+9detWh/PVbSk7aNAg/vvf/7pimYoOpNvuPMI7MFHQ5u+Y2N/6NtXRiXq23U5siLUF5s8SeuOj\n9eDzfdlNXaZQKBRO6b7i0YGJgraH98T+1i17R/s9sgrK8ff2JEhv3e34eXvys4RerDt8HmNV6/tO\nKxQKRbcVj4gA287D/bke2QYjQsD4fiEAZBo6Nlw302AkOliPEJcrCM9Pjqa8ysz6o0236lQoFApH\ndFvxCPTR4uWp6RizVaGRiAAdQb5e9PTz6nizVYGR2JD6fp3RcUHEhej5TJmuFApFG+i24iGEsCYK\ndkCJkmyDkegg68M7JljfoWYrc22kV0wD8RBCcGtyND+eNXD2UuckLioUiq6LW8VDCDFDCHFKCJEm\nhFjiZMwkIcQhIcRxIUSHZgtZe5l3jMPcFtEVF+LboVnmuSXWSK8YBxFlc0dFoRHwxX61+1AoFK3D\nbeIhhPAA3gZuAIYAtwkhhjQY0wN4B7hZSjkUuNVd63FEWAckClZUm8ktqSA62Aew5lpcKKnokPwS\nuFwOJTbYt9G5iEAd1w4I5Yv9OZgtKudD4RrWrFnDiRMn7J+feuopNm/eDMAbb7xRr5hi3dLsLWHr\n1q3s2rXL/nn58uX25MWW4Ky0vG1uVV6+5bhz5zEGSJNSnpFSVgErgFkNxvwCWCWlzAKQUl5043oa\nYd15VLo1We5ckTVM1/bmHxuiR8rL4bvuJtsepus4l2V+cjR5JZV8n6pKYCtcQ0PxePbZZ5k6dSrQ\nWDxaS0PxWLRoEXfccUeLrm1vaXlQ5eXr4k7xiATq2kNyao/VZQAQJITYKoTYL4Ro2W+BiwgP0GGq\nNlNS4b5wVVuYrs1sFVO7A8jqoIirzAIjnhpBr0DHdaymDA4n2NdL5XwonDJ79mySkpIYOnQo7777\nrv24n58fS5cuZcSIEYwdO5a8vDx27drF2rVreeyxx0hMTCQ9Pd3+9r5s2TLOnz/P5MmTmTx5sn2e\nhnMA5OfnM3fuXEaPHs3o0aPZuXMnGRkZLF++nNdff53ExES2b9/O008/zV/+8hcA0tLSmDp1KiNG\njGDUqFGkp6fX+z7aW1oeVHn5unR2hrknkARMAXyAH4QQu6WUp+sOEkIsBBYCxMTEuOzmtkTBiyUV\nBPq4J+M7x1A/u9u2A+gop3mmwUhkkA+eHo7fE7w8NcxOjOTj3RkYyqsI9vXqkHUpWk/uCy9QmeLa\nkuzegwcR8cQTTY55//33CQ4OxmQyMXr0aObOnUtISAjl5eWMHTuW559/nscff5z33nuPJ598kptv\nvpmZM2cyb968evM88MADvPbaa3z33Xf2MiXO5njwwQdZvHgxEydOJCsri+uvv56UlBQWLVqEn58f\njz76KADffvutff7bb7+dJUuWMGfOHCoqKrBY6vetaa60PKDKy7cCd4rHOSC6zueo2mN1yQEKpJTl\nQLkQ4ntgBFBPPKSU7wLvAiQnJ7vMxhTufzlRMD7c31XT1iPLYMTLU0Oon/VeIb5e+Hl7dph4ZBUY\nHTrL6zJ/dBTv7zzLmoPn+N+JfTpkXYquw7Jly1i9ejUA2dnZpKamEhISgpeXl/3BmpSUxKZNm1o9\nt7M5Nm/eXM/0VVJSUq+kSUNKS0s5d+4cc+bMAUCna1vF6KuvvpqvvvrK/vnOO++sd16Vl7+MO8Vj\nLxAvhOiDVTQWYPVx1OVL4C0hhCfgBVwFvO7GNdXjci9z9yUKZhtMRAf5oNFYE/SEELXhuh1jtsoy\nGBkR3avJMYMiAkiICuSzfdn8ekJcvWRCxZVDczsEd7B161Y2b97MDz/8gF6vZ9KkSVRUWINMtFqt\n/XfFw8ODmprWm3+dzWGxWNi9e3ebRcARzZWW//HHH1s0jyovb8VtPg8pZQ1wP7ABSAE+k1IeF0Is\nEkIsqh2TAvwXOAL8CPxdStlhtZjDA9xf3yqrNru7LrEh+g7p61FsrKbYVO0w0qohtyZHczK3lGPn\nSty+LkXXobi4mKCgIPR6PSdPnmT37t3NXuPv709paWmrz9Vl+vTpvPnmm/bPhw4davJ6f39/oqKi\nWLNmDQCVlZWN/BiuKC0Pqry8DbfmeUgp10spB0gp+0kpn689tlxKubzOmFeklEOklMOklG+4cz0N\n0Wk9CPTRui1RUEpJtqGx2SgmRE+OweT28FhbGZSWVA2+eURvvD01KuNcUY8ZM2ZQU1PD4MGDWbJk\nCWPHjm32mgULFvDKK68wcuTIRk7rhQsXMmPGjHoOc0csW7aMffv2kZCQwJAhQ1i+3PrIuOmmm1i9\nerXdYV6Xjz/+mGXLlpGQkMD48ePJza1fescVpeVt86jy8iC6Wk+H5ORkuW/fPpfNN/31bcSF+PLu\nHckum9NGkbGKxGc3sfTGwdx1TV/78X/vyeKJ1UfZueQ6Inv4uPy+NtYdPs/v/nOQbx68msG9Apod\n/+CKg3x38iI/Lp2KTuvhtnUpWk5KSgqDBw/u7GUouhiOfm+EEPullC570HXb8iQ2wt2YZW7L5XBk\ntgL39zPPahDp1Rzzk6Mpqahhw3FVLFGhUDRNtxcPW6KgO7ic41F/d9FR4bpZBUZ6+nnh692yuIhx\nfUOI7OHD5/ty3LouhULR9en24hEeoCO/rNIt/gdb+9eGO49egT5oPYTbxSPTUN7iXQeARiO4NTmK\nnemX7Jnpis6nq5mWFZ1LR/2+KPEI1GG2SArKXL/7yDIY6aHXNmo566ERRAfp3Z5lnm0w2bsHtpR5\nSVEArDygdh9XAjqdjoKCAiUgihYhpaSgoMClIc7O6OwM806nbqJgWIBrf+COIq1sxIS4tzR7ZY2Z\n88WmVu08AKKC9Ezo15PP9+XwwHXx9vwURecQFRVFTk7OT7Y+ksL16HQ6oqKi3H6fbi8e7kwUzDYY\nGdo70OG52GA9+zMKkVK6JSkvp9CElC13ltfl1uQoHlxxiB/OFDChf8/mL1C4Da1WS58+KutfceWh\nzFa1uw1X9zI3WyTnikxOcyxiQnwprayh0Fjt0vvasHUrdFZNtymuHxpBgM5T5XwoFAqndHvx6Onn\njUZYiyO6ktySCqrNslGklY242od6hpvCde1hum0QD53Wg1mJkXxzLJdiN4mbQqHo2nR78fDQCEL9\nvV2eZZ7dTI6FbUfgrn7mmQVGfLQe9oKMrWV+cjRVNRbWHjnv4pUpFIqfAt1ePMCa6+Fqs5U9xyPI\nsXhEBekRwn25Hlm1Ybpt9acMiwxgUIS/6vOhUCgcosQDq9/joosd5jkGIxoBvZ2UH9FpPYgI0Nnr\nT7mazAJjm0xWNoQQzE+O5khOMSkXVLFEhUJRHyUeWMXDHTuPXoE+eHk6/xHHBOvdYraSUpJlMBLb\nhkirusweGYnWQ6iMc4VC0QglHljDdYtN1VRUm102Z3ahyamz3EZciC8ZbhCPi6WVVNZY2hRpVZdg\nXy+mDQln9cEcqmoszV+gUCi6DUo8gLDaREFXFkjMMhid+jtsxITouVRWSXmla3uo2/woLSnF3hy3\nJkdTaKzm25S8ds+lUCh+OijxwPWJgqYqM/mllc0m6NkjrlxcR8o2X2tLkzjimvhQIgJ0KudDoVDU\nw63iIYSYIYQ4JYRIE0IscXB+khCiWAhxqPbPU+5cjzNcnSiY46QgYkNsHf5cHXGVVVCORuCSXiEe\nGsHcpEi2nc53W9MshULR9XCbeAghPIC3gRuAIcBtQoghDoZul1Im1v551l3raQqbeLgqUdBZNd2G\nxNh3Hq6NuMpsgbO+NdyaFI1FqmKJCoXiMu7ceYwB0qSUZ6SUVcAKYJYb79dmAnSe6LQal71ZZxU4\n7uPRkEAfLT30WpfvPDILjO12ltclrqcvY/oE8/m+bFXdVaFQAO4Vj0igrqE8p/ZYQ8YLIY4IIb4R\nQgx143qcIoRwaaJgdqEJnVbTouzu2BBfl4tHtsG14gFwa1IUGQVGDucUu3RehULRNelsh/kBIEZK\nmQC8CaxxNEgIsVAIsU8Isc9dpanDXJgoaIu0akl2d2yw3qWJgmWVNRSUVxET3H5neV1GxwUDcDqv\n1KXzKhSKrok7xeMcEF3nc1TtMTtSyhIpZVnt1+sBrRCiUQ1wKeW7UspkKWVyaGioWxbr0p1HE308\nGhIboud8UQXVZtfkUdj6orelFHtTRAb54KERbqvFpVAouhbuFI+9QLwQoo8QwgtYAKytO0AIESFq\nX8+FEGNq11PgxjU5JSJQR15JRbtt+lJKsg3GFudYxATrreXbC03tuq+NbEPbS7E3hdZDQ+8eOjJV\ne1qFQoEbm0FJKWuEEPcDGwAP4H0p5XEhxKLa88uBecA9QogawAQskJ3kkQ3z96ayxkKxqZoeeq82\nz1NorKa8ytxi8bDlYmQajMT1bL+pyeY/aU9dK2fEBvuS5aYS8gqFomvh1k6Ctaao9Q2OLa/z9VvA\nW+5cQ0upmyjYHvG4XE23ZTkWtr4eVnNT+01ymQYjQQ76pruCmBA93xy94PJ5FQpF16OzHeZXDK5K\nFMxuZROmUH9vfLQeLou4yipoub+ltcQE6yk0VlNSoRpEKRTdHSUetUQE2HYe7ROP5vp4NEQIQUyw\n3nXiYTAS44KyJI6wVelVTnOFQqHEo5ZQW3HEdiYK5hQaCfH1wte75RbBmBC9S7LMq80WzhWZ2l2K\n3RkxdhObEg+ForujxKMWndaDIL2WvNL27zyiWvnwjg3Wk2UwYrG0L1bgfJEJs0W6xVkOdZ37ymmu\nUHR3lHjUITxAR25x+xIFsw2mVvscYkP0VFRbuFjavnvbI63ctPPw8/YkxNdLma0UCoUSj7qEB+ja\n5fOoMVs4X2RqcaSVDfsbfTvDYLPclONRl5gQ1/lnFApF10WJRx0i2ikeF4orqLHINu08gHYn4GUZ\njHh5agj317VrnqaIqTWxKRSK7o0SjzqEB3hzqaySmjaWCmlpKfaG9O7hmtIfmQXlRAf5oNE0X1Or\nrcQG67lQbFJtaRWKbo4SjzqEB+qwSLhUVtWm6+05Hq0UD62HhsgePi7YeZhc0j2wKWJCfLHIyw2v\nFApF90SJRx1s5p62JgpmG0x4aAS9AltvNooN0ber9IeUkqyCcrc5y224ysSmUCi6Nko86nC5REnb\nxCPLYKR3Dx2eHq3/scaG6Mloh9mqoLyK8iqzW53loBIFFQqFFSUedQgLqE0UbOvOo9DY4szyhsQG\n+1JsqqbY2LbSH+4O07Xh6nIqCoWia6LEow49fb3x0Ii2i0cr+ng0xJ693cYEPHeVYm+IrZyKirhS\nKLo3SjzqoNEIwvy925QoaKyq4VJZVasjrWzEtrP0R2aBESEgqo07n9YQHeyacioKhaLrosSjAW1N\nFMw2WJs5tVU8bDuWtr7RZxrKiQjQodN6tOn61hAbYt15dFLrFYVCcQWgxKMB4QHebRSP1vXxaIje\ny5NQf28yLrXdbNVW4WqWqnLI3AVHPgNzjcvKqSgUiq5Ls6VfhRAaYJ6U8rPWTi6EmAH8FWsnwb9L\nKV9yMm408APWToJftPY+riQiQMeu9NZ3ws1qY45HXeJC9G0Ogc0sMHLtABf0d6+ugLxjcP6g9c+5\nA3DpFMjapECtDzHBV9nvaeuDolAouhfNioeU0iKEeBxolXgIITyAt4FpQA6wVwixVkp5wsG4l4GN\nrZnfXYQH6iitqMFYVYPe6/KPR0pJWlEa8UHxDq/LLjSi9/Ig2LftXQhjgn3ZmXap1deZqsxcLK1s\nvbPcXA0XT1wWifMHrZ8tNdbz+p4QOQqGzIJeI+DzOyFrN7FJ1wHWjPYxfYJbvV6FQtH1aWnTic1C\niEeBTwG7XUVKaWjimjFAmpTyDIAQYgUwCzjRYNzvgJXA6JYu2p3YEgXzSirp0/Pyj2dt+lqe3Pkk\nH9/wMYlhiY2us0VaCdH20iCxIXpWHqigotrcKt+FvQFVc7ueS2mQ8+Nlocg9CuZa05OuB/QeCeMf\nsP4dOQoCIqHu99N7JGT/SOQUHzTisqlOoVB0P1oqHj+v/fu+Osck0LeJayKB7Dqfc4Cr6g4QQkQC\nc4DJXCHiUTdRsE/Py6U+Pjtt3XhtyNjgRDxM7fY52HYO2QYj8eH+Lb7ucjXdJkqT5OyHv1t3DHj5\nQa9EGHPj1YPOAAAgAElEQVSXVSR6j4SgPvWFwhHRY2DPcrxkFb0C219ORaFQdF1a6vP4pZRypxvu\n/wbw+1rTWFNrWAgsBIiJiXHDMi4T7iBRMLUwlSP5R/D28GZT5iYeG/0YGnE51kBKSZbByIT+Pdt1\nb5u/JLOgdeJhK+XeZAfBnB+tf//2W+g9CjRtiJWIvgp2LYMLh4lVpdkVim5Ns08QKaUFeKsNc58D\nout8jqo9VpdkYIUQIgOYB7wjhJjtYA3vSimTpZTJoaEucAo3QbiDXuarUleh1Wh5aNRD5BnzOJJ/\npN41BeVVmKrNRAe3LdLKhm3nkNHKGldZBiP+Ok966LXOB+UeA99QiEpum3CAdecBkL3HHq6rUCi6\nJy19inwrhJgrWmfQ3wvECyH6CCG8gAXA2roDpJR9pJRxUso44AvgXinlmlbcw+X4eXui9/KwJwpW\nmitZd2YdU2KmMKv/LLQaLZsyN9W7xhWRVgBBei3+Os9WP5QzC1rgb8k7BuHD2rU+/MKs5q3sPcQE\n+2Ior6K0om3lVBQKRdempeJxN9Zoq0ohRIkQolQIUdLUBVLKGuB+YAOQAnwmpTwuhFgkhFjUrlW7\nESGEtSlUbS/zzZmbKa4s5pb4W/D38md87/FsytxUL0Euu6UO6xbcuy3moGyDselIK3MNXEyB8KHt\nWh9gNV1l/0hs7S5Lma4Uiu5JS8UjELgT+LOUMgAYijUEt0mklOullAOklP2klM/XHlsupVzuYOyd\nnZ3jYSMswJu8Yqt4rExdSaRfJFf1svr6p8dN50L5BY5dOmYffzlBsP1JerHBvq3aeZgtkuxCIzHB\nTTjLDenWqKqI4e1eH9FjoPwi/bTWkGJlulIouictFY+3gbHAbbWfS2mbH6TTqMnPp/jLL5FVzTd6\nigjQkVtSQWZJJntz9zI3fq7dQX5t1LV4ajzrma6yDSZ6+nnj49X+0iAxIXpyCo0t7mZ4odhEtVk2\nvfPIPWr9u71mK7DuPIBoo1U8lXgoFN2TlorHVVLK+4AKACllIdD2bLhOwLh/P+d/v4SK06nNjg0P\n1HGxpJJVp1fhITyY1X+W/VygdyBje41lY+ZGu+kqy2Akpp3OchuxwXqqzZILxS0rkdIif0veMdBo\noeeA9i8wbDB4+aPP3U+QXqvMVgpFN6Wl4lFdmwkuAYQQoUCXamKtG2Y12VQcO9rs2HB/HVXmatak\nfcnVUVcTpg+rd3567HTOlZ3jhMGa75hd6Lq6UraIq5Y+lLNa0scj7ziEDgRPF+i9xgOikqxO8xBf\nVV1XoeimtFQ8lgGrgTAhxPPADuAFt63KDWgje+MRHIzpSPPiERGow9M/BUNlAfPi5zU6f13MdXgK\nTzZlbKLabOF8kcllTZhiW9nXI9NgROsh6N2jiZ1PrgsireoSfRXkHWdAD+UwVyi6Ky0SDynlJ8Dj\nwIvABWC2lPJzdy7M1Qgh8Bk+nIqjR5odGx7gjbbHXgK9ejIhckKj84HegYzpNYaNmRs5X2jCIl3j\nLAerv8XLU9PiNq9ZBUaigvR4aJyE6RoNUHreNZFWNqLHgLQwRnuG80Umqmq61CZU0QZMNSb+36H/\nx4WyC529FMUVQouzxaSUJ6WUb0sp35JSprhzUe5ClzCcyrR0zGVNv9ULzyI8fE8zPGAqnhrHSfjT\nY6eTXZrNrhzrTsZVZiuNRhAd5NNys1VzpdjzaqPCIly484hMBgSDa05ikXCuyOS6uRVXHBU1FTyw\n5QHeOfwOy480CpRUdFO6VT8Pn4QEkJKK48ebHLcz7xsAIj0mOR1zXcx1eAgPtuZsBmh3dnldYkN8\nW1w3KrOgvOmyJLm14uFKs5VPDwgbTFS5VThVxNVPl0pzJQ999xB7LuyhX2A/NmRswFSjXhYU3Uw8\ndMOsD9CmTFdmi5m1Z9agqRiA0RTgdFyQLojkiGSOFn2PpwZ6BbpOPGKC9WQWlDfbqa/IWEVJRU3T\nYbp5x8A3zJod7kqixxBw6SACC1mtLKei6BpUmatY/N1idp7fyWNJf6Ts/M2UV5ezJWtLZy9NcQXQ\nrcTDMygIbUxMk07zXed3kVueS5DlanuioDOmx06n1HyB8J5Fzn0ObSAuRI+xysylsqZzUmymrWbN\nVq40WdmIvgpNZQlDPC8op/lPkGpzNY9sfYTt57bz1LinOJ02mLSsUDzMwaxJ+7Kzl6e4AuhW4gHg\nM2wYpibCdVemriRYF0ycz2hym2lHe13MdSAFPkFNm8Faiy1ct7kw2Mul2J2Ih7kGLp50rbPcRm2y\n4BS/DFWa/SdGtaWaR7c9ytacrTx51ZMkBMzgkz2Z9A8LwGhIZM+F3eSW53b2MhWdTLcTD13CcGrO\nX6AmP7/RuUumS2zL3sbN/W6mV4AfeSVN9+ju6dMTUdmXcs8DLl1jTMjl0uxN0WyCYEGqtSxJuAvK\nkjQkuC/oQ7jKM63FkWGKK59qSzW///73bMnewpIxS5g/cD7PfXUCf52W/9w1llAxAYnkq/SvOnup\nik6m24mHT0ICAKajxxqdW5O2hhpZwy3xtxAeoKOgvJLqJsqElFXWYCoaRpnlHOlF6S5bY1SQD0I0\nLx6ZBeWE+nvXa5dbj7zaHZE7zFZCQPRVDKxOIctgbNY/o7jyqbHU8Iftf7D2rEl+jNsH3863KRfZ\nkXaJh6bGE+rvze+uHkuNMY4VJ1erf/NuTrcTD93gweDhgamB01xKyarUVSSFJ9EnsA8RgTqkhPxS\n57uPbIORmtKhCAQbM1zXgt3b04PegT72Jk/OyKptfeuU3KPWsiQhjvuut5voMfSszEJXXUh+WdO7\nNMWVjdli5okdT7AhYwOPJD3CHUPvoKrGwvPrU+gX6ssvx8YCcMuoKHyrriLPlNWor42ie9HtxEPj\n44P3gAFUNHCa783dS3ZpNnPj5wKXOwo25ffIMhiRNQEM6JHAxkzXiQdY/RjN+RKyCoxNh+nmHYfQ\nQa4pS+KIKGtzqJGaK990lVVgZM47O3llw8kmXwi6I2aLmT/u/CPfnP2GB0c9yJ3D7gTgw10ZnL1U\nzpMzh6D1sD4qvDw1LBw1B2nx5O+HulSesMLFdDvxAPAZPhzT0aP1tt1fpH6Bv5c/02KtleZtHQUv\nNiEetlLsM+Kmk1aUxpniMy5bY2yIvskHcmWNmQslFXb/iEPcFWllo/dIpMaTJM3pKzriqsZs4aFP\nD3L8fAnvbE1nwstbeGL1Uc5eUiHGFmnhT7v+xLoz6/jdyN/x2+G/BeBSWSXLvk1l0sBQJg+sH+Z9\nx9hBeJgS2H5+E1Xm5qtUK36adE/xSBiOpaSE6sxMAIoqiticuZmZfWei87SKhk08cpsI1802GPH3\n9mRm/+kAbMrY5HRsa4kJ9qWgvIqyyhon9zYhZRORVuUFUHrBPZFWNrz0yIgEkjSpV3TE1dvfpXMg\nq4hX5iXw7cPXMndUFF/sz+G6V7dyz7/2czi7qLOX2ClYpIVnf3iWL9O/5N4R97IwYaH93GubTmOs\nNvPkz4Y0uk6n9eCmvjdjFkb+efDrjlyy4grCreIhhJghhDglhEgTQixxcH6WEOKIEOKQEGKfEGKi\nO9djQzfcGn1kc5p/deYrqi3VdpMVQLDeC62HIK8pn0ehiahgPRG+ESSGJjZqT9se7AUSnfg9spuL\ntMpzQ2a5AzTRV5GoSSfnUrFb79NWDmYVsmxLKrMTezMrMZK+oX68eMtwdvx+Mvdc248daZeY9fZO\nFrz7A1tPXew2TmApJX/e/WdWpq5kYcJCFo243Nwz5UIJK37M4ldjY+kf5ufw+sevnQk1gXx87Iro\n36boBNwmHrUl3N8GbgCGALcJIRq+xnwLjJBSJgL/C/zdXeupi3e/fggfH0xHjyClZGXqSob3HM7A\n4IH2MRqNIMxf12SiYN0+HtPjpnOq8BSZJZkuWaNNFJyZg2yi4rSDYAeJB9Fj0FGF5mLj6LXOpryy\nhsWfHiIiQMczs+r/HML8dTw+YxA//GEKS28cTMYlI3f+cy83/HU7qw/mNBll19WRUvL8nuf5/PTn\n/GbYb7g/8X6EEPZzz647QYCPloemOg+0CPDxJjHoOoo4xo4zZztq6YorCHfuPMYAaVLKM1LKKmAF\nMKvuACllmbz8qudLbb8QdyM8PdENHULFkaMczj9MWlFavV2HjfAAb3sv84ZIKck2GO3VdKfGTAVw\n2e4jtplcj0yDEb2XBz39nDjDc4+BXzj4hbpkPU6pTRYMLzrs3vu0gee+OkGmwchr80cQ6KN1OMbP\n25O7runL949P5i+3jsBskSz+9DCTXtnK+zvOUu7EbNhVkVLy8t6X+fTUp9w59E4eHPWgXTgANp7I\n44czBTw8bQA99E0HWjw28XaEsPB/O/7t7mUrrkDcKR6RQHadzzm1x+ohhJgjhDgJfI1199Eh+AxP\noOLECValfI6Ppw8z+sxoNCY8QOfU55FfWklljcXusO7l14uEngkuC9n112kJ8fVymmWeXRumW/c/\nfj3yXNzDwxmBkZR6hzOoJsWpf6Yz+O+xXFbszWbRtf24qm9Is+O9PDXMS4piw0PX8I//SaZ3Dx3P\nfnWC8S9t4bWNp7j0EwhFllLyyr5X+CTlE345+Jc8nPRwvd+fyhozL6xPIT7Mj1+MiWl2voSwgfTU\n9ifNtJVTuSXuXLriCqTTHeZSytVSykHAbOA5R2OEEAtrfSL78h1khrcFn4ThyKoqUvb+lxv73Iiv\ntrH5JzxA5zTLPLuwtq5UnT4e0+Omk2JIIbs02+E1rSUmRN+E2aqJHA9zNeS7qSyJA0pCRzFKk3rF\nhOteLKngD6uOMCwygMVTW9d6V6MRTBkczueLxrPynnFc1SeYZVvSmPDSFp5cc5RiU7WbVu1+3jjw\nBh+f+JhfDPoFj49+vNGLxz93ZpBZYOSPM4fg6dGyR8Mvh83FQ5fL/2351h1LVlzBuFM8zgHRdT5H\n1R5ziJTye6CvEKKng3PvSimTpZTJoaGuMcPohlszzaOzKxyarMDaUbCsssbhG3W2wVqWum5Rwqmx\nLjZdBTsWD4tFkmUwOo+0upQK5iqIcENZEgeIqKuIFAXk5bguy76tWCySR784gqnazBs/H4mXZ9t/\nxZNig3n3jmQ2P3wtsxMjWfFjNvf8a3+X9Iccv3Sc94+9z7wB81gyZkkj4cgvreStLWlMGRTGNQNa\n/n9s3sCb0ODJzrwNZKjQ526Fk7oWLmEvEC+E6INVNBYAv6g7QAjRH0iXUkohxCjAGyhw45rsaCN7\nU+7rwagCf4b1dGzesSUK5pVU4BdaP+rEVlcqKuhyKfZIv0iGhgxlU8Ym/ndY+y1wMSG+rD18nsoa\nM96eHvbjF+0mM2fO8tqyJB1htgICBkyA3VCTuQfGjGzTHNXmagorCymsKMRQYaCwopDCykIKTAX2\n42Zp5uGkh+kT2MfpPB/9kMH3p/N5bvYwp5FCraV/mB8vz0tgbJgXi9en88c1x3jxluHOTYZXIJ+f\ntppnG5qqbLy68RSVNWaW/mxwq+YN9A5kYuS1bKvZzdvfneKVW0e5asmKKxy3iYeUskYIcT+wAfAA\n3pdSHhdCLKo9vxyYC9whhKgGTMDPZQfFSp40nORUhIUheVqnDwFbrkdeSQX9GohHtsFIeIA3Oq1H\nvePT46bz+v7XOV92nt5+vdu1xthgPRYJOYWmevdvtiBi3lHw8IKebipL0gC/2JGY8MInb7/TMWaL\nmXVn1pFZktlIIAwVBkqrSh1epxEaenj3IFgXzIXyCzz+/eP8+8Z/o/Vo7AA/nVfKC9+c5LpBYfzy\nquZt9q2hKjubQffN55lf/54/7YW+ob4svKafS+/hLsqry1l/dj3Xx12Pv5d/o/PHzhXz6b5sfjOh\nD31DWy+48wfO4ftz3/Jl6rc8VDSYyB6u622juHJx584DKeV6YH2DY8vrfP0y8LI71+CMlakrEZGe\njNp+CXNZOR5+jn0eYBWPhmTVibSqy7SYaby+/3U2ZW7if4b+T7vWaDNLZRUY64mHLUzXaWmS3GMQ\nOhAcPGDdgoeWdM8BRBQ7j7j6z8n/8PLel/EQHvTw7kGQLohgXTCDgwcTpAuyfvYOth8P1lm/DvAK\nwENjFehvs77loe8e4u1Db/NQ0kP15q+sMfPgikP4e3vy8twEl+8KjHv2QHU108vOsGf4jbz4zUli\nQ3y5fmiES+/jDr45+w2mGpND86yUkme/OkGQ3ovfTWnby8b4yPH08A6mIGA/f9uWzrOzOmbHq+hc\n3CoeVyqmGhNfn/maXyQmw/c/UHHiOL5jxjQad1k8GjvNcwpNXNUnuNHx6IBoBgcPZmPmRheIh1XQ\nGiYKZhmMeGgEkUFO3vDyjkO/ye26d2s5H5DAdYZPodoE2vrrKq4s5v8d/n+M7TWWv037GxrRNj/E\nlJgp3BJ/C+8fe5+JkRNJjki2n3t142lSLpTwj/9JJtTfu13fiyOMhw4BYNq/n9f+8ATninbz0IpD\nfL5oHMMiA11+P1ey8vRK+vfoz4jQEY3O/fdYLj+eNfDn2cOchjM3h1aj5aZ+P+Nflf9mxf6T3D+5\nP2G1/3cUP106PdqqM9iYsZGy6jLGXfcrACqOOm4O5eftiZ+3Z6Nw3aoaC+eLrdnljpgeN50j+Ufa\n3TCnp58Xei+PRqU/MguM9O6hsxerq0f5JSjL7TB/h42y0CQ8MVOT09h0tfzwcsqqy3hs9GNtFg4b\nvx/9e6L8o1i6Y6nd1LUr7RLvbT/D7VfFMGVweLvmd4bpoFU8Kk+fRmsq5707kgj29eI3H+5tsoSN\nqzBW1bB8WzrHzrUuk/+k4STHCo4xb8C8Rruximozz69PYVCEPwtGRzuZoWXM6jcLiRl8D/LedtfV\neFNcuXRL8ViZupK4gDiSBk5CGx3dZFva8ABvLjZIFDxfZK0r5czn4KqEQSEEMcGNCyQ2WYrdllnu\nzoKIDtDEWHdupak76x3PKM5gxckVzOk/hwFBrQubdYReq+eFiS+Qa8zlxT0vUmys5pHPD9MnxLfV\nzt6WYi4qoio9Ha+rRoOUGA8cIMxfx9//J5myihp+8+FetyYTZhuM3PLOLl765iQz39zB7/5zsNly\n/Ta+OP0FXhovZvad2ejcP3acJafQ1KrQXGcMDB7IoOBB9Ox1lH/tzsJQrgom/tTpduKRXpTOwYsH\nuSX+FoQQtRV2nfclcJQoaHNYRzsxG8UFxjEgaIBLQnZjQ/RkODBbOS1LkttBZUkaENErknRLL8yZ\ne+odf33/63h5eHH/yPtddq/EsEQWJixk3Zl13LXq7+SXVvLGgkTnTbHaiemw1ZfzTMwBajTw5coX\n+fD4h2i8L7DstkRSLpTw0KeHMFtcH+uxK/0SN7+1g3NFJt65fRT3Te7HphO5THl1G099eazJ8vI2\n8+z0uOkEetc3rV0sqeCd79KYNiScCf0bRce3iZv73UypPEuV5jzv71AlS37qdDvxWJW6Ck+NJzf3\nuxmwFkl01pYWIMJBoqAtQbCpcujTY6dz8OJB8srz2rXe2BBfsgtNWGofTKUV1RjKq5zneOQdA78I\n8HXNA6GlxIboOWCJx+/iAagNmNubu5ct2Vv47fDf0tPHtetZmLCQSJ8BnKj6J3dNDiEhqodL56/L\n8W2rMQvQj06muE9PQk5f5C/7/sK8dfN49sjPGZX8NVvPfcUfv/7eZfeUUvLBzrP86h8/EuLnzdr7\nJ3Lj8F48dv0gtj02mfmjo/lkTxbXvvIdr208RWlF4+TFDRkbKKsuc+gof2XDKarMFpbe6Lrd2o19\nbsRTeDKg30k+3JXRpRMqFc3TrcSjylzF2vS1TI6eTIiPtWSFT0L9CrsNCQ/UcbG0wv7wBuubv5eH\nhnB/507BaXHWviCbsza3a82xIXqqaiz2plS2XY/TSCt39/BwQri/jsNiILrqQjCcwSItvLL3FSJ8\nI/jVkF+5/H65RdXknJ6DRmMm1fJ3LNI9iXu7L+wma9cmLkb68NoN7xA/6WbizpvZOPMrnpvwHON7\nj6fAnIKu90rWGe7nmn9P55kfnmFDxgYKKwrbdM/KGjOPf3GEp9edYPLAMFbfO54+PS/vNMMDdLww\nZzibFl/D5IFhLNuSxrW1tbgqa8z2cStP15pnw5PqzX80p5gvDuTwvxP6ENfTyQ62DYT4hDAxaiLl\n2h8prazko10ZLptbceXRrcRjS9YWiiqL6r2J2drSVhxz7PcI9/em2iwxGC/bcHMMJqKCfNBonIeD\n9g3sS/8e/dttuooNtkVcWUXD5v+IdiQe5mrIP+WwLEl2aTaPbXuMrJKsdq3HGRqN4EJAbTRP9h6+\nOvMVKYYUHhr1kL1HiqswWyQPf3YIUR3KfQmL2ZO7m09SPnHpPQCO5h/loc2/I/68pP/En6HX6tEn\nJUN1NQHpeczuP5sXr36RLbdu4YuZq4iy3MYlQxBfpa/n0W2Pcu2n1zJ/3Xxe2/cau87twlRjavae\neSUV/Pxvu/l8fw4PTInn3V8l4a9zHAXVN9SPt28fxZf3TWBQhD/PfnWCKa9uY9WBHE4ZUjmUf6iR\no9wamnucYL0X913X32U/Kxuz+s2isKqApEH5/GPnT6+wZGdyMreEqporp7pBtxKPL1K/oLdvb8b1\nHmc/ptHr8Y6Pd+o0jwhsnOuRZTA6jbSqy/TY6RzIO0C+se31uOy5HrUFEm2RVw7NVpdOW8uShNcv\nS2K2mFm6Yyn/zfgvi7cupqLGPdFBoucAyvDFmLWLvx74K8NChnFDnxtcfp/l29LZm1HIs7OHsnDk\n7VwbdS1v7H+D1MJUl90jvSide769hyFF/nhXSQKTxwKgH2XNoDftvxxVJoRgYEg8n972GFHV91KR\n9ideGvse9yXeh6/Wl49TPubuzXdz9Yqr2XV+l9N7Hsgq5KY3d3A6r5TlvxzFw9MGNPmCYmNEdA8+\n+e1VfPS/Ywj00fLwZ4f5n8/fxEN4NnKUf330AnszCnn0+oEEOBElR5jLyjlz000UrV7T5Lhroq4h\n0DuQoPAjFBmr+WSPa1oUuILvc77nTFHXjASrqrEw951dPP/1ic5eip1uIx7ZpdnsubCHOfFzGoWL\nOmpLayPMQaJgduHlPh5NMS12GhLJt1ltLxrXK1CHp0bYdx6ZBUaCfb0cv43ay5LU33l8kvIJBy8e\nZG78XFILU3l+z/NtXk9TRPf046CM58O8H7hovOiS0NyGHMkp4vVNp5mZ0IvZiZEIIXh6/NP4efmx\nZPsSl7RFPVd2joUbF6LVaHnS17pL9UlMBMCjRw+84+Mx7msckhyg0/KP/xmNt6eW/1tr4tb+v+af\nM/7JzgU7WT51ORG+Ebz040vUWBq/jX+2L5sFf9uNt1bDqnvHM2NYr1atWQjBNQNCWXf/RF79+VCM\nXnuoKBrCPR+d5kCW1XxWUW3mxfUnGdwrgPnJrQvNNf64h8rUNC489RTGAwedjvPy8OKGuBs4eGkH\n4/rreff7s1RUm52O7ygOXTzEfd/ex6wvZ3H/t/ezL3dfl2r8tS/TQHmVmYnxbm6x0Aq6jXicNJzE\nV+vL7P6zG53TJQzHUlxMdVZjk05Eg0TBkopqiozVDrPLG9KvRz/6BPZpl+nK00NDVJCPXTyyDUbH\nJiuA3MZlSc4Wn2XZwWVMip7En8b9iYUJC1mTtoZVqavavCZnxATr2SZjeN+zgmlRkxgV7sI6R4WZ\nmP9+Pa//ey2h/t48P/tybamePj15ZvwznC48zVsH32rXbS6ZLrFw40JMZhN/m/Y3vFMy8QwNRRug\ngbevgrPf45OchOngQWRNYxGIDtbz7h3JXCiu4O6P91FZY0av1TMhcgIPJz3M2eKz9X721WYLT689\nzuNfHGF0nyDW3jeRQREBbV6/RiPwDjyG1Bi5feh8zuSXccs7u1j40T6e/zqFc0Umnpo5BI8W7Gjq\nUr5jB0KvR9urFzkPPkB1nvNAkNn9Z1NlqWLk4EwulVXy6V7XVJluK1JK/nrgr4ToQrh3xL0cyT/C\nrzf8mtvX387GjI2YLZ0vbs2x7VQ+Wg/BuH7NtxfoKLqNeEyLncbW+VuJ8G1cTsInwVph15HpKtTf\nGyEu9zJvtv1rHYQQTI+dzr68fRSY2l7vMSbEl0y72aq8aWd56CB7WRKzxcyTO59E56njT+P+hBCC\ne0bcw7he43h+9/OkFKS0eU2OiA3Rs7FnCWYhWBzm4o7CG57AI2c300tX8er8EQTq6++8JkVP4tYB\nt/LB8Q/Ym7u3TbcoqSph0aZF5JvyeWfKOwwIGoDp4EF8Ro5EHPzYWuZ+3YPoExOwGI1UnDzlcJ5R\nMUG8eusI9mYU8oeVl3e0k6MnMypsFG8fepvy6nIM5VXc8Y8f+WBXBr+Z2IcPfz2GIN+mGzC1hJWn\nVxLlF8XS625i22OTeXjaAHalF/Dx7kxmDI1o0wOobMdOfMeMIeqtN7GUG8n53QNYKh2HCQ8JGUK/\nwH4cKd7MmLhglm9L71Rb/Q8XfmBf3j7uSriLexLvYeO8jfxx7B8prizmkW2PMHP1TP5z8j8t8kl1\nFttO5zM6Lhg/7yunKEi3EQ/AqePWu18/hE7nMN9D66EhxPdyoqBNPJy+/TdgWuw0LNLSLtNVXG1f\nj2qzhfNFFU2E6R6vl9/x4YkPOZJ/hCfGPGEPlfXQePDSNS/RQ9eDh7c+TEmV65r41HjkYAg4wy+K\ny4i+5ELbcvoWOPkVBdKfudrdjI9yXH7k0eRHiQmI4YkdT7T6+zLVmPjdt78jvTid1ye9TmJYIjX5\n+VTn5OAzIgEOfAQ9YsBwBr3ZarIy7d/ndL6bRvTm4WkDWHXwHG9/lwZYXyYeTX4UQ4WBV35Yzs1v\n7WB/ViGv3jrCJYl6YE3K3Je3j7kD5qIRGny9PXlgSjzbHpvEkz8bzHOzWx+JV5WVRXVWFr4TJ6Ib\nMIDeL79ExZEj5D79jEPTjxCCm/vfzOH8w9w6zpsLxRWsOpDT7u+tLUgpWXZgGb18e3HrgFsB63Ng\n/mwJs18AACAASURBVMD5rJ29ltcnvU6wTzAv7HmB6V9M5+1Db2OoMHTKWp1xodjEydxSrm1FqfyO\noFuJhzOsbWmHUuHEaR4e4F1n59G4j0dTDAgaQGxAbLtMVzHBekorajh+vgSzRTre9ZTlQ1mePUw3\nvSidtw6+xdSYqY2c1sG6YF699lVyy3NZumOpS2y/Ukr+k/42mPXcUh0C2Xuav6glmKvhmyUUekdx\nf81DeFlMcPQLh0P1Wj0vTnyRfGM+L+x5ocW3qLZU88jWRzh48SAvXv0iEyInAJfrWfmEmaE4G6Y9\nB8PmoT3+Ltpe4Q79HnX53XX9mTMykr9sPM1XR84DMDx0OAk9ruWL9E+okkV8fvc45iZFtXitzbEy\ndSWewrOReTbEz5vfXt23TXW/yndaqwb4ThgPQMC0afS8916KV6+m8ON/ObxmZt+ZaISGC+YdjIgK\n5J2t6dR0Qh+ULdlbOF5wnHtG3IOXR/1dnYfGg6mxU/nXDf/iwxkfkhiWyPLDy5n+xXSe++E5t0Um\ntpZtp6wBN5MGhnXySuqjxKMWn+HDqUhJQVY3TmyqmyiYZTASoPNscRE5m+nqwpE9pN4yh+q8i61e\nm61A4vbT1l8ih+KRVyt84cOosdSwdMdS/LR+PDn2SYcVZhPDEnkk+RG2Zm/ln8f/2eo1NWRr9lb2\n5e1FV3YDFzyHQs4+cIUt+cd34dIp/qL5NSJ2AoQNhf3O1zs8dDh3j7ibr898zTdnv2l2eou08OSO\nJ9l+bjt/HPdHZsRdbkdsOngIodWiK/4OfENh4I1w/fPg6Y2+pwnj/v1NCq8QgpfmDmd0XBCPfHaY\n/ZmF/N9/T7Jr3xiEsDBp7AFGRLsuubHaXM3a9LVcG32tS5Myy3buRBsZiVdcnP1Yz/vvw++668h7\n+WXKdzd+UQjThzGu9zjWnVnHvZP7kWUwsq5WQDsKs8XMWwffIi4gjpv63eR0nBCCUeGjePO6N/ly\n9pfM7DuT1Wmrmbl6Jou/W8zh/MNIKVl3+HyntCPedjqfXoE6BoS7pj+Nq1DiUYtPwnBkZSWVqY3D\nPcMCdPZoq+xCY5OZ5Y6YFjuNn+2poebESUo3/LfVa7OZqbanXqr97CCxq04DqH8e+yfHC46zdOxS\nezKkI24ffDvXx13PXw/8tc1+ArA+tF7d/yp9AvvQz2cK+ywDoKoULrbTp1J2Eba+RHnMdXxSOIjr\nh/WCpDvhwmE47zzi567hd5EQmsBzu59rsjillJIX97zI+rPreXDUg3azhg3ToUPoBg9Ak74BEm8H\nTy/wj4DrnsRHl4XZYKDqbEaT34K3pwd/+1Uy4QE65v/tB97Zms6CkYncPvgXbMha59Lw4i3ZWzBU\nGJg3YJ7L5pTV1Rh/2I3vxIn1XkKERkPv/3sZr7g4zj30EFU5jZuEzuo3i9zyXAJ6ZDAowp+3tqS5\npYSLM9afXU9aURr3j7wfT03LfAV9A/vy9Pin2TB3A78Z/hv25O7hl+t/ybX/upWH1n3CS984r4Pn\nDkzVVew4e5LBfc/xScon7fp/6mqUeNSia8JpHhGgo6C8isoas9M+Hk0RLyKYkFJbXuS771q9NttO\n40BWId6eGsIcmR5yj4F/L05VXuKdw+8wI24G18dd3+S8QgieGf8MMf4xPLbtsTbno3x2+jMySzJ5\nNPlRYoMD2FIeZz3RXtPV5meg2sSaiPsBwbQh4ZAwHzx9YP8HTi/z1Hjy0sSX7DswZ9nnbx96mxWn\nVnDn0Dv5zbDf1Dtnqaqi4tgxfMIFSDOMuuPyydG/RT/E2s3QuPv/s3fe8U1V/R9/n7RJm3RPumhL\nyy67bGTKEmWJkyWCG9wD9ZHH/fPBgYiggAgCyhJEUVFAFLBFdpllzw66S/fO/f1x0tKRtklJoUje\nrxcv05t7b+6JSc493/H5RNQ6DHcHDYsndaFFIyfeHdWG/xvdlic7PI6D2oFZ+2fVeryprD21Fl8H\nX3r49qh9ZxPJO3wYfU5OWciqPDaOjgTM/RylpITYadPQ51YU8OzfuD9Oaid+Pvcz0wY05WxyDhsO\nVetEbVGKSor44uAXtHJvxaCgQWYf76Xz4tlOz7Jx9CZCbcaSWpCIrvEyNuU+TP/V/Rm/cTyv7HiF\nOQfmsPbUWnbG7+RS5iWKSsyXZFEUheTcZPYm7GXtqbV8su8Tnv7zaUb8OIIeK7tB4P/Ym/8RM/fO\nZFvMNrPPX1/Ua+peCDEU+AzpJLhIUZT/VXp+HDAdEEAW8KSiKNU7CtUjan9/bNzcyDtyGLcH7q/w\nnI+L/LFOyiwgNj2PQWbKfmesX4+6GKJCVXTcs5eSrCxsnKo6ulWHvdoGH2d7EjLzaeLpYLxxLPEY\nRd5hvBH5Bs4aZ17v9rpJ53ZQOzCr3yzGbRzHyzteZtHgRSbfpUFFr47e/r054nGWNTmu6D28UcXs\ngS5Taj+JMWL3wcFvodezfH/KnnYBdvi5agEttLlb5j0Gvwd2xt/Hxs6NebXrq7y5802WRy+v4q2y\nPHo5Cw4v4O5mdxu1Zi2IjkYpLERLNDTpAx7lXANVNmjGz8Fm1VjyfluO29ja5Veaejuy8dneZX+7\n2LnweLvH+Xjfx/wT/0+FxtW6EJMVw67Lu3iqw1Nl5lmWIDsiAmxscOje3ejzdk2a4P/Jx8Q8/gSX\n33gDv08+KXsv7W3tGdJkCL+e+5VX73md9o1deefnaHqFeta738f6M+uJzY5l3u3z6txrlF9Uwkur\nT3DweDteGnI3ebaH+Wr3PzQPs6VYlcrh5MOy1Fe5Gp4VCLx13vg7+uPn6Ie/o3/ZP18HXzIKM7iQ\neYGLmRe5mHGx7HFu8dWJV6PSEOgcSIhLCLqidkSdVfP1uGG09grFzc7tmt8bS1Fvk4cQwgaYBwwC\nYoG9QogNiqKUb5E8D/RVFCVdCHEHsBDoVl/XVBNCCOzbtTWaNC/9oB+Jy6CwWG9Sd3kpil5P+qpV\niA5h/N7hIh3PZnPi91WE3fuoWdcX6KEjITPfeL6juBCST7CokR8n0o4zu99s3OxN/5A1c2vGjO4z\neD3ideYcmMMLnV8w+diFhxeSWZDJS51fKpOQB0G2Vyec67ry0Oth48vg6ENih2c4uHU3Lw9pcfX5\n8Elw8Ds4uk4+robRTUezPWY7nx34jO6+3WnhLs+x4ewGPtz7IQMDBzKj+wyjOaFcg3+HVhcP4e9U\neV407oyumR+50efh8mHwbWf2MB9s+SArT6xk1v5ZrPZdfU0NletPr0clVIxuOrrO5zBGTuROtO3a\nYeNcfe+JY58+eD3/PMmzZmHXqhWej179bI8MHcnaU2v5M+YPZt03kGGf/c30dYdZPKlLvXnA5xfn\ns+DQAjp6d6S3f+/aDzBCTkExjy7bx86zqbw7MowJPYLJLQxl0SZHfAob87bBLbFYX0xSbhJx2XFl\n/+Kz44nLjmNf4j5+PfcrCkYq0hD4OfoR7BxMR++OBDkHEewcTJBLED46n7IbgGGf/U17N1tua9yl\n7m9IPVGfK4+uwBlFUc4BCCFWASOBsslDUZTyWg27AMuVndQBbZu2pOz4u4otbWmj4N4LsoTPlB6P\nUnIiIymKicHvuY95o2somavHcHDNHBK7hTIgcIDJ5wly17HnfJrxfEvKKY7bwsKsk9wZcie3B91u\n8nlLGR46nINJB1lybAntvdtze2Dt57iUeYkVJ1Zwd7O7y36YS/Mz8U7tcL7wu6wCczSzxPDgdxB/\nAEYvZPOZbACGhJVb7QV0Ae/WMnRVw+QhhODNnm9y909381rEa6y8cyU743by38j/0s23GzP7zKx2\nlZV38CBqVw1qDzdoWdULA0A75H6yPplD0apnUT+7FVTm/fhrbDQ80/EZpv89nV/O/VKm9GwuRfoi\n1p9ZT2//3kb7mOpKcXo6+UeO4Dltaq37ejz6CAUnjpM861PsW7TAsU8fANp7tSfIOYifzvzEqKGj\neO2Olrz1czSr9sbwYFfL+syXsvrkapLykpjZZ2adJqiM3CImfbOHQzFX+OTe9owJD6AkO4eEBx9k\nbr6KrbEdKertj9rdDVuVLX6Ofvg5+tGFqj/wRSVFJOQmEJcdx+XsyzjbORPsHExjp8ZVqr8qk5SZ\nT/TlTF4Z2qLG/W4U9Znz8AfKt5bGGrZVxxSg9vKYekTbri0oCvnRxypsb+Rsjx2FiFObEOir9fEw\nRvqKldh4eOA8aBDNvVriPmAgHc8qPL/1WZYeW2pymWzpj7KxBsGiy4f5j5cHbnYuvNb1NZOvrTLT\nu04nzCOMGREziMmsvSv40/2folapK3h1lAo5nlIbpL5j95h3EXlX4I+3oHE3aHcfm44lEuLpUMHD\nHSHkpBEfBfEHazydu7077/R6h9Ppp3lp20u8tP0lWnu0Zk7/OdV+eRVFIe/AfrQumdBhLNgaL2/V\ndZeNkLmHoyFqmXnjNDC0yVDCPML4POrzOmuO7YjdQUpeilHp9Wshd9cuUBQce/WqdV8hBL7vvYdd\nixbEvfgShRculG0fHjKcfYn7iM2KZWKPYHo19eDdX6KrmJxZguzCbBYdWUQvv14VbIpNJSW7gAe+\n2sXRuAy+GNeprIw6bckSCk6fxqs4h4d2reJ0797EPPEkGb/+WiXXUx61jZrGTo3p7tud0c1Gc3vg\n7YS6htY6cYCssgLo29wLRVEoOHPGaGHCjaJBJMyFEP2Rk8f0ap5/TAixTwixL7ka3w1LYN9WCgpW\ntqV106mZqfma/2a9zX0226v3Dq9EYWwc2du24XrvPQiN/LB4DhyKLreEsSXhfLzvY97b9Z5RraPK\nBBoqrIxVWn15bh2nNRre7PFWFdMfc9DYaPikn4xZ1yaguC9hH39c+qOKV4eLTo2LVs2B4kBQqc1P\nmm+fCbmpcMeHZOQVs+tcKoPDfKreQba7D2zt4cDSWk/ZJ6AP97e4n22x22js1Jgvbv8Cnbr61WNx\nfDzFySloPfOh06Rq97Nv1Qqh05KXHwxb3pQWwGaiEipe7PwiCTkJfHvceM8EAEX5sHuhXMlVYt2p\ndXhrvekdULcQTXVkR0SgcnYu+17UhkqnI2DuXISNDTFTp1GSLVeNI0JHIBD8fO5nVCrBR/e0x0Yl\nePF7yxtoLY9ezpWCKzzd6Wmzj72ckcd9C/7hfEo2ix7qUqYvVpyaStqSJTgNGULwr7/yzIAXONNn\nOPknThD/4kucuq03cS+/QvaOHUZL/euCUlLCse17GB+7E+f/e4PTPXtx7q7hpK9YYZHzW4L6nDzi\ngPLqawGGbRUQQrQDFgEjFUUxquGhKMpCRVE6K4rS2cur/rosbd3dUQcEVPH2EIdWMUr1N9mKPS+q\n12GnmCa+d2XNGhACt/vuK9vmcNttYGvLpPQwJreZzJpTa5i2dRrZhdk1nqtvcy8m9QymW4h7he1H\nU46yOPsUI4s19A0yPQxWHf6O/nzQ+wNOpp+sttFOr+j5aJ/06pjYemKV54M8dJxNLwG/DhBjxsoj\n6QTsXiBXFX4d+PNkIsV6pWLIqhStG4SNhsPfQ0HN7x3Ai51f5Pnw51k4eCGu9jX3VuQeOACArm0b\n8KxetlzY2qLr0IHcDDcozJYTSB3o4tOFfgH9+PrI18a7m/Mz4bt74LeX4aenysy2AC5nXyYiLoJR\nzUaZVehQG4qikBMRiUOPHggb0xPwmgB//GfPpvDCBeJfmY6i1+Pr6EtXn65sOLMBRVHwc9Xy9ogw\n9l5It6jfeXp+OkujlzIwcCBhHlVtCWriYmoO93z5D0mZBSyb3K1CN3fKl/PRFxTg9dyzuOo0+HZu\nzychgwnd+geBy5bictddZO/YQcxjj3O6T18S3nmH3ANRZjXfKoWF5EZFkbLwKy499hinunXn/q/e\nYNy+Hyg4eRLHfv3wff993MeNNWtc9Ul9Th57gWZCiCZCCA3wALCh/A5CiEDgB2CCoiin6vFaTEbb\nri35h8vJlKSchl9f5Ji6LY8WvYg3abJxrRb0hYVcWbsWx/79Ufv5lW23cXJC16Uz2du383z487zV\n4y12X97NhN8mEJ9dfROVi1bNWyPCKlitFpQU8EbEG3iUKLziHl7tsebSJ6APj7Z9lPVn1rP+9Poq\nz/967leiU6N5ttOzRiVfGrvrpGlV424Qd0Am9GtDUeC3V2T11IAZAGw6moi3kx3tq3MJDJ8k+0mO\n1S7yqLXVMrnNZLx1tXfp5u3YiLDVYzek9koxbXg4BecvUdL+MVkddvGfWo8xxvPhz5NXnMeCQwsq\nPpGTAkuHw6V/oNVwOL25Qof9+jPy/8/dze6u0+tWR+HZsxQnJuJwW+0hq8o4dO9Go+nTyf7zT1Lm\nfQHAiKYjiM2OZU+CvJkY3dGfO9r4MGvzKY5ftoxEzuKji8ktyjXb8vhUYhb3zv+H3MJiVjzaja5N\nrt6gFcbEkL56Na5jxmDXRJZn39HGh4upuZxMysGha1d833mbZn/vIOCLeTj06M6VdT9wcexYzg4c\nRNKsT432junz8sjZtZvkufO4OOlhTnbtxsUHx5I8axZFcfEU9LmdmeFjufDlKppu2YzfB/+H65i7\nUfvXFPm/vtTb5KEoSjEwDdgEHAfWKIpyTAjxhBDiCcNu/wU8gC+EEAeFENWLBV0n7Nu2oyg+nuKU\nFBkq+P5hUNvzXcAM/tGHcdyxG/z9CeTV7BKXtWkzJWlpuD34YJXnnPr1o/DMWQpjYhjTfAxfDvqS\nxJxExv46liPJpjchfXHwC85mnOXt5GScfTqYPdaamNphKt18uvH+7vc5kXaibHtecR6zD8wmzCOM\nYU2GGT02yF1HXHoeJf5doKQAEqr3iC/j+M9wfjsMeAMcPMgvKmH7qWQGhzWq3tOicTcpBFlDz0dd\nyNu/B62ngmhTVYG5MrrwzqAo5Or6gktj+PUFKaliJiGuIYxpNoY1J2XPDABXYmDxUCnI+MAKuHcp\n+HeG36dDTiol+hJ+OP0DPf164u9o2R+V7AjZv2JKvsMYbhPG4zJ6NCnz5pG5ZQsDAwfirfPmlR2v\ncCr9FEII3hvVBmetmudXH6zggFgXknKTWHliJcNDhxPqGlr7AQaOxGZw/wI54a9+vEcVO+Pkz+Yg\nbGzwnHq1aGBQ60YIAb8dvdqAqtJocBowAP9Zs2gWGYnfzP+hCQkh9euvOTd8BOdGjiJl/gKSPpnF\nhQfHcrJrNy5NmkTKvHmUZGbget+9+M/5jGaREYT++gt/DH2YHYGd6Nmt9TW9L/VJveY8FEXZqChK\nc0VRQhVFed+wbb6iKPMNjx9RFMVNUZQOhn/mZ7gszFVb2iOwZYaU/Rj1JfbuMgJ3oOmzMowQMbvG\n86SvXIk6KBCHnlXr9x379wcg29Aw2N23O98O+xZ7W3smb5pskg7WoeRDfHPsG8Y06s5tefkWt561\nUdkws89MXDQuFQQUlx5bWqtXR5CHjmK9QoKToXy1trxHUR5s+o+UHgl/GJDd9HlFJQwJq6F6qDRx\nHrdflstaAH3yRfIv56Bt2xrUtfciaNu3A7WavENH4Y6ZkBQNu76s02s/2UHqL3124DPpCLl4iOyy\nn/AjNB8CKhsY8bn8/P3+KpHxkSTmJjKmuWUT5SBLdDUhIRVWzeYghMDnrTexb9eOy9NfxeZCHF8P\n/hpbYcuUTVM4mXYSD0c7/nd3W04kZPHZH9fWab/w8EJKlBKebP+kycfsOZ/Gg1/tQqex5fsnetC8\nUcWeofzoaDJ/+QX3iRNRN7q6YvV0tKNLsDubjhpXL7BxdMBl5EgCv1pIs+3baPTGG6i0WpJnzyZ1\nyRLQ6/F4aCKNF8yn+e5dhPzwAz6vv47z4MHYekhFiG2nkukU6FZFPboh0SAS5g0J+1atQKUi/88f\nZHiqxzRoPqSsUVAb2F4ma3fPh0zjYab8kyfJO3AAtwceRBgp39QEBqIJDa3QbR7iGsJ3w76jhXsL\nXtj2AouPLq42ZppfnM8bEW/QSNeIlxwMZXyNLO9b7qH14ON+H3M5+zIzImaQmJPI4qOLGRQ0qIov\ndnkCDRVX5wqdpBJtzG4URSH77wiUEiN3mJGfQcYlGPYh2Miw3KZjCTjZ29I9pBb58Hb3g42dSYlz\nU8jbMA8UgXbgvbXvDKi0WrStW5O7/wC0vBOa3wHb/gcZ5qvIemo9ebjNw2y5uIWo5cPkCubhjRBU\n7gakUWvo/QIcWcPaqC/xsPegX+N+Zr9WTegLCsjduxeHOq46SlHZ2RHw+RyEg46YqdNojBtLhi7B\nzsaOKZunEJ0azcDWjbi/c2Pmbz/L/ot1U7ONyYph3al1jGk2hgAn06r9t59KZuLi3Xg727H2yR5G\nC1GSZn2KysUFj0eqhi/vaOPDycQsziXXnG+z9fTEffw4gletpOmO7bTYs5vg1avwfuklHPv2Ndo/\nk5JdwOHYjAanolsZ6+RRCZVOh11IMHmRm8C3A9wuk6CNDL0ege466P+6FP3b9j+j50hfuRJhZ4fr\n6OrDHk79+5G7d19ZRQrIH+tFgxcxNHgon+7/lLf/eZsifdUQyOdRn3Mh8wJv93wbx+RT4OQHOvcq\n+1mCjt4deT78ef6M+ZOHfn+IIn0Rz3d6vsZjrlrnGvIeMXvI2vgbMY8+SvqqVRV3vnIJIj6FsLsh\nWJa+Fpfo2Xo8kdtbeqOuTaZc5w5ho+DwGijMqfM4AVAU8nbIanFtH+MhOWNoO4eTd/Qo+vx8ufpQ\n9PD7q3W6hIm6JniV6PnERYcy+XfjK8reL5Lk1ZwdqUcY2eQO1CrL3p3m7d+Pkp+PYx3yHZVRN2pE\nwGdzKLp8mbgXXiTArhFLhi5BZ6vjkc2PcCzlGG/c1Qo/Vy0vrDlUJ8/zLw9+iY3KhsfaPWbS/r8f\nvcwjS/fSxNORNY/3wNelavVkzq7d5ERE4PnYY0Z/4EtXxL8fq147rTJqb29Uutp7xP4+3TBVdCtj\nnTwqU1KMVpdIfooK5Z7FUgwPGed8c3hrOga6gVuwlN2IWg7JFfP8JdnZZGz4Gec778TGtfqqHsf+\n/aG4mJyIitpI9rb2zOwzk0fbPsq60+t48o8nK3hTHEg8wPLo5dzf4n4paZFwtIrtrKWZ0HoCg4IG\nEZcdx7iW42jsXLOFqY+zPRoblazjb9wNsi6TtmQRAGmLl1R04Nv0HxAqGPxu2aa9F9JJzy2qOWRV\nnvBJUJAJx6om983iwt/kxWSj8fPE1s30Dn1deGcoKiLv8GFwC4K+L8sczqnN5r1+9AZ0q8YzrVjL\nIVvYknXW+H62dvzY7g5KhGBMouU9wrMjIhFqNboululq1nXqiM9/Z5ATGcnJDh3Jv3M8c7935skf\nC/n1tfFcWDmP2c2KKIyJ4YOfzRMePJN+hl/O/cLYlmNNKob44UAsU1dE0dbfhVWPdsfTsWoPj6Io\nJH3yCba+vriNH2f0PH6uWtoHuFQburoWtp1MxtNRQ5hf3R0lrwfWyaMy2z7A3i6OkkJBUc7VOzqd\nxpaHezW5at/Z+yVQ6+DPitIVGT/9hJKbazRRXh5t+/bYuLiU5T3KoxIqnun0DO/2epf9ifuZsHEC\nsVmx5BblMiNyBn6OfrwQ/oKsYko5afF8R2WEELzb611e7foqT3aoPaasUgkC3A3WuY27kpeiJu/o\ncRz69KYoLo7M3zfJHc9tg+MbZBjG5Wq4YdOxBOxsVfRtYeKyPbAHeDa/5sS5sm8Jeal2aLtWFQGs\nCV2njoC8Ywegx9Pyeja+BIUmNsLtXwrfPwR+HRk57neaujZl9oHZRoX29IqeH5L20E3tQeD+FXDJ\nQt4pBnIiI9GGh5t0l2wqbvfeS+MF8/GcNhWH227Dzs6BrslODIvMx3bmAhxef47Fmz/ggTfGcfS2\nvlwcP4H46a+S/Plcrvywnpzdeyi6fLnKeecenItOrWNym8k1vr6iKCzdeYEX1hyiWxN3lk/pVm0+\nIWvTZvKPHMFr2jRUdtX7nwxp48Oh2AzirljOgbBEr7DjVDJ9mnlVXyjSQLBOHuU5tw3+/gRtz8GA\ncYXdMhy9oOcz8g4zRsokK4pC+sqV2Ldti7ZtzT/owtYWh759yN6+w3geAOkFvXDQQlLyUhi3cRyv\n/f0al7Iu8W6vd2WTW8pJ0BfXS76jMg5qB8a1GoeD2ogcvBGC3HVcTMsF7zDSTrugsrfFf9YsNKGh\npC5ahFJcCL9Nl6u4HlcbuhRFYUt0Ir2beVUoS66R0sR57F65EqsLOSkU7tlISYFAG25e3YaNqyt2\nzZpdNYey1cCdn8CVixBhgmpuxGz4+RkIHQAT1mPj4MkL4S8QkxXDmlNrquy+K34XcdlxjOn8rJx0\nNzwNxZbxmShKSpJ9BRYIWVXGsW9fvKZOxe//3ido2VJa/PkXHv/8wfsvN+Z/E3TkvvwoWzoN4x/n\nYIpL9OTs2UPKF19w+fXXufTQQ5zpP4DEjz4qO9/RlKNsvbSVh8IeqrF353xKDhMX7+HNDccY2Mqb\nxZO64FCNnatSXEzy7NlomobiMmpkjeMZalgZW3L1cSQug/TcItNvnG4g1smjlOwk+OEx8GyO3UNz\nEPb25Buxpa1Aj6nSJOiPt2S55t69FJ45W+uqoxSn/v0puXKFvIPVS2x08enCt8O+xUHtwJ8xfzK2\n5Vi6+BjCCeU8PBoaQR4OXErNoTAxicxLGlzDNNg4OuIxZQoFJ06Qs3iGLEEd8kGFqqajcZnEXclj\nsLHGwJpo/yDYaOqeOD+0krxEeaen69jR7MO1ncPJi4q6GpJr0kcm8yNmy14hYygKbJ4Bf7wJbe6B\nB1aCRk7Ot/nfRjffbsw/NL+Kpe7a02txtXPl9tA74a7Z8iZix8dmX7MxcnZKublrTZabiq+LP7Me\n/Ja01n48ab+KgNeHM7PD/Xw6/EWa/fUnLQ8dJHTzJgKXLMZl1CjSvl5M+mo5oX4e9Tludm5GG1VB\nquLO2nySIZ/u4OClK7w1vDXzx4djr66+6fHKuh8ovHAB7+efr7U5MsTLkRaNnMzKe9TGtpNJWxi2\neAAAIABJREFUCAG9m1knj5sDvR7WPwH5GXDvEoTOBfuwsJpXHgB2jtDnFbgYAWf+IH3lSlQuLjgP\nu6Pm4wyUdptnb9tW435NXJrw3bDveK3razwX/tzVJxKOyEojj+q7oG8Uge46cgpLiF/6LSBw978I\nhTm43HUnto28SF2xHkJvhxYV36vN0QmoBAw0U/YenTu0HgmHVpseKipFUWD/N+TlB6BydkYTEmLe\n8ci8hz43l/wTJ69uHPyeDG3++kKFrnAASophwzTYOQe6PAp3f1WWX4OrfucZBRl8feTrsu0peSn8\ndekvRoSOkPpIzQYaJqlZV28mroGciEhsPD2xa3H9xPi8dd4sHrIYXwdfPjr0Mvfels8vhy+z4VA8\nQqNBExiIQ48e+L73Lg69e5Pw7rsc+PUbdsbvZErbKUZXw3+dSGLQp9uZ8+cZhrX1YeuLfZnUq0mN\nPvH6vDxS5s5F27EjjgNMU2sY0saHvRfSSM6yzMpv+6lk2ge44u5Qu/bVjcY6eQD88zmc3QpDPyhL\nPmvbtCE/Orp2rZrwSeAWTNFPb5K15Q9cR49GZW+aV4GNkxO6zp1NMohys3djbKuxaG3LVYYkHgXv\nVmXlrQ2JIA8ddsUF5P6wFqce7VHriiDuAEKjwb2LK7kJNuQFTpYhp3JsOpZA1ybudfvyhD8MBRkQ\n/aN5x12MhNQz5KXr0LZvb7S8ujZ0nWXpct7+cn2ujt5w+ww4v0PKx5dSlC/zG1HfQt9XYdhHRhV5\nW7q3ZHjocL6N/pbL2TLev+HsBoqV4ooiiEM+AHsXGb66ButfRa8nJzISh5496vQeXAteOi++HvI1\n/o7+/JH2AS2CLzPjx6MkZFzVVxO2tvjP+kTa4b7xMW1y3Lm/RUXvnbgreTy+fB8Pf7MXjY2KFY92\nY/YDHU3yD0lbtpzi5GS8X3rRZDXeO9r4oCiwJTrRrPEaIz2nkIMxVxp8iW4p1skjZi9sfUfetRoa\n1ADsa7ClrYCtBgbM4Mqui1BcXMVIqjac+l/tNjebxGMNMmQFcvK4PWY/IisL90celxtjdkPcftzs\n/kalU5O6dlOFY86n5HAqMZvBresoKx7UEzyamZ843/8NJcKFgtgUtB3r1qmv9vFB7e9/Ne9RSufJ\n4NcRNr0uV7alOlUnfoGhM6H/a1Um0PI83fFphBB8HvU5iqLww+kf6OTdiRDXcqsjBw+440PZLLl7\nfp2uHyD/+HFK0tNxvO22Op/jWvDUerJo8CIaOzcmxWE+RZrjvLLucIV+JxsnJy6/+TCFooTp3xdj\nmyWT1YXFer7cdpaBn2xnx6kUpg9tyW/P9qFnqGle7iVXrpC6aBGO/fqhCzdd6qeljxNBHjqLhK52\nnE5GUaDfTZDvgFt98si7Ausmyz6J4XMqfIm1NdjSVkZpMYIrF1xxCBBoAnzNuoTK3eYmk5UIOcn1\nXmlVV/xd7Bl5NoKMwKZoe/QFzxZwaRdsfAWVmxdu48aT9cdWCs6dLztmk+ELaHa+o5TSxHnMbkiM\nrnV3AHLTIHoDeU79QVHqlO8oRdc5nNz9+ys2d6ps4M5ZMqf223SpU3VxpwxTdX+i+pMZ8HHwYXyr\n8fx87meWRS/jYuZF4x7lbcZAsyHw53uQfqFO158TEQmAQ0/zqs0siYfWg68Hf00Tl2A0/kuJjPub\nb3dfKnter+j5LGElS8f7YJeSRdwzz/LPiQSGzfmbmb+foHczT7a80Icn+4WisTX95y1l4Vfos7Px\neqHmHqbKCCEYGubDzjMpZORdm6Lu9lPJuOnUVSRSGiq37uShKPDzs7JL/J7FoK34P0wdEICNqyt5\nR2ufPLK2b6c4W49bcCrsW2LWZZR2m9eW96hCouG66rnHo66U7NlFYHYS+7sMkSGAxl3hzBaI2wcD\n38Z90hSERkPq4qvx/M3HEmjj70yAmR7xFTA3cX5oFZQUkFcQBCoV9m3NdwQsRRseTklaGoXnL1R8\nwr8TdHkEDq2URQIPrpQqBSYype0U3Ozc+HjfxzhpnIx7cgsBd80CYSM/12YoupaSExmJXatW2Hqa\ndrdeX7jZu/H1kK9p7tYUh8bf8sFf6zifIhtAN1/czIm0E9w56kUcZ7xF7t697H5mOgVFxSye1JmF\nEzub/fkpio8n/dtvcRk5Evvmzc2+3qFtfCjWK2w9XvfQld5Qotu7mdfVdoAGzq07eez/RsbGB8wA\nIxaPNdnSVubKypXY+vrg2KsL7PhQhibMwKl/P3IqdZvXSgOutAJIW7qMTAcXdpTaszY2uAsHdIV2\n92Pr4YHrmDFk/LSBosREkjLzOXDpCkPqGrIqxcEDWo2QP9RFtdTfGxLlBHQh73Qcds2bV3CQNBdd\nZ1nim7vfiL7ngDegwziY+JPUqTIDJ40Tj7eXob+7Qu4yqmQMyLLdgW/KkvOD5vk+6HNyyI2KwrHX\njVt1lMfFzoWvBn9FM7dm2Pgu4/F135BfVMi8qHk0dW1K8uVWDDvmwOqWgxhyaQ/fu51jQMu6rViT\n584DRcHrafPUeEtpH+CKj7M9v19Dye6x+ExSsgtvmpAV3KqTR+IxKR8ROkD2alSDtm07Cs6cQZ9T\nvexFwfnz5Oz8B7f770cMflsaGf0z16zLcezXD4qKqnSb10jCUXD2rzdZkmuh4MwZciIiON19KOcy\nDHLsTQeCf7jsfzAkY90nPwx6PWlLl7HZkHAcbGpXeU2ET5L5heifat7v0i5IOYnSYQJ5hw7VOd9R\niqZJE2zc3cmrnPcAubId9QUEdq/Tue9rcR9TO0zlkbaP1Lxj5ymyaXLT6zJUZiI5e/ZAUZGsAGwg\nuNi5sOSORQQ4NCPebiFjf3yeC5kXuBI3gLd/PkGHQFcmLnwfpzuGkv7pp2T98YfZr1Fw5gwZP/6I\n29ixdZY7V6kEQ8Iasf1UMrmF5surAGw/Jf9f9blJkuVwK04ehTlSZt3OGUYvqNF32r5tG9DryY+u\nPn5+ZdVqUKtxHTNG/ji2Hgk755r1xdV26FBtt3m1JNa/LEldSVv+LcLOjpwhw0nOKiCvsAScfeHR\nP8H3alhIExCA89ChXFm1iu37zxLsoaN5I8cazmwiwbeBe2jtifP934DGiQL79uhzcq4p3wFytaoL\n70TufiOTxzWiVql5ov0TtUtwqFQyf1eUCxtfNvn8ORGRCK0WbadO13illsVZ48yaEYtxUYVwOncH\nJXkB5F5pxbyxnVg2uSsh3k74ffAB9m3bEvfyKzV+V42R9OlsVDodHk88fk3XObSNLwXFeradrJvT\n6baTybT1dzEql9JQufUmj99fhZRTcPdCWUpZA1qD/WZ1SXN9Xh5X1q/HedAgbEsdDgf8F4rzYcdH\nRo8xhind5hUoLpBjaIAhq+L0dDJ++gmXESPwDZLFA5fSqu+78Hj0EfS5uXhu/YUhxuxm60Jp4vzS\nP9Kd0Bi5aVILq9195B2VvRnaDtfuiaIND6coNpaixGsv3awzXs2h7ysyLHviV5MOyYmMRNe1CypN\nw+svcLJzYs2IxbgW9WeY79NsfbEfd7bzLfusqOztaTxvLjaursQ8+RRFiabduOUeOED21q14PDLF\nLC0zY3QJdsPdQVOn0FVGbhEHLqXfVCEruNUmjyNr4cAyqaUU2r/W3W09PFD7+0tvDyNkbtyIPjMT\nt7HlOso9m0KniTJxnma6xWZZt/mhQ7XvnGyQJWmAlVZXvl+Lkp+P+8QJBLnLxOXF1OrDfvYtW5LT\nvgvDz/7N4KYWrDLpMFZ6qFeXOD+8RhpVhU8iLyoKGw8P1I1rFnw0BZ1B2iR33w32Nev1nLy5+PVF\nGcKrgcLYOAovXMCxV8MJWVXG39WNvx+Zw8cjh+FoRFrE1suLxl9+QUlWFrFPPYU+r+Z8lxQ/nIWN\nlyfuE413qJuDrY2KQa0a8eeJJLONrSLOpKBXuGn6O0qp18lDCDFUCHFSCHFGCFFFo1oI0VII8Y8Q\nokAI8VJ9Xgtp5+Dn52Titt9rJh9mX9mW1oCiKKR/twK7Zs3QVq4L7zsdVLbw5/smv05Zt7kpoatE\ng35TA1t5KEVFpH/3HQ49e2LXrJmUr6fmlQfApnaDcCvIpsm+bZa7GAdPadt6cEXVxHlpotyvE/i2\nI/dgFNqOHSyy6rFv1RKVTndVJPFGYaOWxlHZibDlvzXumhNpKNGtBz0rQHbT714AHzeXlWAFWfXy\nMvYtW+L/8cfkR0cTP/1VFL2+2n2zt20jb/9+vKZOtZgA5NC2PmQXFBN5JsWs47adTMLZ3pYOjW+O\nEt1S6m3yEELYAPOAO4DWwINCiMqeimnAM4BlhHlqIvWsTFqOWSS/WCaiLbWlTU2tsD3/yBHyo6Nx\nG/tg1R8dZ1/o/iQcXQuXTVhJYF63OYnHwNZexvUbEJmbNlOcmIj7Q/JOzlWnxsneVqrrVkN+UQnf\n5XuSEtCUtCWV5NqvlfBJkH8FojdU3B67F5KPQ/gkilNTKbp4CZ0FQlYgQ5DaDh2qNgveCPw7Qfen\n5ER5ofpijJyICGz9fNEYPLotytm/YP5t0p/ewVuqB3/Zs8bruRacBvTH+5VXyNq8meTP5hjdRykp\nIXnWLDRBQTJXaSF6hnrgZGdrVuhKURS2n0qmd3OvGqVTGiL1ebVdgTOKopxTFKUQWAVUkKlUFCVJ\nUZS9wLV115hCs0Hw9AHpbGcGFWxpy5G+YiUqnQ7n4SOMH9jrWbB3hT/eNvm1TO42TzjS4GRJFEUh\nbelSNMHBOPTuDcgEcpCHQV23GiLPpJBbpMd2wiSKYmLI2mymB0ZNBPcG95CqifP934DGEdqMKQsT\naq8xWV4ebedwCk6fpiSj5nDRdaH/f6Ry8YanjZYuK8XF5OzahWOvXpbJN5WSdh5WjYPlo2Ty/v7v\n4Im/YfLvshflm7vg99drL6euA+6THsL1vvtIXbCAKz9WlarJ2PAzBafP4PXcswi15Yy07GxtGNDK\nmy3RiRSXVL/qKc/xy1kkZRXcdCErqN/Jwx8o/ysYa9h247A1Pxlo37q1tKUtlzQvTk8nc+NGXEaN\nrL4vQOsKfV6Smlnntpv0Wle7zbdVv5OiNMhKq7yDB8k/cgS3iRMq6CIFuUt13erYdCwBJztbOj04\nEk2TJqQsWlSt/a7ZqFTQ6SG4tFPmiUCqChz9AdreA3aO5EVFgVqNfZjl3k9deGepsnzggMXOWWc0\nOhj+mQzbGnG+zDt8BH1WFg6WyncUZMsbpnld5apjwAyYugda3SULGQK7w5OR0kxt1zyY3xtiLbtK\nE0LgM+MNdN27c3nGfyvkn/QFBSR/Pgf7sDCchpjXb2MKQ8N8SM8tYs950yx1txlKdPtZJ4/6QQjx\nmBBinxBiX3Jy3Urh6opKp8OuadMKK4+MH9ajFBbi+sADNR/c5VFwDiiTbK+Nq93mNYSushNlL0mj\ntiaO4PqQtmwZKmdnXEdW9EAI9NARm55Hib7q+ItL9PxxPIn+Lb2x09jiMWUyBdHHyYncabkL6zBO\nJs73GxLnR76H4jwZ0gJyo6Kwb9XKZDFLU9C2bwdq9Y3Pe5QS0g86joedn8Pfs2QI10BORASoVDj0\nqFv/SRl6vezW/zxcKvyGjYan98kbKHWl91bjIPt9JqyXq5KvB0lZleLCa7uGcgi1moDPZqPx9yd2\n2tMUXpISJ+krV1Icf1mKH9aD+GPfFl7Yq1Uma11tP5lMa19nk4QbGxr1OXnEAeXLVwIM28xGUZSF\niqJ0VhSls5fX9Z+hS5PmiqKg6PWkr1qFrnPn2qUM1PZS+C7+QO0NawYc+/Wtudu81OyoAa08iuLj\nydq8Bdd770HlUHElFuSuo1ivEG/EbW3fxXTScgrL7GadR4zA1tub1EWLLHdxjl7Q8k44tEKq2e5b\nAr7twa8jSmEh+UeOorvG5sDKqOzt0YaFNYy8RymD34OALrD1bfi8E8zrBn+8Tc5fm9G2bYuNi0vd\nzx27HxYPhvWPy3zflC2yFN7Zr+bjQgfAkzulpPyOj+CrAXU38zKCjYsLjed/CYoiS3jj40n9cj4O\nPXvi0KOHxV6nPDqNLX2be7HpWAJ6IzdM5cnKL2L/xfSbwvjJGPU5eewFmgkhmgghNMADwIZajmmQ\naNu2oyQjg6KYGHIiIymKialYnlsT7R8Er5bw57tgxFK0Mk79+9fcbV5aaWVmmW5xWhpJs2dTlGR6\n86KppK+QUhju46r6PQd6VF9xtflYIppydrMqjQb3hx4id9euasuj60T4JMhLl6ZLScfKVh35J0+i\nFBRYNN9Riq5zOHnHjqHPz6995+uB1g2mbILnjkgFXkdvSv6aQ96J0zjYHIRfXoAzW827+89KgB+f\ngkUDIP0ijPwCHvlT6piZfF2uMPpLeGAFZCfAwn5ydVRimcIJTXAw/p/PofDSJc6NGk1JRgZeL75g\nkXNXx9A2PiRmFhAVc6XG/SLPpFKsV27KkBXU4+ShKEoxMA3YBBwH1iiKckwI8YQQ4gkAIYSPECIW\neAF4QwgRK4RocK7vZUnzw0dIX7ESG09PnAYONO1glQ3c/iaknoGo5bW/Vm3d5olHZShMa3pTk76w\nkNip00idv4CL48bXTf69unPn5pK+5nucBg1C7Vf1TjOwrNej4uShKAqbjiVwW1PPCnX7rvffh8rJ\nidSvLLj6aNIX3JpIuXK1g3TtA5nvwDLNgZXRhodDURF5h2pxo7zeuAZCt8fhoZ/J6bYQFIFDx1ZS\nC+zbu+GjUFg7WfqPVKfRVlwgHRI/D5f9Mr2ehaf3Q8dxNSo21EjLO+Gp3dBymFwdLRkKKWfqPs5y\nOHTtiu9bb6HPzMR52B1o65rfMnFCG9CyEWobUaYSXR3bTyXJfF/QtTUo3ijqNeehKMpGRVGaK4oS\nqijK+4Zt8xVFmW94nKAoSoCiKM6KorgaHpunKngdsGvaFGFvT9amTWRv347rPWMQ5nTitrgDGneH\nbTOvJm6rodZu8wTzkuWKopAwYwZ5UVF4TpuGPjOTC2PHkn/ylOnXXwMZP/2EPjOz2kYrXxctahvB\nxbSKSfNj8dJudkgl+XUbR0fcHnyQrC1bKDh/HougUkH4Q/Jx2zFgL+9PcqOisPX1Re1jAT2tSug6\ndQIhjIskNhBy9h5E5eSE9vnv4ZVz8OBqKa9zbrucQD4MgeV3w96vIfOyzNud/A2+6C5XccG3wdTd\nMOidsvf0mnDwgHuXwpivpXXv/Ntg13yZT7lGXMfcTdDyZfi8827dTpB6Fma1lLL6teQvXbRqeoZ6\n8vvRhGqLPxRFYdvJZHo19URtaolu5BxItsz31hLcFAnzG41Qq7Fv3ZqsLVsAcLvPdDlteQIBQ96H\ngkxZhbLifunpUM0Hq9pu81JZEjNCVqlfLSLjpw14PvM0XtOmEvTdtwiVDRcnTCD3QJR546iEoteT\ntmw59u3aVSsqaKMSNHbTcanSymNzdGK1drPuEycg1GrSFpsnb18jHSdK29seT5dtyjt4yOL5jlJs\nXFywa9bMuEhiA0BRFLIjI3Ho3h1hawtqLbQYCiPnwkunYPIm6TeSdk7a6M5qCXM6wMoHZAPsuHUw\ndjV4WLjXSAhZCffULjk5/T4dlo2AK5dqP7YWdF261E01uTAHVo+Xkja758OuL2o9ZGgbHy6l5RJ9\n2fi98OmkbC5n5JsuSbLnK9gyA6KWmXPl9Yp18jCRUp0rx/79jYZnaiWgMzx3VHa3x+6FJXfAooEy\nkV7JOrTabvPkE6CUmNxZnrllC8mzZuF85514PvkkIFdRQd99h42bK5emTCH777o3a+VERFB4/jzu\nEyfW2CMQ6KGrkvPYfCyBzkHueBgRgrP19MTl7tFk/Pij5XI0Dh4w4Qep+wQUXb5M8eXLaDtYPt9R\niq5zOHkHD1q28dFCFJ4/T3H8ZeMquiobWVI7+D14Jkr+kA+YIUN/Q/5PJrmbmRi2rSvOvjDueyny\nGB8FX/SU0kKWKuM2FUWRPTLJJ2DcGin3v+k/cPyXGg8b1LoRKgGbqmkY3HZSfq5NSpYf/0WKXLYY\nBgNN7xurb6yTh4mUJlVNTpQbw8ED+r0qJ5FhH0NuCqyZCHM7y9CAoWGq2m7zBNNlSfKjo4l/ZTr2\n7dvh+/57FX7cNQH+BH/3HZrgYGKeeorMjRvrNJy0pcuw9fbGecjgGvcLcpcrj9Il/MXUHE4kZNXo\nGOgxeTJKSQnpy+rnTivv4EGAa5ZhrwlteDj63Fzyj1cjzngDKXMN7FWLJIkQsiG1z0sw8UfoMdUs\nhYZrQggZbnxyJ/h1kD/iG1++vhPIri9k/mfADGkrMHqB7Nz/4VE5qVWDp6MdXYLdqy3Z3XYymRaN\nnPB10db8+pd2w7opUrF7zNdyYm8gWCcPE3EaPIjg79fgWNuXzRQ0Ouj6qOx4v3ep7ET/9QX4tA1s\n/xBy04x3m5fKktQSKihKSiLmqanYuLrSeO5coz0Mtp6eBC1birZ9O+JefIn0VavMGkLB6dPkREbi\nNm5crV26gR4OZBUUk54rq81KE4lDavDu0AQG4jx0COmrVlOSZXktpNyoKIS9PfYtW1r83KXUaA51\ng8mOjEATHIwm4Mb27ZqEWxBM3ADdp8Ler4w2O9YLFyJg8wxoeRfcZrCn1ejggZWg84AVD0BGbLWH\nD23jw6nEbM4mVyy7zykoZu+FtNpDVimnYeX9suR57Gr52g0I6+RhIkKlKgtdWQyVDYSNkj4Xk36V\ndxd/vQ+fhuEodgOVus0TDbIkNdx96PPziZ06jZLMTBp/Me+qVLwRbJycCPzqKxz79CHhrbdJmb/A\n5O7utGXLEXZ2uN53b637BlZS1918LJHWvs40dq/5y+DxyCPos7PNnthMIe/gIbRt2lhUnqIy6kaN\nUAcENJxmQQP6wkJy9+ytfdXRkFCpZN6wwzjY/j+ZSK9PMuLg+0lS3mbUl3IVVIpTIxi7RjY4rri/\nWqHH0pujylpXO8+mUlSi1CxJkpUI346RUi7j10mhzwaGdfJoCAghk4Pj1sj4cthoNBfWoHEuJnv1\nPLk8VhRDpVX1IStFUbj8+uvkHz2K/4czsW/VqtaXVmm1BMz9HOfhw0mePZukmR/WqEYKBs+ODRtw\nGTnSJB+EoHK9HslZBey/lF5jyKoU+9atcejZk7Sly9AXFNS6v6no8/PJj46u15BVKbrwcHL3H7Cc\n5IoFyDtwACUvr/5UdOsLIWQOpMWdMpF+eE39vE5xgQwnF+XBA98ZryRr1BruWwpJx6W5nJEyXj9X\nLe0bu1Yp2d1+KgmdxobOwdW4gBZkw4r7ICdZ/ia4h1hiVBbHOnk0NLxbSbvS5w7j2CWMnPNXKJnX\nHxYPhbw08Kl+9ZMy7wsyN/6G1wvPm96Hgqwm85v5P9zGjyftm2+4/J83akzyXlm9BqWgAPeJE0w6\nf/lejy3RiShKzSGr8ng8+gglKSlk/Ghah74p5B89CsXF9dIcWBlt53BK0tIotFTZsQXIiYgAtRqH\nrmY08zUUbGzhnsVS9PLHJ+HUJsu/xu+vQtw++T30alH9fqEDpMzKmS1yMjNygzA0zIfDsRnEpsuC\nkdIS3Z6hnmhsjfz8lhTJFU/CYbj3GxmNaKBYJ4+GirMfTg+/AXpBjvdDV0sVq/kwZW7cSMrcubiM\nGoXHI7X4XBtBqFQ0+s/reE6bRsb69cQ++5zRu32lsJD0FStw6NULu6ZNTTq3vdqGRs52XErLZdOx\nBALddbT0cTLpWF337ti3aUPq4q9Nc1k0gbJkeT00B1amwZhDlSM7cie6jh2rSMncNKjtZUd6ozC5\nQrj4j+XOfWA57FsszbRaj6x9/84PQ8+nYe8i2PVllaeHtpE3SZuOSWfJs8k5xKbnGc93KAr88pyc\njO76FJpbXrjRklgnjwZMWbd5vBaePQRPRMqS30rkHT5M/Guvo+3UCZ933q6ztLYQAq9pU2n0n/+Q\nvXUrMY89XkVjK3PTZoqTkso8O0wlyN2BY/GZ7DybwuDWjUy+RiEEHo88QtHFS2V9NtdKbtRBNEFB\n2LpXEzawIJomwdh4eDSYvEdxSgoFx4/fXPkOY9g7y14TlwCZd0iwgJxN3AHpvBjST1ZXmcrAd2RS\nfdPrcKJi5WITTwda+jiVlexuPyWFXY3mO7bPhKhvoc8rZRI6DRnr5NGAKes237EDRdgYbQ4sSkgg\nZupUbD09CZj7uUU8qN0njMfvow/J3bePSw9NojhNykuXeXY0aWK8P6AGAj10HL+cSVGJwpA25nV0\nOw0aiCYoiNSvrl2uXVEU8qKirsuqA+Tkp+vUqcGIJObslIrFN12+wxiOXlKZ185RdsKbYftchZxU\nuYpx9IYxi83zylGp4O6vZDnxuilVSniHhPmw92IaSVn5bDuZRKiXQ9VikQPLYNsHsiCg/+t1H8d1\nxDp5NHCc+vWjJD3dqLe5PjeXmKeeQsnNI+DLLyx6J+0yfDgB8+ZScOYMF8eNp+jyZfKiDpJ/9Kjs\nADdTw6jUz9zTUUOnQPO0fISNDe5TJpN/7BgJb75VxdXRHIouXaIkLe265DtK0XUOpygujqIE0x3m\n6ovsiAhs3N1NKqa4KXANlBOIvhiWjZIyKuaiL4F1kyE7Ce5bJvuxzEWjgwdXgda9Sgnv0DY+stfw\nYDy7z6fRr4V3xWNPbZYW2aG3S+8VS5py1SPWyaOB49C7t9Fuc0WvJ376qxScOIn/rE9ql4evA079\n+hH49SKKk5O5MHYcyZ9+isrFBZeRJsSCK1GqrjuodSNsVOZ/OVxHj8ZtwgSurFvH2cFDSFmwsE6K\ntVebA6/f5KEty3vc2NWHoteTs/MfHHr2rBcvixuGVwsYtxZyUqS4Y166ecf/+S6c2wZ3zZINgHXF\nyUdWRxXmVCjhbenjRLCHjjlbT1NYrK8Ysoo7AN8/JPM39y29fg2YFuBf9An6d1LabZ69bVuF7cmf\nzSFryxYaTX8Fx7596+31dZ07E7R8GUphIbl79+J2372odOY3K4X5uWCjEoxoX7emNKHrkAhwAAAM\n20lEQVRW4/Of1wn5+Wd03buT/OmnnL1jGBk//VRraXF5cqOiUDk4YNf0+vm/27dsgUqnu+HNggUn\nT1KSkvLvCFlVJiBcltWmnpE/3IXVu1dWIPoniPgUwh+WhlnXSqMwuO8bWcK7djKUFCOEYEgbHzLz\ni9GqbejaxBAhSDsvS3J1nlKKxc60IpKGgnXyuAlw6t+PgtNnKIyVS+GMn34idcECXO+9F7dq1Gwt\niX2rVgR/9y2u99+P+6RJdTpHU29HDr05mB6hdQgJlMMupAmN580lcNlSbN3diZ/+KhfuuZecXbtN\nOj4v6iDa9u0RNtdP5kHY2qLt2PGGiyTmRBokSXr2vKHXUW+E9ocxi6R23JqJtXuTJJ+UfiT+neGO\nmZa7jqYDYdhHcHozbHoNkCW7AD1CPbBX28gcy7djZLht/Dq5arnJsE4eNwFl3uZ//kXugSguvzED\nXdeu+Mx4o86VVeaiCQ7G9+23sPWo+49/ed+Oa8Wha1eCv1+D30cfUnwlnUuTJhHz5FMUnD1b7TEl\n2dkUnDp1XUNWpeg6h1Nw+jQlV2o2CKpPsiMisWvRArW3d+0736y0HinLXM/8AT8+Ub2ce36mVMpV\na2Wew7aqQOc10WUK9JgGexbCrvm0D3BlVAc/JvYIgsJcqU6cESvzJF6WDzlfDyz3bbZSb2gCA9GE\nhHDlx/UUz5+PrZ8v/p/NNs9T5F+IUKlwGT4cp0GDSFu+nNQFCzk3YiSu992L17RpVSa6/MOHQVFu\nyOShDZf9ORcnPYx9q1bYNQ1FExKCXWgoan9/i6+E9Hl5FF64QMG5cxSev0DhuXPk7t+P+wTTGjtv\nasqcI9+SpmnDPq6YhFYU2WCYehYe2gAu9aTvNegdGZr6/VVUbkHMfuAOmZxfM1Guju5bJtWLb1Lq\ndfIQQgwFPgNsgEWKovyv0vPC8PwwIBeYpCjKgfq8ppsVx/79SPt6MSonJxp/+aVJsiC3Cip7ezwf\nfRTXMWNImTuP9NWrydzwMx6PPYb7QxPLhCFzo6JACLTt2133a9R16oTHI1PIO3KU7Ii/yVi/vuw5\nYWeHpkkT7EJD0YSGYBfaFLvQEDSBgTXeICiKQklKCgXnzlN4/pycKM6dp/DcOYri46/uKATqgAAc\ne/XC9d576nOYDYdez0FuKuz8XIoYli9/jfgUTvwi5eWDzSs5NwuVDYz5CpYMg7VTYPJvsiT3xC/S\nCrj1iPp77euAqC/NHSGEDXAKGATEIj3NH1QUJbrcPsOAp5GTRzfgM0VRutV03s6dOyv7GlC37vUi\n/+RJYqdOw+fttyyj7PsvpuDceZI++YTsrVux9fXF+7lncR4+nJjHHqc4MZGQnzfc6EukJCODgrPn\nKDx3loIzZyk4d5bCs+coiou7upOtLZrAQDmRhISiCQykODWVwnPnKDgvJwp9OcVhodVi16QJmpAQ\nNCFNsAsJQdMkBE1wECo7C4dlbgYUBTZMk413Q2dKc6uzf8pcQ9hoKXF+PcK+mZdh0e2QdwWKcqDn\nMzC4jo6G14AQYr+iKFW7jOt6vnqcPHoAbymKMsTw92sAiqJ8UG6fBcA2RVFWGv4+CfRTFKXaYu1b\ndfKwYj45u/eQ9OGH5B87hn3r1hRevIjznXfi+07DMdSpjD43l4Lz5yk8e7bC5FJ46RIY5FlsGzWS\nk0OTEDlRNAnGLiQE20aN/l0luJagpFiWwp74RRopRX4mk9OP/AGa6yjPknBUrkCaD5GeIDfg/5Ol\nJ4/6DFv5A+XMKIhFri5q28cfqEOnjxUrFXHoJpPqmb/+StKsT9Hn5KALv4Y6/uuASqdDGxaGNqyi\nT71SWEhRfDw2np7YODreoKu7CbGxlSuM7+6Rvut2LnD/t9d34gCpDvHiCZmgv0maAGvjpkiYCyEe\nAx4DCAwMvMFXY+VmonxSPeeff3Ds3ftGX1KdEBoNmuDgG30ZNyelQoq/vwbt7rW877qpNDAzp2ul\nPiePOKBxub8DDNvM3QdFURYCC0GGrSx7mVZuBVT29jgZSp6t3ILYO8OoeTf6Kv5V1GfgbS/QTAjR\nRAihAR4AKmcqNwAThaQ7kFFTvsOKFStWrDQM6m3loShKsRBiGrAJWaq7WFGUY0KIJwzPzwc2Iiut\nziBLdR+ur+uxYsWKFSuWo15zHoqibEROEOW3zS/3WAGm1uc1WLFixYoVy2Ot67NixYoVK2ZjnTys\nWLFixYrZWCcPK1asWLFiNtbJw4oVK1asmI118rBixYoVK2ZTb9pW9YUQIhm4aMYhnkBKPV3OzYB1\n/Nbx38rjB+t7UDr+IEVRvGrb2VRuusnDXIQQ+ywpBnazYR2/dfy38vjB+h7U1/itYSsrVqxYsWI2\n1snDihUrVqyYza0weSy80Rdwg7GO/9bmVh8/WN+Dehn/vz7nYcWKFStWLM+tsPKwYsWKFSsW5l89\neQghhgohTgohzgghXr3R12MJhBCNhRB/CSGihRDHhBDPGra7CyG2CCFOG/7rVu6Y1wzvwUkhxJBy\n28OFEEcMz80R4uaxOBNC2AghooQQvxj+vmXGL4RwFUKsFUKcEEIcF0L0uMXG/7zhs39UCLFSCGH/\nbx+/EGKxECJJCHG03DaLjVkIYSeEWG3YvlsIEVzrRSmK8q/8h5SBPwuEABrgEND6Rl+XBcblC3Qy\nPHYCTgGtgQ+BVw3bXwVmGh63NozdDmhieE9sDM/tAboDAvgNuONGj8+M9+EFYAXwi+HvW2b8wFLg\nEcNjDeB6q4wfaVN9HtAa/l4DTPq3jx/oA3QCjpbbZrExA08B8w2PHwBW13pNN/pNqcc3uwewqdzf\nrwGv3ejrqodx/gQMAk4CvoZtvsBJY+NG+qv0MOxzotz2B4EFN3o8Jo45ANgKDCg3edwS4wdcDD+e\notL2W2X8/kAM4I60lPgFGHwrjB8IrjR5WGzMpfsYHtsimwpFTdfzbw5blX7ISok1bPvXYFhadgR2\nA42Uqy6MCUAjw+Pq3gd/w+PK228GZgOvAPpy226V8TcBkoElhrDdIiGEA7fI+BVFiQM+Bi4Bl5Hu\no5u5RcZfCUuOuewYRVGKgQzAo6YX/zdPHv9qhBCOwDrgOUVRMss/p8jbh39lGZ0Q4i4gSVGU/dXt\n828eP/KusBPwpaIoHYEcZMiijH/z+A1x/ZHISdQPcBBCjC+/z795/NVxI8b8b5484oDG5f4OMGy7\n6RFCqJETx3eKovxg2JwohPA1PO8LJBm2V/c+xBkeV97e0OkFjBBCXABWAQOEEN9y64w/FohVFGW3\n4e+1yMnkVhn/QOC8oijJiqIUAT8APbl1xl8eS4657BghhC0yPJpa04v/myePvUAzIUQTIf6/vbt5\njeoK4zj+/aGURCjSdicuVBAXsVaKigsXil2FbFpHKFgawb+g3YkQUFwoSKGQlNIgBkqxRRfWhYio\nf0C1JaTiC1oqCOLCglV8Q+RxcY7OZUwcD72TSTO/Dwxc7pm585y88HBe7n30DmkR6GSXY/rP8u6I\nw8CViPim0nQSGM7Hw6S1kJfnP8+7KZYDK4Hf8nD3vqSN+ZpfVj4zZ0XE7ohYGhHLSL/T8xHxBb3T\n/zvALUmr8qmtwGV6pP+k6aqNkhbluLcCV+id/lfV2efqtRqk/6s3j2S6vQjU4QWmQdJupL+APd2O\np6Y+bSINT6eAyfwaJM1PngOuA2eB9yuf2ZN/Bteo7CgB1gGXctsobRbI5toL2Exzwbxn+g+sBS7m\nv4ETwHs91v+9wNUc+4+kXUXzuv/AUdIazzPS6HNXnX0G+oBjwA3SjqwV7WLyHeZmZlZsPk9bmZlZ\nhzh5mJlZMScPMzMr5uRhZmbFnDzMzKyYk4eZmRVz8jArJGmnpCVvaN8n6ZM215iQ1Kg/OrPZsbDb\nAZj9D+0k3Wh1u7VB0oKIGJn1iMxmmUceZqQnFOfCSuO50NAZSf3TvK9Bukv3J0mTkvol3ZR0UNIf\nwPbqqELSiKQLuXDRD9MVHJJ0QKm415SkQx3vrFkNnDzMmlYCYxExANwDtrW+ISKOkx4NsiMi1kbE\n49z0T0R8HBE/t3xkNCLWR8RqoB8YqjZK+gD4FBiIiDXA/nq7ZNYZTh5mTX9HxGQ+/p1UfOdt/TLD\n+S25rOefpOJVAy3t/wJPgMOSPgMeFXynWdc4eZg1Pa0cP6dsTfBh6wlJfcB3QCMiPgTGSQ+geyVS\n4Z0NpEerDwGnC2M26wonD7NyD0j149t5mSju5uJdr+2uyucXR8Qp4Cvgo9qiNOsg77YyKzcBfC/p\nMak29LQi4p6kcdLOrDukGjOt3gV+zaMUAV/XH65Z/fxIdjMzK+ZpKzMzK+ZpK7MZSBoj1Uyv+jYi\njnQjHrO5xNNWZmZWzNNWZmZWzMnDzMyKOXmYmVkxJw8zMyvm5GFmZsVeABM19Fuf73rUAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113af7e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "mc_trails = performance_df2[performance_df2[\"method\"]==\"MC\"][\"n_trials\"]\n",
    "mc_err = performance_df2[performance_df2[\"method\"]==\"MC\"][\"err\"]\n",
    "ax.plot(mc_trails,mc_err,label=\"regular MC\")\n",
    "\n",
    "anti_mc_trails = performance_df2[performance_df2[\"method\"]==\"antithetic MC\"][\"n_trials\"]\n",
    "anti_mc_err = performance_df2[performance_df2[\"method\"]==\"antithetic MC\"][\"err\"]\n",
    "ax.plot(anti_mc_trails,anti_mc_err,label=\"antithetic regular MC\")\n",
    "\n",
    "ohmc_trails = performance_df2[performance_df2[\"method\"]==\"OHMC\"][\"n_trials\"]\n",
    "ohmc_err = performance_df2[performance_df2[\"method\"]==\"OHMC\"][\"err\"]\n",
    "ax.plot(ohmc_trails,ohmc_err,label=\"OHMC\")\n",
    "\n",
    "anti_ohmc_trails = performance_df2[performance_df2[\"method\"]==\"antithetic OHMC\"][\"n_trials\"]\n",
    "anti_ohmc_err = performance_df2[performance_df2[\"method\"]==\"antithetic OHMC\"][\"err\"]\n",
    "ax.plot(anti_ohmc_trails,anti_ohmc_err,label=\"antithetic OHMC\")\n",
    "\n",
    "ax.legend()\n",
    "plt.xlabel(\"n_trials\")\n",
    "plt.ylabel(\"err\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Runtime Efficiency Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X94VdWd7/H3F4QCgmglFCsgqCi/BEyCopYpFlSqFaF2\nNNQ7jr32AVodLU/LLa0OUms7dtJbMFon0NZH7XSkthUaldaWqQgXZJofRQVEiaAmIBKxoGD4EfK9\nf5zD9pxwQk5I9vmVz+t58pi99jp7r2wP53v2Wmt/l7k7IiIiAJ3S3QAREckcCgoiIhJQUBARkYCC\ngoiIBBQUREQkoKAgIiIBBQUREQkoKIiISEBBQUREAieluwGt1adPHx80aFC6myEiklUqKyvfc/e8\nluplXVAYNGgQFRUV6W6GiEhWMbO3kqmn7iMREQkoKIiISEBBQUREAlk3piDSkRw+fJja2loOHDiQ\n7qZIlujWrRv9+/enS5cuJ/R6BQWRDFZbW0uvXr0YNGgQZpbu5kiGc3d2795NbW0tgwcPPqFjqPtI\nJIMdOHCA008/XQFBkmJmnH766W26sww1KJjZZDN7zcyqzWxugv29zexpM3vJzDaa2VfCbI9INlJA\nkNZo6/sltKBgZp2BnwKfB4YD081seJNqtwGb3H00MAH4v2bWNaw2iYjI8YV5p3ARUO3uW939ELAE\nuK5JHQd6WSS09QTeBxpCbJOIZJCePXu26fUTJkxg4MCBxK41P3Xq1Ljjvv7661x99dUMGTKE/Px8\nbrjhBt599902nTeXhRkUzgRqYrZro2WxHgKGATuAV4A73b2x6YHMbIaZVZhZRV1dXZsaVV62iJ3z\nz6Xxnt7snH8u5WWL2nQ8kY7C3WlsPOafZ9rPf+qpp7JmzRoA9uzZwzvvvBPsO3DgANdccw1f+9rX\n2LJlC1VVVXz961+nrZ8juSzdA81XAeuBTwNjgIfM7JSmldx9sbsXunthXl6LqTuaVV62iJGVd9OP\nOjoZ9KOOkZV3KzBIzlj2t+1cdv9fGDz3WS67/y8s+9v2Nh3vzTff5Pzzz+fmm29m5MiR1NTU8Kc/\n/YlLLrmE/Px8/vEf/5F9+/YBsHz5coYOHUpBQQF33HEHX/jCFwCYP38+P/7xj4Njjhw5kjfffDPu\nPPv27WPixInk5+dzwQUX8Pvf/77Z8zdVVFTEkiVLAHjqqaf44he/GOz7r//6Ly655BKuvfbaoGzC\nhAmMHDmyTdcll4UZFLYDA2K2+0fLYn0FeMojqoFtwNCwGjSgqpjudiiurLsdYkBVcVinFEmZZX/b\nzneeeoXte+pxYPueer7z1CttDgxbtmzh61//Ohs3buTkk0/mvvvuY8WKFVRVVVFYWMhPfvITDhw4\nwMyZM/nDH/5AZWVlq7+Jd+vWjaVLl1JVVcXzzz/PN7/5zaBLKPb8Z5111jGvnThxIqtWreLIkSMs\nWbKEG2+8Mdi3YcMGCgoK2vT3dzRhBoVyYIiZDY4OHhcBZU3qvA1MBDCzTwHnA1vDalBfT/xG7evv\nhXVKkZQpfu416g8fiSurP3yE4udea9NxzzrrLMaNGwfAunXr2LRpE5dddhljxozhscce46233mLz\n5s2cffbZwdz46dOnt+oc7s53v/tdRo0axaRJk9i+fXvQ7x97/kQ6d+7MZz7zGZYsWUJ9fT3Kotw2\noT285u4NZnY78BzQGXjE3Tea2azo/lLg+8CjZvYKYMC33cP7hN5lefTj2MCwy/rQL6yTiqTIjj31\nrSpP1sknnxz87u5cccUVPPHEE3F11q9f3+zrTzrppLixgERz6H/1q19RV1dHZWUlXbp0YdCgQUG9\n2PM3p6ioiGnTpjF//vy48hEjRvDCCy+0+Hr5WKhjCu6+3N3Pc/dz3P0H0bLSaEDA3Xe4+5XufoG7\nj3T3/wyzPTX5c6j3+Bmv9d6Vmvw5YZ5WJCU+fWr3VpWfiHHjxrFmzRqqq6sB2L9/P6+//jrnn38+\nW7duDcYKfv3rXwevGTRoEFVVVQBUVVWxbdu2Y467d+9e+vbtS5cuXXj++ed5662ksjwHxo8fz3e+\n851j7lC+/OUvs3btWp599tmgbNWqVWzYsKFVx+9I0j3QnFJjp8xkQ8F97CSPRjd2kseGgvsYO2Vm\nupsm0mZzrjqf7l06x5V179KZOVed327nyMvL49FHH2X69OmMGjWKSy65hM2bN9O9e3cefvhhJk+e\nTEFBAb169aJ3794AXH/99bz//vuMGDGChx56iPPOO++Y4950001UVFRwwQUX8PjjjzN0aOuGFs2M\nb33rW/Tp0yeuvHv37jzzzDM8+OCDDBkyhOHDh/Pwww/Tlgkruc5i5/dmg8LCQtciO9JRvPrqqwwb\nNizp+sv+tp3i515jx556Pn1qd+ZcdT5TL2w6Ezwc+/bto2fPnrg7t912G0OGDGH27NkpObfES/S+\nMbNKdy9s6bVKiCeSQ6ZeeGbKgkBTP/vZz3jsscc4dOgQF154ITNn6g48GykoiEi7mD17tu4MckCH\nGlMQEZHjU1AQEZGAgoKIiAQUFEREJKCgICLtatmyZWzatCnYnjdvHitWrABg4cKFfPTRR8G+1qbO\nXrlyJWvXrg22S0tLefzxx9vY4vaRK2nAFRREpF01DQr33nsvkyZNAo4NCq3VNCjMmjWLm2++udXH\nURrw5ikoiOSSl5+EBSNh/qmR/778ZJsPOXXqVAoKChgxYgSLFy8Oynv27Mldd93F6NGjGTduHO++\n+y5r166lrKyMOXPmMGbMGN544w1uueUWfvvb31JSUsKOHTu4/PLLufzyy4PjND0GQF1dHddffz1j\nx45l7NixrFmzhjfffJPS0lIWLFjAmDFjWL16dVxa7urqaiZNmsTo0aPJz8/njTfeiPs7lAY8Se6e\nVT8FBQUu0lFs2rQp+cov/dr9vk+533PKxz/3fSpS3ga7d+92d/ePPvrIR4wY4e+99567uwNeVlbm\n7u5z5szx73//++7u/s///M/+m9/8Jnh97PZZZ53ldXV1wb7mjjF9+nRfvXq1u7u/9dZbPnToUHd3\nv+eee7y4uDh4fez2RRdd5E89FcnEX19f7/v374/7O7Zt2+Zm5i+++KK7u9fV1fn48eN937597u5+\n//33+/e+9z2vr6/3/v37+9atW93dvaioyK+55pqE5x8xYoRv27bN3d1PPvlkd3c/fPiw7927NzjH\nOeec442Njcecv6nPfvazvm7dOr/gggu8oaHBr7jiCt+2bVtw3NmzZ/vChQsTvrapRO8boMKT+IzV\nw2siueK/74XDTTKiHq6PlI+64YQPW1JSwtKlSwGoqalhy5YtnH766XTt2jX4Bl1QUMCf//znVh+7\nuWOsWLEirgvqgw8+CL7FJ/Lhhx+yfft2pk2bBkTWZ0ikuTTgAIcOHQpyOTVNAx57h9QSj6YBX7Vq\nFZ06dcq6NOAKCiK5Ym9t68qTsHLlSlasWMGLL75Ijx49mDBhQpDSukuXLkSWV498mDU0tH559eaO\n0djYyLp165r9cD9RSgPeMo0piOSK3v1bV56EvXv3ctppp9GjRw82b97MunXrWnxNr169+PDDD1u9\nL9aVV17Jgw8+GGwf/aBu7vW9evWif//+LFu2DICDBw+2OKCtNOCJKSiI5IqJ86BLk7UTunSPlJ+g\nyZMn09DQwLBhw5g7d+5xuz6OKioqori4mAsvvPCYwd4ZM2YwefLkuIHmREpKSqioqGDUqFEMHz6c\n0tJSAK699lqWLl0aDDTH+uUvf0lJSQmjRo3i0ksvZefOncc9h9KAN3N+DzF1tplNBh4gsvLaz939\n/ib75wA3RTdPAoYBee7+fnPHVOps6Uhamzqbl5+MjCHsrY3cIUyc16bxhI4q29OAZ2TqbDPrDPwU\nuAKoBcrNrMzdg9Ejdy8GiqP1rwVmHy8giEgLRt2gINAOOnIa8DAHmi8Cqt19K4CZLQGuAzY1U386\n8EQz+0REUqYjpwEPc0zhTCD26YzaaNkxzKwHMBn4XTP7Z5hZhZlVtPfTeyIi8rFMGWi+FljTXNeR\nuy9290J3L9TaqiIi4QkzKGwHBsRs94+WJVKEuo5ERNIuzKBQDgwxs8Fm1pXIB39Z00pm1hv4LPD7\nENsiIiJJCC0ouHsDcDvwHPAq8KS7bzSzWWY2K6bqNOBP7r4/rLaISNvU1tZy3XXXMWTIEM455xzu\nvPNODh06xMqVK4M0FUcdTYAHmZMOWpIX6piCuy939/Pc/Rx3/0G0rNTdS2PqPOruRWG2IyzlZYvY\nOf9cGu/pzc7551JetijdTRJpd+7OF7/4RaZOncqWLVt4/fXX2bdvH3fddVdSr8+EdNCSvEwZaM46\n5WWLGFl5N/2oo5NBP+oYWXm3AoOk1bNbn+XK317JqMdGceVvr+TZrc+2/KIW/OUvf6Fbt2585Stf\nASI5ihYsWMAjjzyS1NoIGZEOWpKmoNCMlu4CBlQV090OxZV1t0MMqCpOZTNFAs9ufZb5a+fzzv53\ncJx39r/D/LXz2xwYNm7cSEFBQVzZKaecwsCBA6murmb16tWMGTMm+Ckrix86nDhxIqtWreLIkSMs\nWbKEG2+8Mdi3YcOGY44t6aWgkEAydwF9PfHtbV9/L1XNFInzQNUDHDgSn7XzwJEDPFD1QKjnHT9+\nPOvXrw9+pkyZErc/E9JBS/IUFBJI5i5glyV+XmKX9UlYLhK2nfsTJ4BrrjxZw4cPp7KyMq7sgw8+\n4O233+bcc89N6hhFRUXccccd3HBDfAqOESNGHHNsSS8FhQSSuQuoyZ9DvXeN21/vXanJnxNq20Sa\n0+/kfq0qT9bEiRP56KOPePzxxwE4cuQI3/zmN7nlllvo0aNHUsdIdzpoSZ6CQgLJ3AWMnTKTDQX3\nsZM8Gt3YSR4bCu5j7JSOkzhLMsud+XfSrXP8ojTdOnfjzvw723RcM2Pp0qX85je/YciQIZx33nl0\n69aNH/7wh606RjrTQUvyQk2dHYZUpM4+OqYQ24VU7131oS8p19rU2c9ufZYHqh5g5/6d9Du5H3fm\n38k1Z18TYgslE2Vk6uxsNnbKTMqJjC309ffYZX2oKZijgCAZ75qzr1EQkDZRUGjG2CkzIRoE+kV/\nRERyncYUREQkoKAgIiIBBQUREQkoKIiISEBBQUTa1bJly9i06eOl2OfNm8eKFSsAWLhwYVwSvdgU\n2slYuXIla9euDbZLS0uDh+qS0VwK8KPHVhpwBQURaWdNg8K9997LpEmTgGODQms1DQqzZs3i5ptv\nTuq1bU0BDh0jDbiCgkgO2fv002z53EReHTacLZ+byN6nn27zMadOnUpBQQEjRoxg8eLFQXnPnj25\n6667GD16NOPGjePdd99l7dq1lJWVMWfOHMaMGcMbb7wRfNsuKSlhx44dXH755Vx++eXBcZoeA6Cu\nro7rr7+esWPHMnbsWNasWcObb75JaWkpCxYsYMyYMaxevZr58+fz4x//GIDq6momTZrE6NGjyc/P\n54033oj7O9qaAhw6RhrwUIOCmU02s9fMrNrM5jZTZ4KZrTezjWb2QpjtEclle59+mnf+dR4NO3aA\nOw07dvDOv85rc2B45JFHqKyspKKigpKSEnbv3g3A/v37GTduHC+99BL/8A//wM9+9jMuvfRSpkyZ\nQnFxMevXr+ecc84JjnPHHXfw6U9/mueff57nn3++2WMA3HnnncyePZvy8nJ+97vf8dWvfpVBgwYx\na9YsZs+ezfr16xk/fnxcO2+66SZuu+02XnrpJdauXcsZZ5wRt7+lFOCA0oAT4sNrZtYZ+ClwBVAL\nlJtZmbtviqlzKvAwMNnd3zazvmG1RyTX7VqwED8QnzrbDxxg14KF9I759tpaJSUlLF26FICamhq2\nbNnC6aefTteuXYM++IKCAv785z+3+tjNHWPFihVxXVAffPAB+/bta/Y4H374Idu3b2fatGkAdOvW\nrdm6xzN+/HieeeaZYPuWW26J298R0oCH+UTzRUC1u28FMLMlwHXAppg6Xwaecve3Adx9V4jtEclp\nDTH928mUJ2PlypWsWLGCF198kR49ejBhwgQORANPly5dMDMg8mHZ0NDQ6uM3d4zGxkbWrVt3wh/u\niQwfPjwYND4qNgX4X//616SOU1RUxLRp05g/f35c+YgRI3jhhezv7Aiz++hMoCZmuzZaFus84DQz\nW2lmlWaWcMTIzGaYWYWZVWTboI1IqpzUpLukpfJk7N27l9NOO40ePXqwefNm1q1b1+JrevXqxYcf\nftjqfbGuvPJKHnzwwWB7/fr1x319r1696N+/P8uWLQPg4MGDx4wTtEcKcMj9NODpHmg+CSgArgGu\nAv7VzM5rWsndF7t7obsXKqWuSGJ9Z38Da/LN2rp1o+/sb5zwMSdPnkxDQwPDhg1j7ty5jBs3rsXX\nFBUVUVxczIUXXnjMYO+MGTOYPHly3EBzIiUlJVRUVDBq1CiGDx9OaWkpANdeey1Lly4NBppj/fKX\nv6SkpIRRo0Zx6aWXsnNn/OJC7ZEC/OhxcjkNeGips83sEmC+u18V3f4OgLv/W0yduUB3d78nuv0L\n4I/u/pvmjpuK1NkimaK1qbP3Pv00uxYspOGddzjpjDPoO/sbbRpPkOyUqamzy4EhZjYY2A4UERlD\niPV74CEzOwnoClwMLAixTSI5rfe11yoISJuEFhTcvcHMbgeeAzoDj7j7RjObFd1f6u6vmtkfgZeB\nRuDn7p5dHXAiIjkk1PUU3H05sLxJWWmT7WKgOMx2iGQzdw9m6Ii0pK1DAukeaBaR4+jWrRu7d+9u\n8z906Rjcnd27d7dpKq9WXhPJYP3796e2tjbr8udI+nTr1o3+/fuf8OsVFEQyWJcuXRg8eHC6myEd\niLqPREQkoKAgIiIBBQUREQkoKIiISEBBQUREAgoKIiISUFAQEZGAgoKIiAQUFEREJKCgICIiAQUF\nEREJKCiIiEgg1KBgZpPN7DUzq44uvdl0/wQz22tm66M/88Jsj4iIHF9oWVLNrDPwU+AKoBYoN7My\nd9/UpOpqd/9CWO1Ih/KyRQyoKqav17HL8qjJn8PYKTPT3SwRkRaFeadwEVDt7lvd/RCwBLguxPNl\nhPKyRYysvJt+1NHJoB91jKy8m/KyRelumohIi8IMCmcCNTHbtdGypi41s5fN7A9mNiLE9qTEgKpi\nutuhuLLudogBVVpxVEQyX7oX2akCBrr7PjO7GlgGDGlaycxmADMABg4cmNoWtlJfr4MEy+n29fdS\n3xgRkVYK805hOzAgZrt/tCzg7h+4+77o78uBLmbWp+mB3H2xuxe6e2FeXl6ITW67XZa4fbuO/bNE\nRDJOmEGhHBhiZoPNrCtQBJTFVjCzfmZm0d8virZnd4htCl1N/hzqvWtcWb13pSZ/TppaJCKSvNC6\nj9y9wcxuB54DOgOPuPtGM5sV3V8KfAn4mpk1APVAkbt7WG1KhbFTZlIO0dlH77HL+lBToNlHIpId\nLNs+gwsLC72ioiLdzRARySpmVunuhS3V0xPNIiISUFAQEZGAgoKIiAQUFEREJKCgICIiAQUFEREJ\nKCiIiEhAQUFERAIKCiIiElBQEBGRgIKCiIgEFBRERCSgoCAiIoEWg4KZdTKzG1LRGBERSa8Wg4K7\nNwL/JwVtERGRNEu2+2iFmX3LzAaY2SeP/oTasiyy8hf3subikWwcOow1F49k5S/uTXeTREROSLJB\n4UbgNmAVUBn9aXGlGzObbGavmVm1mc09Tr2xZtZgZl9Ksj0ZY+Uv7uXUhU/wyb1H6AR8cu8RTl34\nhAKDiGSlpMYUgP/l7oOb/Jzdwus6Az8FPg8MB6ab2fBm6v0I+NMJ/QVp1mXxk3zicHzZJw5HykVE\nsk2yYwoPncCxLwKq3X2rux8ClgDXJaj3L8DvgF0ncI60O3XvkVaVi4hksmS7j/7bzK43M2vFsc8E\namK2a6NlATM7E5gG/EcrjptR9vTu3KpyEZFMlmxQmAk8CRw0sw/M7EMz+6Adzr8Q+Hb0bqRZZjbD\nzCrMrKKurq4dTtt+Ds+4gYNd4ssOdomUi4hkm5OSrNcbuAkY7O73mtlA4IwWXrMdGBCz3T9aFqsQ\nWBK9AekDXG1mDe6+LLaSuy8GFgMUFhZ6km1OiQm3zmMlkTGEU/ceYU/vzhyecQMTbp2X7qaJiLSa\nubf8GWtm/wE0Ap9z92FmdhrwJ3cfe5zXnAS8DkwkEgzKgS+7+8Zm6j8KPOPuvz1eWwoLC72iosWJ\nTyIiEsPMKt29sKV6yd4pXOzu+Wb2NwB3/7uZdT3eC9y9wcxuB54DOgOPuPtGM5sV3V+a5LlFRCRF\nkg0Kh6NTRx3AzPKI3Dkcl7svB5Y3KUsYDNz9liTbIiIiIUl2oLkEWAr0NbMfAP8P+GForRIRkbRI\n6k7B3X9lZpVExgcMmOrur4baMhERSblku49w983A5hDbIiIpsPIX92q2nDRL6ymIdCDK1SUtUVAQ\n6UCUq0taknT3kUgmKS9bxICqYvp6Hbssj5r8OYydMjPdzcp4ytUlLdGdgmSd8rJFjKy8m37U0cmg\nH3WMrLyb8rJF6W5axlOuLmmJgoJknQFVxXS3Q3Fl3e0QA6qK09Si7KFcXdISBQXJOn09cVLEvv5e\niluSfSbcOo8935jO+7070wi837sze74xXbOPJKAxBck6uyyPfhwbGHZZH/qloT3ZZsKt80BBQJqh\nOwXJOjX5c6j3+NRb9d6Vmvw5aWqRSO7QnYK0m1TNCBo7ZSblED3Xe+yyPtQUaPaRSHtIKnV2JlHq\n7Mx0dEZQ7ABwvXdlQ8F9+rAWyQDJps5W95G0C80IEskNCgrSLjQjSCQ3KChIu9hlec2U90lxS0Sk\nLUINCmY22cxeM7NqM5ubYP91Zvayma03swoz+0yY7ZHwaEaQSG4IbfZRdKW2nwJXALVAuZmVufum\nmGr/DZS5u5vZKOBJYGhYbZLwaEaQSG4Ic0rqRUC1u28FMLMlwHVAEBTcfV9M/ZOJLvcp2WnslJkQ\nDQL9oj9KXCdavyG7hNl9dCZQE7NdGy2LY2bTzGwz8Czwv0Nsj6SYEteJ1m/IPmkfaHb3pe4+FJgK\nfD9RHTObER1zqKirSzzLRTJPrk1TXfmLe1lz8Ug2Dh3GmotH6oMtCVq/IfuEGRS2AwNitvtHyxJy\n91XA2WbHTldx98XuXujuhXl5iWe5SOq19CGZS9NU9Y33xGj9huwTZlAoB4aY2WAz6woUAWWxFczs\nXDOz6O/5wCeA3SG2SdpJMh+SuTRNVd94T4zWb8g+oQUFd28AbgeeA14FnnT3jWY2y8xmRatdD2ww\ns/VEZird6NmWd6ODSuZDMpemqeob74nR+g3ZJ9SEeO6+HFjepKw05vcfAT8Ksw0SjmQ+JNM1TTWM\nGU97enfmkwn+Zn3jPb4Jt85jJWj2URZRllQ5Icl+SCaaphqmuMR80RlPvSvvpvxoW07Q4Rk3cHDh\nE3F3R/rGe6xmp58qCGSNtM8+kuyUqd0CYc140oplLdNgfG5Q6mw5YZn4UFLjPb3pZAnK3ej0vT2p\nb1AHsubikQnvHt/v3ZnL/mdDGloksZJNna3uIzlhmdgtoKU600eD8blB3UdZRg9QNa+8bBGf4ABN\nb36zdcZTU5n+/17TT3ODgkIWUZ9t844OMJ/Gh1i0+8gd/k7PnFj9LRv+32fqOFNrZHrgTQUFhRZk\n0ptED1A1L9EAsxkcpHvWBwTIjv/32T4Ynw2BNxU0pnAcR98kR/8xfnLvEQ4ufIKVkJY3ekfvsz3e\n8wd9vQ4SDDBnY0qNRLLl/30mjjMl67iBN0v/phOhO4XjyLRvZ831zf79lNz/39hSxtVcSqmRiPrr\nw5ctgTdsuf9p0gaZ9iY5POMGDja5tzt4Ehy48KOcT0fd0vMHuZRSI5Fc6K/PdAq8EQoKx5Fpb5IJ\nt85jx2WNvH8KkT7bU2DPZfuZdMb7WZuOOlktZVwdO2UmGwruYyd5NLqxk7ycGGA+Ktv767OBAm+E\nxhSOIxNTG0zut5NOVx9bnit9581J5vmD9kipkckrxWVzf302UJ6mCAWF48jEN0muPZyV7IdwTf4c\neh/NaRRV712pKZjTbn93WHmTpHXS+aR8JgXedF0HpbnIMnEfXFH13jUru0pa+7d8HECiGVfb6Vv8\n0eN+yuuCZxxi7SSPfvOr23weaVnTGX8QuTvvaF1lYVyHZNNcKCikSHtG/bA+HFNt5/xzE971pPJD\nOFFgakp5k9qmNe995U+KCOM6KPdRBmnv5x1SnY46LJnwbEGiWU1NZWvXXCZo7Xs/02b8pUs6r0Oo\ns4/MbLKZvWZm1WY2N8H+m8zsZTN7xczWmtnoMNuTLpn2vEOmyIRnC5qb1XRULk1rTYfWvvczbcZf\nuqTzOoQWFMysM5ElNj8PDAemm9nwJtW2AZ919wuA7wOLw2pPOunbT2Lt+WzBiaYjaS4wuZNz01rT\nobXvfU0LjUjndQjzTuEioNrdt7r7IWAJcF1sBXdf6+5/j26uA/qH2J42KS9bxM7559J4T292zj+3\nVQ+L6dtPYu31bEFbctY0F5gqCv6dfvOrFRDaqLXvfT2PEZHO6xDaQLOZfQmY7O5fjW7/E3Cxu9/e\nTP1vAUOP1m9OOgaa2zrjRzMqwtXWQblcGbjPRHrvZ46sGmg2s8uBW4HPNLN/BjADYODAgSlsWcRx\nUywk8eGRic875JK2ds/lysB9qiUzq0jv/ewT5p3CJcB8d78quv0dAHf/tyb1RgFLgc+7++stHTcd\ndwpa4jGzaRpj6ukOIPske6cQ5phCOTDEzAabWVegCCiLrWBmA4GngH9KJiCkSybMkpHmaXAy9TSj\nLneFFhTcvQG4HXgOeBV40t03mtksM5sVrTYPOB142MzWm1nKn0pLZgA51zNwZjsNTqaeZtTlrg79\nRHNrBpCzcTAynTlkJLepyy77ZNVAc7q0ZgA52wYjM23VOMktmZhBWNpHh15PoaUc/dlMfb4dT1ue\npWktddnlrg59p5Bpaajbs7tHfb4dSzrSfmdSmmlpPx36TiGTBpDb8lRuInqKumNpablSkWR16KCQ\nSUs4tndciwS/AAAIw0lEQVR3j6Zpdiy53BUqqdWhu48gcwaQ27u7R0+SdiyZ1hUq2avDB4VMsad3\n54RT/NrS3aM+344jFcuVSsfQobuPMom6e6QtMqkrVLJbh354LdPoYTPJVnrvZj6t0SwiKaHkeNkh\nExLiiUgH0FEflDzR1f4ynQaaRULwca6sOnZZXlbkyjpRHfFByVxOI6M7BZF2dvTp4n7U0Sn6dPHI\nyrtDTTuRTh3xQclcvjtSUBBpZx3t6eKOOHMul++OFBRE2llHe7q4IybHy+W7I40piLSzjvh0cUd7\nUDKXU4eHeqdgZpPN7DUzqzazuQn2DzWzF83soJl9K8y2iKRKJiValHDk8t1RaM8pmFln4HXgCqCW\nyJrN0919U0ydvsBZwFTg7+7+45aOq+cUJBtk40p9ktsyYeW1i4Bqd98abdAS4DogCAruvgvYZWbX\nhNgOkZTLlESLIq0VZvfRmUBNzHZttExERDJUVsw+MrMZZlZhZhV1dYlndoiISNuFGRS2AwNitvtH\ny1rN3Re7e6G7F+bl5bVL40RE5FhhBoVyYIiZDTazrkARUBbi+UREpI1CG2h29wYzux14DugMPOLu\nG81sVnR/qZn1AyqAU4BGM/sGMNzdPwirXSIi0rxQH15z9+XA8iZlpTG/7yTSrSQiIhkgKwaaRUQk\nNRQUREQkoKAgIiIBBQUREQkoKIiISEBBQUREAgoKIiISUFAQEZGAgoKIiAS0HGcG+HhBljp2WZ4W\nZBGRtFFQSLPyskWMrLyb7nYIDPpRR+/KuykHBYYoBU2R1FH3UZoNqCqOBIQY3e0QA6qK09SizHI0\naPajjk7RoDmy8m7Kyxalu2kiOUlBIc36euJFg/r6eyluSWZS0BRJLQWFNNtliRcN2mV9UtySzKSg\nKZJaCgppVpM/h3rvGldW712pyZ+TphZlFgVNkdRSUEizsVNmsqHgPnaSR6MbO8ljQ8F9GkiNUtAU\nSS1z9/AObjYZeIDIyms/d/f7m+y36P6rgY+AW9y96njHLCws9IqKipBaLJno49lH77HL+mj2kcgJ\nMLNKdy9sqV5oU1LNrDPwU+AKoBYoN7Myd98UU+3zwJDoz8XAf0T/KxIYO2UmRINAv+iPiIQjzO6j\ni4Bqd9/q7oeAJcB1TepcBzzuEeuAU83sjBDbJCIixxFmUDgTqInZro2WtbYOZjbDzCrMrKKuLvFs\nFBERabusGGh298XuXujuhXl5iWejiIhI24UZFLYDA2K2+0fLWltHRERSJMygUA4MMbPBZtYVKALK\nmtQpA262iHHAXnd/J8Q2iYjIcYQ2+8jdG8zsduA5IlNSH3H3jWY2K7q/FFhOZDpqNZEpqV8Jqz0i\nItKyULOkuvtyIh/8sWWlMb87cFuYbRARkeRlxUCziIikhoKCiIgEQk1zEQYzqwPeasMh+gBKsRmh\naxFP1yOerke8bL8eZ7l7i3P6sy4otJWZVSST/6Mj0LWIp+sRT9cjXke5Huo+EhGRgIKCiIgEOmJQ\nWJzuBmQQXYt4uh7xdD3idYjr0eHGFEREpHkd8U5BRESakZNBwcwmm9lrZlZtZnMT7DczK4nuf9nM\n8tPRzlRJ4nrcFL0Or5jZWjMbnY52pkpL1yOm3lgzazCzL6WyfamWzPUwswlmtt7MNprZC6luYyol\n8e+lt5k9bWYvRa9HbqXncfec+iGSZ+kN4GygK/ASMLxJnauBPwAGjAP+J93tTvP1uBQ4Lfr75zv6\n9Yip9xciaVq+lO52p/n9cSqwCRgY3e6b7nan+Xp8F/hR9Pc84H2ga7rb3l4/uXinoBXf4rV4Pdx9\nrbv/Pbq5jkgK81yVzPsD4F+A3wG7Utm4NEjmenwZeMrd3wZw91y+JslcDwd6RdeY70kkKDSktpnh\nycWg0G4rvuWI1v6ttxK5i8pVLV4PMzsTmEZkzfBcl8z74zzgNDNbaWaVZnZzylqXeslcj4eAYcAO\n4BXgTndvTE3zwhdqllTJLmZ2OZGg8Jl0tyXNFgLfdvfGyJfBDu8koACYCHQHXjSzde7+enqblTZX\nAeuBzwHnAH82s9Xu/kF6m9U+cjEoaMW3eEn9rWY2Cvg58Hl3352itqVDMtejEFgSDQh9gKvNrMHd\nl6WmiSmVzPWoBXa7+35gv5mtAkYDuRgUkrkeXwHu98igQrWZbQOGAn9NTRPDlYvdR1rxLV6L18PM\nBgJPAf/UAb79tXg93H2wuw9y90HAb4Gv52hAgOT+vfwe+IyZnWRmPYCLgVdT3M5USeZ6vE3krgkz\n+xRwPrA1pa0MUc7dKbhWfIuT5PWYB5wOPBz9dtzgOZr4K8nr0WEkcz3c/VUz+yPwMtAI/NzdN6Sv\n1eFJ8v3xfeBRM3uFyAzGb7t7NmdPjaMnmkVEJJCL3UciInKCFBRERCSgoCAiIgEFBRERCSgoiIhI\nQEFBpI3MbJCZfTlmu9DMStLZJpETpSmpIjGiSc6sNblszGwC8C13/0JoDRNJEd0pSIcX/ab/mpk9\nDmwAjsTs+5KZPRr9/dHoOhxrzWxrzDoL9wPjo+sNzI6uPfBM9DXzzewxM1ttZm+Z2RfN7N+ja1f8\n0cy6ROsVmNkL0YRzz+Vw1l7JcAoKIhFDgIfdfQSw/zj1ziCSMPALRIIBwFxgtbuPcfcFCV5zDpHk\naVOA/wSed/cLgHrgmmhgeJDIug0FwCPAD9rhbxJptZxLcyFygt6Krq3RkmXRrqVN0bw3yfiDux+O\npkXoDPwxWv4KMIhI7pyRRLJtEq2Tq7m4JMMpKIhExN4dxA60dWtS72DM78nm1T4IEE3Ffdg/Hshr\nJPJv0ICN7n5JK9orEgp1H4kc610zG2ZmnYgsttOSD4FebTjfa0CemV0CYGZdzGxEG44ncsIUFESO\nNRd4BlhLct04LwNHogu5z27tyaLLPn4J+JGZvURkAZdLW3sckfagKakiIhLQnYKIiAQUFEREJKCg\nICIiAQUFEREJKCiIiEhAQUFERAIKCiIiElBQEBGRwP8HyQBwh3LWOT8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11745b208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "mc_runtime = performance_df2[performance_df2[\"method\"]==\"MC\"][\"runtime\"]\n",
    "mc_err = performance_df2[performance_df2[\"method\"]==\"MC\"][\"err\"]\n",
    "ax.scatter(mc_runtime,mc_err,label=\"regular MC\")\n",
    "\n",
    "anti_mc_runtime = performance_df2[performance_df2[\"method\"]==\"antithetic MC\"][\"runtime\"]\n",
    "anti_mc_err = performance_df2[performance_df2[\"method\"]==\"antithetic MC\"][\"err\"]\n",
    "ax.scatter(mc_runtime,mc_err,label=\"antithetic regular MC\")\n",
    "\n",
    "ohmc_runtime = performance_df2[performance_df2[\"method\"]==\"OHMC\"][\"runtime\"]\n",
    "ohmc_err = performance_df2[performance_df2[\"method\"]==\"OHMC\"][\"err\"]\n",
    "ax.scatter(ohmc_runtime,ohmc_err,label=\"OHMC\")\n",
    "\n",
    "anti_ohmc_runtime = performance_df2[performance_df2[\"method\"]==\"antithetic OHMC\"][\"runtime\"]\n",
    "anti_ohmc_err = performance_df2[performance_df2[\"method\"]==\"antithetic OHMC\"][\"err\"]\n",
    "ax.scatter(ohmc_runtime,ohmc_err,label=\"antithetic OHMC\")\n",
    "\n",
    "ax.legend()\n",
    "plt.xlabel(\"runtime\")\n",
    "plt.ylabel(\"err\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4 Conclusion\n",
    "\n",
    "From the test, one can see OHMC reduce the variance considerably, which make it more appreciable when number of trails is limited. For example, if we can only make trails less than 1000, OHMC is definite better than regular MC. However the runtime efficiency is similar based our observation. It seems slightly stable than regular MC, but still not enough to affirm that. As for the antithetic variates' version, there is no promotion in time efficiency and the variance reduction effect is valid only when number of trajectories is small enough and not gauranteed. Therefore, the antithetic variates method is not effective enough, the delta-based method is neither time efficient nor practical. At last, we highly recommand the optimal hedged Monte Carlo simulation as it is efficient, accurate and compatible to most of the stochastic process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Boundary Scheme Test\n",
    "\n",
    "Five kinds of boundary scheme:\n",
    "\n",
    "* absorption\n",
    "* reflection\n",
    "* Higham and Mao\n",
    "* partial truncation\n",
    "* full truncation\n",
    "\n",
    "We used Euler discretization schemes as presented in Table 1 of the paper \n",
    "\n",
    "Lord, Roger, Remmert Koekkoek, and Dick Van Dijk. A comparison of biased simulation schemes for stochastic volatility models. Quantitative Finance 10.2 (2010): 177-194.\n",
    "\n",
    "https://www.ssoar.info/ssoar/bitstream/handle/document/22127/ssoar-2010-02-lord_et_al-a_comparison_of_biased_simulation.pdf?sequence=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def boundary_scheme_comparison(n_trials,n_steps,boundarySchemes):\n",
    "    simulation_start = timeit.default_timer()\n",
    "    mc = MonteCarlo(S0=stock_price,K=strike,T=time_to_maturity,r=risk_free_rate,q=dividend,sigma=volatility,\n",
    "                kappa=kappa,theta=theta,xi=xi,rho=rho,V0=V0,underlying_process=\"Heston model\")\n",
    "    simulation_end = timeit.default_timer()\n",
    "    theoretical_value = 6.8061\n",
    "    \n",
    "    se_estimate_list = []\n",
    "    rmse_list = []\n",
    "    abs_err_list = []\n",
    "    runtime_list = []\n",
    "    \n",
    "    for boundaryScheme in boundarySchemes:\n",
    "        # regular MC\n",
    "        rmc_start = timeit.default_timer()\n",
    "        mc.simulate(n_trials,n_steps,boundaryScheme=boundaryScheme)\n",
    "        rmc_price = mc.MCPricer('c')\n",
    "        rmc_end = timeit.default_timer()\n",
    "        \n",
    "        rmc_se_estimate = mc.standard_error()\n",
    "        rmc_rmse = np.sqrt(np.mean((mc.value_results - theoretical_value)**2))\n",
    "        rmc_runtime = rmc_end-rmc_start + simulation_end - simulation_start\n",
    "        rmc_err = np.abs(rmc_price - theoretical_value)\n",
    "        se_estimate_list.append(rmc_se_estimate)\n",
    "        rmse_list.append(rmc_rmse)\n",
    "        abs_err_list.append(rmc_err)\n",
    "        runtime_list.append(rmc_runtime)\n",
    "        \n",
    "    result = {\"boundary scheme\": boundarySchemes,\n",
    "              \"n_trials\": [n_trials]*len(boundarySchemes),\n",
    "              \"n_steps\": [n_steps]*len(boundarySchemes),\n",
    "              \"runtime\": runtime_list,\n",
    "              \"abs err\": abs_err_list,\n",
    "              \"root mean square err\": rmse_list,\n",
    "              \"standard err estimate\": se_estimate_list\n",
    "             }\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_trials: 100; n_steps: 50\n",
      "n_trials: 100; n_steps: 100\n",
      "n_trials: 100; n_steps: 150\n",
      "n_trials: 600; n_steps: 50\n",
      "n_trials: 600; n_steps: 100\n",
      "n_trials: 600; n_steps: 150\n"
     ]
    }
   ],
   "source": [
    "n_trials_list = np.arange(100,1000,500)\n",
    "n_steps_list = np.arange(50,200,50)\n",
    "boundarySchemes = [\"absorption\", \"reflection\", \"Higham and Mao\", \"partial truncation\", \"full truncation\"]\n",
    "boundary_scheme_df = pd.DataFrame()\n",
    "for n_trials in n_trials_list:\n",
    "    for n_steps in n_steps_list:\n",
    "        print(\"n_trials: {}; n_steps: {}\".format(n_trials,n_steps))\n",
    "        new_df = pd.DataFrame(boundary_scheme_comparison(n_trials,n_steps,boundarySchemes))\n",
    "        boundary_scheme_df = pd.concat((boundary_scheme_df,new_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>abs err</th>\n",
       "      <th>root mean square err</th>\n",
       "      <th>runtime</th>\n",
       "      <th>standard err estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>boundary scheme</th>\n",
       "      <th>n_steps</th>\n",
       "      <th>n_trials</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">Higham and Mao</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.486407</td>\n",
       "      <td>8.224467</td>\n",
       "      <td>0.004541</td>\n",
       "      <td>0.825143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.315152</td>\n",
       "      <td>7.818877</td>\n",
       "      <td>0.006758</td>\n",
       "      <td>0.319211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.527316</td>\n",
       "      <td>8.275765</td>\n",
       "      <td>0.004912</td>\n",
       "      <td>0.830055</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.429531</td>\n",
       "      <td>8.367403</td>\n",
       "      <td>0.014096</td>\n",
       "      <td>0.341432</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>1.570269</td>\n",
       "      <td>8.525049</td>\n",
       "      <td>0.007029</td>\n",
       "      <td>0.842140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.388945</td>\n",
       "      <td>7.282411</td>\n",
       "      <td>0.021307</td>\n",
       "      <td>0.297127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">absorption</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.756648</td>\n",
       "      <td>8.064365</td>\n",
       "      <td>0.007919</td>\n",
       "      <td>0.806924</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.094931</td>\n",
       "      <td>7.756724</td>\n",
       "      <td>0.012094</td>\n",
       "      <td>0.316907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>1.206100</td>\n",
       "      <td>8.070553</td>\n",
       "      <td>0.006275</td>\n",
       "      <td>0.802012</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.186890</td>\n",
       "      <td>7.652998</td>\n",
       "      <td>0.021262</td>\n",
       "      <td>0.312600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.678147</td>\n",
       "      <td>7.121772</td>\n",
       "      <td>0.009942</td>\n",
       "      <td>0.712513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.580559</td>\n",
       "      <td>8.236398</td>\n",
       "      <td>0.025793</td>\n",
       "      <td>0.335693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">full truncation</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.344331</td>\n",
       "      <td>7.646292</td>\n",
       "      <td>0.005200</td>\n",
       "      <td>0.767702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.617695</td>\n",
       "      <td>8.117524</td>\n",
       "      <td>0.007140</td>\n",
       "      <td>0.330711</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>1.490339</td>\n",
       "      <td>7.788039</td>\n",
       "      <td>0.004902</td>\n",
       "      <td>0.768262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.569383</td>\n",
       "      <td>7.878562</td>\n",
       "      <td>0.014538</td>\n",
       "      <td>0.321068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.089846</td>\n",
       "      <td>7.374268</td>\n",
       "      <td>0.008343</td>\n",
       "      <td>0.741087</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.036406</td>\n",
       "      <td>7.561167</td>\n",
       "      <td>0.021863</td>\n",
       "      <td>0.308937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">partial truncation</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.774601</td>\n",
       "      <td>7.282134</td>\n",
       "      <td>0.005450</td>\n",
       "      <td>0.727730</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.412367</td>\n",
       "      <td>7.321595</td>\n",
       "      <td>0.007170</td>\n",
       "      <td>0.298677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.001891</td>\n",
       "      <td>8.496781</td>\n",
       "      <td>0.004603</td>\n",
       "      <td>0.853959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.046981</td>\n",
       "      <td>7.703554</td>\n",
       "      <td>0.016748</td>\n",
       "      <td>0.314753</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>0.865988</td>\n",
       "      <td>8.340313</td>\n",
       "      <td>0.007377</td>\n",
       "      <td>0.833702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.192895</td>\n",
       "      <td>8.065819</td>\n",
       "      <td>0.021324</td>\n",
       "      <td>0.329466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">reflection</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">50</th>\n",
       "      <th>100</th>\n",
       "      <td>0.040811</td>\n",
       "      <td>7.655841</td>\n",
       "      <td>0.005360</td>\n",
       "      <td>0.769430</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.746866</td>\n",
       "      <td>8.913519</td>\n",
       "      <td>0.009944</td>\n",
       "      <td>0.362916</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">100</th>\n",
       "      <th>100</th>\n",
       "      <td>0.343496</td>\n",
       "      <td>9.255441</td>\n",
       "      <td>0.004634</td>\n",
       "      <td>0.929566</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.100272</td>\n",
       "      <td>7.797641</td>\n",
       "      <td>0.014512</td>\n",
       "      <td>0.318577</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">150</th>\n",
       "      <th>100</th>\n",
       "      <td>1.140397</td>\n",
       "      <td>7.690254</td>\n",
       "      <td>0.007680</td>\n",
       "      <td>0.764354</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>600</th>\n",
       "      <td>0.264661</td>\n",
       "      <td>7.908352</td>\n",
       "      <td>0.020880</td>\n",
       "      <td>0.322946</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      abs err  root mean square err   runtime  \\\n",
       "boundary scheme    n_steps n_trials                                             \n",
       "Higham and Mao     50      100       0.486407              8.224467  0.004541   \n",
       "                           600       0.315152              7.818877  0.006758   \n",
       "                   100     100       0.527316              8.275765  0.004912   \n",
       "                           600       0.429531              8.367403  0.014096   \n",
       "                   150     100       1.570269              8.525049  0.007029   \n",
       "                           600       0.388945              7.282411  0.021307   \n",
       "absorption         50      100       0.756648              8.064365  0.007919   \n",
       "                           600       0.094931              7.756724  0.012094   \n",
       "                   100     100       1.206100              8.070553  0.006275   \n",
       "                           600       0.186890              7.652998  0.021262   \n",
       "                   150     100       0.678147              7.121772  0.009942   \n",
       "                           600       0.580559              8.236398  0.025793   \n",
       "full truncation    50      100       0.344331              7.646292  0.005200   \n",
       "                           600       0.617695              8.117524  0.007140   \n",
       "                   100     100       1.490339              7.788039  0.004902   \n",
       "                           600       0.569383              7.878562  0.014538   \n",
       "                   150     100       0.089846              7.374268  0.008343   \n",
       "                           600       0.036406              7.561167  0.021863   \n",
       "partial truncation 50      100       0.774601              7.282134  0.005450   \n",
       "                           600       0.412367              7.321595  0.007170   \n",
       "                   100     100       0.001891              8.496781  0.004603   \n",
       "                           600       0.046981              7.703554  0.016748   \n",
       "                   150     100       0.865988              8.340313  0.007377   \n",
       "                           600       0.192895              8.065819  0.021324   \n",
       "reflection         50      100       0.040811              7.655841  0.005360   \n",
       "                           600       0.746866              8.913519  0.009944   \n",
       "                   100     100       0.343496              9.255441  0.004634   \n",
       "                           600       0.100272              7.797641  0.014512   \n",
       "                   150     100       1.140397              7.690254  0.007680   \n",
       "                           600       0.264661              7.908352  0.020880   \n",
       "\n",
       "                                     standard err estimate  \n",
       "boundary scheme    n_steps n_trials                         \n",
       "Higham and Mao     50      100                    0.825143  \n",
       "                           600                    0.319211  \n",
       "                   100     100                    0.830055  \n",
       "                           600                    0.341432  \n",
       "                   150     100                    0.842140  \n",
       "                           600                    0.297127  \n",
       "absorption         50      100                    0.806924  \n",
       "                           600                    0.316907  \n",
       "                   100     100                    0.802012  \n",
       "                           600                    0.312600  \n",
       "                   150     100                    0.712513  \n",
       "                           600                    0.335693  \n",
       "full truncation    50      100                    0.767702  \n",
       "                           600                    0.330711  \n",
       "                   100     100                    0.768262  \n",
       "                           600                    0.321068  \n",
       "                   150     100                    0.741087  \n",
       "                           600                    0.308937  \n",
       "partial truncation 50      100                    0.727730  \n",
       "                           600                    0.298677  \n",
       "                   100     100                    0.853959  \n",
       "                           600                    0.314753  \n",
       "                   150     100                    0.833702  \n",
       "                           600                    0.329466  \n",
       "reflection         50      100                    0.769430  \n",
       "                           600                    0.362916  \n",
       "                   100     100                    0.929566  \n",
       "                           600                    0.318577  \n",
       "                   150     100                    0.764354  \n",
       "                           600                    0.322946  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted_boundary_scheme_df = boundary_scheme_df.set_index(['boundary scheme','n_steps','n_trials']).sort_index()\n",
    "sorted_boundary_scheme_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average the above form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>abs err</th>\n",
       "      <th>n_steps</th>\n",
       "      <th>n_trials</th>\n",
       "      <th>root mean square err</th>\n",
       "      <th>runtime</th>\n",
       "      <th>standard err estimate</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>boundary scheme</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Higham and Mao</th>\n",
       "      <td>0.619603</td>\n",
       "      <td>100.0</td>\n",
       "      <td>350.0</td>\n",
       "      <td>8.082329</td>\n",
       "      <td>0.009774</td>\n",
       "      <td>0.575851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>absorption</th>\n",
       "      <td>0.583879</td>\n",
       "      <td>100.0</td>\n",
       "      <td>350.0</td>\n",
       "      <td>7.817135</td>\n",
       "      <td>0.013881</td>\n",
       "      <td>0.547775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>full truncation</th>\n",
       "      <td>0.524667</td>\n",
       "      <td>100.0</td>\n",
       "      <td>350.0</td>\n",
       "      <td>7.727642</td>\n",
       "      <td>0.010331</td>\n",
       "      <td>0.539628</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>partial truncation</th>\n",
       "      <td>0.382454</td>\n",
       "      <td>100.0</td>\n",
       "      <td>350.0</td>\n",
       "      <td>7.868366</td>\n",
       "      <td>0.010445</td>\n",
       "      <td>0.559714</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>reflection</th>\n",
       "      <td>0.439417</td>\n",
       "      <td>100.0</td>\n",
       "      <td>350.0</td>\n",
       "      <td>8.203508</td>\n",
       "      <td>0.010502</td>\n",
       "      <td>0.577965</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     abs err  n_steps  n_trials  root mean square err  \\\n",
       "boundary scheme                                                         \n",
       "Higham and Mao      0.619603    100.0     350.0              8.082329   \n",
       "absorption          0.583879    100.0     350.0              7.817135   \n",
       "full truncation     0.524667    100.0     350.0              7.727642   \n",
       "partial truncation  0.382454    100.0     350.0              7.868366   \n",
       "reflection          0.439417    100.0     350.0              8.203508   \n",
       "\n",
       "                     runtime  standard err estimate  \n",
       "boundary scheme                                      \n",
       "Higham and Mao      0.009774               0.575851  \n",
       "absorption          0.013881               0.547775  \n",
       "full truncation     0.010331               0.539628  \n",
       "partial truncation  0.010445               0.559714  \n",
       "reflection          0.010502               0.577965  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boundary_scheme_df.groupby(\"boundary scheme\").apply(lambda x: x.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the absolute error, Higham and Mao's scheme works the best and the reflection scheme performs the worst."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}