Portfolio_selection_Random_Forest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 Portfolio selection using Random Forest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Goal: Develop and train a predictive model that can reliably predict future changes in the stock prices, and then use it for construction of stock portfolio.  The objective is to 'beat the market'. <br>\n",
    "Data: stocks obtained from Wikipedia by scrapping. Stock prices and fundamentals are obtained using API 'yahoofinancials' (https://pypi.org/project/yahoofinancials/). S&P100 index and stock constituents are considered. Planned expansion - to S&P500, and then to Russell 3000 stocks. <br>\n",
    "Training is performed using 65% of stocks data. Testing (and back-testing) performed for the remaining 35% of the time period. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ++++++++++++++++++++++++++++++++++++++++++++ <br>\n",
    " By: Yuriy Podvysotskiy <br>\n",
    " Date: November 4, 2019 <br>\n",
    " ++++++++++++++++++++++++++++++++++++++++++++"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages:\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import yahoofinancials as yf\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from bs4 import BeautifulSoup\n",
    "import requests\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>Ticker</th>\n",
       "      <th>Name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>AAPL</td>\n",
       "      <td>Apple Inc.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>ABBV</td>\n",
       "      <td>AbbVie Inc.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>ABT</td>\n",
       "      <td>Abbott Laboratories</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>ACN</td>\n",
       "      <td>Accenture</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>ADBE</td>\n",
       "      <td>Adobe Inc.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Ticker                 Name\n",
       "0   AAPL           Apple Inc.\n",
       "1   ABBV          AbbVie Inc.\n",
       "2    ABT  Abbott Laboratories\n",
       "3    ACN            Accenture\n",
       "4   ADBE           Adobe Inc."
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Scrap stocks names in S&P 100 index:\n",
    "url = 'https://en.wikipedia.org/wiki/S%26P_100'\n",
    "textnames = requests.get(url)\n",
    "\n",
    "soup = BeautifulSoup(textnames.text, 'lxml')\n",
    "tables = soup.select('table')\n",
    "tables_data = tables[2].select('tr')\n",
    "\n",
    "col_names = [\"Ticker\", \"Name\"]\n",
    "stocks = pd.DataFrame(columns=col_names)\n",
    "\n",
    "for i in range(1, len(tables_data)):\n",
    "    a1 = tables_data[i].select('td')[0].getText().replace(\"\\n\",\"\")\n",
    "    a2 = tables_data[i].select('td')[1].getText().replace(\"\\n\",\"\")\n",
    "    \n",
    "    stocks = stocks.append(pd.DataFrame([[a1, a2]], columns = col_names))\n",
    "\n",
    "stocks = stocks.reset_index(drop = True)\n",
    "stocks.head(5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "issue with:  BRK.B\n"
     ]
    }
   ],
   "source": [
    "# Scrap stock prices:\n",
    "yahoo_financials_stocks = yf.YahooFinancials(stocks['Ticker'])\n",
    "\n",
    "# 1 Read and get prices:\n",
    "weekly_stock_prices = yahoo_financials_stocks.get_historical_price_data('2015-01-01', '2019-11-01', 'weekly')\n",
    "\n",
    "## A. Some exploration:\n",
    "weekly_stock_prices['AAPL'].keys()\n",
    "pd.DataFrame(weekly_stock_prices['AAPL']['prices'])\n",
    "\n",
    "## B. Make up the panel:\n",
    "my_stocks = list(weekly_stock_prices.keys())\n",
    "df_prices = pd.DataFrame(columns = ['date', 'ticker', 'formatted_date', 'adjclose', 'volume'])\n",
    "\n",
    "for tick in my_stocks:\n",
    "    # tick = my_stocks[0]\n",
    "    try: \n",
    "        res = pd.DataFrame(weekly_stock_prices[tick]['prices'])\n",
    "    except:\n",
    "        print(\"issue with: \", tick)\n",
    "        continue \n",
    "    # res.columns.T\n",
    "    res = res[['date', 'formatted_date', 'adjclose', 'volume']]\n",
    "    res['ticker'] = tick\n",
    "    df_prices = df_prices.append(res, sort = False)\n",
    "\n",
    "df_prices = df_prices.fillna(method = 'ffill')\n",
    "df_prices.to_csv(\"SP100_prices.csv\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doing Income for : AAPL\n",
      "Doing Balance for : AAPL\n",
      "Doing Income for : ABBV\n",
      "Doing Balance for : ABBV\n",
      "Doing Income for : ABT\n",
      "Doing Balance for : ABT\n",
      "Doing Income for : ACN\n",
      "Doing Balance for : ACN\n",
      "Doing Income for : ADBE\n",
      "Doing Balance for : ADBE\n",
      "Doing Income for : AGN\n",
      "Doing Balance for : AGN\n",
      "Doing Income for : AIG\n",
      "Doing Balance for : AIG\n",
      "Doing Income for : ALL\n",
      "Doing Balance for : ALL\n",
      "Doing Income for : AMGN\n",
      "Doing Balance for : AMGN\n",
      "Doing Income for : AMZN\n",
      "Doing Balance for : AMZN\n",
      "Doing Income for : AXP\n",
      "Doing Balance for : AXP\n",
      "Doing Income for : BA\n",
      "Doing Balance for : BA\n",
      "Doing Income for : BAC\n",
      "Doing Balance for : BAC\n",
      "Doing Income for : BIIB\n",
      "Doing Balance for : BIIB\n",
      "Doing Income for : BK\n",
      "Doing Balance for : BK\n",
      "Doing Income for : BKNG\n",
      "Doing Balance for : BKNG\n",
      "Doing Income for : BLK\n",
      "Doing Balance for : BLK\n",
      "Doing Income for : BMY\n",
      "Doing Balance for : BMY\n",
      "Doing Income for : BRK.B\n",
      "Doing Balance for : BRK.B\n",
      "Doing Income for : C\n",
      "Doing Balance for : C\n",
      "Doing Income for : CAT\n",
      "Doing Balance for : CAT\n",
      "Doing Income for : CELG\n",
      "Doing Balance for : CELG\n",
      "Doing Income for : CHTR\n",
      "Doing Balance for : CHTR\n",
      "Doing Income for : CL\n",
      "Doing Balance for : CL\n",
      "Doing Income for : CMCSA\n",
      "Doing Balance for : CMCSA\n",
      "Doing Income for : COF\n",
      "Doing Balance for : COF\n",
      "Doing Income for : COP\n",
      "Doing Balance for : COP\n",
      "Doing Income for : COST\n",
      "Doing Balance for : COST\n",
      "Doing Income for : CSCO\n",
      "Doing Balance for : CSCO\n",
      "Doing Income for : CVS\n",
      "Doing Balance for : CVS\n",
      "Doing Income for : CVX\n",
      "Doing Balance for : CVX\n",
      "Doing Income for : DD\n",
      "Doing Balance for : DD\n",
      "Doing Income for : DHR\n",
      "Doing Balance for : DHR\n",
      "Doing Income for : DIS\n",
      "Doing Balance for : DIS\n",
      "Doing Income for : DOW\n",
      "Doing Balance for : DOW\n",
      "issue with Balance for:  DOW\n",
      "issue with Balance for:  DOW\n",
      "issue with Balance for:  DOW\n",
      "Doing Income for : DUK\n",
      "Doing Balance for : DUK\n",
      "Doing Income for : EMR\n",
      "Doing Balance for : EMR\n",
      "Doing Income for : EXC\n",
      "Doing Balance for : EXC\n",
      "Doing Income for : F\n",
      "Doing Balance for : F\n",
      "Doing Income for : FB\n",
      "Doing Balance for : FB\n",
      "Doing Income for : FDX\n",
      "Doing Balance for : FDX\n",
      "Doing Income for : GD\n",
      "Doing Balance for : GD\n",
      "Doing Income for : GE\n",
      "Doing Balance for : GE\n",
      "Doing Income for : GILD\n",
      "Doing Balance for : GILD\n",
      "Doing Income for : GM\n",
      "Doing Balance for : GM\n",
      "Doing Income for : GOOG\n",
      "Doing Balance for : GOOG\n",
      "Doing Income for : GOOGL\n",
      "Doing Balance for : GOOGL\n",
      "Doing Income for : GS\n",
      "Doing Balance for : GS\n",
      "Doing Income for : HD\n",
      "Doing Balance for : HD\n",
      "Doing Income for : HON\n",
      "Doing Balance for : HON\n",
      "Doing Income for : IBM\n",
      "Doing Balance for : IBM\n",
      "Doing Income for : INTC\n",
      "Doing Balance for : INTC\n",
      "Doing Income for : JNJ\n",
      "Doing Balance for : JNJ\n",
      "Doing Income for : JPM\n",
      "Doing Balance for : JPM\n",
      "Doing Income for : KHC\n",
      "Doing Balance for : KHC\n",
      "Doing Income for : KMI\n",
      "Doing Balance for : KMI\n",
      "Doing Income for : KO\n",
      "Doing Balance for : KO\n",
      "Doing Income for : LLY\n",
      "Doing Balance for : LLY\n",
      "Doing Income for : LMT\n",
      "Doing Balance for : LMT\n",
      "Doing Income for : LOW\n",
      "Doing Balance for : LOW\n",
      "Doing Income for : MA\n",
      "Doing Balance for : MA\n",
      "Doing Income for : MCD\n",
      "Doing Balance for : MCD\n",
      "Doing Income for : MDLZ\n",
      "Doing Balance for : MDLZ\n",
      "Doing Income for : MDT\n",
      "Doing Balance for : MDT\n",
      "Doing Income for : MET\n",
      "Doing Balance for : MET\n",
      "Doing Income for : MMM\n",
      "Doing Balance for : MMM\n",
      "Doing Income for : MO\n",
      "Doing Balance for : MO\n",
      "Doing Income for : MRK\n",
      "Doing Balance for : MRK\n",
      "Doing Income for : MS\n",
      "Doing Balance for : MS\n",
      "Doing Income for : MSFT\n",
      "Doing Balance for : MSFT\n",
      "Doing Income for : NEE\n",
      "Doing Balance for : NEE\n",
      "Doing Income for : NFLX\n",
      "Doing Balance for : NFLX\n",
      "Doing Income for : NKE\n",
      "Doing Balance for : NKE\n",
      "Doing Income for : NVDA\n",
      "Doing Balance for : NVDA\n",
      "Doing Income for : ORCL\n",
      "Doing Balance for : ORCL\n",
      "Doing Income for : OXY\n",
      "Doing Balance for : OXY\n",
      "Doing Income for : PEP\n",
      "Doing Balance for : PEP\n",
      "Doing Income for : PFE\n",
      "Doing Balance for : PFE\n",
      "Doing Income for : PG\n",
      "Doing Balance for : PG\n",
      "Doing Income for : PM\n",
      "Doing Balance for : PM\n",
      "Doing Income for : PYPL\n",
      "Doing Balance for : PYPL\n",
      "Doing Income for : QCOM\n",
      "Doing Balance for : QCOM\n",
      "Doing Income for : RTN\n",
      "Doing Balance for : RTN\n",
      "Doing Income for : SBUX\n",
      "Doing Balance for : SBUX\n",
      "Doing Income for : SLB\n",
      "Doing Balance for : SLB\n",
      "Doing Income for : SO\n",
      "Doing Balance for : SO\n",
      "Doing Income for : SPG\n",
      "Doing Balance for : SPG\n",
      "Doing Income for : T\n",
      "Doing Balance for : T\n",
      "Doing Income for : TGT\n",
      "Doing Balance for : TGT\n",
      "Doing Income for : TXN\n",
      "Doing Balance for : TXN\n",
      "Doing Income for : UNH\n",
      "Doing Balance for : UNH\n",
      "Doing Income for : UNP\n",
      "Doing Balance for : UNP\n",
      "Doing Income for : UPS\n",
      "Doing Balance for : UPS\n",
      "Doing Income for : USB\n",
      "Doing Balance for : USB\n",
      "Doing Income for : UTX\n",
      "Doing Balance for : UTX\n",
      "Doing Income for : V\n",
      "Doing Balance for : V\n",
      "Doing Income for : VZ\n",
      "Doing Balance for : VZ\n",
      "Doing Income for : WBA\n",
      "Doing Balance for : WBA\n",
      "Doing Income for : WFC\n",
      "Doing Balance for : WFC\n",
      "Doing Income for : WMT\n",
      "Doing Balance for : WMT\n",
      "Doing Income for : XOM\n",
      "Doing Balance for : XOM\n"
     ]
    }
   ],
   "source": [
    "# Scrap financials:\n",
    "df_income = pd.DataFrame(columns = ['date', 'ticker', 'totalRevenue', 'ebit', 'netIncome'])\n",
    "df_balance = pd.DataFrame(columns = ['date', 'ticker', 'totalAssets', 'totalLiab'])\n",
    "\n",
    "# Progress bar:\n",
    "# https://stackoverflow.com/questions/3002085/python-to-print-out-status-bar-and-percentage\n",
    "\n",
    "for tick in stocks['Ticker']:\n",
    "\n",
    "    yahoo_financials_stocks = yf.YahooFinancials(tick)\n",
    "    try:\n",
    "        fin_data = yahoo_financials_stocks.get_financial_stmts('annual', ['income', 'balance'])\n",
    "    except:\n",
    "        print(\"issue with: \", tick)\n",
    "        continue\n",
    "    \n",
    "    # 2.1 Income Statement:\n",
    "    print(\"Doing Income for :\", tick)\n",
    "    n_dates = len(fin_data['incomeStatementHistory'][tick])\n",
    "    \n",
    "    for i_date in range(n_dates):\n",
    "        try:\n",
    "            date0 = list(fin_data['incomeStatementHistory'][tick][i_date].keys())[0]\n",
    "            ni0 = fin_data['incomeStatementHistory'][tick][i_date][date0]['netIncome']\n",
    "            ebit0 = fin_data['incomeStatementHistory'][tick][i_date][date0]['ebit']\n",
    "            rev0 = fin_data['incomeStatementHistory'][tick][i_date][date0]['totalRevenue']\n",
    "        except:\n",
    "            print(\"issue with income for: \", tick)\n",
    "            continue\n",
    "        df_income= df_income.append({\n",
    "                'date': date0,\n",
    "                'ticker': tick,\n",
    "                'totalRevenue': rev0,\n",
    "                'ebit': ebit0,\n",
    "                'netIncome': ni0\n",
    "                }, ignore_index = True)\n",
    "    \n",
    "    # 2.2 Balance sheet:\n",
    "    print(\"Doing Balance for :\", tick)\n",
    "    n_dates = len(fin_data['balanceSheetHistory'][tick])\n",
    "    \n",
    "    for i_date in range(n_dates):\n",
    "        try:\n",
    "            date0 = list(fin_data['balanceSheetHistory'][tick][i_date].keys())[0]\n",
    "            assets0 = fin_data['balanceSheetHistory'][tick][i_date][date0]['totalAssets']\n",
    "            liabs0 = fin_data['balanceSheetHistory'][tick][i_date][date0]['totalLiab']\n",
    "        except:\n",
    "            print(\"issue with Balance for: \", tick)\n",
    "            continue\n",
    "        df_balance= df_balance.append({\n",
    "                'date': date0,\n",
    "                'ticker': tick,\n",
    "                'totalAssets': assets0,\n",
    "                'totalLiab': liabs0\n",
    "                }, ignore_index = True)\n",
    "\n",
    "    \n",
    "df_fs = df_income.merge(df_balance, how = 'outer', on = ['date', 'ticker']) \n",
    "df_fs = df_fs.fillna(method = 'ffill')\n",
    "df_fs['ebitm'] = df_fs['ebit']/df_fs['totalRevenue']\n",
    "df_fs['nim'] = df_fs['netIncome']/df_fs['totalRevenue']\n",
    "df_fs['date'] = pd.to_datetime(df_fs['date'])\n",
    "df_fs = df_fs.sort_values(by = ['ticker', 'date'])\n",
    "df_fs['revg'] = df_fs.groupby(['ticker'])['totalRevenue'].pct_change()\n",
    "\n",
    "df_fs = df_fs.rename(columns={\"totalAssets\":\"assets\"})\n",
    "df_fs['lever'] = df_fs['totalLiab']/df_fs['assets']\n",
    "df_fs = df_fs[['date', 'ticker', 'assets', 'ebitm', 'nim', 'lever', 'revg']]\n",
    "df_fs = df_fs.reset_index(drop = True)\n",
    "df_fs['year'] = df_fs['date'].map(lambda x: x.year)\n",
    "\n",
    "df_fs.to_csv(\"SP100_financials.csv\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge data (prices + financials) together: \n",
    "df_prices = df_prices.drop(columns = ['date'])\n",
    "df_prices = df_prices.rename(columns = {\"formatted_date\" : \"date\"})\n",
    "df_prices['date'] = pd.to_datetime(df_prices['date'])\n",
    "df_prices['year'] = df_prices['date'].map(lambda x: x.year)\n",
    "\n",
    "\n",
    "df_prices['year-1'] = df_prices['year'] - 1\n",
    "df_panel = df_prices.merge(df_fs, how = 'inner', left_on=['ticker', 'year-1'], \n",
    "                           right_on = ['ticker', 'year'])\n",
    "df_panel = df_panel.drop(columns = ['date_y', 'year_y', 'year-1'])\n",
    "df_panel = df_panel.rename(columns = {\"date_x\" : \"date\",\n",
    "                                      \"year_x\": \"year\"})\n",
    "df_panel = df_panel.sort_values(by = ['ticker', 'date'])\n",
    "df_panel['return'] = df_panel.groupby(['ticker'])['adjclose'].pct_change()\n",
    "df_panel.isnull().sum()\n",
    "df_panel = df_panel.dropna()\n",
    "\n",
    "# Save: \n",
    "df_panel.to_csv(\"SP100_panel.csv\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get market index:\n",
    "SP100 = yf.YahooFinancials('^OEX').get_historical_price_data('2015-01-01', '2019-11-01', 'weekly')\n",
    "df_mkt = pd.DataFrame(SP100['^OEX']['prices'])\n",
    "df_mkt = df_mkt[['date', 'formatted_date', 'adjclose']]\n",
    "df_mkt = df_mkt.drop(columns = ['date'])\n",
    "df_mkt = df_mkt.rename(columns = {\"formatted_date\" : \"date\"})\n",
    "df_mkt['date'] = pd.to_datetime(df_mkt['date'])\n",
    "df_mkt = df_mkt.sort_values(by = ['date'])\n",
    "df_mkt = df_mkt.reset_index(drop = True)\n",
    "df_mkt['mkt_ret'] = df_mkt['adjclose'].pct_change()\n",
    "df_mkt = df_mkt[['date', 'mkt_ret']]\n",
    "df_mkt.to_csv(\"SP100.csv\", index = False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add more features & prepare the data:\n",
    "df_panel['ret5w'] = df_panel.groupby(['ticker'])['return'].rolling(5).mean().values\n",
    "df_panel['vol5w'] = df_panel.groupby(['ticker'])['return'].rolling(5).std().values\n",
    "\n",
    "df_panel.isnull().sum()\n",
    "df_panel = df_panel.dropna()\n",
    "df_panel = df_panel.reset_index(drop = True)\n",
    "df_panel = df_panel.sort_values(by = ['ticker', 'date'])\n",
    "df_panel['fret'] = df_panel.groupby(['ticker'])['return'].shift(-1)\n",
    "\n",
    "df_panel.isnull().sum()\n",
    "df_panel = df_panel.dropna()\n",
    "df_panel = df_panel.reset_index(drop = True)\n",
    "\n",
    "df_X = df_panel\n",
    "df_X = df_X.drop(columns = ['ticker', 'year', 'date', 'fret'])\n",
    "\n",
    "df_y = df_panel['fret']\n",
    "df_y = pd.DataFrame(df_y)\n",
    "\n",
    "df_y.loc[:,'chg'] = 0\n",
    "df_y.loc[df_y['fret'] > 0,'chg'] = 1\n",
    "df_y = df_y.drop(columns = ['fret'])\n",
    "df_y = df_y.reset_index(drop = True)\n",
    "len(df_panel)\n",
    "len(df_X)\n",
    "len(df_y)\n",
    "\n",
    "n_size = len(df_panel)\n",
    "n_test = np.floor(n_size*.6)\n",
    "\n",
    "X_train = df_X.loc[0:n_test,:]\n",
    "X_test = df_X.loc[n_test:n_size, :]\n",
    "\n",
    "y_train = df_y.loc[0:n_test].values.reshape(-1,1)\n",
    "y_test = df_y.loc[n_test:n_size].values.reshape(-1,1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 4 folds for each of 84 candidates, totalling 336 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Using backend LokyBackend with 8 concurrent workers.\n",
      "[Parallel(n_jobs=-1)]: Done  88 tasks      | elapsed:    3.7s\n",
      "[Parallel(n_jobs=-1)]: Done 283 tasks      | elapsed:   18.5s\n",
      "[Parallel(n_jobs=-1)]: Done 336 out of 336 | elapsed:   25.3s finished\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>0</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>Accuracy</td>\n",
       "      <td>0.549623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Recall</td>\n",
       "      <td>0.863848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>Precision</td>\n",
       "      <td>0.567461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>F1</td>\n",
       "      <td>0.684968</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  0\n",
       "Accuracy   0.549623\n",
       "Recall     0.863848\n",
       "Precision  0.567461\n",
       "F1         0.684968"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train and test the model:\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import *\n",
    "\n",
    "#X_train.columns.T\n",
    "\n",
    "params = {'n_estimators':[10, 20, 30, 40, 50, 100], \n",
    "          'max_depth':[4, 5, 6, 7, 8, 9, 10],\n",
    "          'criterion':['gini', 'entropy']}\n",
    "mod_forest = GridSearchCV(RandomForestClassifier(random_state=0), params, cv=4, verbose = 1, n_jobs = -1)\n",
    "mod_forest.fit(X_train, y_train)\n",
    "mod_forest_best = mod_forest.best_estimator_.fit(X_train, y_train)\n",
    "\n",
    "def get_performance(res0_0, res1_1, res1_0, res0_1):\n",
    "    accuracy = (res0_0 + res1_1)/(res0_0 + res1_0+res0_1+res1_1)\n",
    "    recall = res1_1/(res1_1+res1_0)\n",
    "    precision = res1_1/(res1_1+res0_1)\n",
    "    f1 = 2/(1/precision + 1/recall)\n",
    "    return([accuracy, recall, precision, f1])\n",
    "    \n",
    "y_Forest = mod_forest_best.predict(X_test)\n",
    "res = confusion_matrix(y_test, y_Forest)\n",
    "res_table = pd.DataFrame(get_performance(res[0,0], res[1,1], res[1,0], res[0,1]),\n",
    "                          index = ['Accuracy', 'Recall', 'Precision', 'F1'])\n",
    "res_table\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAke portfolio and Evaluate:\n",
    "# ********************** Make portfolio *****************************\n",
    "df_panel = df_panel.loc[n_test:n_size,]\n",
    "df_panel['prob_top'] = 0\n",
    "df_panel['prob'] = mod_forest_best.predict(df_panel[['adjclose', 'volume', 'assets', 'ebitm', 'nim', 'lever', 'revg',\n",
    "       'return', 'ret5w', 'vol5w']])\n",
    "\n",
    "\n",
    "def backtest_strategy(df_panel, df_mkt, prob_quantile):\n",
    "    quantiles = df_panel.groupby(['date'])['prob'].quantile(prob_quantile)\n",
    "    quantiles = pd.DataFrame(quantiles)\n",
    "    quantiles.columns = ['top_quantile']\n",
    "\n",
    "    quantiles['date'] = quantiles.index\n",
    "    quantiles = quantiles.reset_index(drop = True)\n",
    "\n",
    "    df_panel = df_panel.merge(quantiles, how = 'left', on = ['date'])\n",
    "    df_panel['top_quantile'].describe()\n",
    "\n",
    "    df_panel.loc[df_panel['prob'] >= df_panel['top_quantile'],'prob_top'] = 1\n",
    "\n",
    "    df_panel['chosen'] = df_panel['prob_top']\n",
    "    df_panel[['prob_top', 'chosen']].describe()\n",
    "\n",
    "    df_panel = df_panel.sort_values(by = ['ticker', 'date'])\n",
    "    df_panel = df_panel.reset_index(drop = True)\n",
    "\n",
    "    df_panel['chosen1'] = df_panel.groupby(['ticker'])['chosen'].shift(1)\n",
    "\n",
    "\n",
    "    ## 2.3) calculate equally-weighted returns\n",
    "    n_sums = df_panel.groupby(['date'])['chosen1'].sum()\n",
    "    n_sums = pd.DataFrame(n_sums)\n",
    "    n_sums.columns = ['n_sum']\n",
    "\n",
    "    n_sums['date'] = n_sums.index\n",
    "    n_sums = n_sums.reset_index(drop = True)\n",
    "\n",
    "    df_panel = df_panel.merge(n_sums, how = 'left', on = ['date'])\n",
    "    df_panel['n_sum'].describe()\n",
    "\n",
    "    df_panel['portf_ret'] = df_panel['return'] * df_panel['chosen1'] / df_panel['n_sum']\n",
    "    df_panel['portf_ret'] = df_panel['portf_ret'].fillna(0)\n",
    "\n",
    "    df_panel = df_panel[['date', 'portf_ret']]\n",
    "    df_panel = df_panel.dropna()\n",
    "\n",
    "    df_portf = df_panel.groupby(['date']).sum()\n",
    "    df_portf['date'] = df_portf.index\n",
    "    df_portf['date'] = pd.to_datetime(df_portf['date'])\n",
    "    df_portf = df_portf.reset_index(drop = True)\n",
    "    df_portf = df_portf.merge(df_mkt, how = 'left', on = ['date'])\n",
    "    df_portf = df_portf.fillna(method = 'ffill')\n",
    "    df_portf = df_portf.sort_values(by = ['date'])\n",
    "    df_portf = df_portf.reset_index(drop = True)\n",
    "\n",
    "    df_portf['portf_cumret'] = df_portf['portf_ret'].cumsum()\n",
    "    df_portf['mkt_cumret'] = df_portf['mkt_ret'].cumsum()\n",
    "\n",
    "    return(df_portf)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x238a2137160>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD2CAYAAADbPoDqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd1hVR9rAf0PviIAFkKKgUQFFMWqMvcQ0ExPTE2PWxPRsNm3TvrTdTXbTq5uY1UTTNdUYayyxGxF7AQERkN57udz5/hiqUi4KXMr8nuc+cM+ZM+e9lPPOvFVIKdFoNBpN98XC3AJoNBqNxrxoRaDRaDTdHK0INBqNppujFYFGo9F0c7Qi0Gg0mm6OlbkFaCkeHh7S39/f3GJoNBpNp2Lfvn2ZUkrPhs51OkXg7+9PRESEucXQaDSaToUQ4nRj57RpSKPRaLo5WhFoNBpNN0crAo1Go+nmdDofQUNUVFSQlJREaWmpuUXpNNjZ2eHj44O1tbW5RdFoNGamSyiCpKQknJ2d8ff3RwhhbnE6PFJKsrKySEpKIiAgwNziaDQaM9MlTEOlpaW4u7trJWAiQgjc3d31Dkqj0QBdRBEAWgm0EP3z0mjajj9PZbM5Kp0yQyWFZQbS8zv2oqtLmIY0Go2moxCXUchNi3YhJYwLdCersJzc4gq2PjUZG6uOufbumFJ1QiwtLRk+fDjBwcHccMMNFBcXt+j6V199td77999/n8GDB3Pbbbc1es3nn3/OQw89BMDHH3/MsmXLWi64RqNpVT7bEY+1hQVPX34RO2KyOJFaQGp+KWuPpppbtEbRiqCVsLe358CBAxw5cgQbGxs+/vhjk66TUmI0Gs9RBAsXLmT16tV89dVXJs1z3333MXfu3BbLrdFoWoe0/FJeW32c7/clcc1wL+6bOIDX54Ty5GWD8Hd3YNnOeHOL2ChaEbQB48ePJyYmBoC3336b4OBggoODeffddwGIj49n8ODBPPDAA4wYMYL58+dTUlLC8OHDue2227jvvvuIi4tj1qxZvPPOO2RnZ3PttdcSGhrKmDFjOHTo0Dn3fOmll3jzzTcBOHDgAGPGjCE0NJTZs2eTk5PTfh9eo+mmvLkuik+2xlFeaWT+eBWNd2N4Px6cHMjtY/yIOJ3Dug66K+hyPoKXfz3KseT8Vp1ziJcLL1491KSxBoOBNWvWMHPmTPbt28dnn33Gnj17kFIyevRoJk6ciJubG1FRUXz22WcsXLgQgBUrVnDgwIGaedauXcvmzZvx8PDg4YcfJiwsjJ9//plNmzYxd+7cemPPZu7cuXzwwQdMnDiRF154gZdffrlGCWk0mtbnTG4JP+0/w+1jfHlkahC9nO3qnb99jB+/Hkzmb98dYMuTk845b270jqCVqF7Rh4eH4+vry/z589m+fTuzZ8/G0dERJycnrrvuOrZt2waAn58fY8aMMWnu7du3c8cddwAwZcoUsrKyyMvLa3BsXl4eubm5TJw4EYA777yTrVu3tsIn1Gg0jbF42ykA7p8U2OBD3s7akheuHkJxeSWHEhv+3zUnXW5HYOrKvbWp9hHURUrZ6HhHR0eT525oHh3+qdGYh4pKIxHxOeSVVFBplEwc5MnKg8nMGNob7x72jV4X6OkMQGxGIdPo3V7imkSb7QiEEEuEEOlCiCPNjBslhKgUQsxpK1nMxYQJE/j5558pLi6mqKiIn376ifHjxzc41tramoqKikbnqXYab9myBQ8PD1xcXBoc6+rqipubW83O44svvqjZHWg0mguj3GDkvi/2ccunu7nvy308+HUksz/aQWZhGZcN7dPkta4O1rg72hCXUdSiexoqjXzzZwIl5ZXnLXdKXkmT59tyR/A58CHQaEyjEMIS+A+wrg3lMBsjRoxg3rx5XHzxxQDcfffdhIWFER8ff87YBQsWEBoayogRI86JFHrppZe46667CA0NxcHBgaVLlzZ536VLl3LfffdRXFxM//79+eyzz1rtM2k03Zm3N0Sz8UQ6z15xEeMCPdh2MpN/rzmBtaVg8kW9mr1+gKcTcZmFLbrn78fTeObHw+SXVHDvxAEtlrmwzMBdn+1tcoxoynxxoQgh/IFVUsrgRs4/ClQAo6rGfd/cnOHh4fLsxjTHjx9n8ODBFyxvd0P/3DQa04lOK+CK97Zx3QhvXp8zDFBm26e+P4SVpQWvXRfS7BxP/3CIDcfS2Pd/002+72PfHeDH/Wfwd3dg0+OTsLBomVn4qe8P8kPkGeJeu3KflDK8oTFmcxYLIbyB2UCzAfdCiAVCiAghRERGRkbbC6fRaDR1kFLy/M9HcLKz4unLaxdPQgjeuGGYSUoAoL+nI1lF5eQWl5s03lBpZOOJdDycbInPKmZXXFaLZd8SlcHVoX2bHGPOqKF3gb9LKZs1fEkpF0kpw6WU4Z6eDbbc1Gg0mmaJSi0gLb+ULVHpvPzrUcoMptndf4w8w5+nsnl65kX0dLQ57/sP8HQCINYEP0FpRSV//+EweSUVvHj1ENwcrPn7D4c4lJRr8v0yCspILygj2Nu1yXHmjBoKB76tin7xAK4QQhiklD+bUSaNRtNFKSmvZM5/d+Jga0lucQVlBiMZBWW8f3MYFhaCL3ef5pOtsUy9qDcvzaoffbhwSwyhPq7cGN7vgmToX6UI4jIKGenn1uTYV1Yd44fIJB6dFsRVoX3x6mHPw19H8ui3B1j/twlYWTa/jj+WonKqhno1rQjMtiOQUgZIKf2llP7A98ADWgloNJq2Yt3RVArKDGQVluPVw54HJw9g1aEUPtwcw6nMIl5ceZTcogq+2nOa9ILaaqGJ2cXEZhRx7XDvFtvnz6afmz02lhZEpxU0OS4+s4jv9iYyd6wfj04biBCCkX5uvHD1EOIyi1h5MNmk+x1NVjkLQ/o2HGVYTZvtCIQQ3wCTAA8hRBLwImANIKU0rRCPRqPRXACGSiNFZZW4Oljz/b4k+vW058v5o3G1t8bV3prk3FLe3hDN6sMp2FhasHjeKG78ZBff7Enkr9OCANgSlQ7ApEEXbpa2srQg1MeViNNNl335ZGssVhaChyYH1js+Y0gfhvR14a310Qzu68LgZh7wx5Lz8XGzx9Wh6U6EbbYjkFLeIqXsK6W0llL6SCkXSyk/bkgJSCnnmRIxpNFoNKYipeThb/Yz9t8bWbz9FDtiM7l+hA9+7o70cLBBCMFr14UwPsiDE6kFzBvnz8UBPZk0yJMv95zGUGkEYHNUBn7uDgR4mJ4E2hQXB/TkcFIexeWGesdPZxWRnl9KYZmBXw4kM2uYF71c6mcpW1gI/nFtMGUGZebKLCxr8l7HUvKb3Q2ALjHRavzrX/9i6NChhIaGMnz4cPbs2QOoP8YFCxYwZMgQQkJC2LVrV73r/P39CQkJYdiwYcyYMYPU1HOLUn344YcEBgYihCAzM7PmuJSSRx55hMDAQEJDQ4mMjKw5t3TpUoKCgggKCmo270Cj6YqsPJjMmiOp2FhZ8I9Vx/DuYc+tF/vWG2NnbcmSeaP46NYR/HWq2gHccrEvGQVl7IzNoqS8kp2xmUwe1KvVsvlHBfTEYJTsT6h1+v55KpuJb2zh4lc3cs2H2ykur+Tms2StZqSfG98uGEtxRSWf7TjV6H3SC0o5lVnUrH8AumCJCXOwa9cuVq1aRWRkJLa2tmRmZlJersLDtm/fzsmTJzl69CglJSUUFJxrG6wuLvfss8/y6quv8v7779c7P27cOK666iomTZpU7/iaNWs4efIkJ0+eZM+ePdx///3s2bOH7OxsXn75ZSIiIpRtceRIZs2ahZtb084pjaarIKXkg00xDO7rwtd3j+a7iESuH+GDp7PtOWOtLS24sk545cSBnjjbWvHrwWROpOZTWmGsd/5CCfdzw0LAH9EZjAv0AGBrdAaWFoIbw/vxzZ8JBPVyYoRvj0bnCOzlxBXBfVm28zQLJgzA1b7W9PPz/jP8389HcLazwtbKgitCms54Br0jaBVSUlLw8PDA1lb9kXl4eODl5QWAjY0NaWlpVFRU4ODgQO/ejdcYmTBhQk356rqEhYXh7+9/zvFffvmFuXPnIoRgzJgx5ObmkpKSwrp165g+fTo9e/bEzc2N6dOns3bt2tb5sBpNJ2BHTBYx6YXcfWkAbo423DdxQINKoCHsrC2ZPqQ3vxxI5u0N0UwY6Mko/56tJpuznTWXh/Tly92nyShQpp09p7II9nLhX9cG87dpA3n5mqHN7kDunzSAgjIDv1Y5jo8m5xGTXsjXfyZgbWWBvY0l79w4nKDezs3K1PV2BGuehtTDrTtnnxC4/N+Nnp4xYwavvPIKAwcOZNq0adx000019X169+5Nfn4+8+bN46uvvmryl7tq1SpCQkxLTAE4c+YM/frVhrP5+Phw5syZRo9rNN2BU5lFvPt7NB5ONlw17PxW8g9OCcTKUpBfYuCJywa2soTw+PSBrD2SykebY3j68os4mJjHvHH+WFiIGid1cwz1csHHzb4qYcyLWxbtxiihqNzAX6cG8eg00+XWO4JWwMnJiX379rFo0SI8PT256aab+PzzzwGYM2cOGzduxMHBgb/97W8APPDAA/z2228110+ePJnhw4eTn5/PM888Y/J9G6tKqquVarorqXmlXP3BdiITcnhkahC2VpbnNc8ATydenzOMj+8YSWCv5lfULaW/pxM3hvvw1Z7TvL/xJOWVRsb0b9muQwjBpEGe7IzN5LU1xykoM+Bka4WUcGVIyxRg19sRNLFyb0ssLS2ZNGkSkyZNIiQkhKVLl3LFFVeQmZnJoEGD+OSTT7j++utrbPdvvPFGzbXVPoKW4uPjQ2JiYs37pKQkvLy88PHxYcuWLfWOn+1f0Gi6Iv/87RjllUY2PDaxJou3o/LI1CB+jDzDwi2xjA/yYEJQy8NTJw3sxZe7E/h2byK3j/HlwcmBHE7KM8kcVBe9I2gFoqKiOHnyZM37AwcO4Ofnh6enJ1JKNm/ejKWlJYsWLeK9995jxIgRLepH0BizZs1i2bJlSCnZvXs3rq6u9O3bl8suu4z169eTk5NDTk4O69ev57LLLrvg+2k0HZmdMZmsOpTCA5MGdHglANDX1Z4nZgxifJAHC28bYVKm8NmMH+jBvRP78/ldo/jHNcH0dbVnRjPlsBtEStmpXiNHjpRnc+zYsXOOtScRERFy7NixcvDgwTIkJETOnj1bZmRkSCml3Lt3rxw7dqwcNmyYHDt2rPzqq6/ksGHD5IoVK6SUUvr5+dWMbYz33ntPent7S0tLS9m3b185f/58KaWURqNRPvDAA7J///4yODhY7t27t+aaxYsXywEDBsgBAwbIJUuWNDivuX9uGk1rUVZRKae+tUVe+p+NsqTcYG5xOiRAhGzkudqmZajbAl2GuvXQPzdNV+G7vQn8/YfDfDo3nOlDOlb3r46CEKLjlaHWaDSa1kBKyeLtp7iojzPTBjffHEZzLl3PWazRaLoFxeUGvv0zkSNn8ohOK+SNOaE6Ou486TKKQEqp/whaQGczCWo0ddlwTLVvzCwsw8pC0NvFllnDvcwtVqelSygCOzs7srKycHd318rABKSUZGVlYWdn1/xgjaaDUVJeyYNfR9Lfw5GPbx9BUG9nDJXG884Z0HQRReDj40NSUhK6jaXp2NnZ4ePjY24xNJoWcyQ5j3KDkSdmDCK8FUs/dGe6hCKwtrYmICDA3GJoNJo2oLDMgKHSSA8H1SLyYKKq2hnar/mqmhrTaLOoISHEEiFEuhDiSCPnbxNCHKp67RRCDGsrWTQaTefloa8jmfLWHyRmFwNwIDEX7x729HLWps3Woi3DRz8HZjZx/hQwUUoZCvwDWNSGsmg0mk5IYnYxW6IyyC4q57HlBwClCIb3a7xEs6bltGWHsq1AdhPnd0opq/u17Qa0wVqj0dTjh8gkhIB5l/izNz6HiPhsknJKGKbNQq1KR0komw+saeykEGKBECJCCBGhHcIaTden0ij55I9Yvtx9mnEDPLhv4gCEgBd+OQrAMB+9I2hNzK4IhBCTUYrg742NkVIuklKGSynDPT0vvIG0RqPpuFQaJU+uOMhra07g7mjLYzMG0sfVjrH93TmWkk9gLydG+ulue62JWaOGhBChwP+Ay6WUWeaURaPRdAw+23GKH/ef4fHpA3l4am2TlutH+LAzNovnrhx8XpU6NY1jNkUghPAFfgTukFJGm0sOjUbTccgrqeCDTTFMGOhZTwkAXDfCm2H9XNukUUx3p80UgRDiG2AS4CGESAJeBKwBpJQfAy8A7sDCqmxgQ2OV8TQaTfdg+d5E8koq+PvMQeecE0JoJdBGtJkikFLe0sz5u4G72+r+Go2m87HqUDKhPq4M9dJRQe2JNrRpNJoOQUJWMQeT8lrcb1dz4WhFoNF0ItLyS0nPLzW3GG3Cb4dTALhCK4J2RysCjaaTYKg0MufjnYx/fTM/7U+qdy4mvZDotAIqKo1mkq5lJOeWsO5oas37ojIDX+yKZ4RvD/r1dDCfYN2ULlF0TqPpDqw5kkpidgkeTja8vSGaa4d7I4Rg3dFU7v1iHwCezra8cNUQrh7WcWvzH03O48r3twOw6fGJ9Pd04t3fo0nOK+WDW8PMLF33RO8INJpOgJSSRVvj6O/hyF+nDSQxu4TYjEIyC8t4YsVBQn1ceeemYXg42fJ/vxyh3NBxdwZLtsfXfL87LpvE7GI+3xnPjeE+jPTTZaXNgVYEGk0n4Js/Ezl8Jo+7x/dn6kWqL+/G4+l88kcsRWUG3r5xOLPDfHjqskHkFlewOSrdzBI3TEWlkd+Pp3FdmDe9nG3ZHZfFm+ujsLQQPDb93JBRTftgkmlICOEN+NUdX1VUTqPRtDFp+aW8uPIIEwd6ctOoflhaCAb3deGn/WeIzyri2jBvAns5ATA+yAMPJxte/OUoJ1IKeGhKIJYWHadr3+64LPJKKpgZ3IcKo2TziXQKygzcP2kAfVx1WWlz0eyOQAjxH2AH8DzwZNXriTaWS6PRVLHvdA4VlZLHpg+seahPG9yLE6kFSAmPTKnNwLWytODqYV6k5pfyzu/RHErKNZfYDbL6cCr21pZMGOjJ6ICeFJQZsLWy4J7x/c0tWrfGlB3BtcAgKWVZWwuj0WjO5WBSLtaWgov61mbVLpjQnyF9XQj2dj0nyubByYEIBEt2nGJ/Qi5hvqYXaMsoKMPT2bbVZK/LjphMlkckMjvMGztrS8b0dwfghnAfejratMk9NaZhio8gjqrSEBqNpv05nJTH4L4u9ZqzO9tZc3lI3wZDLT2cbHnh6iF4udqxP9H0HcG2kxlc/Orvje4iKiqN/HLgTItDVIvLDdz/5T5u+98e/N0dePHqIQAE9nJi4W0jeHLGRS2aT9P6mLIjKAYOCCE2AjW7AinlI20mlUajAcBolBxOymPW8JaHg4b5urE/Iaf5gVWsOZKKlCqxK7SBev9f7DrNK6uOYaiUXD/S9D5SL608ytqjqTw8JZD5lwbgbFe7rtTJYx0DU3YEK1GtJHcC++q8NBpNG/Da6uN8sfs0APFZRRSUGc6rEUuYbw+SckpIL2g+E1lKyR9RqunThmNp55wvLjewcEsMQL1EsOaISS9keUQS900cwOMzBtU0oNd0LJrcEQghLIHpUsrb20kejaZLY6g0sj0mkxF+bry9Ppox/XsyM7h2Vbw7LotPtsZhIWBQb2eSc0sACPFpeRG2MF+lPA4k5DJjaJ8mx55ML+RMbglD+rpwLCWf4yn5DO7rUnM+Ij6HzMJyBvZ2YuvJDErKK7G3sWxiRsW+06pb7ZwW7CA07U+TOwIpZSXgKYTQalyjuUCklDzz42HmfbaXUf/8nc93xvP8z0coKa+sOf/Guih6u9ji4+bAU98fZOXBZHo4WBNUFR7aEoZ6uWJtKer5CQrLDHy15zRGo6w3dvMJlXfwxg2hONtZ8Z+1J+qdT8lTCmn+pQGUVhj5I9q0lrH7E3Jxtbemv4dji+XXtB+m+AjigR1CiJVAUfVBKeXbbSWURtMVWXUohRX7krh1tC9nckrwdLbl+31JPPfzYfq42FFcXsm+0zn8a3Ywfj0duX3xHuKzinlw8oDz6shlZ23JkL4u9fwEqw4m89xPRxjU25lw/9os3j2nshng6chQL1f+OjWIf/52nC1R6UwapJLXUvJKEQJmDfPmX78d54/odGYGN73LAKqilnpQ1XNE00ExRREkV70sAJO7QgghlgBXAelSyuAGzgvgPeAKlEN6npQy0tT5NV0PKdUqtSs+NKSULNwSS39PR/55TTAWVfkAidnF/Bh5pmbcDSN9uGWULxYWguvCvFl1OIW5Y/3P+75hvm4sj0jEUGnEytKCmPRCAA4l5dUoAiklh5LymDhQ9QOfO9afz3bE8+m2uFpFkFuKh5Mt9jaWDPVy5URqAYeSctkTl809ExrOAcgvrSA6vUA7hDsBzSoCKeXL5zn358CHwLJGzl8OBFW9RgP/rfqq6aY8/cNhjqbk8dX8Mbg6dK2I5T+iMzieks/rc0JrlADAmzcMIzIhh6mDe5OcW0Kgp1PN+deuD+Gv04Lo7XL+Gbdhvj34fGc80WmFDPFyIS5TbeqPnMmrGZOaX0pmYRmhVX4IGysLbh7Vj7c2RJOQVYyvuwMp+aX0rcr8DertxE+RZ3h/40l+P57O7WP8GvQXHErMQ8paX4Wm42JKZvFmIcSms1/NXVdVgiK7iSHXAMukYjfQQwihlw7dlIj4bL6LSOTImXye+uGgucVpdRZvP0VvF1uuHe5d73i/ng5cM9wbJ1srBvZ2rqckbK0s8XO/MNt6WD+VTBZZZR6KzajaEdRRBIeS1Pd1HdJzwn2wEPBjVbnrlNySWkXQy4mCMgNbqqKMTmXWWIzrEXE6GyFgWD+tCDo6ppiG6paTsAOuBwytcG9vILHO+6SqYylnDxRCLAAWAPj6+rbCrTUdCaNR8o/fjtPbxZbRAe7siMk0t0gt5lhyPvsTcxjm04OL+jjXs+nHZxax7WQmf5s2EBur9q3z2K+nPR5ONkQm5HBDuA+J2cXYW1sSm1FIQWmFCu/cm4ilhWBInSihvq72+Hs4EpVaAEBqXinjAj0AavoGG6oczqcyixji5cLZ7IrNYqiXC672XWt31xUxxTR0ds7ADiHEH61w74YMwbKBY0gpFwGLAMLDwxsco+m8/HoomYOJubx5wzCSc0tYeTCZMkNlvUzajkpOUTmRCTk89PV+SipU9I+7ow3f3Tu2phDcN3sTsLQQ3Hxxv3aXTwhRlViWS0JWMUYJN4b7sHTXaZ7+4TBrjqRglDBjSG/srOv/vH17OpCQXUxBaQUFZYaaonBBvdXnsrWyoMxgJK5ql1GXkvJK9ifkMm+cf5t/Rs2F06wiEELULRBuAYwEmg8XaJ4koO5/hg/KKa3pRlQaJa+vjSLY24Xrwrz5fp8yRaTnl3X4TlX7Tudw66e7KTMY6e/pyPs3hxGbUcjzPx/h1dXHWTJvFEajZOWBZCYN9LwgW/+FMMLXjQ3H0tgbr8xDc0b2o8xg5Nu9ibjYWfH7YxPp1YBsfj0d2BefQ2qeSkirNg25O9rg4WTD8H49OJacX+N3qMu+0zmUVxoZO8C9DT+ZprUwxTS0D7VSFyiT0ClgfivceyXwkBDiW5STOE9KeY5ZSNO1OZSUy5ncEp6aOQgLC1Gz6kzNL+3QiuBMbgkLlkXQx9WOv00byISBnvR0tCHY25W0/FJeXX2CnTGZWFoIUvJKefpy89XTGemn/ATf71OW2ABPR567cjBxmUXcPKpfg0oAlP+ioMzAsZR8QJmLQO0ylswbhYeTLU9+f7BBRbDtZAZWFoJR/rrRTGfAFEUwWEpZL0ddCNFseUIhxDfAJMBDCJEEvEhV8Top5cfAalToaAwqfPSuFkmu6RJsO5mJEDA+SIUuVq86U/I6doP2F385QklFJcvvG8sAz/rJXnPH+rNkezxvbYjmoj7O2FtbMm1wbzNJCqE+rlgIiEzIpY+LHU626t9++b1jm7zOt0oR/3lKxXz0rdMvoLoWUX8PJ34+cAYpZb2w3w3H0hg7wL3mXpqOjSm/pZ3AiLOO7WrgWD2klLc0c14CD5pwf00XZtvJDIK9XGvKEPeu3hFUZbJ2RLZGZ/D78XSeveKic5QAqESuh6cG8txPRziYmMvlIX1xNOMD0c7aEncnWzIKyrhxlOl+Cl93pQj2VCmChkxbvj0dKCg1kF9iqAn5jUkvJC6zSPsHOhGN/nUKIfqgonjshRBh1Dp3XYCOu2fXdBoKSiuITMjl3joJSc62VjjaWJKa13HbXyzbdRoPJ1vmXRLQ6JgbRvZj5YFkPJxtef7Kwe0oXcMsvG0EybklzGpBU/vqHUFMeiFernYNRjz5uClzUVJuMa4OKvy0umidOXdBmpbR1DLlMmAeyolbt5xEPvBsG8qk6SbsiMmk0ihrzEKg7M+9Xe1Ize+YO4K0/FI2R6Vzz/j+TYaC2lhZ8F0zppf25Hxs9Q42VvRwsCa3uIJLgzwaHONdrQhyShjqpRTB2qOphHi74tXD/vwF1rQrjSoCKeVSYKkQ4nop5Q/tKJOmm/D9vjN4ONkS7l+/g1ZfV7uaSJWWkF9awf+2xnHbGL9Wj9DJK6ngyJk8lkckUmmU3NwCE0tnJre4AoDJVaUmzsbHTe0aknKU4j6anMfBxNwOsQvSmI4phssdQojFgJeU8nIhxBBgrJRycRvLpunCVK+sF0zoj/VZBdV6u9ixOzarxXO+vvYEX+5OYHNUBivuG3tOXPz5IKVk0dY43t4QTZlBdeaaM9IH/25STdPV3pq8ksZ3BG4O1jjYWHKmShF8ufs0dtYW3DCyeyjKroIpiuCzqtdzVe+jge8ArQg0580Xu05jlJKbws99YPR1tSOtoIxKo6xp1t4URqNkzZFUvtqTQJhvD/Yn5LL5RDqXt0Kxsy3RGby25gTTh/Rm7lg/gr1ccetG/XWX3zuWmPTCel3F6iKEwLuHPUk5xQBsOpHO9CF9ulytqK6OKfnuHlLK5YARQEppACrbVCpNlya/tIKlu+KZObRPgyvrPq72VBolWYWmOYz/vfYED34diV9PBz6dG44QEJ12brZrSzFUGnn1twn8zHAAACAASURBVOP4uzvw0a0jGB/k2a2UAMCgPs5cGdq0QvVxs+dMbgm5xeWk5ZcR3EC5CU3HxhRFUCSEcKeq/IMQYgyQ1/QlGk3jLNsZT0GpgQcnBzZ43qsqhDQpt3GHcUaBUhLpBaV8vjOeq4d5sfbRCXg42eLdw56YBsoenE25wUhMekGj53+ITOJkeiFPX35Ru9cI6kz4uDmQlFNSo3wH9jG5Wr2mg2CKaegxVBbwACHEDsATmNOmUmm6LEfO5PH+phimD+lNsHfD7RerdwnxmUWM8K3vSC4qM/DEioOsOZLKg5MHUFJuxFBp5LHpA2t8AkG9nGrq7jfFiyuPsiIikZ3PTKGXc33nsqHSyEebYwn1ceWyZto8dne83ezJK6lg32lVwmJgb60IOhvN9Sy2QFUcnQgMQuUSREkpK9pBNk0X5I11UbjaW/Pv60IaHdPPzQFLC1GvvHFafinJuSWs2JfE2qOpXDLAnY82xwIwO8ybgDompsBeTuyIzWrSx3AwMZdv9yYgJWw6ns7NF9evarvhWBoJ2cU8f+XILtkopzWpziXYHJWOk61VzY5O03loUhFIKY1CiLeklGOBo+0kk6YLcyI1n/FBHrg7NV6lxMbKAh83+3o1bJ78/hBbq/rk3jM+gGevGMzXfyZga2XJVWfZsAN7OVFuMJKUU9xoPf9XVx/H3dEWG0vBhmNp5yiCA0m52FhaMOWihsMmNbX0c6stRTFCt6XslJhi+FwvhLhe6N+u5gLJL60gLb+MoF7Nmw4CPBw5laEUgZSSA1WNVYb5uPL4jEEIIbhttB9zRvqcEyZaXf65MfNQbEYhe05lc/f4AGYM7cP2mExKK+rHP8SmF+Lv4XBevYK7G0O9XGrqEGmzUOfElL/yx4AVQJkQIl8IUSCEyG9juTRdkOoHc/WDuikCPBw5lVmElJKE7GLySw3889pgfn5wXLP5AYGezvXudzYrIpKwtBBcN8Kb0QE9KTMYaxqwVBObUWSSnBqwsrSo8aP07GZRVV2FZhWBlNJZSmkhpbSRUrpUvdfxYRqScorJNDHEEyAmzXRF0N/DkZKKSlLzS4moqqMf7u9mktnB1cEaT2dbjqc0vF5ZfTiFCUEe9HK2qymLcDS5dmyZoZLTWUUNFpTTNMzDUwIZH+TBDQ3khWg6PnrfqzkvKo2Smz7ZzRMrTO8vHJNRiI2VBf3cmq9BM7yq1+7qw6lEnM7G2c6KgSaYlKq5NNCDzVEZVFQa6x1PySshIbuYS6vqG/XraY+znRVHk2sjouMzVScvvSMwHXcnW76YP7qe017TedCKQHNe7IzN5ExuCbtis86xr++Oy+K/W2LPaWF4Mq2A/h6OJtndQ3xcGR3Qk0+3xrE7LpuRfm71Grs3x+XBfcgrqWBXbBaLt59iwuubKSmvrKmtPzpAFWETQvXqrbsjqG7wrncEmu5CmyoCIcRMIUSUECJGCPF0A+d9hRCbhRD7hRCHhBBXtKU8mtZjRYRqKVlmMNaYbkA5dh9ffpD/rD3B7f/bQ7mhdkUek1HYolX2/ZMGkJpfyqnMIq5sYbmICQM9cbSx5K31Uby7IZqE7GJWHjzDnlPZONtaMbhOo/ahXq6cSM2nsqoZe7Vvob+nXt1qugcmKQIhxKVCiLuqvvcUQjReiL32GkvgI+ByYAhwS1XBuro8DyyXUoYBNwMLWyK8xjzklVSw7mgq14V5Y20p2BaTUXMuKq2AM7klXB7ch+S8UlbsSyQ2o5C4jEKSckpapAgmDvTknvEBvHvT8Bbbnu2sLXnt+lBiM4ooLDfg5WrH0p2n2ROXRbi/W738gqFeLpRW1DZhj0kvxLuHPQ42uruWpntgSvP6F4FwVELZZ6h2k18C45q59GIgRkoZVzXPt8A1wLE6YySq0Q2AK7p5fadg1aFkygxG7rzEn4TsYj75I46Nx9NZ+dA4Nh5PB+DlWUNJyC7muZ+OAGAhoIe9NdeP8DH5PkIInrvy7LWD6cwa5sUofzcSs0uISS/k2Z8OA6p5e12qM5w3nUinj6sdf0RnMC5QN13XdB9MWfLMBsKASAApZbIQwhSvnTeQWOd9EqpJfV1eQuUpPAw4AtMamkgIsQBYAODr69vQEE07siIiiYG9nQj1cWVmcB8iTucQk17I78fT2Xg8jVAfV3q52PHcFYN5a0M0lwf34VhyPneM9Wv3hvR9Xe3p62pPqI8rR5LzKC4zcM3w+l26BvZ2YuJAT95YF8Xao6nklVRw/8SG6yBpNF0RUxRBuZRSCiGqi86ZajhtyLMnz3p/C/C5lPItIcRY4AshRLCUsl6oh5RyEbAIIDw8/Ow5NO3I3vhsDiTm8twVgxFC8JdxAcwO8+bK97ezePspDiXl8siUIAAuCfTgksCG69i3N3bWlrw6u+GyFkIIPrg1jGd+OExkQg43hvsQ4tNwHSSNpitiiiJYLoT4BOghhLgH+AvwqQnXJQF19+A+nGv6mQ/MBJBS7hJC2AEeQLoJ82vamZLySh5ffpB+Pe25+WL1q7WwELg72TJruBeLtsZha2XRogbpHQUXO2s+um2EucXQaMyCKQllbwLfAz+g/AQvSCk/MGHuvUCQECJACGGDcgavPGtMAjAVQAgxGFXgLgNNh2Td0VQSsot5dXbIOY1Kqs0td48PwFv3qtVoOhWmOIv/BqyQUm5oycRSSoMQ4iFgHWAJLJFSHhVCvAJESClXAo8Dn1bdQwLzpJTa9NMBkVKy8mAyXq52jBtwrrlnqJcrvz50KYP76lozGk1nwxTTkAuwTgiRDXwLfC+lTDNlcinlamD1WcdeqPP9MZqPPtKYkZyicp7+8RBHzuSTll/K/EsDGk3s0nZ1jaZzYopp6GUp5VDgQcAL+EMI8XubS6YxO9lF5cxeuINNJ9LJLirHYJRcPcyr+Qs1Gk2noiUZM+lAKpAF6CLt3YDHlh8gOa+Ur+8ZQ2lFJTtjsxiq+9FqNF0OU3wE9wM3oVpUfg/cU2XS0XRhErOL2RKVwRMzBjLKX9XlGV9VqE2j0XQtTNkR+AGPSikPtLUwmo7DqkMpAFwz3NvMkmg0mramUUUghHCRUuYDr1e971n3vJQyu41l07QzSTnFbD6RTn6pgSXbTzG8X492zwTWaDTtT1M7gq+Bq4B9qNDOuqEiEujfhnJp2pn1R1N55Nv9lFbUJnXfP2mAGSXSaDTtRaOKQEp5VdXXZiuNajo/H22OwauHPYvvHIWFgF8OJHNTJ8wQ1mg0LafZ8FEhxEZTjmk6L4nZxRxMyuOGkf0I8HDEz92RR6YGnZM9rNFouiZN+QjsAAfAQwjhRq1pyAWVT6DpIqw5ohzDLW3+otFougZN+QjuBR5FPfT3UasI8lENZzSdHCklN3y8i0NJeQR7u+Drrh3DGk13pCkfwXvAe0KIh00sMqfpZJxMLyTitGoz+dBkXX9fo+muNJtHIKX8QAgRjGo3aVfn+LK2FEzT9vwRpQq97nh6iq4YqtF0Y0xtVTkJpQhWo3oQbwe0IujkbIlOJ6iXk1YCGk03x5Tm9XNQPQNSpZR3AcMA2zaVStPmFJUZ2Hsqh0mDdNkIjaa7Y4oiKKlqHWkQQrigis/pZLJOzqKtcZRXGpk8SNcP1Gi6O6YoggghRA9Ue8p9qCb2f5oyuRBiphAiSggRI4R4upExNwohjgkhjgohvjZZcs15s+90Nu9tPMn1I3wYO8Dd3OJoNJq2JnZzk6dNcRY/UPXtx0KItYCLlPJQc9cJISxRYabTUf2L9wohVtatXCqECAKeAcZJKXOEEHp52g6siEjC0caSf1w7FCEabjKj0Wi6CBWlsPLhJoc0lVDWaCdvIcQIKWVkM7e/GIiRUsZVXfMtcA1Qt4T1PcBHUsocACmlblrfRhw5k4e/hyPWloLVh1O4bGgfHGxa0o5Co9GcN6V5cHIDDLkWLNvx/64oC3a8A3mJTQ5rSqK3mjgngSnNiOAN1L17EjD6rDEDAYQQO1B9jV+SUq49eyIhxAJgAYCvr28zt9XURUrJI98e4NeDyYwP8uC20b7klxq4JkyXl9Zo2oWSHPhiNiTvVw/kS/92fvPseB8SdoHRAC5eMP0VsGukPayhDH66D479ArISLroKVUe0YZpKKJt8ftLW0JDN4ezG9FZAECo81QfYJoQIllLmniXLImARQHh4uG5u3wJOZRbx68FkhvfrwbaTmRxIyKWXsy3jtG9Ao2ldDOVw8GvIiYfCDOg3CgbPgu/ugLSj4BUGW/4NgdOgT0jL5s6Igg0vQI9+YN8TYjfB6V1w7X/BZ+S54w+vgKM/wuj7IPRG6BsGt5yHIqhGCDG3oeMmJJQlAXXLV/oAyQ2M2S2lrABOCSGiUIphb3NyaUxje0wmAO/eNJx//nacA4k5LJobjpWlKXECGo3GZLa/DVteAwsrtVI/8CX8+iggYfYn0H8yfDoZvroR7tkELi2o7bX9HbC2h3s2g6MHxP0BP98P/5sKd60Gv0tqx0oJuxZC72CY+W8wwQ9oirFqVJ3v7VA5BZE0n1C2FwgSQgQAZ4CbgVvPGvMzcAvwuRDCA2UqijNBJo2JbDuZSb+e9vh7OPLJHSOpqDRiZ21pbrE0mq6FoRwilqjV/q0r1MM3aS8cXwkO7jDsZjXu1u/gf9OU0rjiDdPmzk2EQ8th9L1KCQD0nwgP7II3BynzT11FELcF0o/CNQtNUgJgWtRQPXezEMIV+MKE6wxCiIeAdSj7/xIp5VEhxCtAhJRyZdW5GUKIY0Al8KSUMsskyTXNUlFpZHdsFlcNU8ViLS0ElhZaCWg0rc6JX6EwDS7+ECyqdtv9LlavuvQJAb9xcGqb6XMf+k7Z+UffW/+4nSv4jVUP/roc/AbsekDIHJNvcT72gWKU+aZZpJSrpZQDpZQDpJT/qjr2QpUSQCoek1IOkVKGSCm/PQ95NI0QeTqHgjID44M8zC2KRtO1iVwGbv5qR9Ac/pdCxnHlR2gOKdVuwPcSNf/ZBEyEjBOQr0rJYyiDqDXKOWxlegEIUxrT/CqEWFn1WgVEAb+YfAeN2Vi0NQ5Xe2utCDSaajJj4OhPEPkFxPzeOnMWZ6sVfvD1tbuBpvAfr76e3tH82FNbITMKQm9o+Hz/Sepr9Wc5uR7K8mHotc3PXQdTfARv1vneAJyWUia16C6adicyIYeNJ9J5auYg3WmsO5ObAM5e7Ru73lE5tBx+eRAqy6sOCHhgN/S6qOVzVVbAxpfBwUOtwmUlXHSladd6DQcbZ4jdWP+BbTTC4mlqxzD9FSgrVIlgbgEQelPDc/UdpnYKR76HAVPgt8fV+ICJLfo4zaovKeUfUso/gP3AcaBYCNGzRXfRtDtvr4/G3dGGO8f6m1sUjTmQUsWdvxsCEYvNLY15kRL+eAN+vAd8LoZ7tykFYOOoonzOhx3vwc4P4PcXYcurStl6NZqDWx9Laxg0E46vUgqlmtPb4cw+NXfsJvjj30qRX7tQydoQQkDIDWrn8PVNUF4Mt3wDVjYt+jimmIYWCCHSgENABKreUESL7qJpV3bFZrE9JpP7Jw3A0VavBLsUFaVgrGx+3Pa3YcP/qe+TuvG/a2UF/PIQbP6nWlXf8SP0DYVeg2HsQ3DsZ4g6J4e1aTJPwh//UVnCj51QsfxzlpgcoQOoa0uyIb6O03j/V2DrAh6DYPk8+PN/Sua6EUENEXIDSCOkHYar3lafrYWY8pR4Ehgqpcxs8eyadmPlwWTKKioZH+TJv1Yfo7eLLbeP8TO3WJrWpLJChR4KAfNW1c8qlVLZvmN+Vw+qpD8h5EYozoL0Y43P2ZUxGuH7v6gQzol/h0nP1H9Yj38MTvymzEX37wTn3qbNu/lVsLRV4Z9OvWD42VHxJhA4DWyc1O9swBQozVdhoMNuVnItngGFhTDxqebn8hwEvmOVEglpxJfQDKYoglhUpJCmg3IwMZdHvtkPqL9zKWHhbSN0vkBXY+//1KoPAV9eDzd/A06eUJCqzB6ntoJjL/VgGPdXmPwcbP4X7N6mlIhlN/EVSQnR62DXh2rFPeOfcEkDRdesbOH6/8GiifDLAyr+vzlnb+ZJ9fC+9FGlBM4XazsV2XPkJ5j8vJLTUKIUQQ9fmL8eck6D+wDT5rtzlfp6nkUkTVEEzwA7hRB7gLLqg1LKR87rjppW59NtcTjbWrFs/sUsj0jC2c6KK0JakLWo6fgUZ8Pm19TqMfwv8MM98P5wsHVWq0kkXPUOjJhX/2HWa6hyjmbFnp9TtLMhJXx9o4qecfFWmbWj72t8fK+LlKJY/YQypc34p6oNZKxUSvZstr8DVnYw5sELl3XiU6oMxPrnQVio0hE+Vfm7PXzVy1QuMBjAlKs/ATYBhwHjBd1N0yqUG4xYCLCytGB/Qg5rjqQy/9IAwnzdCPN1M7d4mrYgcimU5alokj4h8Je16pjRABbWSjn0CT73ump7cfqx7qEITvymlMCkZ2H846Y9IEfdrVb6uz5UYaVleWDtCI8dBfs6/085p+Hgt3DxgoaVREtxHwCXPqacwhbWKoLITAmfpigCg5TysTaXRFOPtUdS+PVQCr2d7Zg/PgAvV7ua3gHzl+4lJa+UiwN68vWeBFztrblrnL95Bda0HunH1QPIuY96X1mhHIcBE2qLlXkNV6/m8BykHjIpByH4uraTuSMgpXLi9uxvuhIAZU65/D/g5qeKu9m5ws734cRqCLutdtzO99WDuiEz0/ky8SlI3ANxmyFweuvN20JM+UltrioD/Sv1TUPZbSZVNyWvpIKHv9mPlYVg04l0+rjYkVlYxpIdp+jjYsc3C8aQXVTGtpPKbx+TXsj8SwN4ZEoQrg7dxP7b1ZASKoprwwNzE+HTqarK5H071MNs878gPwmufrfl81vZKoWRuKd15e4oRK+Dokz1wI7dBKmHVI2dlppKhICxVeaeasf7sV9qFUFFCRz8DoLngGsrlnC3sFQRR3sXw+CrW2/eFmLKT6vaJf5MnWMS3be41VkRkcjW6AwcbCz5y7gAnr3iIpJzS1l3NJWFW2K4Z1kEfVzs6OFgzbs3qdXgJN1zuHPzy0MQvRbu36F2AGufBkOpKhuw93/KnLP9HRhxJwSd54qx32j481OV+NRU2YGSHBXDPvHvqtJlR6aiRFX2PFRVlaYwDdKOqJ1UC2rsNIgQMOQa2POJUjKOHqpsQ3lBbfG41sShJ0x8svXnbQGmFJ0LaA9BujtGo2TZrtOM8ndj+b1ja8xAvu4O3DOhP0O9XJi75E9i0gv569SgrqkAijJVJydTIyU6Oyd+U6WKQWWEjrkfTqyCqS9C/HaVrOTUS5k6Zv77/O/jO0bZv5MPqKqUkcvg9h/VA6gu+79USqffaBh0+fnfr60xlMOXc1SJholPQ9ZJleULMOqeFtXYaZSwO9TPLOIz9fDf/rZyPleXh+hitGU/Ak0L2ByVTkJ2MU9eNqjBPsKXBHqw7e+TSc4tIcS7hxkkbGMyomHZNVCQDANnKufnwMvMLVXbUZytVrR9QlRy0aZ/qO5Tjp5KIYTdAZ+MV5mlN30FNg7nf69+VY0Bl14NlVXW3eTIcwuknfit6tyBjq0INr2isnBnf6Ie0kajUpa7/wsj72yde/S6CAZMhd0LVQaxrFSJY6bUEuqEmPKpRtV5jQdeAma1oUzdjopKI6+vjcLHzZ6ZwX0aHdfX1Z6Rfj2xsepif4x5SfD5lWCsUI645AMqBDDtqLklaztWP6kyS6/9r3Jsht6kkr9G36vMMk6ecNsKuOw102vYNEZ1vHtlGQydrb7PiK4/piiz1o+QcuDC7teWxO9QD+bw+bVmGgsLmPI8PJ3Y8s5fTTH2AfU7cu4D922HIV33sddm/Qg0pvNjZBJRaQV8csdIrLtb57BKAyyfq2y+d/+uVmJhc+GjUZByCHoPNbeErc+xX1SRsMnP1T64Zn0AQTOqestW0Sek9R5sl72q/A5Xvavq12dG1T8fvVaVKeg1VPXW7QgYjSrOfs/HqobOiDtU+KaLj4r3P5vWXq0PmAq3Llex/Web0boYbdqPQAgxUwgRJYSIEUI83cS4OUIIKYQIPw95Oj174rLp5WzLjCEmprh3JQ5+owptXf1ubZx7zwAV8phxwryytQXF2bDqMVU1sm4Tcytb5eS0tmub+459UCkbC0tVyybzZP3z1QlYYbcrx2t1fXtzUZgBSy6DH+ZDSa4qorb2abVbmfL8hZnKTEUIZZ7s4koATPMR/Ept03kLYAiw3ITrLIGPgOmo3sR7hRArpZTHzhrnDDwCdNH4tuY5mpzPUC+XBn0DXZqKElX90XukquVejaU1uAeqmO6uxp6PoTgT5v5svpIPngNVjHw1UkLCblW6uDqzNXF3rRmpvSlIg6VXqVDaa/8LoTer1X5GlEqMG3yNeeTqwpiyI3gTeKvq9RowQUrZ6Oq+DhcDMVLKOCllOfAt0NBv8B/A60CpaSJ3LUorKonJKGSol2vzg7saB76G/DMw7aVza6R4Dup6O4LyIvhzEQy6snVt2S3FY5BSRsVVqUA5p9QuwHcMeIWp4mVntz/MiILX+0NCG6/XyouVqTAvCW7/QRV0qzb5eA5SyqmLOmzNSaM7AiFEINC7qhdB3ePjhRC2UsrYZub2BhLrvE8CRp81VxjQT0q5SgjxRBOyLAAWAPj6tqD+RicgOq2ASqNkiJeLuUVpX6RUoXl9QhsOyfO8SNnSK0o6bkx7eRF8djnknQFLG5UENniWsvV7Djx3fOQXKlb/0kfbX9a6VJedSD2smqBXP9x9x6hELP/x5yqCoz8pZ/a6Z+Dujedd3KxRMmNg6xtw/FeoKILrF4P/uNa9h6ZRmlKt7wIFDRwvqTrXHA39pciak0JYAO8Ajzc3kZRykZQyXEoZ7unZCjU+zICUkuJyA1LKesePJecDMKRvN1MECbtUJc3wuxp+qHgOAmTH3hXsWqhKNwy8DAKnqNXs+ufgv2OVo7sulRUqLt137LkNzdsb76oGKmeq+hTEbgRbV/CsUhD9J0FOPGSfqr0meh1YOyh/zr7PWleeilL4crbKoQi9AeatvvCkME2LaEoR+EspD519UEoZAfibMHcS0K/Oex8guc57ZyAY2CKEiAfGACu7qsN40dY4hrywjqlv/0FybknN8WMp+TjZWuHbsx2cXx2F4mz46V7lnGysfnrNw2qf6fMm/qnqz+cmNj/2QilMhx3vqrIA1y6Eaz6C+7fDw5GAgEPf1R9/9CfIS1Tloc2NvZvywZyJVArr8PeqlEK1yaX/JPW1eldQmK7yDsY9Cv0nw9pnVD2k1uLPRSpf4uav4Or39E7ADDSlCJoKXzBlr74XCBJCBAghbICbgZXVJ6WUeVJKDymlv5TSH9gNzKpSNF2OTSfS8e5hT0Z+GXcs3sOhpFze3hDNb4dSGNzXGQuLbuQo/vURVUP/xi9UGeWG6OGnkqtM6a4lpcqIXXIZHPkBtr7euvKCCmVM3Fv7APzzU2UamvpS/XHuAyBwKhz9WV1TTeQy6DkAgjpIkpz3SPWzXfesUgx1G6B4BIFrP1VWAeDkBvV10EyVxGXjCGv+3jpylObDtjdVwbX+k1pnTk2LaSpqaK8Q4h4p5ad1Dwoh5qPaVTaJlNIghHgIWAdYAkuklEeFEK8AEVLKlU3P0PkwGiWLtsWxIyaTHg427IvPZvG8UQzwdOJAYi63j/FjxpDe/OXzvcz6cAdCQLifG/dO6CYlFQBiNio78NQXwGdk4+OEUBEsSXsbPh+9XmXj5iUqRVCaC0OvU7XiD3yjSg+0RnGwSoNa+e/9HxSkAEJlPZ9YpXwBHoHnXjP0OhWXH79N2eAL01U5hPFPdBxHp/dItWspTFXlK+qWWxZCVSrd9REUZamdgWMv5c8RQu1qNrygdmveTfwOTeHQd6qsyORnmh+raTOaUgSPAj8JIW6j9sEfDtgAJsWVSSlXA6vPOvZCI2MnmTJnR+ajzTG8tSGa/h6OHE3Op7DUwLJdp7kh3Icyg5FwPzdG93dn8xOT+PVQCpcGejCoTyMr4q5I+gn46T5wC1D9YpvDJxyiVitTkqWNcm4WpqmV6qFv1Qo7eI56OLkHqjozeYlqV7DuWbhx6YXJW1YIK+ZBzAZVjmH6P1QLyIglqg/AqPkNXzdklrr/ro+U43XvYpWsNfTaC5OnNal+gNu6qHIWZxNyoypAd/RHpYx9R9f6ckbeBdvegk3/VDWLztdxbKxUCtZrxIUrFM0F0agikFKmAZcIISajbPkAv0kpN7WLZJ2QdcdS6xWNe3z5QX49mExfV2VlG+mvVl29XOyYf2k3q+WXm6hiw4WlytY0pTCYT5VT9bMrVJhpmXKsIyxgwlMw4UmVaFQXNz9l5tj0D9WUfNDMxuevKFFhkQ3V9a80KCUQu1Fl44bfpY6H3qDKYCRFqB1BQ1jbq+YlW16Fj8ep2He/S6HXkOY/c3vRJ0TV3R85D2ydzj3fe6hKetvxnlKu1Z8fwM5FNX5Z+3elKOrmgJhKpQF+WqCCAa5ffN4fQ9M6mFJiYjOwuR1k6dTkl1ZwLDmfh6cE1SSG3RDuww+RSby9IRo/dwd6ObdR1mhHJTdRrejzk1XUSUUp3LOp4dDKhvAKU18zjquHTciNytzj4t10tuclj6juXRFLmlYE6/8PIhbD346Ci1f9c38uUjuBq96p/xAE09oIjrpbPUSLs5VdvXrn0lGwsoWH9tU3CdWl2gT0/V/U++pEs2ouvkftytb/nzKFtfSzrXtG7dymvaQjhDoAF9boUlNDRHw2Rgmj+9c+oEYH9GSEbw8iE3IZF+hhRunamZJc2PhK/VaKLl5w/aemKwFQK9U7fgan3tC7BatpKxtlkjm5XvkPGnpI5acoB640KmU16u7ac5UGVXXSb5zyB5wPju7wwC71oLXroKHBzbVbHHIt9PynatHY96xdk4UlDLsF1jylfCDOJpZHKS9W+QJ/CPOLzQAAGcRJREFULlLmwbplNjRmQyuCVmJPXDY2lhaMqNMzWAjB8nvHkpxbSi+XVqiR3hmoKFGmnIzj6iE69iFw8z//1fCAyed3nVcYHPhKZaj26Hfu+Z0fKCXl1EeVX66rCI6vVOaQy/9zfveuxs3vwq43NxaWKpwz9UjDtX08qpR6ZpRpiiAzBr69BTKjlRKZ9lJrSqu5ALQiaAXySir4IzqDYf1csbOu33zaytICX/dulCOw9hnV/OTWFTCwERt6e1A3D+FsRVBRopqwDJ2tTE27PqrtRAUqNLRnfxjYgWvytxcBE9SrIWoUQXTjY6pJP6FKjSPhjp9gwJRWFVNzYXSQWLbOS6VRcs2H2zmRWsA1w1uxl2lnJH67yjq95BHzKgGA3sHKJJUcee65qDVQlqcqbQ67Re0M9ldVVs9PhoSd6nhHCfXsqLh4gY3TuZVMzyY/RTUdsrCEv6zXSqADoncEF8j+hBzis4r5z/Uh3DSqa9VBahHZp1TIpIsPTH7W3NJUNW0PU41MQCV/bfwHZMeqMFRnL7WKtbBU/oSIJTD2YZVlC8o+rmkaIVTyWVNVYo1GFR1Ulq9qFDWUd6ExO1oRXCAbjqdhbSm4PKSvuUUxD2nHlMMwfpt6P+ezjlMkbsAUlWVcnK162u5bCn1DVUG7sNuVEgDlH1hxJ/yrt9od9A5umVO7O+MxqPZ33xAn18GprSr6qiUOf027ohXBBbLhWBpj+rvjYmem2vLmJCsWllW175v6ogrx7EgO0sCp8Me/Yf3zynE87q8w/ZVzx110lYpekRJcfVQ9HY1p9BqswkiLsxsO6d29UO0Swxpsfa7pIGhFcAHEZhQSl1HEnWP9zS1K+1NWCF/NUeGXd63tmCtorxEqaerAV8qxOfm5hsdZWukIlvOlOiP4TCQETat/LvWI2g1Me1n9jDUdFu0NM5GiMgOrDiXz9oZoft5/BoA1h1U7v+lducVk6hH45SH45lYwlKljUqqiY9mn4MZlHVMJgHr49J+kvp/yf6ZlM2tahtdwQDRcJTZymar9NELvBjo6Wk2bQJmhkpnvbSUxW5WPtrQQBHu78uvBFEb5u+HVo4PYxFuTSgOseVI5Ua3swVCiEoGmPA9b/g0HvlQlHvwvNbekTTP2IVXJdPDV5paka2LrrHpHnK0IjEaVjxE4rVv0/O3saEVgAlujM0nMLuH160OZMNCTqW9t4c4lf3Imt4RXrhlqbvFaH2Ol6hdw5HsY84Cq3bP2GVVorKwQ9vwXht/WuKmlI9HvYvM3gunqeI9UJUTqZnEn7VXVWnX0VadAKwIT+O1QMj0crJk9whtrSwtenDWU5XsTGeLlwjXDOnnuQOJetXI79Ycy/Yycpwq8Hfle2c2rSwBc8aYqnrbnv6oY2VXvdqzaORrz4T1C+WFyE2qDBfZ8rHaSAztI/wVNk2hF0AQVlUa2x2Sy4VgaVw/zwtpSuVRuDO/HjeENlC3oTBgrYctrytxjYa1WzRbWsPZpdX7YrfXrwNg6qWzhTa+oTlVnV/3UdF+qq8QeX6kqs57eqaqSTvz/9s48vqrq2uPflRAGgYCIAZkckDoxFRF5KoiCiNU6FBWcnrYodX7VUqdSfA4tan3WOtWhWpVW0Vq1VhG0CqhFhVi1akVFlEFBoARBIAzJen+sfckxhCE399ybm7O+n8/95Aw797fXuWefdfa09hX1N86S8y1idQQiMgz4LbYwze9V9YZq5y8FzgY2AkuBH6nqvDjzVBsmv7+Yix59G4CjezageQKVFfDXC+DdR6H36XDUDVUrhS350NYT7jly8/9r2c6WZHScKLv0hK6D7aWi92k2XLe4Y/1YltPZLmIbNSQihcCdwFHAvsApIlJ9RsnbQF9V7Qk8AcSwxmD6fLZsNQD/d1IvDmko0UMrK+Dp88wJDLoSjr/z28tFluxjweJqCjLmOFti6PU2e/vhY63jeMClfg/lEXEOH+0HzFHVuaq6HpgIHBdNoKpTVXVN2H0DW+C+3rCwbA0lLZswfP9Om9YYyDtU4ct3oOxz23/uUlse8PCxMOiKnGbNaUC029cm6y1+z+Zu9Dol1zlyakGcTUMdgQWR/YXAgVtJPwp4vqYTIjIaGA3QpUv24vksLFtLpx3zdGioqq2b+/Iv4av37FjLXWwkx8E/saGfjpNJ+p8PqxZbCI/GzXOdG6cWxOkIanqF1hoTipyOrYd8aE3nVfVe4F6Avn371vgdcbCwbC29O7fOllxmee03Fl9np24W56ViI3z+CjRtbQvHO06mEYGh1+U6F04axOkIFgLRoTWdgC+rJxKRIcDPgUNVdV2M+akVFZXKlyvWckw+dhJ/+DdzAt2H2zKJhSEO0oGjc5svx3HqJXE6gllANxHZHfgCGAmcGk0gIt8F7gGGqeqSGPNSaxavLGdjpdJpxxx0eFVWwKJ3Ye40mDvV2l03rgetsFmcXQ6CXiOq1vSNsnoZPHOxxdk57q4qJ+A4jrMFYnMEqrpRRC4EpmDDRx9Q1Q9E5FqgVFWfAX4NtAD+HDpj56vqsXHlqTYsXG592FntIyj/Gib9zNbaXVtmx9p1t9mZqTbXRe/a4i9vT4Cz/26jfCo2wr+fNqcxdzqsWwXH3wVFTbOXd8dx8pZY5xGo6iRgUrVj4yLbQzb7p3rCwjKLK9S5TZZqBBvKLbDbgjeh58kWCnmPQ6FFyeZpV34J9w6CCT+AI66x0A9LZ9tC6Z3721q7JftkJ9+O4+Q9PrN4C6QcQYfWWXqrnnw5zHsNht8PPU7cetriDnD6kzDhBHjyHGjdBU6eYHH1fXlFx3FqiTuCLbCgbA3tipvQpFHhthPXRNk8G6q5tswW7VhbBmvD3/WrrQ1/zyGWZvZz8NaDttbvtpxAivbdYdQU+OApW2Graav08uk4TuJpUI6gslJZ+s06WjUrommRPcDXb6xk9uKVbKioZP9dty8c7tdrNvDCB4vpv8dOtc+EKkz9pU23r44UWPNNYWOb1JWioMje5msbzbPNHjDgp7XPo+M4ToQG5QhumvIRd0//lMaFBdw6sjdH7teeEfe+ztvzVwAwcXT/7Xq43/3Kp6xat5FLjkhjwZUpV9nyfL1OhR7D7cHfrI39bVJc1XSzfC58OhWatYZuR1pQN8dxnBzQYBzBnCXf8PtX53L43iWUf/0VVz31HtM+WsLb81dw2bC9eHTmfM7/0z9p07wxTYsK2G+XVvTv2obOO+5Ak0aFdO9YjIiwfmMlj81awJH7tmefXbYzcuKyT2ykz7wZMPtZOPBcGHbD1sM0t9nDPo7jODmmwTiCO6fOoWlRIb/ps4Tip87gATmB60qHM6BbW87bfQnHrHiNyxYPpnXzpqxev5Hn31/EY6VVETB6dWrFPWf05Z0FZSxfvZ4RB4S5cJUVFqt/9iSLt77mPyEe/y12fsUCuHuAreBV3BH6/RiO/JXH6nccJ29oEI5g6ap1PPevRYzuswOtJv8QCooYVfEXTj22P43Xf4A8OJ4uWsnEE/ayiVjYzOHZi1ey7Jv1LCxbw/hJszn1vjfYqUVjSlo2YUC3trD0I3h0pDXjFDWHtnvaEnyl91tQrc4HwIvjAIULZtpkL8dxnDyjQTiCiTPns76ikh+XP2Ajcs55CV66jmYvhMBqPU6yh/rL18M+x0Dj5hQWCPt1qBpps1e7lpz1h1nMXbaa8wd1pVFhgS3Qvma5Denc+2goamZLNd7aHabfaNE7UwtwuBNwHCdPyXtHMHHmfG6fOodzuiym5SdPw4Ax0L4HjPyTLZDRtLU9sOfNgAePhucvs5m6H0+GDWuhfU/o3I++u/VhxpWH8/HiVXTv2ArmvGQzdY8c/+0hnU1a2MpdL46D5Z/asM2DLsrdBXAcx6kjopq1YJ4ZoW/fvlpaWgrA2vUV9LnuBa5uNYkRax5BWrSH81/f8vJ4f78GXgtt+0XNLWzD6hDiqMfJMORqaBWWRLh/qI3xv7AUGjX59vdUbID7DrMYQAddZItyOI7j1GNE5C1V7VvTubyuEUz7aAk/1/sY+c1L0P1E+N6vt75G6mFXWfNOyT7Qbag94FcthtIH4NVbLF5PnzOh62EW6uHIX23uBMACuR3/O3jpOovB7jiOk8fkdY3g+gefZuznZ1J5wGgKvndT3UbqrFgAr94Mb/8RKjdCo2bw0w9t/L/jOE6e0yBrBOUbKtjrs4fZUNCYokGX1324ZuvO8P3fwsDLbD3f4o7uBBzHSQR56wimzHyfY3mF/3QdTvvmGVxYvlVHGDgmc9/nOI5Tz8nLUJWqyrxXH6GJbKDdkAtznR3HcZy8JlZHICLDROQjEZkjIlfUcL6JiDwWzr8pIrttz/dO+3gp/dZMZ2WLPZB23TOdbcdxnEQRmyMQkULgTuAoYF/gFBHZt1qyUUCZqu4J/Aa4cXu++53JD9GvYDbNv3uSh3JwHMepI3HWCPoBc1R1rqquByYCx1VLcxzwUNh+AhgssvUne/my+VxSdj0rmu9O4QE/zHimHcdxkkacjqAjsCCyvzAcqzGNqm4EvgY2ixMtIqNFpFRESheua8r/ci6NLpgBxbvElHXHcZzkEKcjqOnNvvqkhe1Jg6req6p9VbVvp5KduHDMtRTvkMVF5R3HcRowcTqChUDnyH4n4MstpRGRRkArYPnWvrRpUSFtW9Qw29dxHMdJizgdwSygm4jsLiKNgZHAM9XSPAOcGbZPBF7WfJvq7DiOk+fENqFMVTeKyIXAFKAQeEBVPxCRa4FSVX0GuB+YICJzsJrAyLjy4ziO49RMrDOLVXUSMKnasXGR7XLgpDjz4DiO42ydvJxZ7DiO42QOdwSO4zgJxx2B4zhOwsm79QhEZCkwL8Nf2xZYluHvrM+6SdVOos1J1U6a7vZo76qqO9d0Iu8cQRyISOmWFmxoiLpJ1U6izUnVTppuXbW9achxHCfhuCNwHMdJOO4IjHsTpptU7STanFTtpOnWSdv7CBzHcRKO1wgcx3ESjjsCx3GchJMIR7CtVc8aKkm0O4k2QzLtdpszRyIcAdAaNq15kFVE5FQR6RW2s33jNo3kI+uFRkRycX+1CNqF2RYWkWNFpGu2dQOb7E3QAzIn5bohlukG7QhEpJWIvABMhk3LYWZLe4iIvArcCnw36GelZ15EhorIDOAOETkty9rHisil2dCKaIqIlIjINOD3AKpakUX9ISLyOhZWPavrp4rI0SLyInCLiAyE7PzWInK8iFwXt84WtHNSrhtymW7QjgAoB8qA7iJyEsT7phgeSM1E5HFgLHA98ASwQ9zakTzsDFwL3AQ8AowQkSvDudh+bxFpJCKXA7cBN4tIb1WtzIbNoUCUh09PETkq5ClOe0VEWojI37DfeizwBrBr3NqRPOwG/BK4HfgQGC0iZ8elH2wuDBo3A1eIyIBM62wHWSvXiSnTqtogP1hVuR1wCXAMsDhyTmLWPi6yfTrwepZsFqA7cE/k2L7Af4C2cdsOHI9VXX8CvJnF37og2HkDcFy2rnfQHhHZvgB4PIvag4E7wnZTYBDwLrBjnL910GkJnANMy5a9QTsn5bqhl+kGUyMQkYtF5D4R+ZGIiFrTwErgaFV9FviXiIwTke6qqplsX4tonwOgqn8NxwuBz4APRKTz1r6jDtpnisgRQVeBb4CDRKRNOPZv4M/YW2OmtS8WkRtE5ORw6DlVLVfVW4ESETk1pCuKSXc4gKpWYuthfwf4B7BIRM4VkW6Z1K2mfVLQfiwcLwRWAAtEJJZFtUXkRBE5MHJoITBcRJqE6z4NmAH8IsO6qfv77HBouqquUtX7gOYiMiqki6MWkpNynbgynU1vHqPXPAurlg8DpgNXAV2BEuD6kOZHwEZsmUyAohi194ic74Gt39wywzbviFVRFwH/Agoj5x4GJlRL+yawe4a0BXsj+we21vSH4TqURNKcAHyRYZu3pNsG6AtcHdKNAVYDfwv7jWLU3jmS5iBgdgz3d0m4t74EngYKqv3Wt0by2CvcF+1iur+vBLpGzh8FfECohWTY7pyU6ySW6YZSIxgM3Kiqk4GfYtXkk4C1wFGhY+li4GWqQlhnqoOpunZjrOoIgKq+F/KR0fWYVbUMeAHYB3gLGBc5fSEwTEQOCPursSaD9RnSVuAwYKyqPoE9IHsBR0bSPAV8LCJjwDraYtLtDRwBLAYGisgk4IfYA3tu+Nc6dxxvxeZhkTQzgIUicmxd9appLwH+GrQWAT+OnL4WOEZE9gt5LAdWYW+RmaCmsnVaJG/PU9U/0TJVU4pROxvlOnFlOq8dQaQq+jbWXoiqlmLV492BQ4AXgZmq2ltVhwKDRGT3UGji0H4D6CAiB4d0gv24TTNYbU19z8OqugK4C/iBiOwa8rESuAb4hYiciXVydSeNh0P1PEfsLgUGBL3JwMfAfiKyVyT5ecBNIrIY6BiT7kdAT8whLABmqep+WCEdJCIda/tbp2Hz3iFdMTAb2FAbve3Uvh34N3YvHS0iu4R8zMFGK90lIodgD6wSoDLdPFTTrX5/v07k/g5cDowHPgHa10V3G9qxluuklOmayCtHICIHS2SctlrbMNjbX4GE4XNYVfULrENrnKqOjXxNF1X9LEbt97G3tg4hnWIFc3W6N2kN2hr+loe/s4DnsREkqTR3YMPc9sdGspyoql+nId8suhOxew7QUkR6hP3pQCvsmiMivYH7gL8AfVT1oZh0XwGKgaXAuap6dUi/HDhYVb+opW5ttFM2twjpVgKdsM7MdKlRW1U3qA2TnIE5m/+JpBmPOYNRwF7AKFVdWxvR0Pa96YG0vfe3iOyJPbSexn7nWrdb10I7o+U6XZszVKara2ezTG9OJtqX4v4AfTAPvA672VLHC8LfNli78J2EdjXgbuBnYbuQSLtqFrR/l9IO+40zrC3V7QG6YG8u+2EPoj1Ttqep3R97iP8BGBqxrVH4uyc2SufSyLFnsIcxQGegRw50C0lzFEVdtcN+0wxrS9SeYN9A4EnM6ZRQNUqo1vcZ8F+Ywx4DFEeOp+zb6v0N7Ax0S9PmdLTrXK7ranO613ob2pvZQobL9NY+9bpGICJFInIPFl71NmAKNnQNESnUKg++CngVa8u7OYxSaY0NsUJVKyJps6G9Y0o76NeqHW87tFVtjH4zEUm9jc4HngLew95Ui1O210Y7aAzC3vKexJpeTgd2FJECDZN31JokZmEPxyvCv64jtNWq6gK1ttRs6X6esldDacm2dkhTnmFtVVUVkSZio4MqVPUV7O34fey3bhu0a3ufDQTuwNrYOwBXisjQ8F2ptvat3t+qulRVP0nD5nS161qu62xzSFvrtvltaFfEWaa3SaY9SyY/WJX7NKBZ2D8L86aNImmuwYZS7Y3N6nwQ86L3UAfPmQfaV2MPjp5h/xTsIXwTdRw5gXXA3RO2OwGPEnnTxSbV3A/sFmx/Buvcuoc0a1651M0D7WuBCcBuYf9cYAlwY11+a6xmc1vYbgM8gNVI2kVszvj9nUvtem7ztcRUpreZtzi/PM2L1R/4TtiWaudGAXenzmGdhI/w7eFsBaQ5rCvPtfuT5lCyqHbY7w0sx5zNV8C0cNOOwIZIPkKopob0LYDW+aLbALSHRPfroHs08BDQIezfFrRGY3MyYrm/s6md5zanXaZrnddsiGznRWsNPIdVy8YCzcPxTe3hWJX8K2qYOUnd3kTzWbsub0fVtVtEzvXDHkbDw/4orFbSq65250q3AWin299To254+NyO9UP9BWuG+Bkwptr/Z9zmuLXz3OaM9wFs61Of+giaY23hF4XtTQG01NrOCrC22CnAoalzYMO+tJZthQ1Iuy7thdW1N8WNUdWZWEdganz2y9gNXhbRTtfuXOnmu3a6v/WW7rGPsXHy44E/q+oJWL/DYal/jPH+jls7n23OWsDEFDl1BCLy3yJyqIgUqw3zuxd4HJsUc6CIpIaoSfhhUiFYy1PH4VvDvlw7c9pNsKGK54d/HYy1a6aGt9W2ky4nuknV3oZuv5Suqq5X1amqOjH86/7YsEVisjk27STanCmy7gjE2EVEpgJnYp2ivxORtmrxUtYAf8d66Q8HewMOo2W+wZpM+qeOu3bGtQcHjXVYh2gLEXkF67i6UG2Wa73WTap2OvdY5H8PEZG3sFrKs7WxN5faSbQ5FjSL7VBUjcn9DvDHsN0Iazd7slraS7Ae/FbADpHjafWeu3attFtTNWKpGZE4K/VdN6nadbjHUn1SHYDvZfn+rpN2Em2O65MdEbtAv8KGux0KfB94KHJesJl7h0aOtcBm0c3EOko7uHbWtGcF7Y75optU7QzdY51yeH/XWjuJNsf9iV/ALta72My8c7CQAMOA+UC/SLrzgKmR/RFYQKX7iES1dO36q51Em/16u83Z0I77E7+AtYGdEdm/K1yos4C3wrECLFjV41RNmjkOGOja+aOdRJv9ervN2dCO+xO/gC3p1oSqNrXTgPFh+x3gorDdF3jUtfNXO4k2+/V2m7OhHfcn9lFDqrpGVddp1djYI7BIkWBx4/cRkWexafX/hM1D8bp2fmgn0eZcarvN2dPNtXbsZMvjEKLrYWNmU1H09sRGTBxCmp10rl3/tJNos19vtzkb2nF9sjmPoBIoApYBPYPn/AVQqaqvaXpx4127fmon0eZcarvNybA5PrLpdbAJUZXAa9jiGa7dQLWTaLNfb7c5Xz8SjMoKItIJOAO4RW0mpWs3UO0k2pxLbbc5GTbHRVYdgeM4jlP/qE/RRx3HcZwc4I7AcRwn4bgjcBzHSTjuCBzHcRKOOwLHcZyE447AcRwn4bgjcBzHSTj/D5Eej8r6L11yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run portfolio and see performance visually:\n",
    "df_portf = backtest_strategy(df_panel, df_mkt, prob_quantile = 0.8)\n",
    "\n",
    "ax = df_portf.plot(x = \"date\", y = ['portf_cumret', 'mkt_cumret'])\n",
    "ax.set_ylabel(\"Cumulative return\")\n",
    "ax.set_xlabel(\"\")\n",
    "ax.legend([\"Portfolio\", \"S&P 100\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comment: \n",
      "\n",
      " The model is helpful as it produced return that is higher than that of the stock market. \n",
      " Moreover, risk-adjusted return (Alpha) is 15.34% on annualized basis. \n",
      " Also, Sharpe ratio (reward-to-variability ratio) is almost twice higher for the suggested strategy\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>Portfolio</th>\n",
       "      <th>Market</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>return</td>\n",
       "      <td>0.2521</td>\n",
       "      <td>0.1099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>st dev</td>\n",
       "      <td>0.1891</td>\n",
       "      <td>0.1205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>sharpe</td>\n",
       "      <td>1.2275</td>\n",
       "      <td>0.7463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>alpha</td>\n",
       "      <td>0.1534</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>beta</td>\n",
       "      <td>0.8978</td>\n",
       "      <td>1.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Portfolio  Market\n",
       "return    0.2521  0.1099\n",
       "st dev    0.1891  0.1205\n",
       "sharpe    1.2275  0.7463\n",
       "alpha     0.1534  0.0000\n",
       "beta      0.8978  1.0000"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Portfolio assessment:\n",
    "# Assessment:\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "df_portf = df_portf.fillna(0)\n",
    "portf_Y = df_portf['portf_ret']\n",
    "portf_X = df_portf['mkt_ret']\n",
    "\n",
    "reg_p = LinearRegression().fit(portf_X.values.reshape(-1, 1), portf_Y.values.reshape(-1, 1))\n",
    "\n",
    "df_perf1 = pd.DataFrame( columns = ['Portfolio', 'Market'], \n",
    "                       index = ['return', 'st dev', 'sharpe', 'alpha', 'beta'])\n",
    "\n",
    "df_perf1.loc['return', 'Portfolio'] = \"{:.4f}\".format(df_portf['portf_ret'].mean()*51)\n",
    "df_perf1.loc['st dev', 'Portfolio'] = \"{:.4f}\".format(df_portf['portf_ret'].std()*(51**0.5))\n",
    "df_perf1.loc['sharpe', 'Portfolio'] = \"{:.4f}\".format((df_portf['portf_ret'].mean()*51 - .02)/(df_portf['portf_ret'].std()*(51**0.5)))\n",
    "df_perf1.loc['alpha', 'Portfolio'] = \"{:.4f}\".format(reg_p.intercept_[0]*51)\n",
    "df_perf1.loc['beta', 'Portfolio'] = \"{:.4f}\".format(reg_p.coef_[0][0])\n",
    "\n",
    "df_perf1.loc['return', 'Market'] = \"{:.4f}\".format(df_portf['mkt_ret'].mean()*51)\n",
    "df_perf1.loc['st dev', 'Market'] = \"{:.4f}\".format(df_portf['mkt_ret'].std()*(51**0.5))\n",
    "df_perf1.loc['sharpe', 'Market'] = \"{:.4f}\".format((df_portf['mkt_ret'].mean()*51 - .02)/ (df_portf['mkt_ret'].std()*(51**0.5)))\n",
    "df_perf1.loc['alpha', 'Market'] = \"{:.4f}\".format(0)\n",
    "df_perf1.loc['beta', 'Market'] = \"{:.4f}\".format(1)\n",
    "\n",
    "print(\"Comment: \\n\\n\", \"The model is helpful as it produced return that is higher than that of the stock market. \\n\", \n",
    "      \"Moreover, risk-adjusted return (Alpha) is 15.34% on annualized basis. \\n\", \n",
    "     \"Also, Sharpe ratio (reward-to-variability ratio) is almost twice higher for the suggested strategy\")\n",
    "\n",
    "df_perf1\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}