Part 3a: Picking Pairs with Hurst Exponent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "72a0b261",
   "metadata": {},
   "source": [
    "## Use the Hurst exponent to evaluate optimal pairs from a basket of assets\n",
    "* assets in the basket are from sp500\n",
    "* use to asses whether time series speed of diffusion from its initial value is slower than geo random walk\n",
    "* speed of diffusion:\n",
    "    * Var(tau) = <|z(t+tau)-z(t)|^2>\n",
    "        * can be aprox <|z(t + τ) − z(t)|^2> ∼ tau^(2H)\n",
    "        * H is Hurst expo\n",
    "        * Geo random walk has H=0.5\n",
    "        * As H decreases to 0, speed of diffusion reduces so price mean-reverting\n",
    "        * As H increases to 1, price is increasingly trending\n",
    "* use spread to calculate H, rank pairs by H value (low being optimal for mean reversion strat)\n",
    "* as described in (A Machine Learning based Pairs Trading Investment Strategy (Simão Moraes Sarmento, Nuno Horta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "562d1d02",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import yfinance as yf\n",
    "import itertools\n",
    "import math\n",
    "from co import *\n",
    "import ndtest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfda1df1",
   "metadata": {},
   "source": [
    "## Calculate the Hurst exponent of a time series\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4bde4e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_hurst_exponent(time_series, max_lag=20):\n",
    "    lags = range(2, max_lag)\n",
    "    \n",
    "    # variances of the lagged differences\n",
    "    tau = [np.std(np.subtract(time_series[lag:], time_series[:-lag])) for lag in lags]\n",
    "\n",
    "    # calculate the slope of the log plot -> the Hurst Exponent\n",
    "    reg = np.polyfit(np.log(lags), np.log(tau), 1)\n",
    "\n",
    "    return reg[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c5bb6fd",
   "metadata": {},
   "source": [
    "## Asset list and dates or calibration and trading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d73a4561",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Asset list\n",
    "tickers = pd.read_csv('constituents_csv.csv', index_col=1)\n",
    "asset_list=tickers['Symbol'].tolist()\n",
    "lst=asset_list[0:50]\n",
    "\n",
    "#dates for model calibration\n",
    "start_data='2019-08-19'\n",
    "end_data='2022-08-19'\n",
    "\n",
    "#dates for trading with algo\n",
    "start_trade='2022-01-01'\n",
    "end_trade='2022-06-01'\n",
    "\n",
    "index_comp='SPY'\n",
    "\n",
    "#number of pairs traded\n",
    "n=10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10785b5d",
   "metadata": {},
   "source": [
    "## Calculate Hurst exponent for the spread of every possible pair combo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed46cc0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "pair_order_list = itertools.combinations(lst,2)\n",
    "pairs=list(pair_order_list)\n",
    "pairs[0][1]\n",
    "lst=[]\n",
    "for i in range(0,len(pairs)):\n",
    "    asset1=pairs[i][0]\n",
    "    asset2=pairs[i][1]\n",
    "    \n",
    "    #adjusted close and spread of asset pairs\n",
    "    X = yf.download(asset1, start = start_data, end=end_data)\n",
    "    X.reset_index(inplace=True)\n",
    "    np.log(X['Adj Close']) - np.log(X['Adj Close'].shift(1))\n",
    "    date=X['Date']\n",
    "    date=date.to_frame()\n",
    "\n",
    "    X=X['Adj Close']\n",
    "    X_df=X.to_frame()\n",
    "    X_df=X_df.rename(columns={'Adj Close': asset1}).dropna()\n",
    "\n",
    "    Y = yf.download(asset2, start = start_data, end=end_data)\n",
    "    Y.reset_index(inplace=True)\n",
    "    np.log(Y['Adj Close']) - np.log(Y['Adj Close'].shift(1))\n",
    "\n",
    "    Y=Y['Adj Close']\n",
    "    Y_df=Y.to_frame()\n",
    "    Y_df=Y_df.rename(columns={'Adj Close': asset2}).dropna()\n",
    "\n",
    "    pairs_df=X_df.join(Y_df).join(date).set_index('Date').dropna()\n",
    "\n",
    "    #calculate spread of adjusted close price between asset pair\n",
    "    pairs_df['spread'] = pairs_df[asset1] - pairs_df[asset2]\n",
    "    \n",
    "    hurst_exp = get_hurst_exponent(pairs_df['spread'].values, 500)\n",
    "    lst.append(tuple((pairs[i],hurst_exp)))\n",
    "    \n",
    "lst.sort(key=lambda a: a[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "4631da28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hurst exponent for all possible pairs:\n"
     ]
    },
    {
     "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>Hurst expo</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ADM-ABC</th>\n",
       "      <td>0.030749</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AOS-ADI</th>\n",
       "      <td>0.052112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>APD-AMZN</th>\n",
       "      <td>0.053331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ARE-GOOG</th>\n",
       "      <td>0.058664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ARE-GOOGL</th>\n",
       "      <td>0.058985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AFL-AMZN</th>\n",
       "      <td>0.472588</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ALGN-ANTM</th>\n",
       "      <td>0.472766</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AMZN-AIG</th>\n",
       "      <td>0.473629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AMZN-APA</th>\n",
       "      <td>0.481570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AMZN-ANET</th>\n",
       "      <td>0.511900</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1225 rows × 1 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           Hurst expo\n",
       "ADM-ABC      0.030749\n",
       "AOS-ADI      0.052112\n",
       "APD-AMZN     0.053331\n",
       "ARE-GOOG     0.058664\n",
       "ARE-GOOGL    0.058985\n",
       "...               ...\n",
       "AFL-AMZN     0.472588\n",
       "ALGN-ANTM    0.472766\n",
       "AMZN-AIG     0.473629\n",
       "AMZN-APA     0.481570\n",
       "AMZN-ANET    0.511900\n",
       "\n",
       "[1225 rows x 1 columns]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = list(zip(*lst))\n",
    "pairs_lst=[]\n",
    "asset1_lst=[]\n",
    "asset2_lst=[]\n",
    "h_expo=[]\n",
    "possible_pairs=[]\n",
    "for i in range(0,len(lst)): \n",
    "    asset1_lst.append(res[0][i][0])\n",
    "    asset2_lst.append(res[0][i][1])\n",
    "    h_expo.append(res[1][i])\n",
    "    possible_pairs.append(asset1_lst[i]+'-'+asset2_lst[i])\n",
    "hurst_df=pd.DataFrame(possible_pairs,h_expo)\n",
    "hurst_df = pd.DataFrame(h_expo, index =possible_pairs,columns =['Hurst expo'])\n",
    "print(\"Hurst exponent for all possible pairs:\")\n",
    "\n",
    "hurst_df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9371ac9e",
   "metadata": {},
   "source": [
    "## Find optimal unique pairs with lowest Hurst exponent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91e3f05f",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_stocks = []\n",
    "selected_pairs = []\n",
    "opt_asset1=[]\n",
    "opt_asset2=[]\n",
    "for i in range(0,len(asset1_lst)):\n",
    "    s1=asset1_lst[i]\n",
    "    s2=asset2_lst[i]\n",
    "  \n",
    "    if (s1 not in selected_stocks) and (s2 not in selected_stocks):\n",
    "        selected_stocks.append(s1)\n",
    "        selected_stocks.append(s2)\n",
    "        pair=s1+' and '+s2\n",
    "        selected_pairs.append(pair)\n",
    "    \n",
    "    if len(selected_pairs) == math.comb(len(lst),2):\n",
    "        break\n",
    "        \n",
    "\n",
    "print('Optimal Pairs with lowest Hurst expo:',selected_pairs)\n",
    "opt_asset1=selected_stocks[0:len(selected_stocks)-1:2]\n",
    "opt_asset2=selected_stocks[1:len(selected_stocks):2]\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c00be043",
   "metadata": {},
   "source": [
    "## Fit copulas to each optimal pair using AIC and BIC and test if good fit with K-S test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2654e8de",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=pd.DataFrame()\n",
    "best_cop_lst=[]\n",
    "for i in range(0,len(opt_asset1)):\n",
    "    #assets (yfinance ticker)\n",
    "    asset1=opt_asset1[i]\n",
    "    asset2=opt_asset2[i]\n",
    "    pairs=selected_pairs[i] #insert name of asset 1 and asset 2\n",
    "\n",
    "    #fetch log returns of pair (to model calibration)\n",
    "\n",
    "    X = yf.download(asset1, start = start_data, end=end_data)\n",
    "    X.reset_index(inplace=True)\n",
    "    X.tail()\n",
    "    X['Log_Returns'] = np.log(X['Adj Close']) - np.log(X['Adj Close'].shift(1))\n",
    "    x = X['Log_Returns'].drop(labels=0, axis=0).dropna()\n",
    "    date=X['Date'].drop(labels=0,axis=0)\n",
    "    date=date.to_frame()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    Y = yf.download(asset2, start = start_data, end=end_data)\n",
    "    Y.reset_index(inplace=True)\n",
    "    Y.tail()\n",
    "    Y['Log_Returns'] = np.log(Y['Adj Close']) - np.log(Y['Adj Close'].shift(1))\n",
    "    y = Y['Log_Returns'].drop(labels=0, axis=0).dropna()\n",
    "\n",
    "    #fit distribution to element in pair + probability integral transform to get marginal uniform\n",
    "\n",
    "    param1 = stats.t.fit(x)\n",
    "    dist1=stats.t(*param1)\n",
    "    u=dist1.cdf(x)\n",
    "    param2=stats.t.fit(y)\n",
    "    dist2=stats.t(*param2)\n",
    "    v=dist2.cdf(y)\n",
    "\n",
    "\n",
    "\n",
    "    #put returns in a data frame with corresponding dates\n",
    "    x_df=x.to_frame()\n",
    "    y_df=y.to_frame()\n",
    "\n",
    "    a=y_df.rename(columns={'Log_Returns': asset2})\n",
    "    b=x_df.rename(columns={'Log_Returns': asset1})\n",
    "\n",
    "    p =a.join(b)\n",
    "    x=p.join(date)\n",
    "\n",
    "    returns_df=x.set_index('Date')\n",
    "\n",
    "\n",
    "    #fit bivariate copula family to data \n",
    "    #calculate Akaike information criterion, Bayesion information criterion and KS value to help determine good fit\n",
    "\n",
    "    pair=[pairs]\n",
    "    copulas_df = pd.DataFrame(index=pair, columns=['copula', 'parameter', 'aic', 'bic', 'KS_pvalue'])\n",
    "\n",
    "\n",
    "    best_aic = np.inf\n",
    "    best_copula= None\n",
    "\n",
    "\n",
    "    copulas=[GaussianCopula(),GumbelCopula(), ClaytonCopula(),FrankCopula()]\n",
    "    for copula in copulas:\n",
    "        copula.fit(u,v)\n",
    "\n",
    "        L = copula.log_likelihood(u,v)\n",
    "        aic = 2 * copula.num_params - 2 * L\n",
    "        if aic < best_aic:\n",
    "            best_aic = aic\n",
    "            best_bic = copula.num_params * np.log(len(u)) - 2 * L\n",
    "            best_copula = copula\n",
    "\n",
    "            # calculate KS-pvalue\n",
    "            smp = copula.sample(size=len(u)) # generate sample from fit copula\n",
    "            s_u = smp[:,0]\n",
    "            s_v = smp[:,1]\n",
    "            ks_pval = ndtest.ks2d2s(u,v,s_u,s_v)\n",
    "\n",
    "            if isinstance(copula, ArchimedeanCopula):\n",
    "                best_param = copula.alpha\n",
    "            else:\n",
    "                best_param = copula.rho\n",
    "    \n",
    "    best_cop_lst.append(best_copula)\n",
    "    copulas_df.loc[pair] = [best_copula.name, best_param, best_aic, best_bic, ks_pval]\n",
    "    df=df.append(copulas_df)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6ca66feb",
   "metadata": {},
   "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>copula</th>\n",
       "      <th>parameter</th>\n",
       "      <th>aic</th>\n",
       "      <th>bic</th>\n",
       "      <th>KS_pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ADM and ABC</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.425854</td>\n",
       "      <td>-451.892147</td>\n",
       "      <td>-447.264105</td>\n",
       "      <td>0.566356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AOS and ADI</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.505171</td>\n",
       "      <td>-502.235208</td>\n",
       "      <td>-497.607166</td>\n",
       "      <td>0.186951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>APD and AMZN</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.379266</td>\n",
       "      <td>-420.098778</td>\n",
       "      <td>-415.470736</td>\n",
       "      <td>0.282734</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ARE and GOOG</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.447343</td>\n",
       "      <td>-465.743336</td>\n",
       "      <td>-461.115295</td>\n",
       "      <td>0.222887</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAPL and AMAT</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.60703</td>\n",
       "      <td>-573.252947</td>\n",
       "      <td>-568.624906</td>\n",
       "      <td>0.207337</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AON and ANET</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.407004</td>\n",
       "      <td>-439.371301</td>\n",
       "      <td>-434.74326</td>\n",
       "      <td>0.259305</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAP and ALB</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.354038</td>\n",
       "      <td>-401.771051</td>\n",
       "      <td>-397.14301</td>\n",
       "      <td>0.323573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ABT and AWK</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.512812</td>\n",
       "      <td>-507.073452</td>\n",
       "      <td>-502.445411</td>\n",
       "      <td>0.226867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MO and AEE</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.343429</td>\n",
       "      <td>-393.5404</td>\n",
       "      <td>-388.912359</td>\n",
       "      <td>0.260405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MMM and A</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.428641</td>\n",
       "      <td>-453.716618</td>\n",
       "      <td>-449.088576</td>\n",
       "      <td>0.514289</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAL and AJG</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.340555</td>\n",
       "      <td>-391.294957</td>\n",
       "      <td>-386.666916</td>\n",
       "      <td>0.152344</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ABMD and ANSS</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.562151</td>\n",
       "      <td>-539.679757</td>\n",
       "      <td>-535.051715</td>\n",
       "      <td>0.239223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ABBV and AEP</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.185299</td>\n",
       "      <td>-244.863478</td>\n",
       "      <td>-240.235436</td>\n",
       "      <td>0.257857</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AKAM and ALLE</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.379113</td>\n",
       "      <td>-419.994158</td>\n",
       "      <td>-415.366117</td>\n",
       "      <td>0.434909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AMCR and AMT</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.338721</td>\n",
       "      <td>-390.01003</td>\n",
       "      <td>-385.381988</td>\n",
       "      <td>0.273636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LNT and ALL</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.421197</td>\n",
       "      <td>-448.803086</td>\n",
       "      <td>-444.175045</td>\n",
       "      <td>0.291062</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AMD and APTV</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.428866</td>\n",
       "      <td>-453.855946</td>\n",
       "      <td>-449.227905</td>\n",
       "      <td>0.194612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACN and AMP</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.63078</td>\n",
       "      <td>-593.403364</td>\n",
       "      <td>-588.775322</td>\n",
       "      <td>0.649818</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GOOGL and AME</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.525197</td>\n",
       "      <td>-514.873319</td>\n",
       "      <td>-510.245277</td>\n",
       "      <td>0.557398</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AES and AMGN</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.245192</td>\n",
       "      <td>-307.79248</td>\n",
       "      <td>-303.164438</td>\n",
       "      <td>0.823554</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AXP and AIG</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.742064</td>\n",
       "      <td>-732.813147</td>\n",
       "      <td>-728.185106</td>\n",
       "      <td>0.279836</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AFL and ALK</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.608809</td>\n",
       "      <td>-574.69476</td>\n",
       "      <td>-570.066719</td>\n",
       "      <td>0.227348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ANTM and APA</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.271938</td>\n",
       "      <td>-333.357909</td>\n",
       "      <td>-328.729868</td>\n",
       "      <td>0.668182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ADBE and ALGN</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.541779</td>\n",
       "      <td>-525.829389</td>\n",
       "      <td>-521.201348</td>\n",
       "      <td>0.499007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ATVI and APH</th>\n",
       "      <td>Gaussian</td>\n",
       "      <td>0.313396</td>\n",
       "      <td>-369.54967</td>\n",
       "      <td>-364.921629</td>\n",
       "      <td>0.411017</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 copula parameter         aic         bic KS_pvalue\n",
       "ADM and ABC    Gaussian  0.425854 -451.892147 -447.264105  0.566356\n",
       "AOS and ADI    Gaussian  0.505171 -502.235208 -497.607166  0.186951\n",
       "APD and AMZN   Gaussian  0.379266 -420.098778 -415.470736  0.282734\n",
       "ARE and GOOG   Gaussian  0.447343 -465.743336 -461.115295  0.222887\n",
       "AAPL and AMAT  Gaussian   0.60703 -573.252947 -568.624906  0.207337\n",
       "AON and ANET   Gaussian  0.407004 -439.371301  -434.74326  0.259305\n",
       "AAP and ALB    Gaussian  0.354038 -401.771051  -397.14301  0.323573\n",
       "ABT and AWK    Gaussian  0.512812 -507.073452 -502.445411  0.226867\n",
       "MO and AEE     Gaussian  0.343429   -393.5404 -388.912359  0.260405\n",
       "MMM and A      Gaussian  0.428641 -453.716618 -449.088576  0.514289\n",
       "AAL and AJG    Gaussian  0.340555 -391.294957 -386.666916  0.152344\n",
       "ABMD and ANSS  Gaussian  0.562151 -539.679757 -535.051715  0.239223\n",
       "ABBV and AEP   Gaussian  0.185299 -244.863478 -240.235436  0.257857\n",
       "AKAM and ALLE  Gaussian  0.379113 -419.994158 -415.366117  0.434909\n",
       "AMCR and AMT   Gaussian  0.338721  -390.01003 -385.381988  0.273636\n",
       "LNT and ALL    Gaussian  0.421197 -448.803086 -444.175045  0.291062\n",
       "AMD and APTV   Gaussian  0.428866 -453.855946 -449.227905  0.194612\n",
       "ACN and AMP    Gaussian   0.63078 -593.403364 -588.775322  0.649818\n",
       "GOOGL and AME  Gaussian  0.525197 -514.873319 -510.245277  0.557398\n",
       "AES and AMGN   Gaussian  0.245192  -307.79248 -303.164438  0.823554\n",
       "AXP and AIG    Gaussian  0.742064 -732.813147 -728.185106  0.279836\n",
       "AFL and ALK    Gaussian  0.608809  -574.69476 -570.066719  0.227348\n",
       "ANTM and APA   Gaussian  0.271938 -333.357909 -328.729868  0.668182\n",
       "ADBE and ALGN  Gaussian  0.541779 -525.829389 -521.201348  0.499007\n",
       "ATVI and APH   Gaussian  0.313396  -369.54967 -364.921629  0.411017"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d177ecb",
   "metadata": {},
   "source": [
    "## Algo returns for specific trading period\n",
    "*refer to (Pairs Trading with Copulas.ipynb) for further explanation of algo/\n",
    "* select n pairs from optimal pairs with the lowest Hurst expo\n",
    "##### specific trades can be seen by printing out the positions DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4df878c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dataframe of the returns for the specific trading period for each pair\n",
    "algo_returns = {}\n",
    "cl = 0.99\n",
    "\n",
    "for q in range(0,n):\n",
    "    X_trade = yf.download(opt_asset1[q], start =start_trade, end=end_trade)\n",
    "    X_trade.reset_index(inplace=True)\n",
    "    X_trade['Log_Returns'] = np.log(X_trade['Adj Close']) - np.log(X_trade['Adj Close'].shift(1))\n",
    "    x_trade = X_trade['Log_Returns'].drop(labels=0, axis=0).dropna()\n",
    "    date=X_trade['Date'].drop(labels=0,axis=0)\n",
    "    date=date.to_frame()\n",
    "    \n",
    "\n",
    "    Y_trade = yf.download(opt_asset2[q], start = start_trade, end=end_trade)\n",
    "    Y_trade.reset_index(inplace=True)\n",
    "    Y_trade.tail()\n",
    "    Y_trade['Log_Returns'] = np.log(Y_trade['Adj Close']) - np.log(Y_trade['Adj Close'].shift(1))\n",
    "    y_trade = Y_trade['Log_Returns'].drop(labels=0, axis=0).dropna()\n",
    "\n",
    "    #put returns in a data frame with corresponding dates\n",
    "    x2_df=x_trade.to_frame()\n",
    "    y2_df=y_trade.to_frame()\n",
    "    \n",
    "    \n",
    "    a=x2_df.rename(columns={'Log_Returns': opt_asset1[q]})\n",
    "    m=y2_df.rename(columns={'Log_Returns': opt_asset2[q]})\n",
    "    \n",
    "\n",
    "    p =m.join(a)\n",
    "    x=p.join(date)\n",
    "\n",
    "    returns_trade_period=x.set_index('Date')\n",
    "    \n",
    "    # calculate conditional probabilities\n",
    "    prob1=[]\n",
    "    prob2=[]\n",
    "\n",
    "    for u,v in zip(dist1.cdf(returns_trade_period[opt_asset1[q]]), dist2.cdf(returns_trade_period[opt_asset2[q]])):\n",
    "        prob1.append(best_cop_lst[q].cdf_u_given_v(u,v))\n",
    "        prob2.append(best_cop_lst[q].cdf_v_given_u(u,v))\n",
    "\n",
    "\n",
    "    probs_trade= pd.DataFrame(np.vstack([prob1,prob2]).T,index=returns_trade_period.index, columns=[opt_asset1[q],opt_asset2[q]])\n",
    "\n",
    "\n",
    "\n",
    "    #calculate positions\n",
    "    positions=pd.DataFrame(index=probs_trade.index, columns = probs_trade.columns)\n",
    "    long = False\n",
    "    short = False\n",
    "\n",
    "\n",
    "    for t in positions.index:    \n",
    "        # if long position is open\n",
    "        if long:\n",
    "            if (probs_trade.loc[t][opt_asset1[q]] > 0.5) or (probs_trade.loc[t][opt_asset2[q]] < 0.5):\n",
    "                positions.loc[t] = [0,0]\n",
    "                long = False\n",
    "            else:\n",
    "                positions.loc[t] = [1,-1]\n",
    "\n",
    "        # if short position is open\n",
    "        elif short:\n",
    "            if (probs_trade.loc[t][opt_asset1[q]] < 0.5) or (probs_trade.loc[t][opt_asset2[q]] > 0.5):\n",
    "                positions.loc[t] = [0,0]\n",
    "                short = False\n",
    "            else:\n",
    "                positions.loc[t] = [-1,1]\n",
    "\n",
    "        # if no positions are open\n",
    "        else:\n",
    "            if (probs_trade.loc[t][opt_asset1[q]] < (1-cl)) and (probs_trade.loc[t][opt_asset2[q]] > cl):\n",
    "                # open long position\n",
    "                positions.loc[t] = [1,-1]\n",
    "                long = True\n",
    "            elif (probs_trade.loc[t][opt_asset1[q]] > cl) and (probs_trade.loc[t][opt_asset1[q]] < (1-cl)):\n",
    "                # open short positions\n",
    "                positions.loc[t] = [-1,1]\n",
    "                short = True\n",
    "            else:\n",
    "                positions.loc[t] = [0,0]\n",
    "\n",
    "\n",
    "    # calculate returns\n",
    "    algo_ret = (returns_trade_period * positions.shift()).sum(axis=1)\n",
    "    algo_returns[selected_pairs[q]] = algo_ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6305e8c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "rets = pd.DataFrame.from_dict(algo_returns)\n",
    "rets = np.exp(rets) - 1 # convert log-returns to simple returns\n",
    "total_ret = rets.sum(axis=1) / len(rets.columns) * 2 # double capital (from short positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bc05f24",
   "metadata": {},
   "outputs": [],
   "source": [
    "index_price = yf.download(index_comp, start=start_trade, end=end_trade)\n",
    "index_price = index_price['Adj Close']\n",
    "index_ret = index_price.pct_change().dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a08e1b7d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1071eb820>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAAFlCAYAAABMeyuWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB6nElEQVR4nO3dd3iUZdbH8e+dHnqvSeg99NAUaaIUFRUsIHZX17Luu666666uuu661V1d14qKvay9IIKCgKj03kOHUEPvpN3vH/dEAgQySWbmmUx+n+vK9STPPOVE4zhz5tznGGstIiIiIiIiIiLBEOV1ACIiIiIiIiISuZR4EBEREREREZGgUeJBRERERERERIJGiQcRERERERERCRolHkREREREREQkaJR4EBEREREREZGgifE6gOKoVauWbdy4sddhiIiIiIiIiEgB8+bN22WtrV3YY2Uq8dC4cWPmzp3rdRgiIiIiIiIiUoAxZuOZHtNSCxEREREREREJGiUeRERERERERCRolHgQERERERERkaBR4kFEREREREREgkaJBxEREREREREJGiUeRERERERERCRolHgQERERERERkaBR4kFEREREREREgkaJBxEREREREREJGiUeRERERERERCRolHgQERERERERkaCJ8ToAERERERGR8mTrvqOs2nEwZPdrVKMCTWtXCtn9RE6lxIOIiIiIiEgI3fTqnJAmHoyBYR0bcM/AljSuVTFk9xXJp8SDiIiIiIhIiKzZeZBVOw5yV/9mDGxTN+j3s8A3y3fw2g8bGLd4G1d2TeLu81vQsFpi0O8tkk+JBxERERERkRCZuGwHANf1bEy9qgkhuWeXlOrcdG5jnpuylndmbeLj+Vu4pkcKd/VvTu3K8SGJQco3NZcUEREREREJkQlLt9M5pVrIkg756lRO4NFh7Zhyfz+Gd2nImzM30ucfU/j7hJXsO5IV0lik/FHFg4iIiIiISAhk7D3Cki37+d2Q1p7F0LBaIn8b0YGf923GU5PSeWHaWt6asZEh7esRF6PPpcPFyG4ppDas6nUYAaPEg4iIiIiISAjkL7MY1K6ex5FAk1oV+c/IztzRrxlPfpPO5BU7vQ5JChjQug6gxIOIiIiIiIgUw8Sl22ldr3JYTZZoXa8KL16X5nUYEuFUSyMiIiIiIhJkmQePM2fjHganel/tIBJqSjyIiIiIiIgE2dfLt2MtSjxIuaTEg4iIiIiISJBNWLqdxjUr0KpuZa9DEQk5JR5ERERERESCaP+RbGas3c2g1HoYY7wORyTklHgQEREREREJoskrd5CTZxkcBtMsRLygxIOIiIiIiEgQTVi6nfpVE+iYVM3rUEQ8ocSDiIiIiIhIkBzJymFaeiaD2tUjKkrLLKR88ivxYIwZa4zZaYxZeobHjTHmaWPMGmPMYmNMlwKP3WCMWe37uqHA/q7GmCW+c542WuwkIiIiIiIRZtqqTI7n5DFIyyykHIvx87jXgGeAN87w+BCghe+rB/A80MMYUwN4BEgDLDDPGPO5tXav75hbgVnAeGAw8FXJfg0ROZu8PEtWbp7XYZwkJsoQE62iKxEREYlsE5Ztp0bFOLo1ru51KCKe8SvxYK39zhjT+CyHXAq8Ya21wExjTDVjTH2gH/CNtXYPgDHmG2CwMWYqUMVaO9O3/w3gMiIs8ZC+4yCfTPqevXH1sUZvsCR4svPyOHI8l8NZORzJyuXwcbc9kpXD4eO5HM3O9TrE01RNjOX1m7vTKbma16GIiIiIBMXxnFy+XbGToe3r6wMXKdf8rXgoSkNgc4GfM3z7zrY/o5D9pzHG3AbcBpCSkhKgcENkxTh+u/p2bon5G8uiWnodjUSw6ChDxfhoKsTFUDE+mhoVK1ApPoYKcdFU9G3jYqIwhM+KprdnbeRnr8/ls1+cS8NqiV6HIyIiIhJwP67dzcHjOQxO1TILKd8ClXgIGmvtGGAMQFpamvU4nGJp2X0QTIvmlZ47YeAdXocjElbOb1OHEc/9yC2vzeHDO86hUnzYPx2JiIiIFMvEpdupFB/DOc1reh2KiKcCVe+zBUgu8HOSb9/Z9icVsj+yVKgBjc6BVRG1gkQkIFrWrcwzo7uQvuMgv3pvAbl5ZSqvKCIiInJWuXmWr5fvYEDrOsTHRHsdjoinApV4+By43jfdoiew31q7DZgIXGiMqW6MqQ5cCEz0PXbAGNPTN83ieuCzAMUSXloNhcwVsGed15GIhJ2+LWvz6LB2TFqxk79PWOl1OCIiIiIBM2fDHvYcztIyCxH8H6f5LjADaGWMyTDG3GKMud0Yc7vvkPHAOmAN8BJwJ4CvqeSfgDm+r8fyG036jnnZd85aIqyx5E9aD3VbVT2IFOr6Xo25oVcjxny3jvdmb/I6HBEREZGAmLB0O/ExUfRtWdvrUEQ85+9Ui1FFPG6Bu87w2FhgbCH75wKp/ty/TKveGOq0c4mHXoX+IxIp9/5wcVvW7z7CQ58uJaVmBc5pVsvrkERERMRj1lqWbzvA3sPZXodykrYNqlCjYtxZj7HWMnHZdvq0rE1F9bESCf/mkhGh1RD4/kk4ssf1fRCRk8RER/HMNZ0Z8dyP3PHWfD658xya1q7kdVgiIiLiAWst09IzeXbKGuZs2Ot1OKepHB/D3ec354ZzGp+xd8PijP1s23+M+y5sFeLoRMKTEg+h0HooTH8CVn8DHa/2OhqRsFQlIZZXbujGZc/9wC2vz+WTO8+hWoWzf5ogIiIikSMvz/LNih088+0almzZT/2qCTx6SVvaNazqdWg/OZ6dx9gf1vOX8St5e9YmHhzahgva1sW1rTvhq6XbiYkynN+mjkeRioQXJR5CoX5nqFQPVn2pxIPIWaTUrMCL13Vl9EuzuOOt+bx+c3fiYgLVA1dERETCUW6eZdzirTw7ZQ3pOw7RqGYF/j6iPZd3TgrL1wG9W9RiWnomfxq3nNvenMe5zWvyh4vb0rpeFcBVbExYuo1ezWrqQxQRHyUeQiEqyi23WPIB5ByHmHivIxIJW90a1+BvI9rz6/cX8YdPl/K3Ee1P+xShOKy1rNt1mKNZuQGMUqRkjIH4mGjiY6JIiI0mPjaKhJhoYqNNqf7ORUTKoqycPD5dsIXnp61l/a7DtKhTif+M7MRF7esTEx1+CYeC+raszTn/dx5vz9zIk5NWM/Q/07mmRwq/vqAVmQePs2H3EX52XlOvwxQJG0o8hEqroTDvVVg/HVoM9DoakbA2vEsS6zIP88yUNcTHRvHoJe2Iiir+m7K8PMufv1zB2B/WByFKkcBxCQmXjIiNjkIpCBEpD45m5XLweA6pDavwwrVduLBtvRL9/94rsdFR3HhuEy7t1JCnJqXz1qxNfLZwK23rV8EYuLBtXa9DFAkbSjyESpM+EFsRVo1X4kHED/de2JKs3DzGfLeO/UezeeLKjsQW49OP7Nw8fvPhYj5ZsIVre6bQt6XWWIr3cvMsx3NyOZ6T576y3ffHfNvj2blk5eZ5HaaISEhEGcPANnXp16p2ma76ql4xjj9emsrono3407jlTF+9i7RG1alTJcHr0ETChhIPoRKbAM0HuLGaF/3LfbwlImdkjOF3Q1pTNTGWf05cxcFjOTw3ugsJsYV3jy7oaFYud70zn29X7uS+C1tyV//mZfoFjYiIiIS/lnUr88bN3Zm1fg/1qyrpIFJQeC+eijStLoKDW2HbQq8jESkTjDHc1b85f74slSmrdnL92NkcOHb2Wd77j2Rz7SuzmLJqJ49fnsovBrRQ0kFERERCwhhDz6Y1aVSzotehiIQVJR5CqcWFYKJg5XivIxEpU67t2Yj/jOzM/I17uealmew+dLzQ43YcOMZVL85gScZ+nr2mC6N7NApxpCIiIiIiciolHkKpYk1I6eWWW4hIsQzr2ICXrk9j9Y5DXPniDLbsO3rS4+t3HWbE8z+SsfcIr97UjaHt63sUqYiIiIiIFKTEQ6i1GgI7lsDejV5HIlLm9G9dhzdv6UHmgeNc+fyPrM08BMDSLfu58oUfOZKVy7u39eTc5rU8jlRERERERPIp8RBqrYa6bfoEb+MQKaO6N6nBu7f15HhOHle9MIO3Zm5k5JiZxMdE88HtveiQVM3rEEVEREREpAAlHkKtZjOo1QpWful1JCJlVmrDqnxwey8SYqN56NOl1K+awId39KJZ7UpehyYiIiIiIqfQOE0vtB4KP/4Xju6DxGpeRyNSJjWtXYkPbu/F/+Zs5sZzGlO9YpzXIYmIiIiISCFU8eCFVkMhLwfWTPI6kpLLOQ7Weh2FlHMNqiVyzwUtlXQQEREREQljSjx4oWEaVKwNq8rgWM0dy+CTO+AvDWHm815HIyIiIiIiImFOiQcvREVBy8Gw+hvIyfI6mqJZC2u/hTcvh+fPgeWfQeV6brlIbrbX0YmIiIiIiEgYU+LBK60vguMHYOMPXkdyZjlZsPBdeKG3SzrsWAbnPwz3LIWhT8DBrbDic6+jFBERERERkTCm5pJeadIXYhLdcotm/b2O5mRH98G8V2HWi3BwG9RuA5c+B+2vgJh4d0yLC6FGU5j5AqSO8DRcERERERERCV+qePBKXAVoNgBWfRVeTRr3bYKnO8OkR6FWSxj9Edw5AzqPPpF0ALdcpPvPIWM2bJnnWbgiIiIiIiIS3pR48FKrIbB/M2xf4nUkjrXw+S8hNwtu/RZu+BxaDARjCj++0zUQV9lVPYiIiIiIiIgUQokHL7UcDBhX9RAO5r8B66bABX+Ehl2LPj6hCnS+FpZ9Age3Bz8+ERERERERKXOUePBSpdqQ3ANWfel1JLBvM0x8EBqfB11v9v+87rdCXg7MeSV4sYmIiIiIiEiZpcSD11oNgW2LYNo/YcU4yEwP/YhKa+GL/wObB5c+4/o3+KtmM2g5COaOhexjwYtRREREREREyiRNtfBa6nCY+wpM+fOJfVExUL2Ja+5Yq4Xb1mkNDbqcud9CaSx4C9ZOdiMyqzcu/vk9bof0y2DpR64JpYiIiIiIiIiPEg9eq5YCv1oCxw7A7jWwazXsSvd9rYbVX0OerwKiaX8Y9rQ7J1D2b3FLLBr1hrRbSnaNpv3cyM1Zz7uGk8FIjoiIiIiIiEiZpMRDuEioAg27uK+CcnNg30ZY/Q18+yd4rhcMfNQlCYqzJKIw1sK4X7nExqX/Lfn1jIEeP3fX2vgjND63dHGJiIiIiIhIxFCPh3AXHeP6KPS8He6cAcndYfx98PolsHtt6a696F1XUXH+I1Cjaemu1eFqSKzuqh5EREREREREfJR4KEuqpcC1H8OwZ2D7Enj+XJjxLOTlFv9aB7bCVw9AyjnQ/bbSxxZXAbrcACu/hL0bS389ERERERERiQh+JR6MMYONMauMMWuMMQ8U8ngjY8xkY8xiY8xUY0ySb39/Y8zCAl/HjDGX+R57zRizvsBjnQL5i0UsY6DLdXDXTGjaFyb+HsYOdtMw/GUtjLsHcrOKP8XibLrfChiY81JgriciIiIiIiJlXpHvOI0x0cCzwBCgLTDKGNP2lMOeAN6w1nYAHgP+CmCtnWKt7WSt7QQMAI4AXxc47/78x621C0v7y5QrVRrAqPdg+EuwezW80Bum/wsOZRZ97uL3IX0CnP+wW8YRKFWToM0lMP8NyDocuOuKiIiIiIhImeXPR93dgTXW2nXW2izgPeDSU45pC3zr+35KIY8DXAF8Za09UtJg5RTGQIer4M5Z0PJCmPwYPNEcnkyF/10H3z8J66bBsf0nzjm4Hb76DST3dA0hA63nHe5+i94N/LVFRERERESkzPFnqkVDYHOBnzOAHqccswgYDvwHuByobIypaa3dXeCYkcC/TznvcWPMw8Bk4AFr7fHiBC8+levC1W/B5tmweRZsmQ9b58OKz08cU7M5NOgC+zdDzjG49FmIig58LMk9oH4nmPUidL05cMs4REREREREpEwK1DjN+4BnjDE3At8BW4CfOh4aY+oD7YGJBc75HbAdiAPGAL/FLdM4iTHmNuA2gJSUlACFG6GSu7uvfEf2uATE1gWwZQFsmA4Ht8Ggv0Kt5sGJwRhX9fDJz2Hdt9B8YHDuIyIiIiIiImWCP4mHLUBygZ+TfPt+Yq3diqt4wBhTCRhhrd1X4JCrgE+stdkFztnm+/a4MeZVXPLiNNbaMbjEBGlpadaPeCVfhRrujX/BN/9H90FiteDet93l8PUfXNWDEg8iIiIiIiLlmj+JhzlAC2NME1zCYSRwTcEDjDG1gD3W2jxcJcPYU64xyre/4Dn1rbXbjDEGuAxYWqLfQIon2EkHgJh4SLsZpv0NFv3P7Tu2zyU9ju2Do3tPfH/sAPS5F1JHBD8uERERERERCbkiEw/W2hxjzC9wyySigbHW2mXGmMeAudbaz4F+wF+NMRa31OKu/PONMY1xFRPTTrn028aY2oABFgK3l/q3kfCRdjP88BR8ctvJ++Mqu+RHQjW3PbILfnhaiQcREREREZEIZawtO6sX0tLS7Ny5c70OQ/y1a7WbcJFQDRKrQ0JViD4l1/XjM/D1g3D3/MCO9hQREREREZGQMcbMs9amFfaYRg5I8NRqAUlprpFlxZqnJx0A2l3mtss+DmloIiIiIiIiEhpKPIi3qiZBSi9YqsSDiIiIiIhIJFLiQbzXbjjsXA47V3gdiYiIiDcObIWsI15HISIiEhRKPIj32l4KJkpVDyIiUj4d2QPP9oSXB8Lh3V5HIyIiEnBKPIj3KteFxr1dn4cy1OxUREQkIGY8C8cPwO418OalLhEhIiISQZR4kPCQOsK94Nq+2OtIREREQufIHpj1oqv+G/UOZK6CNy+Do3u9jkxERCRglHiQ8NBmGETFaLmFiIiULzOfg6yD0Pc30HwgXP027FgObw53I6lFREQigBIPEh4q1ICm/bTcQkREyo+je121Q5thULed29fyQrj6Tdi+BN4aAccOeBujiIhIACjxIOEjdQTs2wRb5nkdiYiISPDNfN71duj725P3txoCV74KWxfA21fC8UPexCciIhIgSjxI+Gh9EUTHabmFiIhEvqP7YOYL0OYSqJd6+uNtLoERr0DGHHjnKsg6HPIQRUREAkWJBwkfCVWh+QWw7BPIy/M6GhERkeCZ+Twc3396tUNB7S6D4WNg0wx452rIOhKy8ERERAJJiQcJL6nD4eBW2DzT60hERESC4+g+l3hofTHUa3/2Y9tfAZe/CBu+h/dGQfbRkIQoIiISSEo8SHhpORhiEmHpR15HIiIiEhyzXvRVO/zGv+M7XAWXPQfrpsHXDwU3NhERkSBQ4kHCS3wlaDkIln8GuTleRyMiIhJYx/bDzGeh1UVQv6P/53W6BjqNhoXvwvGDwYtPREQkCJR4kPCTOgIOZ8KG6V5HIiIiElizXnTJB3+rHQrqeiNkH4YlHwY8LBERkWBS4kHCT4sLIK4SLNN0CxERiSDHDsCMZ6HlEGjQqfjnJ6VBnXYw77VARyYiIhJUSjxI+IlNhFZDYfnnkJPldTQiIiKBMftFOLYP+p1lksXZGOOqHrYthK0LAxeXiIhIkCnxIOEpdYR7cbZuqteRiIiIlN5P1Q6DoUHnkl+nw5UQkwDzXw9cbCIiIkGmxIOEp2YDIKGqpluIiEhkmD0Gju6FviWsdsiXWB3aXQ6LP4DjhwIT26ny8uCTO+CrUsYqIiLio8SDhKeYOGhzCaz8ErKPeR2NiIhIyR0/CDOegRYXQsMupb9elxsg62DweiHNeh4WvQOzXoCdK4JzDxERKVeUeJDw1W64e2G15huvI5F81sLRfV5HUTwHtsEHN8KhTK8jEZHyavZLvmqHBwJzvZSeUKsVzAvCcotti+CbR1zlYWxFmP6vwN9DRETKHSUeJHw16QsVasJSTbcIG5P/CE+murXKZcWMZ2DZJ7Dwba8jEZHyaPda+OE/0PwCSOoamGvmN5ncMhe2LwnMNQGyDsOHt0DFWjD8Zeh2i1vyuHtt4O4hIiLlkhIPEr6iY6DtpZA+wb0YEm9tXehePGcdhI0/eB2Nf7IOw4I33fcazyoioXZ4N7x9BZgoGPL3wF6740iIjgts1cNXv4Xda2D4GKhYE865291j+r8Ddw8RESmXlHiQ8JY6ArKPuOSDeCc3Bz6/GyrWdt3U103zOiL/LH4fju13/UK2LdKndiISOtnH4L1RsH8LjHoPajYL7PUr1HDJ+cXvQ9aR0l9v2ScuUdv7HmjSx+2rVMdVVix+D/ZuLP09RESk3FLiQcJbSi+okgQ/PuO6bIs3Zj4L2xfD0H9Ccg9YXwYSD9a6LvL12sNg3yeNWrYjIqGQlwef/Bw2z4LhL0JKj+Dcp+uNcHw/LP+0dNfZtwk+/z9omAb9f3/yY+f80lVs/PBU6e4hIiLlmhIPEt6iomHgI7B1Pix8y+toyqfda2HKX6D1xdBmGDTtCzuXw6GdXkd2dhu+d3F2/zlUbeiSWBrPKiKhMOkRlwy44E9u9GWwNDoXajaHea+V/Bq5OfDRrWDzYMTLEB178uNVG0Kn0bDgLTiwtVThiohI+aXEg4S/9le6N42THnVdwSV0rIVx97g1vkP/6RqaNennHlv/nZeRFW32i5BYA9pf4X5OHQGZK2DHcm/jEpHINudl+PFp6PYz1yMhmIxxozU3zyr52Mvv/gmbZ8LF/4YaTQo/pvevIC8Xfni6xKGKiEj5psSDhD9j3Jveo3vdJ+8SOgvfccsqBj4KVRq4fQ06QXxVWDfVw8CKsG8zrPwSulwPsYluX9tLXbmwmkyKSLCsmgDj74eWg90SL2OCf89O10BUbMmaTG78Eb77B3QYCR2uOvNx1Ru7ZpbzXgv/ajcREQlLSjxI2VCvvfv0aM7LgR0dJmd2aCdM/L2rNul604n9UdHQuHd493mY+4rbdrvlxL5KdaDxeW65hbXexCUikWvrAvjwJqjXAa4Y6yYzhULFWtDmYlj0rmto6a+je90Si2opcNETRR/f+9eQexxmPFvyWEVEpNzyK/FgjBlsjFlljFljjHmgkMcbGWMmG2MWG2OmGmOSCjyWa4xZ6Pv6vMD+JsaYWb5r/s8YExeYX0kiVv/fQ2J192mS3jgG34QH3ESRS/4DUac8VTTt65qR7VnvTWxnk33UffLXaqh7QV1Q6gjYs85NuBARCZR9m+Cdq6FCLbjmfYirGNr7d70Rju2DFZ8XdaRjLXzxf3BoO4wYC/GViz6nVnNoN9x9AHBkT2miFRGRcqjIxIMxJhp4FhgCtAVGGWPannLYE8Ab1toOwGPAXws8dtRa28n3NazA/r8DT1prmwN7gQIfTYoUIrE6nP8IbJoBSz7wOprIlj7RVQacdx/UbnX64036um04Vj0s+RCO7oEePz/9sTaXQFSMmkyKSOAc3QdvX+mqDUZ/AJXrhj6Gxn3ccgh/mkzm5cHM52D5Z9D/QUjq6v99zrsXsg7BrBdKGqmIiJRT/lQ8dAfWWGvXWWuzgPeAS085pi3wre/7KYU8fhJjjAEGAB/6dr0OXOZnzFKedb4OGnSBrx+CYwe8jiYyHT8I434Ntdu4ee6Fqd0KKtWDdWGWeLDWNZWs09YtqzhVhRrQbICbV6+qGREpraP74L3RbvrPyLegTmtv4oiKck0mN/4Au1YXfkz2MVcN9lwPt4yu2QA491fFu0/dtm7C0awX4Nj+UoctIiLlhz+Jh4bA5gI/Z/j2FbQIGO77/nKgsjGmpu/nBGPMXGPMTGPMZb59NYF91tqcs1xT5HRRUTD0CTi0wzXEksCb/Cc4sAWGPQ0xZ1gBZQw06eMmW+TlhTa+s9k00/UA6X7rmZu6tRsO+zdDxpzQxiYikWX7UhjTz02EuPwF95zopU6jXUXXqVUPR/a4yRVPtYcvfgkxCTDiFbjmg9OX0fmjz/0u6TD7pYCELSIi5UOgmkveB/Q1xiwA+gJbgFzfY42stWnANcBTxphmxbmwMeY2X+JibmZmZoDClTItqaurfJj5PGSu8jqayLJ5Nswe4964J3c/+7FN+8KRXbAzjMZTzn4REqpCh6vPfEzroRAdr+UWIlJyiz+Alwe6njI3fnlibK+XKteFVkPcNKKc47B3I3z1W3gyFb79s2vSfP1n8PPvXLwlbX7ZoBO0uNA1mcw6HNBfQUREIpc/iYctQHKBn5N8+35ird1qrR1ure0MPOjbt8+33eLbrgOmAp2B3UA1Y0zMma5Z4NpjrLVp1tq02rVr+/lrScQb+Khr3vXVb1QyHyg5WfD5L93YzPMfLvr4cOvzcGArLP/cJaXO1tgtoSq0uACWferm0ouI+Cs3272Z//hn0KCzexOf0tPrqE7oeqPrcfPaxfB0Z9cIsu0wuP0HuO5jaNovMCM++9zv7jP31dJfS0REygV/Eg9zgBa+KRRxwEjgpLbJxphaxpj8a/0OGOvbX90YE59/DHAusNxaa3G9IPI/IrgB+Ky0v4yUIxVrQf+HYN1UWPGF19GUfdbCNw9D5gq46N/+dTivlgw1moZPn4e5Y8HmubGrRUkd7rq5b/wx+HGJSGQ4uB1ev8T1N+hxB9zwuTeNJM+m6QCo2Rx2roBed8L/LXbLQOqlBvY+yd1d8vnHp13Vh4iISBGKTDz4+jD8ApgIrADet9YuM8Y8ZozJn1LRD1hljEkH6gKP+/a3AeYaYxbhEg1/s9bm12X/Fvi1MWYNrufDKwH6naS8SLsZ6rRzTbKyjngdTdk2/V8w63nocTu0Guz/eU36umZmudnBi80fOcfduuaWg6FGk6KPbzkYYivAso+DHpqIRIBNM+HFPm4U74hXYMjfIDrW66hOFxUFP5sM966AC/8MVYPYPqvP/a7f0oK3gncPERGJGH71eLDWjrfWtrTWNrPWPu7b97C19nPf9x9aa1v4jvmZtfa4b/+P1tr21tqOvu0rBa65zlrb3Vrb3Fp7Zf45In6LjoGh/3SNAr9/0utoyq45r8C3f4L2V8GgvxZ9fEFN+7rRalsXBCc2fy37BA5nQo/b/Ds+rqJLPiz/zPukiYiEL2th1ovw2kXueeNnk8Kjn8PZJFbzr2qttBr3hpRzYNo/3HQPERGRswhUc0kRbzQ+F9pfCT/8B/as8zqasmfpx/DlvdBiEFz2XPE7nDf2dXH3ernFrBehVkto2t//c1JHwJHd4dOjQkTCi7Xw2V2ul1DzC+DWKVC3nddRhQ9jYPBfXZPhyX/0OhoREQlzSjxI2XfBn1zJ68SHvI6kbFkzCT6+zTVGu/K1kpUNV6zpOqV7+eY9Yy5snQ/dbyte07TmAyG+Ciz9JHixiUjZtWE6LHwbzvkljHzHVRLIyRp0ckv05r4KmzWiWEREzkyJByn7qtSH834Nq76E9d95HU3ZsHk2/O86qN0aRr0HcRVKfq0mfWHzLO/6bMx6EeIqQ8eRxTsvNgFaX+Sak+ZopZeInGL2S5BYHfr/vvjVYOVJ/9+7aUjjfqWlayIickb6P6lEhp53QtVk12hSIxLPbsdyePtKqFzPjVcr7ad4TftBbhZsnhmI6Irn4HbX36Hz6JKtaU4dAcf3w9pvAx+biHgj64hbQrZ7bcmvsX8LrPzSjeeNTQxcbJEovjIM+TvsWAozn/c6GhERCVNKPEhkiE2EgY/C9iWw6D2vowlfe9bDm5e7f17XfQKV6pT+mim9ICom9H0e8vLg0zvd8orufjaVPFXTfu4TzaUfBTQ0EfHQrBdgzssw6dGSX2Peq77xvLcELKyI1vpiaDUUpv4V9m3yOhoREQlDSjxI5EgdAQ3TYPJjkHXY62jCz8EdLumQcwyu/RiqNw7MdeMrQVK30Pd5mP4vWDvZfdJWs1nJrhEdC22GwaqvNJJVJBIc2QPfP+XG5a74HHauKP41fhrPOyhwz5ORzhgY8g/AwPjfuMacIiIiBSjxIJHDGBj0Fzi0HX542utowsvRffDWcDi0E0Z/CHXbBvb6TfrC1oVwdG9gr3sm66bClMfdCNCuN5XuWqnD3UjQ1V8HJDQR8dD0f8HxAzD6A4itCNP/XfxrLP/cjeftdmvg44tk1ZKh/+8g/StYOc7/8w5lwhuXwRMtYdIfYX9G0EIUERHvKPEgkSWlB7S73I3XPLDV62jCx5THIXMlXP0mJHcL/PWb9gUsbPg+8Nc+1YGt8NHPoHYruPjJ4k2yKEzj86BiHS23ECnr9m2C2WOg0zXQuDd0uxmWflj8Xg9zXoIaTaHZgODEGcl63AF127uqh+MHiz5+6wIY0w82zXQTkn54Cp7qAO9fDxt/VOWEiEgEUeJBIs/AR8HmwuQ/eR1JeMjNhiUfuiUFzc8Pzj0aprnS5mD3ecjNhg9vdssirnrDLfMorahoaHeZq3jw54WyiISnKX8FDPT7nfu5190QFQvfP+n/NbYtdlN6uv1MkyxKIjoGLnkKDm6Dbx8/+7EL34VXBoGJglsmwrUfwS8XQq+73P9LXh0CL54H89+E7KOhiF5ERIJI/1eVyFO9MfS8Axa968r/y7u1U+DoHuhwVfDuERMHjc4Jfp+HyX+ETTPgkv+4iodAaTfc9b5Inxi4a4pI6OxY5p7ze9zmSv4BKteFrje4/fs2+3edOS9BTKKrmpCSSUqDtJth9ouF/z84Nxu+egA+vR2Su8NtU6F+R/dY9UZw4Z/g1yvc83xeLnz+C/h3W9csVJWMIiJllhIPEpnOuxcq1ICvH1Kp5pL33eSGZkGqdsjXpC/sSvfvhWHWYdenoTijT1eMgx//C2m3QIcrSxxmoZJ7uOUWK74I7HVFJDQm/RHiq0DvX5+8/5xfAsYtvyvK0b2w+AP3/JJYPShhlhvnPwwVa8O4X538PH94l2tyPOt5Nwb7uk+hYs3Tz4+rAF1vhDt+hBvGucT2D/+Bly+A3JwQ/RIiIhJISjxIZEqoCv1/Dxumw6rxXkfjnazDbhZ928tcVUIwNe3rtuu/O/txO5a7Nb1vXArP9nDLQPLyzn7OnvVudGaDzjD4rwEJ9yRRUdD6Ilj9DWQfC/z1RSR4NnwPqyfCefe4hHNB1ZKh0yiY/wYc3H726yx4G3KOqqlkICRWc82ety5wo00Bti1yz/2bZ8PlL7rn8uiYs1/HGGhyHox8G654FQ5kwNpvgx29iIgEgRIPErm63Ai1WsHXf4CcLK+j8cbK8ZB9BNoHuEKgMHXbQ2KNM/d5sNa9+H+pv5uyMegvbpzlR7fAC+e6TvKFVadkH3ONxoyBK1+HmPjgxN/mYsg+DOumBOf6IhJ41sI3j0DlBtDj9sKP6X0P5GW7iqkzyctzb5CTe0L9DsGJtbxJHeEadE7+E8x83vVzsHlw8wToOLL412s11P0/ZtG7gY9VRESCTokHiVzRMXDhn2HPWpj7itfReGPJB1AlCVJ6Bf9eUVHuk6n1005PIBw/CB/fBp/f7ZY13P69ayB2+w8w4hW35vf96+DFPrBqwsnnT/gtbF/sPiGr3ih48TfuA/FV3ZIOESkbVnwBW+a6MY6xiYUfU6OpS77OfRUO7y78mLXfwt710F3VDgFjDFz0L5f0mfAANOwCt01z25KIiXPJjJVfwrH9gY1VRESCTokHiWwtLoCm/WHq3+DIHq+jCa3Du2HtZGg/InTd2Zv0hQNbTh5ft32JK69d+iH0fxCu+8Q1fQMXV/sr4M6ZcNkLcPwAvHs1vDwQ1kyGRe/BvNfcJ5atBgc39pg4aHmhW5qjNcQi4S83xzWcrdUKOhbRDLL3r13116znC3989hjX56XNsMDHWZ7VaAqXPusmjVz/GVSqXbrrdRwFucdh+WeBiU9EREJGiQeJbMbAoMfdG9rv/ul1NKG1/BPIywnNMot8Tfu57fqprmph7lh46Xw4fghu+AL6/saNrzxVdIxbh/2LuXDJ03BoB7w1HD65HRr1hv4PhSb+Npe4CSCbZoTmfiJScgvehN1rYOAjRfcKqNMa2g6DWS+6pV4F7d3gxul2vTH4vXDKo/ZXQL8H3NK60mrYBWq2cElpEREpU5R4kMhXtx10vs59orVrjdfRhM7iD6B2G6ibGrp71mjqlnas+go+vBnG3QONz3VLKxr3Lvr86Fg3/u7ueTD0CWg1BK54peg3FYHSfCDEJGi6hUi4yzrsKtmSe7i1//447z6XhJ7z0sn757wCJgrSbgp8nBJYxrj+EBt/cAkjEREpM5R4kPKh/4PuDeW0v3sdSWjs3QibZ7pPmowJ3X2NcdMt1kyC5Z+6kWqjPyp+eW1MvFtrPepdqFwvKKEWKq6ia4a28ktvxrAePwQZc93ykm8edtM8ROR0M5+HQ9th4B/9f46r3wFaDoYZz7n/1gCyj7rKidYXQZUGwYtXAqfD1W67+H1v4wik8j72W0TKhRB9jCjiscp1XVOqpR+5KQmxCV5HFFxLP3Lb9leE/t6dr4XMla6xZ6NzQn//0mp9sevzsHVByZugFSUv1yUVdiyFncthxzL3tfeUREP2URhazpYIiRTl8G744T+u0qFRMRvnnncfvDIQ5r0K59ztniuP7lVTybKkWjI0Ps9Nt+hzf2iT68GwehK8OxIq1ISqSe73q5oEVfO3vu8Tq5f931VEyjUlHqT8aDMM5r8O66YGv1Gh15Z84EqQqzcO/b0bnQO3luE5662GgImGleOCk3g4vNs129y/yf1soqBGM/dpbKdr3NKgOm1dF/j0CTDkH3qxKeWDtbD6G9ebJq4ixFXybSue+DkmDqb/C7IOuYqq4kru5nrR/Phf6PYzmP0S1G7t3shK2dFxJHx2F2TMgeTuXkdTOjOfc0mFFgNhfwZsW+xGYeceP/m4Kklw0/jgTncSEQkiJR6k/GjSB+KruPX7kZx42LHMfYo+9AmvIymbKtRwfSlWfFGyNzZFWfSOSzoM+Yd7wVy7deFjAFsOcomHzFWuMZ5IpFv4tnszeTZRsS4x0Xk01GlTsvv0uR9eu8j1oNm20D1XKrlXtrQZBl/e55pMluXEw96NbpRr39+6kbD5rIXDu2D/ZpeM2LcRJj3qelUNetyzcEVESkOJByk/YuLc+t5VX0Luf0LXsDDUFr/vPrFvd7nXkZRdrS+Br+6HzHSo3TJw17UW5r3uqlF6/Pzsx7YY5LbpE5R4kMh37ABM+iMkdYeh/3DNI7MOu8qGrCMFvj8Medlwzv+V/F6NzoWUXq5UP66y+/RcypaEKtDmYrdUZvBfXV+gsmjBmy7p1fnak/cb43ojVap9ovJuyzx3fP8HIa5C6GMVESklNZeU8qXNJW4978YfvI4kOPLy3AuxZgOgYi2voym7Wl/ktisDPN1i44+wezV0uaHoY6s2hLrt3Zg/kUj33T/gcCYM+Ts06Oym4LQc5HrzdLkOet4Ofe5zozMv/HPxG9YWZIy7FrikQ3zlwPwOElodR8KxfZA+0etISiY3Bxa85aYpVUsu+vjut8Gx/W4ppYhIGaTEg5Qvzc+HmMTIHZe4eaYrzexwldeRlG1VG0LDrrBiXGCvO/91iK/qfzVKy0GwaaZLlolEql1rYOYLbvlEsBq6nqrZ+TD8Jej/+9DcTwKvST+oVM8ttwi0o3tdP6hgWv01HNwGXW/07/iUXm489uyXNAVDRMokJR6kfImr6Bo4rRznqgMizZIPXGLF37n2cmatL4at89362kA4sgeWfeqSQv6WybYcBDYX1kwOTAwi4Wji71yfk/MfCd09jXH/LVaoEbp7SmBFx7jJTasnuqa9gfTpnfDGpW68cbDMe80lTvKX1RXFGNcQdccSl5AWESljlHiQ8qfNMPcpw5Z5XkcSWDlZsOwTaD0U4it5HU3Z1+YSt135ZWCut/h/rkt5Vz+WWeRr2NWNWCurpcQiRUn/2n3y2/c3UKmO19FIWdNxlGs2mj9COhDWT3cjlQGm/CVw1y1ofwas+cb1dihOv6kOV7mqudljghOXiEgQKfEg5U+LC11n9BWfeR1JYK391pWHttcyi4Co1QJqtQrMspz8ppINu0K99v6fFxUNzS9wL1Dzcksfh0g4ycly1Q41W0D3IpqtihSmXqrrhbM4QMst8vLg6weharJr4rh2cnCqCxa8BTbP9S8pjriKLlmx4nM4uD3wcYmIBJESD1L+JFaDpn3dG8pIWie55H03C7zZAK8jiRxtLnENIUtbxrt5NmSu8H8tb0EtB7mEUsac0sUgEm5mvQC71/imEsR5HY2UVR1HugrGzPTSX2vx/2DbIrfsp9ddULF24Kse8nJh/pvu/9XVGxf//G63uCqPea8FNi4RkSBT4kHKpzaXwN4NsGOp15EExvFDsHK8a1qoF/CB0+Zi12Mh/avSXWfeaxBXCdoNL/65zQa48ahabiGR5OAOmPYPt769xQVeRyNlWfsrwUSVvuoh6whMfgwadHHTVOIqQu97YP002PB9YGIF17PnQIZ/040KU7OZq4SbO9ZVDYmIlBFKPEj51OoiwETOdItV4yHnqHsBJoFTv5MruS3NdIuj+1zvjfZXlqz3RmI1181ciQeJJJMfg5xjrtpBpDQq13VTShb9r3RNo2c8Awe3wqC/QJTv5XHazVCprqt6CFSF5PzXXSVFaZpAd78NDu0I/MhnEZEg8ivxYIwZbIxZZYxZY4x5oJDHGxljJhtjFhtjphpjknz7OxljZhhjlvkeu7rAOa8ZY9YbYxb6vjoF7LcSKUql2tDonMhJPCx+371BTu7pdSSRxRg33WLtt66qpCSWfOCSQiVZZpGv5SDYuQz2bS75NUTCxZZ5sPAt6HmH+/RWpLQ6jnRVBBtLWJlwcDt8/5SrhmzU68T+2EQ4717Y+AOs/670cR7cDqu+gk7XlK46sflAt0xj9kulj0lEJESKTDwYY6KBZ4EhQFtglDGm7SmHPQG8Ya3tADwG5H+EcQS43lrbDhgMPGWMqVbgvPuttZ18XwtL9ZuIFFebS2DncjdDviw7vMu9MU4dceJTGgmcNhe7aRRrvin+uda6ZRb1O0KDTiWPoaVv3Nrqr0t+DZFwkJcHX/0WKtaBPvd7HY1EilZDIa6yq3ooiSmPQ24WDPzj6Y91uQEqNwhM1cOCt9zyvZIus8gXFQXdboVNM2Db4tJdS0QkRPx5l9IdWGOtXWetzQLeAy495Zi2wLe+76fkP26tTbfWrvZ9vxXYCdQOROAipdb6Yrct66WKyz5xL2S0zCI4Unq5kZYlWW6xZb7rI1KaageAWi3dp1tabiFl3ZL3XaPUgY9CQhWvo5FIEVcB2l0Kyz91vRqKY8cylxDofmvhFTixCdDnXtg80yX5SyovD+a/AY3PC0ylT+fREJMIc1T1ICJlgz+Jh4ZAwfreDN++ghYB+V3TLgcqG2NqFjzAGNMdiAPWFtj9uG8JxpPGmPjCbm6Muc0YM9cYMzczM9OPcEX8VC0ZGnQu+8st0idAzeZurJgEXlQ0tBriqg2K28hr3qsQWxFSryhdDMa4JnzrpxX/RbVIuDh+EL55xI2V7TjK62gk0nQcBVmHYOWXxTvv64cgvsrZK3A6X+eWM5am6mH9VNi3sfSJ6HyJ1aHDVbD4AziyJzDXFBEJokDVZd8H9DXGLAD6AluAn4bOG2PqA28CN1lr8zv//A5oDXQDagC/LezC1tox1to0a21a7doqlpAAa3OJW2+8P8PrSEomJ8uNe2za3+tIIlubYXD8QPHW+B47AEs/gtThgflkt+Ug14xvw/TSX0vEC9P/BYe2w5B/aFmYBF7KOVA1xS1v8zdJvHqSq2Lo+xuoUOPMx8XEQ5/7YMtcWF2CZXfg4kqs4V53BEr3W10PoYVvB+6aIiJB4s//+bcAyQV+TvLt+4m1dqu1dri1tjPwoG/fPgBjTBXgS+BBa+3MAudss85x4FXckg6R0GozzG2L+wlJuMiYA9lHoGlfryOJbE36unGYKz73/5ylH7p/N11vCkwMjXu76gktt5CyaM96mPGs+1Q6Kc3raCQSRUVBj9tcg8nnzyk6QZCb46odqjdx/RKK0mk0VGvk+kEUt+rh0E73OqPTNS6JESj12ruEy5yXSzfRQ0QkBPxJPMwBWhhjmhhj4oCRwEmvvo0xtYwx+df6HTDWtz8O+ATXePLDU86p79sa4DJgaSl+D5GSqdUCarcuu8st1k9z88sb9/Y6ksgWmwAtLnBjS/Nyiz4e3KdbdVOhYZfAxBATD836u8RDoMa6iYTKtL+756rzH/E6Eolk59wN13wAWHj7Cnj7yjM3kF7wJmSugAv+6N+EiehYVxmxbaGbTFEcC9+BvBzocn3xzvNH91th7wZYMynw1xYRCaAiEw/W2hzgF8BEYAXwvrV2mTHmMWOM7+Ni+gGrjDHpQF3gcd/+q4A+wI2FjM182xizBFgC1AL+HKDfSaR42gxzo7IO7/I6kuJbNw3qd3JrPSW4Wl8MhzNh6l+LHq25dQFsW+TW8hoTuBhaXOhGxu1cHrhrigRb5ipY/D/3BqlKfa+jkUjX8kK4YwZc8CfYOAOe6+kqG44dOHHM8YOuciGl14nKR390GOkqJKb8xf8KA2th/uuuMqF2q+L9Lv5ocwlUqgezxwT+2iIiAeTXIktr7XhrbUtrbTNr7eO+fQ9baz/3ff+htbaF75if+ZZPYK19y1obW2Bk5k9jM621A6y17a21qdbaa621RbySFwmSNpeAzSt7yy2OH3LrTbXMIjRaXwStLoLv/gn/6Qg//heyjxZ+7LzXXbfxQE8aaXGh26ZPCOx1RYJpyl8gtgKc+yuvI5HyIiYOzv0l3D0POlztnq//29VNr8jLg++fconkCx8vXnI4Ogb6PQA7lsBKPycdbZgOe9ZB11KO0DxjTLGQdrMb+bx7bdHHi4h4RN2dROq1d+s2/VlucWgnfPxzmPxY8OMqysYfXelmEyUeQiI2EUa9A7d84yaIfP2QS0DMGgM5x08cd/wQLPkA2l0OidUCG0OV+lC/I6R/HdjrigTLtsVuxGHPO6BiLa+jkfKmcl247Fm49Vuo3gg+uwteHgAznnHThpK6Fv+aqVdAzRau+s2fqod5r0NCVWh76iT6AOp6I0TFwpxXgncPEZFSUuJBxBhX9bBuKhzbX/gx1sLi9+HZ7rD4PZj+b9jucVuS9dMgOh5SenobR3mT3B2u/wxuHO/GmH51PzzdBea+CrnZbpJF1qHAjUw7VcvBkDFb49OkbJjyF/emq9cvvI5EyrOGXeHmr+HyMXBgG2BgYAn7jeRXPexc7pJqp7LWjT0+sBW2zHdNiTuMdMnrYKlcF9oOcxUdWYeDdx8RkVIwtgw1KUtLS7Nz5871OgyJRJtmwdgLYfjL0OGU8viD22HcPa6xYFI3GPQXeOsKaHwujHrXm3gBnj/Xjf+6oYw2xowE1rqE1ZTH3YSRao1cA72YBLhzRmD7O+TLmOc+sRv+kpvhLhKuMubCy+fDgIegz/1eRyPiZB12idtqyUUfeyZ5uW5yxrH9UKet2x7b57ZH90Fe9snH3/Ej1G1XmqiLtm4avDEMRr7jlgaKiHjAGDPPWlvo+KqYUAcjEpaSukGlur5PJnyJB2th0bsw4QFXSn/hn6HnnRAV7daOfvsn2DzbfQIeaocyYcdSGPCH0N9bTjDGTZpo2s+NbpvyuOt4PuSfwUk6ADToDBVruz4PSjxIOPv2z1ChJvS43etIRE6Iq+i+SiMqGgY9DhMfdAmHhKpQLcVtE6u5bYJvW71R8JMOAMk9IDoONs1U4kFEwpISDyLg5n+3vtglGrKOwNG9MO5XsPprSO4Jlz4LtZqfOL7H7TDrBdfr4YYvgvcm80w2fOe2TfuF9r5SOGNcJ/UWF8D2JW6MZrBERbkmkyvHuTn00XoalzC04XtYN8UlbOMrex2NSOA1H+i+wkVsgktMb5rpdSQiIoVSjweRfG0ugewjMP5+N35r/XQY/De4afzJSQeA+Epw3n2uW/W6qaGPdd00iK/iRmlK+DAG6ndwyYFganGhK+ndPCu49xEpCWvh28fdiL9uP/M6GpHyI6WnG+d8polLIiIeUuJBJF/j3q40cuFbbtLFHT+4TuxR0YUfn3YTVElyVQ+h7pWyfpqLV592l0/NBkBUDKye6HUkxZd97OQpIMG04XtY+G7o//ss79Z+C5t+hD73BbehnoicLKWX6y+xdYHXkYiInEaJB5F80bEw7GkY9gzcMA5qNjv78THxrrP11vmw8svQxAiwd4P70jKL8iuhCjQ6B9LLYOLhjWHwzxausmjHsuDdx1o3Ou/T2+GDG+DYgeDdS06w1vV2qJoCXa73OhqR8iW5h9tumuFtHCIihVDiQaSgtpdCl+v8L5XvOMqNVPz2z67LdSism+a2TfqG5n4SnloOhsyVsHej15H4b99mtzykakOY95rrCv/ywOCMgMuY6xJ0zS+AFePgpQGwc0Vg7yGnW/WVS8b2/Y1LzopI6FSoAbVauUldIiJhRokHkdKIjoH+D0LmCljyYWjuuX6aWztdu1Vo7ifhqcUgt5071ts4iiN/aciVr8O9q9xo2mP7XWXCv1rDl/e65pyBsOR9iI6HK16BGz5393lpQOj+Oy2P8vLcZJcaTV1SVkRCL6UnbJ7p/nsUEQkjSjyIlFbby1xPiKl/gZys4N7LWlj/HTTpE/pJGhJeajWHDiPhh6fg+6e8jsY/6ROhehOo1cJ9MtfrLrhrNtw0AVoNhflvwgu9XYJg/Xclv09uNiz9GFoNduPsGveGn38H9TvCR7fA+N8E/7/V8mj5J27Mb7/fq/+MiFdSerpEa+ZKryMRETmJEg8ipRUVBQMedmXdC94M7r12LofDmervIM6lz0LqCJj0CPzwtNfRnF3WYbdMqOXgk5NmxkCjXjD8Rbh3JQz+OxzKhI9vcwmEklg3FY7sgvZXndhXpb4bfdvzLpj9Irx2EezfUqpfSQrIzYEpf4XabSB1uNfRiJRfKT3dVn0eRCTMKPEgEggtLoDknjDtH8EdY5U/urOp+jsI7lPly8dAu8vhmz/Aj894HdGZrf8Oco9Dy0FnPqZCDeh5O1z0BBzcBivHlexei993lQ4tLjh5f3QsDP4LXPmaS+K92OdEzxQpnSXvw+7V0P/3Z54EJCLBV70JVKqrccsiEnaUeBAJBGNg4CNwaDvMfil491k3DWo0g6pJwbuHlC3RMTD8Zbfk5+sHYcazXkdUuPQJEFcZGp1b9LHNB0L1xjBrTPHvk3XYTZlpe9mZmxu2uxxu/RYq1IQ3Lwvff2ZlRU4WTP0b1OsAbS7xOhqR8s0YV/WgigcRCTNKPIgESqNz3Bum7/8dnNF9udmw8QdVO8jpomNgxMvQZhhM/D3MfN7riE5mrevv0HwAxMQVfXxUNHT7GWz6sfjNJleOh+zD0OGqsx9Xu5VLPrS+yP0zW/5Z8e4jJ6waD/s2uka76j0j4r3knrBvExzY6nUkIiI/UeJBJJAGPARH98KMIJS8b5kPWYfU30EKFx0LV4x1nzhPeABmvuB1RCdsW+SWTrQc7P85na+FmMTiVxAteR+qJEHKOUUfG18JRrwCSd3g0zthp5qxlUj6BEis7hKvIuK9n/o8zPQ2DhGRApR4EAmkBp2h7aWudPvwrsBee/00wEDj8wJ7XYkc0bFwxavQ+mKY8NuSLVUIhvSJgIHmFxR56E8Sq7uqhcXvu2SePw7vgjWTof0I1/TVHzHxcNUbEFsB3rvGdYMX/+Xlwuqv3b9bTbIQCQ/1OkBsRSUeRCSsKPEgEmj9H4TsI/DVb+DAtsBdd91UqN/BNeATOZP85EOri+Cr+4Pbc8Rf6RMgKQ0q1S7eed1vhZyjsOAt/45f9gnY3JOnWfijSgOXfNi30U3TyMsr3vnl2ZZ5cGS3G10qIuEhOgaSuqrPg4iEFSUeRAKtdis4525Y+hE82RbeuRpWjCv5aEBwDfM2z4Ym6u8gfoiJc5MbWg2F8fe5qgGvHNwBW+effZrFmdRr75ZMzHnZfbJelMXvQ522UC+1+Pdq1AsG/dUlSb77R/HPL69WfQUmGpqd73UkIlJQSi/YsRSOH/Q6EhERQIkHkeC44DG4ez6c+yvYuhD+Nxr+3Ra+/gPsWl38622aAXnZ6u8g/ouJgytfd03GvvotHNnjTRyrv3bblkNKdn6P22DvBlgz6ezH7VkPGbOh/ZUluw+4CouO18DUv8KqCSW/TnmSPtE11k2s5nUkIlJQSk+weZAxx+tIREQAJR5EgqdmMzdi855lMOo918BuxrPwTBqMHQwL3naVDP5YNw2i49wnGCL+iomDi//t+hZM/qM3MaRPcM0e67Yr2fmtL4bK9WHWi2c/bsmHbtv+ipLdB9xEhov/DfU7wse3wq41Jb9WebBvE+xcVrJqFhEJrqRuYKLU50FEwoYSDyLBFh0DrYbAqHfg1ytg4KNwaCd8dic80w02zSr6GuunQVJ3iKsQ9HAlwtRtBz3vgHmvQ8a80N47+xisneLemJZ0zGJ0LKTdDGsnnzkRYK2bZpFyDlRLKXm8ALGJcPVb7r7/G60y5bNJn+i2xZlWIiKhEV8Z6qYq8SAiYUOJB5FQqlwXet8Dd8+DG75wVQyvDYUfn3FvngpzZA9sWwxN1d9BSqjvb6FSXfjy1/71SgiUjd9D9uHSvzHtcgNExbpeD4XZtgh2pUOHUiyzKKhaimvQuSsdPrvrzP9tlnfpE6FGM6jVwutIRKQwKb0gY27pekyJiASIEg8iXjAGmvSBn09zb8q+fhD+dy0c3Xf6seu/A6waS0rJJVSBQY/DtoUw77XQ3Td9IsQkQpNSjoCtXBfaXQ4L34bjh05/fMkHLjHR9rLS3aegpn1dr5bln8EPTwXuupEi67B7blK1g0j4Sunhkr/bl3gdiYiIEg8inkqo6sq6B/3FrYUf09d9elvQ+mkQVxkadvEmRokMqSOg8Xkw+TE4vCv497PW/U037eeWL5RW99vg+AFY/N7J+/NyXX+HFhcEftRsr1+4f26TH4M1kwN77bJu3VTIPa7+DiLhLLmn22q5hYiEASUeRLxmDPS6C24c78ohX74A5o49Ud69bho0PtetORcpKWNg6BOQdQgmPRL8++1c4ZoPtgrQJ+JJaVC/E8x+6eSlDxumw6HtpZtmcSbGwLD/Qu3W8MWvtOSioPQJEF9FDW9FwlnVhm7p2GYlHkTEe0o8iISLlB7w8+nQuDeMuwc+vg0yV8GetVpmIYFRp7VLci14y7+mpqWR7htH2eLCwFzPGFf1kLnSJRvyLf7AVQS1KuG4zqLEVYSuN8H+TXBgS3DuUdbk5UH619BsgJucIiLhK7mnq3hQ4lREPKbEg0g4qVgTRn8I/R+CpR/CS+e7/WosKYHS5zdQpSGMvxdyc4J3n/SJbixllQaBu2bqCEisAbPHuJ+zj8GKz6HNJYFZznEmyd3cdvPs4N2jLNm+yFWZBCvZIyKBk9ITDu2Aveu9jkREyjklHkTCTVQU9L0frvsUYhPcm8Q6bb2OSiJFfCXXU2T7Epj7SnDucXg3ZMwOfOPB2AToegOs/BL2bXZVFccPBG6axZnUTXVNMjPmBvc+ZUX6RMBA8wu8jkREipK/HCrYVW4iIkXwK/FgjBlsjFlljFljjHmgkMcbGWMmG2MWG2OmGmOSCjx2gzFmte/rhgL7uxpjlviu+bQxJR3yLhKhmvaFX8yBW752ZeYigdL2UmjaH779MxzaGfjrr5kENi84jQfTbnbbuWPdNItKdYO/FCk6Fhp0dskUcQmf5O6uQktEwlvt1q6R9aYZXkciIuVckYkHY0w08CwwBGgLjDLGnPrx6xPAG9baDsBjwF9959YAHgF6AN2BR4wx1X3nPA/cCrTwfWkml8ipEqtD1aSijxMpjvxGk9lH4es/BP766RNcQqB+58Bfu1oKtBoK816F1V+75RdR0YG/z6mS0tzEmZzjwb9XODuwDbYu0DQLkbIiKgqSe2iyhb/WToE967yOQiQi+VPx0B1YY61dZ63NAt4DLj3lmLbAt77vpxR4fBDwjbV2j7V2L/ANMNgYUx+oYq2daa21wBvAZaX7VURExG+1msO5v3TjKTf8ELjr5ma70ZMtLnQveIOh+21wdC/kZgVnmkVhkru7+21bHJr7havVX7ttoJfRiEjwpPSEXavgyB6vIwlvudnw7iiY8hevIxGJSP68KmwIbC7wc4ZvX0GLgOG+7y8HKhtjap7l3Ia+7892TRERCabz7oOqyTD+PveCKxA2zYDj+4P7xrRJH6jdBmq1dEsgQiHJ12CyvC+3SJ/o/mbUd0ak7Mjv87DZjz4PezfAuF/DocyghhSWti+BnKNuKyIBF6iPo+4D+hpjFgB9gS1AbiAubIy5zRgz1xgzNzOzHD4JiogES1wFGPw32Lkc5r4amGumT4ToOGjaLzDXK4wxMPoD9xWq/ieV60HVFMiYE5r7haPsY7Buiksqqe+MSNnRoDNExRbd52HrAnj5Atd4eOlHoYktnOQ/v+9a7Z7vRCSg/Ek8bAGSC/yc5Nv3E2vtVmvtcGttZ+BB3759Zzl3i+/7M16zwLXHWGvTrLVptWvX9iNcERHxW+uL3Kf5c18JzJz39AnQ+Dw3PSOYqiVD9cbBvcepktJgczlOPGz4HrKPaJmFSFkTm+iSD2ebbJH+Nbx6EcQkQIVa/lVHBMKqr+DovtDcqyj5iQebC5krvI1FJAL5k3iYA7QwxjQxxsQBI4HPCx5gjKlljMm/1u+Asb7vJwIXGmOq+5pKXghMtNZuAw4YY3r6pllcD3wWgN9HRESKwxjodA1kroRtC0t3rV1rYPeayH1jmtwdDmTAga1eR+KN9K8gtgI07u11JCJSXCk9Yev8wj/Jn/c6vDvS9f752TduOVsoqrt2rXb3nTu26GNDYfNsqNPOfb99qbexiESgIhMP1toc4Be4JMIK4H1r7TJjzGPGmGG+w/oBq4wx6UBd4HHfuXuAP+GSF3OAx3z7AO4EXgbWAGuBrwL1S4mISDG0uxyi42Hhu6W7TvoEt43UiQc/9Xkoh1UP1rplNE37Q2yC19GISHGl9HQNcrcuOLHPWvj2cfjil9CsP9z4pVtWltwD9m+G/YUWIwfOuqluuyMM3uQf2gn7NkKHq1yCNRxiEokwMf4cZK0dD4w/Zd/DBb7/EPjwDOeO5UQFRMH9c4HU4gQrIiJBkFgdWg2BpR/ChX+GmLiSXSd9gms6WL1RYOMLF/U6uARNxhxoe+pwpwi3c7l7I9Lnfq8jEZGSSO7htptmQKNekJMFX/wfLHoHOl8HFz8J0bG+Y7u7bcZsqHp58GJaO8VtdywL3j38lZ9QTunp/j+migeRgAvSrDMRESlTOl0DR3bDmm9Kdn7mKtgw3VVPRKqYOKjfsXz2eYj0ahaRSFexlpsEtHkWHDsA71zlkg79fg/D/nsi6QBQrz3EJLqlB8GSm+P+n2Gi3JKLnOPBu5c/Ns92DTjrd3S//44lgel7JCI/UeJBRESg2QCoWBsWlXC5xcznXTVA2s2BjSvcJHd3vTBysryOJLTSJ7rmdJXreR2JiJRUcg9X8fDqUPem/9LnoN9vT59SEx0LDbsEt8Hk1gVw/IBrcGxzXfLaSxlzXcIhNhHqpcKx/bA/w9uYRCKMEg8iIuJeaLa/ClZNgCN7ij6+oMO7YdF70PFq96laJEtKg5xj7tOwYEj/GsYODu4njcV1eLeLJ1KbhoqUFym93Bvqvevhmveh8+gzH5vcHbYtguyjwYll3VTAQM873c87lwfnPv7IzXGNN/OXmNRt77bq8yASUEo8iIiI03Ek5GUXf377vLGQc/TEC8hIlpS/9nlucK6/4A33ieTYQTD5sfCorFj9NWC1zEKkrGs9FDqNhpvGQ/Pzz35scg/Iy4GtC4MTy7qpUL+Da9obHedtn4edy9yo4PwGwnXbuq36PIgElBIPIiLi1O8AdVNd9YK/co7D7Jeg2flQp03wYgsXVRtC5QbBqUjIy4MN37vGlZ2ugen/gpcGeN94LX0CVKoH9Tp6G4eIlE5idbjsOdfHoCj5SdZgLLfIOuyu27Sfq7ar3crb57n85/P8xEN8ZajeJHiVbSLllBIPIiJyQseRsGWua/blj6Ufw6Ed0KscVDvkS+4WnJGaO5fD0b3Qcghc+iyMfBcObYcx/eD7pyAvN/D3LEpOFqz91lU7ROklg0i5UbEm1GgWnCTrxhmuuq5pP/dznXbeLrXImAOV6kK1lBP76qWq4kEkwPQqQkRETmh/lesy7k+TSWth5rNQu7WreCgvkrq5ee+Hdgb2uhumu23j3m7beijcOdO96Z/0iGsIt2ddYO9ZlE0zXAM49XcQKX+Se7jKhEBPd1g3xTUjTunlfq7bDg5uK35/oUDJmOOe1ws22azb3j3fZh32JiaRCKTEg4iInFC5rksiLPqfK/0/mw3fw/YlrrfDqV3RI9lPfR4CXPWw4Xuo3hiqJZ/YV7EWXPUmXD4Gdq6A53vD3LGhG/O26iuISYCmfUNzPxEJH8nd4ciuwCc8102DlB5uggSc6KngxXKLw77fL3+ZRb56qYCFHR5WYohEGCUeRETkZJ1GwYGME5/An8mMZ6FCTehwVWjiChf1O7p574EsQc7v75Bf7VCQMW5iyJ0/uqka4+6Bd0e5TuzBZC2s+tKVQ8dVDO69RCT8JPdw20A+1x3KdL0T8pdZgFtqAd4st8hvFJw/0SJf3VS33b44tPGIRDAlHkRE5GSthkJ81bMvt9i91jUdTLvlxKdW5UVsgmvEGcjJFjuWwrF90LjPmY+pmgTXfQrnPwLpX8Hi/wXu/oXZuRz2bXJ/DyJS/tRuDfFVICOAiYf109y2YOKhcj1IrOHN+MqM2RAVA/U7nby/Wor7/6BGaooEjBIPIiJysthEaHcZLP8cjh8q/JiZz7tu5N1+FtLQwkZSNzf3PVBVB6f2dziTqCjofY97kTztb8Edt7lyPGDU30GkvIqKcs91gax4WDcFEqqe/EbfGNfnwYtlDZtnu+qGuAon78+PSQ0mRQJGiQcRETldx1GQfRhWfHH6Y0f2wMK3of2VridEeZTUzc193xmgNckbvocaTd24zqIYAwP+4KoR5r8emPsXZtV4t7SjvP47FhG3BGHHMjh2oPTXshbWToUmfSAq+uTH6rR1fWyK6i0USHm5sGX+6f0d8tVLdb97KGMSiWBKPIiIyOlSerpGh4Utt5j/unvT3fOOkIcVNvJfqAbik8C8XNjwQ9HVDgU1P991hP/un5B1pPQxnOrANlfR0WpI4K8tImVHcnfAujHLpbVnnesfVHCZRb667Vyye9+G0t/HXzuXu3ue2t/hp5hS3eN714cuJpEIpsSDiIiczhhX9bD+O9ifcWJ/bjbMGgNN+kK99t7F57VqKW7ueyD6PGxfAsf3n72/w6nyqx4O7YA5L5U+hlOlf+W2rS4K/LVFpOxomAaYwCRZ101x26b9T3+srq/BZCiXW+RPJjpbxQOoz4NIgCjxICIihetwNWBPbmK47FM4uBV63eVVVOHBGPdiNRBN1/zt73CqxudCswHw/ZOBKYMuaOV4qN4EarcK7HVFpGxJqOKSAptnlf5a66ZC1WS3rOxUtVsDJrQjNTfPgQq1XHVfYeq0BROlPg8iAaLEg4iIFK5GE0g5Bxa959bmWgszn4WaLaD5BV5H572kbq50+PDu0l1nw/dQszlUqV/8cwc8BEf3wsznShdDQccPus7zrS9yCRYRKd+Su7vqrtL0OsjLdRV0TfsV/rwSX8klAALVN8cfGbPd73am57nYRPfcrIoHkYBQ4kFERM6s40jYle4acG2aAVsXuN4OUfrfx0/lufnluiWRmwMbfyx+tUO+hl2h9cXw4zOu6WcgrP0WcrPU30FEnOQecPwAZK4s+TW2LYRj+wvv75CvbrvQVTwc2QO717gGumdTr70qHkQCRK8cRUTkzNpdBjEJrsnkjGchsbrr/SDQoDOY6NIlHrYvdi/oG59X8msMeAiyDsEPT5X8GgWtHO/+PSf3DMz1RKRs+6mZbimWW6yb6rZN+p75mLrtXBVZ9tGS38df+f15ks7QWPKnmFJh/yY4ui/oIYlEOiUeRETkzBKqupL7Re/Byi8h7ebT552XV3EVXPOx0vR5KGl/h4LqtHGjTWeNgYPbS34dcBUYqydCi0EQHVO6a4lIZKjR1PVCKE2DyXVToW57qFT7zMfUaQs2r3SVFf7KmOP6NzTscvbj8psoh7L3hEiEUuJBRETOruM1kHUQomKg261eRxNekrq7ZSh5uSU7f/10qNUSKtcrXRz9HnDLI6b/q3TX2TzT9YxoPbR01xGRyGGMW25R0oqHrCOwaSY0PUu1A7jqAgjNm/yM2a7CIq6inzFpuYVIaSnxICIiZ9e0nxsf2fHqkjVAjGRJ3dwyh50rin9ubo7rm1Gaaod8NZtB52th7quwb1PJr7PqK4iOc9MyRETyJXeHPWvh8K7in7t5pkuMFjZGs6AaTSAmMfgjNfNyIWNe0csswCWFK9R0Y49FpFSUeBARkbOLjoE7foSLnvQ6kvCTXIoGk9sWuqRFafo7FNT3N+6TyWl/L9n51rrlNE36QnzlwMQkIpEhuYfbluS5bu0UiIqFRr3OflxUtBvhG+zqgsxVroovv3fF2Rjjqh5U8SBSako8iIhI0eIrQ0yc11GEn+pN3KdhJXkxHoj+DgVVTYK0W2Dhu7BrTfHPz1wFe9drmoWInK5BJ7fcriTLLdZNdYmLopY1gHuTvzPIFQ/5fXmS/ah4ANfnYecKV6UmIiWmxIOIiEhJGePKdUuSeFg/HWq3hkp1AhfPeb+GmHiY+pfin7vqS7dV4kFEThWbCPU7Fr/B5OHdbnrP2cZoFlS3LRzOhEM7ix2i3zLmQGIN1zTTr5hSIeeYG78pIiWmxIOIiEhpJKXBrnQ3F95fudmu2Vqgqh3yVaoDPW6HpR8Vf/b8qq/ciNAqDQIbk4hEhuQesGWee/7y1/ppbut34qGd2wazweTmOW6ZhTH+HV9PDSZFAkGJBxERkdLIL9fdMt//c7YugOzDgevvUNC5v4T4qjDlcf/PObjDzbVvdVHg4xGRyJDc3X3yv32x/+esmwrxVVxS0x91fImHYC23OLoXdq060Z/HH7VauR4VajApUipKPIiIiJRGgy5uHnxGMUqQA93foaDE6nDu3bBqPMx5xb9z0icAVsssROTM8qdAbC7G0rJ1U12CNTrGv+Mr1YaKtYNX8bBlntv601gyX0xcaJpeikQ4JR5ERERKI76S+5SuOH0e1k+HOm2hYq3gxHTur6DFIPjyXlj2adHHrxrvRqbmlzmLiJyqakOomux/g8k962HfRv+XWeSr2y54iYfNc1yiuGHXYsaUWvzlayJyEiUeRERESqtRL9jwg1uuUJScLPfCPRjVDvmiY+HK19ya7I9vdZ86nknWYfd4q6H+r3kWkfIpqZv/DSbzn3ea9S/ePeq0g8yVkJdbvPP8kTHHJX2LOzK4Xioc2g6HdwU+JpFyQokHERGR0upzP1SpD+9cBbvXnv3YrfMh+0hw+jsUFFcBrnkParaA90afuQfF2ilu3XarocGNR0TKvuQecCAD9mcUfey6qVClIdRsXrx71G3nnpP2rCtRiGeUl+eSw0lpxT+3rq/BpPo8iJSYX4kHY8xgY8wqY8waY8wDhTyeYoyZYoxZYIxZbIwZ6ts/2hizsMBXnjGmk++xqb5r5j8WwHliIiIiIVSpDoz+CKyFt0ac/VOx/P4Ojc4NflyJ1eHaj6BCDXj7CthVyDi4VV+5ZpSNzgl+PCJStuU30y2q6uHwLjfRomm/4ldS1W3rtoFebrErHY7vP9GrojjqtffFpOUWIiVVZOLBGBMNPAsMAdoCo4wxbU857CHgfWttZ2Ak8ByAtfZta20na20n4DpgvbV2YYHzRuc/bq0N4sBeERGRIKvVHK55Hw5uc5UPWYcLP279dPfpWcWaoYmrSn247lO3rvnNy+DA1hOP5eW6xpItLnDLM0REzqZee4hJLDzxYK3rofDxbfDvNm6CRLvhxb9H7dbu+SrQiYf8PjzJJUg8VKwFleqpz4NIKfhT8dAdWGOtXWetzQLeAy495RgLVPF9XxXYyulG+c4VERGJTMnd4IqxblzmhzdDbs7Jj+ccdy/Yg9nfoTA1m8HoD+HoPnhzOBzZ4/ZnzIEju6C1llmIiB+iY11jxoJTfLKOwPw34MU+8MpAWDkeut4Id82GFgOLf4/YRKjRLPAjNTNmQ0I1d+2SqNdeFQ8ipeBP4qEhsLnAzxm+fQU9ClxrjMkAxgN3F3Kdq4F3T9n3qm+ZxR+MKbwOyxhzmzFmrjFmbmZmph/hioiIeKj1RTD0n66SYPx97lPAfFvmQc7R4Pd3KEyDTjDqHdizFt652lVkrPzSzadvXoI3ByJSPiV3h22LXEXCxAdddcPnd0NuNlz0L7h3hXsOrN2q5Peo2zawFQ/ZR2Hjj645ZlQJW9zVS4XMVa5BsIgUW6CaS44CXrPWJgFDgTeNMT9d2xjTAzhirS2YJhxtrW0PnOf7uq6wC1trx1hr06y1abVr1w5QuCIiIkHU7WfQ+x6Y9ypM/9eJ/Ru+B4x3/RSa9HEVGVvmwvs3uMRD496QUNWbeESk7EnuDnk58Pw5MOsF18fhxi/hzhnuua+4EyMKUzcV9q6H44dKf62Mea4aY/caSC3B0o+CMeVlw65VpY9JpByK8eOYLUBygZ+TfPsKugUYDGCtnWGMSQBqAfl9G0ZySrWDtXaLb3vQGPMObknHG8X9BURERMLSgIdh/xb49k+us3unUbD+O/epWYUa3sXV5hK4+Cn44pfu5x63exeLiJQ9jXtDyyGuiqrLDa6PTKDV8bWTy1xZsikU4CoTpv0dvn8SKteD6z6BZgNKHlN+g8ntS098LyJ+8yfxMAdoYYxpgks4jASuOeWYTcD5wGvGmDZAApAJ4Kt8uApX1YBvXwxQzVq7yxgTC1wMTCrl7yIiIhI+oqLg0mfh0A74/BduwkTGHEi72evIoOsNrvHbjGehzcVeRyMiZUl8ZTeqN5h+mmyxtGSJh+1L4JPb3fmdRsOgv0BitdLFVKMZxCSoz4NICRW51MJamwP8ApgIrMBNr1hmjHnMGDPMd9i9wK3GmEW4yoYbrf1pUWsfYLO1tuAw3nhgojFmMbAQl9B4KRC/kIiISNiIiYOr33Rd2t+7xs2m96K/Q2F6/wruXQVVGngdiYjIyao1htiKsKOYDSZzc+C7f8KY/nBoJ4x6Dy57rvRJB4DoGKjTBrYvLv21pPzasx6e7gJzX/U6kpDzp+IBa+14XNPIgvseLvD9cqDQgeTW2qlAz1P2HQa6FjNWERGRsiehKoz+AF4e6EZtNurldUQnlLTJmohIMEVFuTf5xZlskbnKVTlsne/GeF70r8Ava6ub6nrjWAuF98UXObsVn7smz+N+BTtXuGqcaL/ekpd5esUhIiISbFUauOZrV73hllyIiMjZ1W3nljUUnAxUGGthxnPwwnmwdwNc8Spc+WpweunUaw9H97gkskhJrJkEtdtAz7tg9ovwzpVu1HU5oMSDiIhIKNRo4ho7iohI0eq2c71oDm4/8zG5Oa5R7sTfQfPz4a5ZpZtcUWRMqW67XX0epASOH4KNM6DFBTD4LzDsv7B+uquI3L3W6+iCTokHEREREREJL/mTLXYsK/zx44fgvVEw/w047z4Y+Q5UqhPcmOq288W0JLj3kci0YbobydriAvdzl+vh+s/gyG54aQCsm+ppeMGmxIOIiIiIiISX/Df5OwtJPBzaCa9f7MrWL34Szv9DaHouJFaDWi1h1hjYpiaTUkxrJrmmqckF2h82PhdumwKV68Obw2F25M5bUOJBRERERETCS4Ua7s3YqZMtdq2BVy6AnStdlUOoRxRf+TpERcOrQ2Htt6G9t5Rd1sLqb6BpXzfxqqDqjeGWr10lxPj74Mt7ITfbkzCDSYkHEREREREJP3XbnbzUYvMcl3Q4ftA17G01xIOY2sLPJkH1RvD2lbDwndDHIGXP7rWwb6PrRVKYhCoukXbOL2HOy/DWCNfjJIIo8SAiIiIiIuGnTlvYtcp9+rvyS7e8IqEq3PINJHX1Lq4qDeCm8dDoXPj0Dpj2j6Knb0j5tmaS2zY7Q+IBXCXNhX+CS5+DTTMirudD+RgaKiIiIiIiZUvddpCbBd88DLNegAad4Zr3oWItryNzCZDRH8Lnd8OUx2H/ZrjoSYjW2yspxJpJULO5m3BVlM6jocl5UC0l+HGFkP7LEBERERGR8JPfYHLmc9ByMFwxFuIqehtTQTFxcPkLUDUJpj8BB7bBla9BfCWvI5Nwkn0UNnwPXW/w/5wISzqAllqIiIiIiEg4qtUSarWCbrfC1W+HV9IhnzFuqsbFT8HayfDaUDi4w+uoJJxs/BFyjkLzgV5H4ilVPIiIiIiISPiJiYdfzPY6Cv+k3eSmcHx4E7wyEK79BGo19zoqCQdrJkNMAjTu7XUknlLFg4iIiIiISGm1Ggw3jnNTN75+yOtoJFysmeQakcYmeh2Jp5R4EBERERERCYSGXaH1RbB5liZdCOzb5CazlPNlFqDEg4iIiIiISOA07ApH98DeDV5HIl7LH6OpxIMSDyIiIiIiIgHTMM1tt8zzNg7x3prJUDUFarXwOhLPKfEgIiIiIiISKHXaQkyiEg/lXU4WrJsGzc9300/KOSUeREREREREAiU6Bhp0goy5XkciXsqYDVkHtczCR4kHERERERGRQGrYFbYtgtxsryMRr6yZBFEx0KSP15GEBSUeREREREREAqlhV8g9DjuWeh2JeGXNJEjpBQlVvI4kLCjxICIiIiIiEkgNu7qt+jyUTwe3w/Ylrr+DAEo8iIiIiIiIBFa1FKhYGzLKWeJh+1KY97rXUXhvzWS3VX+HnyjxICIiIiIiEkjGuKqH8lbxMPmP8MX/QdYRryPx1ppJUKku1E31OpKwocSDiIiIiIhIoDVMg13pcGy/15GExpE9sPZbwLrfu7zKy3X/HJoP1BjNApR4EBERERERCbSGXQALWxd4HUlorPgC8nLc9ztXeBuLl7bMh2P71N/hFEo8iIiIiIiIBFrDLm6bMdfbOEJl2cdQvTFEx0FmOU48rJkEJgqa9vc6krCixIOIiIiIiEigJVaHms3dJ+CR7lAmrP8OUq+Ami3Kd8XDmklumU2FGl5HElaUeBAREREREQmGhmmwZS5Y63UkwbXiM7B5kDoc6rSBnSu9jsgbR/a4hqKaZnEaJR5ERERERESCoWFXOLQDDmzxOpLgWvox1G4Nddq6xMP+TXDsgNdRhV5+c00lHk6jxIOIiIiIiEgwJHV120ju83BgG2z8EdoNd1Mc6rRx+zNXeRuXF9ZMgsQa0KCT15GEHb8SD8aYwcaYVcaYNcaYBwp5PMUYM8UYs8AYs9gYM9S3v7Ex5qgxZqHv64UC53Q1xizxXfNpYzRrREREREREIkjdVNdsccs8ryMJnuWfAtYts4ACiYdy1uchLw/WTIZmAyAq2utowk6RiQdjTDTwLDAEaAuMMsa0PeWwh4D3rbWdgZHAcwUeW2ut7eT7ur3A/ueBW4EWvq/BJf81REREREREwkxMPNRrH9mJh6UfQ932UKuF+7laY4hJLH8NJrcvgsM7tcziDPypeOgOrLHWrrPWZgHvAZeecowFqvi+rwpsPdsFjTH1gSrW2pnWWgu8AVxWnMBFRERERETCXsM02LoAcnO8jiTw9m2CjNmQevmJfVFRULsV7FzuXVxeWPYJRMVAiwu9jiQs+ZN4aAhsLvBzhm9fQY8C1xpjMoDxwN0FHmviW4IxzRhzXoFrZhRxTRERERERkbKtYVfIPgKZETjpYdknbttu+Mn767QtX5Mt8vJgyUfQ7HyoWNPraMJSoJpLjgJes9YmAUOBN40xUcA2IMW3BOPXwDvGmCpnuc5pjDG3GWPmGmPmZmZmBihcERERERGREEhKc9tIXG6x9GNo0AVqNDl5f53WcGi7Gy9ZHmyaAQcyoP2VXkcStvxJPGwBkgv8nOTbV9AtwPsA1toZQAJQy1p73Fq727d/HrAWaOk7P6mIa+I7b4y1Ns1am1a7dm0/whUREREREQkTNZpCQjXYEmGTLXavhW0LTzSVLKiOryVgJFZ5FGbJBxBbAVoP9TqSsOVP4mEO0MIY08QYE4drHvn5KcdsAs4HMMa0wSUeMo0xtX3NKTHGNMU1kVxnrd0GHDDG9PRNs7ge+Cwgv5GIiIiIiEi4MMYtt9gy3+tIAuunZRaXn/5Y/mSL8tDnISfLTfZofRHEVfQ6mrBVZOLBWpsD/AKYCKzATa9YZox5zBgzzHfYvcCtxphFwLvAjb6mkX2AxcaYhcCHwO3W2vx6mzuBl4E1uEqIrwL3a4mIiIiIiISJhl3dm/Csw15HEjjLPoHkHlA16fTHqjSE+CrlY7LF2slwdC+0v8rrSMJajD8HWWvH45pGFtz3cIHvlwPnFnLeR8BHZ7jmXCC1OMGKiIiIiIiUOUlpYPNg60JofNrbprIncxXsWAqD/17448ZA7dblo8Hkkg8gsQY06+91JGEtUM0lRUREREREpDANu7ptpDSYXPoxYKDdZWc+pk5rV+VhbaiiCr3jB2HleLfcJDrW62jCmhIPIiIiIiIiwVSxFlRrFBkNJq2FZR9D495Qud6Zj6vTFo7ugcMRPJlw5XjIOQodtMyiKEo8iIiIiIiIBFukNJjcsQx2pRfeVLKg8tBgcskHUDUFkrp7HUnYU+JBREREREQk2JLSYP9mOLjD60hOd2QPZKb7d+yyj8FEQ9tLz35c7fzEQ4Q2mDyUCWu/hfYjIEpvq4uif0IiIiIiIiLBFs59HsbdA892g49/Dge2nfk4a2HpR9Ckj1s+cjaV6rimi5GaeFj+KdhcTbPwkxIPIiIiIiIiwVa/o6sUCLc+DznHYfU3UKulq2b4b1eY/m+3/1RbF8DeDZA6oujrGuOWW0Rq4mHx+1CnHdRt63UkZYISDyIiIiIiIsEWmwh124VfxcOG6ZB9GC58HO6a5cZCTv4jPNvDNU8sOJVi2ccQFQttLvbv2nXaQObKyJtssWc9ZMyGDld6HUmZocSDiIiIiIhIKCSluQaTeXleR3LCqgkQkwhNzoMaTWHk23DdpxATD++NgreGw05f8mDZp9BsACRW9+/addrA8QNwYEswf4PQW/qR2/pT+SGAEg8iIiIiIiKh0bCreyO+e43XkTjWQvpEV+UQm3hif7P+cPv3MPjvrkLj+XPggxtcc8zU4f5fPxIbTFrrplmk9IJqKV5HU2Yo8SAiIiIiIhIKDdPcNlz6POxcDvs3QctBpz8WHQs9b4e750PXG2DFFxAdD62G+n/9OhGYeNix1C0faa9lFsUR43UAIiIiIiIi5UKtFhBX2VURdLrG62ggfYLbtigk8ZCvYi24+Eno9jM4th8Sqvh//Qo1oFLdyEo8LPkAomKg7WVeR1KmKPEgIiIiIiISClHR0LAzZIRJxcOqCVC/E1SpX/SxdduV7B512kBmhCQe8vJgyUfQ7HyoWNPraMoULbUQEREREREJlYZdXbl+9rHCHz+yBxa8Be9cDVP/Hrw4Du+CjDnQakjw7gFQp61rThlODTVLatMMOJABHa7yOpIyRxUPIiIiIiIiodIwDfJyYPtiSO7u9h3KhJXjYMXnsP4793hMIqyZBF2u968iobhWfw3Ywvs7BFLt1pBzFPZtcFMzyrIl70NsheAnayKQEg8iIiIiIiKh0rCr26ZPhG2LYPlnsPEHsHlQvQmccze0vRQSqsLTXWDOS3D+w4GPI30CVK7vlloEU522brtzZdlOPORkuXGirS+CuIpeR1PmKPEgIiIiIiISKlXqQ5WGMP0J93PNFnDevS7ZUDcVjDlxbJuLYe5YOO8+iKsQuBhysmDNt240ZsH7BUPtVm67czm0LsZEjHCzdjIc2wfttcyiJJR4EBERERERCaVBj0NmOrQd5pYinOnNf8+73BjLRe9Ct1sCd/+NP0DWwdAsGUioAlWT3QjKsmzx+1ChJjTr73UkZZISDyIiIiIiIqHU7nL/jkvpCQ06w8znoetNEBWg2QDpEyAmAZr0Dcz1ilKnTdkeqXn8IKz6CjqPhuhYr6MpkzTVQkREREREJBwZ46oedq92jSYDwVr3JrpJ38Au3zib2q1hVzrk5oTmfoE29W+uQWaHq72OpMxS4kFERERERCRctbsMKjeAmc8G5nqZq2DfxuBPsyioTlvIzYI960J3T4Bti91Y0n2bSn6NxR/AjGeg+20nppBIsSnxICIiIiIiEq6iY6HHbbBuKmxfWvrrpU9w25aDS38tf9Vp47Y7l4fungDT/+V+3zcvdyNLi2vrQvj8F9DoXBj0l4CHV56U+R4P2dnZZGRkcOzYMa9DCYqEhASSkpKIjdVaIhERERGRcqnrjTDtH67Xw2WlrHxInwD12kPVhgEJzS+1WgImtA0mD++GlV9C036waRa8PQJuGOeaXfrjUCa8Nxoq1IIrX1dvh1Iq84mHjIwMKleuTOPGjTHBHgUTYtZadu/eTUZGBk2aNPE6HBERERER8UJideh0Dcx/AwY+ApXqlOw6R/bA5lluPGcoxVWAGk1CW/Gw5APIy3aVCge2wrsj4d1RcO2HEJt49nNzs+GDG+HILrh5IlSqHZKQI1mZX2px7NgxatasGXFJBwBjDDVr1ozYag4REREREfFTjztcn4Q5r5T8Gqu/AZsX2mUW+WqHcLKFtbDgTTcRpG47aHEBXP6iGyP64c1FN7mc+CBs/B6G/RcadApJyJGuzCcegIhMOuSL5N9NRERERET8VKu5SxjMeRmyS/jBZPoEqFjHvSEPtTptYPdayDke/HttWwg7lkLna0/sa38FDP0nrBoPn98NeXmFn7vgLZj9IvT6BXS4KvixlhMRkXgIR40bN2bXrl1ehyEiIiIiIpGi112u/H/J+8U/Nzcb1kyGlhdClAdvA+u0AZsLu1YH/17z34SYBEi94uT93W+F/g/Confg64dcZURBGXNh3D2uL8TAPwY/znJEiQcREREREZGyoPF5ULc9zHju9DfNRdk0A47vh5ZDghNbUfInWwS7wWT2UVjyIbQZBonVTn+8z/1u2crMZ2H6Eyf2H9wB/7sWKteHK16F6DLfDjGsKPEQAJdddhldu3alXbt2jBkz5rTH//SnP9GqVSt69+7NqFGjeOIJ9we+cOFCevbsSYcOHbj88svZu3dvqEMXEREREZGywhjodSdkroB1U4p37qoJEB3nPs33Qs0WEBUT/AaTK75wCZYu1xX+uDGu4WSHkfDtn93SlZzj8P51cGw/jHwHKtQIbozlUESlcf74xTKWbz0Q0Gu2bVCFRy5pd9Zjxo4dS40aNTh69CjdunVjxIgRPz02Z84cPvroIxYtWkR2djZdunSha9euAFx//fX897//pW/fvjz88MP88Y9/5Kmnngpo/CIiIiIiEkFSR8A3j7iqh2YD/D8vfQI06QPxlYIX29nExEHN5sFvMLngTajWCBr1PvMxUVFw6TMu0fDlfbDkIzft48rXoF5qcOMrp/yqeDDGDDbGrDLGrDHGPFDI4ynGmCnGmAXGmMXGmKG+/RcYY+YZY5b4tgMKnDPVd82Fvq8SzoTx3tNPP03Hjh3p2bMnmzdvZvXqE+uWfvjhBy699FISEhKoXLkyl1xyCQD79+9n37599O3bF4AbbriB7777zpP4RURERESkjIiJd70K1nwDmav8O2fXatiz1ptpFgXVbh3cxMOe9bD+O+h8XdF9LKJj4cpXodE5sOlH6P1raHd58GIr54qseDDGRAPPAhcAGcAcY8zn1tqCNTIPAe9ba583xrQFxgONgV3AJdbarcaYVGAi0LDAeaOttXMD86tQZGVCMEydOpVJkyYxY8YMKlSoQL9+/TT+UkREREREgiftZvjuCZj5HFzyn6KPT5/gti0HBTeuotRpC8s/g6wjEFch8Ndf+A5goNMo/46PTYRr/gdrp0DriwIfj/zEn4qH7sAaa+06a20W8B5w6SnHWKCK7/uqwFYAa+0Ca+1W3/5lQKIxJr70YYeP/fv3U716dSpUqMDKlSuZOXPmSY+fe+65fPHFFxw7doxDhw4xbtw4AKpWrUr16tWZPn06AG+++eZP1Q8iIiIiIiJnVLEWdBwJi96Dw7uLPn7VBKjTDqqlBD+2s6nTGrCwy89KjeLIy3WJh+bnQ9Uk/8+Lrwxth0FUdOBjkp/4k3hoCGwu8HMGJ1ctADwKXGuMycBVO9xdyHVGAPOttQUHt77qW2bxB2OM8T/s8DF48GBycnJo06YNDzzwAD179jzp8W7dujFs2DA6dOjAkCFDaN++PVWrVgXg9ddf5/7776dDhw4sXLiQhx9+2ItfQUREREREypqed0LOMZjyOGxfAlmHCz/u6F430aKVx8sswFU8QHCWW6ybAgcyoPO1gb+2lFqgmkuOAl6z1v7LGNMLeNMYk2qtzQMwxrQD/g5cWOCc0dbaLcaYysBHwHXAG6de2BhzG3AbQEqKxxm6QsTHx/PVV1+dtn/Dhg0/fX/ffffx6KOPcuTIEfr06fNTc8lOnTqdViEhIiIiIiJSpDqtoe2lMPcV9wVQuQHUbOa+ajRzzRz3rAWb631/B4DqTSA6PjiTLRa8BYk1oNXQwF9bSs2fxMMWILnAz0m+fQXdAgwGsNbOMMYkALWAncaYJOAT4Hpr7dr8E6y1W3zbg8aYd3BLOk5LPFhrxwBjANLS0oo5rDY83HbbbSxfvpxjx45xww030KVLF69DEhERERGRsm7EWDhvuUsu7F4Du9e6r+Wfw9E9J46rUAsadvUuznzRMVCrJexcGdjrHtkDK7+EtFtc800JO/4kHuYALYwxTXAJh5HANaccswk4H3jNGNMGSAAyjTHVgC+BB6y1P+QfbIyJAapZa3cZY2KBi4FJpf1lwtU777zjdQgiIiIiIhJpomOgfgf3daqje2H3OpeQqNEkfHoY1GkDG76HvLyiJ0/4a/H7kJulZRZhrMh/09baHOAXuIkUK3DTK5YZYx4zxgzzHXYvcKsxZhHwLnCjtdb6zmsOPHzK2Mx4YKIxZjGwEJfQeCnAv5uIiIiIiEj5lFgdkrpCx6shubvX0ZzQtC8c3AofXO+mW5SWtbDgTWjQGeqllv56EhR+9Xiw1o7HNY0suO/hAt8vB84t5Lw/A38+w2XDoNZHREREREREQqbTaDi2HyY+CPsvglHvQeW6Jb/etoWwYylc9K+AhSiBF6DaFhEREREREZEiGAO97oKRb0PmSnj5fNhRimaT89+EmARIvSJwMUrAKfEgIiIiIiIiodX6IrhpPORmw9hBsGZy8a+RfRSWfAhthkFitYCHKIGjxEMAPP7447Rr144OHTrQqVMnZs2aRb9+/WjVqhUdO3bk3HPPZdWqVfTu3fuk0ZsffPABgweHwVgbERERERGRUGvQGW6dDNVS4O0rYe7Y4p2/Yhwc3w9drgtOfBIwfvV4kDObMWMG48aNY/78+cTHx7Nr1y6ysrIAePvtt0lLS2PMmDHcf//9vPDCC1x55ZX079+fnJwcfv/73zNhwgSPfwMRERERERGPVE2CmyfAhzfDuHvcONAL/uTfxIsFb0C1RtCod/DjlFJR4qGUtm3bRq1atYiPd/Nia9Wqddoxffr04amnniI1NZVLLrmEv//97xw+fJjrr7+eZs2ahTpkERERERGR8BFfGUa+CxN/BzOegb0bYPhLEFfBPW4t5BxzTSmPHXDbA1tg/XfQ/6HAjeWUoImsxMNXD8D2JYG9Zr32MORvZ3z4wgsv5LHHHqNly5YMHDiQq6++mr59+550zBdffEH79u0BeOSRR+jSpQtxcXHMnTs3sLGKiIiIiIiURdExMPSfUKMZTHgA/tsVYhNOJBrysgs5Jw46jQp9rFJskZV48EClSpWYN28e06dPZ8qUKVx99dX87W8uUTF69GgSExNp3Lgx//3vfwGoWLEiV199NZUqVfqpSkJERERERESAnrdDjaYw71WITYSEqhBfxW0TqkBCtRM/V23olmpI2IusxMNZKhOCKTo6mn79+tGvXz/at2/P66+/Dpzo8XCqqKgoolQOJCIiIiIicrqWF7oviRh691tKq1atYvXq1T/9vHDhQho1auRhRCIiIiIiIiLhI7IqHjxw6NAh7r77bvbt20dMTAzNmzdnzJgxXHHFFV6HJiIiIiIiIuI5JR5KqWvXrvz444+n7Z86deoZz3n00UeDF5CIiIiIiIhIGNFSCxEREREREREJGiUeRERERERERCRolHgQERERERERkaCJiMSDtdbrEIImkn83ERERERERiXxlPvGQkJDA7t27I/INurWW3bt3k5CQ4HUoIiIiIiIiIiVS5qdaJCUlkZGRQWZmptehBEVCQgJJSUlehyEiIiIiIiJSImU+8RAbG0uTJk28DkNEREREREREClHml1qIiIiIiIiISPhS4kFEREREREREgkaJBxEREREREREJGlOWpkEYYzKBjV7HUUy1gF1eByESQPqblkijv2mJNPqblkijv2mJNJH6N93IWlu7sAfKVOKhLDLGzLXWpnkdh0ig6G9aIo3+piXS6G9aIo3+piXSlMe/aS21EBEREREREZGgUeJBRERERERERIJGiYfgG+N1ACIBpr9piTT6m5ZIo79piTT6m5ZIU+7+ptXjQURERERERESCRhUPIiIiIiIiIhI0SjwEkTFmsDFmlTFmjTHmAa/jESkuY0yyMWaKMWa5MWaZMeb/fPtrGGO+Mcas9m2rex2riL+MMdHGmAXGmHG+n5sYY2b5nqv/Z4yJ8zpGEX8ZY6oZYz40xqw0xqwwxvTSc7SUZcaYe3yvOZYaY941xiToeVrKEmPMWGPMTmPM0gL7Cn1eNs7Tvr/txcaYLt5FHlxKPASJMSYaeBYYArQFRhlj2noblUix5QD3WmvbAj2Bu3x/xw8Ak621LYDJvp9Fyor/A1YU+PnvwJPW2ubAXuAWT6ISKZn/ABOsta2Bjri/bT1HS5lkjGkI/BJIs9amAtHASPQ8LWXLa8DgU/ad6Xl5CNDC93Ub8HyIYgw5JR6Cpzuwxlq7zlqbBbwHXOpxTCLFYq3dZq2d7/v+IO4FbUPc3/LrvsNeBy7zJECRYjLGJAEXAS/7fjbAAOBD3yH6e5YywxhTFegDvAJgrc2y1u5Dz9FStsUAicaYGKACsA09T0sZYq39Dthzyu4zPS9fCrxhnZlANWNM/ZAEGmJKPARPQ2BzgZ8zfPtEyiRjTGOgMzALqGut3eZ7aDtQ16u4RIrpKeA3QJ7v55rAPmttju9nPVdLWdIEyARe9S0fetkYUxE9R0sZZa3dAjwBbMIlHPYD89DztJR9Z3peLjfvGZV4EJEiGWMqAR8Bv7LWHij4mHWjcTQeR8KeMeZiYKe1dp7XsYgESAzQBXjeWtsZOMwpyyr0HC1liW/d+6W4pFoDoCKnl6yLlGnl9XlZiYfg2QIkF/g5ybdPpEwxxsTikg5vW2s/9u3ekV8G5tvu9Co+kWI4FxhmjNmAW/42ALc+vpqvpBf0XC1lSwaQYa2d5fv5Q1wiQs/RUlYNBNZbazOttdnAx7jnbj1PS1l3puflcvOeUYmH4JkDtPB14Y3DNcb53OOYRIrFt/79FWCFtfbfBR76HLjB9/0NwGehjk2kuKy1v7PWJllrG+Oek7+11o4GpgBX+A7T37OUGdba7cBmY0wr367zgeXoOVrKrk1AT2NMBd9rkPy/aT1PS1l3puflz4HrfdMtegL7CyzJiCjGVXpIMBhjhuLWE0cDY621j3sbkUjxGGN6A9OBJZxYE/97XJ+H94EUYCNwlbX21CY6ImHLGNMPuM9ae7ExpimuAqIGsAC41lp73MPwRPxmjOmEa5YaB6wDbsJ9sKTnaCmTjDF/BK7GTdZaAPwMt+Zdz9NSJhhj3gX6AbWAHcAjwKcU8rzsS7A9g1tSdAS4yVo714Owg06JBxEREREREREJGi21EBEREREREZGgUeJBRERERERERIJGiQcRERERERERCRolHkREREREREQkaJR4EBEREREREZGgUeJBRERERERERIJGiQcRERERERERCRolHkREREREREQkaP4f9Ju9ffBPhIcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1296x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(18,6))\n",
    "plt.plot(np.nancumprod(total_ret + 1), label='algo')\n",
    "plt.plot(np.nancumprod(index_ret + 1), label=index_comp)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48490d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "def metrics_calc(cumret):\n",
    "    total_return = (cumret[-1] - cumret[0])/cumret[0]\n",
    "    apr = (1+total_return)**(252/len(cumret)) - 1\n",
    "    rets = pd.DataFrame(cumret).pct_change()\n",
    "    sharpe = np.sqrt(252) * np.nanmean(rets) / np.nanstd(rets)\n",
    "    \n",
    "    # maxdd and maxddd\n",
    "    highwatermark=np.zeros(cumret.shape)\n",
    "    drawdown=np.zeros(cumret.shape)\n",
    "    drawdownduration=np.zeros(cumret.shape)\n",
    "    for t in np.arange(1, cumret.shape[0]):\n",
    "        highwatermark[t]=np.maximum(highwatermark[t-1], cumret[t])\n",
    "        drawdown[t]=cumret[t]/highwatermark[t]-1\n",
    "        if drawdown[t]==0:\n",
    "            drawdownduration[t]=0\n",
    "        else:\n",
    "            drawdownduration[t]=drawdownduration[t-1]+1\n",
    "    maxDD=np.min(drawdown)\n",
    "    maxDDD=np.max(drawdownduration)\n",
    "    \n",
    "    return total_return, apr, sharpe, maxDD, maxDDD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "084a9877",
   "metadata": {},
   "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>Total return</th>\n",
       "      <th>APR</th>\n",
       "      <th>Sharpe</th>\n",
       "      <th>Max Drawdown</th>\n",
       "      <th>Max Drawdown Duration</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Algo</th>\n",
       "      <td>0.007799</td>\n",
       "      <td>0.019378</td>\n",
       "      <td>0.577382</td>\n",
       "      <td>-0.012455</td>\n",
       "      <td>63.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SPY</th>\n",
       "      <td>-0.132629</td>\n",
       "      <td>-0.296394</td>\n",
       "      <td>-1.330918</td>\n",
       "      <td>-0.170587</td>\n",
       "      <td>95.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Total return       APR    Sharpe Max Drawdown Max Drawdown Duration\n",
       "Algo     0.007799  0.019378  0.577382    -0.012455                  63.0\n",
       "SPY     -0.132629 -0.296394 -1.330918    -0.170587                  95.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metrics = pd.DataFrame(index=['Algo', index_comp], columns=['Total return', 'APR', 'Sharpe', 'Max Drawdown', 'Max Drawdown Duration'])\n",
    "metrics.loc['Algo'] = metrics_calc(np.nancumprod(total_ret + 1))\n",
    "metrics.loc[index_comp] = metrics_calc(np.nancumprod(index_ret + 1))\n",
    "metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18b70033",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}