nb/gauss-mix-kurtosis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Mixture and Leptokurtotic Assets\n",
    "\n",
    "Mixture of Gaussian distributions is an intuitive way to create \n",
    "distributions with skewness and kurtosis by adding \n",
    "two or more normal distributions.\n",
    "Such mixtures can synthesize *leptokurtosis* (\"fat tails\").\n",
    "\n",
    "Given the component variances of a GM(n) model it is easy,\n",
    "given the associated probabilities, to compute the moments\n",
    "of the resulting mixture. Our problem here is the *inverse*:\n",
    "from observable statistics, deduce the model parameters.\n",
    "This has been difficult due to solving tedious simultaneous\n",
    "equations for the so-called *method of moments problem*.\n",
    "However, with some simplications, the problem becomes\n",
    "tractable analytically and can provide insightful\n",
    "interpretations with concrete probabilities.\n",
    "\n",
    "We illustrate how a GM(2) model can synthesize Gaussian\n",
    "risk-equivalence for leptokurtotic financial assets.\n",
    "Both analytical and numerical solutions are provided.\n",
    "\n",
    "Appendix 1 confirms our analytical results by simulation.\n",
    "Also, we gain an experimental understanding of how kurtosis\n",
    "itself is distributed under small-sample conditions.\n",
    "\n",
    "Appendix 2 visualizes non-Gaussian distributions\n",
    "through quantile-quantile (Q-Q) probability plots.\n",
    "\n",
    "Appendix 3 points to Python package *scikit-learn*\n",
    "to fit multivariate Gaussian mixtures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "VIEWING: *if you encounter this message,*\n",
    "[Math Processing Error], *it means\n",
    "GitHub is timing out on rendering LaTeX equations.\n",
    "Please view this notebook at*\n",
    "[Jupyter](http://nbviewer.jupyter.org/github/rsvp/fecon235/blob/master/nb/gauss-mix-kurtosis.ipynb)\n",
    "instead. Issue [#3](https://github.com/rsvp/fecon235/issues/3)\n",
    "has been raised in this regard.\n",
    "\n",
    "*Dependencies:*\n",
    "\n",
    "- Repository: https://github.com/rsvp/fecon235\n",
    "- Python: matplotlib, pandas, sympy\n",
    "\n",
    "*CHANGE LOG*\n",
    "\n",
    "    2017-06-02  Visual demo of plotqq() for non-Gaussian distributions.\n",
    "                   Refactor gm2_* code to module ys_gauss_mix.\n",
    "    2017-05-08  Edit sympy code. Simulation confirms theory.\n",
    "    2017-05-05  Solve GM(2) for 2nd and 4th moments.\n",
    "    2017-05-04  Minor edits.\n",
    "    2017-05-03  First rough draft."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from fecon235.fecon235 import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ::  Python 2.7.13\n",
      " ::  IPython 5.1.0\n",
      " ::  jupyter_core 4.2.1\n",
      " ::  notebook 4.1.0\n",
      " ::  matplotlib 1.5.1\n",
      " ::  numpy 1.11.0\n",
      " ::  scipy 0.17.0\n",
      " ::  sympy 1.0\n",
      " ::  pandas 0.19.2\n",
      " ::  pandas_datareader 0.2.1\n",
      " ::  Repository: fecon235 v5.17.0221 develop\n",
      " ::  Timestamp: 2017-06-03T12:24:26Z\n",
      " ::  $pwd: /media/yaya/virt15h/virt/dbx/Dropbox/ipy/fecon235/nb\n"
     ]
    }
   ],
   "source": [
    "#  PREAMBLE-p6.15.1223d :: Settings and system details\n",
    "from __future__ import absolute_import, print_function, division\n",
    "system.specs()\n",
    "pwd = system.getpwd()   # present working directory as variable.\n",
    "print(\" ::  $pwd:\", pwd)\n",
    "#  If a module is modified, automatically reload it:\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "#       Use 0 to disable this feature.\n",
    "\n",
    "#  Notebook DISPLAY options:\n",
    "#      Represent pandas DataFrames as text; not HTML representation:\n",
    "import pandas as pd\n",
    "pd.set_option( 'display.notebook_repr_html', False )\n",
    "from IPython.display import HTML # useful for snippets\n",
    "#  e.g. HTML('<iframe src=http://en.mobile.wikipedia.org/?useformat=mobile width=700 height=350></iframe>')\n",
    "from IPython.display import Image \n",
    "#  e.g. Image(filename='holt-winters-equations.png', embed=True) # url= also works\n",
    "from IPython.display import YouTubeVideo\n",
    "#  e.g. YouTubeVideo('1j_HxD4iLn8', start='43', width=600, height=400)\n",
    "from IPython.core import page\n",
    "get_ipython().set_hook('show_in_pager', page.as_hook(page.display_page), 0)\n",
    "#  Or equivalently in config file: \"InteractiveShell.display_page = True\", \n",
    "#  which will display results in secondary notebook pager frame in a cell.\n",
    "\n",
    "#  Generate PLOTS inside notebook, \"inline\" generates static png:\n",
    "%matplotlib inline   \n",
    "#          \"notebook\" argument allows interactive zoom and resize."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian mixture in general and its four moments\n",
    "\n",
    "We build upon the cumulative distribution function (cdf) $\\Phi$ of the standard Gaussian $N(0, 1)$.\n",
    "\n",
    "For $0\\leq p_i \\leq 1$ such that $\\sum_{i=1}^n{p_i} = 1$, \n",
    "we define **Gaussian mixture GM(n)** by the following cdf \n",
    "and the associated probability density function (pdf):\n",
    "\n",
    "$$\\begin{aligned}\n",
    "F_n(x) &= \\sum_{i=1}^{n}{p_i}\\Phi (\\frac{x-\\mu_i}{\\sigma_i}) \\\\\n",
    "f_n(x) &= \\sum_{i=1}^{n}{p_i}\\phi (x;\\mu_i,\\sigma_i^2) \\\\\n",
    "\\phi (x;\\mu_i,\\sigma_i^2) &= \\frac{1}{\\sqrt{2\\pi}\\sigma_i}e^{-{(x-\\mu_i)^2}/{(2\\sigma_i^2)}}  \n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *intuitive idea* boils down to having n jars, each containing real numbers\n",
    "which are normally distributed, specified only by their mean and variance. \n",
    "We then randomly pick $x$ from the i-th jar with probability $p_i$.\n",
    "This is essentially how we will later simulate GM(n)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition 1:** The centralized moments of a Gaussian mixture GM(n) are:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mu &= \\sum_{i=1}^{n}{p_i}\\mu_i \\\\\n",
    "\\sigma^2 &= \\sum_{i=1}^{n}{p_i}(\\sigma_i^2+\\mu_i^2)-\\mu^2 \\\\\n",
    "skewness =: s &= \\frac{1}{\\sigma^3}\\sum_{i=1}^{n}{p_i}(\\mu_i-\\mu)[3\\sigma_i^2+(\\mu_i-\\mu)^2] \\\\\n",
    "kurtosis =: \\kappa &= \\frac{1}{\\sigma^4}\\sum_{i=1}^{n}{p_i}[3\\sigma_i^4+6(\\mu_i-\\mu)^2\\sigma_i^2+(\\mu_i-\\mu)^4]\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The history behind Proposition 1 goes back to Pearson (1894).\n",
    "The algebraic expressions used in solving the ***method of moments***\n",
    "in general are very tedious, see Cohen (1967) per Wang (2015, section 4.1).\n",
    "\n",
    "For the standard Gaussian distribution, $s=0$ and $\\kappa=3$.\n",
    "The importance of Proposition 1 is that mixture GM(n) can be a\n",
    "distribution with values for skewness and kurtosis which are non-standard.\n",
    "\n",
    "Returns on financial assets are said to be \"non-Gaussian.\"\n",
    "They are known to be \"fat-tailed\", or more technically,\n",
    "***leptokurtotic***: $\\kappa > 3$.\n",
    "This is where our mixture model can be helpful for risk management."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moments of the zero-mean GM(2)\n",
    "\n",
    "We focus our attention to a mixture built with two Gaussians,\n",
    "and examine its *even* moments, $\\sigma^2$ and $\\kappa$\n",
    "by setting the component means $\\mu_i$ to zero.\n",
    "We can later account for a Gaussian mixture GM(2) with non-zero mean \n",
    "since the first *odd* moment of a distribution is merely an additive shift,\n",
    "without any deformation to shape."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Corollary 1.1:** For $p+q=1$ and $\\mu_1=\\mu_2=0$, \n",
    "the moments of this Gaussian mixture, which we shall call **zero-mean GM(2)**,\n",
    "are given by:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mu &= 0 \\\\ \n",
    "\\sigma^2 &= p\\sigma_1^2 + q\\sigma_2^2 \\\\ \n",
    "skewness =: s &= 0 \\\\\n",
    "kurtosis =: \\kappa &= \\frac{3}{\\sigma^4}[p\\sigma_1^4 + q\\sigma_2^4]\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Corollary 1.2:** Combining the equations for even moments,\n",
    "the following is a necessary condition for the **zero-mean GM(2)**: \n",
    "\n",
    "$$\\rm K := \\frac{\\kappa}{3} = \\frac{p\\sigma_1^4 + q\\sigma_2^4}{(p\\sigma_1^2 + q\\sigma_2^2)^2}$$\n",
    "\n",
    "Note: Our defined $\\rm K$ should not be confused with\n",
    "Fischer \"*excess kurtosis*\" := $\\kappa-3$\n",
    "(which is sometimes erroneously reported as the kurtosis statistic).\n",
    "Our $\\kappa$ is also known as the *Pearson* kurtosis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structured zero-mean GM(2)\n",
    "\n",
    "Given data set {$x$}, we can compute its variance and kurtosis.\n",
    "Thus, $\\sigma$ and $\\rm K$ are observable.\n",
    "Corollary 1.2 says that GM(2) parameters can indeed synthesize $\\rm K$,\n",
    "but its messiness suggests that we impose further structure for clarity.\n",
    "\n",
    "A sensible requirement is a strict ordering,\n",
    "$\\sigma_1 < \\sigma < \\sigma_2$\n",
    "in consideration of the weighted equation for $\\sigma^2$.\n",
    "We can impose this ordering by constrained constants $a$ and $b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemma 2.1:** Let $\\sigma_1 = a\\sigma$ and $b\\sigma = \\sigma_2$\n",
    "such that $0<a<1<b$, then for zero-mean GM(2)\n",
    "the following probability satisfies the condition for *kurtosis*:\n",
    "\n",
    "$$ p = \\frac{\\rm K - b^4}{a^4 - b^4}$$\n",
    "\n",
    "*PROOF*. Substitute the two equalities for $\\sigma$ into the kurtosis equation\n",
    "of Corollary 1.1. Also note that by construction: $q=1-p$.    $\\square$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemma 2.2:** Let $\\sigma_1 = a\\sigma$ and $b\\sigma = \\sigma_2$\n",
    "such that $0<a<1<b$, then for zero-mean GM(2)\n",
    "the following probability satisfies the condition for *variance*:\n",
    "\n",
    "$$ p = \\frac{1 - b^2}{a^2 - b^2}$$\n",
    "\n",
    "*PROOF*. Substitute the two equalities for $\\sigma$ into the $\\sigma^2$ equation\n",
    "of Corollary 1.1. Also note that by construction: $q=1-p$.    $\\square$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition 2:** Let $\\sigma_1 = a\\sigma$ and $b\\sigma = \\sigma_2$\n",
    "such that $0<a<1<b$, then for zero-mean GM(2)\n",
    "the following satisfies conditions for both *kurtosis* and *variance*:\n",
    "\n",
    "$$\\frac{\\rm K - b^4}{a^4 - b^4} = \\frac{1 - b^2}{a^2 - b^2}$$\n",
    "\n",
    "*PROOF*. Immediate from Lemma 2.1 and 2.2.    $\\square$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic strategy and conclusion\n",
    "\n",
    "Proposition 2 shows how constants $a$ and $b$ are jointly constrained\n",
    "by the moments of the Gaussian mixture.\n",
    "\n",
    "Our **strategy** will be to choose $b>1$, then use the kurtosis $\\kappa$\n",
    "to compute $a$. Probability $p$ can then be computed using\n",
    "Lemma 2.2, for which $q=1-p$ follows.\n",
    "Since $\\sigma_1 = a\\sigma$ and $b\\sigma = \\sigma_2$\n",
    "the components of our zero-mean GM(2) are fully resolved,\n",
    "given specific $\\sigma$.\n",
    "This strategy will satisfy the necessary condition derived in Corollary 1.2.\n",
    "\n",
    "From two observable statistics $\\sigma$ and $\\kappa$,\n",
    "we have *deduced*, not fitted, the parameters of a zero-mean GM(2).\n",
    "An observable mean $\\mu$ which is non-zero can be added to\n",
    "our model to obtain a zero-skew GM(2).\n",
    "\n",
    "Fitting a general Gaussian mixture is computational very intensive,\n",
    "and for returns on financial assets, the results are often\n",
    "not robust (see Appendix 3). Our analytical method of moments\n",
    "allows us to constrain parameters in a way which has an optimal\n",
    "interpretation in Hilbert space, which for finance\n",
    "boils down to the Markowitz mean-variance framework.\n",
    "(An alternative model of leptokurtosis, the jump-diffusion process,\n",
    "does not have such an ideal mathematical interpretation.)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerically solving Proposition 2\n",
    "\n",
    "Given $b$ we could analytically solve for $a$ rewriting the solution for quadratic equations.\n",
    "Instead we shall demostrate the use of Python's sympy package (imported as *sym*) for symbolic mathematics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mgm2_strategy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkurtosis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0mgm2_strategy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkurtosis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m'''Use sympy to solve for \"a\" in Proposition 2 numerically.\u001b[0m\n",
       "\u001b[0;34m    >>> round( gm2_strategy(7, 2), 4 )\u001b[0m\n",
       "\u001b[0;34m    0.7454\u001b[0m\n",
       "\u001b[0;34m    '''\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  sym.init_printing(use_unicode=True)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  #        ^required if symbolic output is desired.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msymbols\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mK\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkurtosis\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m3.0\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  Use equation from Prop. 2\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mLHS\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mRHS\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0ma_solved\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolveset\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLHS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mRHS\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdomain\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mReals\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  ... expect negative and positve real solutions\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#      provided in sympy's FiniteSet format.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mif\u001b[0m \u001b[0ma_solved\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEmptySet\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#               ^when no feasible solution was found in domain.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#                Do not accept imaginary solutions :-)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0msystem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdie\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Argument b could not treat kurtosis; INCREASE b.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#     ^dies when kurtosis > 12 and b=2, for example.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#           SPX returns since 1957 have kurtosis around 31.6\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#           which is very high, requiring b>3.4 for feasiblity.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#  But a>0 by construction, so extract the positive real number:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0ma_positive\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma_solved\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0ma_positive\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Dropbox/ipy/fecon235/lib/ys_gauss_mix.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm2_strategy??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In writing the function `gm2_strategy`, we find that\n",
    "even a numerical solution (as opposed to a symbolic one)\n",
    "constrained to the real domain (as opposed to the complex) is messy.\n",
    "Surprisingly, $b$ is not entirely a *free* parameter\n",
    "since its particular value may be inadmissible when kurtosis is large.\n",
    "The function will work when $b$ is adjusted slightly higher\n",
    "when kurtosis is large enough to generate a fatal error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.745355992499930"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Example of deducing constant \"a\"\n",
    "#  using default value for \"b\":\n",
    "gm2_strategy(kurtosis=7, b=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " !!  FATAL 1: __main__.py Argument b could not treat kurtosis; INCREASE b.\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "1",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\u001b[0;31m:\u001b[0m 1\n"
     ]
    }
   ],
   "source": [
    "#  Example of deducing constant \"a\", INTENTIONALLY fatal:\n",
    "\n",
    "#  gm2_strategy(kurtosis=15, b=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.487950036474267"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Example of deducing constant \"a\"\n",
    "#  and re-adjustment of b given high kurtosis:\n",
    "gm2_strategy(kurtosis=15, b=2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical solution for GM(2)\n",
    "\n",
    "So for suitably chosen b, we can implement our strategy\n",
    "to completely specify our GM(2) model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mgm2_main\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkurtosis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0mgm2_main\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkurtosis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m'''Compute specs for GM(2) given observable statistics and b.'''\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgm2_strategy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkurtosis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  Probability p as given in Lemma 2.2:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#   The returned parameters can then be used in simulations, \u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#   e.g. simug_mix() in lib/yi_simulation.py module,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#   or for synthetic assets in risk management:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Dropbox/ipy/fecon235/lib/ys_gauss_mix.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gm2_main??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constants a, b: [0.745355992499930, 2]\n",
      "GM(2), p:  0.870967741935484\n",
      "GM(2), sigma1: 0.0968962790249909\n",
      "GM(2), q:  0.129032258064516\n",
      "GM(2), sigma2: 0.26\n"
     ]
    }
   ],
   "source": [
    "#  Demo of zero-mean GM(2) parameters using gm2_main():\n",
    "gm2_print(kurtosis=7, sigma=0.13, b=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical application in finance\n",
    "\n",
    "Suppose we are examining the returns of a financial asset XYZ\n",
    "and find evidence of a leptokurtotic distribution where $\\kappa=7$.\n",
    "Clearly this is non-Gaussian (though kurtosis estimates\n",
    "are not known to be stable due to their fourth power).\n",
    "The annualized volatility was observed to be 13%, \n",
    "so $\\sigma=0.13$ and we wish to use the GM(2) model\n",
    "to manage risk.\n",
    "\n",
    "We establish two synthetic assets on our books to represent XYZ\n",
    "by our choice of constrained constant $b$:\n",
    "\n",
    "- XYZ1 with volatility $a\\sigma$ where $a=0.7454$ is deduced,\n",
    "- XYZ2 with volatility $b\\sigma$ where $b=2.0000$\n",
    "\n",
    "The idea here is that XYZ1 models XYZ behaviour\n",
    "during normal trading (with probability $p=0.8710$),\n",
    "while XYZ2 models XYZ in an extremely volatile environment\n",
    "(with probability $q=0.1290$).\n",
    "\n",
    "Note that if kurtosis should increase, then upon recomputation, \n",
    "$q$ will increase because more \"fat-tail\" risk has emerged.\n",
    "\n",
    "We can still use our structured zero-mean GM(2) with\n",
    "*non-zero* mean returns because the mean is an *odd* moment\n",
    "merely providing an additive shift -- without changing\n",
    "the *shape* of the distribution. The pressing issue was \n",
    "actually getting suitable probability weights $p$ and $q$.\n",
    "\n",
    "If \\$100,000 of XYZ was on our books at 13% volatility, we use the\n",
    "probability weights to replace it with: \n",
    "- \\$87,097 of XYZ1 at 9.69% volatility, and\n",
    "- \\$12,903 of XYZ2 at 26% volatility.\n",
    "\n",
    "Important to emphasize that synthetic assets XYZ1 and XYZ2\n",
    "are Gaussian assets, so we are free to use the \n",
    "optimal linear mean-variance framework,\n",
    "yet the leptokurtotic risks are approximately covered.\n",
    "GM models are excellent extensions for\n",
    "VaR (Value at Risk) evaluations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "***Historical note:*** The specifics of the XYZ asset were actually *stylized*\n",
    "after SPX (S&P 500 index). The volatility from 2012 to 2017 was approximately\n",
    "12.5%, and during the Great Recession from 2007 to 2012 the volatility was 26.5%.\n",
    "So $b=2$ is not unreasonable since the volatility was roughly double. \n",
    "Round figure of 7 for kurtosis is reasonable, but when ten years is considered\n",
    "from 2007 to 2017 the overall kurtosis is about 13.5 since\n",
    "the extreme returns during the Great Recession is combined\n",
    "with compressed returns of the latter period.\n",
    "\n",
    "The analytics of leptokurtosis demonstrates that there is a\n",
    "small chance of double volatility lurking in financial returns every day."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "***Note on option pricing***: Models for option pricing usually assume\n",
    "logarithmic returns follow a Gaussian stochastic process, which by\n",
    "assumption sets Pearson kurtosis rigidly fixed to 3.\n",
    "By forecasting kurtosis during the life of an option,\n",
    "we can apply Gaussian mixture techniques to re-evaluate the standard model.\n",
    "The details will covered in a separate notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Appendix 1: Confirm analytical results using simulation\n",
    "\n",
    "We now use our simulation module `nb/yi_simulation.py`\n",
    "with generic GM(2) specifications to see if our \"observed\" statistics\n",
    "$\\kappa$ and $\\sigma$ can be reproduced. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0msimug_mix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msigma1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m256\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0msimug_mix\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0msigma1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m256\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m'''Simulate array from zero-mean Gaussian mixture GM(2).'''\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#     Mathematical details in nb/gauss-mix-kurtosis.ipynb\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  Pre-populate an array of shape (N,) with the FIRST Gaussian,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#  so that most work is done quickly and memory efficient...\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0marr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msimug\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0msigma1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mN\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;31m#     ... except for some random replacements:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#                p = 1-q = probability drawing from FIRST Gaussian.\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#  So with probability q, replace an element of arr\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;31m#  with a float from the SECOND Gaussian:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mmaybe\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mq\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0marr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrandog\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0msigma2\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0marr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Dropbox/ipy/fecon235/lib/yi_simulation.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  Here is the function for simulating zero-mean GM(2):\n",
    "simug_mix??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Y\n",
      "count  2560.000000\n",
      "mean      0.005134\n",
      "std       0.131963\n",
      "min      -0.674101\n",
      "25%      -0.068051\n",
      "50%       0.006167\n",
      "75%       0.079101\n",
      "max       0.668929\n",
      "kurtosis  6.00261\n"
     ]
    }
   ],
   "source": [
    "#  Use XYZ example to generate sample array of GM(2) mixture:\n",
    "mix = simug_mix(sigma1=0.096896, sigma2=0.26, q=0.12903, N=2560)\n",
    "\n",
    "#  Show observed stats:\n",
    "stat(mix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By repeatedly executing the previous block,\n",
    "we see that sigma (standard deviation) is about 0.13\n",
    "as expected from our analytics,\n",
    "but the kurtosis can fluctuate. In small samples,\n",
    "the fourth power can show up looking like outliers.\n",
    "The count=2560 corresponds to ten years of daily data.\n",
    "\n",
    "So *what is kurtosis on average here?*\n",
    "This is a typical question asked of small-sample distributions.\n",
    "Let's collect some data from repeatedly running a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Y\n",
      "count  1000.000000\n",
      "mean      6.944915\n",
      "std       0.708962\n",
      "min       5.102336\n",
      "25%       6.455984\n",
      "50%       6.877510\n",
      "75%       7.365816\n",
      "max       9.790339\n",
      "kurtosis  3.711716\n"
     ]
    }
   ],
   "source": [
    "runs = 1000\n",
    "kurt_samples = np.zeros(runs)\n",
    "\n",
    "for i in range(runs):\n",
    "    mix = simug_mix(sigma1=0.096896, sigma2=0.26, q=0.12903, N=2560)\n",
    "    kurt_samples[i] = kurtfun(mix)\n",
    "    #                 kurtfun yields Pearson kurtosis \"kappa\"\n",
    "\n",
    "stat(kurt_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***On average, kurtosis from our experiment is about 7 as expected\n",
    "from the mathematical propositions.***\n",
    "In theory, kurtosis should converge to 7 as we increase N,\n",
    "the number of points per run, and as we increase the number\n",
    "of simulated runs.\n",
    "\n",
    "It is very interesting how much kurtosis can fluctuate depending\n",
    "on the particular sample: its std is about 0.7, so we can expect\n",
    "to observe kurtosis roughly between 5.6 to 8.4 when\n",
    "the true kurtosis is 7."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Appendix 2: Visualizing a non-Gaussian distribution\n",
    "\n",
    "It is helpful to *visualize* the mathematical results above.\n",
    "The most useful plot for non-Gaussian distributions is\n",
    "the **Quantile-Quantile** [Q-Q probability plot](https://en.wikipedia.org/wiki/Q–Q_plot).\n",
    "\n",
    "The default setting for our *plotqq()* procedure will\n",
    "graphically compare the distribution of observed data\n",
    "against a theoretical Normal distribution.\n",
    "Details are covered in module `lib/yi_plot.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEZCAYAAABmTgnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX6wPHPI7ihJqZi7mQuaaaCLVqmlC1WaotlRd2y\nkiwrrW7dut2ude3eyn4tat3UwGu2mOXSZlYq7oZWIpp7Zrgj4i64wvP7Yw404gCDDgwzPO/Xi5dz\nznznnO/D4Dxzvs853yOqijHGGFOYCv7ugDHGmLLPkoUxxpgiWbIwxhhTJEsWxhhjimTJwhhjTJEs\nWRhjjCmSJQsTEETkPhFZ4O9+GFNeWbIox0TkDxG5qpT3WV9Etpzmy726KEhE5ojIA6e5D1NMIjJO\nRIaexusqiUiCiKSKyH4RSRaRHm7PNxWRHBE5ICIHnX//kW8b0SIyz3l+h4g87ouYzKlC/d0BU+7c\nAHzn706UBSISoqrZ/u6HH4UCm4ErVHWLiNwIfC4ibVV1s9NGgZrq4ephEamN629pMDAZqAw0Kp2u\nlz92ZFFOiciHQBPgG+cb29Nu3+T6ichmEdktIgNE5CIRWS4ie0TkHbdt3CciC0XkHRHZJyKrvThS\nuQGYXkCfckTkcRH5XUTSReT1Qvp/mYj8JCJ7RWSJiHR21v8buAJ414lrpLP+GhFZ47R/R0Tm5h59\niEgFEXlDRHaJyAYRGej0xeP/D+eI7K/O72SviHwqIpXcno8Tkd9EJENEvhSR+vliHCgi64H1buse\nEZH1zjfsoSLSTEQWOb/XiSLi1Rc7EakiIm8639b3ish8EansPNdbRFY67+NsETk/X7+auS3nHS2I\nSDcR2SIiT4nIThHZJiL9cmMF7gb+5vy+v/KmnwCqmqWqQ1V1i7P8LfAH0NE9JAr+nHoK+F5VJ6rq\nCVXNVNV13u7fFJOq2k85/cH1H/NKt+WmQA7wHlAJuBo4DEwFagMNgJ24vgkC3AccBwYBIUBfYB8Q\nXsD+QoFdQLUCns8BEoGauL4hrgMecNvXfOdxLWAPEIvrg+ROZ7mW8/yc3Nc5y7WBA8AtTj+fAI65\nbfthYLUTXzgwG8gGKhTye1sM1HParwYecp67yomxPVARGAnMyxfjD87rKrut+wKoBrQGjgAznfej\nBrAK+IuX7+l/nf6fg+uDtpPTj5bAIad/IcAzwG9AqPO6bKCZ23bGAUOdx92c9/lF57XXA5m4vvGf\n1Nbt9d8Ae533Jf+/XxfQ93rO31tLt7/HbGALriOQ/wG13donAsOBRbj+Lr8CGvv7/1Ww/tiRhZF8\ny4rrP/4xVZ2F60PhU1XdrarbgQVAlFv7nao6UlWzVfVzXB/wNxawr65AiqpmFtKf11R1v6puxfVB\ncJeHNjcC61V1gqrmqOpEYC3Qq4Bt3gCsVNUvnH4Ox/Xhkut2YLiqblfVfcCrhfQv1whV3em0/wbo\n4KyPBcaq6nJVPQ78HegsIk3cXvuKqu5T1aNu64ap65vxGmAlMENVN6nqQVxDLe6/c49ERID7gUGq\nmqYui51+9AWmqepsdQ19vQFUBS7LfXkRmz8GvOz8/r7DlXhaFdRYVXupai1VPdvDv7099D0U+BgY\np6rrndUZwMW4kkZHXInzE7eXNQLuBR4HGgOpwKdFxGFOkyUL40m62+PDnPzBehio7ra8Ld9rN+H6\nhu5JgUNQbrZ6sa0GznP599uwgG02wPXt1N2WQp7Pv21P3H8nWfz5Ozmpb05i3J2vb+4x5irO77wg\ndXCN22/08Fz+fimumAv6neW3W1Vz3JbdYz4jTpL7GDiK64M/t4+ZqprsfCHYBTwGXCsi1Zwmh4Ev\nnDbHgH8Bl4lIDV/0y5zMkkX55osph/N/2DQBthfQ1ptk0diLbW0HIj3sNzdx5Y9rh/N8QfvZkW+5\naRF9LMx299c7H2y1OTlBlNRUzxm4hrDOK6pfjsZu/coCwtyeO6cY+/VUfJ7udgZT/p9v8zUfiyvR\n3apFF/yVPz+3VnjYt02jXUIsWZRvaUCzfOuKGo7IL8IpSoeKyO3A+XhICCISCVTSoguQz4hIuIg0\nxnWWy0QPbaYDLUTkThEJEZE7cI31T3Oe38nJcX0LtBGRm532g3GNj+f6HBgkIg1FpBbwbFFBF+JT\n4H4RaecUll8BFqtTxPUFpxjdNf9652jhf8Bb4jpFuYKIdBKRirhivFFErnTeq6dxJZYk5+XLgFjn\nNT1w1Sm8lf/3jareoKo1VPUsDz95w5QiMhrX30xv5+jAPc5LRKSluNQGRgBznKE5cNVKbnF+1xWB\nfwIL3Z43PmTJonx7Dfinc3bMU866or6p5V9eArTA9a32ZaCPqu71sK8bKfqoAlxFyqVAMq5awP/y\nN1DVPUBP4Glnv08DNzrrwfWhcru4zuYarqq7cdUlhjntz8NVFM0Vj6vovBz4BZhSRB8L/Paqqom4\nPrSm4jrSORdXAb6w13r97dhJogeAXwto8rTz3M+4hr9ew1WoXw/cA7yLqwB/I9BLVU84r3sC6I2r\nCH0XroJ7Ydz7OBa4wPk7mlrE69xjaQI8hKves9PtSCS3TtUM+N6JdwWu5Bab1wHVOcDzuP6ucr/4\nxGJKhLi+jPixA65vMcNxJa6xqjrMQ5sY4G1cZ3XsUtUrS7WTxiMRuQ94UFVP+Zbroe23wDuq+n0h\nbXKA5qrqaczdp0RkDvCRqp6SjESkKa5x/4r5xun9TkTuBtqo6j+KbGyMD/n1ojznPPZ3ge64xlR/\nFpGvVHWtW5uauE4HvFZVt4lIHf/01pyhOc5PoCjucFypUNVPim5ljO/5exjqEuA35xTB47jGp2/K\n1yYWmKKq2wBUNaOU+2h8QFXfyHeqqMdmpdIZ7/ZlhVJj3Ph7uo+GnHzK4lZcCcRdS6CiM2xQHRip\nqh+VUv9MIVR1PDDeh9sL8dW2vNhXgVeaq+omXBefGWMc/k4W3ggFonFdeVoNSBKRJFXd4N9uGWNM\n+eHvZLGNk89/b8SpF3ltBTJU9QhwRETm45pK4ZRkISI2dGCMMcWkqkXW6Pxds/gZaC6uCewq4TrF\n8Ot8bb4Cujjnx4cBlwJrCtqgr+ZBKWs/L774ot/7YPFZfBZf8P14y69HFqqaLSKPATP489TZNSIy\nwPW0vq+qa0XkB1znWWcD76vqaj922y9SU1P93YUSZfEFNosv+Pl7GAp1nXffKt+6MfmW38A18Zkx\nxhg/8PcwlPFSv379/N2FEmXxBTaLL/j5/QpuXxIRDaZ4jDGmpIkIGgAFbuOluXPn+rsLJcriC2wW\nX/CzZGGMMaZINgxljDHlmA1DGWOM8RlLFgEi2MdMLb7AZvEFP0sWxhhjimQ1C2OMCRJJiVnMnbKb\nmD616dw9rOgXYDULY4wpV5ISs1gcO5LBo1qxOHYkSYlZPt2+JYsAEexjphZfYLP4/G/ulN0MSB9K\nGIcZkD6UeVN9e584SxbGGBMEYvrUZkzEELKoypiIIXS71bd3oLaahTHGBImkxCzmTc2g2611fF6z\nsGRhjDHlmBW4g0wgjJmeCYsvsFl8wc+ShTHGmCLZMJQxxpRjNgxljDHGZyxZBIhgHzO1+AKbxRf8\nLFkYY4wpktUsjDGmHLOahTHGGJ+xZBEggn3M1OILbBZf8PN7shCRHiKyVkTWi8izhbS7WESOi8it\npdk/Y4wxfq5ZiEgFYD3QHdgO/AzcqaprPbSbCRwG/qeqUwvYntUsjDGmGAKlZnEJ8JuqblLV48BE\n4CYP7R4HJgPppdk5Y4wxLv5OFg2BLW7LW511eUSkAXCzqo4Cisx+wSrYx0wtvsBm8QW/UH93wAvD\nAfdaRqEJo1+/fkRGRgIQHh5Ohw4diImJAf58w23Zlm3Zlsvrcu7j1NRUisPfNYtOwEuq2sNZfg5Q\nVR3m1mZj7kOgDpAJPKSqX3vYntUsjDGmGALifhYiEgKsw1Xg3gH8BNylqmsKaD8O+MYK3MYY4xsB\nUeBW1WzgMWAGsAqYqKprRGSAiDzk6SWl2sEyxP0QMhhZfIHN4gt+fq9ZqOr3QKt868YU0PaBUumU\nMcaYk9jcUMYYE0RUFRHvTxz1dhjK70cWxhhjzoyqsmjLIuKT4xGED27+wOf78Pd1FsZLwT5mavEF\nNovvzCUlZvHqwC0kJWZ5/Zpdmbt488c3afNeG/p/3Z92Ee34v2v+r0T6Z0cWxhjjZ0mJWSyOHcng\n9KGMmTIEJgyic/cwj21zNIfZf8wmPjmeHzb8wE3n38T7Pd+nS5MuxRp+Ki6rWRhjTAlLSsxi7pTd\nxPSp7TEJvDpwC4NHtSKMw2RRlZED1/Lcf5uc1Gb7we2MWzaOscvGclbls4iLjuPudncTXiX8jPoW\nEKfOGmNMsMs7ahjVisWxIz0OM8X0qc2YiCFkUZUxEUPodmsdAE7knGDa+mncNPEm2r7Xls37N/PZ\nbZ+xbMAyHr3k0TNOFMVhySJA2JhwYLP4AtuZxDd3ym4GpA8ljMMMSB/KvKkZp7Tp3D2MThMGMXLg\nWjpNGET9jukMmTOEyOGR/Hv+v+ndsjebn9zMmF5juLjhxSU63FQQSxbGGFOCCjpqyK9jTCjNH/2J\nl3bcwkXvX8T+I/v57u7vWNx/MQ9GP0j1StVLuecns5qFMcaUsKTELOZNzaDbrXVOqVms372ehOQE\nxi8fT+s6rYmLjuPW1rdStWLVUulbQMwN5WuWLIwxgeDw8cNMWTOF+OR41mWs47729/Fg9IO0rN2y\n1PtiBe4gY2PCgc3iC2y+iu/Xnb8y6LtBNH67MR+v+JhBlwxi85ObGXbNML8kiuKw6yyMMaYEHTp2\niIkrJxKfHM/2g9t5oMMD/PLQL0SGR/q7a8Viw1DGGONjqsov238hPjmeSasn0a1pN+Ki4+jRvAch\nFUL83b2T2NxQxhhTyvYd2cfHKz4mITmBA0cP0D+6P6sGrqJBjQb+7toZs5pFgLAx4cBm8QW2wuJT\nVRZsWsC9X9xL5PBIFm5eyJvXvsmGQRt4/orngyJRgB1ZGGPMadmVuYvxy8eTkJyAiBAXHcdb171F\nnTDP11EEOqtZGGOMl3I0h8SNicQnxzPj9xncfP7NxEXHcVnjy/xyVbUv2HUWxhjjI9sObGNcimsS\nv/Aq4cRFxxF7YWypzs1UUuw6iyBTnseEg4HFF3hO5Jzgm3Xf0PvT3pz/9PlsPbCVSbdPIvmhZAZe\nPDAoEkVxWM3CGGPc/LH3D8YuG8u4lHE0qdmEuOg4HqnzCNdfc32hrytqGvJAZ8NQxphy71j2Mb5a\n+xXxyfEk70jmnnb30D+6P20j2npsnz8x5E5DPiB9KGMihtCpkJsXlTVWszDGmCKsy1hHfHI8H634\niDZ12+RN4lcltEqBr/GUGOZO2V3kzYvKKqtZBJlgHBN2Z/EFtkCK7/Dxw3y0/CO6jutKtw+6EVoh\nlIX3L2TOfXOIvTDWY6Jwj8/T/Sm8nYY8kPm9ZiEiPYDhuBLXWFUdlu/5WOBZZ/Eg8Iiq/lq6vTTG\nBLrlactJSE5gwsoJXNLwEp7o9AS9WvaiYkjFYm0npk9txkwZkndkkTft+IRBjJwa63Ea8mDg12Eo\nEakArAe6A9uBn4E7VXWtW5tOwBpV3e8klpdUtVMB27NhKGNMnoNHD+ZN4rfj0A4ejHqQB6IeoElN\n74eIPBWuC7s/RaAJiJqFkwheVNXrneXnAM1/dOHWPhz4VVUbF/C8JQtjyjlV5eftPxO/NJ7JayYT\nExlDXHQc1513XYGT+LknBOCkx4FauPZWoNQsGgJb3Ja3OusK0h/4rkR7VEYF0pjw6bD4AltZiG/v\n4b28s+Qd2o9uT+yUWM47+zxWD1zNF3d8wQ0tbig0USyOHcngUa2Y3OdT5t82gsGjWrE4diQTR6Qx\nIH0oPxVy/+zywu81C2+JyJXA/UCXwtr169ePyMhIAMLDw+nQoQMxMTHAn3/QtmzLthwcy6pKSLMQ\n4pPj+eK7L7ik4SUMv3c4MZExzJ83n3VL11E/pn5e+1VLj3Dg9wuI6VOb5JT5pMw/QJhexKtOQti7\nfy4vM4UwDtM6/UXm745kTMQQWqe/yDPh99Cm2WqgSZmJ/3SWcx+npqZSHGVhGOolVe3hLHschhKR\ndsAUoIeq/l7I9mwYyphyID0znfEp40lYlkCIhBAXHcdf2v+l0En83E95/WfNd4iQdB7f93KBj3OH\nnYCgqU94Eig1ixBgHa4C9w7gJ+AuVV3j1qYJkAj8RVUXF7E9SxbGBKkczWHWxlnEJ8cz8/eZ3NL6\nFuKi4+jcqHORk/glJWbx1uBUxq+6iDAO8zIv8FfezLsu4tneq2nYqELeKa/BnBzy8zZZoKp+/QF6\n4EoYvwHPOesGAA85j+OB3UAysAz4qZBtabCaM2eOv7tQoiy+wFaS8W3dv1WHzh2qTd9uqlGjo/S9\nn97TfYf3ef36H2dl6lsRr+psuunL8oJmUlWfqpmgr4W/oplU1bciXtUfZ2UWuo1gfv+cz80iP6v9\nXrNQ1e+BVvnWjXF7HAfElXa/jDH+cyLnBNN/m058cjyLNi/ijgvuYOodU4muH+31NnLPcErbepxX\nnYvoULiv7S88NTwSgJFT7y43RxBnyqb7MMaUGRv3bmRssmsSv8jwSOKi4+h7QV+qVapW5Gtzk0Pd\nZtVImX+Axgs+LbAOYcnhTwFRs/A1SxbGBJ6jJ47y1TrXJH4paSncc6FrEr8LIi7w6vVJiVl8MiKd\nxgs+5ZJ9P7BIrkCVAmsSlihOZskiyMydOzfvFLhgZPEFttOJb23GWuKXuibxaxvRlrjoOG5pfUuR\nk/h9MiKdEHK4sGt43hHEsX1Z/JU3GcFgBjOC5bQnic48zGifHE0E8/vnbbLwe83CGFN+ZB3PYvLq\nycQnx7Nhzwb6te/Hjw/+SPOzmxf4mtwEsXdXNrV/nUOzzNVEkcyiaVdQT+Fx3mQ57RnNw3Qiibfk\nrzylbzK5Zn+e77qSOwafY0cTPmBHFsaYEpeSlkL80ngmrprIpQ0vJS46jp4te+ZN4ueeEOrUlbyj\nhtwEEZm5mkxcdYuCjiD+WfMdsrteSduu4WT8cciGnLxkw1DGGL86ePQgn678lPjkeHYe2pk3iV/j\nmq6p3fIfMYRm7qcR24gimUncTnN+PylBLKc9k7gtr80iuYKn9M28JGFHEKenRJKFiNQCGqvqijPp\nXEkJ5mQRzGOmYPEFutz4VJUl25aQkJzAlDVTuDLySuKi47j2vGtPmpsp4Y09JD8/hebH1+QlBPjz\nqOEEoackiIcZzXPVRrL7wiupW1dK9QgimN8/n9UsRGQu0NtpuxRIF5FFqvrUGffSGBMUDhw9wMgl\nI4lPjufIiSP0j+rPmkfX8MeSs5j79m62NdufN6y0f79Sad2vROVs5mFG5yUEIK/uMInbGc3DPMxo\nPq32IEsvvJnn6z7GXaccPZztn4DLoSKPLERkmapGiUh/XEcVL4rIClVtVzpd9F4wH1kYU9aoKvM2\nzSMhOYFp66dxQ4sbaLDhDn4ddR6VwiqhVODi3z+j67GZTOJ2qnAEgG005GFGM42eNGT7SUcMQN5R\nw/IFBwjRHBteKmE+G4YSkV+Ba4HxwD9U9WdLFsaUXzsP7WT88vEkJCdQKaQScdFxtNp3KyOfOkbT\nNd9T1S0pjOOBvGGlXF1YQApRdGAZr8nz1OsUySMvN7KE4Ce+vJ/FUOAH4HcnUTTDNY+TKUXu0wsH\nI4uvbMvOyeaHDT9w2+e3cf5/z2dtxlqebjSGbt9/w8wHbuS9G0dweM0mGpBGLfZRi32uaxwYQCeS\n2ME57CWcvYSzjGg6sIxRlZ7g9tcv4sMfW5b5RBHo758vFFmzUNVJwCS35Y1An5LslDGmbNh6YCv/\nmjKGz34bR8UjtWjwW18uWPY3UkNqUmnbfMKyk6lGQ6qxkx58zkraUoWjACwjmiiSGVXpCVr268T+\ntCPsSc9mJ7C5rvCkDS8FFG+GoVoCo4B6qtrWubdEb1X9d2l0sDhsGMqYM3ci5wTfrv+W12eM5ue0\nJFqs6EDU0mgkrT01OEgVjrKNhrRlFeAaVsqtP1Qhiw+5jxpNaxFRvyJ164rVHMo4X9Ys5gHPAGNU\nNcpZt1JV2/qkpz5kycKY07dx70b+NXUUn2/4kNDdDWn9cyfqr+rMRcfXA5BCezqwHPgzQSiSd93D\nsJDnibjE6g+Bxpc1izBV/SnfuhOn1y1zuoJ9zNTi8495M/dy3X3vUueRy2n12kUkTVrDTWP7c0P8\n32icchVPHB/LPmqyl3DCyGQH9fLqDj2ZxuqQC5nbegD/7fwQL/7QJSDqD6ejrL5/pcmbuaEyROQ8\nQAFE5DZcd7UzxgSol19PYkzSaHa3mEo4jYlaeCk11zyMZlflApbThdGMox+LuJyeTOM1niWn6bns\nDQslc/cxNtQOYXPzirzoDDHNnbuDzjHBlyTMn7wZhmoGvA9cBuwF/gDuUdXUEu9dMdkwlDGeJSVm\n8cHITSw5/g0bGk5Ea/5O/ZRr6ZUcwea9V9KB5XRhAePoR3UO0ZyNRJHMiwylUr1w7ny6Ef2ftgvg\ngpHPp/sQkWpABVU9eKadKymWLIz5U1JiFiNeSGNVxjJyWr3Dhra/cM6WptRIvoPb1h9hZ05DqnKE\nXdShBgdPShDZ4WdzXutKVn8oB3xZ4B7iab2qDj3NvpWYYE4WwTw3DVh8vpA7MV/qhhNs2bGHkPPi\n2drxBw6FHaHTstZUXnYXlx3YRm12nXSK62Eqs4yOhITXOO0EYe9f4PLl/Swy3R5XAXoCa063Y8YY\n30l4Yw8f/l8aWZkQnbWAzQ23sq7jArbd9DN1N15Ij9ldqPp7Z6rpMXZRlR3UoxqZtGVl3imuF7Sr\nxNt2eqspQrGnKBeRysAPqhpTIj06A8F8ZGFMroQ39jD53TR27Quh+f6lVKu6hSXtNpDR8SuyQkKJ\nTO7OSym/81XmQ1TnkE+PIEzwKbH7WTjTlP+sqgXf2spPLFmYYJaUmMVzD2ynzebvOUh1fo3cw56O\nU8lo8QuN1l9Ih+SLqJx6KcepxIWs4nIWMYSXOF6tNhFNKtGseUW7QM6cwtcTCeY2CgHqAkNV9d0z\n7qWPBXOyCOYxU7D4CpPwxh7i/51G4/0r2VYthBod3uPH6NXUPHGMiOQbiVpxAeGHXSPKuUcP2dVq\nUrP6iVI7i8nev8Dly5pFT7fHJ4Cdquqzi/JEpAcwHNcFgmNVdZiHNiOB63HVT/qpaoqv9m9MWeJe\ng6hUUTmqoTQ78DOh560m6brZ7D13GU3XRNPli/5cvjWb2uwhgauoUKUKDc51HT1Y/cGUhAKPLESk\n0K8jqrrnjHcuUgFYD3QHtgM/A3eq6lq3NtcDj6nqjSJyKTBCVTsVsL2gPbIwwSn39Na0jUc4eNRV\ngwCozR7SzjpGctQK0qO+pUnmUcKTb+XqlXVocPQQH3IvOTVrERGezc2PNbBrIMxp88WRxVJcw0+e\nNqJAs9Psm7tLgN9UdROAiEwEbgLWurW5CfgQQFWXiEhNEamnqjt9sH9j/CIpMYuXH99B2JqlRJBO\nRWqwmSYcq1CBtJbJLIv+gk2N0ui0sgktJz5J17RMpwbRjwrhtej/j8aWIEypKjBZqOq5pbD/hsAW\nt+WtuBJIYW22OevKVbII5jFTKB/xbfilHZPfTWN/VggNdy1HqY5Sgfrs5PdaUCX6Kb7vkEr1PfVo\nl9yJjpOu4ILjmzibzSTwINOr9+WBF8tmkigP718wx+cNb2oWuWdAtcB1nQUAqjq/pDp1Jvr160dk\nZCQA4eHhdOjQIe9Nzp0MzJZtubSWv/3sAEumNmdzRirtcr4gm6qE0xWlAiEVprO+yWaWdF1HRr0d\nNJ/TkNvjO7DywOtkVAxjZ6WfWFerFQ2qRjPwsQY0v2gFrhl3yk58thx4y7mPU1NTKQ5vzobqDwwG\nGgEpQCcgSVWvKtaePG+7E/CSqvZwlp8D1L3ILSKjgTmq+pmzvBbo5mkYymoWxl9yaw+pa45w4sSf\nxencGkQa9bma2aTQntC6KayI/pnUdj8SkVafBsnXcWxtX6rWrkt41RNWgzClytenzl4MLFbVDiJy\nPvCKqt7qg06GAOtwFbh3AD8Bd6nqGrc2NwCPOgXuTsBwK3CbssA9QTTev5IKzsz9tdnDQacGUYfd\nAHSp+B3fXHCUldFLyKyVTttlF1NnWU/S93amUZua/G1kEzuDyfiFL0+dPaKqR0QEEamsqmtFpJUP\n+oiqZovIY8AM/jx1do2IDHA9re+r6nQRuUFENuA6dfZ+X+w70MwN8jHTQIjPfe6l7ZuO0SxrJRGk\nc4im5FCBHCoBUJ+d7OQc7uAzvqzfmI3Rc5h21nxaSn3CFw2k6c4bCKsUQu3zQ3khSK6iDoT370wE\ne3ze8CZZbBWRcOBLYKaI7AU2+aoDqvo90CrfujH5lh/z1f6MKa7ci+La7V/IESpTjYrUIzyvOK0I\nh6jOMSoCkFr5LLZe+A3/iZ7BkbBDVE6+m3rzr6de2DU21bcJWMWa7kNEugE1ge9V9ViJ9eo02TCU\n8bUXBqTx2/tzOEBNLiOJFNoDEMFODlE9757UirKi8SGWRSdzoPVcam7pROONfWm442JuebRsnsFk\nDPhgGEpEpgMTgC9V9RCAqs7zXReNKduSErNYmLCe2lQkkj/ybi16jIoISl0y2Fi1Osntf2fPpd9Q\nsWo2d5/fnyG3JRBRLcLf3TfGpwq7B/cY4EbgDxH5XERuEZFKpdQvk4/7aW/BqKzFl5SYxfOxG+mb\nM4FQjiHO9GjVOUiqNCWxdRgf3ZXI9KeeofVfdvDtE2PZ88pG3r3vHx4TRVmLz9csvuBX2EV5XwFf\niUgY0Au4FxglIt8BE1R1Zin10ZgS5X5WUwVRtGIobXfNowoNOUIY3ZjPWO7n6FmZVOo5nfR2b3P2\nWdUYEh3H3e0+5uyqNsRkgl9xaxbtgPFAO1UNKbFenSarWZjicJ/NNYQTeae8ZlCXy0iiCwv4Wm4g\no3kyiy5aR0abVPq260NcdByXNLwEkSKHeY0p83x26qyI1AP6AncC9YHPgX5n2kFj/CF3Vtf9++H8\noyk0oiIkkI8JAAAfmklEQVRKBU5QKe+U10j+YFNNZWHUUVKihnD8UGN6nvso/33qQWpUruHvEIzx\niwJrFiISJyKzgWRcU308o6rNVPU5VV1eaj00QPCPmZZGfC8MSGPmMz/QNn02Zx9N4wSVeIz/Uo1D\nhHKMrRXqkNF6HtPvfo9PBrzNtrCKtPjsLV5vuYQP//7EGSUKe/8CW7DH543Cjiw6A68CiaqaU0r9\nMcYn3C+gO7T7GEc1hIq70qhLReqzk7asZB5dWcTlXHv2KP4efR5pHWYStq8JDX+/j4tnX01kZHUe\n+TA4Lpoz5kwV+7aqZZnVLMo390J1+wMLydLK1OAgB6lBFmF0J5F5dKU2e2gaup61rVcxOXo/R+tt\n5MqzYhkx4DHOr3O+v8MwplT5croPY8ocT0cOjXYtp64z/UYjtpBCe1qygZ2cw2O8SwpRtKj7Je93\nDGFvux+om9WG57o8wdN9bqdyaGV/h2RMmWbJIkAE+9w0RcWX8MYeJr+bRsWwUA5lCi03zyGLyjRx\nO3LIcZt+I/cCuh3Uo1LFPbxzQXXWdfwfG2uG0FHu5aPnlnJurdK4ZYt38QU6iy/4FXYFd4nfVtWY\nouTeUa7pmu9pwxF2UYdKhHk8chhHP3ZQj6ocAeBw/VXM6Pgbe9rOJjy9I5Hr/s3w3ncw4Bm7utqY\n4irsHtx/8OdtVZvguuuKAOHA5lK6k16xWM0ieOTWH8J/TmRd9nlcxRwAUmjPY7zLNHqSTt28mkQr\n1nM5i3i+8tOkXbSYPRd9AWEHuKPFA7zY5yEantXQzxEZUzb58n4W8cAXqjrdWb4euFlVB/ikpz5k\nySKwudchGq2bw96cGrRlFbXZxUraUoWj7KJOXmJ4jWfJjDiXULI52GQF6a0msrvFXHq2uY7+0f25\nutnVVJDCZrQxxvj05keqemFR68qCYE4WwThmmpsc9u7KZuP2xbTflkVY9gHW05JL+YkuLGAaPWnI\ndqqQxYfcR42mtagcFgpZR+j+aGXoMo2EZQmcyDlBXHQc97a/t0xO4heM7587iy9w+fJsqO0i8gLw\nsbN8N7D9TDpnyif35LB/v9Jo/RzCsg9QD9jOfho6k/U1ZjN7CWcZ0fRkGq/yPOd0juTtlxtx6VVV\nmPPHHOKT4/nXhu/pldaLUTeO4oomV9j0G8aUIG+OLM4GXgS64qphzAeGlsUCdzAfWQSigpIDwDYa\n0pZVeW1zjyKUPz/wD1OZzaHNueXVS7hxwFE+SPmAhGUJVK9UnbjoOO6+8G5qVa1V6nEZE0x8Ngzl\ntsFqqpp5xj0rQZYsyo6EN/aQ/PwUqhz3Ljk0YhtRJDMs5HkqtWxKjZqh1K6bQ5O//MrC7A+Zt2ke\nt7W+jbiOcVzc4GI7ijDGR7xNFkVW/0TkMhFZDaxxltuLyHs+6KMphkCZmyYpMYs7O29k+rPzqH98\nM7XYRy328TCj2UdN9hJ+0hDT6pAL2dC6J9Nbd+SrXv/jxR+68M6PVWj+yod8ceVVTNo/jJ4te7L5\nic3E944P2NleA+X9O10WX/DzpmbxNnAd8DWAqi4Xka4l2isTkJISs5jc51N0fzUeZ/RJRw65ySH3\nyGFnzVA21xVeHHwOnbuHMStxC/vrL+NfyfH8/P7PxLaN5dvYb2lXr52fozLGgHc1iyWqeqmILFPV\nKGfdclVtXyo9LAYbhvKvwb03UuebcXRhASlE0YFlvCbPU/l817BS3brCHU5yyPXb7t9ISE5g/PLx\ntKrTirjoOPq07kPVilX9GIkx5Ycvz4baIiKXASoiFYHBOENSpvxJSsxi7pTd1G1WjZT5B9i7K5s6\ndYUa9cM4/n1i3hBTFMmMqvQEt//nIvo/ffJkAEdOHGHK6inEJ8ezJmMN97a7l3n95tGqTis/RWWM\nKYo3RxZ1gBHA1biu4J4BDFbV3SXfveIJ5iOLsnCed8Ibe0j75yguPzKTSdxOFY7kFabf5TE+4l6W\n054PuZtNkVfyz4RzTzqKWJm+kvil8Xzy6yd0bNCRuOg4erfqTaWQSmUivpJk8QW2YI7PJ0cWIhIC\n/EVV7/ZZz/7cdi3gM6ApkAr0VdX9+do0Aj4E6gE5QLyqjvR1X0zRkhKzmPXPBfzvyH8YwWDqkwbA\nw4xmBIN5jHcZwwAGMIaGVXpyr5MoDh07xGcrPyNhWQJb9m/h/g7388tDvxAZHunfgIwxxeLNkcXP\nqnqxz3csMgzYraqvi8izQC1VfS5fm3OAc1Q1RUSqA0uBm1R1bQHbDNojC39KSszircGpDFw1MK8W\nkf/IYpFcweW6gPerPsFV/7qCqLv/IH5pPJNWT+KKplcQFx1Hj+Y9CK1gEx0bU5b4crqPt4GKuI4C\n8q6zUNXkM+zgWqCbqu50ksJcVS30zjMi8iXwjqomFvC8JQsfcx96ck8I5/2lE/vTjrAnPZu6dYW2\nXcPZkrqNzE6JJGaOZ/+R/TwY9SD3R91PgxoN/B2GMaYA3iYLVLXQH2COh5/ZRb3Oi+3uKWzZQ/tI\nXMNV1Qtpo8Fqzpw5pb7PH2dl6h1VvtRMqqqCzqab3tZ2lf44KzOvTU5Oji7YtEDv/eJerflqTe07\nqa/O/H2mZudkF2tf/oivNFl8gS2Y43M+N4v8zC5yTEBVryx2qnKIyExc9Ya8VbimDHnB064K2U51\nYDKuwvqhwvbZr18/IiMjAQgPD6dDhw55hancC2ts2bvl/42cxqVHXsmrRXxc8Txi+m2mc/c2ZGRl\nMGTcEL5d/y1VW1QlLjqOWyrfQnjVcGKalY3+27It2/Kpy7mPU1NTKQ5vhqHqAa8ADVT1ehFpA3RW\n1bHF2tOp210DxOifw1BzVLW1h3ahwDTgO1UdUcQ2tah4TNHc53Rqt/ZzLtn3A+9XfYIr/9WFZn1S\niE+O54cNP3DT+TfRP6o/XZp0Ccirqo0xvq1ZfAeMA/6hqu2dD+9leoZTlDsF7j2qOqygArfT7kMg\nQ1Wf8mKblizOUO5V2A33r+JhRvNctZFkXXcBob1mMPPAh9SoVMM1iV+7uwmvEu7v7hpjzpDP5oYC\n6qjq57hOXUVVTwDZZ9g/gGHANSKyDugOvAYgIvVFZJrz+HJcU6JfJSLLRCRZRHr4YN8Bx/0QsqTk\nnvUUvj+V/hVGMbvlYTb2GsAnbXuQU3sbE/tMZNmAZTx6yaM+TxSlEZ8/WXyBLdjj84Y35zFmikht\nnJqCiHQC9hf+kqKpa4rzqz2s3wH0dB4vAkLOdF+maEmJWSyOHcmtx77kzStr8WZUKC0PCE3W3My0\n+0fR/dqyd0MhY0zp8WYYKhp4B2gLrATqArep6oqS717x2DBU8eRO3VHr3IqMn/Up1SKeJqVBDl1/\nbcjOHe9yceN2p8zlZIwJLj6bG0pVk0WkG9AK19lM61T1uA/6aPwoKTGLrx/5J39Efs60antp3iCU\nRsm9WDfxBz48+zE6TbjWkoQxJk+BNQsRuTX3B+iNK1m0BHo560wp8tWY6eHjh/l4xcfEzr6KsTe9\nxXYasmjcCVI+2M/AFft4+PyldJowqNQTRbCPCVt8gS3Y4/NGYUcWvZx/I4DLgNnO8pXAj8DUEuyX\n8bFfd/5KfHI8E36dwEUNLiKu3SAqPrGRi9Jm8Y08QwveJCWiB08Nj7QjCmPMKbypWcwA7nMKz4hI\nfeADVb2uFPpXLFazONmhY4eYuHIi8cnxbD+4nQc6PMD9UfezY2lE3jTjGX8cos651cn44xDdbq1j\nicKYcsaX11mscb9YTkQqAKs8XUDnb5YsXNO3/LL9F+KTXZP4dWvaLW8Sv5/mHOWTEek0XvApj+97\nmTERQ/wy5GSMKTt8eZ1Fooj8ICL9RKQf8C0w60w7aIqnqDHTfUf28e5P7xI1Joo7Jt9BZHgkqwau\n4ss7v+TGljfy05yjLI4dSb1vxvL4vpcJ4zAD0ocyb2pG6QRQhGAfE7b4Aluwx+cNb86GekxEbgFy\n77v9vqp+UbLdMt5QVRZuXkh8cjxfr/uaHs178Oa1b3LluVdSQVzfA3JPj03bepxX04eynPaM5mEe\nZjRjIobQ7dY6fo7CGBMICh2Gcm5+NOtMJhMsTeVlGGpX5i7GLx9PQnICIkJcdBz3tr+XOmF/fvDn\nzu+UO+T0z5rvECHpeY+zu15p11AYY3xas0gEbtV8d7Eri4I5WeRoDokbE4lPjmfG7zO4+fybiYuO\n47LGl50yiV/u1diH0jP5K28SxmGyqMqzvVfTsFEFK2QbY/L4smZxCPhVRMaKyMjcnzPvovHGtgPb\n+M/8/9Dw8Yb8bdbfiImMIfWJVD64+QMub3K5x9le507ZzYD0oVzNLEbzMFlUZUzEEGIHRfDcf5uU\nyUQR7GPCFl9gC/b4vOHN3FBTsWsqStWJnBN899t3xCfHs3DzQvpe0JeXYl7ioT4PeTUVeEyf2oyZ\nMoQB6UOZXLM/z3ddaUNOxpgz4s0wVBWgubO4QVWPlHivTlOgD0Ol7ktlbPJYxqWMo3HNxsRFx9H3\ngr5Ur1S92NtKSsxi3tQMG3IyxhTqjGsWzn0rXgEeADbhmheqMX/e26LMzQ8ViMniWPYxvlr7FfHJ\n8STvSOaedvfQP7o/bSPaFntbuWc+xfSpbQnCGOMVX9Qs/g84GzhXVTuqajRwHhAOvOGbbpZf6zLW\n8fSMp2n0ViPe++U9+nXox9antjK8x3CPiaKoMdPcovbgUa1YHDuSpMSsEup5yQj2MWGLL7AFe3ze\nKCxZ9ATiVPVg7gpVPQA8AtxQ0h0LRoePH+aj5R/RdVxXun3QjRAJYdEDi5hz3xxiL4ylSmiV0952\nblG7rF1sZ4wJDoUNQ61X1ZbFfc6fyuow1PK05SQkJzBh5QQuaXgJcdFx9GrZi4ohFX22j9wjiwHp\nQ20aD2OM13xRs/gSmKqqH+Zbfw/QV1V7+6SnPlSWksXBowfzJvHbcWgHD3R4gAeiHqBpeFOf7se9\nTgFYUdsYUyy+qFk8CjwqInNF5E3nZx4wCNdQlMlHVflp20/EfR1Hk+FNmL5hOi/FvETq4FT+deW/\nzihReBozzV+nAMrsdRRFCfYxYYsvsAV7fN4o8DoLVd0GXCoiVwEXOKunq2piqfQsgOw9vJePV3xM\nfHI8Wcez6B/dn9UDV1O/Rv0S3e/cKbsZ7FanGDk1ls7dm5ToPo0x5VOR11kEktIchlJVFmxeQHxy\nPN+s+4brW1xPXHQcMZExeZP4lZTcoae6zapx8P/etzqFMea0+WxuqEBSGskiPTOd8SnjSViWQIiE\nEBcdx1/a/+WkSfx8LX9dwr2QXeOZh+zGRcaY0+bLuaHKvRzNYcbvM7h90u20fKclqzNWM+6mcawa\nuIonOz9Z4olicexIokY1Z3HsSCaOSDvpFNmMPw4FbJ3CXbCPCVt8gS3Y4/OG35KFiNQSkRkiss65\nuVLNQtpWEJFkEfm6NPu47cA2Xp73Ms1GNOO5Wc9xVeRVbHpiE+NuGudxtteSkHv9RBWOMSB9KDlS\ngTERQ/ImB7T7URhjSoPfhqFEZBiwW1VfF5FngVqq+lwBbZ8EOgJnFXbKri+GoU7knGD6b9OJT45n\n0eZF3HHBHcR1jCO6fvQZbfd0ebp+AuwUWWOMb5T5moWIrAW6qepOETkHmKuq53to1wjXfFT/AZ4q\nqWTxx94/GLvMNYlf05pN8ybxq1ap2mltz5dsUkBjTEkJhJpFhKruBFDVNCCigHZvA88APs9qR08c\n5fNVn3PNR9dwScIlZB7LZMY9M/jxwR+5P+r+MpEoADp3D6PT7RuDOlEE+5iwxRfYgj0+b3hzP4vT\nJiIzgXruq3B96L/gofkpyUBEbgR2qmqKiMQ4ry9Uv379iIyMBCA8PJwOHToQExMD/PmGx8TE8MqC\nV3j9k9c5N/xc/nb337il9S0sXriYXat35aUt9/a2bMu2bMvBsJz7ODU1leLw5zDUGiDGbRhqjqq2\nztfmFeAe4ARQFaiBawqSewvYptfDUNN/m07L2i1pfnbzohsbY0yQCoSaxTBgj6oOK6rA7bTvBvy1\npAvcxhhTngRCzWIYcI2IrAO6A68BiEh9EZnmx36VSe6HkMHI4gtsFl/wK9GaRWFUdQ9wtYf1O3Dd\nSyP/+nnAvFLomjHGmHxsuo8yxm6NaowpTYEwDGXyCfRboxpjgpclizKksFujBvuYqcUX2Cy+4GfJ\nogyJ6VPb5n0yxpRJVrPwI0/1CZvawxhTmsr8dRYlIZCShacJAi05GGNKmxW4y4ikxCxeHbjllGJ1\nYfUJT4J9zNTiC2wWX/CzZFGCCju7yeoTxphAYsNQJejVgVsYPKoVYRwmi6qMHLiW5/7bJO95q08Y\nY/zNahZlgNUljDFlndUsyoDO3cPoNGEQIweuPeNEEexjphZfYLP4gp/f5oYqLzp3D6Nz9yZFNzTG\nmDLMhqGMMaYcs2EoY4wxPmPJIkAE+5ipxRfYLL7gZ8nCGGNMkaxmYYwx5ZjVLIwxxviMJYsAEexj\nphZfYLP4gp8lC2OMMUWymoUxxpRjVrMwxhjjM5YsAkSwj5lafIHN4gt+fksWIlJLRGaIyDoR+UFE\nahbQrqaITBKRNSKySkQuLe2+GmNMeee3moWIDAN2q+rrIvIsUEtVn/PQ7gNgnqqOE5FQIExVDxSw\nTatZGGNMMZT5+1mIyFqgm6ruFJFzgLmqen6+NmcBy1T1PC+3acnCGGOKIRAK3BGquhNAVdOACA9t\nzgUyRGSciCSLyPsiUrVUe1lGBPuYqcUX2Cy+4Fei97MQkZlAPfdVgAIveGju6ZAgFIgGHlXVX0Rk\nOPAc8GJB++zXrx+RkZEAhIeH06FDB2JiYoA/33BbtmVbtuXyupz7ODU1leLw5zDUGiDGbRhqjqq2\nztemHpCkqs2c5S7As6raq4Bt2jCUMcYUQyAMQ30N9HMe3wd8lb+BM0y1RURaOqu6A6tLpXfGGGPy\n+DNZDAOuEZF1uJLAawAiUl9Eprm1GwR8IiIpQHvglVLvaRngfggZjCy+wGbxBT+/3YNbVfcAV3tY\nvwPo6ba8HLi4FLtmjDEmH5sbyhhjyrFAqFkYY4wJEJYsAkSwj5lafIHN4gt+liyMMX71yy+/MG/e\nPF5//XV/d8UUwpJFgMi9sCZYWXyBzZv4Xn31VVq0aMHYsWMZPnw4jzzyCIcPH2bp0qV06tSJjIwM\nMjMzT2v/L7/8Ml9//TWvvOL5ZMmcnBxeeeUVPv30U+Lj4/PWqypPPfVUoe1ycnL48ccfT3lteWPJ\nwhhTKi6++GJuvfVWHnzwQZ544gnS0tKYNWsWAwYMoGLFiuTk5FCtWrVibzcxMRGA3r17c/z4cRYu\nXHhKm08//ZQmTZpw1113sWHDBrZs2cLevXsZPnw48+fPL7Dd5s2bPb62PLJkcQaSErN4deAWkhKz\nSnxfwT5mavEFNm/iW7JkSd4RSHp6Onv27OHyyy8HYPLkyfz973/nxIkTxd73okWLiIqKAiAqKorZ\ns2d7bNOoUSMAmjZtyoIFC6hVqxZPPvkkZ511VqHtFi1aREZGxknryiNLFqcpKTGLxbEjGTyqFYtj\nR5ZKwjAmkP3yyy8cOXKEUaNG8fbbb/P9999z9tlnM2HCBGbMmMHf//53KlQo/kdSenp63hFJ9erV\nSUtLO6VNjRo18hKRqrJt27a859xPt/fUrkaNGmRnZ3t8bXnit4vyAt3cKbsZnD6UMA4zIH0oI6fG\n0rl7kxLbn415BzaLD/bs2cMtt9wCQLdu3ahcuTIAsbGxxMbGntJ+9erVzJw5E5FTLwG47777qFnT\ndb+0nJwcQkJCAMjOzs577O6ee+5hwYIFXH311axYsYKWLVvmPee+fU/tctcBp7y2PLFkcZpi+tRm\nzJQhDEgfypiIIXS7tY6/u2RMmbV582bOOeeck5aPHj1K1aoF33GgTZs2tGnTpsht16tXL68wfuDA\nAerWrXtKmwsvvJDdu3fz3Xff0bBhQ9q2betxW57aefvaYGfDUKepc/cwOk0YxMiBa+k0YRCdu4eV\n6P5szDuwlff4lixZQvv27QE4duwYO3bsoGrVqqSnpxf4mtWrVzNixIhTfkaOHMm+ffvy2nXp0oUV\nK1YA8NNPP9GpUycANm3alNdmxowZbNy4keuvv57du3fTvXv3vOfch6E8tZsxYwbTp0/3+NryxI4s\nzkDn7mElOvRkTDCYP38+o0ePplGjRuzatYu6devSq1cvJk2aROvWrYmI8HTfM++PLK666iq+++47\nJk+ejIhw7bXXsm/fPmJjY1m0aBEALVq0YM2aNYwaNYq+ffsSGhpKZmYm8fHxrF27luHDh/PQQw95\nbNeiRQu+/vrrk9aVR+V2bqikxCzmTtlNTJ/aJX5UYIwxZZXNDVUIO5PJGGOKp1wmi7lTdjPA7Uym\neVMz/N2lIpX3Me9AZ/EFtmCPzxvlMlnE9KnNmIghZFHVzmQyxhgvlOuaxbypGXS7tY7VLIwx5Za3\nNYtymyyMMcZYgTvoBPuYqcUX2Cy+4GfJwhhjTJFsGMoYY8oxG4YyxhjjM35LFiJSS0RmiMg6EflB\nRGoW0O5JEVkpIitE5BMRqVTafS0Lgn3M1OILbBZf8PPnkcVzwCxVbQXMBv6ev4GINAAeB6JVtR2u\nuazuLNVelhEpKSn+7kKJsvgCm8UX/PyZLG4CxjuPxwM3F9AuBKgmIqFAGLC9FPpW5rjPshmMLL7A\nZvEFP38miwhV3QmgqmnAKVNPqup24E1gM7AN2Keqs0q1l8YYY0p2inIRmQnUc18FKPCCh+annMYk\nIuG4jkCaAvuBySISq6oTSqC7ZVpqaqq/u1CiLL7AZvEFP7+dOisia4AYVd0pIucAc1S1db42twHX\nqWqcs/wX4FJVfayAbdp5s8YYU0zenDrrz7t4fA30A4YB9wFfeWizGegkIlWAo0B34OeCNuhNwMYY\nY4rPn0cWZwOfA42BTUBfVd0nIvWBeFXt6bR7EdcZUMeBZUB/VT3ul04bY0w5FVRXcBtjjCkZQXUF\nt4gMFZHlIrJMRL53aiFBQ0ReF5E1IpIiIlNE5Cx/98mXROQ25wLMbBGJ9nd/fEFEeojIWhFZLyLP\n+rs/viYiY0Vkp4is8HdffE1EGonIbBFZJSK/isggf/fJl0SksogscT4vf3VGcQpuH0xHFiJSXVUP\nOY8fB9qo6iN+7pbPiMjVwGxVzRGR1wBV1VMuZgxUItIKyAHGAE+rarKfu3RGRKQCsB5XrW07rnrb\nnaq61q8d8yER6QIcAj50LpwNGs6XzXNUNUVEqgNLgZuC7P0LU9UsEQkBFgGDVPUnT22D6sgiN1E4\nquH64AkaqjpLVXNjWgw08md/fE1V16nqb7hOsQ4GlwC/qeomp842Edep4EFDVRcCe/3dj5Kgqmmq\nmuI8PgSsARr6t1e+papZzsPKuE54KvDoIaiSBYCI/FtENgOxwBB/96cEPQB85+9OmEI1BLa4LW8l\nyD5sygsRiQQ6AEv82xPfEpEKIrIMSANmqmqBZ5sGXLIQkZnOpIK5P786//YCUNUXVLUJ8AmueaUC\nSlHxOW3+ARwPxIsTvYnPmLLEGYKaDAzON3oR8FQ1R1WjcI1SXCoibQpq68/rLE6Lql7jZdMJwHTg\npZLrje8VFZ+I9ANuAK4qlQ75WDHev2CwDWjittzIWWcChDMn3WTgI1X1dC1YUFDVAyIyB+gBrPbU\nJuCOLAojIs3dFm/GNcYYNESkB/AM0FtVj/q7PyUsGOoWPwPNRaSpM7X+nbguRg02QnC8X578D1it\nqiP83RFfE5E6ubeGEJGqwDVAgcX7YDsbajLQEldhexPwsKru8G+vfEdEfgMqAbudVYtVdaAfu+RT\nInIz8A5QB9gHpKjq9f7t1ZlxEvwIXF/Mxqrqa37ukk+JyAQgBqgN7AReVNVxfu2Uj4jI5cB84Fdc\nhV8FnlfV7/3aMR8RkQtxzfhdwfn5TFX/U2D7YEoWxhhjSkZQDUMZY4wpGZYsjDHGFMmShTHGmCJZ\nsjDGGFMkSxbGGGOKZMnCGGNMkSxZGONGRBqKyJfOlOIbRGSkiFT08T66iUhnt+UBInKP83iciNzq\ny/0Z4wuWLIw52VRgqqq2BFoAYcD/+XgfMcBluQuqOkZVP/bxPozxKUsWxjhE5CrgsKp+CK6bhQBP\nAveKyKMi8o5b229EpKvz+D0R+Sn/DWRE5A8ReUlEljo35WopIk2Bh4EnRCRZRC4XkRdF5CkP/YkW\nkbki8rOIfCci9Zz1g5wb8qQ4V1AbU+ICbiJBY0rQBbhucJNHVQ+KSCoQQsFz/T/v3D++ApAoIlNU\ndaXzXLqqdhSRR3Dd0OkhERkNHFTVtyDvplYncSawewfXPGC7RaQv8ArwIPAsEKmqx4Ptbomm7LJk\nYcyZu1NE4nD9fzoHaAPkJosvnH+XArcUY5utgLbATBERXKMA253nlgMTRORL4Msz7LsxXrFkYcyf\nVgO3ua9wvrnXwzV5Y0u3p6o4z0cCfwU6OtM8j8t9zpE7O3A2xfv/JsBKVb3cw3M3Al2B3sA/RKSt\n2x0UjSkRVrMwxqGqiUBVtzOTQoA3cA0HpQJR4tIY1y1TAc7CdQ/qg05NwZtZcg86ryvMOqCuiHRy\n+hLqdmOaJqo6D3jO2U51L0M05rRZsjDmZLcAt4vIeiADyFbV11R1EfAHsAoYjlPbUNUVQAque6d8\nDCx021ZBNY5vgFtyC9z52qmz3eO4jnKGiUgKsAzo7NQyPhaR5U4fRqjqAR/EbUyhbIpyYwrgfKv/\nFLhFVVP83R9j/MmShTHGmCLZMJQxxpgiWbIwxhhTJEsWxhhjimTJwhhjTJEsWRhjjCmSJQtjjDFF\nsmRhjDGmSP8Pl8wXRqXRQ/sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa9aaf80c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                Y\n",
      "count  256.000000\n",
      "mean    -0.022010\n",
      "std      0.142689\n",
      "min     -0.720196\n",
      "25%     -0.094655\n",
      "50%     -0.011186\n",
      "75%      0.063604\n",
      "max      0.511783\n",
      "kurtosis  7.51103\n"
     ]
    }
   ],
   "source": [
    "#  Generate GM(2) mixture, with unconstrained mean mu, and use\n",
    "#  parameters from XYZ example to generate a small sample array:\n",
    "mu = 0\n",
    "sigma = 0.13\n",
    "b = 2  # By construction, b > 1.\n",
    "mumix = mu + simug_mix(sigma1=0.096896, sigma2=b*sigma, q=0.12903, N=256)\n",
    "\n",
    "plotqq( mumix )\n",
    "stat( mumix )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the observed values truly came from a normal distribution,\n",
    "most of the plotted points should be along the solid green line.\n",
    "Points far from the line are usually identified as *outliers*, but\n",
    "for leptokurtotic distributions such points in the tails are to be expected.\n",
    "Statistical tests for normality are essentially testing\n",
    "deviations from the green line ($R^2$ shown within the plot\n",
    "derives from linear regression and should not be\n",
    "considered a serious test statistic).\n",
    "\n",
    "The last cell simulates a small-sample (N=256, say, one year of daily data)\n",
    "where kurtosis is exactly 7 in theory.\n",
    "***We highly recommend re-running the cell using different parameter values***,\n",
    "including unconstrained mean mu, **to visually understand\n",
    "how Q-Q plots can appear so vastly different in practice\n",
    "even when the generating dynamics are fully known**.\n",
    "When applying *plotqq()* to actual returns of financial assets,\n",
    "one should be able to recognize plots where a GM(2) model\n",
    "is appropriate.\n",
    "\n",
    "The *stat(mix)* output illustrates how much statistical estimates of a\n",
    "Gaussian mixture sample can differ from their theoretical values,\n",
    "especially when the sample size is small.\n",
    "Such variability is noteworthy when we base portfolio decisions\n",
    "on estimated statistics mistaken for being stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "##  Appendix 3: Fitting Gaussian mixtures using scikit-learn\n",
    "\n",
    "This notebook imposed structure on GM(2) to discern its\n",
    "parametric relationships in theory. From two observable statistics\n",
    "we deduced a reasonable shape for a Gaussian mixture\n",
    "stylized by actual financial returns.\n",
    "\n",
    "Given several assets, our plain $\\sigma$ must be generalized\n",
    "to a covariance matrix. Python via *scikit-learn* offers\n",
    "code to fit unconstrained mixtures.\n",
    "For starters, please see to http://scikit-learn.org/stable/modules/mixture.html\n",
    "\n",
    "Tiago Ramalho wrote a\n",
    "[quick introduction](http://www.nehalemlabs.net/prototype/blog/2014/04/03/quick-introduction-to-gaussian-mixture-models-with-python)\n",
    "to Gaussian mixture models with Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Cohen, A.C. (1967) Estimation in mixtures of two normal distributions.\n",
    "*Technometrika* 9:15–28.\n",
    "\n",
    "Pearson, K. (1894) Contributions to the theory of mathematical evolution.\n",
    "*Philosophical Trans. Roy. Soc. London* A185:71–110.\n",
    "\n",
    "Wang, Jin (2015), [Multivariate Mixtures of Normal Distributions: Properties, Random Vector Generation, Fitting, and as Models of Market Daily Changes](https://www.researchgate.net/profile/Jin_Wang40/publication/274695842_Multivariate_Mixtures_of_Normal_Distributions_Properties_Random_Vector_Generation_Fitting_and_as_Models_of_Market_Daily_Changes/links/5525314c0cf22e181e73ee4f.pdf), Journal on Computing, Vol. 27, No. 2, Spring 2015, pp. 193–203. Cf. earlier version, Wang (2000), [Modeling and Generating Daily Changes in Market Variables Using A Multivariate Mixture of Normal Distributions](http://ww2.valdosta.edu/~jwang/paper/MixNormal.pdf).\n",
    "\n",
    "Thanks for corrected [skew equation](http://onlyvix.blogspot.com/2011/05/mixture-of-normal-formulas-for-skew-and.html) to @onlyvixx, and jwang@valdosta.edu."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}