QuantitativeMethods_II.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline \n",
    "import math\n",
    "from scipy.stats import norm\n",
    "import scipy\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Common Probability Distributions\n",
    "A probability distribution specifies the probabilities of the possible outcomes of a random variable. Examples of random variables:\n",
    "- return on a stock\n",
    "- earnings per share \n",
    "\n",
    "Use cases examples:\n",
    "- the Black–Scholes–Merton option pricing model\n",
    "- the binomial option pricing model (price up or price down)\n",
    "- the capital asset pricing model\n",
    "\n",
    "We can view a probability distribution in two ways:\n",
    "#### Probability Function\n",
    "$P(X = x)$ is the probability that a random variable $X$ takes on the value $x$. (Note that capital X represents the random variable and lowercase $x$ represents a specific value that the random variable may take.) \n",
    "- For a discrete random variable, the shorthand notation for the probability function is $p(x) = P(X = x)$. \n",
    "- For continuous random variables, the probability function is denoted $f(x)$ and called the probability density function (pdf), or just the density.\n",
    "\n",
    "#### Cumulative Distribution Function (cdf)\n",
    "It gives the probability that a random variable X is less than or equal to a particular value x, P(X ≤ x). For both discrete and continuous random variables, the shorthand notation is $F(x) = P(X ≤ x).$ \n",
    "\n",
    "\n",
    "The cdf has two other characteristic properties:\n",
    "- The cdf lies between 0 and 1 for any $x$: $0 ≤ F(x) ≤ 1$.\n",
    "- As we increase $x$, the cdf either increases or remains constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete Random Variables\n",
    "The mathematical operation that corresponds to finding the area under the curve of a pdf $f(x)$ from $a$ to $b$ is the integral of $f(x)$ from $a$ to $b$:\n",
    "\n",
    "$$P(a\\le X\\le b)=\\int^b_af(x)dx$$\n",
    "<div style=\"text-align: right\"> Equation 2 </div>\n",
    "\n",
    "*Note: For any continuous random variable $X$, $P(a ≤ X ≤ b) = P(a < X ≤ b) = P(a ≤ X < b) = P(a < X < b)$, because the probabilities at the endpoints $a$ and $b$ are $0$.*\n",
    "### The Uniform Distribution \n",
    "Probability distribution that describes equally likely outcomes. <br>\n",
    "#### Discrete Uniform Distribution\n",
    "The pdf for a uniform random variable is:\n",
    "\\begin{equation}\n",
    "  f(x) =\n",
    "   \\begin{cases}\n",
    "      \\frac{1}{b−a} & \\text{for a<x<b}\\\\\n",
    "      0 & \\text{otherwise}\n",
    "   \\end{cases}       \n",
    "\\end{equation}\n",
    "\n",
    "### The Binomial Distribution（二项分布）\n",
    "When we make probability statements about a record of successes and failures, or about anything with binary outcomes, we often use the binomial distribution.<br>\n",
    "\n",
    "\n",
    "$X$~$B(n,p)$ read as “X has a binomial distribution with parameters n and p.”<br>\n",
    "\n",
    "The building block of the binomial distribution is the **Bernoulli random variable**. A trial (an event that may repeat) that produces one of two outcomes. Such a trial is a `Bernoulli trial`. If we let $Y= 1$ when the outcome is success and $Y = 0$ when the outcome is failure, then the probability function of the `Bernoulli random variable` $Y$ is\n",
    "$$p(1) = P(Y = 1) = p$$\n",
    "$$p(0) = P(Y = 0) = 1 – p$$\n",
    "where $p$ is the probability that the trial is a success. \n",
    "\n",
    "For $X$ distributed $B(n, p)$, the probability of $x$ successes in $n$ trials is given by\n",
    "$$p(x)=P(X=x)=\\binom{n}xp^x(1−p)^{n−x}=\\frac{n!}{(n−x)!x!}p^x(1−p)^{n−x}$$\n",
    "<div style=\"text-align: right\"> Equation 1 </div>\n",
    "\n",
    "When modeling stock price movements, the time periods can be made extremely small. Then the number of periods, $n$, can be large. The probability of exactly $x$ up moves in $n$ trials is:\n",
    "\n",
    "$$p(x)=P(X=x)=\\binom{n}xp^x(1−p)^{n−x}=\\frac{n!}{(n−x)!x!}p^x(1−p)^{n−x}$$\n",
    "\n",
    "Given independent trials, any sequence with $x$ up moves must occur with probability $p^x(1 − p)^{n–x}$. The binomial distribution is only symmetric if $p=0.50$. The mean of a binomial distribution is $np$ and the variance is $np(1−p)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean and Variance of Binomial Random Variables\n",
    "Because a single Bernoulli random variable, Y ~ B(1, p), takes on the value 1 with probability p and the value 0 with probability 1 − p, its mean or weighted-average outcome is p. Its variance is p(1 − p).9 A general binomial random variable, B(n, p), is the sum of n Bernoulli random variables, and so the mean of a B(n, p) random variable is np. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "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>X~B(n,p)</th>\n",
       "      <th>Mean</th>\n",
       "      <th>Variance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Bernoulli, B(1, p)</td>\n",
       "      <td>p</td>\n",
       "      <td>p(1-p)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Binomial, B(n, p)</td>\n",
       "      <td>np</td>\n",
       "      <td>np(1-p)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             X~B(n,p) Mean Variance\n",
       "0  Bernoulli, B(1, p)    p   p(1-p)\n",
       "1   Binomial, B(n, p)   np  np(1-p)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "des = {'X~B(n,p)': ['Bernoulli, B(1, p)', 'Binomial, B(n, p)'],\n",
    "      'Mean': ['p', 'np'], \n",
    "      'Variance': ['p(1-p)', 'np(1-p)']}\n",
    "df = pd.DataFrame(data=des)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example \n",
    "Consider a random variable distributed B(n = 5, p = 0.50):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>Number of Up Moves, x</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Number of Up Moves, x\n",
       "0                      0\n",
       "1                      1\n",
       "2                      2\n",
       "3                      3\n",
       "4                      4\n",
       "5                      5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {'Number of Up Moves, x':[0,1,2,3,4,5]}\n",
    "df = pd.DataFrame(data=data)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For convenience of calculation, turning the column into a list\n",
    "x = pd.to_numeric(df['Number of Up Moves, x'].to_list())\n",
    "x\n",
    "p = 0.50\n",
    "n = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 5.0, 10.0, 10.0, 5.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "# Counting number of ways to choose X from n\n",
    "results = []\n",
    "numerator = math.factorial(5) # why 5 instead of 6 (i.e., len(x))?\n",
    "\n",
    "for i in x:\n",
    "    denominator = math.factorial(5-i) * math.factorial(i)\n",
    "    result = numerator/denominator\n",
    "    results.append(result)\n",
    "print(results)\n",
    "\n",
    "# Append a column to represent the results\n",
    "df['Number of Possible Ways to Reach x Up Moves'] = results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating P of each way\n",
    "results2 = []\n",
    "for i in x: \n",
    "    result = p**i * (1-p)**(n-i)\n",
    "    results2.append(result)\n",
    "\n",
    "# Append the list to a new column    \n",
    "df['Probability for Each Way'] = results2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating the probability for each value of the number of up moves from the first column (using an alternative approach)\n",
    "df['Probability for x, p(x)'] = df['Number of Possible Ways to Reach x Up Moves'] * (1/sum(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.03125]\n",
      "[0.03125, 0.15625]\n",
      "[0.03125, 0.15625, 0.3125]\n",
      "[0.03125, 0.15625, 0.3125, 0.3125]\n",
      "[0.03125, 0.15625, 0.3125, 0.3125, 0.15625]\n",
      "[0.03125, 0.15625, 0.3125, 0.3125, 0.15625, 0.03125]\n"
     ]
    }
   ],
   "source": [
    "# Cumulating the entries in the fourth column, the cumulative distribution function\n",
    "cdf = []\n",
    "pdf = df['Probability for x, p(x)'].to_list()\n",
    "\n",
    "# Calculate the cdf\n",
    "for i in range(len(pdf)):\n",
    "    result3 = pdf[0:i+1]\n",
    "    print(result3)\n",
    "    cdf.append(sum(result3))\n",
    "    \n",
    "# Append the list to a new column\n",
    "df['F(x) = P(X ≤ x)'] = cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "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>Number of Up Moves, x</th>\n",
       "      <th>Number of Possible Ways to Reach x Up Moves</th>\n",
       "      <th>Probability for Each Way</th>\n",
       "      <th>Probability for x, p(x)</th>\n",
       "      <th>F(x) = P(X ≤ x)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.03125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.15625</td>\n",
       "      <td>0.18750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.31250</td>\n",
       "      <td>0.50000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.31250</td>\n",
       "      <td>0.81250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.15625</td>\n",
       "      <td>0.96875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>0.03125</td>\n",
       "      <td>1.00000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Number of Up Moves, x  Number of Possible Ways to Reach x Up Moves  \\\n",
       "0                      0                                          1.0   \n",
       "1                      1                                          5.0   \n",
       "2                      2                                         10.0   \n",
       "3                      3                                         10.0   \n",
       "4                      4                                          5.0   \n",
       "5                      5                                          1.0   \n",
       "\n",
       "   Probability for Each Way  Probability for x, p(x)  F(x) = P(X ≤ x)  \n",
       "0                   0.03125                  0.03125          0.03125  \n",
       "1                   0.03125                  0.15625          0.18750  \n",
       "2                   0.03125                  0.31250          0.50000  \n",
       "3                   0.03125                  0.31250          0.81250  \n",
       "4                   0.03125                  0.15625          0.96875  \n",
       "5                   0.03125                  0.03125          1.00000  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXiU9bn/8fedfYOEkLAlgQQSkrAFNEWBorIIAVrAU0Rt3Wp7eWylLriw+TvuGyBYe/QoRz1tj7WAWwuVVUXADQiSsIWEJCxJWBK2kADZv78/MvZECjJAZp6ZZ+7XdeUyM/NM+Iyaj+PM3PcjxhiUUkrZl5/VAZRSSrmWFr1SStmcFr1SStmcFr1SStmcFr1SStlcgNUBzhYTE2MSExOtjqGUUl5l8+bNR4wxsee6zeOKPjExkezsbKtjKKWUVxGRfee7TV+6UUopm9OiV0opm9OiV0opm9OiV0opm9OiV0opm9OiV0opm9OiV0opm/O4z9Er9UOqaur5y4b9nK5tsDqKW3WPjWBC/y6IiNVRlBfSolde5YXlu/jLhv34Ut99d8qIyLBAhqV2sDaM8kpa9MprFJZXsXBTCXcOTuSJ8b2tjuM2dQ1NjJq/lueX5TE0OYYAf33FVV0c/TdGeY0Xlu8iLNCf+0akWB3FrYIC/JiWlUbB4Wre31xqdRzlhbTolVf4qugIn+SV89thyUSHB1kdx+2y+nTiym7teGl1Aad87P0Jdfm06JXHa2oyPLcsj7ioUH45JNHqOJYQEWaNS6eiqpYF64qtjqO8jBa98nh/zy1je9lJHhmdSkigv9VxLHNF13aM69eZBeuKOXyyxuo4yos4VfQikiUi+SJSKCLTz3H7PSKyTURyROQLEenV4rYZjvvli8jo1gyv7K+mvpE5K/LpGxfJ+IwuVsex3LTRaTQ0NTFvVYHVUZQXuWDRi4g/8CowBugF3NKyyB3eNcb0Ncb0B2YD8xz37QXcDPQGsoDXHD9PKae8/eUeDlTWMHNsOn5+PvSZyvPo2j6MOwYlsnhzCXkHT1odR3kJZ57RDwQKjTHFxpg6YCEwoeUBxpiW/8aFA45P/jIBWGiMqTXG7AEKHT9PqQs6Wl3La2uKGJnekUE92lsdx2NMGZ5M25BAnl++y+ooyks4U/RxQEmLy6WO675HRO4VkSKan9Hfd5H3vVtEskUku6KiwtnsyuZ+/+luztQ3Mn1MmtVRPEpUWBC/G57MuoIK1hbo74u6sFZ7M9YY86oxpgcwDXjsIu+7wBiTaYzJjI095ykPlY8pqqjmLxv28/OBXUnuEGF1HI9z26BudI0O4/lleTQ2mQvfQfk0Z4q+DEhocTnecd35LAQmXuJ9lQKah6NCA/25f6RvDUc5KzjAn2lZaew6VMUHOkSlLsCZot8EpIhIkogE0fzm6pKWB4hIy9/GccBux/dLgJtFJFhEkoAUYOPlx1Z2tqH4KKt3HuY31/UgJiLY6jgea2zfTgzoGsXcVfmcrtMhKnV+Fyx6Y0wDMAVYCeQBi40xO0TkKREZ7zhsiojsEJEcYCpwh+O+O4DFwE5gBXCvMabRBY9D2cR3w1GdI0O4a0iS1XE8mogwa2w65VW1vLl+j9VxlAdzaqmZMWYZsOys6/6jxff3/8B9nwWevdSAyrcs3XqA3NJKXroxg9Ag/STuhWQmRjOmTydeX1vEzQMT6NAmxOpIygPpZKzyGDX1jcxekU+vzm25YcC/fDhLnce0rDTqGpqYv3r3hQ9WPkmLXnmMP321l7ITZ3hsnA5HXYzEmHBuG9SNRZv2U3C4yuo4ygNp0SuPcOxUHf+5ppDhaR0YnBxjdRyvc9/wFMKDA3h+WZ7VUZQH0qJXHuGVT3dzqraBGTocdUnahTcPUa3Jr+CL3UesjqM8jBa9styeI6d455t93DywKykd21gdx2vdPiiR+HahPKtDVOosWvTKci8u30VwgB8P6HDUZQkJ9OfRrDTyDp7koy06l6j+jxa9stSmvcdYseMQ91zbQz8a2Ap+2q8zGQlRzF2Zz5k6HVlRzbTolWWMMTz7cR4d2wbz66HdrY5jC98NUR06WcNbX+iZqFQzLXplmY+3HSSn5AQPjUrV4ahWNDApmlG9OvJfnxdRUVVrdRzlAbTolSVqGxp5ccUu0jq14WdXxFsdx3amj0mjtqGJ33+qZ6JSWvTKIv/79T5Kjp1h1rh0/HU4qtV1j43gF1d15a8bSygs1yEqX6dFr9zuxOk6Xvl0N9f2jGVoip5/wFXuG5FCWKA/L+iZqHyeFr1yuz98Vkh1bQMzx6ZbHcXW2kcE89thyXySV85XRTpE5cu06JVb7Tt6ij9/vZfJmQmkdtLhKFf75ZBE4qJCeW5ZHk06ROWztOiVW81ekU+Anx9Tr+9pdRSfEBLozyOjU9ledpK/5+oQla/Soldus3nfcT7edpB/v7Y7HdrqcJS7jM/oQt+4SOasyKemXoeofJEWvXKL5uGonXRoE8zd1+hwlDv5+Qkzx6ZzoLKGt7/UM1H5Ii165RbLtx/i2/0neGhUT8KCnDqxmWpFg3q0Z2R6R15bU8TRah2i8jVa9Mrl6hqaeHHFLlI7tmHSlQlWx/FZ08ekcaa+kd9/qmei8jVa9Mrl3vlmH/uOnmbG2DQdjrJQcocIbhmYwLsb9lNUUW11HOVGWvTKpSpP1/PKZ7sZmhLDtT11OMpqD4zsSUigPy/qEJVP0aJXLvXq54VUnqlnxph0RPTZvNViIoL5zXU9WLXzMBuKj1odR7mJFr1ymZJjp/njl3uZdEU8vbq0tTqOcrhrSBKdI0N0iMqHaNErl5m9Mh8/P3hoVKrVUVQLoUH+PDwqldzSSpZuPWB1HOUGWvTKJbbsP87S3APcPbQ7nSJ1OMrT3DAgjl6d2zJbh6h8gha9anXGGJ5blkdMRDB3X9vD6jjqHPz8hMfGpVN24gx/+mqv1XGUizlV9CKSJSL5IlIoItPPcftUEdkpIltF5FMR6dbitkYRyXF8LWnN8MozrdxxmE17jzP1+p5EBOtwlKcanBzD8LQO/OeaQo6dqrM6jnKhCxa9iPgDrwJjgF7ALSLS66zDtgCZxph+wPvA7Ba3nTHG9Hd8jW+l3MpD1TU08cLyPFI6RDA5U88c5elmjEnjVG0Dr+gQla0584x+IFBojCk2xtQBC4EJLQ8wxqwxxpx2XPwG0N9wH/Xuhn3sdQxHBfjrK4OeLqVjG276UVfe+WYfe46csjqOchFnfhPjgJIWl0sd153Pr4DlLS6HiEi2iHwjIhPPdQcRudtxTHZFRYUTkZQnOllTz+8/3c3gHu0ZltrB6jjKSQ9en0JQgB+zV+gQlV216lMuEbkVyATmtLi6mzEmE/g58LKI/Mu7c8aYBcaYTGNMZmysTk96q9fWFHHiTD0zx+pwlDfp0CaEe67twfLth8jee8zqOMoFnCn6MqDlJqp4x3XfIyIjgVnAeGPMP9fjGWPKHH8tBj4HBlxGXuWhSo+f5u0v93DDgDj6xEVaHUddpF8PTaJj22Ce+TgPY3SIym6cKfpNQIqIJIlIEHAz8L1Pz4jIAOANmku+vMX17UQk2PF9DDAE2Nla4ZXnmLsyHwEe1uEorxQWFMBDo1LJKTnBx9sOWh1HtbILFr0xpgGYAqwE8oDFxpgdIvKUiHz3KZo5QATw3lkfo0wHskUkF1gDvGCM0aK3ma2lJ/hbzgF+PTSJLlGhVsdRl+hnV8ST1qkNL67YRW2DDlHZiVMfcjbGLAOWnXXdf7T4fuR57vcV0PdyAirP1nzmqDzahwdxjw5HeTV/P2HWuHRue2sj//v1Pn49VM8EZhf6+Td1WT7JK2fDnmM8cH1P2oQEWh1HXaahKbFc2zOWVz7dzYnTOkRlF1r06pLVNzbx/PI8esSGc/OP9MxRdjFzbDrVtQ384bNCq6OoVqJFry7Zwo37Ka44xYwx6QTqcJRtpHZqw+TMBP789V72HdUhKjvQ3051SU7W1DP/k91clRTNiHQdjrKbqdf3JMDPj9kr8q2OolqBFr26JK9/XsSxU3XMGqfDUXbUoW0Id1/TnY+3HWTzvuNWx1GXSYteXbQDJ87w1hd7mNi/C/3io6yOo1zk7mu6E9smmOeW6RCVt9OiVxdt7qp8DPDwaB2OsrPw4AAeur4nm/cdZ8X2Q1bHUZdBi15dlO1llXy0pYy7hiQR3y7M6jjKxW7MTCC1YxteWLGLuoYmq+OoS6RFr5z23XBUVGggvx2mw1G+wN9PmDE2jX1HT/PON/usjqMukRa9ctqa/HK+Lj7KAyN70laHo3zGtT1jGZoSwyuf7abydL3VcdQl0KJXTmlobOK5ZbtIignn51d1tTqOciMRYcaYdCrP1PPq5zpE5Y206JVTFmWXUFhezfQxaToc5YN6dWnLpCvi+eOXeyk5dvrCd1AeRX9j1QVV1zYwf3UBAxOjGdWro9VxlEUeGpWKnx/MXqlDVN5Gi15d0BtrizhSXcdMHY7yaZ0iQ7h7aHeW5h5gy34dovImWvTqBx2qrOG/1xfz04wu9E/Q4Shfd/e1PYiJCNIhKi+jRa9+0Eur8mlqgkd1OEoBEcEBPHh9TzbtPc6qnYetjqOcpEWvzmvngZO8/20pdw5JJCFah6NUs5syE0juEMELy3dR36hDVN5Ai16dkzGG55blERkayL3XJVsdR3mQAH8/Zo5NY8+RU7y7Yb/VcZQTtOjVOa0tqOCLwiPcNzyFyDAdjlLfNyy1A4N7tOflTwo4WaNDVJ5Oi179i+bhqDy6tQ/j1qu7WR1HeSARYebYdE6cqee1NUVWx1EXoEWv/sX7m0spOFzN9Kw0ggL0XxF1bn3iIrlhQBxvf7mH0uM6ROXJ9LdYfc+p2gZeWl3Ald3akdWnk9VxlId7eFQqAszVISqPpkWvvmfBumIqqmr1zFHKKV2iQvn10CT+lnOAraUnrI6jzkOLXv3T4ZM1LFhXzLh+nbmiazur4ygvcc+1PWgfHsSzH+sQlafSolf/NG9VAQ1NTUwbnWZ1FOVF2oQE8sDIFDbsOcYneeVWx1HnoEWvANh16CSLN5dw+6BEurbX4Sh1cW4e2JXuseE8vzxPh6g8kFNFLyJZIpIvIoUiMv0ct08VkZ0islVEPhWRbi1uu0NEdju+7mjN8Kr1PL9sF22CA/jdcB2OUhcv0N+PGWPSKa44xcJNJVbHUWe5YNGLiD/wKjAG6AXcIiK9zjpsC5BpjOkHvA/Mdtw3GngcuAoYCDwuIvrir4dZV1DB2oIK7huRQlRYkNVxlJcamd6Bq5KieXl1AVU6ROVRnHlGPxAoNMYUG2PqgIXAhJYHGGPWGGO++yDtN0C84/vRwGpjzDFjzHFgNZDVOtFVa2hsal51kBAdym2DdDhKXToRYda4dI6equP1tTpE5UmcKfo4oOX/i5U6rjufXwHLL+a+InK3iGSLSHZFRYUTkVRr+eDbUnYdqmJaVhrBAf5Wx1Ferl98FBP7d+HN9Xs4cOKM1XGUQ6u+GSsitwKZwJyLuZ8xZoExJtMYkxkbG9uakdQPOF3XwEur8umfEMW4vp2tjqNs4uHRqRhg7iodovIUzhR9GZDQ4nK847rvEZGRwCxgvDGm9mLuq6zx5vo9HD5Zy2M6HKVaUXy7MO4aksRHW8rYXlZpdRyFc0W/CUgRkSQRCQJuBpa0PEBEBgBv0FzyLT9IuxIYJSLtHG/CjnJcpyxWXlXD62uLGNOnE5mJ0VbHUTbz22E9iAoN1CEqD3HBojfGNABTaC7oPGCxMWaHiDwlIuMdh80BIoD3RCRHRJY47nsMeJrm/1hsAp5yXKcsNn/1buoampiWpcNRqvW1DQnkgZE9+br4KGvydYjKauJp/7XNzMw02dnZVsewtYLDVWS9vI7bByXyxPjeVsdRNlXf2MSo+evw9xNW3D+UAH+dz3QlEdlsjMk81236d94HPb8sj/DgAO4bkWJ1FGVjgf5+TMtKo7C8mkXZOkRlJS16H/Nl4RHW5FcwZVgy0eE6HKVca3TvjvwosR3zV++murbB6jg+S4vehzQ1GZ79OI+4qFDuGJxodRzlA747E9WR6loW6BCVZbTofchHW8rYefAkj2alEhKow1HKPQZ0bcdPM7qwYH0xhyprrI7jk7TofcSZukbmrsonIz6Sn/brYnUc5WMeHZ1KUxO8pENUltCi9xFvf7mHg5U1zBybjp+fDkcp90qIDuPOIYm8/20pOw+ctDqOz9Gi9wEVVbW8tqaQUb06clX39lbHUT7q3uuSiQwN5LllOkTlblr0PuD3nxZQ29DE9DE6HKWsExkWyH3DU/ii8AhrC3R5oTtp0dtcYXkVf91Ywi+u6kr32Air4ygfd+vV3ejWPoznluXRoGeichstept7YfkuwgL9dThKeYSgAD+mZ6VRcLia9zeXWh3HZ2jR29hXRUf4JK+c3wzrQfuIYKvjKAVAVp9OXNmtHS+tLuCUDlG5hRa9TTU5zhzVJTKEu4YkWR1HqX/6boiqoqqWBeuKrY7jE7TobWpJ7gG2l53kER2OUh7oym7tGNe3MwvWFVN+UoeoXE2L3oZq6huZszKfPnFtmZDxQ2d9VMo6j2al0tDUxLzVBVZHsT0tehv6ny/3UnbijA5HKY/WrX04tw9KZHF2CbsO6RCVK2nR28zR6ubhqJHpHRjcI8bqOEr9oN8NTyYiOIDnl+2yOoqtadHbzCuf7uZ0faMORymvEBUWxH0jUlhbUME6HaJyGS16GymqqOYvG/Zzy8AEkju0sTqOUk65bVA3EqJDeW5ZHo1NuhrBFbTobeTF5bsICfTngZE9rY6ilNOCA/yZlpXGrkNVfPCtDlG5gha9TWwoPsqqnYf5zXU9iNHhKOVlxvXtTP+EKF5alc/pOh2iam1a9Dbw3XBUp7Y6HKW8k4jw2Lh0Dp+s5c31e6yOYzta9DawdOsBcksreXh0KqFBOhylvFNmYjRZvTvx+toiyqt0iKo1adF7uZr6RmavyKdX57bcMECHo5R3mzYmjbqGJl7+ZLfVUWxFi97L/fnr5uGoWePS8dfhKOXlkmLCufXqbizcuJ/dh6usjmMbWvRe7Gh1LX/4rJBhqbEMSdbhKGUP941IITw4gGc+1jNRtRYtei9ljGHmR9uorW9i5th0q+Mo1Wqiw4N4YGRP1hZU6M76VuJU0YtIlojki0ihiEw/x+3XiMi3ItIgIpPOuq1RRHIcX0taK7iv++DbMlbuOMzDo3uS0lGHo5S93Dk4kauSonly6U5Kjp22Oo7Xu2DRi4g/8CowBugF3CIivc46bD9wJ/DuOX7EGWNMf8fX+MvMq4CSY6d5YskOBiZF86sfd7c6jlKtzt9PeGlyBgAPvZerE7OXyZln9AOBQmNMsTGmDlgITGh5gDFmrzFmK6AngXSxxibDQ+/lAvDSjRn6Bqyyrfh2YTz+015s3HOMt77QE5RcDmeKPg4oaXG51HGds0JEJFtEvhGRiec6QETudhyTXVGhi41+yFtfFLNxzzEe/2kvEqLDrI6jlEtNujKe0b07MndlAXkHdZXxpXLHm7HdjDGZwM+Bl0Wkx9kHGGMWGGMyjTGZsbGxbojknXYdOsnclQWM6tWRSVfGWx1HKZcTEZ67oS9tQwN5cFEOtQ2NVkfySs4UfRmQ0OJyvOM6pxhjyhx/LQY+BwZcRD7lUNvQyIOLcmkbGsDz/9YXEX3JRvmG9hHBvPizvuw6VMX81TpIdSmcKfpNQIqIJIlIEHAz4NSnZ0SknYgEO76PAYYAOy81rC+bv3o3eQdP8sK/9aO9Li1TPmZEekduGZjAG+uK2LjnmNVxvM4Fi94Y0wBMAVYCecBiY8wOEXlKRMYDiMiPRKQUuBF4Q0R2OO6eDmSLSC6wBnjBGKNFf5E27T3GG+uKuPlHCYzs1dHqOEpZ4rFxvUhoF8ZD7+VQXasbLi+GeNrkWWZmpsnOzrY6hseorm1gzO/XIQjL7h9KRHCA1ZGUskz23mNMfuNrbrwygRcn9bM6jkcRkc2O90P/hU7Gerinl+6k7PgZ5k3O0JJXPi8zMZp7ru3BouwSVu04ZHUcr6FF78FW7zzMouwS/v3aHmQmRlsdRymP8MDInvTq3JYZH27jSHWt1XG8gha9hzpSXcv0D7aS3rktD+qpAZX6p6AAP+bf1J+q2gamf7BNF585QYveAxljmPHhNqpqGnj5pv4EBeg/JqVaSu3UhkdHp/JJ3mHey9bFZxeiDeKB3ttcyuqdh3lkdCqpnXRhmVLncteQJAZ1b8+TS3fo4rML0KL3MCXHTvPkkh1c3T2aX/1Yz/+q1Pn4+QlzJ2fgJ8LUxTm6+OwHaNF7kMYmw0OLc/ETYe6NGfjpwjKlflBcVChPjO/Npr3H+e/1uvjsfLToPcib64vZuPcYj4/vTXw7XVimlDP+7Yo4xvTpxEur8tl5QBefnYsWvYfIO3iSl1YVkNW7Ez+7Qk/yrZSzRIRnb+hLZGgQUxfr4rNz0aL3AM0Ly3JoGxrIc7qwTKmLFh0exOxJzYvP5q0qsDqOx9Gi9wDzVhWw61AVsyf1JTo8yOo4Snml4Wkd+flVXVmwvphvio9aHcejaNFbbEPxURasL+aWgV0ZnqYLy5S6HLPGptM1OoyHFudSVVNvdRyPoUVvoaqaeh56L5eu0WE8Ni7d6jhKeb3w4ADmTe7PwcozPLVUF+V+R4veQk8t3cmBE80Ly8J1YZlSreLKbu347XXJvLe5lJW6+AzQorfMqh2HeG9zKb+5rgdXdtOFZUq1pvtGpNC7S/Pis4oqXXymRW+BiqpaZny4jd5d2nL/CF1YplRrCwrw4+Wb+lNd28CMD7f6/OIzLXo3a15YtpWq2gbm68IypVwmpWMbpmWl8UleOYs2lVgdx1LaMm62OLuET/LKeXR0Kj076sIypVzpl4MTGdyjPU//Yyf7j/ru4jMtejfaf/Q0Ty3dyaDu7blriC4sU8rV/PyEOY69Ub68+EyL3k0amwxTF+c0LyybrAvLlHKXuKhQnprQm+x9x3ljXZHVcSyhRe8mC9YVk73vOE9O6E1cVKjVcZTyKRP7xzGub2fmry5gx4FKq+O4nRa9G+w4UMm81fmM7duJGwbowjKl3E1EeGZiH9qFBfHgohxq6n1r8ZkWvYvV1DcydVEuUWFBPDtRF5YpZZV24UG8OKkfBYereWlVvtVx3EqL3sXmrS4g/3AVsyf1o50uLFPKUsNSO3Dr1V1584s9fF3kO4vPtOhd6Jvio/z3+mJ+cVVXhqV2sDqOUgqYOTadxPbhPPxeLid9ZPGZFr2LVNXU89DiXLpFhzFLF5Yp5THCggJ4aXIGByvP8OQS31h85lTRi0iWiOSLSKGITD/H7deIyLci0iAik8667Q4R2e34uqO1gnu6J5fu5GDlGebd1J+wIF1YppQnuaJrO6YMS+aDb0tZsf2g1XFc7oJFLyL+wKvAGKAXcIuI9DrrsP3AncC7Z903GngcuAoYCDwuIu0uP7ZnW7H9EO9vLuXeYclc0dX2D1cpr/S7ESn0jYtkxofbKK+qsTqOSznzjH4gUGiMKTbG1AELgQktDzDG7DXGbAWazrrvaGC1MeaYMeY4sBrIaoXcHqu8qoaZH22jT1xb7huRYnUcpdR5BPr7Mf+mDE7XNTL9g222XnzmTNHHAS03ApU6rnPG5dzX6xhjmP7BNk7VNjB/cn8C/fUtEKU8WXKHNkwfk8Znu8r560b7Lj7ziCYSkbtFJFtEsisqKqyOc8kWbirhs13lTMtKI0UXlinlFe4YlMiQ5PY88/FO9h09ZXUcl3Cm6MuAhBaX4x3XOcOp+xpjFhhjMo0xmbGxsU7+aM+y7+gpnv7HToYkt+fOwYlWx1FKOcnPT5h7YwYBfsKDi3JoaDz7FWjv50zRbwJSRCRJRIKAm4ElTv78lcAoEWnneBN2lOM6W2leWJaLv58wZ5IuLFPK23SODOXpiX34dv8J3lhXbHWcVnfBojfGNABTaC7oPGCxMWaHiDwlIuMBRORHIlIK3Ai8ISI7HPc9BjxN838sNgFPOa6zldfXFrF533GentCHLrqwTCmvND6jC+P6NS8+215mr8Vn4mnvNGdmZprs7GyrYzhte1klN7z2JaN6deI/fz5Ad9ko5cVOnK5j1Px1RIYGsvR3PyYk0N/qSE4Tkc3GmMxz3eYRb8Z6q5r6Rh5clEO7sCCemdhHS14pLxcVFsScGzPYXV7NnJX2WXymRX8Z5q7MZ3d5tS4sU8pGru0Zy+2DuvHWF3v4qvCI1XFahRb9Jfqq6AhvfrGHW6/uynW6sEwpW5kxJp3uMc2LzyrPeP/iMy36S3Cypp6HF+eSFBPOzLG6sEwpuwkN8mfeTf05XFXLk0t2WB3nsmnRX4InluzgcFUt8yZn6MIypWyqf0IUU4Yl8+GWMpZt8+7FZ1r0F2n5toN8+G0Z9w5LZoAuLFPK1qYMT6ZffCQzP9pG+UnvXXymRX8Ryk82LyzrFx/J74YnWx1HKeVigf5+zJvcnzN1jTz6wVavXXymRe8kYwzTPtjK6bpG5unCMqV8RnKHCGaOTefz/Are3bjf6jiXRNvKSe9u3M+a/ApmjEkjuUOE1XGUUm5029XdGJoSwzP/yGPPEe9bfKZF74S9R07xzD/y+HFyDLcPSrQ6jlLKzfwce6yCAvyYutj7Fp9p0V9AQ2MTDy7OIdBfmHNjP11YppSP6hQZwtMT+7Bl/wn+6/Miq+NcFC36C3h9bRFb9p/g6Yl96BypC8uU8mXjM7rw04wu/P7T3Wwr9Z7FZ1r0P2B7WSUvf7Kbn/TrzIT+tj0xllLqIjw9oTcxEcE8uDiHmvpGq+M4RYv+PGrqG3lgUQ7tI5oXlimlFHy3+KwfheXVvLhil9VxnKJFfx6zV+RTWF7N7EkZRIXpwjKl1P8ZmhLLnYMT+Z8v9/KlFyw+06I/hy8Lj/D2l3u4fVA3ru3pnac2VEq51rSsNLrHesfiMy36s1Seqefh93LpHhPOjDG6sEwpdW6hQf7Mn9yf8qpaHv/7dqvj/CAt+rM8sWQH5VW1zL+pP6FB3nN2GaWU+2UkRHHf8BT+lnOAf2w9YHWc89Kib+HjrQf5aEsZU3WH3PUAAAhZSURBVIYlk5EQZXUcpZQXuHdYDzISopj10XYOe+jiMy16h/KTNcz62zYy4iOZogvLlFJOCvD3Y/7kDGobGnnkfc9cfKZFT/PCskfe30pNfSPzbtKFZUqpi9M9NoJZY9NZV1DBOxs8b/GZNhrwzob9rC2oYMaYdHrE6sIypdTFu/XqblzTM5ZnP95JcUW11XG+x+eLfs+RUzz3cR5DU2K47epuVsdRSnkpEWHOpH4EB/jz4OJcj1p85tNF39DYxIOLcggK8GPOpAxdWKaUuiwd24bw7A19yC05wWsetPjMp4v+tc+LyClpXljWKTLE6jhKKRv4Sb8uTOjfvPhsa+kJq+MAPlz020oreeXT3YzP6ML4jC5Wx1FK2chT4/vQoU0wDy7yjMVnPln0zQvLthATEczTE3RhmVKqdUWGBTL3xgyKKk7xwnLrF585VfQikiUi+SJSKCLTz3F7sIgscty+QUQSHdcnisgZEclxfL3euvEvzQvLd1FUcYo5N/YjMizQ6jhKKRsakhzDnYMT+eNXe1m/u8LSLBcsehHxB14FxgC9gFtEpNdZh/0KOG6MSQbmAy+2uK3IGNPf8XVPK+W+ZF/sPsIfv9rLnYMTGZqiC8uUUq4zfUwaPWLDeeS9rVSetm7xmTPP6AcChcaYYmNMHbAQmHDWMROAPzm+fx8YISIe9xGWytOOhWWx4UzLSrM6jlLK5kIC/Xn5pgEcqa7l/1m4+MyZoo8DSlpcLnVcd85jjDENQCXQ3nFbkohsEZG1IjL0XH+AiNwtItkikl1R4br/xfmPJds5Ul3Ly7qwTCnlJn3jI7l/RApLcg+wNNeaxWeufjP2INDVGDMAmAq8KyJtzz7IGLPAGJNpjMmMjXXNyylLcw/w95wD/G54Cv3idWGZUsp9fnNdDwZ0jeKxv23nUKX7F585U/RlQEKLy/GO6855jIgEAJHAUWNMrTHmKIAxZjNQBPS83NAX61BlDY/9bTsZCVHcO6yHu/94pZSPC/D3Y97k/tQ1NPHI+7luX3zmTNFvAlJEJElEgoCbgSVnHbMEuMPx/STgM2OMEZFYx5u5iEh3IAUobp3ozjHG8OgHW6ltaGT+5AwCdGGZUsoCSTHhzBqXzvrdR/jfb/a59c++YOs5XnOfAqwE8oDFxpgdIvKUiIx3HPYW0F5ECml+iea7j2BeA2wVkRya36S9xxhzrLUfxA9555t9rCuoYNbYdLrrwjKllIV+cVVXrkuN5blleRS5cfGZeNru5MzMTJOdnd0qP6uooppxr6xnYFJ7/vTLH+GBHwRSSvmY8pM1jHp5Hd2iw3j/N4NbbS26iGw2xmSe6zbbvo5R39jE1EU5hAT6M2dSPy15pZRH6NA2hOdu6EtuaSWvril0y59p26J/dU0huaWVPDOxDx3b6sIypZTnGNu3MzcMiOMPnxWSU+L6xWe2LPrckhP84bNCJvbvwk/66cIypZTneWJ8bzq2CWbqohzO1Ll28Zntiv5MXSMPLs6hQ5tgntSFZUopDxUZ2rz4rPjIKV5YnufSP8t2Rf/C8jyKK04x98YMIkN1YZlSynMNTo7hriFJ/Onr5k8Huoqtin797gr+9PU+fjkkkSHJMVbHUUqpC3o0K5WUDhE88n4uJ07XueTPsE3Rnzhdx8Pv5ZLcIUIXlimlvEZIoD/zb+rP0eo6/t/fd7jkzwhwyU+1QEOToW9cFPePSCEkUBeWKaW8R5+4SB4alcqZ+kaamkyrn7/a1gNTSinlK3xyYEoppVQzLXqllLI5LXqllLI5LXqllLI5LXqllLI5LXqllLI5LXqllLI5LXqllLI5jxuYEpEK4HJOqBgDHGmlON7C1x6zrz1e0MfsKy7nMXczxsSe6waPK/rLJSLZ55sOsytfe8y+9nhBH7OvcNVj1pdulFLK5rTolVLK5uxY9AusDmABX3vMvvZ4QR+zr3DJY7bda/RKKaW+z47P6JVSSrWgRa+UUjZnm6IXkSwRyReRQhGZbnUeVxORt0WkXES2W53FXUQkQUTWiMhOEdkhIvdbncnVRCRERDaKSK7jMT9pdSZ3EBF/EdkiIv+wOou7iMheEdkmIjki0qpnX7LFa/Qi4g8UANcDpcAm4BZjzE5Lg7mQiFwDVAN/Nsb0sTqPO4hIZ6CzMeZbEWkDbAYm2vyfswDhxphqEQkEvgDuN8Z8Y3E0lxKRqUAm0NYY8xOr87iDiOwFMo0xrT4kZpdn9AOBQmNMsTGmDlgITLA4k0sZY9YBx6zO4U7GmIPGmG8d31cBeUCctalcyzSrdlwMdHx5/7OzHyAi8cA44E2rs9iFXYo+DihpcbkUmxeArxORRGAAsMHaJK7neBkjBygHVhtj7P6YXwYeBZqsDuJmBlglIptF5O7W/MF2KXrlQ0QkAvgAeMAYc9LqPK5mjGk0xvQH4oGBImLbl+pE5CdAuTFms9VZLPBjY8wVwBjgXsfLs63CLkVfBiS0uBzvuE7ZjON16g+AvxhjPrQ6jzsZY04Aa4Asq7O40BBgvOP16oXAcBF5x9pI7mGMKXP8tRz4iOaXpFuFXYp+E5AiIkkiEgTcDCyxOJNqZY43Jt8C8owx86zO4w4iEisiUY7vQ2n+wMEua1O5jjFmhjEm3hiTSPPv8WfGmFstjuVyIhLu+IABIhIOjAJa7RN1tih6Y0wDMAVYSfMbdIuNMTusTeVaIvJX4GsgVURKReRXVmdygyHAbTQ/y8txfI21OpSLdQbWiMhWmp/QrDbG+MxHDn1IR+ALEckFNgIfG2NWtNYPt8XHK5VSSp2fLZ7RK6WUOj8teqWUsjkteqWUsjkteqWUsjkteqWUsjkteqWUsjkteqWUsrn/D0ClKLl9TW50AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Observe how the change of p(x) influence the skewness of the curve\n",
    "plt.plot(df['Probability for x, p(x)']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "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>Block Brokers</th>\n",
       "      <th>Profitable Trades</th>\n",
       "      <th>Losing Trades</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BB001</td>\n",
       "      <td>3</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BB002</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Block Brokers  Profitable Trades  Losing Trades\n",
       "0         BB001                  3              9\n",
       "1         BB002                  5              3"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {'Block Brokers': ['BB001', 'BB002'], \n",
    "       'Profitable Trades': [3, 5], \n",
    "       'Losing Trades': [9, 3]}\n",
    "df = pd.DataFrame(data=data)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are paying a fair price on average in your trades with a broker, what should be the probability of a profitable trade?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under the assumption that the prices of trades were fair,\n",
    "1. calculate the probability of three or fewer profitable trades with broker BB001.\n",
    "2. calculate the probability of five or more profitable trades with broker BB002."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: using binomial distribution to evaluate the performance of an investment manager<br>\n",
    "\n",
    "Expected to keep portfolio return within a band of 75 basis points (bps) of the benchmark’s return, on a quarterly basis, portfolio return is within the 75 bps band 90 percent of the time.\n",
    "\n",
    "Question: The manager meets the objective in six out of eight quarters. How does the manager’s record precisely relate to your expectation of a 90 percent success rate and the sample size, 8 observations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of known variables\n",
    "n = 8\n",
    "p = 0.90"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "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>Number of Success (&gt;=75 bps)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Number of Success (>=75 bps)\n",
       "0                             0\n",
       "1                             1\n",
       "2                             2\n",
       "3                             3\n",
       "4                             4\n",
       "5                             5\n",
       "6                             6\n",
       "7                             7\n",
       "8                             8"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# find the probability that, given an assumed true or underlying success rate of 90 percent, performance could be as bad as or worse than that delivered.\n",
    "data = {'Number of Success (>=75 bps)':[0,1,2,3,4,5,6,7,8]}\n",
    "df = pd.DataFrame(data=data)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 8.0, 28.0, 56.0, 70.0, 56.0, 28.0, 8.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "results = []\n",
    "numerator = math.factorial(n) \n",
    "x = df['Number of Success (>=75 bps)'].to_list()\n",
    "# Counting the number of ways to reach the corresponding sucesses\n",
    "for i in x:\n",
    "    denominator = math.factorial(n-i) * math.factorial(i)\n",
    "    result = numerator/denominator\n",
    "    results.append(result)\n",
    "print(results)\n",
    "\n",
    "# Append a column to represent the results\n",
    "df['Number of Possible Ways to Reach x'] = results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Probability of each way \n",
    "results = []\n",
    "for i in x:\n",
    "    p_x = p**i\n",
    "    p_y = (1-p)**(n-i)\n",
    "    result = p_x * p_y\n",
    "    results.append(result)\n",
    "df['Probability for Each Way'] = results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \n",
    "df['Probability for x, p(x)'] = df['Number of Possible Ways to Reach x'] * df['Probability for Each Way']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "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>Number of Success (&gt;=75 bps)</th>\n",
       "      <th>Number of Possible Ways to Reach x</th>\n",
       "      <th>Probability for Each Way</th>\n",
       "      <th>Probability for x, p(x)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.000000e-08</td>\n",
       "      <td>1.000000e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>8.0</td>\n",
       "      <td>9.000000e-08</td>\n",
       "      <td>7.200000e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>28.0</td>\n",
       "      <td>8.100000e-07</td>\n",
       "      <td>2.268000e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>56.0</td>\n",
       "      <td>7.290000e-06</td>\n",
       "      <td>4.082400e-04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>70.0</td>\n",
       "      <td>6.561000e-05</td>\n",
       "      <td>4.592700e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>56.0</td>\n",
       "      <td>5.904900e-04</td>\n",
       "      <td>3.306744e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>28.0</td>\n",
       "      <td>5.314410e-03</td>\n",
       "      <td>1.488035e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>7</td>\n",
       "      <td>8.0</td>\n",
       "      <td>4.782969e-02</td>\n",
       "      <td>3.826375e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>8</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.304672e-01</td>\n",
       "      <td>4.304672e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Number of Success (>=75 bps)  Number of Possible Ways to Reach x  \\\n",
       "0                             0                                 1.0   \n",
       "1                             1                                 8.0   \n",
       "2                             2                                28.0   \n",
       "3                             3                                56.0   \n",
       "4                             4                                70.0   \n",
       "5                             5                                56.0   \n",
       "6                             6                                28.0   \n",
       "7                             7                                 8.0   \n",
       "8                             8                                 1.0   \n",
       "\n",
       "   Probability for Each Way  Probability for x, p(x)  \n",
       "0              1.000000e-08             1.000000e-08  \n",
       "1              9.000000e-08             7.200000e-07  \n",
       "2              8.100000e-07             2.268000e-05  \n",
       "3              7.290000e-06             4.082400e-04  \n",
       "4              6.561000e-05             4.592700e-03  \n",
       "5              5.904900e-04             3.306744e-02  \n",
       "6              5.314410e-03             1.488035e-01  \n",
       "7              4.782969e-02             3.826375e-01  \n",
       "8              4.304672e-01             4.304672e-01  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.18689526999999992\n",
      "There is a moderate 18.7 percent probability that the manager would show the record he did (or a worse record) if he had the skill to meet your expectations 90 percent of the time.\n"
     ]
    }
   ],
   "source": [
    "# Solving for P(x>0.75), i.e. number of success of 6\n",
    "cdf = df['Probability for x, p(x)'].to_list()\n",
    "print(sum(cdf[0:7]))\n",
    "print('There is a moderate 18.7 percent probability that the manager would show the record he did (or a worse record) if he had the skill to meet your expectations 90 percent of the time.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: You can use other evaluation concepts such as tracking error or tracking risk, defined as the standard deviation of return differences between a portfolio and its benchmark, to assess the manager’s performance.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuous Random Variables\n",
    "Use cases: \n",
    "- modern portfolio theory \n",
    "- risk management technologies\n",
    "- an approximate model for returns\n",
    "\n",
    "### Continuous Uniform Distribution\n",
    "For a continuous uniform random variable between a and b,\n",
    "\\begin{equation}\n",
    "  f(x) =\n",
    "   \\begin{cases}\n",
    "      \\frac{1}{b−a} & \\text{for a<x<b}\\\\\n",
    "      0 & \\text{otherwise}\n",
    "   \\end{cases}       \n",
    "\\end{equation}\n",
    "\n",
    "and \n",
    "\n",
    "\\begin{equation}\n",
    "  f(x) =\n",
    "   \\begin{cases}\n",
    "      0 & \\text{for $x\\le a$}\\\\\n",
    "      \\frac{x-a}{b-a} & \\text{for a<x<b}\\\\\n",
    "      1 & \\text{for $x \\geqslant b$}\n",
    "   \\end{cases}       \n",
    "\\end{equation}\n",
    "\n",
    "For a continuous distribution, the probability of any fixed value is zero because there is an infinite number of outcomes.\n",
    "### The Normal Distribution \n",
    "The normal distribution is popular for quantitative work. This is partly because of the `central limit theorem` – the sum of a large number of independent random variables is approximately normally distributed.\n",
    "\n",
    "**Characteristics**: \n",
    "- The normal distribution can be completely described by two parameters – mean ($μ$) and variance ($σ^2$). Defined as $X∼N(μ,σ^2)$ *(read “X follows a normal distribution with mean $μ$ and variance $σ^2$”)*. \n",
    "- It is symmetric (skewness = 0) and has a **kurtosis of 3**. \n",
    "- The mean, median, and mode are all the same. \n",
    "- Also, a linear combination of two or more normal random variables is (also) normally distributed.\n",
    "\n",
    "#### Multivariate Distribution\n",
    "A multivariate normal distribution for the returns on $n$ stocks is completely defined by three lists of parameters:\n",
    "- the list of the mean returns on the individual securities (n means in total);\n",
    "- the list of the securities’ variances of return (n variances in total); and\n",
    "- the list of all the distinct pairwise return correlations: n(n − 1)/2 distinct correlations in total.\n",
    "\n",
    "**In order to specify the normal distribution for portfolio return, we need the means, variances, and the distinct pairwise correlations of the component securities.**\n",
    "\n",
    "$$f(x)=\\frac{1}{σ\\sqrt{2π}}exp\\big(\\frac{−(x−μ)^2}{2σ^2}\\big)\\ \\ \\  for−∞<x<+∞$$\n",
    "<div style=\"text-align: right\"> Equation 3 </div>\n",
    "\n",
    "The normal density with $μ = 0$ and $σ = 1$ is called the `standard normal distribution` (or `unit normal distribution`). \n",
    "##### Normal Distribution as an Approximate Model for Returns\n",
    "The normal distribution is used to model many security returns. However, the fit is not good for distributions with fat tails (i.e., high kurtosis) or asymmetry. \n",
    "- The normal distribution is a closer fit for quarterly and yearly holding period returns on a diversified equity portfolio than it is for daily or weekly returns.\n",
    "- Normal distribution tends to underestimate the probability of extreme returns. (Since actual return distribution is more fat-tailed).\n",
    "- Option returns are skewed. Cautious in using the normal distribution to model the returns on portfolios containing **significant positions in options**.\n",
    "- The normal distribution is not appropriate for asset prices with a `floor` of zero, even though it is a reasonable fit for asset returns.\n",
    "\n",
    "Having established that the normal distribution is the appropriate model for a variable of interest, we can use it to make the following probability statements:\n",
    "- Approx. 50\\% of all observations fall in the interval μ ± (2/3)σ.\n",
    "- Approx. 68% of all observations fall in the interval μ ± σ.\n",
    "- Approx. 95% of all observations fall in the interval μ ± 2σ.\n",
    "- Approx. 99% of all observations fall in the interval μ ± 3σ.\n",
    "\n",
    "In most cases, the population mean and standard deviation must be estimated from samples. This creates an unlimited number of normal distributions.\n",
    "#### Standardizing a Random Variable\n",
    "It is often useful to standardize the normal random variable $X$ by subtracting the mean and dividing by the standard deviation. This results in a `standard normal random variable` $Z$ that is normally distributed with a mean of 0 and a standard deviation of 1.\n",
    "\n",
    "If we have $X$ ~ $N(μ,σ^2)$ (read “$X$ follows the normal distribution with parameters $μ$ and $σ^2$”), we standardize it using the formula\n",
    "\n",
    "$$Z = \\frac{(X – μ)}{σ}$$\n",
    "<div style=\"text-align: right\"> Equation 4 </div>\n",
    "\n",
    "Tables can be used to determine `cumulative normal distribution` values for $Z$ (i.e., $P(Z≤x)=N(x)$) (`z-score`), N(x) is a conventional notation for the cdf of a standard normal variable (from the z-score table):\n",
    "\n",
    "$$N(−x)=P(Z≥x)=1−N(x)$$\n",
    "$$N(x)=P(Z≤x)=P(Z≥−x)$$\n",
    "\n",
    "Some of the most frequently used values are:\n",
    "- The 90th percentile point is 1.282\n",
    "- The 95th percentile point is 1.650\n",
    "- The 99th percentile point is 2.327\n",
    "\n",
    "The probability that we will observe a value as small as or smaller than 9.5 for $X ~ N(5,1.5)$ is exactly the same as the probability that we will observe a value as small as or smaller than 3 for $Z ~ N(0,1)$. We can answer all probability questions about X using standardized values and probability tables for Z.\n",
    "\n",
    "\\[Revisit Example 8\\]\n",
    "### Applications of the Normal Distribution\n",
    "`Mean–variance analysis` generally considers risk symmetrically in the sense that standard deviation captures variability both above and below the mean. An alternative approach evaluates only `downside risk`. \n",
    "#### Shortfall Risk and Selecting an Optimal Portfolio Using Roy’s Safety-first Criterion\n",
    "`Safety-first rules` focus on `shortfall risk`, the risk that portfolio value will fall below some minimum acceptable level (over some time horizon).\n",
    "\n",
    "$$SFRatio = \\frac{[E(R_P) – R_L]}{σ_P}$$\n",
    "\n",
    "The quantity $E(R_P) − R_L$ is the distance from the mean return to the shortfall level. Dividing this distance by $σ_P$ gives the distance in units of standard deviation.\n",
    "\n",
    "For a portfolio with a given safety-first ratio, the probability that its return will be less than $R_L$ is $N(−SFRatio)$, and the `safety-first optimal portfolio` has the lowest such probability:\n",
    "\n",
    "$$P(R_P<R_L)=N(−SFRatio)$$\n",
    "\n",
    "The SFRatio is just the `Sharpe ratio` with $R_L$ substituted for $R_F$.\n",
    "\n",
    "[Table of the Standard Normal Cumulative Distribution Function](https://courses.cs.washington.edu/courses/cse312/18sp/sections/08/phi_table.pdf)\n",
    "\n",
    "Risk at financial institutions is often measured with `stress testing` and `Value at Risk (VaR)`:\n",
    "- Stress testing uses very unfavorable events or scenarios (`scenario analysis`) to manage the risk. \n",
    "- VaR measures the `minimum` losses expected over a specified time period at a given probability.\n",
    "\n",
    "### Lognormal Distribution and Continuous Compounding\n",
    "The `lognormal distribution` is also widely used in financial modeling, including by the `Black-Scholes-Merton model`. If a random variable $Y$ follows a lognormal distribution, then $ln(Y)$ is normally distributed. The lognormal distribution is bounded below by $0$ and is skewed to the right. This makes it a good choice for asset prices.\n",
    "- Mean ($μ_L$) of a lognormal random variable = \n",
    "$$exp(μ + 0.50σ^2)=e^{μ+0.5σ^2}$$\n",
    "\n",
    "- Variance ($σ_L^2$) of a lognormal random variable = \n",
    "$$exp(2μ + σ^2) × [exp(σ^2) − 1]=e^{2μ+σ^2}\\Big(e^{σ^2}−1\\Big)$$\n",
    "\n",
    "where $μ$ and $σ^2$ are the mean and variance of the associated normal distribution.\n",
    "#### Relationship between the Distribution of Stock Return and Stock Price \n",
    "(theoretical foundation for using the lognormal distribution to model prices)\n",
    "- if a stock’s continuously compounded return is normally distributed, then the future stock price is necessarily lognormally distributed. \n",
    "- stock price may be well described by the lognormal distribution even when continuously compounded returns do not follow a normal distribution. \n",
    "\n",
    "#### Price Relative\n",
    "The `price relative`, is an ending price over a beginning price; it is equal to 1 plus the holding period return on the stock from $t = 0$ to $t = 1$:\n",
    "\n",
    "$$\\frac{S_1}{S_0} = 1 + R_{0,1}$$\n",
    "\n",
    "$$\\frac{S_{t+1}}{S_t} = 1 + R_{t,t+1}$$ \n",
    "\n",
    "where $R_{t,t+1}$ is the rate of return from $t$ to $t + 1$.\n",
    "#### Continuously Compounded Return\n",
    "The continuously compounded return associated with a holding period is the natural logarithm of 1 plus that holding period return, or equivalently, the natural logarithm of the ending price over the beginning price (the price relative).\n",
    "\n",
    "The `continuously compounded return` from $t$ to $t + 1$ is\n",
    "\n",
    "$$r_{t,t+1} = ln\\Big(\\frac{S_{t+1}}{S_t}\\Big) = ln(1 + R_{t,t+1})$$\n",
    "<div style=\"text-align: right\"> Equation 5 </div>\n",
    "\n",
    "The continuously compounded return is smaller than the associated `holding period return`. \n",
    "\n",
    "Continuously compounded return to time $T$ is the sum of the one-period continuously compounded returns:\n",
    "\n",
    "$$r_{0,T} = r_{T–1,T} + r_{T–2,T–1} + …r_{0,1}$$\n",
    "<div style=\"text-align: right\"> Equation 6 </div>\n",
    "\n",
    "Then, the continuously compounded return to the investment horizon $T$ is \n",
    "\n",
    "$$r_{0,T}=ln\\Big(\\frac{S_T}{S_0}\\Big)$$\n",
    "\n",
    "This implies \n",
    "\n",
    "$$S_T=S_0e^{r_{0,T}}$$\n",
    "\n",
    "#### Relation between the Mean and Variance of Continuously Compounded Returns Associated with Different Time Horizons\n",
    "A key assumption in many investment applications is that returns are **`independently and identically distributed (IID)`**. Independence captures the proposition that investors cannot predict future returns using past returns (i.e., weak-form market efficiency). Identical distribution captures the assumption of stationarity.\n",
    "\n",
    "Assume that the one-period continuously compounded returns (such as $r_{0,1}$) are IID random variables with mean $μ$ and variance $σ^2$ (but making no normality or other distributional assumption). Then\n",
    "\n",
    "(Assuming the period returns are independently and identically distributed yields the following:)\n",
    "\n",
    "$$E(r_{0,T}) = E(r_{T–1,T}) + E(r_{T–2,T–1}) + … + E(r_{0,1}) = μT$$\n",
    "<div style=\"text-align: right\"> Equation 7 </div>\n",
    "\n",
    "(we add up $μ$ for a total of $T$ times) and\n",
    "\n",
    "$$σ^2(r_{0,T}) = σ^2T$$  \n",
    "<div style=\"text-align: right\"> Equation 8 </div>\n",
    "\n",
    "This is true even if the individual time period returns are not normally distributed.\n",
    "\n",
    "we can model future stock price $S_{T}$ as a lognormal random variable because $r_{0,T}$ should be at least approximately normal. \n",
    "\n",
    "$$S_T = S_0exp(r_{0,T}) $$\n",
    "#### Volatility\n",
    "Volatility is the standard deviation of the continuously compounded returns on the underlying asset.\n",
    "\n",
    "In practice, we very often estimate volatility using a historical series of continuously compounded daily returns. We gather a set of daily holding period returns and then use Equation 5 to convert them into continuously compounded daily returns. We then compute the standard deviation of the continuously compounded daily returns and annualize that number using Equation 8."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**\n",
    "Suppose an investor views any return below a level of $R_L$ as unacceptable. Roy’s safety-first criterion states that the optimal portfolio minimizes the probability that portfolio return, $R_P$, falls below the threshold level, $R_L$. In symbols, the investor’s objective is to choose a portfolio that minimizes $P(R_P < R_L)$.\n",
    "\n",
    "When portfolio returns are normally distributed, we can calculate $P(R_P < R_L)$ using the number of standard deviations that $R_L$ lies below the expected portfolio return, $E(R_P)$. The portfolio for which $E(R_P) − R_L$ is largest relative to standard deviation minimizes $P(R_P < R_L)$. Therefore, if returns are normally distributed, the safety-first optimal portfolio maximizes the safety-first ratio (SFRatio)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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>Mean and Std (%)</th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Expected annual return</td>\n",
       "      <td>25</td>\n",
       "      <td>11</td>\n",
       "      <td>14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Standard deviation of return</td>\n",
       "      <td>27</td>\n",
       "      <td>8</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Mean and Std (%)   A   B   C\n",
       "0        Expected annual return  25  11  14\n",
       "1  Standard deviation of return  27   8  20"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {'Mean and Std (%)':['Expected annual return', 'Standard deviation of return'], \n",
    "       'A':[25, 27], \n",
    "       'B':[11, 8], \n",
    "       'C':[14, 20]}\n",
    "df = pd.DataFrame(data=data)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Known variables\n",
    "principal = 800000 \n",
    "min_goal = 30000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.749995312505859"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solving for the R_L shortfall level \n",
    "rl = 100* (min_goal / (principal + 1)) # min_goal = rl * (principal + 1)\n",
    "rl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7870372106479311 0.9062505859367675 0.512500234374707\n"
     ]
    }
   ],
   "source": [
    "# Comparing the SFRatio of the three allocations \n",
    "sf_a = (df['A'][0] - rl) / df['A'][1]\n",
    "sf_b = (df['B'][0] - rl) / df['B'][1]\n",
    "sf_c = (df['C'][0] - rl) / df['C'][1]\n",
    "print(sf_a, sf_b, sf_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1824016168069592"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Allocation B's probability of lower than R_L returns: N(−SFRatio)\n",
    "#N(−0.91) = 1 − N(0.91) = 1 − 0.8186 = 0.1814 \n",
    "1 - norm.cdf(0.9062505859367675)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example: estimation of volatility**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>Day</th>\n",
       "      <th>Closing Price (IDR)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Monday</td>\n",
       "      <td>6.95</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Tuesday</td>\n",
       "      <td>7.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Wednesday</td>\n",
       "      <td>6.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Thursday</td>\n",
       "      <td>6.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Friday</td>\n",
       "      <td>6.35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Day  Closing Price (IDR)\n",
       "0     Monday                 6.95\n",
       "1    Tuesday                 7.00\n",
       "2  Wednesday                 6.85\n",
       "3   Thursday                 6.60\n",
       "4     Friday                 6.35"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {'Day':['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday'], \n",
    "       'Closing Price (IDR)':[6.950, 7.000, 6.850, 6.600, 6.350]}\n",
    "df = pd.DataFrame(data=data)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate the volatility of Astra shares. (Annualize volatility based on 250 days in a year.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify the probability distribution for Astra share prices if continuously compounded daily returns follow the normal distribution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Simulation\n",
    "Monte Carlo simulation uses computers to generate many random samples from a specified probability distribution. (An approach to estimating a probability distribution of outcomes to examine **what might happen if particular risks are faced**.) \n",
    "- it provides only statistical estimates, not exact results. \n",
    "- generates a large number of random samples from a specified probability distribution or distributions to represent the role of risk in the system.\n",
    "- requires the identification of risk factors associated with the problem and specify probability distributions for them. \n",
    "\n",
    "Steps of Monte Carlo simulation:\n",
    "1. Specify the quantities of interest in terms of underlying variables.\n",
    "2. Split the time horizon into subperiods.\n",
    "3. Specify the distributional assumptions for the underlying risk factors.\n",
    "4. Use a computer program to draw K random values of each risk factor.\n",
    "5. Calculate the underlying variables based on the random values drawn.\n",
    "6. Compute the quantities of interest.\n",
    "7. Repeat steps 4 to 6 for N trials. The Monte Carlo estimate is the mean quantity of interest over the N trials.\n",
    "\n",
    "### Application Examples\n",
    "- valuing complex securities (e.g. securities with embedded options)\n",
    "- planning\n",
    "- develop estimates of VaR\n",
    "\n",
    "### Comparison/Alternatives\n",
    "A drawback of this approach is that any risk not represented in the time period selected (for example, a stock market crash) will not be reflected in the simulation.\n",
    "\n",
    "`Historical simulation` (or `back simulation`) is an alternative to Monte Carlo simulation. It assumes the past applies to the future. However, unlike the Monte Carlo method, historical simulation does not let itself to \"what if\" analyses because it only reflects the tendencies that appear in the data.\n",
    "\n",
    "`Analytical methods`, such as the Black-Scholes-Merton model provide more `cause-and-effect analysis` than the Monte Carlo method, which provides probability-based estimates rather than exact results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling and Estimation\n",
    "## Sampling\n",
    "    During the exam, you are expected to know the common values for a standard normal distribution (i.e., 85%, 90%, 95%, 97.5%, and 99%). However, if a question requires information from any other distribution table, the relevant portion of the table will be provided with the question.\n",
    "### Simple Random Sampling \n",
    "A sampling plan is needed to set the rules used to select a sample. In a simple random sample, each population element has an equal probability of being selected. It is also often called just a random sample.\n",
    "\n",
    "Simple random sampling requires randomness. This could be done by assigning each member of the population a random number and using a computer program or table of random digits to choose the members.\n",
    "\n",
    "Another approach is `systematic sampling`. This method chooses every $k_{th}$ member until the desired sample size is reached.\n",
    "\n",
    "- The `sampling error` is the difference between the `sample statistic` and the `population parameter` (e.g., the sample mean and the population mean). \n",
    "- The `sampling distribution` of a statistic is the distribution of all statistic values calculated from the same sample size from the same population.\n",
    "\n",
    "### Stratified Random Sampling\n",
    "In a `stratified random sampling`, the population is first divided into subgroups (strata) based on some criteria. Simple random samples are drawn from each subgroup in proportion to the subgroup's relative size to the entire population. This method results in less variance than estimates derived from simple random sampling.\n",
    "\n",
    "Stratified sampling is commonly used to create portfolios that are meant to track a `bond index`. \n",
    "1. First, the entire population of bonds in the index is divided into subgroups based on factors such as maturity, sector, credit quality, etc. \n",
    "2. The manager then selects a sampling of bonds from within each subgroup.\n",
    "\n",
    "### Time-Series and Cross-Sectional Data\n",
    "Assumptions: the sampled returns come from the same population. It is important that the data come from the same population. For time-series data, this means the time period should not be too long because fundamental changes could have occurred (e.g., economic recession). \n",
    "\n",
    "Cross-sectional data are data on certain characteristics of individuals or groups at a single point in time (e.g., year-end P/E ratios for all S&P 500 companies). Assumptions: we must be sure that all the data do come from the same underlying population. (for example, analysts often summarize company-level data by industry). \n",
    "### Distribution of the Sample Mean\n",
    "The `sample mean` is a random variable with a probability distribution called the statistic's sampling distribution.\n",
    "\n",
    "**The `Central Limit Theorem`**: Given a population described by any probability distribution having mean $μ$ and **finite** variance $σ^2$, the sampling distribution of the sample mean $\\bar{X}$ computed from samples of size $n$ from this population will be approximately normal with mean $μ$ (the population mean) and variance $\\frac{σ^2}{n}$ (the population variance divided by $n$) when the sample size $n$ is large (usually greater than 30).\n",
    "\n",
    "**The `Standard Error of the Sample Mean (SEM)`**: For sample mean $\\bar{X}$ calculated from a sample generated by a population with standard deviation $σ$, the standard error of the sample mean is given by one of two expressions:\n",
    "\n",
    "$$σ\\bar{X}=\\frac{σ}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 1 </div>\n",
    "\n",
    "- when we know $σ$, the population standard deviation, or by\n",
    "\n",
    "$$s\\bar{X}=\\frac{s}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 2 </div>\n",
    "\n",
    "- when we do not know the population standard deviation and need to use the sample standard deviation, $s$, to estimate it.\n",
    "\n",
    "In practice, we almost always need to use Equation 2. The estimate of $s$ is given by the square root of the sample variance, $s^2$, calculated as $s^2=\\frac{\\sum_{i=1}^n{(Xi−\\bar{X})^2}}{n−1}$\n",
    "<div style=\"text-align: right\"> Equation 3 </div>\n",
    "\n",
    "To summarize, according to the central limit theorem, when we sample from any distribution, the distribution of the sample mean will have the following properties as long as our sample size is large:\n",
    "- The distribution of the sample mean $\\bar{X}$ will be approximately normal.\n",
    "- The mean of the distribution of $\\bar{X}$ will be equal to the mean of the population $\\mu$ from which the samples are drawn.\n",
    "- The variance of the distribution of $\\bar{X}$ will be equal to the variance $\\sigma^2$ of the population divided by the sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point and Interval Estimates of the Population Mean\n",
    "The two branches of `statistical inference` are `estimation` and `hypothesis testing`. \n",
    "- Estimation seeks to find the value of the parameter.\n",
    "- Hypothesis testing seeks to find if the value of a parameter equals some specific value. \n",
    "\n",
    "### Point Estimators\n",
    "Estimators are the formulas used to calculate the sample statistics. Estimates are the particular values derived from these estimators.\n",
    "- **Unbiased** estimators' expected value (the mean of its sampling distribution) equals the parameter it is intended to estimate.\n",
    "- **Efficient unbiased** estimators have the smallest sampling distribution variance for a given sample size. Repeated samples will generate results that are tightly grouped around the parameter.\n",
    "- **Consistent** estimators get closer to the population parameter's value as the sample size increases.\n",
    "\n",
    "### Confidence Intervals\n",
    "A point estimate is unlikely to exactly equal the population parameter due to sampling error. An `interval estimate` is often more appropriate.\n",
    "\n",
    "A `confidence interval` is a range for which one can assert with a given probability $1 − α$, (called the `degree of confidence`), where α is the `significance level`, that it will contain the parameter it is intended to estimate. This interval is often referred to as the `$100(1 − α)\\%$ confidence interval for the parameter`. \n",
    "\n",
    "    For example, using a 5% significance level creates a 95% confidence interval around the sample mean. We can be 95% confident that the population mean falls somewhere in this interval.\n",
    "#### Construction of Confidence Intervals\n",
    "A $100(1−α)\\%$ confidence interval is calculated by:\n",
    "\n",
    "$$Point\\ estimate ± Reliability\\ factor \\times Standard\\ error$$\n",
    "\n",
    "where\n",
    "\n",
    "`Point estimate` = a point estimate of the parameter (a value of a sample statistic)\n",
    "\n",
    "`Reliability factor` = a number based on the assumed distribution of the point estimate and the degree of confidence (1 − α) for the confidence interval\n",
    "\n",
    "`Standard error` = the standard error of the sample statistic providing the point estimate\n",
    "\n",
    "$z_α$ denotes the point of the standard normal distribution such that $α$ of the probability remains in the right tail.\n",
    "\n",
    "*Note: The reliability factor in this case is based on the standard normal distribution, which has a mean of 0 and a variance of 1*\n",
    "##### Confidence Interval with Known Population Variance\n",
    "A $100(1 − α)\\%$ confidence interval for population mean $μ$ when we are sampling from a normal distribution with known variance $σ^2$ is given by\n",
    "\n",
    "$$\\bar{X}±z_{\\frac{α}{2}}\\times \\frac{σ}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 4</div>\n",
    "\n",
    "This produces a confidence interval with upper and lower bounds with a total of $α$ probability that the population mean is outside the confidence interval. $z_{\\frac{α}{2}}$ is used because $\\frac{α}{2}$ represents what percent would be in each tail.\n",
    "\n",
    "    For example, a 95% confidence interval uses the 1.96 reliability factor for z_{0.025} because there is a 2.5% probability that the population mean is greater than the upper bound and a 2.5% probability that it is less than the lower bound. Other commonly used reliability factors include z_{0.05}=1.65 for a 90% confidence interval and z_{0.005}=2.58 for a 99% confidence interval.\n",
    "##### Confidence Interval with Unknown Population Variance  (the z-Alternative)\n",
    "A $100(1 − α)\\%$ confidence interval for population mean $μ$ when sampling from any distribution with unknown variance and when sample size is large is given by\n",
    "\n",
    "$$\\bar{X}±z_{\\frac{α}{2}}\\times \\frac{s}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 5</div>\n",
    "\n",
    "Reliability Factors for Confidence Intervals Based on the Standard Normal Distribution. We use the following reliability factors when we construct confidence intervals based on the standard normal distribution:\n",
    "- 90 percent confidence intervals: Use $z_{0.05}$ = 1.65\n",
    "- 95 percent confidence intervals: Use $z_{0.025}$ = 1.96\n",
    "- 99 percent confidence intervals: Use $z_{0.005}$ = 2.58\n",
    "\n",
    "##### Confidence Interval with the Population Mean (Population Variance Unknown) (the t-Distribution)\n",
    "It is more appropriate and conservative to use the `t-distribution` for confidence intervals when the population variance is unknown. This is valid even when the **sample size is small**. Since it is more conservative (i.e., the reliability factor is bigger), the confidence interval will be wider.\n",
    "\n",
    "The t-distribution is based on/defined by the `degrees of freedom (df)`. Degrees of freedom relates to the number of independent observations. Since the sample mean must be determined before calculating the sample variance, there are really only $n−1$ independent observations. **If you know $n−1$ observations and the sample mean, then the $n_{th}$ observation can be determined. That is why you must divide by $n−1$ rather than $n$ when calculating the sample variance.** \n",
    "\n",
    "**As the degrees of freedom increase, the t-distribution approaches the standard normal distribution.**\n",
    "\n",
    "If we are sampling from a population with unknown variance and either of the conditions below holds:\n",
    "- the sample is large, or\n",
    "- the sample is small, but the population is normally distributed, or approximately normally distributed,\n",
    "\n",
    "then a $100(1 − α)\\%$ confidence interval for the population mean $μ$ can be given by\n",
    "\n",
    "$$\\bar{X}±t_{\\frac{α}{2}}\\times \\frac{s}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 6 </div>\n",
    "\n",
    "where the degrees of freedom for $t_{\\frac{α}{2}}$ is $n − 1$, $n$ is the sample size.\n",
    "### Selection of Sample Size\n",
    "The reliability factor depends on the chosen statistic ($t$ or $z$) and the chosen degree of confidence. The sample size also affects the width of the confidence interval. Larger sample sizes decrease the width because the sample standard deviation is divided by the square root of $n$.\n",
    "\n",
    "While increasing the sample size improves precision, it is not always recommended. First, it increases the chance of sampling from more than one population. Second, larger samples are more expensive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>Population</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Normally</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Distributed?\\tPopulation</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Known</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Variance?\\tSmall Sample\\tLarge Sample</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Yes\\tYes\\tz\\tz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Yes\\tNo\\tt\\tt or z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>No\\tYes\\tNot Available\\tz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>No\\tNo\\tNot Available\\tt or z</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              Population\n",
       "0                               Normally\n",
       "1               Distributed?\\tPopulation\n",
       "2                                  Known\n",
       "3  Variance?\\tSmall Sample\\tLarge Sample\n",
       "4                         Yes\\tYes\\tz\\tz\n",
       "5                     Yes\\tNo\\tt\\tt or z\n",
       "6              No\\tYes\\tNot Available\\tz\n",
       "7          No\\tNo\\tNot Available\\tt or z"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The following table summarizes which reliability factor to use for various scenarios.\n",
    "df = pd.read_csv('reliability_factors.csv', delimiter='\\t')\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The interval 0.401 and 0.499 includes the population mean is with 90% confidence.\n"
     ]
    }
   ],
   "source": [
    "# Calculate and interpret the 90 percent confidence interval for the population mean of all US equity mutual funds\n",
    "sample_size = 100\n",
    "s = 0.30\n",
    "mean = 0.45 # avg. Sharpe ratio\n",
    "\n",
    "# Step 1. Finding the reliability factor for 90% confidence interval. (Note that Z0.05 = 1.65)\n",
    "z_005 = 1.65\n",
    "\n",
    "# Step 2. Implementing equation 5\n",
    "lower = mean - z_005 * (s / sample_size**0.5)\n",
    "higher = mean + z_005 * (s / sample_size**0.5)\n",
    "\n",
    "#print('The confidence that the population mean falls between interval', str(round(lower, 2)), 'and', str(round(higher, 2)), 'is 90%.')\n",
    "print('The interval',str(round(lower, 3)), 'and', str(round(higher, 3)), 'includes the population mean is with 90% confidence.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example 6\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More on Sampling\n",
    "### Data-Mining Bias \n",
    "Ronald Coase, \"If you torture the data long enough, it will confess.\"\n",
    "\n",
    "Data mining relates to overuse of the same or related data.\n",
    "- **`Out-of-Sample Test`** is a test of a strategy or model using a sample outside the time period on which the strategy or model was developed. If a relationship is genuine, it should be found in other datasets. Ultimately, a model is only successful if it works in the future.\n",
    "- **Intergenerational data mining** involves using information developed by previous researchers using a dataset to guide current research using the same or a related dataset.\n",
    "\n",
    "### Sample Selection Bias\n",
    "When data availability leads to certain assets being excluded from the analysis. \n",
    "- `survivorship bias` sometimes appears when we use both stock price and accounting data. For example, many studies have shown the stocks of companies with low price-to-book ratios tend to outperform in future periods. This could be because companies that fail are excluded from the studies.\n",
    "- A sample can also be biased because of the removal (or `delisting`) of a company’s stock from an exchange.\n",
    "\n",
    "Newer asset classes such as hedge funds may present even greater problems of sample selection bias. \n",
    "### Look-Ahead Bias\n",
    "`Look-ahead bias` occurs if the information is used that would not have been available on the test data. For example, accounting information such as book value will not be available for some time after the end of the period.\n",
    "### Time-Period Bias\n",
    "A test design is subject to time-period bias if it is based on a time period that may make the results time-period specific.\n",
    "\n",
    "Ideally, an analyst should test market anomalies over several business cycles to ensure that results are not period specific. This bias can favor a proposed strategy if the time period chosen was favorable to the strategy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing\n",
    "    During the exam, you are expected to know the common values for a standard normal distribution (i.e., 85%, 90%, 95%, 97.5%, and 99%). However, if a question requires information from any other distribution table, the relevant portion of the table will be provided with the question.\n",
    "    \n",
    "Use cases: \n",
    "- Is the underlying mean return on this mutual fund different from the underlying mean return on its benchmark?\n",
    "- Did the volatility of returns on this stock change after the stock was added to a stock market index?\n",
    "- Are a security’s bid-ask spreads related to the number of dealers making a market in the security?\n",
    "- Do data from a national bond market support a prediction of an economic theory about the term structure of interest rates (the relationship between yield and maturity)?\n",
    "\n",
    "## Introduction\n",
    "Two hypotheses are always stated: \n",
    "- the `null hypothesis`, $H_0$\n",
    "- the `alternative hypothesis`, $H_a$\n",
    "\n",
    "- The null hypothesis is assumed true until tests prove otherwise. \n",
    "- The alternative hypothesis is only accepted if there is sufficient evidence to reject the null hypothesis.\n",
    "\n",
    "### Sequence of Hypothesis Testing\n",
    "1. State the hypothesis.\n",
    "- Identify the appropriate `test statistic` and its probability distribution.\n",
    "- Specify the significance level.\n",
    "- State the decision rule.\n",
    "- Collect data and calculate the test statistic.\n",
    "- Make a `statistical decision`.\n",
    "- Make an `economic decision`.\n",
    "\n",
    "### Stating the Hypotheses\n",
    "The choice of null and alternative hypotheses should be based on the `hoped-for condition`. For example, if an analyst is attempting to show that the mean annual return on European equities has exceeded 8%:\n",
    "- the null hypothesis ($H_0$) should be that the mean return is less than 8%. \n",
    "- the alternative hypothesis ($H_a$) should only be accepted if statistical tests provide sufficient evidence that the mean return is greater than 8%.\n",
    "\n",
    "#### Two-Sided/Two-Tailed Hypothesis Test (Equal/Not Equal to) \n",
    "The null hypothesis is that the population parameter is \"equal to\" a certain value. The alternative hypothesis is that it is \"not equal to\" this value.\n",
    "\n",
    "$$H_0: θ = θ_0 \\ versus\\ H_a: θ ≠ θ_0$$\n",
    "#### One-Sided/One-Tail Hypothesis Test (Greater/Smaller Than)\n",
    "$$H_0: θ ≤ θ_0 \\ versus \\ H_a: θ > θ_0$$\n",
    "\n",
    "$$H_0: θ ≥ θ_0 \\ versus \\ H_a: θ < θ_0$$\n",
    "### Test Statistic and Its Probability Distribution\n",
    "A `test statistic` is a quantity, calculated based on a sample, whose value is the basis for deciding whether or not to reject the null hypothesis.\n",
    "\n",
    "The distribution for the test statistic is assumed to be one of the following:\n",
    "- `t-distribution` (for a t-test)\n",
    "- Standard normal or `z-distribution` (for a z-test)\n",
    "- `Chi-square (χ2)` distribution\n",
    "- `F-distribution` (for an F-test)\n",
    "\n",
    "$$Test\\ statistic=\\frac{Sample\\ statistic−Value\\ of\\ the\\ population\\ parameter\\ under\\ H_0}{Standard\\ error\\ of\\ the\\ sample\\ statistic}$$\n",
    "<div style=\"text-align: right\"> Equation 1 </div>\n",
    "E.g. the test statistic will be calculated by:\n",
    "$\\frac{\\bar{X_{RP}}−\\mu_{0}}{s\\bar{X}}=\\frac{\\bar{X_{RP}}−0}{s/\\sqrt{n}} = \\frac{\\bar{X_{RP}}}{s\\bar{X}}$\n",
    "\n",
    "For a sample mean, $\\bar{X}$, calculated from a sample generated by a population with standard deviation $\\sigma$, the standard error is given by one of two expressions:\n",
    "\n",
    "$\\sigma\\bar{X}=\\frac{\\sigma}{\\sqrt{n}}$\n",
    "<div style=\"text-align: right\"> Equation 2 </div>\n",
    "\n",
    "(when we know $σ$ (the population standard deviation), or)\n",
    "\n",
    "$s\\bar{X}=\\frac{s}{\\sqrt{n}}$\n",
    "<div style=\"text-align: right\"> Equation 3 </div>\n",
    "\n",
    "(when we do not know the population standard deviation $\\mu$ and need to use the sample standard deviation $s$ to estimate it.)\n",
    "\n",
    "The null hypothesis can be rejected or not rejected after the test statistic has been calculated. The decision is based on a comparison that assumes a specific significance level, which establishes how much evidence is required to reject the null hypothesis. The table below summarizes the four possible outcomes.\n",
    "\n",
    "### Significance Level\n",
    "The `level of significance` reflects how much sample evidence we require to reject the null. It is the probability of incorrectly rejecting the null ($\\alpha$).\n",
    "\n",
    "The **`power of a test`** is the probability of correctly rejecting the null—that is, the probability of rejecting the null when it is false.\n",
    "\n",
    "$$power = 1- β$$.\n",
    "\n",
    "The probability of a Type I error in testing a hypothesis, denoted by $α$. Conventional significance levels to conduct hypothesis tests are: \n",
    "- 0.10 (we have some evidence that the null hypothesis is false)\n",
    "- 0.05 (we have strong evidence that the null hypothesis is false)\n",
    "- 0.01 (we have very strong statistical evidence that the null hypothesis is false) <br>\n",
    "\n",
    "The probability of a Type II error is $β$. Reducing the probability of a Type I error by decreasing α will increase the probability of a Type II error.\n",
    "### Stating the Decision Rule (Rejection Point)\n",
    "When we test the null hypothesis, if we find that the calculated value of the test statistic is more extreme than a given value or values determined by the specified level of significance, $α$, we reject the null hypothesis. We say the result is statistically significant. Otherwise, we do not reject the null hypothesis and we say the result is not statistically significant. \n",
    "\n",
    "A `rejection point` (`critical value`) for a test statistic is a value with which the computed test statistic is compared to decide whether to reject or not reject the null hypothesis. In a graph, notice the area in the rejection regions correspond to the significance level.\n",
    "\n",
    "1. one-tailed test\n",
    "rejection points are indicated using the symbol for the test statistic with a subscript indicating the specified probability of a Type I error, $α$; \n",
    "  - For the one-sided test $H_0:θ≤θ_0$, the critical value is $z_α$.\n",
    "  - For the one-sided test $H_0:θ≥θ_0$, the critical value is $−z_α$.\n",
    "2. two-tailed test<br>\n",
    "For a two-sided test $H_0:θ=θ_0$, two rejection points exist. If using the z-distribution, the values would correspond to $±z_{\\frac{α}{2}}$. \n",
    "\n",
    "#### Relationship between Confidence Interval and Hypothesis Testing\n",
    "A two-sided test at an $α$ significance level is the same as a $(1−α)$ confidence interval. If the test statistic falls outside the confidence interval, then reject the null hypothesis.\n",
    "- If the hypothesized population mean, $μ_0$, is less than the lower limit of the $95\\%$ confidence interval based on the sample mean, we must reject the null hypothesis at the $5\\%$ significance level (the test statistic falls in the rejection region to the right).\n",
    "- If the hypothesized population mean is larger than the upper limit of the 95% confidence interval, we reject the null hypothesis at the $5\\%$ level (the test statistic falls in the rejection region to the left). \n",
    "\n",
    "### P-Value (Alternative Approach to Hypothesis Testing)\n",
    "The `p-value` is the smallest level of significance at which the null hypothesis can be rejected. Analysts and researchers often report the p-value (also called the marginal significance level) associated with hypothesis tests.\n",
    "\n",
    "We can use p-values in the hypothesis testing framework presented above as an alternative to using rejection points. **The smaller the p-value, the stronger the evidence in rejecting the null hypothesis.**\n",
    "- If the p-value is less than our specified level of significance, we reject the null hypothesis. \n",
    "- Otherwise, we do not reject the null hypothesis. \n",
    "\n",
    "### Making the Economic or Investment Decision\n",
    "Statistical significance does not necessarily imply economic significance. For example, small differences in portfolio returns can be statistically significant, but transaction costs and other fees could prevent successful implementation.\n",
    "- statistical decision (from step six)\n",
    "- all pertinent economic issues\n",
    "- nonstatistical considerations: such as investors' risk tolerance and financial position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step One Example**: risk premium hypothesis\n",
    "$H_0: \\mu_{RP} \\leq 0\\ \\ population\\ mean\\ risk\\ premium\\ for\\ US\\ equities (\\mu_{RP})\\ is\\ less\\ than\\ or\\ equal\\ to\\ zero$\n",
    "$H_a:\t\\mu_{RP} > 0\\ \\ The\\ population\\ mean\\ risk\\ premium\\ on\\ US\\ stocks\\ is\\ positive.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example** Assume we can conduct a z-test based on the central limit theorem because our US sample has many observations. We can conduct a z-test because we can plausibly assume that the test statistic follows a standard normal distribution. (TBC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example Step 3: specifying the significance level\n",
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwAAAAHOCAYAAAAv/G5jAAAAGXRFWHRTb2Z0d2FyZQBBZG9iZSBJbWFnZVJlYWR5ccllPAAAAxxpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADw/eHBhY2tldCBiZWdpbj0i77u/IiBpZD0iVzVNME1wQ2VoaUh6cmVTek5UY3prYzlkIj8+IDx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IkFkb2JlIFhNUCBDb3JlIDUuNi1jMTQ1IDc5LjE2MzQ5OSwgMjAxOC8wOC8xMy0xNjo0MDoyMiAgICAgICAgIj4gPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4gPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIgeG1sbnM6eG1wPSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvIiB4bWxuczp4bXBNTT0iaHR0cDovL25zLmFkb2JlLmNvbS94YXAvMS4wL21tLyIgeG1sbnM6c3RSZWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20veGFwLzEuMC9zVHlwZS9SZXNvdXJjZVJlZiMiIHhtcDpDcmVhdG9yVG9vbD0iQWRvYmUgSWxsdXN0cmF0b3IgQ1M1LjEiIHhtcE1NOkluc3RhbmNlSUQ9InhtcC5paWQ6MTc1Njc3RjQ1RUMyMTFFOUE0MEJDMzQ2MDExODY0OTkiIHhtcE1NOkRvY3VtZW50SUQ9InhtcC5kaWQ6MTc1Njc3RjU1RUMyMTFFOUE0MEJDMzQ2MDExODY0OTkiPiA8eG1wTU06RGVyaXZlZEZyb20gc3RSZWY6aW5zdGFuY2VJRD0ieG1wLmlpZDoxNzU2NzdGMjVFQzIxMUU5QTQwQkMzNDYwMTE4NjQ5OSIgc3RSZWY6ZG9jdW1lbnRJRD0ieG1wLmRpZDoxNzU2NzdGMzVFQzIxMUU5QTQwQkMzNDYwMTE4NjQ5OSIvPiA8L3JkZjpEZXNjcmlwdGlvbj4gPC9yZGY6UkRGPiA8L3g6eG1wbWV0YT4gPD94cGFja2V0IGVuZD0iciI/Pj0kIGoAAEARSURBVHja7N0/bBxZmiD4UHUZO+gGmlNAAz3bB1QWegaYAWa2soyd23VWKWtnLVHmWEp5d5bEc/bOOVHWrkfJWpMp60xS1mAtUuZaZHm7mBKUBUzNChAwYgGN3h5L9z7yi1IolEkm/yQZmfn7AYEgk5nJzBd/3ve9ePHerQqAhfLbX/9qvaxul6Wfy1rjz0dlOSzLy7Lsvnrz9lCJAdB0SxEALETQ3yurh2UZZsAfgf1+Wb4ty7jx1Hje12UZZHIQf3tWllFJBo6UJAASAIBuB/4R7D8uy6MM5p9nMD8+Z9IQNsrrRkoVAAC6Gfw/Ksu7XIaXSSLKslWW92XZy8QAgBXlCgBA9wL/6LqzXZ104XlalidX0X2n8b6RANwr77mvtAEA4GaD/81GS31/Du8fVwO2838MlTgAANxM4N8vy0F293l0jYnGptIHAIDrDf7rvv4H82j1P+X/DjMJ2LIVAABg/gF4szvO1g19hqHuQACr5WeKAOBGAu9o6d8py78py9++evP22U18jne/+/3hF7/4+Y/lx/9S1t/H77YOAABcbfA/aHT56XXkM9VXIvq2EAAAXF2gXXe52c5Jvrr02fYyMVmzpQAA4PIB9maXb7rNexJex5UJWwsAAC4XXC/E2Ps5HOk7IwMBAMCSB/+Nz1t3U1q39QAAYImD/9bnfteVm5QBuDqfKQKA+QXRZRWB/4NXb96OFuzjb5RlXJZtWxIAAGYI/hd9gi33AwAAwGyB89LMrmumYAAAWLGAuXE/gEnCAACgESj360m+lvC77eUcASYJA1hwbgIGuJoAuVdWe2XZffXm7YMl/Ir3yhLBv5uCARbczxQBwKWD/wiM/64sRxEov/vd7/+wbN8xvtMXv/j5/yg/bpb1t+X3/27LAwCwqgnAyoyZHyMC5XfVFQhgQekCBHC5gPhRdTLW/71Xb96OV+ArP6lOrnToCgQAwMoF//VNv5sr9r0H+b3X7QUAAKxKELyWo+Lsrej31xUIYEHpAgRwMdEFJoLfeyv6/Z/k2izBAAAst8ZkX4MVL4d15QAAwLIHvb3s+qLlu/ppgrADJQGwOHQBAjif6Pozrj50gVl1MelZP0dDAgCA5RFBbnZ56SuNj8pl0w3BAIvDFQCA2YLcXlk9LsuTV2/eHiqRjzytTuYG0C0KAIClSQD0dT+9fIaujgAAsGzB7UBpnJkk7SkJAAAWOahdM+rPzGVlhmCABeAeAIDTxag/0b/dqD9nePXm7X5ZjSr3AgAAsIi0aF+ozHpZZoYFBQBg4YLZ12XZURLnLjfDggJ0mC5AAFOC2LKKAHZDaZzb01y7CgAAwEIE/71swd5UGhdPoFwFAOgmVwAAPhUTfh29evNWAnBxJgcDAKD73Ph7pWVZz5/QUxoAAHQ1aDWZ1dWWZ9xIva0kALpDFyCAD8FqtPoPKmP+X6Uoy6GrAAAAdDEB0FqtXAGWnisAACdBagxZ2au0/s9DDKUaVwEGigIAgC4E/2uG/Zx7Gbu3AqAjXAEA+DBh1VNFMTdxZWXgKgAAADcqJ/16r/X/WsraVQAAAG48KN2Om1SVxLWUdT3HwkBpAABwEwFpPwPSodK4tjJ3FQAAAMHoCpW5qwAAAAhEJV4AADDfIDQmp9pREpIvgFVjGFBgFQPQYXUy6deG0rh+r9683S+rWB4rDQAA5h38r2Xr/7bSuNHtUF8F6CkNgOvlCgCwamLSr7VK6/+NyqsA48pVAAAJAMC8ROt/WT0sy7MSgB4pkRsXswMPXQUAkAAAzEvd2vxUUdy8koSNKlcBAACYh2hlNulXJ7fL0L0AAADMI9DciZt/lUQnt42bsgGukS5AwCoEmIOyWq/c+NtV7gUAAOBKEwAzz3Z/G70ry6aSAJg/VwCAZQ8so+V/UJ20MtNdz8ryMEdqAkACAHBhW2UZ5bjzdFc9MtMjRQEgAQC4kBzxp1dp/e+8nJfBVQAACQDAhYP/CCKj9f9pCS7HSmQh1FcB1hUFgAQA4LzqriRa/xdEXgXYrUwMBiABADiPbP1/WJZnGVSyOCJh65mwDUACAHAe0fUnAv+nimKxZHetUeUqAIAEAGAWOZnUsCxPtP4vrLgZuJcTuAEAwKkJwE5ZXiuJhd+OJm8DmBNXAIBlChoH1ckIMg+UxsKLewEGrgIAAHBaAqDVeLm250FZtpUEwNVyBQBYlmAxWv4HlWE/l0ncCzDM+zoAAOCjBOC11mLbFYCzuQIALEOQOCyrXqX1fxnFVYD1nNsBAAkAIPg/Dgxj3P+nOYY8y2WU60eKAkACANAMDLX+L6GcyyGuAjxUGgASAGDFZet/BIbPTPq11EZlWcuuXgBIAIAVtpXrp4pieWXXrkgCHisNAIAVFUNDluW9VuGV294DpQEAsJoB4XYMEakkVmqbm+gN4AroAgQsYiA4KKthWTaUxkqJG70HJgYDAFi9BEBL8OpuexODAVySKwDAogWAg7KKxbCfqym2+9DEYAAAq5MARAvwjpJY6X3gXVk2lQTAxbgCACxS4Dcsq16l7/+qMzEYgAQAWIHgP7p8xDjwoxwXntUV8z6YGAxAAgAsuUcR9FX6/q+8nPV5VJkYDEACACynbP2PLh/PtP6TohtQz8RgABIAYDk9yvVTRUEoieBhWe1X7gUAAFguMelTWd6X5ZHSoLVvrOe+0VMaAADLE+Rtx9CfSoIp+4eJwQDOSRcgoMvBXb+shpUbf5ku7gVYNzEYgAQAWA5bZdl/9ebtSFEwRb1vDBUFgAQAWGA5ukssWv+ZqjEkqJuBASQAwIKLMd6j9X9fUXCGekjQdUUBIAEAFlDO8Dooy4bS4Cw5N8Ru5SoAAMDCJgBGduG8+8wghwTtKw2A07kCAHQtkIvx/nuVvv+cQ3YVG1euAgAALFTwv1aWd2XZUhpcYP8Z5lUAQ4ICnMIVAKBL6tl+tf5zEXEfwFFjPwJAAgB0VbbaRveNZzm0I5xLY0jQ+0oDQAIAdF90+4kA7qmi4BLqIUGHigIAoKNKsNbLvtuCNq5if9opy56SAJjMFQCgC2LSr/GrN29HioIrEFcBYljQnqIAkAAAHRPjt5fVsDLpF1ekMSToY6UBIAEAuieCtP0StO0qCq5QjCS1bkhQAAkA0CHZ+j+oDPvJ1asTyqGiAADoTgJw4GZN5rh/bZXltZIA+JgrAMBNBWfDsuqX5YHSYE7qIUHXFQUAwM0nAK/Lsq0kmPN+tucqE8DHXAEAbiIoe1RWvUrff+bPkKAAADcc/K+V5V1ZNpUG17TPxdWmLSUBcMIVAOC6Pcr1U0XBNYmrAENDggJIAIBrlgHYwwjIXr15e6REuCajssS+52ZgAAkAcM2iG0YE/lr/uTaZbI4y+QQA4DrETZhleZ/Df8J173/93P8GSgMA4HoCsG2TMnHD++CeoWcBdAECrifwigm/hpVhP7lZzys3AwMAXEsCYDImurIvGoIWWHmuAADzDrgGZRWL1n+6YFSW+4oBkAAAzM/jsuy/evN2X1HQATEnQNyQbkhQQAIAcNVyxJ9BWTaUBl1QEtFxWe1WhgQFJAAAcxGt/6MSdB0qCjokbgYexNC0igKQAABckWz9jwBL3386pSSkcQVgXLkKAEgAAK5U3fo/VhR0UNwLYEhQQAIAcBVymMUIrLT+01Wj3EfdDAwAcMngf81Y6yzIvhqzUx8oCWDVuAIAXLVHuX6qKOi46AbUz5mqASQAAOeV/anjxspnr968PVIidFmOTnVYuRkYkAAAXNhWWSLw1/rPonAzMADARcSY6mV5n8N/wiLtu+5ZAVaKKwDAVYlhP8ev3rwdKQoWTFwFuK8YAAkAwIxyRtVhZdhPFlMkrT1XrwAAZk8AYjjF10oC+zBA97kCAFw2cIohFIdl2VAaLLC4ehVXAUwMBgBwRgKwF4uSYAn25W37MrAKXAEALhMwDcoqFn3/WQbPY3/O/RoAgAkJgNZ/7NMAC8YVAOCigdKg0vrP8ol7WQbuBQAA+DQB2NFSypLu20YEApaaKwDARQKkXllFC+kzpcESqkcEGioKAIBKCyn2cYBF5goAcN7AyLj/rIK4CrBW9vdNRQEArHoCYJQUVmVf3yzLu7KsKQ1gmbgCAJwnIBpURv5hdTwty1FZthQFALCqCcDrGP1HSbBC+/ywLO+z6xvAUnAFAJg1ENosq+gKoe8/K+PVm7ejstqvXAUAJADAigX/vbJ6WJZnJSAaKxFWTHR5GxgWFABYpQRgz5CIrPgxsO2GYGBZuAIAnBX4xIRfg7I8UBqssLrr22NFAQAsc/C/lq2e+j/jeHBDMLAkXAEATrNdnQyDaNhPVl7jhuBtpQEALJ3f/vpXj7K1c6A04Kfjop/HxSOlAQAsU5DT0/UHph4f9QzBPaUBACxLgHMQi5KAU4+RPSUBLCL3AADtwGazrHqVUX/gNHF8DHQFAgAWPfgf6N8MsyfLugIBi+iWIgAymIkJjmKyr/1Xb97eUyIw0zET3YCOyjFzR4kAi0IXIKC2U50M+anrD8ygBP318aIrEACwWGK0HxMcwYWPn7orkOMHAFiI4GU9g/+h0oALH0dGzgIWhi5AsNpBS7RYxqymo5zlFLiY6ArUy1G0ADrtZ4oAVjb4jxsY637/f/vud7//g1KBiynHz5svfvHzfy4//qeyfll+HysVAKBrCcCOIQxhLsfV60ywATpJFyBYzSBls6zWy3Lv1Zu3YyUCVya6AkXwv60oAICuBP9DN/3CXI8xE+oBAJ0JTPrZ7UfrJMz3WNs0tC7QVWYChtUJSOqZfg/NWgrXcszFLMG9snyTk4YBdIJ7AGB1gv8IRsZluadE4FrEseZ+AEACANyICEB6ZXmgJRKuRx5rkQSsux8AALg20d8/+yIPlAbcyDHofgAA4NoCDyP+QDeOxXreDfMDADdOFyBY4uC/Oun68+TVm7cjJQI3KuYHiC5Be4oCAJhH8G+4T3BcAgArFmRoaYTuHZ+65QE37meKAJYquKiH+3xTlv/w7ne//4NSge4ox+ThF7/4eRyn/7msX5Tf3ygV4LqZCAyWL/gPdwz3CZ0+XnfKalCdTBI2ViLAdXITMCxX8N+rjPUPiyBuCo7Af8fIQIAEALiIrQz+o+X/UHF0Plm7ttfRTZmkP8jj1k3BwLXSBQgWP6CM4GFd8N/57RQzwd4ty8uy3C7LxizbK28Wjdd9m+v98rqNKc99nQFlbVSe+0Dpd3q/iMnBDsrydNp2BQD4KPg3oshiBP8ZnNe/D2aZFKr8fT1eVz8v1vn71qREoSwHMfpTY+kp/YXYP4wMBADMFDRsChoWYjutZbD/qPX467PGg8/nbE4JFnutxw37uhzH87rSAObNPQCwmMFCBP2Pq5MbfkdKpNMioIsW/HZ3n8P827RtHF1DehP+VL/PsLU/1MmBewUWUDmOI9GLY3k7tz2ABAD4KNjbFvwvjLu5bo/MFH36104J9tamBIp1AnA794e1TAYHuV+8bl9tYGGSgHpkoD1JADBPn99A8BKVVQQwX5blZTnh7Z7jtYOsTH+sTm6YOpoSHH2dz4mb5fanvFevarSgpZHxmFmQ4H8k+F8Ya63AfaZAPwPBKs9n7XNXM2iM8+BXeX6M5WFZtsrvX7qpdCHdqU6G9I0rAebzgOnx2y/L8nzWwS9miT9b7/191rVHp8Skg9b5eHNRyvGzG9hoMdrBOCumu2f1gW28Ngo1WrmelCU22kG7hSTfK26O6+dz9055/8et5b7gnwUK/o3ssuTyfBQNGOtZ0dTqbkOHrefvZ+XzVf7tkVbkhdzuR5kErGUdpksXfKgH4/y3k3Hgk0yUh1cRfzZG5Lqb59mIJ1+fch7dbsWRXy5SWV53F6AorP1G1hUbYXjWTU9Z+dX9nY8y29utGmMn1ztA+dsfl+VOVoLjfP/+hB1hLU+yzQUE/3TJvUwC6lF9ouK7nX97eUoAea+VLLB4SUBsw6ir3NwN1U8t+FEPPos4MI+TSAK2ZhjxbJb4M54TV91idu6v8r3rLpaT6uTnrRjSFdcpG64/aYSDrNQOznjtTvs5OYzeT+83qc9rDqH3fsIoGtuGx2ORWjxyPzZZ0GJuv3qo1rXW41uTRvOZVvE1hgI9aA4pespr9trnPhay3nzn2IefhlOedC59P2lo5PPEn/mcwYTXxvn2/aTz66KX53VeAbif63ZfqmjN70+7zJmPr095XbidLSZPJ7y8fs248X5R2Ubm9joTi4HDii4HAJWW/0VXt9S3LyPH7+NZuh7mlc+jPF/F62ZtaRor/sWVV7ujZXFdEgAnceSEPvlxnAwuE3/mFYX9KXHkuFUvRwxZz+WysA3K15kA9FuBe+3HKZXjqa9r7ACz9HHdbb3fqPowBN/eaZkj3HDwH60Mu4L/hbablcjdCee25+fYH6IxZCeTwd0ZnrvmRvGlSQLq7gqSAFZZf0IMWQfp/TnEn5NiyLDWOK9HMnCwiPPxXPswoJcY0eDHC2aLT5r/MyrOCKaij1d10rJS3yxn2DwE/8zrnBfbcL3RjSfON+OyPG1s83oiqMGE/aGX+8Oz9v6QY/+/bnSH7GWiYL9Znn1olNtTEsCqO7rEcXSu1+a5NOriJ633iVEo7+V9AvV5duHm7/h8WfeQ3BCnXiqPyz0xzFp1ctf3w2ZlDIJ/rjCA2y3bta4k6kqsPcTjuF25NYY1jsRhY8ol6sN8XXRpjJ/jOfcMH7l8SUBzH3JugLnbyvPu0RnH5WEjjlyY43IpE4BsZduapRLMfrXHQ0nZ1xH8M88koPr0UvJHFUl10j2x/dhZ7xuVzzdKWBIAXFldPKxO7tGa6Rxcnh/P6y3Sd7zOLkAvGwFO09dZgPtTCna/+bxWoPTT+7Y8zqxtPONnG9vd6cgJ50DwD5yRKOoOxCqKeHBSN5t+/u1K4s983tfnnEjx+0UrzOtMAHYbG6qpV53SKtZ47aTXVe3X5g29H80M1xxCb4r428ixxQ0H/0b7ASQBMNmLiNeao+7UAx7k3y4df2bwf78d/M8w0s/xrMQSgMknrHryrrutAo3CftYK1jdbmVp00em1Hov32W8F+nEi/DF3kEG9RGBVdwXKG+baO8/DqnWTBwj+AUkAdEbs8+Pq4wkO62HiR406db05sMs54s962O0XzRgy32vQiFGHzUbl+p7TRRt17dY1BzrHU5tnlhQbJFrrnzULrTHN80fBUAZJEahv5EaLEX7uNAL72GjDaTtN/V4xbmtmi08zWYhLQE+aiQQI/oEFOYes5zkkApwNN3+z5Pt7P2PEuoW+vt+z2Rh8kHHiV3VX8LPiz8b9d9N6ixy/Vz7vIJOOSB5+WX0YpGGhjr1bN7QB44QV2dduu59+bqThlL/1Mtsbt8fBPmNCr3FjJ+g1ssd9gT+Cf2AJgqIIXqKeuyMJYMn397VGHLfb3t8zHuxNapGfFn9mbNib9j+b9wnk8baex9v+Oe43BZzAji8hvjcJHXBVSUDOTHqwqDOTAsAyV9RbGfwPlQZwheeWXiYA7xZtUiIAWOYKelvwD8zxHBM3Ke5lEuA8AwCCf2DFzje6GQLADVTE0SK3ky1yAyUCXNO5p77XaO+MuXCAFXNLEcB8g//qZHSOXnUyOodRp4DrPAfVwybGuehBewQ9YDX9TBHA3CreCPr/TvAP3JR3v/v9my9+8fMY+zwSgc3y81pZ/lt5/A9KB1aXKwAwn+DfuNxA185L9aRhwdUAkAAAV1jJDqqTS+6Cf6Br56foChQ3Bg/Lsp+JwFjJAMDFK9f6prttpQF0+Fw1KMvrPF9tukkYAC5WoW4Zdg9YsPPWoxyhzLwBAHCOCnTNGP/Agp/D6gaM185jsPzcAwCXrDirD8N83nv15u2+UgEW9HwW57HH1cn9AeOyPCnntJGSAQkA8KGybI7088Awn4BEAJAAwPJWkFExxo2++9VJy7+RfoClTwTKsut8B8AqVorbbvYFVikRaJz33hk1CBZfp2YCjhPMF7/4+Y/vfvf7sU1DFyvBsn9Gl59/U5a/ffXm7TOlAiy7UicfleVFOf/FOS96DjysTmYVjnPi9zHbsFJiReKAmEn7P5Xlnxc9Vv2sS4VandxIuZejEWhdoEsHfcygeZC/3jGDJrBqoutPWTbL8sfl1wdlifugDsr58SDnQFFvs+xxwOvqpEvcwrvVwQKOgq27Vmy48YgOJKbR1z8O/KfVyc1w+r8CVD/NfH4/g6I4N0bjyDODIrBE+3gv44DY10cZmy58HHCro4UdQVfcePSoMlU5N5vt1zP6PtDqD3BqvR1JQHQPioApEoDoMuSmYRZ5v97MePQwA//9Zflutzpe8JFtxdWAuMwYLa+bdkeuqSKrW/13M/hXgQHMXnffz3NonE9HZXmhEYUF24e3c/99tozx560F2RCb2apwlMHYvt2TOe1rjzLbt68BXO58upZJQCQDgzyvRjLwXBchOrzPNhsAN5a1B8qtBdooverkaoBWWeaV7f90taksT+1fAFdah9fJQJxnx1mXSwboyj7abADcWPYrVrcWcAOtZ6AWWVp0C3pqt+WKEsv9yv0mAPM+7/arD12EepIBOrA/blcr1gB4a0E3VgT/daZ2mEGbkwbn3YfqG83Hle4+AJIBVjUOiPp/Y5X2uVsLvvHixBGtt4PKEI2cL3l8mA8Zahagu8lABGZuIOaq97XmKH8rGQfcWpINOaw+zB1guEZmCfxjeDr9/AG6ec7uVR/fM3BUJwOVoUW53H61dGP6r2wC0AjwIgkYVvpyI/AHWKagLQK2u5kUhMNGMqCrELPsR5vVko7pv9IJQGMDDypzBzjQTyqMh9WHKbsF/gDLcX6PJOB29aGrUD0D8cuy7Gv8Y0JcuNRj+ksAPs30zB2wegd6vxH4jzPwHwn8AZbynN/LRCASgkEGedHCu99ICJz/V3Pf0DNkFROAxonB3AGrcaCvZ+A/yMD/iZt7AVauLhhUH7oL9fPhurvQvsbAldkPVmpMfwnA6cGhuQOWM7sfZuDfyyTvmRM8AFlHRDJQXx2oE4KoI17m+lDD4FJtc6NDSgAmngjMHXDz22Fw2RNuq5tPPbX8M5f2ALhAQtDsMnR4VXVJ1lVH6qZr27YrO6a/BEB2uAjl/76s7py3lT4P7rqbTz9P2NG/33BwAFwmIeg3koIwzjqmTgj2L/j+e/Eebjqd+3bUy0MCcK4dRv+wBUgA8opBPTFMqLv5yOwBuOo6qpkQxLqXfzpsJQWHM7yXBGC+2yq2jTH9JQAXzv7dId6xBCAP6mEG/r1Kaz8ANxcn9KsPXYfi53isnphs6lUCCcBct0uUqZEeJQBXkvEbI/YGE4BGF5/71YeRfOrWfkkZAF2px+qE4Otc9/JPdUIQ62i42pEAiNckADJKJiQA2Xev2cVnVJYXumQBsCD1Wq+VEPQbf466Lq5gm4/gcmWsx4YE4FoO5Jn7lEXSIAM9dwLwpCxfZtAfB3UE+/WU7k6QACx6sBoxxNbaF1/0jv7pn+o/1aMNfVuZsbguq61SDhtnPMc9mxKAa90pZ7qrvDzvdeVqwazleTez9/pEqF8/AMta7+3d+ff/YXDn3/9NNX71XfX6u+/K+u+P12lcfTxj8XgFy2hqDGXUxvn6XBFMFhlm2fn2M+vcKj/frVrjymaW36tOug1JAKYH/cct/b/+zW+qNz/8EH+6J4MHYFX0fvunx0tV/c1J5H+SEPTe/PAPw9evvhv+4X/9r6gzIwE4bCQESz3aXfm+w4yh+s0YasKY/t8Y+U8CcN1JwPHlprIzPs8s9KD83MxC635+69F1aNUv5zVu5K2D/uov/vKvqt6f/tnxeu2LL6r/9/+KXKmSwQMgIUjROFYSgd74u7/vlfV6JgTNkYaWMSG4n+sYXelpxhHN3hcbxvSXANx0IhAH3Z1GP7QI+Deqj2/0ifsG7qxg0N9rBP2DOuj/87/6V8frf/FHf2QHAoBTxBXyWP7tv7vdTAjW3vzwD+uREMR9BI2EoL6HYH+BY4f16sPEa31j+ksAup4IPC076Siz0xjeq7lzDmJ4qlW4FyCH4apb+XsR5DeDfgDg8glBLRKASAjG3/19nRA8LnVxVX18hWCR4o+txs8R/B9kTHXHPZUSgK4mAcdDhDaG/WqK7PWrJQz46+9at/If9+f/6rd/Wn3zr//3j05SAMDVii6033zx16XO/etmQlCVhGBQ1oNFSghyyPXehD/dMTKSBKDLwXD0S9uZEPwfZ7HLMixoo5U/1sddnepW/gj842QEANx8QhD3DDQTgjc//FAnBIetpOBGu9XkyD6PJ32l6uQ+y3uuAEgAuhgU9zP4753ytDjodhftZp1GX/7bVauVv76JFwDonrorbl1XNxKCfln3S0LwKOv6m04Itk/La8qyF/dYuvlXAtC14D+y1nF1MnHVj2X5ZfXxjcCDXO+U53/T5RtYGiP21AH/cV/+CPi18gPAciYEb/7xh/7r77679oQgJv1qxEzjXEL8z2/z55hN+W7c8Fw+y8iWlADcuGzRvzdjYN2vWuPadiTgHzQC/uOD8Ks/LQH/X54E/PryA8DyJwTHEfjJXAT98au/bycEzbkIxlf4MeI9nxvTXwKwrInCUVcC/xxm66OAv+7W8xd/9a8+GoMYAFgdUyYnqxOCYcYR4+qKZis2AagEgPkE+80W/n7+XDX78cfauPwAwIwJQa8kBMN5JARIALh8wD+oWi38An4A4KoSgsZsxcOyHuZsxc2E4FD3HgkAVx/w96qPW/h/6sP/63/5GwE/ADA3U2YrbicER9XHVwgkBBIAzhnw91sBf68O+Hu//bNc68MPAHQmIfhptmIJgQSAs4P9ujtPv/rQped4GM4/abTuG6UHAJAQIAFYzIC/bt3/Ote9+uA5CfT/N+PwAwBLkxAc/dM/fZQQxO8SAgnAMgf7EdzXffZ/at2vJ906DvZ15wEAllg0an7zxV9X3/zrvz4tIahaCcG+kpMALErAP6g+tO5/1Hc/btbVug8ASAgmJgQxW/GgrAfl98cSAglAV4P9Zsv+TyPzxE5dt+7/yW9+o3UfAEBCIAFYwGB/4o260ZXnJMj/s59a+Q3FCQBwtQnB//zHH35KCN788IOEQAIwl4C/GezrygMAcIMJQSx/8Zd/dfx7jCrUvEIwISE4zKTgSOlJAKYF+xHcR8Bf99sf1DtbPQynrjwAAN0QvS0iGTglIahjvMPq46sEEoJVTQBaN+rGOrr3/DTJVgxbFYG/1n0AgIVOCPpl3S8JwSMJwQolAK3W/Vh/cqOuSbYAACQEq5IQLF0CMEvrfgT8btQFAFjdhODNP/5Qvf7uu/74VUkKvvtupRKChU4AprXu1114tO4DADApIYh7O0/u7/yb48fGr75bmYRgoRIArfsAAMzDjAnBuJUQjCUAVxvs9yp99wEA6E5C0CsJwbCsh4ucEHQmAdC6DwCAhGBJEwCt+wAASAiWOAE4bVZdrfsAACxjQvA/f/ihN/6uJASvvhvGyENdSQiuPAHI1v1msB/Bv9Z9AABWLiH4t//u9vHvMTtxSQQ6kRBcOgEoH7wO8j9p3f/1v/xN1fvTPzOrLgAAKy0av2OZMSH4NhOCwxtPAMqHWqs+9Nm/XTVa9/+kDvbLFzu5/AEAAFwgIagyITiqPr5CcCUJwednBPx16359s26v/sDRjUfrPgAAzC0hWCsJwXpZr19lQnCrvNHmhMfr7jxrzQfrm3bhovb+69/FalSW75UGAEvsfombeuImrsqbH/6hqq8MNNQJwbfnTQD2Gr/38o2OFDNzMijLoX0MgCXXz7purCi4pv3tUDHQSSXhfJ/DwgLAMtd3e1N6WcCN+0wRAACABAAAAJAAAAAAEgAAAEACAAAASAAAAAAJACtivzIHAADLL8ZkHysGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACA2dxSBFyF3/76V72yWn/15u3Tc7ymH68pyy/L8m157WiG/zEoy1pZDsvz95U8S3xMxX4+LMuXZXlZ9vddpQI3fkw+Ksfi5nnrxjyOvy/LqLz+6Iz/EfVcP+s5xz1z8Zki4LKBf1m2yo+vy3L3HK+LwGa7LJEwPCnLw/LY9inPjxPuTlnGkWQI/lmBhPog9/eNOLZOOz6A+Qb+ZXmU9dzjc7wugvi9suzncfx1/J5B/qTnr+dxvxZJhuAf6PSJMdcHZdk7R9LwPoOc9mPDCc/fzvdfU+KsyHG11wz4MwB5nwECcDP13FYch+d43et2nZaPbU947rAs7zJpgLlzBYBLaVzKPDrHyx7ma8eN94mf96tW60peXYig595pl01hiYKNCAAGZXnROs4+OT6Aa63nfjzHcRz1Vq8s49afnpdl2GoAi+dGUvCg/K9DJY4EgGU1rYUjTny9+sRY1hEExWXXZ81kAZbc/SlJdRwffVfCYOHruSqT/PrqQgT/+7r8cJ0+VwTclDjxtVr169aVSADGjUBoP68E9PPk+cTVAFYgcGi3BP7Y+Pu+YoKF0E7Yjxr1XFjP57zI+wzuZv33zNUA5skVAG5C3bVh2Hr8dusEuZ4/D6qTG4Wf5Wv2tIKy7CS5sNDq1vz7rccHrXquHjzj63zNRib5e+4HYJ5cAeAmjPKk+Lic4OIkOM6TYixH0eqR3YAiyN9vDLm227hcGi0lm4oSgA4m8FGPRV03zFHs9jOw/+nKdq77+fwH9WvL8+9VJyMOxZXvO0oTCQBz1+hqc5aNi16ejJbN8n/ipDasTlr916oPNzjWrSa9XLf/x24mALdtLQA6nAQ8KHXdy6yvBll/xSAY40b92WvXc3HPW3ndYfXhagFIAJj7CWvjmv5PtPw/bSQem/njk/x79PtvJgLN5MGGYplFwDCIy/+tJPvr+thQRLAwdeqoOrnqXQ9sEXXag8ZT4nie1OimCyBz5R4Ablz2c4xWkSet0X4O2wlA62+wjOqrYO2goNf4G7BY9Vxc6Y4r7PutWe+jLlubcl+beg4JAAth6uyGObrBtJNidOnZnTC9+vHNUNlqUj8/gqBoGXmmuFlG2eofgf7d1n7ft99DZ4L59mNRV22eMkDFVq7vtR5/knXaowkJ/xOlzbz8TBFwBSfDCEz+j7L8+Re/+Pm37373+//e/Ht57P8rq2FZPy9/O2omBhn8x3Bnn5zoynPH5TUx9OF/LOv/VpZ/kc//Pw2PxjIr+/p/jYCgrP+osd//P8YJhxur5yIgjxt4/7ws/1yOy8NSR/2hccz+x7L6v8vyP8rjh836MevA6Pd/r/marOf+EPVb+fFxWX9fljdl+S/l9xetKwVwpW4pAi55UozWjnZXhXGzK0/d77E+meWJdFidXN7cP2u4w0wwIlmI99w1PCIrdHzVs4numgwPbvRYHLQeOmo2RGW9FsfrKO9VW8t67ijrufEMCcZ6/up4BwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAWBC//fWv1soyUBIArEi91y9LT0nQRZ8pAq5Jvyx7igGAFbFVlqFiQAIAAABIAAAAAAkAAAAgAQAAACQAAACABAAAAJAAAACABAAAAJAAAAAAEgAAAEACAAAASAAAAAAJAAAAIAFgUfz2179ai+WK37OnZAHoaL3XW4T3BAkAc/PqzdujstorJ6/+FZ0Et8pqqGQB6Kio8x5dYfA/KKstxYoEgEVzWJaDchLbvMwJsCwH5cc4qY4UKQAd9SwC9lJn7V225T7rzb2yvFCsSABYxJNheFxOZq/Lsn6Ok1+vLNt5AoyrCPuv3rwdK1IAOmo314MqG7/O2xU2G71eR73Zek+QALAYSsAeVwDqoL1Xlp1sGRnMEPjHCXDY+NNzJQpAh+u8qO/289e1DOJfz5IIZOAfDV57WV+GUXanhWvz+WlBGpzDYeNkFmK/ihNdPP6sThByf7tfTe/nf2SfBGAB6rxmXVUnAg9LHRZ13qgV+MeV8Yet19R0/+Ha3So75XvFwDU4yhMkAKyCcS696uMGMpiH6EJ9Z9Ynf16efEuZcVk5es+kURH2q5MrAMejBZXlq+qklWR9SkLwxy6FAtDxOm9YVttTgv76CsBOWV6WZaM6af0fTnm7+PuhUuWSzhU7fa68uCLrEwL/JyWY38+T5SDW2XfyQfl9IxOGh61EYL0yChAA3XZ/QuAfdd6okSRUWe8dnlHvfVme81SRcp3cBMylZXDfy1/jRHcnLkPVwf8k0cpfls3q5IpA88T3UIkC0OE6L+q7Qf4ara4bpT77qhn8z1DvNZ+7rlSRALCI6paQaP345rTAf8oJMVpFvsnkoW9GRAA6bJjrGLrzq/O23me996BR7/XOM3w2SADoigjYv8mWjQuJS6SRPEQSUbkKAEB33S3Lg1Jn3bvMPWuteu+uYgUWRox5PMsEKDn28fsZ37OnZAHoaL3Xm/F5eznT7yzPHShZrpObgLmUeYzYYyZgADpc743n8J77SpbrpAsQAABIAAAAAAkAAAAgAQAAACQAAACABAAAAJAAAAAAEgAAAEACAAAASAAAAAAJAAAAIAEAAAAJAAAAsLQ+v4l/+ttf/6pfVmuthw9fvXl7dAXvvRXr8l4b1/A9emW1XZaN8v8O7U4AQJfMK+Yq77uWMdCL8l4jMddiuckrALEh93J5WJZHZeMelOVdWTYv8b7ruVyHOKgGuQYA6KJ5xFxrGW/dFnNx3mwudrz3rce247Gy7Fw0I82sdB6fd31KBszZZTdob2sAWOJ6b++SDZqLEHP15vRZe3nlQsw1Jzd9D8CkSzh11531i+xYcUnrKroSTQn++5P+n90IAOi4ecRc4zl91rhKsSbmmp/Pu/aBYuOWnbD+dW1KMF73ZxtP2vni75P6h0UrdH0QTNqJMrOM9z5qvj7/X/Q7ezYl+/0k6Wi8V3yn/Sn/ay0+f+O5h3ZuAGARYq5mLHOBmCvip177fcvjw7J6VJYXYq756dwoQK0d5rC9AeNyVXXS5+xuWV7XN/3m3zbL8rr8utXeefN1g8wqX+cO1nzOVr7u+KaW7BtXdydaz8dv5/8Y5BKXzF5XrSsDeclvO19zPz5TfSkrL2tt5+uG+TkOqpN+ea/ndTkNAOCKYq6PYpnzxFz5vjv5uoh5DvL5dbJR31NwP2Ounphr+Xa8vbo/Wu4Qw+yjdjCl79dHj8eOmH3XBq3n7DWzxTqYbz3nXf1Yvs924+/r+b6bjcfet/vyNZ7X/P+xsx60nreT/6/X+K7v8/sPGgfM+y71F7zqk4x7AABYsRhns2Of58pirkYss3nOmOuglRC033fQjq3EXFfv847slJuZXcaO9mDScFK50Y4v4TQ2/ve5jtfu58/tyzmxk43jvRuXuQ7zf8VOFq+LS0136j+W/79bHt9ovOc0ky4dRbbb7ir0JDPox/n96ktuL+tLVZF552NfOm0CAF2OuVrdh2aNueL3fut/Roz0YzX5HgUx15x0ogtQ2RCxMz7Ijft4yl3egwkb/zAD9+envP3tCa97nq+rd8pJn+npeceZbRww7feq32eWoat6Tk8AwCrEXDl4y+Z5++OLuS7n8w7tkIfZ6h5dcbZy5/xkQ026uWMGa9Ned439v/arKTfYAACsQMx1XV9RzHWGzzq2Q47Karc6uVGjPeZ+ZIZr7T5huUOdNvFXvG7QznCzT1gzw71/zved9r/CpEtKa9XZl7cAAJY65moPxCLmWr0EYFJ2FlnouDoZiad5+WY31zutm1KG1eR+YbUX+X/2Wjtk9A2r73of5wEwaLxv/NxrvdeX9Y48LaPOrHO9dQNMPeTUc7scALCKMVfjfR83e2DkfQnj1nv1xFxLmAA0x2xtbrjsA3av+jAc5yAfH+eOGo8f5F3kMazT7dMuUWWGO8r/9S7vAn9Xfv620d+svvS1l3+PO8q34j6Axlsd5k4WO+mwuXO2EoWNxkFTf6+4vDaqP2fj8a9b5THtAAUAWOiYK9/3SfVh+M+9fN9fNvru1+uHOexoX8x19W7d0I4YG7N9qWe3NflWc+bdozoYz52zniHuRStIr+ohQMvjd1qPR9B+PzPX5zHSz4TP9DB3rBe58xy1/n4/d+JRZq7DxluM6oks8m8P8/Mf5eccNXa6YWOnG2dG/NFjk+7KX/CTT2y3vfK9blUAsPxBd8QjL/Om25v8HPOMuWJY0Sft7zhDzLVefeh6/aId8+Tfb+ff9sVczHTANecBoDsJgHkAAFixeGRzyb/jSo+lv8g+W9Lv9cKmBQCYu11FsHg+X/QvkJe2os/+qPpww8fIpgUAuNKYKyZOjT720fc+utuMzjtnEhKAK9EYy7a+weXeeSeTAADgTKOMt2KJewUeKBIJwE0mAU9tSgCAucZb0cC6qSQW32eKAAAAJAAAAIAEAAAAkAAAAAASAAAAoFtuZBSgnCK613o4xpGNu8uP5j2mbE5tHbMF32tPTw0AsCzEXHQmASg7wKjsELEzvs6H6gkl7pblUfnbflk/KM8bz+kjHDV2fuCGKqU4FyiJuZRtnF+HU857x5X/POdLKf9/Lc/vMUnQhi0CN0fMRRcrqddled96LHbG92U5UEJLta0H7W3NageoeZwPlMZcy3k7y3mvLOuRdMW5tSzvckbPuSUA+T8f2Qqs8PEXx8CmmIsuuul7AMYTMtV6Uq9+2SH7NhEspWGu7y9Yhb5oAe3Leh2X3vOKy53qpCVuK5KCefzTuLpQljsmaYROEXPxk67PBHw0oQJey+Ah1uNJXQiyVXGQv+5O6t+WfeL225e8skKsD4LRhL8P8qDZbzw33mff7gQzi8A/jstokd6YZ3eUKwz+N5ekwj8q3+V5+fFxbgd9coHLxFzx916+fndSVyIxlwRglkr2USNwb+8I/ay0nuXOGJe3H5bnfdN4znbuLA/y+Xu5A4cHWSFu5U50p64g8zk7ZXlRlsiIY0eLy2UPGv3ntvLxJ+X3u/kZYud8XH6/Y4eEmY7xYQad38YxXJY45jdPeX4vj7s43g4n3UTWqBiOppw7etWHqw4fVR51BRctYY3nfVTR5Xkpzj2jrJDG9f9o3WD3SYNDXfHlZ1vP506q6NYafx/nex2d0rgxuoI+u0enbKO6Qh9N+BzN8py2TX6quE/Znp8EFBO2x3rj++pDDB2JufJYjZt7n5XHNvO8EXFT/RZf5XlEzNVBn3VkBxxk/9Tt3NmeZrDerhRiR4udYz8rnCfVyWWrYaPCiZ+fZ8UTFXFULnE5+lbsVLnDTLopbSuf9zQvX4/y/WOHH8SBUZZ7+dzb8bdIMsq6Tj4GdieYyf2sUHYzwLx/RrKwlc+NY3enTvLr80L0s82K4WkemwfNS9n5Hg8bQXg0Cmw1/hY3xm3lz3v5ebbrPrF57jmc8vniObej8sv3Pqi71eR57SCTnH6ub+dnOWg0TNQV7U5WjqM8D+61nrOdlWl8jy+zsjx3F558z7rMn7X/Vpad/C77+bzXWRk3y3Mny/tpVsbvc9nO+zt2siwHE7bndr53fM/7eU/C2oTtUZfZ17kP7Dh0oDsxVzbe9OskPteHmUzcyrhJzCUBmCkoiJ0qugNM6hIwyEopbljZzMvxXzd2jmk7xI8RHJx2s2HjEte3rT/VLVMPW4+/rD9f43N+bXeCMyuefp70x3nsRKXSmxTI5jG7lRVQXZHsVh+PbrOTx2PdOvwkk4H1xns8zHPKfvZ33c3zSC8rnfo4j+P5q1jyOVHRPcrKqW5p+j7fZ5zv3W80ODzNAP5h/t6s+CJJuJcV2IP8jP1WRfukfu8MzPuN52xmuY3yOQ/yfLg1Y9F/mZV+VNgH+do7E7pHRnk+z/sF9huf9XHjs8b/fJHlUpd5+CY+V6vibm7PXgb0sS0O83vey6RmpxFA1ElJL+8jeJABigofuhVz3Z7wvvWVzrMaIcRcN6wTXYCyb9dhVnbR8rM/4dJ2/G2cLW1Nzay1fs3ahH9z2ji3/SmfKyr5o2l/b1mzO8GZHlYftzo/yYrgYfVpX/Tj+wRaFdODOuDMgHLQCECPr/qVx79pnAvisvFRq//+WqOCG2UjQdXqirKRldjdDD6nnbee1OeWZiv5BC9aFWR93tnP/7PWvJyd3V8OG4/F5zhsfY+jTJ56M3YFqoP3wwz+jyZUylEmLyfcDDhofOa1CRV+ld/j8IxtXzWTjrwf4Tipi//Z+ttuqyGnaj8HuPmYK84drfPJWceomEsC8NGGj4ogWoOidSpag76Z8LTBGe8R/cYiI32c49rWQcSDS/QdVdnAFWi0+kTQ+rh9bE8IZHuTzhONY7I35Txw2KpoXk6oxM46H0VFND6rksl+r/HZ4zzzbXW+ca7Xzvge+63vsdF67DzfKa5c7MYN15kEPK4+vSzfr79T6/HNCY0sv7zAuXJapf5t4++HM5YZcMMxV55D+nlOqe+7rPv6X5SY65p81rEdMjZ83cesXQnVmeajCYFF87EY9i5aGB9lsHFnhsmG6kpt0iWlnh0SrsSjTMbvNJfqQ4vSwwsGfGc95/YFP+/4rIA++9Dez64vT6tLTHTT7O8/xeCyG6DVBWp9yucYnJYY5fYaNq4SxHYbXWKGz8PWeRhYgJgrG2Ti6ubLfH2cE76a4SqdmEsCMLGFL3ai/eqkFb9ZQdU3DG7lxBVrebPZdv6t2cc0Lmdt5s0l4xkOgnhOJAnrrZsHe/kZn9lV4NLuTkrG87FxBpVrrUqi3w5Is8W936hEHp7SKBDPGbS7tOS546wbaHunBaX5P4bZR/0yxo0E6aOEoHGzXTznfjtJyIm9ztsqXt8/sD2hvM8qzzpgf5blGufrJzOWQR1Q9KfUAyp9WKyYK85P0T1xN2OumUbqEnOteAKQO1qv8XPTvUYFtZn9Po9vWqs+3Pj2rjoZLeJlI8ivK7PXjVEp3ufsd83s9m5rHTayAtrJu+Mj6IjLYg8aw/3Vgcjdur9v47HeGX2AYWU1hsKcZjeP32Hjsee5rkeEqY+3+42bSPcbgWgdNG82gtn6PfYao/McX7Jut1i3RtypK6Inrc/5deM91hrPrR/rNT7HrOeDuqJ92PiMvTzP7Te+Ry+/R79x3vz6jAp3UoV/lOfY+ubjZqW8n5VyPZrPWgb/zf9xPFpRNrBsnqNPfl2W7e5fcYXGEJ+wmDFXnI/etWKug9b/EHN10K0b2hEn3kjWrEgaFXBdOe03Kulhvn639Zr6Jrfn1ac3q91v7OTNnWbcvEowbTzxfO92y9Vh+zHj0k7d5nHQ7sXQYEpjJSueeujHjSlj9NfjPcdxF63J+9lHdbMVMH40ek2+dq91TI+aLdKNoSenvUf9P0aNIDUqometuQB28jPuZwvVUf7v+r6ElxnMxr4ewfFG6703Wt9pNyu7o0YZNc9ZD1r/f7uVIE28mXdKue62y77xOervU5+79lrntXZ5vqsm3wi82/iO61mG8Z736s/Y2BbRFelF/p/79ffI8+x2fuaNvBm6+diT897PATd47turLnAP0gLFXPXx+6L6uIti3CO0HiOqtd9XzLXiCcAcd/KoJF5MmZQmWrEO7SwSAG5k2089YU+pnA4bQWM/g79xNXlyrOYEWtNm/u5lJfbJRGGNQPirRoA9baKuQX62cet9d3MEovi9nzfctiu+owkNEEetZKZuNZs2m+YgP8PEsmh91kk33X40qlJr20ya3OxwwpWSrUx21lrf5XYGAqPW/54UaEyc2G3CvjDO39emBRAgAbiR71ZfQZzYCBENJpOGA4Z57Iy9vPS0NSn7bU4exM0EgbF9lAQd3Dc37Zszl9Vw0jm2cZ7dVErwIQFY1mMizwXvp8zhsu5c0H2fL8sXyWH74vL942zt388/reXPGzY3wKXElZJxe9zvbA2Mvz1QRLAS4spd9OnfybH7Dxsx13Nd9SQA150EREte9C1tXn4+dHMZMEkGrrfzZ5NMnS0aUuJq6rucSOijSdqca2E11IMJTLt3gO7TH5vrCrQGlXsA6OZ+2aR/+fnLTZnB5ONkae8BYPF9rgiAVaW1SrkBrKLPFAEAAEgAAAAACQAAACABAAAAFoKbgE+RQwTGnAIxk+XTc7wuXvNl9WH2yo3ThsfLWTe/Lsv3ldmKAYDliaV61ck8Ic9njW9yVvD7+WsvX7t7xv+IWOqXEUudJ2aTANDemdZz54v1k3O87njG4bLzPWj8HrMBfjJddiYYO9XJdPcbxtAGAJYolnqUsVQE9M/PEfxHbHQv5mbJ3/dyAsLRhOcPM8F4oAF1dp8t2Y7Wi+mnJ4ztfW6RaZbl3jn//yAz0OaswxuZvT6aEPzHGMExhrYJdACAzgTuGU/1LhlLPW3FRLPYKstuPTFjrp/F4xk7tYP/rUwWBP+rlgBE4J0t7bHcZBea+7mzHjV2/nqK7Iet50Z2u3aBAwMAYG4ycB9HXBXxVbbCX0c8FwnHoCw/tv60nzHTsBn7Zdz3xCzu53djXYAa/bUmmanPfWZ+93MnfdKB2SinZcovY4eOAygvZw1zB9fyDwDMO+baPOXPTyfFItndZpTB/8OM255P6oYz7zgqGnbL/48fb8fnzYcj+B/r779gCUBu5F9m4H7U2EGP+3HNkDjcrU76k93rUBBdf49eKxn5Za7rS1d369/zykV8n8NMCGSxAMBVup0xxjjjlEF10g15dFYMlXHJg3pglLI+yPhrNIf4q36/L1ux31ozjsqkJOLB3fJzdAEa5O/xmfSsmMFNdgGKjfysEfwPMvgfTcsuM0E4yIzvm8j6OtaC/iLX7e4+g1zXScF6/nyYNwv/cT6+d9n+dgAALc8awX89AEkE9jMHyxFvlWUz4q986HX2aLgymWzE51xvxUOT4qj69yf5mTYyQdmyuc92Y1cAmi3ds+6MseOV58bGjktRXzd36C6IxKV8rtu5A8Z3+j6z2H4mLeNGP7pxfa9CHFTl8SdZBqdeAQEAOGd80hxCcyfX5+6GnEH5wwzII17bncPHjQFY4upENIo+zziqTgZe5vrrXL+ov0PGYA8zButUfCgBmG7mnbHRJy12vq3sE/bsrBt/MyOc5SaWjct0w4kW/djx8n/V7xMZcnv4q5ftgzO/S89uCQBctexJMajO2eU4Y666d8OzeXazyXslv6o+tPpHH//oLn3UTjgmxH4vMv6KWEoC0OUE4KI7Y270/XqCibKuJ5kYTXn+tfUJy+9xmN+v3mmfNnbs+PG23Q8AuKZ4K2KtU7taT3jNMAP/4x4a19Wqno3Bu/kZIs47npOp0Uj8bXXSTWhg+M+L+ayLO+N5xvGPnTH70cclo7ip9uC6hqua4fvFDjusPr2ysVtNb+l/YbcEAK4wHpnY1XpavBXxS97sG6+7k/MVjW/o4x9/7ugG3oqjqimx1JGkoMMJwBn9/u9eJFvMm4LjRpCr3km/nPId4p6EnSl/iyRkK4P/dh+5jUxWhq2DMMpiZLcEAK44iK6qTxskp8Vbh3MYbKU3IVaKCVz3ssF0UixV96K404r5Il6KnhWPWy+5XZlfaSY32QVoOzPL40Qgu8VU+fuldrar2FkbN7qEyITjctNowoHTbw77mYH88eNVTmM94fPFzcBxxWI7E6GjfM098wIAAFclZvWtTrpaR3yx1Yi3qurDzbyfxClX9L8HjVgqBnCpWt2P+vkZ4v81u/zEYzHP08vs5TEp1tvIBCKSm+cZR72Y8zwFcHIVozlMVe6Ig/ZU1WcdHKsw9Gd+z/f2GgBWJEbYO2MCLqqfeky044X+OV7fO0/XcUACAAASAFbOZ4oAAAAkAAAAgAQAAACQAAAAABIAAACgW/5/AQYADaAsqI/oY10AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 34,
     "metadata": {
      "image/png": {
       "height": 100,
       "width": 450
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(filename='CFA0129-R-EXH2.png', width = 450, height = 100) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example** Since it's one-tailed, the rejection point will be indicated by $z_{0.05}$, as assumed the test statistic follows a standard normal z-distribution, it can be found that $z_{0.5}$ is 1.645. We reject the null hypothesis if $z > 1.645$. The value of the standard normal distribution such that 5 percent of the outcomes lie to the right is $z_{0.05} = 1.645$. (Let $z$ represent the calculated value of the test statistic.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4\n",
    "rej_point = 1.645"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**<br>\n",
    "Note that the arithmetic mean equity risk premium for US stocks relative to bill returns, $\\bar{X}_{RP}$ was 7.5%, the sample standard deviation was 19.5%. In order to calculate the test statistic: <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.017951205049017305 4.177992496615468\n"
     ]
    }
   ],
   "source": [
    "# Step 5. Know that\n",
    "mean = 0.075 \n",
    "std = 0.195\n",
    "n = 118\n",
    "# Calculate the standard error of the sample mean\n",
    "sesm = std / n**0.5\n",
    "# Calculate the test statistic \n",
    "t_stat = (mean - 0) / sesm\n",
    "print(sesm, t_stat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We reject the null hypothesis in favor of the alternative hypothesis that the risk premium on US stocks is positive.\n"
     ]
    }
   ],
   "source": [
    "# Example step 6\n",
    "if t_stat > rej_point: \n",
    "    print('We reject the null hypothesis in favor of the alternative hypothesis that the risk premium on US stocks is positive.')\n",
    "else: \n",
    "    print('Fail to reject the null hypothesis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4704667741120198e-05"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Using Scipy to obtain the p_value (TBC)\n",
    "p_value = scipy.stats.norm.sf(abs(t_stat)) # z-score\n",
    "p_value # Correct answer is 0.000015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scipy.stats.t.ppf(1-0.05, 117)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example** In the sixth step, we found strong statistical evidence that the US risk premium is positive. The magnitude of the estimated risk premium, 7.5 percent a year, is economically very meaningful as well. Based on these considerations, an investor might decide to commit funds to US equities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis Tests Concerning the Mean\n",
    "### Tests Concerning a Single Mean (T-Test, Z-Test)\n",
    "#### T-Test\n",
    "A `t-test` is a hypothesis test using a statistic (t-statistic) that follows a t-distribution. The t-distribution:\n",
    "- is a probability distribution defined by a **single parameter**, degrees of freedom (df). \n",
    "- is more spread out: It has a $\\sigma > 1$ and more probability for outcomes distant from the mean (it has fatter tails) than the standard normal distribution (z-test). \n",
    "- As the number of df increases with sample size, the spread decreases and the t-distribution approaches the standard normal distribution as a limit.\n",
    "\n",
    "For hypothesis tests concerning the population mean of a normally distributed population with **unknown variance $\\sigma^2$**, and if the population sampled has *unknown variance* and either of the conditions below holds:\n",
    "- the sample is large (>=30), or\n",
    "- the sample is small but the population sampled is normally distributed, or approximately normally distributed,\n",
    "then the test statistic for hypothesis tests concerning a single population mean, $μ$, is\n",
    "\n",
    "$$t_{n−1}=\\frac{\\bar{X}−μ_0}{s/\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 4</div>\n",
    "\n",
    "where\n",
    "\n",
    "$t_{n−1}$ = t-statistic with $n − 1$ degrees of freedom ($n$ is the sample size)\n",
    "\n",
    "$\\bar{X}$ = the sample mean\n",
    "\n",
    "$μ_0$ = the hypothesized value of the population mean\n",
    "\n",
    "$s$ = the sample standard deviation\n",
    "\n",
    "The denominator of the t-statistic is an estimate of the sample mean standard error, $s\\bar{X}=\\frac{s}{\\sqrt{n}}$\n",
    "#### Z-Test\n",
    "Practically speaking, the two tests are similar since the t-distribution approaches the z-distribution as $n$ gets large. \n",
    "\n",
    "1. If the population sampled is normally distributed with **known variance $σ^2$**, then the test statistic for a hypothesis test concerning a single population mean, $\\mu$, is\n",
    "\n",
    "$$z=\\frac{\\bar{X}−μ_0}{\\sigma/\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 5</div>\n",
    "\n",
    "2. If the population sampled has **unknown variance $s^2$** and the sample is large, in place of a t-test, an alternative test statistic (relying on the central limit theorem) is\n",
    "\n",
    "$$z=\\frac{\\bar{X}−μ_0}{s/\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 6</div>\n",
    "\n",
    "where\n",
    "\n",
    "$\\sigma$ = the known population standard deviation\n",
    "\n",
    "$s$ = the sample standard deviation\n",
    "\n",
    "$μ_0$ = the hypothesized value of the population mean\n",
    "\n",
    "The following rejection points are commonly used for a z-test.\n",
    "- $z_{0.005}=2.575$\n",
    "- $z_{0.01}=2.33$\n",
    "- $z_{0.025}=1.96$\n",
    "- $z_{0.05}=1.645$\n",
    "- $z_{0.10}=1.28$\n",
    "\n",
    "### Differences Between Means (with Independent Samples)\n",
    "#### When Variances are Assumed Euqal\n",
    "When we can assume that the two populations are normally distributed and that the unknown population variances are equal, a `t-test` based on `independent random samples` is given by\n",
    "\n",
    "$$t=\\frac{(\\bar{X_1}−\\bar{X_2})−(μ_1−μ_2)}{\\displaystyle\\sqrt{(\\frac{s^2_p}{n_1}+\\frac{s^2_p}{n_2})}}$$\n",
    "<div style=\"text-align: right\"> Equation 7</div>\n",
    "\n",
    "where  \n",
    "- $s^2_p=\\frac{(n_1−1)s^2_1+(n_2−1)s^2_2}{n_1+n_2−2}$, (a pooled estimator of the common variance), \n",
    "- $df = n_1 + n_2 − 2$.\n",
    "\n",
    "#### When Variances are not Assumed Equal\n",
    "When we can assume that the two populations are normally distributed but do not know the population variances and cannot assume that they are equal, an `approximate t-test` based on independent random samples is given by\n",
    "\n",
    "$$t=\\frac{(\\bar{X_1}−\\bar{X_2})−(μ_1−μ_2)}{\\displaystyle\\sqrt{(\\frac{s^2_1}{n_1}+\\frac{s^2_2}{n_2})}}$$\n",
    "<div style=\"text-align: right\"> Equation 8</div>\n",
    "\n",
    "where we use tables of the t-distribution using `“modified” degrees of freedom` computed with the formula\n",
    "$$ df=\\frac{(\\frac{s^2_1}{n_1}+\\frac{s^2_2}{n_2})^2}{\\frac{(\\frac{s^2_1}{n_1})^2}{n_1}+\\frac{(\\frac{s^2_2}{n_2})^2}{n2}}$$\n",
    "<div style=\"text-align: right\"> Equation 9</div>\n",
    "\n",
    "### Differences Between Means (with Dependent Samples) (Paired Comparison T-Test)\n",
    "The prior tests were only valid if the samples were unrelated to each other (i.e., independent). Another method must be used if the independence assumption is not valid.\n",
    "\n",
    "`Paired Observations` refer to observations that are dependent on each other. *(For example, the monthly returns of two funds over the same time period are likely related due to common market factors. A paired comparison test is a statistical test for differences in dependent items.)*\n",
    "\n",
    "**By calculating a standard error based on differences, the t-test presented below takes account of correlation between the observations.**\n",
    "\n",
    "The difference between two random variables taken from dependent samples, denoted di, is calculated. Then, the list of differences is statistically analyzed.\n",
    "\n",
    "Assume $μ_d$ and $μ_{d0}$ represent the observed and hypothesized population mean difference, respectively. The hypotheses are normally formulated in the following manner:\n",
    "- $H_0:μ_d=μ_{d0}$ versus $H_a:μ_d≠μ_{d0}$\n",
    "- $H_0:μ_d≤μ_{d0}$ versus $H_a:μ_d>μ_{d0}$\n",
    "- $H_0:μ_d≥μ_{d0}$ versus $H_a:μ_d<μ_{d0}$\n",
    "\n",
    "For normally distributed populations with unknown population variances, the test statistic is:\n",
    "1. Step 1 Find the **sample mean difference**\n",
    "$$\\bar{d}=\\frac{1}{n}\\sum_{i=1}^n{d_i}$$\n",
    "<div style=\"text-align: right\"> Equation 10</div>\n",
    "where $n$ = number of pairs of observations. \n",
    "\n",
    "2. Step 2 The sample variance, denoted by $s^2_d$, is\n",
    "$$s^2_d=\\frac{\\sum_{i=1}^n{(d_i−\\bar{d})^2}}{n−1}$$\n",
    "<div style=\"text-align: right\"> Equation 11</div>\n",
    "\n",
    "3. Step 3 `Standard error of the mean difference` using the sample standard deviation, $sd$\n",
    "\n",
    "$$s_{\\bar{d}}=\\frac{s_d}{\\sqrt{n}}$$\n",
    "<div style=\"text-align: right\"> Equation 12</div>\n",
    "\n",
    "4. Step 5\n",
    "$$t=\\frac{\\bar{d}−μ_{d_0}}{s\\bar{d}}$$\n",
    "<div style=\"text-align: right\"> Equation 13</div>\n",
    "\n",
    "where $df =n − 1$, \n",
    "\n",
    "$n$ = number of paired observations,\n",
    "\n",
    "$\\bar{d}$ = sample mean difference (as given by Equation 10),\n",
    "\n",
    "$s_{\\bar{d}}$ = standard error of $\\bar{d}$ (as given by Equation 12).\n",
    "\n",
    "Use cases: \n",
    "- the dividend policy of companies before and after a change in the tax law affecting the taxation of dividends. We then have pairs of “before” and “after” observations for the same companies.\n",
    "- whether the mean returns earned by two investment strategies were equal over a study period. The observations here are dependent in the sense that there is one observation for each strategy in each month, and both observations depend on underlying market risk factors. Because the returns to both strategies are likely to be related to some common risk factors, such as the market return, the samples are dependent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**\n",
    "If the level of confidence is 0.10 with 24 samples, in a two-tailed test, the test statistic will be:\n",
    "$t_{\\frac{\\alpha}{2}, 23}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6765508919142629"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Finding the t score at 0.05 level of confidence with 50-1 degree of freedom \n",
    "scipy.stats.t.ppf(1-0.05, 49)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4048917596601207"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Finding the t score at 0.01 level of confidence with 50-1 degree of freedom \n",
    "scipy.stats.t.ppf(1-0.01, 49)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example 4, 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**<br>\n",
    "Quarterly returns for a six-year period for two managed portfolios specializing in precious metals. we want to test the hypothesis that the mean quarterly return on Portfolio A equaled the mean quarterly return on Portfolio B during the six-year period. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the data\n",
    "df = pd.read_csv('metal_prices.csv', delimiter='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Change Quarter column into index\n",
    "df.set_index('Quarter', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deleting first row as it's obsolete \n",
    "df = df.iloc[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data cleanning\n",
    "df['Portfolio A (%)'] = df['Portfolio A (%)'].str.replace('−', '-')\n",
    "df['Portfolio B (%)'] = df['Portfolio B (%)'].str.replace('−', '-')\n",
    "df['Difference'] = df['Difference'].str.replace('−', '-')\n",
    "df = df.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "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>Portfolio A (%)</th>\n",
       "      <th>Portfolio B (%)</th>\n",
       "      <th>Difference</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Quarter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4Q:Year 6</th>\n",
       "      <td>11.40</td>\n",
       "      <td>14.64</td>\n",
       "      <td>-3.24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3Q:Year 6</th>\n",
       "      <td>-2.17</td>\n",
       "      <td>0.44</td>\n",
       "      <td>-2.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2Q:Year 6</th>\n",
       "      <td>10.72</td>\n",
       "      <td>19.51</td>\n",
       "      <td>-8.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1Q:Year 6</th>\n",
       "      <td>38.91</td>\n",
       "      <td>50.40</td>\n",
       "      <td>-11.49</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4Q:Year 5</th>\n",
       "      <td>4.36</td>\n",
       "      <td>1.01</td>\n",
       "      <td>3.35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Portfolio A (%)  Portfolio B (%)  Difference\n",
       "Quarter                                                \n",
       "4Q:Year 6            11.40            14.64       -3.24\n",
       "3Q:Year 6            -2.17             0.44       -2.61\n",
       "2Q:Year 6            10.72            19.51       -8.79\n",
       "1Q:Year 6            38.91            50.40      -11.49\n",
       "4Q:Year 5             4.36             1.01        3.35"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**<br>\n",
    "Because the two portfolios shared essentially the same set of risk factors, their returns were not independent, so a paired comparisons test is appropriate. Let $μ_d$ stand for the population mean value of difference between the returns on the two portfolios during this period. We test $$H_0: μ_d = 0$$ versus $$H_a: μ_d ≠ 0$$ at a 0.05 and 0.10 significance level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6488"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the sample mean difference\n",
    "diff = df['Difference'].sum()/len(df['Difference'])\n",
    "diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3204691085620541\n"
     ]
    }
   ],
   "source": [
    "# Calculate the standard error of the mean difference\n",
    "numerator = 0 \n",
    "for i in range(len(df)):\n",
    "    sqr_diff = df['Difference'][i] **2 \n",
    "    numerator += sqr_diff\n",
    "\n",
    "denominator = len(df) - 1\n",
    "samp_var = (numerator / denominator)**0.5 \n",
    "std_error = samp_var / len(df)**0.5\n",
    "\n",
    "print(std_error) # slightly different with 6.71/(24**0.5) = 1.369673"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.4913405363238835"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate test statistics with n-1 = 23 degree of freedom \n",
    "t_stat = (diff - 0) / (std_error)\n",
    "t_stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0686576104190406\n"
     ]
    }
   ],
   "source": [
    "# Calculate the rejection points at 0.05 significance level\n",
    "t_score = scipy.stats.t.ppf(1-0.025, 23)\n",
    "print(t_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fail to reject null hypothesis.\n"
     ]
    }
   ],
   "source": [
    "# Making statistical decision\n",
    "if t_stat < -t_score: \n",
    "    print('Reject null hypothesis.')\n",
    "elif t_stat > t_score: \n",
    "    print('Reject null hypothesis.')\n",
    "else: \n",
    "    print('Fail to reject null hypothesis.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7138715277470473"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the rejection points at 0.10 significance level\n",
    "scipy.stats.t.ppf(1-0.05, 23)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis Tests Concerning Variance and Correlation\n",
    "### Tests Concerning a Single Variance\n",
    "The hypotheses are formulated similarly for means and variances:\n",
    "- $H_0:σ^2=σ^2_0$ versus $H_a:σ^2≠σ^2_0$\n",
    "- $H_0:σ^2≤σ^2_0$ versus $H_a:σ^2>σ^2_0$\n",
    "- $H_0:σ^2≥σ^2_0$ versus $H_a:σ^2<σ^2_0$\n",
    "\n",
    "The `chi-square`, `χ2`, test statistic is used. It is an `asymmetrical` family of distributions defined with degrees of freedom. It does not take on negative values. The following rejection points are used for the three hypotheses. \n",
    "- $H_a:σ^2≠σ^2_0$ . Reject $H_0$ if the test statistic is greater than $χ^2_{\\frac{α}{2}}$ or less than $χ^2_{1−\\frac{α}{2}}$\n",
    "- $H_a:σ^2>σ^2_0$. Reject $H_0$ if the test statistic is greater than $χ^2_α$\n",
    "- $H_a:σ^2<σ^2_0$. Reject $H_0$ if the test statistic is less than $χ^2_{1−α}$\n",
    "\n",
    "If we have $n$ independent observations from a normally distributed population, the appropriate test statistic is\n",
    "$$χ^2=\\frac{(n−1)s^2}{σ^2_0}$$\n",
    "<div style=\"text-align: right\"> Equation 14</div>\n",
    "where $df=n − 1$,  $s^2=\\frac{\\sum_{i=1}^n{(X_i−\\bar{X})^2}}{n−1}$\n",
    "<div style=\"text-align: right\"> Equation 15</div>\n",
    "\n",
    "    In contrast to the t-test, for example, the chi-square test is sensitive to violations of its assumptions. If the sample is not actually random or if it does not come from a normally distributed population, inferences based on a chi-square test are likely to be faulty.\n",
    "    \n",
    "### Differences Between Variance\n",
    "#### Rejection Points for Hypothesis Tests on the Population Variance\n",
    "- “≠” $H_a$: Reject $H_{0}$ if the test statistic is > the upper $\\frac{α}{2}$ point (denoted $χ^2_α/2$) or < the lower $\\frac{α}{2}$ point (denoted $χ^2_1−α/2$) of the chi-square distribution with $df = n − 1$.\n",
    "- “≥” $H_a$: Reject $H_{0}$ if the test statistic is greater than the upper $α$ point of the chi-square distribution with $df = n − 1$.\n",
    "- “≤” $H_a$: Reject $H_{0}$ if the test statistic is less than the lower $α$ point of the chi-square distribution with $df = n − 1$.\n",
    "\n",
    "### Equality (Inequality) of Two Variances (F-Test)\n",
    "Like the chi-square distribution, the `F-distribution` is a family of `asymmetrical distributions` bounded from below by $0$. **Each F-distribution is defined by two values of degrees of freedom**, called `the numerator and denominator degrees of freedom`. (The `F-test`, like the chi-square test, is not robust to violations of its assumptions.)\n",
    "\n",
    "The hypotheses are formulated similarly for a single variance and for two variances.\n",
    "- $H_0:σ^2_1=σ^2_2$ versus $H_a:σ^2_1≠σ^2_2$\n",
    "- $H_0:σ^2_1≤σ^2_2$ versus $H_a:σ^2_1>σ^2_2$\n",
    "- $H_0:σ^2_1≥σ^2_2$ versus $H_a:σ^2_1<σ^2_2$\n",
    "\n",
    "A test concerning differences between the variances of the two populations is based on the ratio of sample variances\n",
    "$$F=\\frac{s^2_1}{s^2_2}$$\n",
    "<div style=\"text-align: right\"> Equation 16</div>\n",
    "\n",
    "with \n",
    "- $df_1 = n_1 − 1$ `numerator degrees of freedom` and \n",
    "- $df_2 = n_2 − 1$ `denominator degrees of freedom`. \n",
    "\n",
    "Note that $df_1$ and $df_2$ are the divisors used in calculating $s^2_1$ and $s^2_2$, respectively.\n",
    "#### Rejection Points for Hypothesis Tests on the Relative Values of Two Population Variances\n",
    "- $H_a:σ^2_1≠σ^2_2$. Reject $H_0$ if the test statistic is greater than upper $\\frac{α}{2}$ point of the F-distribution\n",
    "- $H_a:σ^2_1>σ^2_2$. Reject $H_0$ if the test statistic is greater than upper $α$ point of the F-distribution\n",
    "- $H_a:σ^2_1<σ^2_2$. Reject $H_0$ if the test statistic is greater than upper $α$ point of the F-distribution\n",
    "\n",
    "### Tests Concerning Correlation\n",
    "A `significance test of a correlation coefficient` allows us to assess whether the relationship between two random variables is the result of chance.\n",
    "- $H_0$: the correlation in the population is 0 ($\\rho = 0$); \n",
    "- $H_a$: the correlation in the population is different from 0 ($ρ ≠ 0$). \n",
    "\n",
    "As long as the two variables are distributed normally, we can test to determine whether the null hypothesis should be rejected using the sample correlation, $r$. The formula for the t-test is\n",
    "\n",
    "$$t=\\frac{r\\sqrt{n−2}}{\\sqrt{1−r^2}}$$\n",
    "<div style=\"text-align: right\"> Equation 17</div>\n",
    "\n",
    "This test statistic has a t-distribution with $n − 2$ degrees of freedom if the null hypothesis is true.\n",
    "\n",
    "The magnitude of $r$ needed to reject the null hypothesis $H_0: ρ = 0$ decreases as sample size $n$ increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**RECAP**\n",
    "\n",
    "- data-mining bias\n",
    "- sample selection bias\n",
    "- look-ahead bias\n",
    "- time-period bias\n",
    "- Common values for a standard normal distribution (85%, 90%, 95%, 97.5% and 99%)\n",
    "- What's a level of significance?\n",
    "- What's the power of test?\n",
    "- When to use Chi-square test? Formula?\n",
    "- When to use F-test? Formula?\n",
    "- Test concerning a single variance?\n",
    "- Rejection rules for all tests. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example 8, 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9893185569368186"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example 10\n",
    "scipy.stats.t.ppf(1-0.025, 82)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Issues: Non-Parametric Inference\n",
    "A test that is not concerned with a parameter, or that makes minimal assumptions about the population from which a sample comes.\n",
    "\n",
    "#### Situations for Using Non-Parametric Testing \n",
    "- when the data we use do not meet distributional assumptions \n",
    "- when the data are given in ranks \n",
    "- when the hypothesis we are addressing does not concern a parameter\n",
    "\n",
    "#### Nonparametric Tests Concerning Correlation: The Spearman Rank Correlation Coefficient"
   ]
  }
 ],
 "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": 4
}