{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Foundations of Computational Economics #34\n",
"\n",
"by Fedor Iskhakov, ANU\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Numerical integration, quadrature\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"\n",
"[https://youtu.be/cc14S679x2M](https://youtu.be/cc14S679x2M)\n",
"\n",
"Description: Gaussian quadrature. Monte Carlo integration."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Integration in economics\n",
"\n",
"- Expected (discounted) utility \n",
"- Expected (discounted) profits \n",
"- Bayesian posterior \n",
"- Likelihood function with unobservables \n",
"- Stochastic elements in (dynamic) economic models \n",
"\n",
"\n",
"*Most integrals can not be evaluated analytically*"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Two main approaches: Monte Carlo and quadrature\n",
"\n",
"1. Based on simulations – **Monte Carlo integration** \n",
"1. Based on the fixed points and weights – **quadrature integration** "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Newton-Cotes formulas\n",
"\n",
"Goal: definite integral $ \\int_a^b f(x) dx $\n",
"\n",
"Idea: Approximate the function with low order polynomial, then integrate\n",
"approximation\n",
"\n",
"1. First order >> Step function approximation \n",
" - Constant, level at midpoint of $ [a,b] $ \n",
"1. Second order >> Linear approximation \n",
" - Trapezoid rule \n",
"1. Third order >> Quadratic approximation \n",
" - Simpson rule "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Composite Newton-Cotes\n",
"\n",
"Preform Newton-Cotes on a grid separately on each interval\n",
"\n",
"- Equally spaced points \n",
"- Newton-Cotes on each sub-interval \n",
"\n",
"\n",
"*Note that the points are placed exogenously*"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Gaussian quadrature\n",
"\n",
"General formula\n",
"\n",
"$$\n",
"\\int_a^b f(x) dx = \\sum_{i=1}^n \\omega_i f(x_i)\n",
"$$\n",
"\n",
"- $ x_i \\in [a,b] $ quadrature nodes \n",
"- $ \\omega_i $ quadrature weights \n",
"\n",
"\n",
"*Note that the points are placed endogenously*"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Quadrature accuracy\n",
"\n",
"Suppose that $ \\{\\phi_k(x)\\}_{k=1,2,\\dots} $ is family of\n",
"polynomials of degree $ k $ *orthogonal* with respect to the\n",
"weighting function $ w(x) $\n",
"\n",
"- let $ q_k $ denote the leading coefficients so that $ \\phi_k(x)=q_k x^k + \\dots $ \n",
"- let $ x_i $, $ i=1,\\dots,n $ be $ n $ roots of $ \\phi_n(x) $ \n",
"- let $ \\omega_i = - \\frac{q_{n+1}/q_n}{\\phi'_n(x_i)\\phi_{n+1}(x_i)}>0 $ \n",
"\n",
"\n",
"Then\n",
"\n",
"- $ a"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### General Monte Carlo integration algorithm\n",
"\n",
"1. sample $ N $ points $ x_1,\\cdots,x_N $ from distribution $ p(x) $ of $ \\tilde{X} $ on $ \\Omega $ \n",
"1. approximate the expectation $ E \\left[ \\frac{f(\\tilde{X})}{p(\\tilde{X})} \\right] $ by the sample average \n",
"\n",
"\n",
"$$\n",
"I_f = \\int_{\\Omega} f(x)\\,dx \\approx \\frac{1}{N} \\sum_{i=1}^{N} \\frac{f(x_i)}{p(x_i)} = Q_f(N)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### General Monte Carlo integration algorithm (simple naive approach)\n",
"\n",
"1. sample $ N $ points $ x_1,\\cdots,x_N $ uniformly on $ \\Omega $ \n",
"1. approximate the expectation $ E \\left[ V f(\\tilde{X})\\right] $ by the sample average \n",
"\n",
"\n",
"$$\n",
"I_f = \\int_{\\Omega} f(x)\\,dx \\approx \\frac{V}{N} \\sum_{i=1}^{N} f(x_i), \\; V =\\int_{\\Omega} \\,dx\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams['figure.figsize'] = [12, 8]\n",
"\n",
"def mc_int_cube(f,ndims=1,N=1000):\n",
" '''Computes the integral of function f on a hypercube of dimension ndims\n",
" using Monte Carlo integration with N uniformly distributed points\n",
" Assume function f uses axis=0 for inputs, and can be vectorized in other axis\n",
" Return: value and standard error\n",
" '''\n",
" # generate uniform numbers on the hypercube\n",
" x = np.random.random(ndims*N).reshape(ndims,N) # uniform random numbers in a matrix\n",
" y = f(x) # function value\n",
" Q = y.mean() # sample average\n",
" seQ = y.std()/np.sqrt(N) # standard error of sample average\n",
" return Q,seQ"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of Monte Carlo samples : 1000\n",
"Estimate (pi_hat) : 3.088\n",
"Standard error (pi_hat) : 0.0530684087\n",
"Approximation error (pi_hat-pi) : -0.0535926536\n"
]
}
],
"source": [
"# pi example from video 33 as two-dim integral\n",
"# Approximate pi using 2-d Monte Carlo integration\n",
"\n",
"N=1000 # Number of Monte carlo Samples\n",
"g = lambda x: (x[0,:]**2 + x[1,:]**2)<1 # indicator function to inegrate\n",
"\n",
"q,se = mc_int_cube(g,ndims=2,N=N)\n",
"pi_hat = 4*q\n",
"se_pi_hat = 4*se\n",
"\n",
"print('Number of Monte Carlo samples : ', N);\n",
"print('Estimate (pi_hat) : ', pi_hat.round(10));\n",
"print('Standard error (pi_hat) : ', se_pi_hat.round(10));\n",
"print('Approximation error (pi_hat-pi) : ', (pi_hat-np.pi).round(10))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Properties of Monte Carlo integral\n",
"\n",
"**Consistency**: Law of large numbers ensures that the sample average converge to the mean\n",
"\n",
"$$\n",
"\\lim _{{N\\to \\infty }}Q_f(N)\n",
"=\\lim _{{N\\to \\infty }}{\\frac{1}{N}}\\sum _{{i=1}}^{N}\\frac{f(x_i)}{p(x_i)}\n",
"=E\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right] = \\int_{\\Omega} f(x)\\,dx = I_f\n",
"$$\n",
"\n",
"**Assymptotic Normality**: By the central limit theorem we have\n",
"\n",
"$$\n",
"\\sqrt{N}\\left(Q_f(N)-I_f \\right)\\ \\xrightarrow {d} \\ N\\left(0,\\sigma ^{2}\\right), \\; \\sigma^2= \\operatorname{Var}\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right]\n",
"$$\n",
"\n",
"The standard error of $ Q_f(N) $ is then given by $ \\sigma_{Q_f(N)}=\\sigma \\big/ \\sqrt{N} $"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Standard error of Monte Carlo integral\n",
"\n",
"Given our estimate $ Q_f(N) $ of $ I_f $, we can obtain an unbiased estimate of $ \\sigma^2= \\operatorname{Var}\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right] $\n",
"\n",
"$$\n",
"{\\hat{\\sigma}}^2_N=\\frac{1}{N-1}\\sum _{i=1}^N \\left(\\frac{f(x_i)}{p(x_i)}-Q_f(N)\\right)^2\n",
"$$\n",
"\n",
"and the estimate of the standard error of $ Q_f(N) $\n",
"\n",
"$$\n",
"{\\hat{\\sigma}}_{Q_f(N)}={\\hat{\\sigma}}_N \\big/ \\sqrt{N}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Convergence of Monte Carlo integrals\n",
"\n",
"The standard error of $ Q_f(N) $:\n",
"\n",
"- is given by $ \\sigma_{Q_f(N)}=\\sigma \\big/ \\sqrt{N} $ \n",
"- can be estimated by $ {\\hat{\\sigma}}_{Q_f(N)}={\\hat{\\sigma}}_N \\big/ \\sqrt{N} $ \n",
"\n",
"\n",
"**Decreases with the standard parametric rate** $ \\sqrt{N} $\n",
"\n",
"- doubling of precision requires 4 time as many random points \n",
"- but does not depend on the dimensionality of the integral, $ \\Omega $ can be high dimensional "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True value (pi) : 3.1415926536\n",
"Average estimate across all runs : 3.13976\n",
"Mean bias : -0.0018326536\n",
"Average std err across all runs : 0.0519299844\n",
"Std dev of bias : 0.0533409261\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHiCAYAAADh4aRaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7jsdV0v8PcnNmgpCMiGEIWlhSX1iNaOLM1jYabhhec8VpolGR2yzDrHOrm9lNnNbZ3T/aKcLNG8YCaPHLfHJJLMShQKL4SG4U4JFARUUFPB7/nj99s6Ltbaa2avmbUW+/t6Pc88a+Z3/cx3frPmPd/5zm+qtRYAAOjBV2x2AQAAsFGEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCL6xTVb2oqn5hTts6vqpuqaqDxtsXVdWPzWPb4/b+X1WdMa/tzbDfX62qj1XVRzZ63+xbVT2sqq7e7DoWaZ7P0Sn29R1V9f6N2Nc0avBnVXVTVb1jP9bflP8ZsEjCL+xDVe2pqs9U1c1V9fGq+oeqempVffG501p7amvtV6bc1sP3tUxr7UOttbu21m6bQ+2/VFV/vmz7j2qtnbPebc9Yx72S/GySk1prX73C/EOq6rVj+7Sqetiy+VVVL6yqG8bLb1RVTcx/QFVdWlWfHv8+YNp1l+3nYeP+X7ds+snj9IvW1xJJVS2N29q2jm0cMj62V1bVp8Z2+9OqWlpvfQeqaZ+j+2N8PL92Yl9/11r7ukXsaz89JMl3J7lna+2UWVfejP8ZsGjCL6ztMa21Q5OckGRXkmcmecm8d7KeQLTFnZDkhtbadftY5m1JfijJSj3DZyU5PcnJSe6f5NFJfjwZgmCS1yf58yRHJDknyevH6ftcdxXXJ/n2qrr7xLQzkvzrPtbZaK9N8tgkP5jkbhnu26VJTp11Q3eUY+6OUucWdUKSPa21T212IbBltNZcXFxWuSTZk+Thy6adkuQLSb5xvP3SJL86Xj8qyRuSfDzJjUn+LsObzJeP63wmyS1Jfj7JUpKW5MwkH0ry1olp28btXZTkBUnekeQTGYLekeO8hyW5eqV6kzwyyeeSfH7c37smtvdj4/WvSPLcJP+e5LokL0tyt3He3jrOGGv7WJLn7KOd7jauf/24veeO23/4eJ+/MNbx0jXa++okD1s27R+SnDVx+8wkbx+vPyLJfySpifkfSvLItdZdYd8PG/f/oiRPG6cdNE77xSQXTSz77UneOT4m70zy7RPzLkryK0n+PsnNSd6c5KiJ2trYFrck+bZx+o8muSLJTUn+KskJq9S4tz3vtY82fMq4rZuTXJXkx1e4j8/M8Ebj5cuPoyT3G+/Dx5NcnuSx69zXs8fjZ0+SJ03Mf+nY1heM6//t5P0e2+lpSa5M8sF9tXuSI8d9PWa8fdckH0jy5BWeo3vr+vkMx/21Gd4gfW+GNzk3Jnn2suf7P47tcW2SP0hyyDjvrWOdnxofzx+YpT3Huv4wye6xDS5O8jXjvEry22ONn0jy7oz/c1Z4HO6R5Pyx9g8k+W8Tx/t/JrltrO/5K6z7IxmO1d8f9/O+JKcuO55/bNH/a11cNvKi5xdm1Fp7R4YXz+9YYfbPjvO2Jzkmwwt/a639cIbg85g2DGv4jYl1/kuGF8jvWWWXT84Qju6R5NYkvzdFjW9K8utJzh33d/IKi/3IePnOJPfJEBj+YNkyD0nydRl6FX+xqu63yi5/P0MAvs94f56c5Cmttb9O8qgk14x1/Mhata/gG5K8a+L2u8Zpe+e9u7U2+Tvt7142f7V1V/Oysf5keEwuT3LN3plVdWSGsPJ7Se6e5LeS7F7WW/yDGYLh0UkOSfJz4/SHjn8PH9vjH6vq9AzHyX/NcNz8XZJXrVLbw5O8o7X24X3Uf12GHu7Dxhp+u6q+aWL+V2cIiydk6Bn/oqo6OMn/zRDYj07y9CSvqKrVPsafZl9HJTkuwxups5dt60kZ3igcleSyJK9Ytv3Tk3xrkpP21e6ttRszPEf+T1UdnSE0XtZae9kqdX91kjuPdf1ikv+T4ZOHb87wvP7FqrrPuOxtSf7HWOO3ZXgu/GSStNb2Pp4nj4/nuZM7mbI9n5jk+Rk+ufhAkl8bpz8iw/Fy3ySHZwjWN6xyf16V4f/OPZI8PsmvV9WprbWXJHlqkn8c63veKut/a4Y3L0cleV6S143tDQck4Rf2zzUZAsRyn09ybIYerM+3YfxfW2G5Sb/UWvtUa+0zq8x/eWvtvW342PIXknz/3i/ErdOTkvxWa+2q1totSZ6V5AnLPmJ+fmvtM621d2UIjrcL0WMtP5DkWa21m1tre5L87yQ/PIcakyGUf2Li9ieS3HUcu7t83t75h06x7opaa/+Q5MgxoDw5QxiedFqSK1trL2+t3dpae1WG3rLHTCzzZ621fx0f09ckeUBW9+NJXtBau6K1dmuGNy0PqKoTVlj27hl6H1fVWtvdWvu3NvjbDMFr8o3aF5I8r7X22RWOuQdlaLNdrbXPtdb+JsMnGU/cz30lyS+M+/rbDOH1+yfm7W6tvbW19tkkz0nybeMY8b1e0Fq7caxzn+3eWntzkr9IcuG47L6Gt3w+ya+11j6f5NUZQt/vjsfv5Rne8Nx/3O6lrbW3j/vck+TFGd7gTWOa9nxda+0d42P/inzpWPl8huP46zN8snFFa+12j/3YXg9J8szW2n+21i5L8ieZ7fl3XZLfGf9nnZvk/RnaEA5Iwi/sn+MyfMS43G9m6L15c1VdVVU7p9jWvnrxls//9yQHZ3ixXq97jNub3Pa2DD3We02Owf10hhfy5Y7K0Lu5fFvHzaHGZPi49rCJ24cluWV8U7F83t75N0+x7r68PMlPZegVP2/ZvOXtltz+/k7TbnudkOR3xy9U7h0uU1m5/W7I8OZqVVX1qKp6e1XdOG7ve/Plx8v1rbX/XGX1eyT5cGvtCxPTVn0sp9jXTe3Lx5r++7iPvb54bI9vwG5cbX6ma/ezk3xjhjcfq/WSJsMY9L1fKt37BuCjE/M/k/Exq6r7VtUbquojVfXJDG9Opn3+TdOeKx4rY1D+gwzDIj5aVWdX1fJjfe8+bmyt3Twxbdbn338se04sf5zggCL8woyq6lsyvLC8bfm8sefoZ1tr98nQI/WMqtr7RaTVAtdaQWyyJ+z4DD1CH8swzvCrJuo6KMPH5tNu95oMwWty27fmy0PAND421rR8W/8x43ZWc3m+vMf55HHa3nn3X9aTe/9l81dbd19enuGj7Te21j69bN7ydkumv78rPSYfzjBW9vCJy1eOPdDL/XWSU6rqnittvKrulOQvk/yvJMe01g5P8sYMYXpfNex1TZJ7TZ7NJKvctyn3dURV3WXZtq6ZuP3FY7uq7prh05TJ+ZO17rPdx+P/xRl66n9i8gwM6/THGXqYT2ytHZZhiMqqnxwsM3V7rqS19nuttW/OMFTnvkn+5yr7OLKqDp2YNuvz77hlz6HljxMcUIRfmFJVHVZVj87wMemft9bes8Iyj66qrx1fSD6ZYbzg3h6mj2YYEzurH6qqk6rqq5L8cpLXjr1W/5rkzlV12ji28LlJ7jSx3keTLC174Z30qiT/o6ruPQaPvWOEb52luLGW1yT5tao6dPy4/hkZzsAwlaq6U1Xdebx5SFXdeeLF+GUZ3kQcV1X3yDCu+qXjvIsytO9Pj9v4qXH630yx7r7u0wczfLT9nBVmvzHJfavqB6tqW1X9QJKTMnycvZbrMww7mDwOXpTkWVX1DUlSVXerqu9bpa6/zvAFsfOq6pvH/R9aw+n3fjRDD/ydxv3cWlWPyjB2dFoXZ3hT9fNVdXANp517TIZjfrlp9/X88fRs35FhfPBfTMz73qp6yHh2jl9JcvE+xjOv1e7PHv/+aIZA/rI5DQ86NMNz+Zaq+vokP7Fs/r6e17O055epqm+pqm8dn9ufype+uPZlxvb6hyQvGJ8398/wRbfl46f35egMz6GDx2PvfhnaGw5Iwi+s7f9W1c0Zeuiek+GLNk9ZZdkTM/TO3ZLhG+J/1Fq7aJz3giTPHT/e/rlV1l/JyzMEto9k+JLOTydJa+0TGXon/yRDL8+nMnzpZa+9IeOGqvqnFbb7p+O235rkgxleXJ8+Q12Tnj7u/6oMPeKvHLc/rfdn+Kj5uAxnO/hMvtTL9+IMXxp6T5L3Zhg3+uIkaa19LsOXop6c4dv0P5rk9HH6PtddS2vtba212/V+jR+nPzpDkL4hw1kDHt1a+9gU2/x0hi80/f14HDyotXZekhcmefX4sfp7M3xJcDWPzxBMzs0whvm9SXYk+evxo++fzvBm5KYMX7w7f5r7O9b3uQynUXtUhh79P8pwxoT3rbDsNPv6yDjvmgxh7KnLtvXKDF+wujHDl82etI/aVm33qvrmDG+4njy+GXthhl7jaYYdreXnMty3mzN8Me7cZfN/Kck54+M5OZ55pvZcwWHj/m7KMAzhhgyhfiVPzHCGlmsyDNN5Xmvtgin2sdfFGf53fSzD8fn4NYaNwB1arT30DQBmM/Zy/nlrbbUhGi/NcEqw525kXXy5qvqRDKcye8hm1wIbRc8vAADdEH4BAOiGYQ8AAHRDzy8AAN0QfgEA6Ma2tReZn6OOOqotLS1t5C4BAOjMpZde+rHW2vaV5m1o+F1aWsoll1yykbsEAKAzVbX859C/yLAHAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbmzb7AIA2LelnbtnWn7PrtMWVAnAHZ+eXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6MZU4beqDq+q11bV+6rqiqr6tqo6sqouqKorx79HLLpYAABYj2l7fn83yZtaa1+f5OQkVyTZmeTC1tqJSS4cbwMAwJa1ZvitqsOSPDTJS5Kktfa51trHkzwuyTnjYuckOX1RRQIAwDxsm2KZ+yS5PsmfVdXJSS5N8jNJjmmtXZskrbVrq+rolVauqrOSnJUkxx9//FyKBliUpZ27Z1p+z67TFrp9AOZrmmEP25J8U5I/bq09MMmnMsMQh9ba2a21Ha21Hdu3b9/PMgEAYP2mCb9XJ7m6tXbxePu1GcLwR6vq2CQZ/163mBIBAGA+1gy/rbWPJPlwVX3dOOnUJP+S5PwkZ4zTzkjy+oVUCAAAczLNmN8keXqSV1TVIUmuSvKUDMH5NVV1ZpIPJfm+xZQIAADzMVX4ba1dlmTHCrNOnW85AACwOH7hDQCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDoxrbNLgCAA9vSzt0zLb9n12kLqgRAzy8AAB0RfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDoxrbNLgBgkZZ27t7sEgDYQvT8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohlOdAayDU6kB3LHo+QUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A2nOgM4wMx6+rU9u05bUCUAW4+eXwAAujFVz29V7Ulyc5LbktzaWttRVUcmOTfJUpI9Sb6/tXbTYsoEAID1m6Xn9ztbaw9ore0Yb+9McmFr7cQkF463AQBgy1rPsIfHJTlnvH5OktPXXw4AACzOtOG3JXlzVV1aVWeN045prV2bJOPfoxdRIAAAzMu0Z3t4cGvtmqo6OskFVfW+aXcwhuWzkuT444/fjxIB2EpmPZsEwFYyVc9va+2a8e91Sc5LckqSj1bVsUky/r1ulXXPbq3taK3t2L59+3yqBgCA/bBm+K2qu1TVoXuvJ3lEkvcmOT/JGeNiZyR5/aKKBACAeZhm2MMxSc6rqr3Lv7K19qaqemeS11TVmUk+lOT7FlcmAACs35rht7V2VZKTV5h+Q5JTF1EUAAAsgl94AwCgG8IvAADdmPZUZwBbgtNsAbAeen4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0I1tm10AAJtraefuzS4BYMPo+QUAoBvCLwAA3RB+AQDohvALAEA3hF8AALrhbA/ApnKmAQA2kp5fAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuTB1+q+qgqvrnqnrDePveVXVxVV1ZVedW1SGLKxMAANZvlp7fn0lyxcTtFyb57dbaiUluSnLmPAsDAIB5myr8VtU9k5yW5E/G25Xku5K8dlzknCSnL6JAAACYl2l7fn8nyc8n+cJ4++5JPt5au3W8fXWS4+ZcGwAAzNWa4beqHp3kutbapZOTV1i0rbL+WVV1SVVdcv311+9nmQAAsH7T9Pw+OMljq2pPkldnGO7wO0kOr6pt4zL3THLNSiu31s5ure1ore3Yvn37HEoGAID9s2b4ba09q7V2z9baUpInJPmb1tqTkrwlyePHxc5I8vqFVQkAAHOwnvP8PjPJM6rqAxnGAL9kPiUBAMBibFt7kS9prV2U5KLx+lVJTpl/SQAAsBh+4Q0AgG7M1PMLsC9LO3dvdgkAsE96fgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdMOpzgDYUjbilHl7dp228H0AW5OeXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAurFtswsANs7Szt0zLb9n12kLqgQ2l+cC9EvPLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBurBl+q+rOVfWOqnpXVV1eVc8fp9+7qi6uqiur6tyqOmTx5QIAwP6bpuf3s0m+q7V2cpIHJHlkVT0oyQuT/HZr7cQkNyU5c3FlAgDA+q0ZftvglvHmweOlJfmuJK8dp5+T5PSFVAgAAHMy1Zjfqjqoqi5Lcl2SC5L8W5KPt9ZuHRe5OslxiykRAADmY9s0C7XWbkvygKo6PMl5Se630mIrrVtVZyU5K0mOP/74/SwTWMnSzt2bXQIA3KHMdLaH1trHk1yU5EFJDq+qveH5nkmuWWWds1trO1prO7Zv376eWgEAYF2mOdvD9rHHN1X1lUkenuSKJG9J8vhxsTOSvH5RRQIAwDxMM+zh2CTnVNVBGcLya1prb6iqf0ny6qr61ST/nOQlC6wTAADWbc3w21p7d5IHrjD9qiSnLKIoAABYBL/wBgBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOjGNL/wBnRqaefuzS4BAOZKzy8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBubNvsAgDgQLO0c/dMy+/ZddqCKgGW0/MLAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbTnUGAGuY9dRlwNal5xcAgG4IvwAAdGPN8FtV96qqt1TVFVV1eVX9zDj9yKq6oKquHP8esfhyAQBg/03T83trkp9trd0vyYOSPK2qTkqyM8mFrbUTk1w43gYAgC1rzfDbWru2tfZP4/Wbk1yR5Lgkj0tyzrjYOUlOX1SRAAAwDzON+a2qpSQPTHJxkmNaa9cmQ0BOcvQq65xVVZdU1SXXX3/9+qoFAIB1mDr8VtVdk/xlkv/eWvvktOu11s5ure1ore3Yvn37/tQIAABzMVX4raqDMwTfV7TWXjdO/mhVHTvOPzbJdYspEQAA5mOasz1UkpckuaK19lsTs85PcsZ4/Ywkr59/eQAAMD/T/MLbg5P8cJL3VNVl47RnJ9mV5DVVdWaSDyX5vsWUCAAA87Fm+G2tvS1JrTL71PmWAwAAi+MX3gAA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbmzb7AJgsyzt3D3T8nt2nbagSgCAjaLnFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IazPcACOaMEsCj+v8D+0fMLAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG5s2+wCgC9Z2rl7s0sANoHnPmwcPb8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALrhVGcAwO3Mevq1PbtOW1AlMF96fgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdWDP8VtWfVtV1VfXeiWlHVtUFVXXl+PeIxZYJAADrN03P70uTPHLZtJ1JLmytnZjkwvE2AABsaWuG39baW5PcuGzy45KcM14/J8npc64LAADmbn/H/B7TWrs2Sca/R8+vJAAAWIyFf+Gtqs6qqkuq6pLrr79+0bsDAIBV7W/4/WhVHZsk49/rVluwtXZ2a21Ha23H9u3b93N3AACwfvsbfs9PcsZ4/Ywkr59POQAAsDjTnOrsVUn+McnXVdXVVXVmkl1Jvruqrkzy3eNtAADY0rattUBr7YmrzDp1zrUAAMBC+YU3AAC6IfwCANCNNYc9AIOlnbs3uwQAYJ30/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG9s2uwAAYPGWdu7e7BJgS9DzCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdMPZHtiyZv1m8p5dpy2oEgDmzf94NoueXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3XCqMw4Ys542BwDoj55fAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRj22YXAADc8S3t3L3ltr9n12kLqIQ7Oj2/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOjGts0uYCta2rl74fvYs+u0hW5/1vuw6HqSjWlXANhrK74WLtL+vM7e0e/z/tDzCwBAN4RfAAC6sa7wW1WPrKr3V9UHqmrnvIoCAIBF2O/wW1UHJfnDJI9KclKSJ1bVSfMqDAAA5m09Pb+nJPlAa+2q1trnkrw6yePmUxYAAMzfesLvcUk+PHH76nEaAABsSes51VmtMK3dbqGqs5KcNd68parev459bqSjknxsURuvFy5qy/tnDvUstL0OQNprNtprNtprdtpsNgdkey3wtXnLttdWyyOjebTXCavNWE/4vTrJvSZu3zPJNcsXaq2dneTsdexnU1TVJa21HZtdxx2F9pqN9pqN9pqN9pqdNpuN9pqN9prNottrPcMe3pnkxKq6d1UdkuQJSc6fT1kAADB/+93z21q7tap+KslfJTkoyZ+21i6fW2UAADBn6/p549baG5O8cU61bDV3uKEam0x7zUZ7zUZ7zUZ7zU6bzUZ7zUZ7zWah7VWt3e47agAAcEDy88YAAHSju/BbVfeqqrdU1RVVdXlV/cwKyxxRVedV1bur6h1V9Y0T8/ZU1Xuq6rKqumRjq994VXXnsQ3eNbbX81dY5k5Vde74M9cXV9XSxLxnjdPfX1Xfs5G1b4b1tFdVLVXVZ8Zj67KqetFG17/Rpmyvh1bVP1XVrVX1+GXzzqiqK8fLGRtX+eaYQ3vdNnF8HfBfUJ6yvZ5RVf8y/r+/sKpOmJjn+Lr9MvtqL8fX7Zd56kRmeNvkL+F6fZy+veb++tha6+qS5Ngk3zRePzTJvyY5adkyv5nkeeP1r09y4cS8PUmO2uz7sYHtVUnuOl4/OMnFSR60bJmfTPKi8foTkpw7Xj8pybuS3CnJvZP8W5KDNvs+beH2Wkry3s2+D1uwvZaS3D/Jy5I8fmL6kUmuGv8eMV4/YrPv01Ztr3HeLZt9H7Zge31nkq8ar//ExPPR8TVDezm+Vm2vwyauPzbJm8brXh9na6+5vj521/PbWru2tfZP4/Wbk1yR2/8y3UlJLhyXeV+Spao6ZkML3SLa4Jbx5sHjZflA8cclOTz33/IAAAQfSURBVGe8/tokp1ZVjdNf3Vr7bGvtg0k+kOFnsQ9Y62yv7kzTXq21Pa21dyf5wrLVvyfJBa21G1trNyW5IMkjF13zZlpne3VnyvZ6S2vt0+PNt2c4Z33i+Jq1vbozZXt9cuLmXSbme32crb3mqrvwO2n8uPmBGd59THpXkv86LnNKhl8J2fsEb0neXFWX1vDrdQe8qjqoqi5Lcl2GF4Pl7fXFn7purd2a5BNJ7p5OfwJ7He2VJPeuqn+uqr+tqu/YsKI30RTttRrH12ztlSR3rqpLqurtVXX6gkrcUmZsrzOT/L/xuuNrtvZKHF8rtldVPa2q/i3JbyT56XGy42u29krm+PrYbfitqrsm+csk/33ZO40k2ZXkiPEBenqSf05y6zjvwa21b0ryqCRPq6qHblTNm6W1dltr7QEZ3gCcUhNjoEer/dT1VD+BfaBZR3tdm+T41toDkzwjySur6rDFVrv5pmiv1Ti+ZmuvZDi+diT5wSS/U1Vfs5Ait5Bp26uqfijJjgzD3hLH16ztlTi+Vmyv1toftta+Jskzkzx3nOz4mq295vr62GX4raqDMwTfV7TWXrd8fmvtk621p4wP0JOTbE/ywXHeNePf65KclwP8Y4pJrbWPJ7kot//o74s/dV1V25LcLcmNmfInsA9Us7bX+PHXDeO6l2YYA3bfDSt4k+2jvVbj+JqtvSb/f101rvvARdS2Fe2rvarq4Umek+SxrbXPjpMdX7O1l+Nr7efjq5Ps7RF3fM3QXvN+fewu/I5jK1+S5IrW2m+tsszhNfxkc5L8WJK3ttY+WVV3qapDx2XukuQRSd67EXVvlqraXlWHj9e/MsnDk7xv2WLnJ9n7TejHJ/mb1lobpz+hhrMb3DvJiUnesTGVb471tNe47kHjuvfJ0F5XbUzlm2PK9lrNXyV5RA1nZzkiw/PxrxZT6dawnvYa2+lO4/Wjkjw4yb8sqtatYJr2qqoHJnlxhiB33cQsx9cM7eX4WrW9Tpy4eVqSK8frXh9naK95vz6u6xfe7qAenOSHk7xnHNaQJM9OcnyStNZelOR+SV5WVbdlePKeOS53TJLzxu8mbUvyytbamzaw9s1wbJJzxoPuK5K8prX2hqr65SSXtNbOz/Bm4uVV9YEMPb5PSJLW2uVV9ZoMbXhrkqe11m7blHuxcfa7vZI8NMkvV9WtSW5L8tTW2o0bfxc21JrtVVXfkuFTliOSPKaqnt9a+4bW2o1V9StJ3jlu65e11+rtleH/2our6gvjurtaawd0OMl0z8ffTHLXJH8x/m//UGvtsY6v2dorjq/V2uunxp7yzye5KWPHh9fH2dorc3599AtvAAB0o7thDwAA9Ev4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBu/H951TSwopBiYAAAAABJRU5ErkJggg==\n",
"text/plain": [
"