Notebooks/GMM/GMMest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Method of Moments (GMM) Estimation\n",
    "### by [Richard W. Evans](https://sites.google.com/site/rickecon/), February 2019\n",
    "The code in this Jupyter notebook was written using Python 3.6. It also uses data file `Econ381totpts.txt` and the image file `MLEplots.png`. This data file and image file are stored along with the Jupyter notebook file `GMMest.ipynb` in the public GitHub repository [https://github.com/rickecon/Notebooks/tree/master/GMM](https://github.com/rickecon/Notebooks/tree/master/GMM).\n",
    "\n",
    "For teaching and learning purposes, it is best to clear all cell output by selecting `Cell` then `All Output` then `Clear` from the menu across the top of this notebook. However, this notebook comes with all the cells' output displayed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. MLE vs. GMM: Strengths and weaknesses\n",
    "A paper by Fuhrer, et al (1995) studies the accuracy and efficiency of the maximum likelihood (ML) estimator versus the generalized method of moments (GMM) estimator in the context of a simple linear-quadratic inventory model. They find that ML has some very nice properties over GMM in small samples when the model is simple. In the spirit of the Fuhrer, et al (1995) paper, we list the strengths and weaknesses of MLE vs. GMM more generally. I recommend you read the introduction to Fuhrer, et al (1995)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. MLE strengths\n",
    "* More statistical significance. In general, MLE provides more statistical significance for parameter estimates than does GMM. This comes from the strong distributional assumptions that are necessary for the ML estimates.\n",
    "* ML estimates are less sensitive to parameter or model normalizations than are GMM estimates.\n",
    "* ML estimates have nice small sample properties. ML estimates have less bias and more efficiency with small data samples than GMM estimates in many cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. MLE weaknesses\n",
    "* MLE requires strong distributional assumptions. For MLE, the data generating process (DGP) must be completely specified. This assumes a lot of knowledge about the DGP. This assumption is likely almost always wrong.\n",
    "* MLE is very difficult in rational expectations models. This is because the consistency of beliefs induces a nonlinearity in the likelihood function that makes it difficult to find the global optimum.\n",
    "* MLE is very difficult in nonlinear models. The likelihood function can become highly nonlinear in MLE even if the model is linear when the data are irregular. This difficulty is multiplied when the model itself is more complicated and nonlinear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. GMM strengths\n",
    "* GMM allows for most flexible identification. GMM estimates can be identified by any set of moments from the data as long as you have at least as many moments as you have parameters to estimate and that those moments are independent enough to identify the parameters. (And the parameters are independent enough of each other to be separately identified.)\n",
    "* Good large sample properties. The GMM estimator is strongly consistent and asymptotically normal. GMM will likely be the best estimator if you have a lot of data.\n",
    "* GMM requires minimal assumptions about the DGP. In GMM, you need not specify the distributions of the error terms in your model of the DGP. This is often a strength, given that most error are not observed and most models are gross approximations of the true DGP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4. GMM weaknesses\n",
    "* GMM estimates are usually less statistically significant than ML estimates. This comes from the minimal distributional assumptions. GMM parameter estimates usually are measured with more error.\n",
    "* GMM estimates can be sensitive to normalizations of the model or parameters.\n",
    "* GMM estimates have bad small sample properties. GMM estimates can have large bias and inefficiency in small samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5. Key questions to answer when deciding between MLE and GMM\n",
    "1. How much data is available for the estimation? Large data samples will make GMM relatively more attractive than MLE because of the nice large sample properties of GMM and fewer required assumptions on the model\n",
    "2. How complex is the model? Linear models or quadratic models are much easier to do using MLE than are more highly nonlinear models. Rational expectations models (macroeconomics) create an even more difficult level of nonlinearity that pushes you toward GMM estimation.\n",
    "3. How comfortable are you making strong distributional assumptions? MLE requires a complete specification of all distributional assumptions of the model DGP. If you think these assumptions are too strong, you should use GMM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. The GMM estimator\n",
    "GMM was first formalized by Hansen (1982). A strength of GMM estimation is that the econometrician can remain completely agnostic as to the distribution of the random variables in the DGP. For identification, the econometrician simply needs at least as many moment conditions from the data as he has parameters to estimate.\n",
    "\n",
    "A *moment* of the data is broadly defined as any statistic that summarizes the data to some degree. A data moment could be as narrow as an individual observation from the data or as broad as the sample average. GMM estimates the parameters of a model or data generating process to make the model moments as close as possible to the corresponding data moments. See Davidson and MacKinnon (2004, ch. 9) for a more detailed treatment of GMM. The estimation methods of linear least squares, nonlinear least squares, generalized least squares, and instrumental variables estimation are all specific cases of the more general GMM estimation method.\n",
    "\n",
    "Let $m(x)$ be an $R\\times 1$ vector of moments from the real world data $x$. And let $x$ be an $N\\times K$ matrix of data with $K$ columns representing $K$ variables and $N$ observations. Let the model DGP be characterized as $F(x,\\theta)$, where $F$ is a vector of equations, each of which is a function of the data $x$ and the $K\\times 1$ parameter vector $\\theta$. Then define $m(x|\\theta)$ as a vector of $R$ moments from the model that correspond to the real-world moment vector $m(x)$. Note that GMM requires both real world data $x$ and deterministic moments from the model $m(x|\\theta)$ in order to estimate $\\hat{\\theta}_{GMM}$. There is also a stochastic way to generate moments from the model, which we discuss later in our section on Simulated Method of Moments (SMM).\n",
    "\n",
    "The GMM approach of estimating the parameter vector $\\hat{\\theta}_{GMM}$ is to choose $\\theta$ to minimize some distance measure of the data moments $m(x)$ from the model moments $m(x|\\theta)$.\n",
    "\n",
    "$$ \\hat{\\theta}_{GMM}=\\theta:\\quad \\min_{\\theta}\\: ||m(x) - m(x|\\theta)|| $$\n",
    "\n",
    "The distance measure $||m(x)-m(x|\\theta)||$ can be any kind of norm. But it is important to recognize that your estimates $\\hat{\\theta}_{GMM}$ will be dependent on what distance measure (norm) you choose. The most widely studied and used distance metric in GMM estimation is the $L^2$ norm or the sum of squared errors in moments. Define the moment error function $e(x|\\theta)$ as either the percent difference in the vector of model moments from the data moments or the simple difference.\n",
    "\n",
    "$$ e(x|\\theta) \\equiv \\frac{m(x) - m(x|\\theta)}{m(x|\\theta)} \\quad\\text{or}\\quad m(x)-m(x|\\theta) $$\n",
    "\n",
    "It is important when possible that the error function $e(x|\\theta)$ be a percent deviation of the moments (given that none of the data moments are 0). This puts all the moments in the same units, which helps make sure that no moments receive unintended weighting simply due to its units. This ensures that the problem is scaled properly and does not suffer from ill conditioning. However, percent deviations become computationally problematic when the data moments are zero or close to zero. In that case, you would use a simple difference.\n",
    "\n",
    "The GMM estimator is the following,\n",
    "\n",
    "$$ \\hat{\\theta}_{GMM}=\\theta:\\quad \\min_{\\theta}\\:e(x|\\theta)^T \\, W \\, e(x|\\theta) $$\n",
    "\n",
    "where $W$ is a $R\\times R$ weighting matrix in the criterion function. For now, think of this weighting matrix as the identity matrix. But we will show in Section 3 a more optimal weighting matrix. We call the quadratic form expression $e(x|\\theta)^T \\, W \\, e(x|\\theta)$ the *criterion function* because it is a strictly positive scalar that is the object of the minimization in the GMM problem statement. The $R\\times R$ weighting matrix $W$ in the criterion function allows the econometrician to control how each moment is weighted in the minimization problem. For example, an $R\\times R$ identity matrix for $W$ would give each moment equal weighting, and the criterion function would be a simply sum of squared percent deviations (errors). Other weighting strategies can be dictated by the nature of the problem or model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. The Weighting Matrix (W)\n",
    "In the GMM criterion function in the problem statement above, some weighting matrices $W$ produce precise estimates while others produce poor estimates with large variances. We want to choose the optimal weighting matrix $W$ with the smallest possible asymptotic variance. This is an efficient of optimal GMM estimator. The optimal weighting matrix is the inverse variance covariance matrix of the moments at the optimal parameter values,\n",
    "\n",
    "$$ W^{opt} \\equiv \\Omega^{-1}(x|\\hat{\\theta}_{GMM}) $$\n",
    "\n",
    "where $\\Omega(x|\\theta)$ is the variance covariance matrix of the moment condition errors $E(x|\\theta)$ from each observation in the data (to be defined below). The intuition for using the inverse variance covariance matrix $\\Omega^{-1}$ as the optimal weighting matrix is the following. You want to downweight moments that have a high variance, and you want to weight more heavily the moments that are generated more precisely.\n",
    "\n",
    "Notice that this definition of the optimal weighting matrix is circular. $W^{opt}$ is a function of the GMM estimates $\\hat{\\theta}_{GMM}$, but the optimal weighting matrix is used in the estimation of $\\hat{\\theta}_{GMM}$. This means that one has to use some kind of iterative fixed point method to find the true optimal weighting matrix $W^{opt}$. Below are some examples of weighting matrices to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. The identity matrix (W = I)\n",
    "Many times, you can get away with just using the identity matrix as your weighting matrix $W = I$. This changes the criterion function to a simple sum of squared error functions such that each moment has the same weight.\n",
    "\n",
    "$$ \\hat{\\theta}_{GMM}=\\theta:\\quad \\min_{\\theta}\\:e(x|\\theta)^T \\, e(x|\\theta) $$\n",
    "\n",
    "If the problem is well conditioned and well identified, then your GMM estimates $\\hat{\\theta}_{GMM}$ will not be greatly affected by this simplest of weighting matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Two-step variance covariance estimator of W\n",
    "The most common method of estimating the optimal weighting matrix for GMM estimates is the two-step variance covariance estimator. The name \"two-step\" refers to the two steps used to get the weighting matrix.\n",
    "\n",
    "The first step is to estimate the GMM parameter vector $\\hat{\\theta}_{1,GMM}$ using the simple identity matrix as the weighting matrix $W = I$.\n",
    "\n",
    "$$ \\hat{\\theta}_{1, GMM}=\\theta:\\quad \\min_{\\theta}\\:e(x|\\theta)^T \\, I \\, e(x|\\theta) $$\n",
    "\n",
    "You then use the errors from each observation in the data $E(x|\\theta)$ that generated the error vector $e(x|\\theta)$, which data matrix is $R\\times N$ and if summed across the rows gives the error vector. When the errors are simple differences the $E(x|\\theta)$ matrix is the following,\n",
    "\n",
    "$$ E(x|\\theta) =\n",
    "  \\begin{bmatrix}\n",
    "    m_1(x|\\theta) - m_1(x_1) & m_1(x|\\theta) - m_1(x_2) & ... & m_1(x|\\theta) - m_1(x_N) \\\\\n",
    "    m_2(x|\\theta) - m_2(x_1) & m_2(x|\\theta) - m_2(x_2) & ... & m_2(x|\\theta) - m_2(x_N) \\\\\n",
    "    \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    m_R(x|\\theta) - m_R(x_1) & m_R(x|\\theta) - m_R(x_2) & ... & m_R(x|\\theta) - m_R(x_N) \\\\\n",
    "  \\end{bmatrix} $$\n",
    "\n",
    "where $m_r(x_i)$ is a function associated with the $r$th moment and the $i$th data observation. When the errors are percent deviations, the $E(x|\\theta)$ matrix is the following,\n",
    "\n",
    "$$ E(x|\\theta) =\n",
    "  \\begin{bmatrix}\n",
    "    \\frac{m_1(x_1) - m_1(x|\\theta)}{m_1(x|\\theta)} & \\frac{m_1(x_2) - m_1(x|\\theta)}{m_1(x|\\theta)} & ... & \\frac{m_1(x_N) - m_1(x|\\theta)}{m_1(x|\\theta)} \\\\\n",
    "    \\frac{m_2(x_1) - m_2(x|\\theta)}{m_2(x|\\theta)} & \\frac{m_2(x_2) - m_2(x|\\theta)}{m_2(x|\\theta)} & ... & \\frac{m_2(x_N) - m_2(x|\\theta)}{m_2(x|\\theta)} \\\\\n",
    "    \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    \\frac{m_R(x_1) - m_R(x|\\theta)}{m_R(x|\\theta)} & \\frac{m_R(x_2) - m_R(x|\\theta)}{m_R(x|\\theta)} & ... & \\frac{m_R(x_N) - m_R(x|\\theta)}{m_R(x|\\theta)} \\\\\n",
    "  \\end{bmatrix} $$\n",
    "\n",
    "where the denominator of the percentage deviation or baseline is the model moment that does not change. We use the $E(x|\\theta)$ data matrix and the Step 1 GMM estimate $e(x|\\hat{\\theta}_{1,GMM})$ to get a new estimate of the variance covariance matrix.\n",
    "\n",
    "$$ \\hat{\\Omega}_2 = \\frac{1}{N}E(x|\\hat{\\theta}_{1,GMM})\\,E(x|\\hat{\\theta}_{1,GMM})^T $$\n",
    "\n",
    "This is simply saying that the $(r,s)$-element of the estimator of the variance-covariance matrix of the moment vector is the following.\n",
    "\n",
    "$$\\hat{\\Omega}_{r,s} = \\frac{1}{N}\\sum_{i=1}^N\\Bigl[m_r(x|\\theta) - m_{r}(x_i)\\Bigr]\\Bigl[ m_s(x|\\theta) - m_s(x_i)\\Bigr] $$\n",
    "\n",
    "The optimal weighting matrix is the inverse of the two-step variance covariance matrix.\n",
    "\n",
    "$$ \\hat{W}^{two-step} \\equiv \\hat{\\Omega}_2^{-1} $$\n",
    "\n",
    "Lastly, re-estimate the GMM estimator using the optimal two-step weighting matrix.\n",
    "\n",
    "$$ \\hat{\\theta}_{2, GMM}=\\theta:\\quad \\min_{\\theta}\\:e(x|\\theta)^T \\, \\hat{W}^{two-step} \\, e(x|\\theta) $$\n",
    "\n",
    "$\\hat{\\theta}_{2, GMM}$ is called the two-step GMM estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Iterated variance covariance estimator of W\n",
    "The truly optimal weighting matrix $W^{opt}$ is the iterated variance-covariance estimator of $W$. This procedure is to just repeat the process described in the two-step GMM estimator until the estimated weighting matrix no longer changes between iterations. Let $i$ index the $i$th iterated GMM estimator,\n",
    "\n",
    "$$ \\hat{\\theta}_{i, GMM}=\\theta:\\quad \\min_{\\theta}\\:e(x|\\theta)^T \\, \\hat{W}_{i} \\, e(x|\\theta) $$\n",
    "\n",
    "and the $i+1$th estimate of the optimal weighting matrix is defined as the following.\n",
    "\n",
    "$$ \\hat{W}_{i+1} \\equiv \\hat{\\Omega}_{i+1} = \\frac{1}{N}E(x|\\hat{\\theta}_{i,GMM})\\,E(x|\\hat{\\theta}_{i,GMM})^T $$\n",
    "\n",
    "The iterated GMM estimator is the $\\hat{\\theta}_{i,GMM}$ such that $\\hat{W}_{i+1}$ is very close to $\\hat{W}_{i}$ for some distance metric (norm).\n",
    "\n",
    "$$ \\hat{\\theta}_{it,GMM} = \\hat{\\theta}_{i,GMM}: \\quad || \\hat{W}_{i+1} - \\hat{W}_{i} || < \\varepsilon $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4. Newey-West consistent estimator of $\\Omega$ and W\n",
    "The Newey-West estimator of the optimal weighting matrix and variance covariance matrix is consistent in the presence of heteroskedasticity and autocorrelation in the data (See Newey and West, 1987). Adda and Cooper (2003, p. 82) have a nice exposition of how to compute the Newey-West weighting matrix $\\hat{W}_{nw}$. The asymptotic representation of the optimal weighting matrix $\\hat{W}^{opt}$ is the following:\n",
    "\n",
    "$$ \\hat{W}^{opt} = \\lim_{N\\rightarrow\\infty}\\frac{1}{N}\\sum_{i=1}^N \\sum_{l=-\\infty}^\\infty E(x_i|\\theta)E(x_{i-l}|\\theta)^T $$\n",
    "\n",
    "The Newey-West consistend estimator of $\\hat{W}^{opt}$ is:\n",
    "\n",
    "$$ \\hat{W}_{nw} = \\Gamma_{0,N} + \\sum_{v=1}^q \\left(1 - \\left[\\frac{v}{q+1}\\right]\\right)\\left(\\Gamma_{v,N} + \\Gamma^T_{v,N}\\right) $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\Gamma_{v,N} = \\frac{1}{N}\\sum_{i=v+1}^N E(x_i|\\theta)E(x_{i-v}|\\theta)^T $$\n",
    "\n",
    "Of course, for autocorrelation, the subscript $i$ and the sample size $N$ can be changed to $t$ and $T$, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. The GMM Variance-Covariance Estimator of the Estimated Parameters\n",
    "The estimated variance-covariance matrix $\\hat{\\Sigma}$ of the estimated parameter vector $\\hat{\\theta}_{GMM}$ is different from the variance-covariance matrix $\\hat{\\Omega}$ of the moment vector $e(x|\\theta)$ from the previous section. $\\hat{\\Omega}$ from the previous section is the $R\\times R$ variance-covariance matrix of the $R$ moment errors used to identify the $K$ parameters $\\theta$ to be estimated. The estimated variance-covariance matrix of the estimated parameter vector $\\hat{\\Sigma}$ is a $K\\times K$ matrix. We say the model is exactly identified if $K = R$. We say the model is overidentified if $K<R$.\n",
    "\n",
    "Similar to the inverse Hessian estimator of the variance-covariance matrix of the maximum likelihood estimator from the [maximum likelihood estimation notebook](https://github.com/rickecon/Notebooks/blob/master/MLE/MLest.ipynb), the GMM variance-covariance matrix is related to the derivative of the criterion function with respect to each parameter. The intuition is that if the second derivative of the criterion function with respect to the parameters is large, there is a lot of curvature around the criterion minimizing estimate. In other words, the parameters of the model are precisely estimated. The inverse of the Hessian matrix will be small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define $R\\times K$ matrix $d(x|\\theta)$ as the Jacobian matrix of derivatives of the $R\\times 1$ error vector $e(x|\\theta)$.\n",
    "\n",
    "\\begin{equation}\n",
    "  d(x|\\theta) \\equiv\n",
    "    \\begin{bmatrix}\n",
    "      \\frac{\\partial e_1(x|\\theta)}{\\partial \\theta_1} & \\frac{\\partial e_1(x|\\theta)}{\\partial \\theta_2} & ... & \\frac{\\partial e_1(x|\\theta)}{\\partial \\theta_K} \\\\\n",
    "      \\frac{\\partial e_2(x|\\theta)}{\\partial \\theta_1} & \\frac{\\partial e_2(x|\\theta)}{\\partial \\theta_2} & ... & \\frac{\\partial e_2(x|\\theta)}{\\partial \\theta_K} \\\\\n",
    "      \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "      \\frac{\\partial e_R(x|\\theta)}{\\partial \\theta_1} & \\frac{\\partial e_R(x|\\theta)}{\\partial \\theta_2} & ... & \\frac{\\partial e_R(x|\\theta)}{\\partial \\theta_K}\n",
    "    \\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "The GMM estimates of the parameter vector $\\hat{\\theta}_{GMM}$ are assymptotically normal. If $\\theta_0$ is the true value of the parameters, then the following holds,\n",
    "\n",
    "\\begin{equation}\n",
    "  \\text{plim}_{N\\rightarrow\\infty}\\sqrt{N}\\left(\\hat{\\theta}_{GMM} - \\theta_0\\right) \\sim \\text{N}\\left(0, \\left[d(x|\\theta)^T W d(x|\\theta)\\right]^{-1}\\right)\n",
    "\\end{equation}\n",
    "\n",
    "where $W$ is the optimal weighting matrix from the GMM criterion function. The GMM estimator for the variance-covariance matrix $\\hat{\\Sigma}_{GMM}$ of the parameter vector $\\hat{\\theta}_{GMM}$ is the following.\n",
    "\n",
    "\\begin{equation}\n",
    "  \\hat{\\Sigma}_{GMM} = \\frac{1}{N}\\left[d(x|\\theta)^T W d(x|\\theta)\\right]^{-1}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the examples below, we will use a finite difference method to compute numerical versions of the Jacobian matrix $d(x|\\theta)$.\n",
    "\n",
    "$$ f'(x_0) = \\lim_{h\\rightarrow 0} \\frac{f(x_0 + h) - f(x_0)}{h} $$\n",
    "\n",
    "The following is a centered second-order numerical approximation of the derivative of a function. (See [BYU ACME numerical differentiation lab](https://github.com/UC-MACSS/persp-model-econ_W19/blob/master/Notes/ACME_NumDiff.pdf) for more details.)\n",
    "\n",
    "$$ f'(x_0) \\approx \\frac{f(x_0 + h) - f(x_0 - h)}{2h} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Examples\n",
    "In this section, we will use GMM to estimate parameters of the models from the [maximum likelihood estimation notebook](https://github.com/rickecon/Notebooks/blob/master/MLE/MLest.ipynb). We will also go through the standard moment conditions in most econometrics textbooks in which the conditional and unconditional expectations provide moments for estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Fitting a truncated normal to intermediate macroeconomics test scores\n",
    "Let's revisit the problem from the MLE notebook of fitting a truncated normal distribution to intermediate macroeconomics test scores. The data are in the text file `Econ381totpts.txt`. Recall that these test scores are between 0 and 450. The figure below shows a histogram of the data, as well as three truncated normal PDF's. The black line is the ML estimate of $\\mu$ and $\\sigma$ of the truncated normal pdf. The red and the green lines are just the PDF's of two \"arbitrarily\" chosen combinations of the truncated normal parameters $\\mu$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages and load the data\n",
    "import numpy as np\n",
    "import numpy.linalg as lin\n",
    "import scipy.stats as sts\n",
    "import scipy.integrate as intgr\n",
    "import scipy.optimize as opt\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "cmap1 = matplotlib.cm.get_cmap('summer')\n",
    "# This next command is specifically for Jupyter Notebook\n",
    "%matplotlib notebook\n",
    "\n",
    "import requests\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlcVFX/B/DPsIPIqoIgCu6CO5rh8riDhkRaiVZm2vIzc8uyMit9slDryeWxsizTfCq13NLKBUsFFTcEUjE0RVERcQXZYTi/P6a5zMAwzAwDMwOf9+s1r673nnvu94w3+XLuOefKhBACRERERNRgWJk6ACIiIiKqW0wAiYiIiBoYJoBEREREDQwTQCIiIqIGhgkgERERUQPDBJCIiIiogWECSERERNTAMAEkIiIiamCYABIRERE1MEwAiYiIiBoYJoBEREREDQwTQCIiIqIGhgkgERERUQPDBJCIiIiogWECSERERNTAMAEkIiIiamCYABIRERE1MEwAiYiIiBoYJoBEREREDQwTQCIiIqIGhgkgERERUQPDBJCIiIiogWECSERERNTAMAEkIiIiamCYABIRERE1MEwAiYiIiBoYJoBEREREDQwTQJIcOHAAMplMr8+sWbNMHbZJ/frrr1iwYAHCw8PRqVMnNGnSBLa2tnB3d0dwcDBee+01pKam6lRXcXEx1qxZgxEjRqB58+awt7eHs7MzOnTogMmTJ+Po0aNazy8rK0NKSgrWrVuHqVOnonfv3rC3t5f+rg4cOGCEFpMlSk1NxbJly/DYY48hICAAjo6OcHJyQkBAAKKiovDrr79CCKFzffHx8ZgwYQL8/f3h4OCA5s2bY8SIEdi4cWO155aWliIxMRFffvklXnjhBXTt2hU2NjbSfXr58mWdYrh//z5iYmLw4Ycf4rHHHoOPj49Ux6BBg3Ruiy5KS0sRExODOXPmYMCAAWjatClsbW3h5uaGnj174vXXX8fFixf1qu/LL7/Ev/71LzRt2hSOjo5o27YtpkyZgpSUlGrPN0bbjd0mskCC6B/79+8XAPT6zJw509Rhm0xJSYlO35Gtra1YtGiR1rrS09NFly5dqq3r1VdfFWVlZRrrWLdundZz9+/fXwvfApm7Z599Vqf7NCwsTNy7d6/a+v79738LKyurKuuJiIgQBQUFVZ6/YMECrXGkpaXp1C5/f/8q6xg4cKCO3071srKyhKenZ7Xfn52dnVi+fHm19d2+fVv06dOnynrs7e3FmjVrtNZR07Ybu01kmWx0SxOpoXn55ZcxderUass1adKkDqIxX66urhg0aBD69OmD1q1bo3nz5nByckJGRgYOHDiAb775BtnZ2Zg7dy7c3NwwZcqUSnWUlpYiPDwcp0+fBgB07doVs2fPRocOHfDgwQMcOnQIn3zyCfLy8rBs2TI0b94cc+bMqVSPUOnBsbW1RefOnVFaWirVSw3T9evXAQAeHh544oknMGjQIPj7+8PGxgaJiYlYunQpUlNTsWfPHkRERODgwYOwstL8cOjrr7/G/PnzAQBt2rTB22+/jS5duiAjIwMrVqzA/v37sXPnTrzwwgv47rvvNNahep86ODige/fuuHXrlt69Tar1eHl5oXfv3vjll1/0qkMXRUVFuHPnDgCge/fuiIyMRJ8+feDl5YXs7Gzs2rULK1euRGFhIWbNmgVHR0e89NJLGuuSy+UYM2YMjh07BgAYM2YMXnzxRXh4eODYsWP44IMPkJWVhZdeegm+vr4ICwvTWE9N227MNpEFM23+SeZEtQdw/vz5pg7HIpSWlmo9funSJeHu7i4AiKZNm2osv3nzZul7DwkJ0Vjm5MmTwtbWVgAQ7u7uoqSkpFKZY8eOiRUrVoj4+HipB2b+/PnsAWzgJk6cKL788ktRWFio8XheXp7o37+/dJ+sX79eY7l79+4JNzc3AUC0bNlS3Lp1S+14aWmpiIiIkOo5ePCgxnp2794tvvjiC5GQkCDdxxMnTtS7B/Djjz8WmzdvFleuXJH2QY9eMF1du3ZNDB8+XMTHx1dZ5ujRo8LR0VEAEK6uriInJ0djubVr10oxTp06tdLxCxcuCBcXFwFAtGvXTuP/50LUvO3GbBNZLiaAJGECWDumTJkifa9nzpypdPzVV1+Vju/YsaPKekaPHi2VO336tE7XZgJIujh9+rR0nzz66KMay3z00UdSmQ0bNmgsc/XqVWFtbS0AiFGjRul8fUMSQE1qIwHU1WuvvSZdf+vWrRrLBAYGSr/E5eXlaSyzaNEiqZ7NmzfrfP3aaLsubSLLxUkgVGuKi4uxevVqhIeHw9fXF/b29mjWrBmCg4Mxbdo0xMXFVTnwPDc3F4sXL0ZISAg8PDxgb2+PFi1a4Iknnqj2UcegQYPUBkNfv34ds2fPRtu2beHo6AhPT0+EhYVh165dxm6yRo0aNZK2CwsLKx0vLi6Wtlu3bl1lPW3atJG2i4qKjBSdYe7fv48PP/wQISEhcHd3h62tLZo2bYrAwECMHj0aq1atQlZWVpXnm8u9ceHCBUybNg3t2rWDk5NTlZMQrl27hrlz56Jnz55wd3eHg4MDWrZsiaioKOzfv1/rNeVyOdatW4ewsDB4e3vDzs4Obm5uaNeuHYYOHYro6GidBv7Xps6dO0vDOap6FLt9+3YAgIuLC8aMGaOxTIsWLTBs2DAAQExMDHJzc2shWvM0ePBgaVvTd3jhwgXp7zkqKgpOTk4a63nuueek7a1btxo3SD1V1yaycKbOQMl8GLMHMDExUQQEBFQ7yFjTb/unTp0SPj4+Ws8bM2ZMlQPNBw4cKP0mHBcXp3Ww88cff1yjdlYnPz9ftG3bVgAQVlZWGh+jrFixQq8eQJlMJrKzs3W6fm30AKakpFT79wNArFy5UuP55nJvbN++XTRq1Kja63799dfSo7CqPs8//7zGx3UPHjwQAwYMqLatjz/+uMZ4W7VqJZWpbcpHj126dKl0rKioSNjY2AhAMVlEm+joaCnmP/74Q6dr14cewC1btkjX/+STTyodX7NmTbU9qErt27cXgOJRu65qo+3VtYksG3sAyehSUlIwYMAApKWlAQBGjx6NTZs24cSJEzh69Ci+/fZbPPPMM2o9Y0rXr1/H0KFDkZGRAZlMhkmTJmHPnj04efIk1q9fj27dugFQ/GY8ceJErXHcuHEDo0ePhrW1NRYvXoxDhw7h+PHjWLp0Kdzc3AAAc+fOxdmzZ43a/pKSEqSnp2Pjxo3o27cv/v77bwDApEmT0Lhx40rlx48fDxcXFwDAkiVLIJfLK5VJTEzEr7/+CgAYN26cVN4UJkyYgIyMDNja2mLq1KnYuXMnTpw4gWPHjmHbtm2YO3cu2rdvr/Fcc7k30tPT8cwzz8DJyQmLFy/G4cOHcfToUaxcuRLOzs5SuW+++QYvvPACCgoK0LlzZ6xcuRKHDh3CqVOnsGXLFjzyyCMAgDVr1uDNN9+sdJ0FCxYgLi4OADBq1Chs2LABhw8fRkJCAnbv3o0lS5ZgwIABkMlkOnzztScxMRE5OTkAgI4dO1Y6fuHCBZSWllZ5XJXq8XPnzhkxSvN28OBBaVvTd6T6Xej6HV69ehV5eXlGilB/1bWJLJypM1AyH6o9gC+//LI4ffp0tZ/i4uJK9fTo0UPq8dL2m+7t27dFfn6+2r4nnnhCiuHrr7+udE5hYaEYPHiwVOa3336rVEbZywNAtGrVSly7dq1Smbi4OCGTyQQAMWPGDF2+Hq3S0tK09vAMGzZM3L9/v8rzt27dKvUy9ejRQ3z77bciPj5exMTEiAULFojGjRsLAKJ79+7ixo0bOsdl7B7AixcvVtvDJ4QQZWVl4u7du5X2m9O94ePjozaIvqL09HTh5OQkAIiJEydWOSD/7bffltqUmpqqdszPz08AEE888USV1xFCiDt37mjcX1c9gKrf7U8//VTp+K5du3TuNT9x4oRU9q233tLp+pbeA5iRkSH9P9qkSZNK964QQkRFRUnxVZxAU9Err7wilf3rr790isHYbdelTWTZmACSxJB1ACv+Y717927pmL5rBGZkZEgDyLU9ZkpLS5MeRz3yyCOVjqv+kNf2SPXhhx+WEq6aqioB9PT0FBs2bKh2trAQQpw9e1Y8//zzUmKq+vHy8hJLly4Vubm5esVl7ATw8OHDUn3Jycl6nWtu90ZVs12VlAPgfXx8tK5rV1JSInx9fQUAMW/ePLVjypnbK1asqKaFmtVFAqg6Cz04OFjjOpM//vijVGbVqlVa60tJSZHKTps2TacYLDkBLCsrE48++qh07ar+rh955BGpjLb7SQgh3njjDansyZMndYrDmG3XtU1k2fgImIxK+ZgSAF599VW9zt2/f7/0+PP555+vspy/vz+GDx8OQPH2Ek2PTAHAzc0N4eHhVdYTHBwMALh06ZJecWri6+uL06dP4/Tp00hMTMQvv/yCadOmIT8/H1OnTsWSJUu0vmmhpKQEP/zwA3bu3Kmx3M2bN7FhwwaTv82jefPm0va6dev0Otec7g07Ozs8+eSTWq/5888/AwAiIiLg4OBQZTkbGxuEhIQAULwhQ5Xy+9q0aRPy8/O1Xk+Ty5cvQyh+Udf7XF389ddfmDRpEgDA0dER69ev1/g4WnXykp2dndY67e3tpe2CggIjRWq+oqOjsWPHDgCKSRPTpk3TWM6SvkNd20SWjQkgaTR//nzpB4+2j7+/v9p5iYmJAICWLVuiVatWel3zzJkz0nafPn20llUez8/PrzKBa9euXZUL2gKKhXEB4MGDB3rFqYly4eXOnTuje/fuCA8Px8qVK3H06FHIZDLMmzevysQlLy8Pw4YNw4cffog7d+7gjTfewLlz51BUVITs7Gzs3bsX/fv3x4kTJxAREYEVK1bUOF5DBQQEYMCAAQCAZcuWISgoCO+99x7++OOPahMcc7s3tCV12dnZ0tjNL7/8stpXIm7evBkAkJmZqVaPcizikSNHEBAQgGnTpmHbtm24detWNS2ufRkZGRg5ciQePHgAmUyGNWvWIDAwUGNZ1e9Kdda6Jqoz1B0dHY0TbC04f/48zpw5o/Fz//59ner4/vvv8e677wJQ/PLxww8/VPlvjqV8h/q0iSwb/1bJqG7fvg1AvadIV3fv3pW2vby8tJb19vbWeJ6qqpZZUFL+o1ZWVqZriHrr2rUrPvjgAwDA2rVrsXfv3kpl5s+fj9jYWACKyQRLlixBx44dYWdnBxcXFwwfPhz79+/H4MGDIYTA7Nmz8eeff9ZazNXZsGGD1OOVkpKChQsXYujQoXBzc8PAgQPxxRdfaFzuxpzuDXd3d611aFvCRpuKSfC7776LyZMnQyaTISsrC5999hnGjBkDLy8vdOnSBfPnz8fNmzcNulZN3L17F6GhodKSNytWrMD48eOrLK86eam6pV1UJy2oTqgxN6GhoejSpYvGj3LJG21+/fVXTJo0CUIIeHl5ISYmRu3eq8gSvkN920SWjQkg1YrantVYW4/EakNkZKS0rewpUhJCYO3atQCA9u3bVzl71cbGBgsXLgSgSFiV55iCr68vjhw5gn379mHq1KkICgqCTCZDSUkJYmNj8fLLL6Nz5844f/68xvPN4d6wtrbWelz10fGsWbOkx/vVffbs2aNWj62tLdasWYMzZ87gnXfeQd++fWFnZwchBM6cOYP3338fbdu2lR4314UHDx5gxIgR0uz3hQsXYvr06VrPadGihbR97do1rWWvXr0qbfv5+dUgUvN14MABPPHEEygpKYG7uzv27t2Ltm3baj3HkO9QJpOpnVebDGkTWTa+C5iMSrmYbEZGht7nKh/JAooxby1btqyyrGqviep55qhp06bS9pUrV9SO3bx5U+ql6tGjh9Z6lGMWAcXYLVMbOnQohg4dCgC4c+cO9u3bh9WrV+OPP/7AxYsXERUVJT32BSzr3vD09JS28/Pz0blzZ4PqUQoMDMTChQuxcOFCFBQU4PDhw/jhhx+wfv165ObmYvz48bh48aJBvaP6KCgoQEREBE6cOAEAmDNnDt55551qz2vfvj2sra0hl8urvfdUj3fq1KlmAdciTQt+6+L48eOIiIhAYWEhnJ2dsWvXLnTt2rXa81Qfr//111/o3r17lWWV36Gfn5/GJZGMzdA2kWVjDyAZVc+ePQEo1lmrmOxUR/WHrPJl6VU5fvw4AMVj3oCAAD2jrFvXr1+Xtis+zrGxKf8dTLnOWlVKSko0nmcOPD09ERUVhd9//x2PPvooACApKQkXLlyQyljSvdG0aVP4+voCAPbt22fUHmdHR0cMGzYM33zzDT7++GMAisSsureY1FRJSQkef/xxaW23KVOm4KOPPtLpXDs7Ozz00EMAFBNdtI1hU9Zvb2+PXr161TBq8/Lnn39ixIgRyM3NhYODA3bu3FntmFSl/v37S9uq6+tVlJmZKfWe9+vXr2YB66AmbSLLxgSQjCoiIkLaXrZsmV7nDho0SHo0t2bNmirLpaenIyYmRjrH3JKhin766Sdpu0uXLmrHPDw8pEWd4+PjtSaBqj80zDnpVfYKAuXj/gDLuzeUieylS5cqPbo3lqq+K2OTy+V46qmnpNcfTpgwAZ9//rledTz22GMAgJycnCpfUXbt2jXs27cPgKJtmhY+t1Tnz59HaGgo7t27B1tbW2zZskV6paAu2rdvL/WI/vjjj1VOmlKdXT969OiahFytmraJLFydLjpDZs1Yr4ILDg7WabHfO3fuaF3sd82aNZXOKSoqEkOHDtVpsd/q1sNSXSPPENu2bRMZGRlayxw8eFA4OzsLAMLGxkbjoq7jx4+X4liwYIHGeu7evSu9SB6A2LNnj04xGnsdwMTERJGYmFjl8bKyMhERESEAxSvrKn4/lnJvCCHEpUuXhL29vQAgPDw8xIkTJ7SW//XXX9XWRrxz5474+eefNa6rp/Txxx9L8Wr6PoyxDmBZWZmYNGmSVM/jjz+u07qUFd25c0e4uroKQLHA+u3bt9WOl5aWSn/3gO6vgRPC/NcBvHLlirSot7W1tcbFsnWh+jq4V155pdLxv//+W3olX5s2bapcfFwTfdturDaR5TLvrhMymaysLLWlN6ri6OiINm3aqO373//+h4ceekga2/TTTz9h3LhxaN26NeRyOf7++2/ExMRg8+bNOH36tNpSMsuWLcPvv/+Oe/fu4YUXXsDhw4cxbtw4eHh44K+//sJ//vMfJCUlAQDGjh2LkSNHGrXd+ti+fTuioqIQHh6OoUOHIigoCG5ubigqKsLFixexc+dO/Pjjj9Is43fffRcdOnSoVM97772Hn3/+Gfn5+ViwYAESEhIwceJEtG7dGoWFhTh69CiWL1+O9PR0AIqeldDQUI0xVVybT/ldAcDu3bvVxj21bdtW7bGULpKSkjBp0iT07t0bERER6NmzJ7y9vVFSUoK0tDSsXbtW6oGLjIysNKbNku6NgIAAfPHFF5g0aRLu3r2Lfv36YcKECRg1ahRatmyJ0tJSXLt2DcePH8fmzZulv3Pl2KmcnBxERkbC398fY8aMQZ8+fdCqVSvY2Njgxo0b2LlzJ77++msAigkCqj2kxvT6669Lk4Y6d+6Mt99+u9pXtGka8+jh4YElS5ZgypQpuHLlCvr06YN58+ahS5cuyMjIwPLly7F//34AitcbDh48WGPdubm5lXpUlUvuAIqJUsrxogDQvXt3jePlkpKS1O5vVZmZmZX+X3jiiScMmlF7584dDBs2TJqY8dprr6Fjx45a/310d3eXhhComjhxIr755hscPnwYn332GTIzM/Hiiy/C3d0dx48fx8KFC5GTkwMrKyusXLmyyh7smrbdmG0iC2bqDJTMhyFvAunWrZvGuk6ePCn9dqnto+m3/VOnTgkfHx+t540ZM6bK1fTrqgdQtddC28fR0VH85z//0VpXTEyMaNKkSbV1DRkyROMr1pT0+bubOHGi3m1eu3atTnX379+/ytebWcK9oWrjxo1Sr4y2j5WVlVqvV3WvB1R+fH19xalTpzRe2xg9gKp16PrR5r333tP4thrl55FHHtH6pgtdvxflp6qnEar//+ryMbRn0ZB/F7X9v3Xr1i3Ru3fvKs+1s7MTq1ev1hpTTdtu7DaRZWIPINWK4OBgpKam4uuvv8b27dtx5swZ3Lt3D56envD19UX//v0xbty4SgtJA4rZsKmpqfj000+xfft2pKamIj8/H02aNMHDDz+M5557rtZ6S/Txn//8B4888gj++OMPnDp1CpmZmcjKyoKVlRU8PDwQFBSEIUOG4Nlnn612duewYcPw119/Yc2aNdi1axfOnj2L+/fvw8bGBt7e3ujduzeeeuopPProo7W+jIo2Tz31FPz9/RETE4O4uDhcu3YNN2/eRGlpKZo1a4aePXti3LhxiIqKqnLxWEu7N6KiohAaGorVq1dj9+7dSElJkcZMeXt7IygoCIMHD8YTTzyhtuxJq1atkJSUhJiYGPzxxx+4dOkSbt68idzcXLi5uSEoKAgRERF46aWXLGqs3L///W+EhYXhs88+Q1xcHG7evAk3Nzd069YNkyZN0rqeIClmwx85cgRfffUVfvjhB5w7dw55eXnw8fHB0KFDMXPmTAQFBZk6TGoAZEJY0IJqRERERFRjnAVMRERE1MAwASQiIiJqYJgAEhERETUwTACJiIiIGhgmgEREREQNDBNAIiIiogaGCSARERFRA8OFoGugrKwMGRkZaNy4sUkX5yUiIiLdCSHw4MED+Pj4VLlofX3HBLAGMjIy1Fb+JyIiIstx9epVtGjRwtRhmAQTwBpQvr7p6tWrcHFxMXE0REREpIucnBz4+flZ1GsYjY0JYA0oH/u6uLgwASQiIrIwDXn4VsN88E1ERETUgDEBJCIiImpgmAASERERNTAcA1jLhBAoLS2FXC43dShEAABbW1tYW1ubOgwiIjIhJoC1qLi4GDdu3EB+fr6pQyGSyGQytGjRAs7OzqYOhYiITIQJYC0pKytDWloarK2t4ePjAzs7uwY924jMgxACt27dwrVr19CuXTv2BBIRNVBMAGtJcXExysrK4OfnBycnJ1OHQyRp2rQpLl++jJKSEiaAREQNFCeB1LKG+ooZMl/siSYiImYnRERERA0ME0AiIiKiBoYJIBEREVEDwwSQKomNjUVERAR8fHwgk8mwfft2U4cEAFi1ahW6du0qvXs5JCQEu3btqlTu888/R0BAABwcHBAcHIy4uDiDytQmY7aFiIhIX0wAqZK8vDx069YNn376qalDUdOiRQssXrwYJ0+exMmTJzFkyBBERkbi7NmzUplNmzZh1qxZmDdvHhITEzFgwACMHDkS6enpepWxlLaQBXrwADh7FsjONnUkRNSQCTJYdna2ACCys7MrHSsoKBApKSmioKDABJEZDwCxbds2vc+Li4sTNjY2orCwUNp36dIlAUBcvnzZaPG5u7uLr7/+WvrzQw89JKZMmaJWpmPHjuKtt97Sq4wm8fHxYsiQIcLT01MAUPvcu3fPJG0xRH25Ny1OUZEQs2cL4eAgBCCEra0QL74oxIMHpo6MqMHR9vO7oWAPIOlt3bp11S4lkpSUhE6dOsHe3l5tn5ubG1q1agUAiI6OhrOzs9ZPVY885XI5Nm7ciLy8PISEhABQrL2YkJCA0NBQtbKhoaE4cuSIzmU0SU5OxqBBg9CtWzfExsZi9+7d8PDwwODBg7Fp0ya4ubkZ3B5D20IWpKQEePRRYOlSoLCwfN9XXwHDhgF5eaaNj4gaHC4EXdd69QIyM+v2mt7ewMmTRqvO1dUVHTp00FomOTkZPXr0UNuXlJSEbt26SX+eMmUKxo4dq7UeX19ftT+fPn0aISEhKCwshLOzM7Zt24bAwEAAwO3btyGXy+Hl5aV2jpeXFzL/+c51KaPJjBkzEBkZiaVLlwIAAgMDMX78eBw7dkxqg77tqWlbyIK8/TawZ49i294eCAsDfv9dkfgdOwZMnQp8+61pYySiBoUJYF3LzASuXzd1FDUyevRojB49WmuZpKQkPPXUU2r7EhMT1RJADw8PeHh46HXtDh06ICkpCffv38eWLVswceJEHDx4UEqcgMoLHQshKu3TpYzSzZs3cejQIfzxxx9q+xs1aqR2jr7tMVZbyMwlJCh6/gDA1laRCA4cCCQnAwMGKMYErl8PjB8PjBhh2liJqMFgAljXvL3r/TXlcjnOnj1bqQfw1KlTaoljdHQ0oqOjtda1a9cuDBgwQPqznZ0d2rZtCwDo1asXTpw4gRUrVuDLL79EkyZNYG1tXamHLCsrS+pJ06VMRQkJCSgrK1NLXpX7e/XqZXB7atoWshBvvAGUlSm2Fy5UJH8A0K0b8N//ApMmlZcbPhzg6/mIqA4wAaxrRnwUa65SU1NRUFAAHx8faV98fDyuX79e40fAFQkhUFRUBECRUAUHByMmJkYt0YyJiUFkZKTOZSoq++eHd0FBAdzc3AAoHt/Gxsbi/fffN1p79G0LWYC4OEDZc9y2LTB7tvrxiROBzz8HTpwATp8Gtm8HHn+87uMkoobHpFNQLFx9nQX84MEDkZiYKBITEwUAsXTpUpGYmCiuXLkihBBi69atokOHDlWe//333wsAYtq0aeL8+fPit99+E23bthUAxJEjRwyOa+7cuSI2NlakpaWJP//8U7z99tvCyspK7N27VyqzceNGYWtrK9asWSNSUlLErFmzRKNGjdRmHutSRlVWVpZwdHQUzzzzjDh37pz45ZdfREBAgJg+fbrJ22IIS743Lc6YMYoZv4AQ69ZpLrN7d3mZvn3rNj6iBoqzgIVgAlgD9TUB3L9/f6VlTgCIiRMnCiGEWLt2rdD2u8OcOXPE8OHDRXh4uLCzsxPdu3cXmzdvFi4uLuLpp582OK7JkyeLVq1aCTs7O9G0aVMxdOhQtYRJ6bPPPpPK9ezZUxw8eFCvMprat3PnTtG+fXtha2srWrduLZYsWSLkcrlZtEVflnxvWpT0dCGsrBSJXfPmimVgNCkrEyIoqDwJPHOmbuMkaoCYAAohE0IIU/Q81gc5OTlwdXVFdnY2XFxc1I4VFhYiLS1NeotDQxIWFoaePXti0aJFpg7FIAsWLMCBAwcbOuXxAAAgAElEQVRw4MABU4dSKxryvVmnPv5YMa4PAObPBxYsqLrsypXAjBmK7blzgWrGkhJRzWj7+d1QcB1AMrrk5GR07drV1GEYbM+ePfjoo49MHQZZuo0by7efeUZ72aio8skf339fPmmEiKiWMAEko8rMzMTNmzctOgGMj4/HQw89ZOowyJJduACcOqXYDg5WTADRplkzQLnod3o6cPhw7cZHRA0eE0AyKm9vbwghEBQUZOpQiExn06by7XHjdDtHtZfw+++NGw8RUQVMAImIjO3HH8u3q1kaSBIZCSjHZP7yi2JKCBFRLWECSERkTNeuKdb0A4CHHgJattTtvEaNgEGDFNvXrwNnztRKeEREABNAIiLj2ru3fFvfV7uNHFm+vWuXceIhItKACSARkTGpJoBhYfqdywSQiOoIE0AiImORy4GYGMW2q6viEbA+2rUD2rRRbB86BOTkGDc+IqJ/MAEkIjKWhATg7l3F9tChgI0Br1tX9gKWlpa/R5iIyMiYABIRGYuy9w8oX9dPX8OHl2/HxtYsHiKiKjABJCIylri48u0hQwyro18/zfURERkRE0AiImOQy4H4eMW2l1f1b/+oiqcnoFxIPTERyM01TnxERCqYABIRGcOZM+WTNvr3B2Qyw+saMEDxX9WkkojIiJgAUiWLFi1C79690bhxYzRr1gyPPfYYUlNTTR2WmkWLFkEmk2HWrFmVjn3++ecICAiAg4MDgoODEafhMZouZWqTv78/ZDJZpc8rr7xiVnGSHg4dKt/u379mdSkTQICPgYmoVjABpEoOHjyIV155BUePHkVMTAxKS0sRGhqKvLw8U4cGADhx4gRWr16Nrl27Vjq2adMmzJo1C/PmzUNiYiIGDBiAkSNHIj09Xa8yddGGGzduSJ+YfyYPPPnkk2YVJ+nh8OHybdVxfIZQTQA5EYSIaoMgg2VnZwsAIjs7u9KxgoICkZKSIgoKCkwQmXFlZWUJAOLgwYM6nxMXFydsbGxEYWGhtO/SpUsCgLh8+bLBsTx48EC0a9dOxMTEiIEDB4qZM2eqHX/ooYfElClT1PZ17NhRvPXWW3qV0SQ+Pl4MGTJEeHp6CgBqn3v37hncJiGEmDlzpmjTpo0oKyurcZzVqU/3plnx8xMCEMLJSYji4prX5++vqM/BwTj1EZFE28/vhoI9gFSt7OxsAICHhwcAYN26dZBVM74pKSkJnTp1gr29vdo+Nzc3tGrVCgAQHR0NZ2dnrZ+KjzxfeeUVhIeHY9iwYZWuWVxcjISEBIRWWH4jNDQUR44c0bmMJsnJyRg0aBC6deuG2NhY7N69Gx4eHhg8eDA2bdoENzc3g9qjjOm7777D5MmTpe/V0DjJRK5eVXwA4OGHAVvbmtf58MOK/xYWAmfP1rw+IiIVBqxSSjXRa3UvZOZm1uk1vZ29cfKlkwadK4TA7Nmz0b9/f3Tu3BkA4Orqig4dOmg9Lzk5GT169FDbl5SUhG7dukl/njJlCsaOHau1Hl9fX2l748aNOHXqFE6cOKGx7O3btyGXy+Hl5aW238vLC5mZmTqX0WTGjBmIjIzE0qVLAQCBgYEYP348jh07JrVB3/Yobd++Hffv38dzzz2nV1vIjKjek8rEraZ69wY2biyvv3t349RLRAQmgHUuMzcT1x9cN3UYOps2bRr+/PNPHFIZ4D569GiMHj1a63lJSUl46qmn1PYlJiaqJYAeHh5Sr2J1rl69ipkzZ2Lv3r1wcHDQWrZi76QQotI+Xcoo3bx5E4cOHcIfFd7K0KhRI7Vz9GmPqjVr1mDkyJHw8fGpdEyfOMmETqr8gtW7t3HqVH2N3PHjwIsvGqdeIiIwAaxz3s7eFnPN6dOnY8eOHYiNjUWLFi10Pk8ul+Ps2bOVegBPnTqlljhGR0cjOjpaa127du3CgAEDkJCQgKysLAQHB6tdJzY2Fp9++imKiorQpEkTWFtbV+ohy8rKknrSdClTUUJCAsrKytSSV+X+Xr16GdQepStXrmDfvn3YunWrWjlD4iQTUk0AVe6JGunRA7CyAsrK1HsYiYiMgAlgHTP0UWxdEkJg+vTp2LZtGw4cOICAgAC9zk9NTUVBQYFaj1Z8fDyuX79u8CPgoUOH4vTp02rHJk2ahI4dO+LNN9+EtbU1rK2tERwcjJiYGLVEMyYmBpGRkQAAOzu7astUVFZWBgAoKCiAm5sbAOD06dOIjY3F+++/b1B7lNauXYtmzZohPDxcbb8hcZKJCFGeAHp5ARoe8xukUSPFgtCnTyvWGMzPB5ycjFM3EZFJp6BYuPo6C/jll18Wrq6u4sCBA+LGjRvSJz8/XwghxNatW0WHDh2qPP/7778XAMS0adPE+fPnxW+//Sbatm0rAIgjR44YLU5Ns4A3btwobG1txZo1a0RKSoqYNWuWaNSokdrMY13KqMrKyhKOjo7imWeeEefOnRO//PKLCAgIENOnT69R/HK5XLRs2VK8+eabGo/rG6euLPneNEsXLypm6wJChIcbt+7Jk8vrPnTIuHUTNWCcBSyE2SSAn332mfD39xf29vaiZ8+eIjY2Vmv5zZs3i06dOgk7OzvRqVMnsXXrVrXjW7ZsEaGhodKSHYmJiWrH79y5I6ZNmybat28vHB0dhZ+fn5g+fbq4f/++zjHX1wQQFZY4UX7Wrl0rhBBi7dq1QtvvDnPmzBHDhw8X4eHhws7OTnTv3l1s3rxZuLi4iKefftpocWpKAIVQ3EutWrUSdnZ2omfPnhqXr9FWRlP7du7cKdq3by9sbW1F69atxZIlS4RcLq9R/Hv27BEARGpqapVldGmLviz53jRLmzaVJ2kLFhi37i++KK972TLj1k3UgDEBNJMEUNnT8dVXX4mUlBQxc+ZM0ahRI3HlyhWN5Y8cOSKsra1FdHS0OHfunIiOjhY2Njbi6NGjUpn169eLf//73+Krr77SmACePn1ajBkzRuzYsUP8/fff4vfffxft2rUTjz/+uM5x19cEsKZCQ0NrvFadKc2fP18MHDjQ1GHUmoZ8b9aKOXPKk7RffjFu3QkJ5XU/84xx6yZqwJgACiETQog6fupcSZ8+fdCzZ0+sWrVK2tepUyc89thjWLRoUaXyUVFRyMnJwa5du6R9I0aMgLu7OzZs2KBW9vLlywgICEBiYiK6V7OMwk8//YRnnnkGeXl5sLGpfnhkTk4OXF1dkZ2dDRcXF7VjhYWFSEtLk17j1ZB4e3tj2bJlGD9+vKlDMUhISAhWrFiBh1RnYdYjDfnerBVDhgD79yu2b9wAvI040auoCHB2BkpLgS5dgD//NF7dRA2Ytp/fDYXJF4I2ZMHb+Pj4SuXDwsJqvECu8kbQJfkjzTIzM3Hz5k2Nr2mzFPHx8fU2+SMjE6I8KWve3LjJHwDY2wOdOim2z51TJIREREZg8kzHkAVvMzMzjb5A7p07d7Bw4UL83//9X5VlioqKUKTyD3BOTo7B16uvvL29YQadykR148YN4M4dxXZt/dLTrZtiJnBpKZCSolgehoiohkzeA6ik74K3xlwgNycnB+Hh4QgMDMT8+fOrLLdo0SK4urpKHz8/P4OuR0T1hOrSRF261M41VIeuJCXVzjWIqMExeQJoyIK33t7eRlsg98GDBxgxYgScnZ2xbds22Gp5h+fcuXORnZ0tfa4q3/1JRA2T6pi82uwBVEpOrp1rEFGDY/IEUHXBW1UxMTHo27evxnNCQkIqld+7d2+V5auSk5OD0NBQ2NnZYceOHdUOiLe3t4eLi4vah4gaMNUEsLZ6AJkAElEtMPkYQACYPXs2JkyYgF69eiEkJASrV69Geno6pkyZAgB49tln4evrK80InjlzJv71r39hyZIliIyMxM8//4x9+/apva/27t27SE9PR0ZGBgDF2ykARe+ht7c3Hjx4gNDQUOTn5+O7775DTk6ONKavadOmsLa2rsuvgIgskfIRsLV1+WQNY2vaFPDxATIyFI+AhQD4PmgiqiGzSACjoqJw584dvP/++7hx4wY6d+6M3377Da1atQIApKenw8qqvLOyb9++2LhxI9555x28++67aNOmDTZt2oQ+ffpIZXbs2IFJkyZJfx43bhwAYP78+ViwYAESEhJw7NgxAEDbtm3V4klLS4O/v39tNZeI6oOSEsWkDADo2FExY7e2dO+uSADv3wfS04F//m0kIjKUWawDaKm4DiBZIt6bRnL2LNC5s2J73DigwhqkRjV3LrB4sWL7l1+ACu+OJiL9cB1AMxgDSERkkepiAohSUFD5trLXkYioBpgAEhEZoi4mgCgxASQiI2MCSERkCNU1AGu7B7BDh/KJH2fP1u61iKhBYAJIRGQIZQ+gqytQ24vCOzkBrVsrtlNSFDOBiYhqgAkgVbJq1Sp07dpVWuswJCQEu3btMnVYWLBgAWQymdrHu8K7V3UpU1painfeeQcBAQFwdHRE69at8f7776OsrKzO2hIbG4uIiAj4+PhAJpNh+/btlcr4+/tXaotMJsMrr7yiVu769et45pln4OnpCScnJ3Tv3h0JCQl11ZSG6f59QLkQfJcudbMsS2Cg4r95eYqZwERENWAWy8CQeWnRogUWL14sLY/z7bffIjIyEomJiQhSHYtkAkFBQdi3b5/0Z03rNVZXZsmSJfjiiy/w7bffIigoCCdPnsSkSZPg6uqKmTNn1l7wKvLy8tCtWzdMmjQJjz/+uMYyJ06cgFwul/585swZDB8+HE8++aS07969e+jXrx8GDx6MXbt2oVmzZrh48SLc3NxqvQ0N2pkz5du1Pf5PKTAQ2LlTsZ2SwqVgiKhGmABSJREREWp//vDDD7Fq1SocPXpU5wTw0KFDGDx4MHJzc2H/z/poaWlpaN26NS5fviyt8agvGxubSj16+paJj49HZGQkwv9ZSsPf3x8bNmzAyZMntdZ79OhRzJs3D8nJybhz547asXv37umVdI0cORIjR47UWqZp06Zqf168eDHatGmDgQMHSvuWLFkCPz8/rF27VtrHNSzrwF9/lW8re+Zqm+r/e2fPAtXcP0RE2vARMGkll8uxceNG5OXlISQkBACwbt06yKp55JWUlIROnTpJyZ9yn5ubm5T8RUdHw9nZWesnLi5Ord4LFy7Ax8cHAQEBGDduHC5dulTp2tWV6d+/P37//XecP38eAJCcnIxDhw7hkUceqbI9ycnJGDRoELp164bY2Fjs3r0bHh4eGDx4MDZt2gQ3NzeD2qOr4uJifPfdd5g8ebLad79jxw706tULTz75JJo1a4YePXrgq6++MugapIdz58q3O3asm2uqJpqcCUxENcQewDrWq1cvZGZm1uk1vb29q+3dquj06dMICQlBYWEhnJ2dsW3bNgT+8wPI1dUVHTp00Hp+cnIyevToobYvKSkJ3VTeazplyhSMHTtWaz2+vr7Sdp8+fbB+/Xq0b98eN2/exAcffIC+ffvi7Nmz8PT01LnMm2++iezsbHTs2BHW1taQy+X48MMPMX78+CrjmDFjBiIjI7F06VIAQGBgIMaPH49jx45JbdC3PfrYvn077t+/j+eee05t/6VLl7Bq1SrMnj0bb7/9No4fP44ZM2bA3t4ezz77rEHXIh2o9gDW1ivgKurUSTHWUAjOBCaimhNksOzsbAFAZGdnVzpWUFAgUlJSREFBgdp+X19fAaBOP76+vnq3raioSFy4cEGcOHFCvPXWW6JJkybi7NmzOp/fq1cvsXTpUrV9ERERYsaMGXrHUpXc3Fzh5eUlPvnkE73KbNiwQbRo0UJs2LBB/Pnnn2L9+vXCw8NDrFu3TmMdmZmZwsrKShw4cEBt/xtvvCF69+5d43YAENu2bdNaJjQ0VIwaNarSfltbWxESEqK2b/r06eLhhx+usq6q7k3SQ+vWQgBCODsLUVZWd9cNCDDNdYnqGW0/vxsK9gDWserGr5nLNe3s7KRJIL169cKJEyewYsUKfPnll9WeK5fLcfbs2Uo9gKdOncLo0aOlP0dHRyM6OlprXbt27cKAAQM0HmvUqBG6dOmCCxcuVHm+pjJz5szBW2+9Jb0fukuXLrhy5QoWLVqEiRMnVqojISEBZWVlar2Xyv29evUyWnuqcuXKFezbtw9bt26tdKx58+ZSz6xSp06dsGXLFr2uQXooKADS0hTbHTvWzQxgpaAgxbVzc4Fr12p/+RkiqreYANYxfR/FmgshBIqKinQqm5qaioKCAvj4+Ej74uPjcf369Ro9Aq6oqKgI586d05pQaSqTn58PKyv14a/W1tZVLgOj3F9QUCBN9Dh9+jRiY2Px/vvvG609VVm7di2aNWsmTVpR1a9fP6SmpqrtO3/+vMGTbEgHFy6Ur8NXV49/lQIDFe8CBhSPgZkAEpGBmABSJW+//TZGjhwJPz8/PHjwABs3bsSBAwewe/duAMC2bdswd+5c/KU6DkpFUlISAGDlypWYMWMG/v77b8yYMQMA1JJIDw8PeHh46BzX66+/joiICLRs2RJZWVn44IMPkJOTo9Zrp0uZiIgIfPjhh2jZsiWCgoKQmJiIpUuXYvLkyRqv26dPHzg6OuKNN97AvHnzcPHiRUyfPh1TpkxB3759DWpPbm4u/v77b+nPaWlpSEpKgoeHB1q2bCntLysrw9q1azFx4kTY2FT+3/XVV19F3759ER0djbFjx+L48eNYvXo1Vq9erVMcZADV+76uJoAoVZwIMmJE3V6fiOoPUz+DtmSGjAG0BJMnTxatWrUSdnZ2omnTpmLo0KFi79690vG1a9cKbbfOnDlzxPDhw0V4eLiws7MT3bt3F5s3bxYuLi7i6aefNjiuqKgo0bx5c2Frayt8fHzEmDFjKo1L1KVMTk6OmDlzpmjZsqVwcHAQrVu3FvPmzRNFRUVVtm/nzp2iffv2wtbWVrRu3VosWbJEyOVyg9uyf/9+jeM1J06cqFZuz549AoBITU2tsq6dO3eKzp07C3t7e9GxY0exevVqrde25HvTLCxYoBiHBwixZUvdXvvEifJrT55ct9cmqkc4BlAImRB8p5ChcnJy4OrqiuzsbLi4uKgdKywsRFpaGgICAuDg4GCiCE0jLCwMPXv2xKJFi0wdikEWLFiAAwcO4MCBA6YOpVY05HvTKMaPBzZuVGynpNTtY+DcXKBxY8V2v37AoUN1d22iekTbz++GgusAktElJyeja9eupg7DYHv27MFHH31k6jDIXCnXALS2Btq0qdtrOzsDyrG1/6xjSURkCI4BJKPKzMzEzZs3LToBjI+PN3UIZK7KygDlpJs2bQA7u7qPoX17ICMDuHULuHcPcHev+xiIyOKxB5CMytvbG0IIk78zmKhWpKcDhYWK7bqeAazUvn35tpYlkIiItGECSESkK1O8Aq4i1bfwVFgCiIhIV0wAiYh0ZcolYJRUewA5DpCIDMQEkIhIV0wAiaieYAJIRKQr1YRL9VFsXQoIUMxArhgPEZEemADWMi6zSOaG92QNKCddeHqabvatrS3QurVi+/x5xcxkIiI9MQGsJba2tgAU750lMifFxcUAFO8/Jj3k5wPXryu227UzbSzK3sf8fMWSMEREeuI6gLXE2toabm5uyMrKAgA4OTlBJpOZOCpq6MrKynDr1i04OTlpfLcwaaHy7maTJ4AVxwG2aGG6WIjIIvEnQC3y9vYGACkJJDIHVlZWaNmyJX8h0ZfqmnvmlgAOGWK6WIjIIjEBrEUymQzNmzdHs2bNUFJSYupwiAAAdnZ2sLLi6A+9mXMCSESkJyaAdcDa2prjrYgsnbkmgFwMmogMwG4AIiJdmFMC6OMDNGqk2GYPIBEZgAkgEZEulAlgs2aAi4tpY5HJynsB09KAf2Z2ExHpigkgEVF1cnOBzEzFdtu2po1FSZkAyuXA5csmDYWILA8TQCKi6pjTEjBKbdqUb1+8aLo4iMgiMQEkIqqOOY3/U1JNAFUTVCIiHTABJCKqjrkngOwBJCI9MQEkIqoOE0AiqmeYABIRVUc1ATSXSSA+PoCDg2KbCSAR6YkJIBFRdZQJoLc30LixaWNRsrICWrdWbF+6BJSVmTYeIrIoTACJiLTJyQGU7/M2l8e/SsrHwEVFwPXrpo2FiCwKE0AiIm3S0sq3VcfdmQOOAyQiAzEBJCLS5tKl8u2AANPFoYnqeEQmgESkByaARETaqPYAKsfcmQuuBUhEBmICSESkjWoCaG49gHwETEQGMpsE8PPPP0dAQAAcHBwQHByMuLg4reW3bNmCwMBA2NvbIzAwENu2bVM7vnXrVoSFhaFJkyaQyWRISkqqVEdRURGmT5+OJk2aoFGjRnj00Udx7do1o7aLiCycOSeArVopZgMDTACJSC9mkQBu2rQJs2bNwrx585CYmIgBAwZg5MiRSE9P11g+Pj4eUVFRmDBhApKTkzFhwgSMHTsWx44dk8rk5eWhX79+WLx4cZXXnTVrFrZt24aNGzfi0KFDyM3NxahRoyCXy43eRiKyUMoE0N5esQyMObGzUySBgCIBFMK08RCRxZAJYfp/Mfr06YOePXti1apV0r5OnTrhsccew6JFiyqVj4qKQk5ODnbt2iXtGzFiBNzd3bFhwwa1spcvX0ZAQAASExPRvXt3aX92djaaNm2K//3vf4iKigIAZGRkwM/PD7/99hvCwsKqjTsnJweurq7Izs6Gi4uL3u0mIjMnBNCoEVBQAHToAPz1l6kjqmz4cGDfPsX27duAp6dp4yGyAPz5bQY9gMXFxUhISEBoaKja/tDQUBw5ckTjOfHx8ZXKh4WFVVlek4SEBJSUlKjV4+Pjg86dO1dZT1FREXJyctQ+RFSP3bypSP4A83v8q8SJIERkAJMngLdv34ZcLoeXl5fafi8vL2RmZmo8JzMzU6/yVdVhZ2cHd3d3netZtGgRXF1dpY+fn5/O1yMiC2TOM4CVOBGEiAxg8gRQSSaTqf1ZCFFpX03K60pbPXPnzkV2drb0uXr1ao2vR0RmzJwngCgxASQiA5g8AWzSpAmsra0r9bplZWVV6uVT8vb21qt8VXUUFxfj3r17Otdjb28PFxcXtQ8R1WOWkAByMWgiMoDJE0A7OzsEBwcjJiZGbX9MTAz69u2r8ZyQkJBK5ffu3VtleU2Cg4Nha2urVs+NGzdw5swZveohonrMEhJA1UfTHANIRDqyMXUAADB79mxMmDABvXr1QkhICFavXo309HRMmTIFAPDss8/C19dXmhE8c+ZM/Otf/8KSJUsQGRmJn3/+Gfv27cOhQ4ekOu/evYv09HRkZGQAAFJTUwEoev68vb3h6uqK559/Hq+99ho8PT3h4eGB119/HV26dMGwYcPq+BsgIrNkCQmgszPg5aWYsMIeQCLSkVkkgFFRUbhz5w7ef/993LhxA507d8Zvv/2GVv+sb5Weng4rq/LOyr59+2Ljxo1455138O6776JNmzbYtGkT+vTpI5XZsWMHJk2aJP153LhxAID58+djwYIFAIBly5bBxsYGY8eORUFBAYYOHYp169bB2tq6DlpNRGZP+R5gV1egwoQxs9KmjSIBzMwE8vIUS9cQEWlhFusAWiquI0RUj5WWAg4OgFwO9OgBnDpl6oiqNnEisH69YvvPP4EuXUwbD5GZ489vMxgDSERklq5eVSR/gPk+/lXiWoBEpCcmgEREmljC+D8l1YkgqnETEVWBCSARkSaWlACqxscEkIh0wASQiEgTJoBEVI8xASQi0kQ5Axgw/wTQ2xuwt1dsMwEkIh0wASQi0kQ1kfL3N1kYOrGyKo8xLQ3g4g5EVA0mgEREmigTwObNAUdH08aiC+VEkIICxZqARERaMAEkIqooP788iTL3x79KHAdIRHpgAkhEVNHly+XbTACJqB5iAkhEVJElzQBWYgJIRHpgAkhEVJElzQBWYgJIRHpgAkhEVJFqAqX6lg1zxgSQiPTABJCIqCJLfATs7g64uiq2mQASUTWYABIRVaRMoGxsgBYtTBuLPpTJano6UFpq2liIyKwZJQHMycnB9u3bce7cOWNUR0RkOkKUJ4AtWwLW1qaNRx/KBFAuB65eNW0sRGTWDEoAx44di08//RQAUFBQgF69emHs2LHo2rUrtmzZYtQAiYjq1L17QE6OYttSHv8qcRwgEenIoAQwNjYWAwYMAABs27YNQgjcv38f//3vf/HBBx8YNUAiojpliTOAlZgAEpGODEoAs7Oz4eHhAQDYvXs3Hn/8cTg5OSE8PBwXLlwwaoBERHXKEmcAK6nGywSQiLSwMeQkPz8/xMfHw8PDA7t378bGjRsBAPfu3YODg4NRAyQiqlOWOANYqR70AKanp+P27dt6ndOkSRO0bNmyliIiqp8MSgBnzZqFp59+Gs7OzmjZsiUGDRoEQPFouEuXLsaMj4iobllyAujvX75tgQlgeno6OnTshMKCfL3Oc3B0Qupf55gEEunBoARw6tSpeOihh3D16lUMHz4cVlaKJ8mtW7fmGEAismyWnAA6OgLe3kBmpkUmgLdv30ZhQT48R70GW08/nc4puXMVd375BLdv32YCSKQHgxJAAOjVqxe6du2KtLQ0tGnTBjY2NggPDzdmbEREdU+ZODk5AU2bmjYWQwQEKBLAzEygoECRFFoYW08/2Hu3NXUYRPWaQZNA8vPz8fzzz8PJyQlBQUFIT08HAMyYMQOLFy82aoBERHWmrAy4fFmx3bo1IJOZNByDqPZaKttCRFSBQQng3LlzkZycjAMHDqhN+hg2bBg2bdpktOCIiOpURgZQXKzYtrTHv0qqcasuaUNEpMKgR8Dbt2/Hpk2b8PDDD0Om8htyYGAgLl68aLTgiIjqlCWP/1OqBzOBiaj2GdQDeOvWLTRr1qzS/ry8PLWEkIjIojABJKIGwqAEsHfv3vj111+lPyuTvq+++gohISHGiYyIqK4xASSiBsKgR8CLFi3CiBEjkJKSgtLSUqxYsQJnz55FfHw8Dh48aOwYiattCA4AACAASURBVIjqRn1IAP38AGtrQC5nAkhEVTKoB7Bv3744cuQI8vPz0aZNG+zduxdeXl6Ij49HcHCwsWMkIqoblvweYCUbG0C5Hh4TQCKqgt49gKWlpfj+++8RFhaGb7/9tjZiIiIyDWXC1KQJ0LixaWOpiYAARVuys4F79wB3d1NHZHb4yjlq6PROAG1sbPDyyy/j3LlztREPEZFpFBUB168rti2190+p4jhAJoBq+Mo5IgPHAPbp0weJiYlo1aqVseMhIjKN9HRACMV2fUoAL10CevY0XSxmiK+cI6rBu4Bfe+01XLt2DcHBwWjUqJHa8a5duxolOCKiOlMfJoAocSawTvjKOWrIDEoAo6KiAChe/aYkk8kghIBMJoNcLjdOdEREdaW+JoB8HRwRaWBQApjG3yiJqL5RnQHcurXp4jAG9gASUTUMSgA59o+I6p361APo5QU4OgIFBUwAiUgjgxJAALh48SKWL1+Oc+fOQSaToVOnTpg5cybatGljzPiIiOqGMlGSycrX0bNUMhng7w+cO6d4BCyEYl89VibKcPTaURy8fBCnMk/hWs41ZBdmw87aDs0aNUMHzw7o17IfhgYMNXWoRGbBoARwz549ePTRR9G9e3f069cPQggcOXIEQUFB2LlzJ4YPH27sOImIapcyAWzRArCzM20sxhAQoEgACwuBzEygeXNTR1QryqxygUHAqH2jcLPwZpXl9lzcg/8e/y+sZdYY4DUACAAERN0FSmRmDEoA33rrLbz66qtYvHhxpf1vvvkmE0AisiwPHgB37ii2Lf3xr1LFcYD1LAEsQyFybH5Ctt9WoBU0Jn9Otk4olhejtKxU2icXchzIPABMBO4WroQnpsNetKvDyInMg0EJ4Llz5/Djjz9W2j958mQsX768xkEREdWp+jT+T6liAti3r+liMbJCqz9x23YZ5Fa3pH3WMmuEtQ3D6I6jEdIiBG092sLexh5CCGTlZSEpMwl7L+7FhjMbcCP3BgCgxOESMsVraCyPgFvJs7CCvamaRFTnDHoXcNOmTZGUlFRpf1JSEpo1a1bjoIiI6lR9mgGs5O9fvl1PJoIIlOGezf9w025eefInrIGjwM9DfsavT/2KF3q+gKBmQbC3USRzMpkMXs5eCGsbhk/CPkH6q+mI7hkNKN8CJyvDA5ufkWn/Gkpk10zTMCITMCgBfPHFF/HSSy9hyZIliIuLw6FDh7B48WL83//9H1566SVjx0hEVLsaQg+ghStDIW7bLUaO7SZAphi7Zy/vgibX3gR2A82ddHvEbWNlgzDfMOBzoPHdRyETikSxxOoybti/hkKr07XWBiJzYlAC+O677+K9997DypUrMXDgQPzrX//Cp59+igULFmDevHkGBfL5558jICAADg4OCA4ORlxcnNbyW7ZsQWBgIOzt7REYGIht27apHRdCYMGCBfDx8YGjoyMGDRqEs2fPqpU5f/48IiMj0aRJE7i4uKBfv37Yv3+/QfETkQVjAmjW5MjGTfu3kG99RLFDWMGtZBK8ij+ETamBT53KgEbZQ+BdtBS2ZYpZ30KWh5t27yLP6pCRIicyXwYlgDKZDK+++iquXbuG7OxsZGdn49q1a5g5cyZkBiw1sGnTJsyaNQvz5s1DYmIiBgwYgJEjRyI9PV1j+fj4eERFRWHChAlITk7GhAkTMHbsWBw7dkwq89FHH2Hp0qX49NNPceLECXh7e2P48OF48OCBVCY8PBylpaX4448/kJCQgO7du2PUqFHIzMzU/0shIstVHxNAd3fA1VWxbcEJoBw5uGn/Doqt/gYAyIQjmhXPh2vp45AZ9iNMjZ1oBe+iT+AgD1bskJXitt0S5Fr/XuO6icyZQf/3pKWl4cKFCwCAxo0bo3HjxgCACxcu4LIBrx1aunQpnn/+ebzwwgvo1KkTli9fDj8/P6xatUpj+eXLl2P48OGYO3cuOnbsiLlz52Lo0KHSBBQhBJYvX4558+ZhzJgx6Ny5M7799lvk5+fjhx9+AKB4Gfjff/+Nt956C127dkW7du2wePFi5OfnV+opJKJ6Tpkg2dvXr9myymT26lWgtFR7WTMkxwNk2b+LEivF34+18IB30X/gWBZs1OtYwRHNit9Fo9J/VrCQCdyxXYE8q8NGvQ6ROTEoAXzuuedw5MiRSvuPHTuG5557Tq+6iouLkZCQgNDQULX9oaGhGq8BKHoAK5YPCwuTyqelpSEzM1OtjL29PQYOHCiV8fT0RKdOnbB+/Xrk5eWhtLQUX375Jby8vBAcrPkfl6KiIuTk5Kh9iMjCCVGeALZqBVjVvFfJbCgTQLkcuGZZExwEinHLbiGKrS4CUCR/XkWLYCdq501UMtjAs2QGGpdG/LOjDLftPkaBVUKtXI/I1Az6ly4xMRH9+vWrtP/hhx/WODtYm9u3b0Mul8PLy0ttv5eXV5WPYjMzM7WWV/5XWxmZTIaYmBgkJiaicePGcHBwwLJly7B79264ublpvO6iRYvg6uoqffz8/PRqKxGZoawsID9fsV1fZgArWeg4QAGB27YrUGSdAgCwEm7wKvoQtsK3Vq8rgwzuJS+iUemwf3aU4pbdYhTLrtTqdYlMweAxgKpj6ZSys7Mhl8sNCqTi2EEhhNbxhLqU11ZGCIGpU6eiWbNmiIuLw/HjxxEZGYlRo0bhxv+zd+dxUZX9/8dfMwMMOyoIuIsrIu4rmlsmanbX7ZKWt2tlmS2SlrdklmlqLil3bt80K21R2/vd5YLelbnvuyYuKKaisu/LzJzfHyMDBCoMM8wAn+fjwcMzZ65zzgfK8c11znVdN28We83w8HDTM4/Jyclcu3atxN+fEMJOVcbn//JU0ACYXm0bGQ47AVApWnyzZ+GolM8v3CrUeOe+gqveOG+iosrkjtMcDOr0crm+EOXFrADYo0cP5s+fXyjs6fV65s+fz0MPPVSqc/n4+KDRaIr09t2+fbtID14ef3//+7b39/cHuG+bX3/9lZ9//pmNGzfSvXt32rdvz8qVK3FxcWHdunXFXler1eLp6VnoSwhRwUkAtC/NIK36VuO2osIn5w20SpNyLUGFBu+cKTgZjGvb69SxJPl+Zua/mELYJ7P+d164cCG//vorzZs3Z/z48YwfP57mzZvzxx9/sGjRolKdy8nJiQ4dOrB9+/ZC+7dv3063e8xcHxISUqR9ZGSkqX1AQAD+/v6F2uTk5LBz505Tm4y7t3zUf3veR61WYzAYSvU9CCEqMAmAduNmxk0YnP+6mm48roauNqlFjTM1c2agVoyPBOW4XIDeNilFCKswKwAGBQVx8uRJhg8fzu3bt0lNTWXMmDH8+eefBAcHl/p8U6ZM4eOPP+aTTz7h3LlzvPbaa8TExDBx4kQAxowZQ3h4uKn95MmTiYyMZMGCBfz5558sWLCAHTt2EBYWBhhv/YaFhTFv3jx++OEHTp8+zbhx43B1dWXkyJGAMURWr16dsWPHcuLECaKionjjjTeIjo5m0KBB5vxYhBAVUWUOgBVoNZAcfQ7Tj0wHF+NrV303PHWD73+QlTkovtTMCTeuNgLQAw7FHbJpTUJYillrAQPUrl2befPmWaSIESNGEB8fz+zZs7l58ybBwcFs3ryZBg2Mo71iYmIK9dR169aNjRs38tZbbzFz5kwaN27Mpk2b6NKli6nNtGnTyMzMZNKkSSQmJtKlSxciIyNNU9b4+PiwdetWZsyYwcMPP0xubi4tW7bkp59+ok2bNhb5voQQFUBlDoCuruDnB7du2X0AfPu3tzmddBoATa433rpXUVH6eWUtzdnQkmq6USQ5rgMVzDw2k6EPDcXH1cfWpQlRJipFUZTSHrR161bc3d1Nz/utWLGCNWvWEBQUxIoVK6hevbrFC7VHKSkpeHl5kZycLM8DClFRNWpkDEdeXpCYCGZMZm/XQkJg/37jdkYGuLjYtp5i7Lu2j4c+fQiDYgA9eMdOwb3GwyU6Njv2IrHrwjhy5Ajt27cv0TFHjx6lQ4cO+I+NQOv/4OcLFQzE8jo5LlEA/DPwn3w//HuzFj4Q9kH+/TbzFvAbb7xhmgPv1KlTTJkyhUcffZTLly8zZcoUixYohBBWo9NB3opDAQGVL/xB4V7Nq/Y3nUlGbgZjfxxrDH8Av4FjTn3bFvU3KtR43fkX3B0I/OOfP/L1ma9tW5QQZWT2SiBBQUGAcU3ef/zjH8ybN4+VK1eyZcsWixYohBBW89dfxkmSofLd/s1j5wNB3vzfm1xIMK4sFVwtGIqf/9/mNHov+CX/9ctbXuZO+h3bFSREGZkVAJ2cnEyjaHfs2GFacaNGjRqyOoYQouKozM//5bHjAHjw+kE+PPAhAM4Ozrzb7l2w50kYzkLfWn0BiMuI49Wtr9q4ICHMZ1YAfOihh5gyZQpz5szh4MGDplGzUVFR1K1b16IFCiGE1UgAtBm9Qc+Lv7yIgvEx9Dl95tDQvaFtiyqBacHTqOFSA4CNpzey+cJmG1ckhHnMCoDLly/HwcGBb7/9llWrVlGnjnF5ni1btjBgwACLFiiEEFZT1QLglSs2K+PvVh1exdGbRwFo7deasK5hNq6oZHycfYjoH2F6/cqWV8jSZdmwIiHMY9Y0MPXr1+fnn38usn/p0qVlLkgIIcrN5cv525VtHeA89eqBWg0Gg930AMamxTLj1xmm1ysfXYmD2uxZycrdqNajWHtsLTuv7uRy4mUW7F7AO73fsXVZQpSKLGwjhKi6CgaigpMmVyaOjpD3aI6dBMDXI18nJdv4vPgzbZ+he/3uNq6odFQqFSseXWEKrfN3z+dSwiUbVyVE6UgAFEJUXXmByN/fLufHs5i828AJCWDjgXq7ru7iy1NfAlDduTrvP/K+TesxV0vfloR1Md62ztZnMzVyqo0rEqJ0JAAKIaqmzEyIjTVuV9bn//LYyUAQg2IoFJTm951PTbeaNqunrN7u9Ta13GsB8NP5n/jj6h82rkiIkitxADx58iQGgz2PzxdCiFIoOCBCAmC5+PrM1xy6YVxLN9g3mOfaP2ezWizBQ+vBew+/Z3r9euTr+RNaC2HnSvzUbbt27bh58ya+vr40atSIQ4cO4e3tbc3ahBDCegoGoco6ACSPHQTAbF024f8LN71e1G8RGrXGYuc/d+6cVdo+yNg2Y4nYH8Gp26c4dOMQm05v4ulWT1vs/EJYS4kDYLVq1YiOjsbX15crV65Ib6AQomIrOAJYegCtbsWhFVxJugJAv0b96N+4v0XOq08zrt88atQoi5yvtDRqDYtDF9P/C+P3E/6/cAa3GIyzg7NN6hGipEocAIcOHUqvXr2oVasWKpWKjh07otEU/9vb5YIfrEIIYY+qwhyAeWwcABMyE5jzxxwAVKhY1G8RKgutu2zITgNFwfuxqTh61yvRMZmXD5O86wuLXB8gtHEooY1DibwUydXkqyw/uJzXu71usfMLYQ0lDoCrV69myJAhXLx4kVdffZUJEybg4eFhzdqEEMJ6qlIArFULtFrIzrZJAHx/9/skZSUBMKbNGNr4t7H4NRy966H1b1Kitrnx1yx+/UX9FrH90nYUFObumsuE9hPwcvay+HWEsJRSzbyZt8rHkSNHmDx5sgRAIUTFlReENJr8efIqK7UaGjSAqCjj960oYKEeuAeJTYtl+cHlAGg1Wub0mVMu1y1vrf1aM7btWD47/hlJWUks3b+UWb1n2bosIe7JrGlgPv30U1P4++uvv7h+/bpFixJCCKvLC4D164NDxVmFwmx5vZwZGXDnTrlddsHuBWTqMgF4seOL1PMq2W3aiujtnm+bJodeun8pCZkJNq5IiHszKwAaDAZmz56Nl5cXDRo0oH79+lSrVo05c+bI4BAhhP1LTITkZON2ZR8BnMcGzwFeT7nOqsOrAHBxcGH6Q9PL5bq2ElA9gGfaPgNASnYKH+z9wMYVCXFvZgXAGTNmsHz5ct5//32OHTvG0aNHmTdvHsuWLWPmzJmWrlEIISyrKo0AzlPw+yw4B6IVzd89n2x9NgAvd34ZP3e/crmuLc3oOQMnjRMA/znwH+6kl19vqxClYVYAXLduHR9//DEvvvgirVu3pk2bNkyaNIk1a9bw2WefWbhEIYSwsKo0ACRPOfcAxiTHsOboGgDcndyZ1n2a1a9pD+p71WdC+wkApOems3DPQhtXJETxzAqACQkJBAYGFtkfGBhIQoI88yCEsHMSAK1+ubl/zCVHnwPAq51fxcfVx+rXtBdv9ngTrUYLGOc/vJV2q8TH6nQ60tLSrFWaECZmBcA2bdqwfPnyIvuXL19OmzaWH94vhBAWVZVWAclTjgHwesp1Pj3+KQCeWk+mdpv6gCMql9oetXmx44sAZOoyidgfUaLjoqOj6dWrF+PHj0dRFGuWKETppoHJs3DhQgYNGsSOHTsICQlBpVKxd+9erl27xubNmy1doxBCWFbBZwCrSgCsUQM8PCA11eoB8IN9H5BryAXg5U4vU8OlhlWvZ49e7/Y6Kw+vJEefw4pDK/j3Q/+mmnO1e7b/8ssvmTRpEikpKQB89tlnjB8/vrzKFVWQWT2AvXr1IioqisGDB5OUlERCQgJDhgzh/Pnz9OjRw9I1CiGEZeUFQHd38KkityZVKmjY0Lh99Sro9Va5THxGPB8d+QgAZwdnJnedbJXr2Ls6nnUY22YsAKk5qaw8tLLYdsnJyYwaNYpRo0aZwl9AQECxj1kJYUlmT35Vu3Zt5s6da8lahBDC+vT6/FGwAQHlNiGyXQgIgFOnIDcXbtyAepafk+/DAx+SkZsBwHPtnsPXzdfi16gopnWfxtpjazEoBpbuX0pY1zBcHV1N7x86dIinnnqq0PKpY8aMYdmyZXh6etqiZFGFmNUDKIQQFdb168YABFXn9m8eKz8HmJqdyrKDywBwUDtU+fVwm9RowvCWwwGIy4hj7dG1gHEu3Q8++IBu3bqZwp+XlxcbNmxg3bp1Ev5EuZAAKISoWqri8395rBwAVx9ZTWJWIgD/avUvGlRrYPFrVDThD4WbthftXcT12Ov84x//4PXXX0en0wHQtWtXjh8/zlNPPWWrMkUVJAFQCFG1VMURwHmsGACzddl8sM+48oUKFf/u/m+Lnr+iau3XmseaPQbAtRPXCG4dXGiw5PTp0/njjz9omPd8phDlRAKgEKJqkR5AIwsHwHUn1nEz7SYAg1sMpkXNFhY9f0U2LWQa/Aash6Q7SQD4+vqybds25s+fj6Ojo20LFFWSWQHw4YcfJikpqcj+lJQUHn744TIXJYQQVlMVl4HLY6UAqDfoWbR3kel1wdueVd3169d5a+xbsBO4O7Vf666tOXHiBKGhoTatTVRtZgXA33//nZycnCL7s7Ky2LVrV5mLEkIIqykYAKvabbeC095YcD3gn6N+5mLCRQD6BvSlY+2OFjt3RbZjxw7atm3LH3/8YdyhAh4Gj+c88Pf3t2ltQpRqGpiTJ0+ats+ePUtsbKzptV6vZ+vWrdSpU8dy1QkhhKXlBcDatcHFxba12EJAAMTFwV9/QU4OODmV+ZRL9i8xbU8JmVLm81V0BoOBefPm8fbbb5tW9KhXrx6aYRqueF1hz197OHj9IJ3rdLZxpaIqK1UAbNu2LSqVCpVKVeytXhcXF5YtW2ax4oQQwqLS0uD2beN2VXv+L09AABw6BIoCMTHQpEmZTnfkxhH+uGrs4Qr0CWRAkwGWqLLCSkhIYPTo0YUGegwcOJDPP/+cH67+wIT/TgBg6f6lbBi6wVZlClG6W8DR0dFcunQJRVE4ePAg0dHRpq/r16+TkpLCM888Y61ahRCibAre9qzKATCPBZ4DXLp/qWk7rEsYalXVHVt4+PBh2rdvbwp/KpWKOXPm8PPPP+Pt7c2/Wv2Lmq41AfjmzDfEJMfYslxRxZWqB7BBA+OcTgaDwSrFCCGEVVXlEcB5CgTAqzt3Eu/tXaLDfHx8qF+/fqF911Ous+nMJgC8XbwZ3Wa05eqsQBRF4f/+7/8ICwszPR/v4+PDhg0beOSRR0ztXBxdmNRpEu/ufBe9omfZgWUsCl10r9MKYVVmLwUXFRXF77//zu3bt4sEwrfffrvMhQkhhMVV5RHAd912cyNvcbav5s7lzRIu6ens4sr5P88VCoErDq1AZzBOZjyx48RCy5xVFenp6bzwwgt8+eWXpn0hISF8/fXX1K1bt0j7SZ0m8f7u98nWZ7P66Gre7vU2HlqP8ixZCMDMALhmzRpefPFFfHx88Pf3R1VgLU2VSiUBUAhhn6QHkDvu7qYA2KJhO/x7jX3gMbnx14j/+QPi4uJMATA9J53/O/x/ADiqHXmp00vWKtlunT9/nqFDh3LmzBnTvsmTJ7Nw4UKc7jG4xtfNl9GtR/PxsY9JyU5h7bG1hHUNK6+ShTAxKwC+9957zJ07l3//W2Z6F0JUIBIAyfH3x4DxAfCG2elo/c0bBLL+xHrTsm9Pt3qaWh61LFdkBfDdd98xbtw40tLSAHB3d2ft2rUMHz78gceGdQ3j42MfA/CfA//h5c4v46A2+4acEGYx62ndxMREnnzySUvXIoQQ1pUXAJ2doYrOw6Y4OXH97nbd5FtmncOgGAoN/nit62sWqKxi0Ov1vPnmmwwbNswU/lq2bMmhQ4dKFP4AWvq2NI2WvpJ0hZ/+/Mlq9QpxL2YFwCeffJLIyEhL1yKEENajKPmjXgMCQF11R6vmjf31yUjGNSez1MdvvrCZCwkXAOjTsA9t/dtasDr7FhYWxvz5802vR44cyYEDBwgMDCzVeQqG5mUHZfo0Uf7M6nNu0qQJM2fOZP/+/bRq1arIOoavvvqqRYoTQgiLiY2FrCzjdhUdAJInGuh5d7tu8i2iajYs1fER+yNM21Vl4ufcROM6x3v37gVAo9GwePFiJk+eXOg5+JJ6pNEjNPduzvn48+y8upNTt07Ryq+VRWsW4n7MCoCrV6/G3d2dnTt3snPnzkLvqVQqCYBCCPsjz/+ZXCmwXa+UAfB83Hn+F/0/AJrUaMKjTR+1aG32KP3P3SRszr/l7ePjw9dff02fPn3MPqdapeblzi/zypZXAFh+cDkf/eOjMtcqREmZFQCjLbiIuBBClAsJgCYFP8HrJt8u1bF5I38BJnaYWKknflYMepJ2fU7K/m9N+5o3b862bdtM8+KWxZg2Ywj/XzhpOWl8ceoL3n/kfaq7VC/zeYUoiTL9zc3JyeH8+fPodLoyF7Jy5UoCAgJwdnamQ4cO7Nq1677tv/vuO4KCgtBqtQQFBfHDDz8Uel9RFGbNmkXt2rVxcXGhd+/ehYbq5/nll1/o0qULLi4u+Pj4MGTIkDJ/L0IIOyQB0KTAT4L6SbH3bPd3mbpMPjvxGQBajZZxbcdZtC57os9M5fa37xYKfwCffPKJRcIfgKfWk3FtxgGQkZvBJ8c+sch5hSgJswJgRkYGzz77LK6urrRs2ZKYGONyNq+++irvv/9+qc+3adMmwsLCmDFjBseOHaNHjx4MHDjQdN6/27dvHyNGjGD06NGcOHGC0aNHM3z4cA4cOGBqs3DhQpYsWcLy5cs5dOgQ/v7+9OvXj9TUVFOb7777jtGjRzN+/HhOnDjBnj17GDlyZKnrF0JUABIATQoHwJslPm77je0kZSUB8FTwU3i7lmwVkYom584VYte/Rlb0UeMOlRqPToMBcHZ2tui1Xu78sml7xaEV6A16i55fiHsxKwCGh4dz4sQJfv/990J/GR555BE2bdpU6vMtWbKEZ599lueee44WLVoQERFBvXr1WLVqVbHtIyIi6NevH+Hh4QQGBhIeHk7fvn2JiDA+mKwoChEREcyYMYMhQ4YQHBzMunXryMjI4KuvvgJAp9MxefJkFi1axMSJE2nWrBnNmzdn2LBhZvxEhBB2r+CjK1V8EMgNIFNjHLzXILHkPYDfXs3vDXux44uWLssupP+5m9jPp6K72zOqdvHE76n3cAvqZZXrNfdpTmjjUACik6LZcnGLVa4jxN+Z9Qzgjz/+yKZNm+jatWuh0U9BQUFcunSpVOfKycnhyJEjTJ8+vdD+0NBQ02irv9u3bx+vvVZ43qn+/fubAmB0dDSxsbGEhoaa3tdqtfTq1Yu9e/fywgsvcPToUa5fv45araZdu3bExsbStm1bFi9eTMuWLYu9bnZ2NtnZ2abXKSkppfpehRA2lNcDWLMmuLvbthYbU4AYd2+aJ8fSIOkmKsWA8qBn+WrBmSTjYzTt/NvRuU7nB14nJiaGuLi4Etd17ty5Ere1NEUxkLx7A8l7N5j2Ofk1puaQGTh4+pIde9Fq136508tEXjJOrbbs4DIea/aY1a4lRB6zAuCdO3fw9fUtsj89Pb3Uw+Hj4uLQ6/X4+fkV2u/n50dsbPG/mcbGxt63fd6fxbW5evUqAJfv/mMwa9YslixZQsOGDfnggw/o1asXUVFR1KhRo8h158+fz7vvvluq708IYQeysuD63emPq/jt3zxXPIwBUKvPxS81gVhPn/sf0Cl/88WOLz7wsz4mJobmgS3IysywQLXWZcjNIv6XpWSc32Pa59ayDzX6v4zaUWv16z/a9FECqgUQnRRN5KVIzsedp7lPc6tfV1RtZgXATp068csvv/DKK8bh63kfBGvWrCEkJMSsQv7+YaIoyn0/YErS/n5tDAYDADNmzGDo0KEAfPrpp9StW5dvvvmGF154ocg1w8PDmTIlf86rlJQU6tWr96BvTQhha1eu5G9LAATgqkd+4GuQdPO+AdCgzoC7U9R5aj0Z2erBz0rHxcWRlZmB92NTcfQu2edk5uXDJO/6okRtLUWXGsed798jx9TDp6J6n/F4dBps1vx+5tCoNbzU6SVe3/46YJwSZtmjMjm0sC6zAuD8+fMZMGAAZ8+eRafT8Z///IczZ86wb9++IvMCPoiPjw8ajaZIb9/t27eL9ODl8ff3v297/7tLPMXGxlKrVq1i2+TtDwoK9x8YvgAAIABJREFUMr2v1Wpp1KjRPQefaLVatFrr/zYohLAwGQBSRHSBwNcg8SYH6t97EuJM90Nwd77/Ma3H4ObkVuLrOHrXK/F6w7nx10p8XkvIjYvhznez0aclAKBycsHn8Wm4Nu70gCMt75l2zzDzt5mmkdZz+87FU+tZ7nWIqsOsQSDdunVjz549ZGRk0LhxYyIjI/Hz82Pfvn106NChVOdycnKiQ4cObN++vdD+7du3061bt2KPCQkJKdI+MjLS1D4gIAB/f/9CbXJycti5c6epTYcOHdBqtZw/f97UJjc3lytXrlhsiL8Qwk5IACzi7z2A96KgkOGRf2t0YseJVq2rPMVvXWYKfxovP/xHLbJJ+AOo7lKdUa1HAZCWk8b6E+ttUoeoOszqAQRo1aoV69ats0gRU6ZMYfTo0XTs2JGQkBBWr15NTEwMEycaP2jGjBlDnTp1TOsvTp48mZ49e7JgwQKeeOIJfvrpJ3bs2MHu3bsB463fsLAw5s2bR9OmTWnatCnz5s3D1dXVNM2Lp6cnEydO5J133qFevXo0aNCARYsWAca1joUQlUjBAFjFRwDnuVIoAN57JHC2+hR6J+Nk0e1rtKelb/GD5CoKRVHIiLo7wFCfC4C2bhA1B89A4+plw8rglc6v8HPUz0zsOJGhLYbatBZR+ZkVADdv3oxGo6F///6F9m/btg2DwcDAgQNLdb4RI0YQHx/P7NmzuXnzJsHBwWzevNnUExcTE4O6wMLt3bp1Y+PGjbz11lvMnDmTxo0bs2nTJrp06WJqM23aNDIzM5k0aRKJiYl06dKFyMhIPDw8TG0WLVqEg4MDo0ePJjMzky5duvDrr79SvbrMxC5EpXKxwAjOJiW7HVnZ/eVWHZ1KjYNioEHijXu2S9VsNm0Pa1ixp8ky5GYTvzmCzKh9pn1uwY/g3f8lVA6O9zmyfLTya0XMazE4qM3umxGixMz6v2z69OnFTvisKArTp08vdQAEmDRpEpMmTSr2vd9//73IvmHDht13zj6VSsWsWbOYNWvWPds4OjqyePFiFi9eXNpyhRAVSV4A1GqhTh3b1mIncjUO3PCsSf3kW8YeQEWBvw160JFAhuZuWEqDh2s9bINKLUOXGs+dH94j5+YF0z6PDo9Tve+EchvsURIS/kR5MesZwAsXLhQaPJEnMDCQixetN1eSEEKUml4PefOTNm4M6sq7dm1pXa1mHAznmZ1OtazUIu+nOUSC6u7KFMfAUW37XjJzZMdeJHb9lPzwd3cSbLfgh+0q/AlRnsz6JPTy8jLNo1fQxYsXcXMr+egwIYSwuuvXISfHuC23fwu5Wj1/loSGiYUHgijoSdNsu/tCBYfLszLLSf9zN7e+/Df6tHgANJ418er+tI2rEsL2zAqAjz/+OGFhYYVW/bh48SJTp07l8ccft1hxQghRZgXvSjRtars67FBeDyAUXRM4U30YvfoOANrMIEgu19LKTFEUkvZuJO6n91F0xhWctHVaUGvMUhw8a9q4OiFsz6wAuGjRItzc3AgMDCQgIICAgABatGiBt7e3PE8nhLAvMgDknmKq+5u2G/ytBzDVIX/wh2tK93KryRIMudnE/XdxoUml3Vr2we+puWjcqtmwMiHsh1lPm3p5ebF37162b9/OiRMncHFxoXXr1vTs2dPS9QkhRNlIALyngj2ADQv0AOaqYslSHwVAY/DDKTOw3Gszly4twbiyx82ou3tUVOs1Fs8uQ+V5PyEKKHUAzM3NJTQ0lI8++ojQ0FBCQ0OtUZcQQljGhfxRnxIACyt0Czgxfy7ANM0WUCkAeOgHoDLvZlG5y7l1idvfzUGfGgeAytEZn8em4trMvCVKhajMSh0AHR0dOX36tPwmJYSoGPJ6AB0dQdbuLiTTyZnbbtXxTU80rQaikEOaw91VlBQH3HX90BFnwypLJuvqSZL3fImSa3zeT+NRE9+hM3Hyk5VfhCiOWb/WjRkzhrVr11q6FiGEsCyDIX8KmEaNQKOxbT12KK8X0Dc9EdecTNI1ezCoUgBw1XdHg30/M6coxp7KpN8/MYU/p1rNqTVmiYQ/Ie7DrGcAc3Jy+Pjjj9m+fTsdO3YsMvXLkiVLLFKcEEKUyc2bkJlp3Jbbv8W6Wr0Wna6fBaB+UiyX3fMHf3joSz+pf3lSdDmkHd9SaJ9bUG+8B76KysHJRlUJUTGYFQBPnz5N+/btAYiKiir0ntwaFkLYDRkA8kBXq+WPBHbPOkG25hwAjob6aA32u+6vPj2R29+/R86N86Z91XqOwbPrk/LvkBAlYFYA/O233yxdhxBCWJ4EwAcqOBn0Nfddpm0P3aOosM8glXP7Mre/nYM+9Y5pX7Xe4/HqMtSGVQlRsZRpaNfFixfZtm0bmXdvseQ9iyGEEHZBAuAD5T0DmOoEJ32NPy+V4oyb3j7X/c24sJ/YL6aZwp/a2R0A5wZtbFmWEBWOWQEwPj6evn370qxZMx599FFu3jSOHnvuueeYOnWqRQsUQgizySogD5TXA/hFa8hyNK7766bvjRpXW5ZVhKIoJO//ljvfz0XJzQLAqVYzvB76l40rE6JiMisAvvbaazg6OhITE4Ora/6HxIgRI9i6davFihNCiDLJC4AODtCggW1rsVNJzh4kaV1Z1Sl/n4fuUdsVVAxFl0v85qUk7fwMMN5pcm3RE7+n55t6AIUQpWPWM4CRkZFs27aNunXrFtrftGlTrl69apHChBCiTBQlPwA2bGgMgaIolYofW3hzyi8DAGddc5wU+5k+RZ+exJ0f5pF9d6QygNdD/8Kr21My2EOIMjDrEzE9Pb1Qz1+euLg4tFptmYsSQogyu3UL0tKM2/L833192i7HtN04qQtpbvdpXI5y7lzh9rez0afcBkDloMV70Gu4BT5k48qEqPjMugXcs2dP1q9fb3qtUqkwGAwsWrSIPn36WKw4IYQwmwwAKRE9yeypZxxQ4Z0BXa77P+CI8pFx8QCxX7xhCn8a9xr4/WuBhD8hLMSsHsBFixbRu3dvDh8+TE5ODtOmTePMmTMkJCSwZ88eS9cohBClJwGwRNIctqNXGwAYfwx8cm7xPxvXlH76V1KP/Je85/2c/JtSc8hbOHh427YwISoRs3oAg4KCOHHiBJ06daJfv36kp6czZMgQjh07RuPGjS1doxBClJ4EwAdSMJCmyV9J44UjEJBw3Xb16HUApB75f5gGewT2wG/kfAl/QliY2U9F16pVi9mzZ1uyFiGEsBwJgA+UpT6KTn0LgNCL0CQBklxsEwD1GcmkHPi20D6v7iPx6v60DPYQwgpK1QOYkZHBSy+9RJ06dfD19WXkyJHExcVZqzYhhDDfhQvGP9Vq4yhgUUSqQ/66v0+dNk6nYosewJw7V7i5fgq6vGtrHPF5fBrVHhop4U8IKylVAHznnXf47LPPGDRoEE899RTbt2/nxRdftFZtQghhHkWBvHXKGzYEmZ2gCL0mkUz1YQA0Bh8aJwYA4J2ZgldmarnVkXHpkHGwR/It0z7vAa/g1qJnudUgRFVUqlvA33//PWvXruWpp54CYNSoUXTv3h29Xo9Go7FKgUIIUWo3b+ZPAdO8uW1rsVMZnntBZRz84a7vz7Xq8XD1FAABiTc47mLdn5uiKKQe+pHE3z4h73k/jZcv+uTbOPrUt+q1hRCl7AG8du0aPXr0ML3u3LkzDg4O3Lhxw+KFCSGE2c6fz9+WAFiUBjLd9xu3FQ3uulAuV69jetvat4EVfS4JW5eR+NtaTIM9mnfHK+Qpq15XCJGvVD2Aer0eJyenwidwcECn01m0KCGEKJOCAbBZM9vVYa8CweBgvM3rauiKA95E1yifAKjPSObOj/PJvnbatM+r21N4PTSS9LM7rXZdIURhpQqAiqIwbty4Qqt9ZGVlMXHiRNzc8qeO//777y1XoRBClJb0AN5fx/xN97vr/hYMgI0SrXNXJycuhjvfzUaXFGvcoXHE59HJuAX1tsr1hBD3VqoAOHbs2CL7Ro0aZbFihBDCIvIGgEClD4AxMTGlmo1hx4kdYBzvgYOhLs6G1gBc8/JDp1LjoBis0gOYefkId35agJJjXHNY7VYN3yEz0dauWP99zp07V6r2Pj4+1K8vzzQK+1OqAPjpp59aqw4hhLCcvB5ANzeoXdu2tVhRTEwMzQNbkJWZUfKDBgBdjZseuoGoME6zotM4EFPNn0aJNwhIvG4cSW2BKVgURSHl8E8k/roWFOOgE0ffRvgOnYmDZ80yn7+86NMSQaUqdaeHs4sr5/88JyFQ2B2zJ4IWQgi7lJ0N0dHG7WbNLBJi7FVcXBxZmRl4PzYVR+96D2xvUGVzu+5bQC4YHHHT9y30fnSNOjRKvIFrbjZ+afHc8vApc40p+zaReWG/6bVLsxB8Bk1F7eRc5nOXJ0N2GihKiX/WALnx14j/+QPi4uIkAAq7IwFQCFG5XLoEBmNPU2W//ZvH0bseWv8Hr3aSqtkGDrkAuKS3Q+PgXuj96Or5vaWNEm6UKQAmJiYCFAp/niHDqdZjFCqVWauQ2oWS/qyFsHcV92+hEEIUpwo9/1caCgppBVb+cE15qEibQiOBE81/DvD06dOFnxnXOOLzj9ep3nNMhQ5/QlQm8jdRCFG5yAjgYuWooshRXzK+uAGOOUVvSV62wFQwP//8MyEhIVy/bjxe7eKJ/8j3ZaSvEHZGAqAQonKRAFisVIct+S8OFd8mugyTQSuKwuLFi3n88cdJy1uFBfAeNKXCjfQVoiqQACiEqFwKBsCmTW1Xhx3Rk0KG5g8AVDotnC6+3S2PGmQ4Gud5bVSKAJidnc348eN54403UBTjyh59+xoHmGjcqpWhciGEtUgAFEJULnnPANauDR4etq3FTqQ57EBR5QCgjQuG3OLbKSo1V+4OBKmfFIuD/sGrPN2+fZuHH36YdevWmfa98847vP/++2UvXAhhNRIAhRCVR0IC5E2KLLd/AVAwkKbJH/zhfKfNfdtfrlEXAAfFQIOkm/dte/LkSTp16sTevXsBcHFxYdOmTcyaNQu1Wv55EcKeyd9QIUTlIWsAF5GlPopObVx6zVnfDk129fu2v1Bgjrsmcdfu2e7HH3+kW7duxMTEAFC7dm127drF8OHDLVC1EMLaJAAKISqPs2fzt1u0sF0ddiTV4RfTtodu0APbXywYAOOLBkBFUZg/fz6DBw8mPT0dgE6dOnHo0CE6dOhggYqFEOVBAqAQovIouE5rUJDt6rATuapYMtWHAdAYauJi6PTAYy765AfApvExhd7Lzs5m9OjRvPnmm6Z9Tz/9NDt37qR2JV5yT4jKSFYCEUJUHgV7ACUAkqbZCirjqFwP/UBUaB54THT1OuhVajSKgSbxfxV67/nnn+f06fwhxO+99x5vvvkmqkq83J4QlZUEQCFE5ZEXAD09jaOAqzCFHNIcIu++cMBd169Ex+U4OHK1mj+NEm/QOP4vVIqBnDtXAEzhz9XVlS+++ILBgwdbo3QhRDmQW8BCiMohLQ2uXjVuBwVBFe+VStfsxqBKAcBV3x0N9x/8UdBFH+MqIS66bBwOfE/C1mWm9+rXr8+ePXsk/AlRwdlNAFy5ciUBAQE4OzvToUMHdu3add/23333HUFBQWi1WoKCgvjhhx8Kva8oCrNmzaJ27dq4uLjQu3dvzpw5U+y5srOzadu2LSqViuPHj1vsexJClKOCI4BlAEjhwR/6Bw/+KOiid11ygVeAizs/A4MegPbt23P48GHatm1ruUKFEDZhFwFw06ZNhIWFMWPGDI4dO0aPHj0YOHCgaXqBv9u3bx8jRoxg9OjRnDhxgtGjRzN8+HAOHDhgarNw4UKWLFnC8uXLOXToEP7+/vTr14/U1NQi55s2bZo8wCxERSfP/5lkqy6SozYGYkdDQ7SG0gXikx41CQWW/23/ypUrqVmzpmWKFELYlF0EwCVLlvDss8/y3HPP0aJFCyIiIqhXrx6rVq0qtn1ERAT9+vUjPDycwMBAwsPD6du3LxEREYCx9y8iIoIZM2YwZMgQgoODWbduHRkZGXz11VeFzrVlyxYiIyNZvHix1b9PIYQVSQA0SXPIn/jZQzcIFSW/HZ5z6zLr927k97uvNSoVnt2eAsDR0dGCVQohbMnmATAnJ4cjR44QGhpaaH9oaKhpdvm/27dvX5H2/fv3N7WPjo4mNja2UButVkuvXr0KnfPWrVtMmDCBzz//HFdX1wfWmp2dTUpKSqEvIYSdkAAIgJ400jU7AVAprrjpe5f42PRzfxD7xRtkpicC4A+s8a6Pa9OuVqhUCGFLNg+AcXFx6PV6/Pz8Cu338/MjNja22GNiY2Pv2z7vz/u1URSFcePGMXHiRDp27FiiWufPn4+Xl5fpq169eg8+SAhRPvICoKsr1K9v21psKF2zA0WVDYC7vi9qXB54jGLQk7hzHXH/byGKznhsW7UDh4GhqXGgKNYsWQhhAzYPgHn+Po+Uoij3nVuqJO3v12bZsmWkpKQQHh5e4hrDw8NJTk42fV27du9lkoQQ5SgrCy5fNm4HBkIVXYdWwUBqodu/jz7wGENWGne+m0PK/m9M+9yC+7KgXjB1AM/sdHwz5W6HEJWNzT8lfXx80Gg0RXr7bt++XaQHL4+/v/992/v7+wPct82vv/7K/v370Wq1ODg40KRJEwA6duzI2LFji72uVqvF09Oz0JcQwg5ERYHBYNyuwrd/s9RH0KlvAOCsb42jcv+7FLnx17j5+VQyLxtXC0GlpnrfCXg/GkZMzQamds2Sb1mtZiGEbdg8ADo5OdGhQwe2b99eaP/27dvp1q1bsceEhIQUaR8ZGWlqHxAQgL+/f6E2OTk57Ny509Tmww8/5MSJExw/fpzjx4+zebPxt+ZNmzYxd+5ci31/QohyIM//AZDi8F/Ttofu8fu2zYo5xc31U9ElXAdA7eyB74g5eHZ8ApVKVWhN4KZJEgCFqGzsYiWQKVOmMHr0aDp27EhISAirV68mJiaGiRMnAjBmzBjq1KnD/PnzAZg8eTI9e/ZkwYIFPPHEE/z000/s2LGD3bt3A8Zbv2FhYcybN4+mTZvStGlT5s2bh6urKyNHjgSMk5kW5O7uDkDjxo2pW7dueX3rQghLkABIruoaWZqjADgY/O657q+iGHtKk35ba9rnWLMhNYe8hWM1f9O+Cz75n5GBSTetUbIQwobsIgCOGDGC+Ph4Zs+ezc2bNwkODmbz5s00aGC8BRETE4O6wDM93bp1Y+PGjbz11lvMnDmTxo0bs2nTJrp06WJqM23aNDIzM5k0aRKJiYl06dKFyMhIPDw8yv37E0JYWcEAWEUngS7U+6d/rNh1f/UZyaQe+L7QPtfAHngPnIzaybnQ/qgCt4ADEyUAClHZ2EUABJg0aRKTJk0q9r3ff/+9yL5hw4YxbNiwe55PpVIxa9YsZs2aVaLrN2zYEEVGuglRMZ06ZfzT2RkaN7ZtLTZgUGeQrvkfACrFudh1f7NjL3Lnh3noU24bd6jUVO/zDB53b/n+XarWjb88a1I35Y70AApRCdlNABRCCLNkZMCFC8btli1BU7Tnq7LLdD/wt6lf3Au9n3YykvjIVaDPNe2rEToJj7YD7nve8zUbUjflDp65WVTdiXWEqJxsPghECCHK5OzZ/HnqWrWybS22oIIMz/y10z10/zBtK7pc4rcuJ37Lh6bw51CtFgBO/k0eeOrzBW4DV8GfrBCVmgRAIUTFdvJk/nbr1rarw1aag94xAQBnfQccFeMgNl3KHWK/mkbaia2mph7tB+HZbUSJT/1nzYambQmAQlQucgtYCFGx5T3/B1WzB7DAKm2ed6d+ybx6grifFmC4O4GzysGJGv1fwj24L2lnfivxqc9LABSi0pIAKISo2AoGwCrWA3g++Tw0NG47GOqi1bch+cC3JP2xHu5O9+Lg5UfNwTNw8mtU6vNfrlGHXLUGR4OeVkDuA48QQlQUEgCFEBVb3i1gX1/jVxWyKXqTads9pR9x/51L5qVDpn3OjTrg89jraFzMm/4qV+PIpRp1CYy7SiBwOlcioBCVhTwDKISouG7dgjt3jNtVrPcvNi2Wzdfvrvsb40TqJ/8tEP5UeHV7Gt9h75gd/vLk3QZ2BLRXrpTpXEII+yE9gEKIiqvgAJAq9vzfsgPLyNXnwgEgMhe9IQ4AtYsnPo9NxaVRB4tc53zNBnDOuO1y8aJFzimEsD0JgEKIiquKDgBJy0lj5e6V8DV3w5lxGhxtnSB8Hp+Gg6ePxa5VcCSwBEAhKg8JgEKIiquKTgHz7oZ3SfpPEiTm7/PsMpRqPUaj0lj2Y/28BEAhKiUJgEKIiiuvB1CthqAg29ZSDhRFYeWqlSx+dTHo7+5UQbU+z+HV6Z9WueZ1z5qkODrjmZuFswRAISoNGQQihKiYdDo4c8a43bQpuLjYth4rS01N5V//+hcvv/SyKfy5N3AHBZzrBVvvwioV56rXBkAbGwsJCda7lhCi3EgAFEJUTFFRkG1c/7ayP/935MgROnTowIYNG/J3doH5K+aXy/VP1aiT/+LYsXK5phDCuiQACiEqpqNH87fbt7ddHVZkMBhYsmQJISEhXLhwwbhTCwyHduPbEeIfUi51nPKum/9CAqAQlYI8AyiEqJgqeQC8ffs248aNY8uWLaZ9Xo28SH4sGWrAG93eQJWrKpdaTtcoEAAL/tyFEBWW9AAKISqmgkGkXTvb1WEFO3bsoE2bNoXC3/iXxpM80hj+6nvVZ1jQsHKrJ6qaP9l5LyQAClEpSAAUQlQ8BkP+rci6dSvNEnC5ubmEh4cTGhpKbGwsAL6+vmzduhWln2K6Z/Na19dw1DiWW106tQbThDtRUZCWVm7XFkJYhwRAIUTFc/kypKQYtyvJ7d/o6Gh69OjB+++/j6IYJ3bu168fJ06cILBLIF+c/AKAas7VeLbds+Ven+nJP0WBEyfK/fpCCMuSACiEqHiOHMnf7mCZJc9sacOGDbRt25YDBw4A4ODgwMKFC9m6dSv+/v4s2rsInUEHwCudX8FDW7b1fc1R6Mav3AYWosKTQSBCiIqnkgwASUxMZNKkSWzcuNG0r1GjRmzYsIHOnTsDcDP1Jh8f/RgAN0c3JneZbJNaC439lZHAQlR40gMohKh4KkEA/PXXX2ndunWh8Ddy5EiOHTtmCn8AS/YtIVtvHILxYscX8Xb1LvdaAU4CikZjfCE9gEJUeBIAhRAVi6LkBxA/P6hVy7b1lFJWVhZTp06lb9++/PXXXwB4eXnx1Vdf8eWXX+Lp6WlqG58Rz6rDqwDQarRMCZlik5oBsoCshg2NL86cyZ+EWwhRIUkAFEJULDEx+cuRtW8PqvKZC88STp48SefOnVmyZIlpX58+fTh16hRPP/10kfYfHviQ9Nx0AJ5t9yy1PGwbdjMCA40bOh2cPm3TWoQQZSMBUAhRsVTA278Gg4HFixfTqVMnTp06BYCTkxMffPABO3bsoF69ekWOSclO4cODHwLgoHZgWvdp5VpzcUwBEODwYdsVIoQoMwmAQoiKpeAI4AoQAGNiYnjkkUd44403yMnJAaBVq1YcPnyYKVOmoFYX/zG88tBKkrKSABjdejQNqjUot5rvJSMoKP/FwYO2K0QIUWYSAIUQFcvdqVIA6NjRdnU8gKIorF69muDgYH777TcAVCoVr7/+OocOHaJVq1b3PDY1O5XFexcDoFapmf7Q9HKp+UEyAgPB4e7kEfv327YYIUSZyDQwQoiKQ6/PD4C1akExt07tQUxMDM899xzbt2837atbty7r16+nT58+Dzz+wwMfEp8ZD8BTwU/RzLuZ1WotDcXZGVq3Nt6GP3fOOBl3gUErQoiKQ3oAhRAVx59/QmqqcbtrV7sbAKIoCh9//DHBwcGFwt8zzzzD6dOnSxT+krKSWLwvv/fvnV7vWK1es3TpYvxTUeDQIdvWIoQwmwRAIUTFUfD2b9eutqujGNeuXWPgwIFMmDCB1LshtU6dOmzevJm1a9fi5eVVovMs3bfU9OzfmDZj7Kb3zyQvAELh/x5CiApFAqAQouIo+NyZnQRARVFYu3YtwcHBbNu2zbR//PjxnD59moEDB5b4XPEZ8SzdvxQwjvx9u+fbFq+3zCQAClEpyDOAQoiKIy8AajR2sQZwTEwMEydOZMuWLaZ9tWvXZs2aNQQHB3P58uVSnW/tlbWk5hh7D59p+wwB1QMsWm9ZnTt3Dpo3p7W7Ow5paeTu2cOpI0eKvRV/7tw5G1QohCgpCYBCiIohNTV/8uHWrcHNzWal6PV6Vq5cyZtvvklaWppp/9ixY1m6dCmpqak0D2xBVmZGyU/qBkwGnMBJ48RbPd+yeN3m0qclgkrFqFGjANgGhAKO8fEM7tiRGJtWJ4QwhwRAIUTFcPiwceAB2PT276lTp5gwYQIHCtz+rFWrFmvWrGHQoEEAREdHk5WZgfdjU3H0LtlI5SSXz8lyMs5x+Hz756nnZT8jnA3ZaaAopu/n9LHNhJ6MBGBAzzH8v4Ci8zFmXj5M8q4vyrtUIUQJSQAUQlQMNn7+Lysri/fee48FCxag0+lM+59//nkWLFhAtWrVihzj6F0PrX+TB55bp7pNltMJALRqLW/2eNNyhVtQ3vdzqlkI3A2AnTKS2VbM95gbf628yxNClIIEQCFExWDDALhz506ef/55oqKiTPuaNWvGmjVr6NmzZ5nPn+TwOaiNoXJEwAibr/n7ICdq5Y9MbnfjvA0rEUKYS0YBCyHsn6LkB8Dq1aFp03K5bGJiIs8//zy9e/c2hT8HBwfeeustTpw4YZHwl6O6RLrmd+OLDHim6TNlPqe1Jbh6cbl6bQBaxV5Em5tt44qEEKUlAVAIYf8uXIDbt43b5TABtKIorFu3jubNm7NmzRrT/q5du3Ls2DHmzJmDs7Nz2a+DQqLjp6C6+2zjH+Dh6FHm85aHg/WCAXAy6Gh3U3oBhahoJAC3GxQ2AAAgAElEQVQKIezfzp352716WfVSp0+fplevXowbN447d+4A4O7uzrJly9i9ezfBwcEWu1aW+ihZmuMAaHK9oQItrHGwXkvTdudrZ2xYiRDCHBIAhRD2rxwCYGpqKq+//jpt27Zl165dpv1Dhw7l7NmzvPzyy2g0GotdT0Fv7P27yz1xEOgtdnqry+sBBOh87bQNKxFCmEMCoBDCvilKfgB0dbX4BNCKovDNN98QGBjIBx98gF5vTGFNmjRh69atfPvtt9SrZ/kpWdI1v5GrvgKAk6EpzultLX4Na/rL05frHjUB6HD9Txz1uTauSAhRGhIAhRD2LToa/vrLuN2tGzg6WuzU58+fp3///gwfPpwbN24A4OzszOzZszl16hT9+/e32LUKMpBJkuPnptfVc8ejqmgfxyqV6Tawiy6b4NhLNi5ICFEaMg2MEMK+Wfj2b0xMDNHR0axZs4aNGzeaevwAunfvzrRp06hbty5nz5417ffx8aF+/fplvnaeZIdv0KviAXDRd8LZ0JpsLgKlW0LN1sutHawXzOCzvwPQ+a/THKsTaNN6hBAlZzcBcOXKlSxatIibN2/SsmVLIiIi6NGjxz3bf/fdd8ycOZNLly7RuHFj5s6dy+DBg03vK4rCu+++y+rVq0lMTKRLly6sWLGCli2Nv7FeuXKFOXPm8OuvvxIbG0vt2rUZNWoUM2bMwMnJyerfrxCihCwYAKOjo2nWvDm63OJvV+7Zs4cnnniiyH5nF1fO/3nOIiEwV3WTFIfvjS8UB6rnPgcUXW6tIvj7QJCPugyzYTVCiNKwiwC4adMmwsLCWLlyJd27d+ejjz5i4MCBnD17ttgP3H379jFixAjmzJnD4MGD+eGHHxg+fDi7d++mS5cuACxcuJAlS5bw2Wef0axZM9577z369evH+fPn8fDw4M8//8RgMPDRRx/RpEkTTp8+zYQJE0hPT2fx4sXl/SMQQtxLXgB0dobOnc0+zW+//cYLL7xQOPxpHHELfhi3lg+jdtQWe1xu/DXif/6AuLg4iwTARMePQWWc9NlT908clTpA0eXWSsLWy61dqlGXOFcvfDKS6fTXWdQGPQa15QbKCCGsxy4C4JIlS3j22Wd57jnjb8IRERFs27aNVatWMX/+/CLtIyIi6NevH+Hh4QCEh4ezc+dOIiIi2LBhA4qiEBERwYwZMxgyZAgA69atw8/Pj6+++ooXXniBAQMGMGDAANM5GzVqxPnz51m1apUEQCHsRUwMXLli3O7aFbTFh7T7uXz5Mm+88Qbff/99of2uLXpSvfc4HDx9LVBoyWSqj5CpMa4hrFFq4KUbUaRNSZePAztYbk2l4mDdljwatRfP7HRa3LnCGb/Gtq1JCFEiNn/qOCcnhyNHjhAaGlpof2hoKHv37i32mH379hVp379/f1P76OhoYmNjC7XRarX06tXrnucESE5OpkaNGvd8Pzs7m5SUlEJfQggrKsPt36SkJKZPn06LFi2KhL8aA16l5uPTyjX8KeSS4Jg/qXS13PGocSm361vL/vqtTNvdrxy3YSVCiNKweQCMi4tDr9fj5+dXaL+fnx+xsbHFHhMbG3vf9nl/luacly5dYtmyZUycOPGetc6fPx8vLy/TlzWmhhBCFBAZmb/du3eJDsnOzmbJkiU0btyYBQsWkJOTAxj//r/99tsAOPk1snSlD5Ts8B06tXE0s1bfAjd973KvwRp2BbQ3bfeMPmrDSoQQpWHzAJhH9belnRRFKbKvtO1Les4bN24wYMAAnnzySdNt6OKEh4eTnJxs+rp2zca3X4SozAyG/ADo7m6cAua+zQ18+eWXNG/enKlTp5KQkACAk5MT06ZNIyoqqtgBHuUhV3WdZIdNxheKmhq5L6LCusvZlZfo6rX5625Paqe/zuCcm2XjioQQJWHzAOjj44NGoynSM3f79u0iPXh5/P3979ve398foETnvHHjBn369CEkJITVq1fft1atVounp2ehLyGElRw/nr/+78MPw31G5+/YsYOOHTsyatQorl69Chh/ARw9ejRRUVEsWLDAZn9fFRQSHFeAyjj4xFP3T5yU8u+BtBqVij8C2gGg1evoGiOrgghREdg8ADo5OdGhQwe2b99eaP/27dvpdo/f+ENCQoq0j4yMNLUPCAjA39+/UJucnBx27txZ6JzXr1+nd+/etG/fnk8//RS12uY/DiFEnm3b8rfvMSHz8ePH6d+/P/369ePYsWOm/aGhoRw9epT169fToEEDa1d6X+maX8nSnARAY/DFSzfSpvVYwx8FbgP3uHLsPi2FEPbCLkYBT5kyhdGjR9OxY0dTT1xMTIzpebwxY8ZQp04d04jgyZMn07NnTxYsWMATTzzBTz/9xI4dO9i9ezdg/M0/LCyMefPm0bRpU5o2bcq8efNwdXVl5Ejjh++NGzfo3bs39evXZ/HixaZF3yG/B1EIYUMFA2CBEfsAZ86cYdasWXz77beF9rdr144FCxbQr1+/8qjwgQzqNBId15pee+e+iBpnG1ZkHXsbtEGvUqNRDPSIlgAoREVgFwFwxIgRxMfHM3v2bG7evElwcDCbN282/eYeExNTqHeuW7dubNy4kbfeeouZM2fSuHFjNm3aZJoDEGDatGlkZmYyadIk00TQkZGReHh4AMYew4sXL3Lx4kXq1q1bqB5FUcrhuxZC3FNKCuzZY9xu0gQaGW+ZRkVF8e6775qme8rToEED5s6dy9NPP21XPfkp3t9jUBlnC3DV9cDF0MnGFVlHirM7x2s1o8ONP2kWH0OtlDtcsHVRQoj7sosACDBp0iQmTZpU7Hu///57kX3Dhg1j2LB7zzqvUqmYNWsWs2bNKvb9cePGMW7cODMqFUJY3W+/gc44WTL9+xMdHc3s2bP5/PPPCy3d5ufnx5tvvskLL7yA1ow5Aq2qBWS5G0fFqhV3qudOsHFB1rUroB0dbvwJwENXjnFBY7k1m4UQlmc3AVAIUXXFxMQQFxdnel3v88+pCVwDph09yjdNmxYKfl5eXrz00kvMmDEDV1fXcqmxNOvu7j+1Hx7Lf10jdyIO3HuO0crgj4D2hO3ZAEDvy0f4tGlXG1dkP0q7ZrOl154WojgSAIUQNhUTE0PzwBZkZWaY9v0OfAl8BuTu21fkmOTkZJYsjeCFF16w+j+UZq3ROxwIMm666ENw1ZdtDeOK4EStZiS4eFIjM4Ve0UfRNupo65Jsztz1nS259rQQ9yIBUAhhU3FxcWRlZuD92FRUGkccDv/Ew9fPYSjQRuWoxTWoN25BvVA7uVp8fd77Ke0avZluR0n2XW+sW++Gd86kSjPn3/3o1Rr+17gzT57egXtOJj3jYtho66JszJz1ncvz/21RtUkAFELYhbST28mOOVlon9bBCW3HJ/DsPBiNi23n3SzJGr061S1StN+ZXnvFPYnGs7q1S7Mb25qF8OTpHQAMir1Q5QNgntKs7yxEeZEAKISwCUVR2LNnD//+978BCoU/byAM+GHsUuJ9bDuPX0kp6IlzXISiSjfuOAnOHm2hCs0Xv6thW9IdnXHLzaL/rUu2n2hWCHFP8vdTCFGudDod33zzDSEhIfTo0YO9e/ea3nN09eID4ArwWK1mFSb8ASQ7fEW2xjgKVp3tBb/YuCAbyHbUsvPupNA1czLpbuN6hBD3JgFQCFEuUlNT+c9//kOzZs0YPnw4Bw4cKPS+Z9cnmd5pMFMAdyCyAo0izVKfJNnha+MLRYP75UGQbduabCWyWYhpe7AN6xBC3J8EQCGEVV2/fp3p06dTr149wsLCiI6ONr3XunVrZs+eDYBr8+4MicrvDdzarPilIO2NnkTinBaDyjgxdTXdKBzTa9m4Ktv5tXEnctUa4G4AlIn1hbBLEgCFEFZx6NAhxowZQ8OGDVmwYAHJycmm9wYMGMD27ds5fvw4gwYNAqB+ahxtb0YBcMa3EZe96xZ7XnuioOOO0/voVQkAOOtb46kbauOqbCvF2Z199VsD0BBoE3/NpvUIIYonAVAIYTFZWVmsX7+eLl260LlzZz7//HN0d1f0cHJyYvz48Zw6dYotW7bwyCOPoFLlT4/yzwJryP6/oJ7lXrs5Eh0/JVtzBgCN4o1PzjRU8rHKz4E9TNtDLx+2YSVCiHuRTyohRJldvXqV8PBw6tWrx9ixYzl48KDpverVqzNjxgyuXr3KJ598QnBwcLHneOLKUdP2LwUChL1K1+wk1eEn4wvFgZrZb6Khmm2LshNbm3cj++5t4H9GH0Vj0D/gCCFEeZNpYESF8/dlw0pCllayPIPBwP/+9z9WrFjBf//7XwwGQ6H327Zty0svvcTIkSMfuFxbKyAo8f+3d+fhTVXpA8e/SZN03+heNhHLvijFgYIKyrSALBVHRWEAd8FBxeXnyLgxyuKKDAooooLDUkYRAcGyKJSlRbYiO1aWArWlpbbpnjTJ+f1RCIQWaIE2Xd7P89wnNyfn3vvekxDe3txzTjoAO8PbcMo3pFIxVHWKrarWvxST5ney9dPtzxuVjsZVtb4u+64P8ty8WBN8I4MyUggqKaDn8d1svDHS2WEJIS4gCaCoUyqaNqwyZGql6yctLY25c+fyxRdfOHToANDr9dx3332MHTuWqKgoh594L+fRC9aXtr/zivWvdoqt68GiySTL9d8oTVk3X09LNF7WvjUeR233beO2DMpIAWDI/vWSAApRy0gCKOqUC6cNk6mVak5paSmrVq1izpw5rFq1qtzVvsaNGzN69Ggef/xxQkNDq7RvTWkp59K4Ep2B5e2uPG/u1UyxBVB8dAfGTfOrFJ/DcSki0/AWVk0OAK7WdgSUNoyp3qpqbXAL/gQaAf1/S+TNkgLy3LycHZYQ4ixJAEWdJFMr1YyUlBS++OIL5s2bR0ZGhsNrGo2G6OhonnjiCWJjY9Hr9Vd1DN+NGwk8ux7fKqpKSUJVPwel19AjVWEly/AupdrjAOhsYQSZX0XD1Z13fWd20TEfeBZws5i5Z/96vo4c5OywhBBnSQIohHCQnZ3NN998w/z589myZUu515s2bcqjjz7KI488QvPm1z5TR8DSpfb1bzpGX/P+qoUGjIFxlLjsBECrvAg2T8AFXycHVrt9TlkCCPDQr6v5ustAqORtAUKI6iUJoBCCkpISVq5cyfz581m5ciWlpaUOr+t0OmJjY3n88ceJjo7GxcXl+hz4t9/wTUoC4IRXIxKbd7o++72OFAr6QYn39rMFOoLMr6JXjZ0bWB2wD9gR1JyuWam0zTrOLX8cJrlxG2eHJYRAEkAhGiybzcbmzZuZP38+33zzDbm5ueXqtG/fnocffpiRI0cSHBx8/YOYMcO++lXr21Ca2jcyVXF4IoSffaK0BJr/DzdbR6fGVJfMj+hB16xUAB7euUISQCFqCUkAhWhAbDYbW7du5X//+x/ffvstaWlp5eqEhYUxbNgwRowYQadOnSrdk7fK8vPhq68AKALiIrpVz3GugVH3LcXhW+3PA0qfxdPW04kR1T3LWtzC68krCSjO4+7Dm5mS9wgZPoFX3lAIUa0kARSinqtM0ufp6cnf/vY3/v73v3PXXXddv594L+fLL8uSQGABkOvqiWv1H7XScnVxGPXnewx7Zw/By+OvToyobirRGVhwc3+eTVqM3mZlZPIPvNfrYWeHJUSDJwmgEPWQ1WolMTGR7777jm+++abCpE+v1xMTE8OwYcOIjY3F09Oz5gI0m+GDD+xPp1+mak1TKIy6+Rj1i88X/gSeN/aCy49nLS7hv10GMPqXJRhsFobtjmdm9wcocJXGFMKZJAEUop4oKChgzZo1LF++nB9++IHs7Oxydc4lfffffz+xsbH4+Tlp6rL58+HUKQBy77iDfRs3UrXRA6uHwkau7ivy9Od7Jnuc7EXRpgS40YmB1XFZXo1Y1q439+9bh19JASOSVzKr+/3ODkuIBk0SQCeqqSnNZOq0mlPTbf3HH3+wYsUKli1bxk8//YTZbC5X51zSd9dddxEZGYm3tzcAR48erdbYLsligXfesT89/cgjsHHj9T3GVVCUckb/EUW687H4m59Ce9qLIhKcGFn9MDPqfu7d/zMuysbj25Yyr8tAigzuzg5LiAZLEkAnqakpzWTqtJpTE21tsVjYtm0bq1evZtWqVezYsaPCel5eXvTr14/BgwczaNAg8vLyas/nYN48SCmbIozevSns5PyhX6wUkGWYiMllX1mB0tKo9Gm8rf0oYL1zg6snjjVqzPK2dzDkwAYCivMYkbySz7rd5+ywhGiwJAF0kpqa0kymTqs51dXWJ0+eZPXq1axevZp169ZVOFwLlE3HNnjwYAYPHsydd96Jq+v5LhVHjx6tHZ+D4mKYMOH884kTr89+r0GpJo0sw0RKtWWzhGiUK4Hml/Gw1b5eyXXdJ1FDiT2QgBbF00nfsLhTDLnuPs4OS4gGSRJAJ6upKc1k6rSac61tXVRUxMaNG+1J38GDBy9Z9+abb2bw4MHExsZyyy23XHHIFqd/Dj7+2H7vH4MGQc+esGuX08Ip0m7jjOEDlKbsyqhW+RJsegNX1dppMdVnRwKb8m3HPjywdx2+pkKeSVzM232ecHZYQjRIkgAKUQts27aNpUuXsn79erZt21ZuJo5z/P39iY6Opm/fvsTExNCkSZMajvQa/PEHvP122bpGA5MmOS2Uss4eCzHqF9rL9LZmBJlfR6/CnBZXQ/Dh7X9n0MFNuFtMjNi1kgU39+doQB36HAtRT0gCKEQNs5WaMP1xCNOJvRT9vg2AMWPGVFhXq9XSrVs3+vbtS79+/ejatWvNjNFXHV56CQoKytafego6Omk2DR/ICZ2FWZ9iL/Kw9iDAPA6tjPNS7U57BzL7L0N4LjEOg83CxDUzGfbgJJkjWIgaJgmgENXMWpiDKe0gprRDZY8ZKWC1XLJ+q1at6N27N9HR0fTp0wd/f/8ajLaarF4NixaVrQcEOO3qn8kvBcaA2f1s8qc0+FlG4mO5Dw2SgNSUWd3v497962lqPE2PE3u4d//PfNehj7PDEqJBkQRQiOus9M8/MKf/Zk/6LLnpV9xmyJAh3HffffTu3Zvw8PAr1q9TcnLgscfOP3/nHWjUqEZDsGIkRz+HwpvO9+h1UQEEml/Ezeb8XsgNTYnejdejxzD32wkAvPHT5yQ268zvzg1LiAZFEkAhrpJSimPHjrFjxw527NjB+vVlyUX2ivcuu53OLwzXph1wa94JrZs3Wd9O4LXXXqNLly41EXbNUgrGjoVzM5HExDgmg9V9eBSFLuvJ0c/Bpsmzl7sWdiZI+woueNdYLMLRhpZdWd72DgYf3IhfSQEfrPqIIXIVUIgaIwmgEJWglOLkyZP2ZO/ckpOTc/kNXfS4hkXg2rht2RLeBhfP87NvmDLq+TWPzz6DhWc7Wvj6wpw5NXavl1lzjBz955S47LGXaSyuqB9M+HV+GJdQSf6c7fXoMdx6cj9hBdnclvor49y8eM3ZQQnRQEgCKMRFSkpKOHDgAHv27LEvv/76a6Vn+HBt3hmPG7vi2rgthpCWaHT6ao64lvrlF3juufPPP/sMmlZuDMJrYSWXXP18ClzWgMZmL/ew3Ibrvg7k7P4UTWe53682MLp788LAF1gQ9xpaFOMPbyERcN7AQEI0HJIAigbLZrNx8uRJ9u3bZ0/y9uzZw2+//YbVar3i9sHBwdx666107dqVrl27otfr6devH/69H5ExF48cKRvn79zUdM8+C0OHVushreSRp/uOfN0PKE2JvdzFFkKj0ifxsHWjwCKzetQ2Sc07M/X24by0aT5aYCHwtz9P8XtD/zckRDWTBFA4VVXnzr3coMiXoqxlY+rNnTuX2bNnc/z4cY4fP05qaiolJSVX2LpMQEAAERERtG/fnrZt29KuXTuCg4MdBl6+mtjqo7Tt2wm47z7csrIAyI+M5Pdhw1CXGPD5WtvNSg55uhXk61agNMX2co1yx9cyFB/LYDQYrukYonp9EjWUdqePcvdvifgCC9d+xv1hrUj1r2cdooSoRSQBFE5ztXPnVkTZrFjysrDkZmDJTceSk05p9klKs09hyc0A4OOPP67iXjWAAiA7O5vs7Gy2bt16zbHWZ2nbtlHUvTuNVVm77Qdu27mT3O7dr//BQiE3cAElbsmguWBYHaXDy9oXv9IHcaEeDKHTEGg0vDjgBQKzjvGXnHSCS/KJWzievz84kSOVnLpQCFE1kgAKp6nq3Lm2UhNFhzZSsGslhfvXU7BnbVmyl5uOxZgJtiv/bGun0eLiHYDON6Rs8Q8vW3yD0WhdKD66A+Om+VWaO/fcNg3W3r0EDR6M4Wzyl+oVwKi+Y3Hz8if0MptVpd1sFFPkkoix9f+gK5Sw/fyLSoeXNRpfywPoVNA1nIhwhmKDG8NuHcLSNTPpCIQVZLN44SuMHvIvdjRp7+zwhKh3JAEUTqcPaIoh5EZshUYseZlY8rKw5mVhyT9T9piXhSU/C1thrn2b/B3LKrVvjd4NrYcvVuNpvG65G7fmndEHNEXvF3bZzhml2SftsVX2fr5z2zRIixfDE09gyM8H4LhXAMNGvM8Zn2Bcr7DpldpNYcWk3U+By08UuWwpu7/vgl90tcoLL0s/vK0DJPGr43IN7twJrG/UhI5/niKwyMiiRf/i332eZP4td8tsIUJcR5IAimpXXFxMeno66enpZGRk2Nd//fVXALKWvI212HjZ2TEuR6N3Q+cXis4/DL1fWNm6Xxj6gCa4eAdQeCCB7B8+xKtTjHTOuN4KCsp6+n75pb1oO/DY3c+R5xN81btVmCnW7qbYZStFLr9g0xjLV8oCH839+HoORYvbVR9L1C7ZwH19n+bLxDh6pu5Bb7Myce0sup/cx4S/PskZT/lZX4jrQRJAcVVMJhNnzpzhzJkzZGVlcfr0aXtid3GyZzRW8J/3BawF2Vc4mgYXr0bofIJQGg3mtIP43jYct+Y3o/cPRevh59AZQ9QAm61sfL9XXjk/yDPwZ9++9F69Gh93nyte+XOggVLDKUy6PZRo91Ci3efQk9deTXngab0dbYo/eYvj8BjVE62nJH/1TZ7Bg5EPvM0/N8zlye1LARh4aBO3HU9mSu9H+LbjX7Fq6+ic2ELUEpIANnBKKZTFjM1UiCopxFqSh+mPwwB8+eWXLFiwwJ7knUv4zpw5Q/7Zn/quB43Bvew+PJ8gXHyC0HkHlj2eXVy8AtC4lH1UC/avJzvtIO4tb5Wrec5gtcKyZWXTuW2/4P47T0+YMYPjHTpQtHo1PlfYjY0izNrfMWl/o6jlFvg/yPb4oMK6GuWKuy0SD2tP3K3d0eJKQaEM51LfWbUuTL7rMX4Na8XENTPxL8nHr6SAd+M/5slt3/Gfng+xss3tzg5TiDpLEsA6rrS0lIKCAvLz8+2PRqORnJwccnNz2b9/PwDGpP+h0bpgMxViKyk8+1iAzVRwyZ9eZ8yYcU2xeXp6EhYW5rCEhoba13Nzcxk6dCghD02RZK62S00tu+L35Zfw+0WzlwweDFOnQsuWcNFQLworFk0GpZqTlGpPUKo5gVl7hFLNKdCUdRapqKOui/LHzRqJhzUKN9vNaKt2PVHUIyvb3s7WZh15/efPuedAAgAt/0xj+ooP+Nf6L5kf1oqZQPHldyOEuEitSQBnzpzJ+++/T3p6Ou3bt2fatGncfvul/7pbsmQJr7/+OkeOHKFly5ZMmjSJIUOG2F9XSvHvf/+b2bNnk5OTQ7du3ZgxYwbt25/vTZaTk8Ozzz7L8uXLARg8eDAff/wxfn5+5Y53tZRSlJSUUFxcTFFREcXFxRQXF7N3714ATKcOYDVmYrOYUBYzylyMzVyMKi0pezQXYystQZmKsRaVTTs2cOBAzGYz+fn5mM8NtHsFxb8lXvO5aDQaAgICCAwMJDAwkKCgIPt6SEiIQ3IXFhaGl5fXZfe36xLjwgnnM1hK6ZKRQheg9ciRcPYPCQedOqHefRdj7+6k5qZy/NAyNh3dBH0hJ/hLbK55ZxO9StzbWQyutk546nvgZuuMTjVBg/ysL8pke/oxbtD/sfDm/ry4aT7dTu4DILTgT15K2cpLwIHl7/JTmzvY1rQ9yeFtKHD1cG7QQtRytSIBXLx4MePGjWPmzJn07NmTzz77jP79+3PgwAGaNWtWrn5SUhJDhw7l7bffZsiQISxdupQHHniAzZs3061bNwDee+89pk6dyty5c2nVqhUTJ04kOjqaw4cP4+1dNgfosGHDOHXqFPHx8QA8+eSTjBgxghUrVlQp/v79+1NaWmpP7i5crjTQcM5Ps6t0LID09PQqb3MxjasnWldPtG7nHr3QunqhdfNE2SwU7FrJ+++/T/fu3e1Jnr+/Py4uct9NfeJtKqSxMZMmxkwaG/+gaW4KjY3HaFT0B0UGK8Y2sES3nz+7w2lPyPSE080DON2sEaddcsjccQ+mX0yOO40CE3sqPiCA0mFQN2CwtcLV1grLbzkYl8zDf+SjciVYXNa2ph0Y+tAUok7s4eGdK+jz+zZ0qmy6v3Y56bRLWgxJYNVoORzUnINBN/B7YDNSAppx0i+EDK8AjG6X/8NUiIaiViSAU6dO5bHHHuPxxx8HYNq0aaxevZpZs2YxZcqUcvWnTZtGdHQ048ePB2D8+PEkJCQwbdo0Fi1ahFKKadOm8eqrr3LvvfcCMG/ePEJCQli4cCFPPfUUBw8eJD4+nq1bt9qTxs8//5yoqCgOHz5M69atKx1/YuK1X12rLI0GfLw88fbyxN3dFQ93t7JHN1c8PNxwd3PFy8sDb08PvLw8KC4q4Lsl3+J7SzSujcLQGQzo9AY02ktcXVFgycvEmAWNvf5Ek7ePbKPizBFVYV2lypcrZauoKlxU9+Spk3S6AXyKNqDL2n9R7Yv2ezZcs+13CluAd+F69Fl7URfXO3+0ss2UwqSOUHQjeBf9hEvW7ouuK11Y91yJwmJySVEAABMGSURBVKyOUNQSvIvWosvcWcExHLcFMKljFN8E3kWrccnc7vCapoJ2ArAUZdO4Fez4+VP+OBRuP77VZsWirFixYVFWLKrs0arKyi2UlVmVDYuyYLGUYrGYsFpKsVjMWCxmSiwllFhNFFtN5BYZ6TwcNK7jseYpLFoLpVorRp3idCBsDIfiSk9ZnA2FV+q4AygX9CocvWqK3tbs7GNT9KqJw8wcBSXry73dQlySRkNS884kNe9McH42A36eQ8yhTURdUMVF2WiXeYx2mcfKbV6iM/CHwYNMoGTtLAq8g8h39SDf1ZN8Vw9MOgOlLjrMLnrMLnpMOj2lWh15p4+QA/ilHcSlOA+FBptGg9JoURqwabRlz9GgNBpsGi2FOenkAAHZJ7n4n1e5j/zZjmzm7FOEAtpi+VFbVC+nJ4Bms5mdO3fyyiuvOJTHxMRcMrFKSkri+eefdyjr27cv06ZNA+DYsWNkZGQQExNjf93V1ZVevXqRmJjIU089RVJSEr6+vvbkD6B79+74+vqSmJhYYQJoMpkwmc5f7SjXu9UF0FPWqvoKnusqWL/4uaGCRX9+XWnASCFGCitsm3JcgccB1lauPkAQ0BKGpU6B1MpvdlUeBPi+8vWDgE4AlRsHEIAQoCNAFa7s2rdZeRXb/Fj5bQBawlPZn5eNf1FdDEBTuOSdUtazSyVp0OBn8MPf1R9/gz/B7sGEuoeichVzp83F9+aRGLxao8HxirGVIqz85lB2bhxAU8bv2MyVm5rvarapyWPJNjWzTSow1T+MCUCb3o/Q01REl6xUIs8cp3VuBtqKNrKYCbSYCQT44zBwuFKx2a37rGr1AX74sMqb/LZ/P3kdO1b9WKJS8vLygIovYjQYysnS0tIUoLZs2eJQPmnSJNWqVasKt9Hr9WrBggUOZQsWLFAGg0EppdSWLVsUoNLS0hzqPPHEEyomJsa+/4iIiHL7joiIUJMnT67wuG+++ea5S1OyyCKLLLLIIksdX06ePFm5ZKUecvoVwHMuHsdNKXXZsd0qU/9KdSra/+WOO378eF544QX789zcXJo3b86JEyfw9fW9ZKzi2uTl5dG0aVNOnjyJj8+VBhgRV0vaueZIW9cMaeeaURfbWSlFfn4+4eHhzg7FaZyeAAYGBuLi4kJGRoZDeWZmJiEhIRVuExoaetn6oaFlM49mZGQQFhZ2yTqnT58ut++srKxLHtfV1RVX1/LDUfj6+taZD31d5uPjI+1cA6Sda460dc2Qdq4Zda2dG/qFmwpvkahJBoOByMhI1q51vEdt7dq19OjRo8JtoqKiytVfs2aNvX6LFi0IDQ11qGM2m0lISLDXiYqKwmg0sm3bNnudX375BaPReMnjCiGEEELUB06/AgjwwgsvMGLECLp27UpUVBSzZ8/mxIkTjB49GoCRI0fSuHFje4/g5557jjvuuIN3332X2NhYli1bxrp169i8eTNQ9tPuuHHjmDx5MhEREURERDB58mQ8PDwYNmwYAG3btqVfv3488cQTfPZZ2U29Tz75JAMHDqxSD2AhhBBCiLrGZcKECROcHUSHDh0ICAhg8uTJfPDBBxQXF/Pf//6Xzp07A/Cf//wHnU7HPffcA0DTpk1p164dU6dOZfLkyZw4cYJZs2YRHR1t32fPnj0xmUy89dZbTJ8+HR8fH+Li4mjatKm9zt13382uXbt47bXX+O9//8sdd9zBnDlzcHOr/NyiLi4u9O7dG52uVuTS9Za0c82Qdq450tY1Q9q5Zkg71z0apRpyH2ghhBBCiIbH6fcACiGEEEKImiUJoBBCCCFEAyMJoBBCCCFEAyMJoBBCCCFEAyMJ4FWaOXMmLVq0wM3NjcjISDZt2uTskOqcjRs3MmjQIMLDw9FoNHz/veOcwEopJkyYQHh4OO7u7vTu3Zv9+/c71MnJyWHEiBH4+vri6+vLiBEjyM3NrcnTqNWmTJnCrbfeire3N8HBwdxzzz0cPuw496nJZOKZZ54hMDAQT09PBg8ezKlTpxzqnDhxgkGDBuHp6UlgYCDPPvssZrO5Jk+lVps1axadOnWyD4QbFRXFjz+enxNa2rh6TJkyxT7s1znS1tfHhAkT0Gg0Dsu5SRZAvp/rA0kAr8LixYsZN24cr776KsnJydx+++3079+fEydOODu0OqWwsJDOnTvzySefVPj6e++9x9SpU/nkk0/Yvn07oaGhREdHk5+fb68zbNgwdu/eTXx8PPHx8ezevZsRI0bU1CnUegkJCfzjH/9g69atrF27FovFQkxMDIWFhfY648aNY+nSpcTFxbF582YKCgoYOHAgVqsVAKvVyoABAygsLGTz5s3ExcWxZMkSXnzxRWedVq3TpEkT3nnnHXbs2MGOHTu46667iI2Ntf+HKG18/W3fvp3Zs2fTqVMnh3Jp6+unffv2pKen25e9e/faX5Pv53rAabMQ12F/+ctf1OjRox3K2rRpo1555RUnRVT3AWrp0qX25zabTYWGhqp33nnHXlZSUqJ8fX3Vp59+qpRS6sCBAwpQW7dutddJSkpSgDp06FDNBV+HZGZmKkAlJCQopZTKzc1Ver1excXF2eukpaUprVar4uPjlVJKrVq1Smm1WpWWlmavs2jRIuXq6qqMRmPNnkAd4u/vr+bMmSNtXA3y8/NVRESEWrt2rerVq5d67rnnlFLyeb6e3nzzTdW5c+cKX5Pv5/pBrgBWkdlsZufOncTExDiUx8TEkJiY6KSo6p9jx46RkZHh0M6urq706tXL3s5JSUn4+vrSrVs3e53u3bvj6+sr78UlGI1GABo1agTAzp07KS0tdWjn8PBwOnTo4NDOHTp0cJg0vW/fvphMJnbu3FmD0dcNVquVuLg4CgsLiYqKkjauBv/4xz8YMGAAf/3rXx3Kpa2vr5SUFMLDw2nRogUPPvggR48eBeT7ub6QIbur6MyZM1itVkJCQhzKQ0JCyMjIcFJU9c+5tqyonVNTU+11goODy20bHBws70UFlFK88MIL3HbbbXTo0AEoa0ODwYC/v79D3Qs/zxkZGeXeB39/fwwGg7TzBfbu3UtUVBQlJSV4eXmxdOlS2rVrx+7du6WNr6O4uDh27drF9u3by70mn+frp1u3bnz99de0atWK06dPM3HiRHr06MH+/fvl+7mekATwKmk0GofnSqlyZeLaXamdK2pzeS8qNnbsWPbs2WOfM/typJ2rrnXr1uzevZvc3FyWLFnCqFGjSEhIuGR9aeOqO3nyJM899xxr1qyp0pSd0tZV179/f/t6x44diYqKomXLlsybN4/u3bsD8v1c18lPwFUUGBiIi4tLub9gMjMzy/01JK7eud5ml2vn0NBQTp8+XW7brKwseS8u8swzz7B8+XLWr19PkyZN7OWhoaGYzWZycnIc6l/czhe/Dzk5OZSWlko7X8BgMHDTTTfRtWtXpkyZQufOnfnPf/4jbXwd7dy5k8zMTCIjI9HpdOh0OhISEpg+fTo6nY6QkBBp62ri6elJx44dSUlJke/nekISwCoyGAxERkaydu1ah/K1a9fSo0cPJ0VV/7Ro0YLQ0FCHdjabzSQkJNjbOSoqCqPRyLZt2+x1fvnlF4xGo7wXZymlGDt2LN999x0///wzLVq0cHg9MjISvV7v0M7p6ens27fPoZ337dtHenq6vc6aNWtwdXUlMjKyZk6kDlJKYTKZpI2voz59+rB37152795tX7p27crw4cPt69LW1cNkMnHw4EHCwsLk+7m+cEbPk7ouLi5O6fV69cUXX6gDBw6ocePGKU9PT3X8+HFnh1an5Ofnq+TkZJWcnKwANXXqVJWcnKxSU1OVUkq98847ytfXV3333Xdq79696qGHHlJhYWEqLy/Pvo9+/fqpTp06qaSkJJWUlKQ6duyoBg4c6KxTqnXGjBmjfH191YYNG1R6erp9KSoqstcZPXq0atKkiVq3bp3atWuXuuuuu1Tnzp2VxWJRSillsVhUhw4dVJ8+fdSuXbvUunXrVJMmTdTYsWOddVq1zvjx49XGjRvVsWPH1J49e9S//vUvpdVq1Zo1a5RS0sbV6cJewEpJW18vL774otqwYYM6evSo2rp1qxo4cKDy9va2/z8n3891nySAV2nGjBmqefPmymAwqC5dutiH1RCVt379egWUW0aNGqWUKhtq4M0331ShoaHK1dVV3XHHHWrv3r0O+8jOzlbDhw9X3t7eytvbWw0fPlzl5OQ44Wxqp4raF1BfffWVvU5xcbEaO3asatSokXJ3d1cDBw5UJ06ccNhPamqqGjBggHJ3d1eNGjVSY8eOVSUlJTV8NrXXo48+av8+CAoKUn369LEnf0pJG1enixNAaevrY+jQoSosLEzp9XoVHh6u7r33XrV//3776/L9XPdplFLKOdcehRBCCCGEM8g9gEIIIYQQDYwkgEIIIYQQDYwkgEIIIYQQDYwkgEIIIYQQDYwkgEIIIYQQDYwkgEIIIYQQDYwkgEIIIYQQDYwkgEKIBumVV16xT2pfnUJDQ/n000+r/ThCCFEVkgAKIZxCo9Fcdnn44YcrtZ8HH3yQBx98sHqDvQZ79+5l1KhRla4fHx+PRqOhpKSkGqMSQjR0OmcHIIRomNLT0+3rixcv5o033uDw4cP2Mnd3d2eEdd0FBQU5OwQhhChHrgAKIZwiNDTUvvj6+qLRaMqVASQnJ9OrVy/c3NwICgri6aefpri4GCj7GXfx4sUsXrzYfuVw69atADz//PNERETg7u5Oy5Yteeutt7BYLJWO79yVuPj4eDp06ICbmxs9evTg4MGDDvXi4uJo27YtBoOBFi1aMH369HLnee4n4JKSEjQaDfPmzWPgwIF4eHjQunVrfvzxRwAOHTpE//79gbIEWKPRMHr0aAAWLVpE+/btcXNzIzAwkJiYGEwmU1WbXQghAEkAhRC1WF5eHn379iU8PJydO3eycOFCVq5cyfPPPw/Aa6+9RmxsLLGxsaSnp5Oenk5kZCQAfn5+fP311xw8eJAPP/yQ6dOnM3PmzCrH8PLLLzN9+nS2bduGj48PsbGxWK1WABITExk+fDijRo1i3759vPrqq7z88svExcVddp9vvvkmo0aNYs+ePdx5550MGzaMvLw8IiIiWLhwIQDHjx8nPT2d9957j9TUVEaMGMHTTz/N4cOH+fnnnxk0aFCVz0UIIeyUEEI42VdffaV8fX3LlU+fPl0FBQWp4uJie9mSJUuUi4uL+vPPP5VSSg0dOlQNHTr0isd46623VM+ePe3P//nPf6pu3bpdsv6PP/6oAPX999/by06fPq0MBoNatmyZUkqpe++9Vw0aNMhhu2eeeUZ16dLF/jwkJETNmjVLKaVUcXGxAtTEiRPtr2dnZytArV+/3uG4F57zli1blEajUenp6Vc8TyGEqAy5AiiEqLUOHjxIZGQkbm5u9rKePXtitVpJSUm57LaLFi2iR48ehISE4OXlxaRJkzhx4kSVY4iKirKvBwcHc+ONN9p/Bj548CA9e/Z0qN+zZ08OHTp02X126tTJvt6oUSMMBgOZmZmXrH/rrbdy22230aZNG4YOHcoXX3yB0Wis8rkIIcQ5kgAKIWotpRQajaZcGVCu/EIJCQmMHDmSe+65h1WrVpGcnMxLL72E2Wy+LnGdO/bl4rscvV5fbn82m+2y9Tds2MCKFSto1aoVH330EW3atOHUqVNXEb0QQkgCKISoxdq1a8eOHTschkRJTExEp9MREREBgMFgsN+Td87mzZtp1aoVL7/8MpGRkURERHD8+PGriuFcpxKAzMxMjh49Sps2bezxbd682aF+YmIibdu2vapjQdn5AOXOSavVcvvtt/P222+TnJyM1Wpl+fLlV30cIUTDJgmgEKLWOjd+3qOPPsr+/ftZu3Ytzz//PI899hh+fn4A3HDDDezevZuUlBTOnDmDxWLhpptu4vfff2fJkiUcOXKEDz/8kJUrV15VDG+88QYbNmywj+fXrFkz7r77bgBeeuklVq5cybvvvktKSgpz5sxh9uzZvPTSS1d9zjfccAMAP/zwA1lZWRQWFrJp0ybeffdddu7cSWpqKt988w05OTnXlGgKIRo2SQCFELWWj48P8fHxpKWlERkZyUMPPcSAAQP46KOP7HXGjBlDs2bNuOWWWwgKCmLHjh088MADjBkzhieffJJbbrmF5ORkxo8ff1UxTJkyhaeffpquXbuSm5vL999/j05XNoRqVFQUCxYsYO7cubRv355Jkybx3nvvXdPA1DfeeCOvvvoqzz77LCEhIbz44ov4+fnx008/0a9fP9q0acPbb7/NjBkzuPPOO6/6OEKIhk2jKnPDihBCNDDx8fH079+f4uJih04oQghRH8gVQCGEEEKIBkYSQCGEEEKIBkZ+AhZCCCGEaGDkCqAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAMjCaAQQgghRAPz/2zpSOpHXKPrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Download and save the data file MLEplots.png\n",
    "url = ('https://raw.githubusercontent.com/rickecon/Notebooks/' +\n",
    "       'master/GMM/images/MLEplots.png')\n",
    "image_file = requests.get(url, allow_redirects=True)\n",
    "open('images/MLEplots.png', 'wb').write(image_file.content)\n",
    "Image(\"images/MLEplots.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1.1. Two moments, identity weighting matrix\n",
    "Let's try estimating the parameters $\\mu$ and $\\sigma$ by GMM. What moments should we use? Let's try the mean and variance of the data. These two statistics of the data are defined by:\n",
    "\n",
    "$$ mean(scores_i) = \\frac{1}{N}\\sum_{i=1}^N scores_i $$\n",
    "\n",
    "$$ var(scores_i) = \\frac{1}{N}\\sum_{i=1}^{N} \\left(scores_i - mean(scores_i)\\right)^2 $$\n",
    "\n",
    "So the data moment vector $m(x)$ for GMM is the following.\n",
    "\n",
    "$$ m(scores_i) \\equiv \\begin{bmatrix} mean(scores_i) \\\\ var(scores_i) \\end{bmatrix} $$\n",
    "\n",
    "And the model moment vector $m(x|\\theta)$ for GMM is the following.\n",
    "\n",
    "$$ m(scores_i|\\mu,\\sigma) \\equiv \\begin{bmatrix} mean(scores_i|\\mu,\\sigma) \\\\ var(scores_i|\\mu,\\sigma) \\end{bmatrix} $$\n",
    "\n",
    "Define the error vector as the vector of percent deviations of the model moments from the data moments.\n",
    "\n",
    "$$ e(scores_i|\\mu,\\sigma) \\equiv \\frac{m(scores_i|\\mu,\\sigma) - m(scores_i)}{m(scores_i)} $$\n",
    "\n",
    "The GMM estimator for this moment vector is the following.\n",
    "\n",
    "$$ (\\hat{\\mu}_{GMM},\\hat{\\sigma}_{GMM}) = (\\mu,\\sigma):\\quad \\min_{\\mu,\\sigma} e(scores_i|\\mu,\\sigma)^T \\, W \\, e(scores_i|\\mu,\\sigma) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something interesting to note here is the $1/N$ weighting on our variance estimator. There is less bias in the estimator of the variance by using the weighting $1/(N-1)$ because one degree of freedom is used in calculating the mean used in the variance calculation. However, in GMM when many moments are used that might have differing degrees of freedom restrictions, it is important to have the same weighting for each moment. So we use $1/N$ in all cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download and save the data file Econ381totpts.txt\n",
    "url = ('https://raw.githubusercontent.com/rickecon/Notebooks/' +\n",
    "       'master/GMM/data/Econ381totpts.txt')\n",
    "data_file = requests.get(url, allow_redirects=True)\n",
    "open('data/Econ381totpts.txt', 'wb').write(data_file.content)\n",
    "\n",
    "# Load the data as a NumPy array\n",
    "pts = np.loadtxt('data/Econ381totpts.txt')\n",
    "#pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function that generates values of a normal pdf\n",
    "def trunc_norm_pdf(xvals, mu, sigma, cut_lb, cut_ub):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    Generate pdf values from the normal pdf with mean mu and standard\n",
    "    deviation sigma. If the cutoff is given, then the PDF values are\n",
    "    inflated upward to reflect the zero probability on values above the\n",
    "    cutoff. If there is no cutoff given, this function does the same\n",
    "    thing as sp.stats.norm.pdf(x, loc=mu, scale=sigma).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, values of the normally distributed random\n",
    "             variable\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    prob_notcut = scalar \n",
    "    pdf_vals = (N,) vector, normal PDF values for mu and sigma\n",
    "               corresponding to xvals data\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: pdf_vals\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    if cut_ub == 'None' and cut_lb == 'None':\n",
    "        prob_notcut = 1.0\n",
    "    elif cut_ub == 'None' and cut_lb != 'None':\n",
    "        prob_notcut = 1.0 - sts.norm.cdf(cut_lb, loc=mu, scale=sigma)\n",
    "    elif cut_ub != 'None' and cut_lb == 'None':\n",
    "        prob_notcut = sts.norm.cdf(cut_ub, loc=mu, scale=sigma)\n",
    "    elif cut_ub != 'None' and cut_lb != 'None':\n",
    "        prob_notcut = (sts.norm.cdf(cut_ub, loc=mu, scale=sigma) -\n",
    "                       sts.norm.cdf(cut_lb, loc=mu, scale=sigma))\n",
    "            \n",
    "    pdf_vals    = ((1/(sigma * np.sqrt(2 * np.pi)) *\n",
    "                    np.exp( - (xvals - mu)**2 / (2 * sigma**2))) /\n",
    "                    prob_notcut)\n",
    "    \n",
    "    return pdf_vals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define a criterion function that takes as inputs the parameters and the estimator for the weighting matrix $\\hat{W}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_moments(xvals):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the two data moments for GMM\n",
    "    (mean(data), variance(data)).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals = (N,) vector, test scores data\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data = scalar, mean value of test scores data\n",
    "    var_data  = scalar > 0, variance of test scores data\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: mean_data, var_data\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mean_data = xvals.mean()\n",
    "    var_data = xvals.var()\n",
    "    \n",
    "    return mean_data, var_data\n",
    "\n",
    "\n",
    "def model_moments(mu, sigma, cut_lb, cut_ub):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the two model moments for GMM\n",
    "    (mean(model data), variance(model data)).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        trunc_norm_pdf()\n",
    "        xfx()\n",
    "        x2fx()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_model = scalar, mean value of test scores from model\n",
    "    m_m_err    = scalar > 0, estimated error in the computation of the\n",
    "                 integral for the mean of the distribution\n",
    "    var_model  = scalar > 0, variance of test scores from model\n",
    "    v_m_err    = scalar > 0, estimated error in the computation of the\n",
    "                 integral for the variance of the distribution\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: mean_model, var_model\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    xfx = lambda x: x * trunc_norm_pdf(x, mu, sigma, cut_lb, cut_ub)\n",
    "    (mean_model, m_m_err) = intgr.quad(xfx, cut_lb, cut_ub)\n",
    "    x2fx = lambda x: ((x - mean_model) ** 2) * trunc_norm_pdf(x, mu, sigma, cut_lb, cut_ub) \n",
    "    (var_model, v_m_err) = intgr.quad(x2fx, cut_lb, cut_ub)\n",
    "    \n",
    "    return mean_model, var_model\n",
    "\n",
    "\n",
    "def err_vec(xvals, mu, sigma, cut_lb, cut_ub, simple):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the vector of moment errors (in percent\n",
    "    deviation from the data moment vector) for GMM.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, test scores data\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    simple = boolean, =True if errors are simple difference, =False if\n",
    "             errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        data_moments()\n",
    "        model_moments()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data  = scalar, mean value of data\n",
    "    var_data   = scalar > 0, variance of data\n",
    "    moms_data  = (2, 1) matrix, column vector of two data moments\n",
    "    mean_model = scalar, mean value from model\n",
    "    var_model  = scalar > 0, variance from model\n",
    "    moms_model = (2, 1) matrix, column vector of two model moments\n",
    "    err_vec    = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: err_vec\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mean_data, var_data = data_moments(xvals)\n",
    "    moms_data = np.array([[mean_data], [var_data]])\n",
    "    mean_model, var_model = model_moments(mu, sigma, cut_lb, cut_ub)\n",
    "    moms_model = np.array([[mean_model], [var_model]])\n",
    "    if simple:\n",
    "        err_vec = moms_model - moms_data\n",
    "    else:\n",
    "        err_vec = (moms_model - moms_data) / moms_data\n",
    "    \n",
    "    return err_vec\n",
    "\n",
    "\n",
    "def criterion(params, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the GMM weighted sum of squared moment errors\n",
    "    criterion function value given parameter values and an estimate of\n",
    "    the weighting matrix.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    params = (2,) vector, ([mu, sigma])\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    args   = length 3 tuple, (xvals, cutoff, W_hat)\n",
    "    xvals  = (N,) vector, values of the truncated normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    W_hat  = (R, R) matrix, estimate of optimal weighting matrix\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        norm_pdf()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    err        = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    crit_val   = scalar > 0, GMM criterion function value\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: crit_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mu, sigma = params\n",
    "    xvals, cut_lb, cut_ub, W = args\n",
    "    err = err_vec(xvals, mu, sigma, cut_lb, cut_ub, simple=False)\n",
    "    crit_val = err.T @ W @ err\n",
    "    \n",
    "    return crit_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can perform the GMM estimation. Let's start with the identity matrix as our estimate for the optimal weighting matrix $W = I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this takes a little time because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_init = 400\n",
    "sig_init = 60\n",
    "params_init = np.array([mu_init, sig_init])\n",
    "W_hat = np.eye(2)\n",
    "gmm_args = (pts, 0.0, 450.0, W_hat)\n",
    "results = opt.minimize(criterion, params_init, args=(gmm_args), tol=1e-14,\n",
    "                       method='L-BFGS-B', bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_GMM1, sig_GMM1 = results.x\n",
    "print('mu_GMM1=', mu_GMM1, ' sig_GMM1=', sig_GMM1)\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_data, var_data = data_moments(pts)\n",
    "mean_model, var_model = model_moments(mu_GMM1, sig_GMM1, 0.0, 450.0)\n",
    "err1 = err_vec(pts, mu_GMM1, sig_GMM1, 0.0, 450.0, False).reshape(2,)\n",
    "print('Mean of points =', mean_data, ', Variance of points =', var_data)\n",
    "print('Mean of model =', mean_model, ', Variance of model =', var_model)\n",
    "print('Error vector=', err1)\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the PDF implied by these results against the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "num_bins = 30\n",
    "count, bins, ignored = plt.hist(pts, num_bins, density=True,\n",
    "                                edgecolor='k')\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel(r'Total points')\n",
    "plt.ylabel(r'Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "\n",
    "# Plot the estimated GMM PDF\n",
    "dist_pts = np.linspace(0, 450, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_GMM1, sig_GMM1, 0.0, 450),\n",
    "         linewidth=2, color='k', label='1: $\\mu_{GMM1}$,$\\sigma_{GMM1}$')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks just like the ML estimate from the previous notebook. Let's see what the criterion function looks like for different values of $\\mu$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this will take a few minutes because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_vals = np.linspace(60, 700, 50)\n",
    "sig_vals = np.linspace(20, 200, 50)\n",
    "# mu_vals = np.linspace(540, 572, 50)\n",
    "# sig_vals = np.linspace(170, 182, 50)\n",
    "crit_vals = np.zeros((50, 50))\n",
    "for mu_ind in range(50):\n",
    "    for sig_ind in range(50):\n",
    "        crit_vals[mu_ind, sig_ind] = \\\n",
    "            criterion(np.array([mu_vals[mu_ind], sig_vals[sig_ind]]),\n",
    "                      pts, 0.0, 450.0, W_hat)\n",
    "\n",
    "mu_mesh, sig_mesh = np.meshgrid(mu_vals, sig_vals)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot_surface(sig_mesh, mu_mesh, crit_vals, rstride=8,\n",
    "                cstride=1, cmap=cmap1)\n",
    "ax.set_title('Criterion function for values of mu and sigma')\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$\\mu$')\n",
    "ax.set_zlabel(r'Crit. func.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the GMM estimator for the variance-covariance matrix $\\hat{\\Sigma}_{GMM}$ of our GMM estimates $\\hat{\\theta}_{GMM}$ using the equation in Section 4 based on the Jacobian $d(x|\\hat{\\theta}_{GMM})$ of the moment error vector $e(x|\\hat{\\theta}_{GMM})$ from the criterion function at the estimated (optimal) parameter values $\\hat{\\theta}_{GMM}$. We first write a function that computes the Jacobian $d(x|\\hat{\\theta}_{GMM})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Jac_err2(xvals, mu, sigma, cut_lb, cut_ub, simple=False):\n",
    "    '''\n",
    "    This function computes the Jacobian matrix of partial derivatives of the R x 1 moment\n",
    "    error vector e(x|theta) with respect to the K parameters theta_i in the K x 1 parameter vector\n",
    "    theta. The resulting matrix is R x K Jacobian.\n",
    "    '''\n",
    "    Jac_err = np.zeros((2, 2))\n",
    "    h_mu = 1e-8 * mu\n",
    "    h_sig = 1e-8 * sigma\n",
    "    Jac_err[:, 0] = \\\n",
    "        ((err_vec(xvals, mu + h_mu, sigma, cut_lb, cut_ub, simple) -\n",
    "          err_vec(xvals, mu - h_mu, sigma, cut_lb, cut_ub, simple)) / (2 * h_mu)).flatten()\n",
    "    Jac_err[:, 1] = \\\n",
    "        ((err_vec(xvals, mu, sigma + h_sig, cut_lb, cut_ub, simple) -\n",
    "          err_vec(xvals, mu, sigma - h_sig, cut_lb, cut_ub, simple)) / (2 * h_sig)).flatten()\n",
    "    \n",
    "    return Jac_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = pts.shape[0]\n",
    "d_err2 = Jac_err2(pts, mu_GMM1, sig_GMM1, 0.0, 450.0, False)\n",
    "print(d_err2)\n",
    "print(W_hat)\n",
    "SigHat2 = (1 / N) * lin.inv(d_err.T @ W_hat @ d_err)\n",
    "print(SigHat2)\n",
    "print('Std. err. mu_hat=', np.sqrt(SigHat2[0, 0]))\n",
    "print('Std. err. sig_hat=', np.sqrt(SigHat2[1, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1.2. Two moments, two-step optimal weighting matrix\n",
    "Similar to the ML problem, it looks like the criterion function is roughly equal for a specific portion increase of $\\mu$ and $\\sigma$ together. That is, with these two moments probably have a correspondence of values of $\\mu$ and $\\sigma$ that give roughly the same criterion function value. This issue has two possible solutions.\n",
    "\n",
    "1. Maybe we need the two-step variance covariance estimator to calculate a \"more\" optimal weighting matrix $W$.\n",
    "2. Maybe our two moments aren't very good moments for fitting the data.\n",
    "\n",
    "Let's first try the two-step weighting matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_Err_mat2(pts, mu, sigma, cut_lb, cut_ub, simple=False):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the R x N matrix of errors from each\n",
    "    observation for each moment. In this function, we have hard coded\n",
    "    R = 2.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, test scores data\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    simple = boolean, =True if errors are simple difference, =False if\n",
    "             errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        model_moments()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    R          = integer = 2, hard coded number of moments\n",
    "    N          = integer >= R, number of data observations\n",
    "    Err_mat    = (R, N) matrix, error by moment and observation data\n",
    "    mean_model = scalar, mean value from model\n",
    "    var_model  = scalar > 0, variance from model\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: Err_mat\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    R = 2\n",
    "    N = len(pts)\n",
    "    Err_mat = np.zeros((R, N))\n",
    "    mean_model, var_model = model_moments(mu, sigma, cut_lb, cut_ub)\n",
    "    if simple:\n",
    "        Err_mat[0, :] = pts - mean_model\n",
    "        Err_mat[1, :] = ((mean_data - pts) ** 2) - var_model\n",
    "    else:\n",
    "        Err_mat[0, :] = (pts - mean_model) / mean_model\n",
    "        Err_mat[1, :] = (((mean_data - pts) ** 2) - var_model) / var_model\n",
    "    \n",
    "    return Err_mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Err_mat = get_Err_mat2(pts, mu_GMM1, sig_GMM1, 0.0, 450.0, False)\n",
    "VCV2 = (1 / pts.shape[0]) * (Err_mat @ Err_mat.T)\n",
    "print(VCV2)\n",
    "W_hat2 = lin.inv(VCV2)\n",
    "print(W_hat2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this takes a little time because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_init = 400\n",
    "sig_init = 60\n",
    "params_init = np.array([mu_GMM1, sig_GMM1])\n",
    "# W_hat3 = np.array([[1. / VCV2[0, 0], 0.], [0., 1. / VCV2[1, 1]]])\n",
    "gmm_args = (pts, 0.0, 450.0, W_hat2)\n",
    "results = opt.minimize(criterion, params_init, args=(gmm_args),\n",
    "                       method='L-BFGS-B', bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_GMM2, sig_GMM2 = results.x\n",
    "print('mu_GMM2=', mu_GMM2, ' sig_GMM2=', sig_GMM2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at how much smaller (more efficient) the estimated standard errors are in this case with the two-step optimal weighting matrix $\\hat{W}_{2step}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_err2_2 = Jac_err2(pts, mu_GMM2, sig_GMM2, 0.0, 450.0, False)\n",
    "print(d_err2_2)\n",
    "print(W_hat2)\n",
    "SigHat2_2 = (1 / N) * lin.inv(d_err2.T @ W_hat2 @ d_err2)\n",
    "print(SigHat2_2)\n",
    "print('Std. err. mu_hat=', np.sqrt(SigHat2_2[0, 0]))\n",
    "print('Std. err. sig_hat=', np.sqrt(SigHat2_2[1, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1.3. Four moments, identity weighting matrix\n",
    "Using a better weighting matrix didn't improve our estimates or fit very much it did improve the standard errors of our estimates. To get the right fit, we might need to choose different moments. Let's try an overidentified model $K<R$, where we estimate $\\mu$ and $\\sigma$ of the truncated normal distribution $K=2$ using the following four moments $R=4$.\n",
    "\n",
    "1. The percent of observations greater than 430 (between 430 and 450)\n",
    "2. The percent of observations between 320 and 430\n",
    "3. The percent of observations between 220 and 320\n",
    "4. The percent of observations less than 220 (between 0 and 220)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_moments4(xvals):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the four data moments for GMM\n",
    "    (binpct_1, binpct_2, binpct_3, binpct_4).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals = (N,) vector, test scores data\n",
    "\n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "\n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    bpct_1_dat = scalar in [0, 1], percent of observations\n",
    "                 0 <= x < 220\n",
    "    bpct_2_dat = scalar in [0, 1], percent of observations\n",
    "                 220 <= x < 320\n",
    "    bpct_3_dat = scalar in [0, 1], percent of observations\n",
    "                 320 <= x < 430\n",
    "    bpct_4_dat = scalar in [0, 1], percent of observations\n",
    "                 430 <= x <= 450\n",
    "\n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "\n",
    "    RETURNS: bpct_1, bpct_2, bpct_3, bpct_4\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    bpct_1_dat = xvals[xvals < 220].shape[0] / xvals.shape[0]\n",
    "    bpct_2_dat = (xvals[(xvals >=220) & (xvals < 320)].shape[0] /\n",
    "                  xvals.shape[0])\n",
    "    bpct_3_dat = (xvals[(xvals >=320) & (xvals < 430)].shape[0] /\n",
    "                  xvals.shape[0])\n",
    "    bpct_4_dat = xvals[xvals >= 430].shape[0] / xvals.shape[0]\n",
    "\n",
    "    return bpct_1_dat, bpct_2_dat, bpct_3_dat, bpct_4_dat\n",
    "\n",
    "\n",
    "def model_moments4(mu, sigma, cut_lb, cut_ub):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the four model moments for GMM\n",
    "    (binpct_1, binpct_2, binpct_3, binpct_4).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        trunc_norm_pdf()\n",
    "        xfx()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    bpct_1_mod = scalar in [0, 1], percent of model observations in\n",
    "                 bin 1\n",
    "    bp_1_err   = scalar > 0, estimated error in the computation of the\n",
    "                 integral for bpct_1_mod\n",
    "    bpct_2_mod = scalar in [0, 1], percent of model observations in\n",
    "                 bin 2\n",
    "    bp_2_err   = scalar > 0, estimated error in the computation of the\n",
    "                 integral for bpct_2_mod\n",
    "    bpct_3_mod = scalar in [0, 1], percent of model observations in\n",
    "                 bin 3\n",
    "    bp_3_err   = scalar > 0, estimated error in the computation of the\n",
    "                 integral for bpct_3_mod\n",
    "    bpct_4_mod = scalar in [0, 1], percent of model observations in\n",
    "                 bin 4\n",
    "    bp_4_err   = scalar > 0, estimated error in the computation of the\n",
    "                 integral for bpct_4_mod\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: bpct_1_mod, bpct_2_mod, bpct_3_mod, bpct_4_mod\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    xfx = lambda x: trunc_norm_pdf(x, mu, sigma, cut_lb, cut_ub)\n",
    "    (bpct_1_mod, bp_1_err) = intgr.quad(xfx, 0.0, 220)\n",
    "    (bpct_2_mod, bp_2_err) = intgr.quad(xfx, 220, 320)\n",
    "    (bpct_3_mod, bp_3_err) = intgr.quad(xfx, 320, 430)\n",
    "    (bpct_4_mod, bp_4_err) = intgr.quad(xfx, 430, 450)\n",
    "    \n",
    "    return bpct_1_mod, bpct_2_mod, bpct_3_mod, bpct_4_mod\n",
    "\n",
    "\n",
    "def err_vec4(xvals, mu, sigma, cut_lb, cut_ub, simple):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the vector of moment errors (in percent\n",
    "    deviation from the data moment vector) for GMM.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, test scores data\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    simple = boolean, =True if errors are simple difference, =False if\n",
    "             errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        data_moments4()\n",
    "        model_moments4()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data  = scalar, mean value of data\n",
    "    var_data   = scalar > 0, variance of data\n",
    "    moms_data  = (2, 1) matrix, column vector of two data moments\n",
    "    mean_model = scalar, mean value from model\n",
    "    var_model  = scalar > 0, variance from model\n",
    "    moms_model = (2, 1) matrix, column vector of two model moments\n",
    "    err_vec    = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: err_vec\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    bpct_1_dat, bpct_2_dat, bpct_3_dat, bpct_4_dat = \\\n",
    "        data_moments4(xvals)\n",
    "    moms_data = np.array([[bpct_1_dat], [bpct_2_dat], [bpct_3_dat],\n",
    "                          [bpct_4_dat]])\n",
    "    bpct_1_mod, bpct_2_mod, bpct_3_mod, bpct_4_mod = \\\n",
    "        model_moments4(mu, sigma, cut_lb, cut_ub)\n",
    "    moms_model = np.array([[bpct_1_mod], [bpct_2_mod], [bpct_3_mod],\n",
    "                          [bpct_4_mod]])\n",
    "    if simple:\n",
    "        err_vec = moms_model - moms_data\n",
    "    else:\n",
    "        err_vec = (moms_model - moms_data) / moms_data\n",
    "    \n",
    "    return err_vec\n",
    "\n",
    "\n",
    "def criterion4(params, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the GMM weighted sum of squared moment errors\n",
    "    criterion function value given parameter values and an estimate of\n",
    "    the weighting matrix.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    params = (2,) vector, ([mu, sigma])\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    args   = length 3 tuple, (xvals, cutoff, W_hat)\n",
    "    xvals  = (N,) vector, values of the truncated normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    W_hat  = (R, R) matrix, estimate of optimal weighting matrix\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        err_vec4()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    err        = (4, 1) matrix, column vector of four moment error\n",
    "                 functions\n",
    "    crit_val   = scalar > 0, GMM criterion function value\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: crit_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mu, sigma = params\n",
    "    xvals, cut_lb, cut_ub, W = args\n",
    "    err = err_vec4(xvals, mu, sigma, cut_lb, cut_ub, simple=False)\n",
    "    crit_val = err.T @ W @ err \n",
    "    \n",
    "    return crit_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what the model moments would be relative to the data moments with the first GMM estimates $mu\\approx 620$ and $\\sigma\\approx 198$ from the previous 2-moment estimation. Let's also look at the resulting criterion function at those values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = np.array([mu_GMM1, sig_GMM1])\n",
    "print('params=', params)\n",
    "print(data_moments4(pts))\n",
    "print(model_moments4(mu_GMM1, sig_GMM1, 0.0, 450))\n",
    "print(criterion4(params, pts, 0.0, 450.0, np.eye(4)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this takes a little time because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_init = 400\n",
    "sig_init = 70\n",
    "params_init = np.array([mu_init, sig_init])\n",
    "W_hat1_4 = np.eye(4)\n",
    "gmm_args = (pts, 0.0, 450.0, W_hat1_4)\n",
    "results_4 = opt.minimize(criterion4, params_init, args=(gmm_args),\n",
    "                       method='L-BFGS-B', bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_GMM1_4, sig_GMM1_4 = results_4.x\n",
    "print('mu_GMM1_4=', mu_GMM1_4, ' sig_GMM1_4=', sig_GMM1_4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = np.array([mu_GMM1_4, sig_GMM1_4])\n",
    "print('params=', params)\n",
    "print(data_moments4(pts))\n",
    "print(model_moments4(mu_GMM1_4, sig_GMM1_4, 0.0, 450))\n",
    "print(criterion4(params, pts, 0.0, 450.0, W_hat1_4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params2 = np.array([350, 98])\n",
    "print('params=', params)\n",
    "print(data_moments4(pts))\n",
    "print(model_moments4(350, 98, 0.0, 450))\n",
    "print(criterion4(params2, pts, 0.0, 450.0, W_hat1_4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "count, bins, ignored = plt.hist(pts, num_bins, density=True,\n",
    "                                edgecolor='k')\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel(r'Total points')\n",
    "plt.ylabel(r'Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "\n",
    "# Plot the estimated GMM PDF\n",
    "dist_pts = np.linspace(0, 450, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_GMM1_4, sig_GMM1_4, 0.0, 450),\n",
    "         linewidth=2, color='k', label='1: $\\mu_{GMM1}$,$\\sigma_{GMM1}$')\n",
    "\n",
    "# Plot the estimated GMM PDF\n",
    "dist_pts = np.linspace(0, 450, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, 400, 70, 0.0, 450),\n",
    "         linewidth=2, color='g', label='1: $\\mu=400$,$\\sigma=70$')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the estimator of the variance-covariance matrix $\\hat{\\Sigma}$ of the GMM parameter estimator by computing the Jacobian of the error vector. In this case, the Jacobian $d(x|\\theta)$ is $R\\times K = 4\\times 2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Jac_err4(xvals, mu, sigma, cut_lb, cut_ub, simple=False):\n",
    "    '''\n",
    "    This function computes the Jacobian matrix of partial derivatives of the R x 1 moment\n",
    "    error vector e(x|theta) with respect to the K parameters theta_i in the K x 1 parameter vector\n",
    "    theta. The resulting matrix is R x K Jacobian.\n",
    "    '''\n",
    "    Jac_err = np.zeros((4, 2))\n",
    "    h_mu = 1e-8 * mu\n",
    "    h_sig = 1e-8 * sigma\n",
    "    Jac_err[:, 0] = \\\n",
    "        ((err_vec4(xvals, mu + h_mu, sigma, cut_lb, cut_ub, simple) -\n",
    "          err_vec4(xvals, mu - h_mu, sigma, cut_lb, cut_ub, simple)) / (2 * h_mu)).flatten()\n",
    "    Jac_err[:, 1] = \\\n",
    "        ((err_vec4(xvals, mu, sigma + h_sig, cut_lb, cut_ub, simple) -\n",
    "          err_vec4(xvals, mu, sigma - h_sig, cut_lb, cut_ub, simple)) / (2 * h_sig)).flatten()\n",
    "    \n",
    "    return Jac_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_err4 = Jac_err4(pts, mu_GMM1_4, sig_GMM1_4, 0.0, 450.0, False)\n",
    "print(d_err4)\n",
    "print(W_hat1_4)\n",
    "SigHat4 = (1 / N) * lin.inv(d_err4.T @ W_hat1_4 @ d_err4)\n",
    "print(SigHat4)\n",
    "print('Std. err. mu_hat=', np.sqrt(SigHat4[0, 0]))\n",
    "print('Std. err. sigma_hat=', np.sqrt(SigHat4[1, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how much tighter the standard errors are here with these four moments than they were in the econometric models of Sections 5.1.1 and 5.1.2 with two moments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1.4. Four moments, two-step optimal weighting matrix\n",
    "Let's see how much things change if we use the two-step estimator for the optimal weighting matrix $W$ instead of the identity matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_Err_mat4(pts, mu, sigma, cut_lb, cut_ub, simple=False):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the R x N matrix of errors from each\n",
    "    observation for each moment. In this function, we have hard coded\n",
    "    R = 4.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N,) vector, test scores data\n",
    "    mu     = scalar, mean of the normally distributed random variable\n",
    "    sigma  = scalar > 0, standard deviation of the normally distributed\n",
    "             random variable\n",
    "    cut_lb = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar lower bound value of distribution. Values below\n",
    "             this value have zero probability\n",
    "    cut_ub = scalar or string, ='None' if no cutoff is given, otherwise\n",
    "             is scalar upper bound value of distribution. Values above\n",
    "             this value have zero probability\n",
    "    simple = boolean, =True if errors are simple difference, =False if\n",
    "             errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        model_moments()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    R          = 2, hard coded number of moments\n",
    "    N          = integer >= R, number of data observations\n",
    "    Err_mat    = (R, N) matrix, error by moment and observation data\n",
    "    mean_model = scalar, mean value from model\n",
    "    var_model  = scalar > 0, variance from model\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: Err_mat\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    R = 4\n",
    "    N = len(pts)\n",
    "    Err_mat = np.zeros((R, N))\n",
    "    pct_1_mod, pct_2_mod, pct_3_mod, pct_4_mod = \\\n",
    "        model_moments4(mu, sigma, cut_lb, cut_ub)\n",
    "    if simple:\n",
    "        pts_in_grp1 = pts < 220\n",
    "        Err_mat[0, :] = pts_in_grp1 - pct_1_mod\n",
    "        pts_in_grp2 = (pts >= 220) & (pts < 320)\n",
    "        Err_mat[1, :] = pts_in_grp2 - pct_2_mod\n",
    "        pts_in_grp3 = (pts >= 320) & (pts < 430)\n",
    "        Err_mat[2, :] = pts_in_grp3 - pct_3_mod\n",
    "        pts_in_grp4 = pts >= 430\n",
    "        Err_mat[3, :] = pts_in_grp4 - pct_4_mod\n",
    "    else:\n",
    "        pts_in_grp1 = pts < 220\n",
    "        Err_mat[0, :] = (pts_in_grp1 - pct_1_mod) / pct_1_mod\n",
    "        pts_in_grp2 = (pts >= 220) & (pts < 320)\n",
    "        Err_mat[1, :] = (pts_in_grp2 - pct_2_mod) / pct_2_mod\n",
    "        pts_in_grp3 = (pts >= 320) & (pts < 430)\n",
    "        Err_mat[2, :] = (pts_in_grp3 - pct_3_mod) / pct_3_mod\n",
    "        pts_in_grp4 = pts >= 430\n",
    "        Err_mat[3, :] = (pts_in_grp4 - pct_4_mod) / pct_4_mod\n",
    "    \n",
    "    return Err_mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Err_mat4 = get_Err_mat4(pts, mu_GMM1_4, sig_GMM1_4, 0.0, 450.0, True)\n",
    "VCV2_4 = (1 / pts.shape[0]) * (Err_mat4 @ Err_mat4.T)\n",
    "print(VCV2_4)\n",
    "# We use the pseudo-inverse command here because the VCV matrix is poorly conditioned\n",
    "W_hat2_4 = lin.pinv(VCV2_4)\n",
    "print(W_hat2_4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this takes a little time because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_init = 600\n",
    "sig_init = 196\n",
    "params_init = np.array([mu_GMM1_4, sig_GMM1_4])\n",
    "gmm_args = (pts, 0.0, 450.0, W_hat2_4)\n",
    "results2_4 = opt.minimize(criterion4, params_init, args=(gmm_args),\n",
    "                          method='L-BFGS-B', bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_GMM2_4, sig_GMM2_4 = results2_4.x\n",
    "print('mu_GMM2_4=', mu_GMM2_4, ' sig_GMM2_4=', sig_GMM2_4)\n",
    "results2_4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the two-step estimator of the optimal weighting matrix creates a small change in the estimated $\\mu$ and $\\sigma$ paramters to $(\\mu=364,\\sigma=111)$ from $(\\mu=362,\\sigma=92)$ in the identity matrix estimation. The criterion function for different values of $\\mu$ and $\\sigma$ here has a clear minimum in a certain area. But it also has some really interesting nonlinearities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this will take a few minutes because the intgr.quad() commands\n",
    "# are a little slow\n",
    "# mu_vals = np.linspace(355, 370, 50)\n",
    "# sig_vals = np.linspace(85, 120, 50)\n",
    "mu_vals = np.linspace(60, 700, 50)\n",
    "sig_vals = np.linspace(20, 200, 50)\n",
    "crit_vals4 = np.zeros((50, 50))\n",
    "for mu_ind in range(50):\n",
    "    for sig_ind in range(50):\n",
    "        crit_vals4[mu_ind, sig_ind] = \\\n",
    "            criterion4(np.array([mu_vals[mu_ind], sig_vals[sig_ind]]),\n",
    "                       pts, 0.0, 450.0, W_hat2_4)\n",
    "\n",
    "mu_mesh, sig_mesh = np.meshgrid(mu_vals, sig_vals)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot_surface(sig_mesh, mu_mesh, crit_vals4, rstride=8,\n",
    "                cstride=1, cmap=cmap1)\n",
    "ax.set_title('Criterion function for values of mu and sigma')\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$\\mu$')\n",
    "ax.set_zlabel(r'Crit. func.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the estimator of the variance-covariance matrix $\\hat{\\Sigma}$ of the GMM parameter estimator by computing the Jacobian of the error vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_err4_2 = Jac_err4(pts, mu_GMM2_4, sig_GMM2_4, 0.0, 450.0, False)\n",
    "print(d_err4_2)\n",
    "print(W_hat2_4)\n",
    "SigHat4_2 = (1 / N) * lin.inv(d_err4_2.T @ W_hat2_4 @ d_err4_2)\n",
    "print(SigHat4_2)\n",
    "print('Std. err. mu_hat=', np.sqrt(SigHat4_2[0, 0]))\n",
    "print('Std. err. sigma_hat=', np.sqrt(SigHat4_2[1, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how much smaller (more efficient) the estimated standard errors are in this case with the two-step optimal weighting matrix $\\hat{W}_{2step}$: $(1.65, 1.49)$ versus $(3.78, 3.24)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. Unconditional and conditional expectations, instruments, and moments\n",
    "Most standard treatments of the generalized method of moments estimator in econometrics textbooks start with this principle and this selection of moments. However, this notebook follows the progression of starting with the most general treatment of GMM and then covering these special cases.\n",
    "\n",
    "In stochastic models, the assumed data generating process might have one or more characterizing equations that involve an unconditional expectation. The unconditional expectation is a strong assumption with many implications on conditional expectations that can create moments for identifying parameters using GMM. In econometric models, these unconditional expectations often show up as an assumption on the error term of one or more of the equations. Note that this is a minimal assumption and does not require knowledge of the distribution of the error term.\n",
    "\n",
    "$$ y_i = \\beta_0 + \\beta_1 x_{1,i} + \\beta_2 x_{2,i} + \\varepsilon_i \\quad\\text{where}\\quad E\\left[\\varepsilon_i\\right] = 0 $$\n",
    "\n",
    "In a macroeconomic model like the Brock and Mirman (1972) model, unconditional expectations show up in two places. The first is in the Euler equation, which is the first equation below. And the second is on the error term in the law of motion for the productivity shock, which is the last equation below.\n",
    "\n",
    "$$ \\quad\\left(c_t\\right)^{-1} = \\beta E\\left[r_{t+1}\\left(c_{t+1}\\right)^{-1}\\right] $$\n",
    "$$ c_t + k_{t+1} = r_{t+1}k_t + w_t $$\n",
    "$$ \\quad w_t = (1 - \\alpha)e^{z_t}k_{t}^\\alpha $$\n",
    "$$ r_t = \\alpha e^{z_t}k_{t}^{\\alpha-1} $$\n",
    "$$ \\qquad\\qquad\\qquad\\qquad\\qquad\\qquad z_{t} = \\rho z_{t-1} + (1 - \\rho)\\mu + \\varepsilon_t \\quad\\text{where}\\quad E[\\varepsilon_t]=0 $$\n",
    "\n",
    "It is valuable to note first that these unconditional expectations imply minimal restrictions on the stochastic distributions in the model. They only imply a restriction on the first moments of those particular parts of the distributions. Furthermore, because they are unconditional distributions (which is a strong assumption), they also imply restrictions on conditional distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\mathcal{I}$ be the set of variables that are in the information set of the model at the time the expectations operator in the model is formed. Let $w\\in\\mathcal{I}$ be the typical element (variable) in the information set. In a cross sectional econometric model, the variables in the information set are $w\\in\\mathcal{I}$ that could possibly be related to the dependent variable $y$ and were determined at the time the expectation was formed. In dynamic models or time series models, variables in the information set include any variables that were determined on or before the period in which the expectation was formed.\n",
    "\n",
    "The following sequence shows how an unconditional expectation can lead to moments that can identify parameters.\n",
    "\n",
    "$$ E[x] = 0 \\Rightarrow E[x|\\mathcal{I}] = 0 \\Rightarrow Cov[x,w] = 0 \\Rightarrow E[xw] = 0 $$\n",
    "\n",
    "The first equation states that the unconditional expectation of $x$ is zero. This implies that the conditional expectation of $x$ given anything else in the information set is also zero. This, in turn, implies that the covariance of $x$ and any element $w$ of the information set is zero so that the expectation of $x$ times $w$ is zero. It is this last equation that generates many of the moments used to identify parameters in GMM. Any variable in the instrument set $w\\in\\mathcal{I}$ can generate a moment condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2.1. Econometric (linear regression) GMM example using unconditional expectation for moments\n",
    "In the econometric model example above, there are three parameters to be estimated $(\\beta_0, \\beta_1, \\beta_2)$. We need at least three independent moments to identify these parameters. The first moment we can use is the assumption that the unconditional expectation of the error terms is zero. This has an empirical counterpart that the sample mean of the error terms should be zero.\n",
    "\n",
    "$$ \\text{(Theoretical moment 1)}\\quad E[\\varepsilon_i] = 0 \\quad\\Rightarrow\\quad E\\left[y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right]=0 $$\n",
    "$$ \\text{(Empirical moment 1)}\\qquad\\qquad\\qquad\\:\\:\\: \\frac{1}{N}\\sum_{i=1}^N\\left[y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right]=0 $$\n",
    "\n",
    "The next two moments are formed from the implications that the unconditional expectation on the errors implies that the expected value of the product of the error and any other variable in the information set is also zero.\n",
    "\n",
    "$$ \\text{(Theoretical moment 2)}\\quad E[\\varepsilon_i x_{1,i}] = 0 \\quad\\Rightarrow\\quad E\\Bigl[\\left(y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right)x_{1,i}\\Bigr]=0 $$\n",
    "$$ \\text{(Empirical moment 2)}\\qquad\\qquad\\qquad\\qquad \\frac{1}{N}\\sum_{i=1}^N\\Bigl[\\left(y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right)x_{1,i}\\Bigr]=0 $$\n",
    "\n",
    "$$ \\text{(Theoretical moment 3)}\\quad E[\\varepsilon_i x_{2,i}] = 0 \\quad\\Rightarrow\\quad E\\Bigl[\\left(y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right)x_{2,i}\\Bigr]=0 $$\n",
    "$$ \\text{(Empirical moment 3)}\\qquad\\qquad\\qquad\\qquad \\frac{1}{N}\\sum_{i=1}^N\\Bigl[\\left(y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i}\\right)x_{2,i}\\Bigr]=0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2.2. Brock and Mirman (1972) estimation by GMM using unconditional expectations for moments\n",
    "Give Brock and Mirman (1972) example. This is currently a problem set.\n",
    "\n",
    "Cite Hansen and Singleton (1982) example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Identification\n",
    "An issue that we saw in the examples from Section 5 is that there is some science as well as some art in choosing moments to identify the parameters in a GMM estimation.\n",
    "\n",
    "* The $\\mu$ and $\\sigma$ parameters were identified more precisely when using the two-step estimator of the optimal weighting matrix instead of the identity matrix.\n",
    "* The overidentified four-moment model of total scores produced much smaller standard errors for both $\\mu$ and $\\sigma$ than did the two-moment model.\n",
    "\n",
    "Suppose the parameter vector $\\theta$ has $K$ elements, or rather, $K$ parameters to be estimated. In order to estimate $\\theta$ by GMM, you must have at least as many moments as parameters to estimate $R\\geq K$. If you have exactly as many moments as parameters to be estimated $R=K$, the model is said to be *exactly identified*. If you have more moments than parameters to be estimated $R>K$, the model is said to be *overidentified*. If you have fewer moments than parameters to be estimated $R<K$, the model is said to be *underidentified*. There are good reasons to overidentify $R>K$ the model in GMM estimation as we saw in the previous example. The main reason is that not all moments are orthogonal. That is, some moments convey roughly the same information about the data and, therefore, do not separately identify any extra parameters. So a good GMM model often is overidentified $R>K$.\n",
    "\n",
    "One last point about GMM regards moment selection and verification of results. The real world has an infinite supply of potential moments that describe some part of the data. Choosing moments to estimate parameters by GMM requires understanding of the model, intuition about its connections to the real world, and artistry. A good GMM estimation will include moments that have some relation to or story about their connection to particular parameters of the model to be estimated. In addition, a good verification of a GMM estimation is to take some moment from the data that was not used in the estimation and see how well the corresponding moment from the estimated model matches that *outside moment*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Linear regression by GMM and relation to OLS\n",
    "The most common method of estimating the parameters of a linear regression is using the ordinary least squares (OLS) estimator. This estimator is just special type of generalized method of moments (GMM) estimator. A simple regression specification in which the dependent variable $y_i$ is a linear function of two independent variables $x_{1,i}$ and $x_{2,i}$ is the following:\n",
    "\n",
    "$$ y_i = \\beta_0 + \\beta_1 x_{1,i} + \\beta_2 x_{2,i} + \\varepsilon_i \\quad\\text{where}\\quad \\varepsilon_i\\sim N\\left(0,\\sigma^2\\right) $$\n",
    "\n",
    "Note that we can solve for the parameters $(\\beta_0,\\beta_1,\\beta_2)$ in a number of ways. And we can do it without making any assumptions about the distribution of the error terms $\\varepsilon_i$.\n",
    "\n",
    "One way we might choose the parameters is to choose $(\\beta_0,\\beta_1,\\beta_2)$ to minimize the distance between the $N$ observations of $y_i$ and the $N$ predicted values for $y_i$ given by $\\beta_0 + \\beta_1 x_{1,i} + \\beta_2 x_{2,i}$. You can think of the $N$ observations of $y_i$ as $N$ data moments. And you can think of the $N$ observations of $\\beta_0 + \\beta_1 x_{1,i} + \\beta_2 x_{2,i}$ as $N$ model moments. The least squares estimator minimizes the sum of squared errors, which is the sum of squared deviations between the $N$ values of $y_i$ and  $\\beta_0 + \\beta_1 x_{1,i} + \\beta_2 x_{2,i}$.\n",
    "\n",
    "$$ \\varepsilon_i = y_i - \\beta_0 - \\beta_1 x_{1,i} - \\beta_2 x_{2,i} $$\n",
    "\n",
    "$$ \\hat{\\theta}_{OLS} = \\theta:\\quad \\min_{\\theta} \\varepsilon^T\\, I \\, \\varepsilon $$\n",
    "\n",
    "Let the $N\\times 1$ vector of $y_i$'s be $Y$. Let the $N\\times 3$ vector of data $(1, x_{1,i}, x_{2,i})$ be $X$. And let the vector of three parameters $(\\beta_0, \\beta_1, \\beta_2)$ be $\\beta$. It can be shown that the OLS estimator for the vector of parameters $\\beta$ is the following.\n",
    "\n",
    "$$ \\hat{\\beta}_{OLS} = (X^T X)^{-1}(X^T Y) $$\n",
    "\n",
    "But you could also just estimate the coefficients using the criterion function in the GMM statement of the problem above. This method is called nonlinear least squares or generalized least squares. Many applications of regression use a weighting matrix in the criterion function that adjusts for issues like heteroskedasticity and autocorrelation. \n",
    "\n",
    "Lastly, many applications use a different distance metric than the weighted sum of squared errors for the difference in moments. Sum of squared errors puts a large penalty on big differences. Sometimes you might want to maximize the sum of absolute errors, which is sometimes called median regression. You could also minimize the maximum absolute difference in the errors, which is even more extreme than the sum of squared errors on penalizing large differences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. References\n",
    "* Adda, Jerome and Russell Cooper, *Dynamic Economics: Quantitative Methods and Applications*, MIT Press (2003).\n",
    "* Brock, William A. and Leonard J. Mirman, \"Optimal Economic Growth and Uncertainty: The Discounted Case,\" *Journal of Economic Theory*, 4:3, pp. 479-513 (June 1972).\n",
    "* Davidson, Russell and James G. MacKinnon, *Econometric Theory and Methods*, Oxford University Press (2004).\n",
    "* Fuhrer, Jeffrey C., George R. Moore, and Scott D. Schuh, \"Estimating the Linear-quadratic Inventory Model: Maximum Likelihood versus Generalized Method of Moments,\" *Journal of Monetary Economics*, 35:1, pp. 115-157 (February 1995).\n",
    "* Hansen, Lars Peter, \"Large Sample Properties of Generalized Method of Moments Estimators,\" *Econometrica*, 50:4, pp. 1029-1054 (July 1982).\n",
    "* Hansen, Lars Peter and Kenneth J. Singleton, \"Generalized Instrumental Variables Estimation of Nonlinear Rational Expectations Models\", *Econometrica*, 50:5, pp. 1269-1286 (September 1982).\n",
    "* Newey, Whitney K. and Kenneth D. West, \"A Simple, Positive, Semi-definite, Heteroskedasticy and Autocorrelation Consistent Covariance Matrix,\" *Econometrica*, 55:3, pp. 703-708 (May 1987)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}