ipynb/pablo_winant-dl_notebook.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "# Deep learning solution method with all-in-one expectation operator\n\n- This companion notebook presents one of the deep learning methods for solving dynamic economic models introduced by Maliar, Maliar and Winant (2019) in the paper \"Will Artificial Intelligence Replace Computational Economists Any Time Soon?\", CEPR working paper DP 14024, (https://lmaliar.ws.gc.cuny.edu/files/2019/10/CEPR-DP14024.pdf).\n- In the paper, we offer a unified approach for casting three fundamental objects of economic dynamics -- lifetime reward, Bellman equation and Euler equation -- into objective functions of the deep learning framework. \n- In this notebook, we illustrate only one of the three approaches -- the Euler residual minimization. \n- We solve a cannonical consumption-saving problem with occasionally binding borrowing constraint and four exogenous stochastic shocks. \n- We parameterize the agent's decision function with a multilayer neural network, and we perform training using stochastic optimization, namely, in each iteration, we train the model on just one or few grid points that are randomly drawn from the state space (instead of conventional fixed grid with a potentially large number of grid points constructed by a discretization of the state-space). \n- Our objective function -- the sum of squared residuals in the Euler equation -- has two types of expectation operators, one is with respect to current state variables (which arises because grid points that are randomly drawn from the state space), and the other is with respect to future state variables (which arises because next-period shocks are randomly drawn from the given distributions).\n- We construct an all-in-one expectation method that merges the two expectation operators into one. Namely, we use two independent random draws for evaluating two terms of a squared residual -- this method eliminates the correlation between the terms and pulls the expectation operator out of the square. Our all-in-one expectation operator allows for efficient parallel calculations and reduces greately the cost of training deep neural networks.\n- In our 5-dimensional problem, we obtain a solution in a mintue or so. Importantly, our DL method is highly scalable and can solve very large models intractable with other methods. In particular, in the paper, we solve Krusell and Smith's (1998) model with 1,000 of heterogeneous agents by constructing decision function of 2,001 state variables."
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## TensorFlow and libraries\n\n[TensorFlow](https://www.tensorflow.org/) is a deep learning library. It can be used to train neural networks, but also to produce computational graphs (aka programs), which can be differentiated automatically, optimized and ran on very scalable architecture.\nVersion 2.0 introduces a new way to build graph, which allows for a more intuitive graph definition, essentially by writing numpy-like code. We can install it using:\n`pip install tensorflow==2.0.0`"
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "metadata": {},
            "outputs": [],
            "source": [
                "import tensorflow as tf"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "This notebook uses several libraries. Any missing library x can be installed using `pip install x`"
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {
                "lines_to_next_cell": 0
            },
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "from math import sqrt\n",
                "from matplotlib import pyplot as plt\n",
                "from tqdm import tqdm as tqdm         # tqdm is a nice library to visualize ongoing loops\n",
                "import datetime\n",
                "# followint lines are used for indicative typing\n",
                "from typing import Tuple\n",
                "class Vector: pass"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## The model\n\nWe consider the following consumption-saving problem: \n\n$$\n\\begin{gather*}\n\\underset{\\left\\{ c_{t},w_{t+1}\\right\\}_{t=0}^{\\infty }}{\\max }E_{0}\\left[\n\\sum_{t=0}^{\\infty }\\exp (\\delta_{t})\\beta ^{t}u\\left( {c_{t}}\\right)\\right]  \\\\\n\\text{s.t. }w_{t+1}=\\left( w_{t}-c_{t}\\right) \\overline{r}\\exp (r_{t+1})+\\exp\n(y_{t+1}), \\\\\nc_{t}\\leq w_{t},\n\\end{gather*}\n$$\n\n\nwhere $c_{t}$ is consumption; $w_{t}$ is the beginning-of-period\ncash-on-hand; $\\beta \\in \\left[ 0,1\\right)$ is a subjective discount\nfactor;  $\\overline{r}\\in \\left( 0,\\frac{1}{\\beta }\\right)$ is a (gross)\nconstant interest rate; and initial condition $\\left( z,w\\right)$ is given.\nThere is an occasionally binding inequality constraint: consumption $c_{t}$\ncannot exceed cash-on-hand $w_{t}$. There are four different exogenous\nstate variables, namely, shocks to the interest rate $r_{t}$, discount factor $\\delta_t$,\ntransitory component of income $q_{t}$ and permanent\ncomponent of income $p_{t}$. The total income is $y_{t}=\\exp(p_{t})\\exp(q_{t})$. All\nexogenous variables  follow AR(1) processes: \n\n\n$$\n\\begin{eqnarray*}\n{y_{t+1}} &=&\\rho_{y}{y_{t}}+\\sigma_{y}\\epsilon^{y}_{t+1}, \\\\\n{p_{t+1}} &=&\\rho_{p}{p_{t}}+\\sigma_{p}\\epsilon^{p}_{t+1}, \\\\\n{r_{t+1}} &=&\\rho_{r}{r_{t}}+\\sigma_{r}\\epsilon^{r}_{t+1}, \\\\\n{\\delta_{t+1}} &=&\\rho_{\\delta }{\\delta_{t}}+\\sigma_{\\delta }\\epsilon^{\\delta}_{t+1},\n\\end{eqnarray*}\n$$\n\n\nwhere $\\epsilon^i_t \\sim \\mathcal{N}\\left( 0,1\\right)$. We assume the Cobb-Douglas utility function $u\\left( {c_{t}}\\right) =\\frac{1}{1-\\gamma }\\left( c_{t}^{1-\\gamma }-1\\right)$."
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "metadata": {
                "lines_to_next_cell": 2
            },
            "outputs": [],
            "source": [
                "# Model parameters\n",
                "\n",
                "β = 0.9 \n",
                "γ = 2.0 \n",
                "# σ = 0.1  \n",
                "# ρ = 0.9\n",
                "σ_r = 0.001\n",
                "ρ_r = 0.2\n",
                "σ_p = 0.0001\n",
                "ρ_p = 0.999\n",
                "σ_q = 0.001\n",
                "ρ_q = 0.9\n",
                "σ_δ = 0.001\n",
                "ρ_δ = 0.2\n",
                "rbar = 1.04"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## Stochastic solution domain\n\nOur DL framework converts a dynamic stochastic model into a static optimization problem. We will solve the static problem on a random grid of points which are drawn from the following 5-dimentional domain:\n\n- $y, p, r, \\delta$ are drawn from their ergodic distribution (recall that for an AR(1) process $z$ with autocorrelation $\\rho$ and conditional standard deviation $\\sigma$, the ergodic distribution is normal with zero mean and standard deviation $\\sigma_z= \\frac{\\sigma}{\\sqrt{1-\\rho^2}}$).\n- $w$ is drawn from a uniform distribution within an interval: $w\\in[w_{\\min}, w_{\\max}]$."
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Standard deviations for ergodic distributions of exogenous state variables\n",
                "σ_e_r = σ_r/(1-ρ_r**2)**0.5\n",
                "σ_e_p = σ_p/(1-ρ_p**2)**0.5\n",
                "σ_e_q = σ_q/(1-ρ_q**2)**0.5\n",
                "σ_e_δ = σ_δ/(1-ρ_δ**2)**0.5\n",
                "\n",
                "# bounds for endogenous state variable\n",
                "wmin = 0.1\n",
                "wmax = 4.0"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "In the recursive form, the solution can be characterized by the Kuhn-Tucker\n(KT) conditions \n\n$$\n\\begin{equation*}\na\\geq 0,\\quad b\\geq 0\\quad and\\quad ab=0,\n\\end{equation*}\n$$\n\nwhere $a$ is the share of wealth that goes to savings and $b$ is the\nLagrange multiplier\n\n$$\n\\begin{eqnarray*}\na &\\equiv &w-c, \\\\\nb &\\equiv &u^{\\prime }(c)-\\beta \\overline{r}E_{\\epsilon }\\left[ u^{\\prime\n}\\left( c^{\\prime }\\right) \\exp \\left( \\delta ^{\\prime }-\\delta +r^{\\prime\n}\\right) \\right] .\n\\end{eqnarray*}\n$$\n\n(In the absence of borrowing constraint $b=0$, the KT conditions lead to the familiar\nEuler equation). \n\nInequality constraints are not directly compatible with the deep learning\nframework developed in the paper, so we reformulate the KT conditions as a set of equations that hold with equality. We use a smooth representation of the KT conditions, called\nthe Fischer-Burmeister (FB) function, which is differentiable \n\n$$\n\\begin{equation*}\n\\Psi_{FB}\\left( a,b\\right) =a+b-\\sqrt{a^{2}+b^{2}}=0.\n\\end{equation*}\n$$\n\nThe restriction $\\Psi_{FB}\\left( a,b\\right) =0$ is also equivalent to the\nKT conditions. \n\nFor numerical treatment, we rewrite the $\\Psi_{FB}$ function in the following\nunit-free form \n\n$$\n\\begin{equation*}\n\\Psi_{FB}\\left( 1-\\zeta ,1-h\\right) =(1-\\zeta)+(1-h)-\\sqrt{(1-\\zeta)^{2}+(1-h)^{2}}=0,\n\\end{equation*}\n$$\n\nwhere $\\zeta $ and $h$ are respectively the consumption share and normalized\nLagrange multiplier \n\n$$\n\\begin{eqnarray*}\n\\zeta  &\\equiv &\\frac{c}{w}, \\\\\nh &\\equiv &\\beta \\overline{r}E_{\\epsilon }\\left[ \\frac{u^{\\prime }\\left(\nc^{\\prime }\\right) }{u^{\\prime }(c)}\\exp \\left( \\delta ^{\\prime }-\\delta\n+r^{\\prime }\\right) \\right] .\n\\end{eqnarray*}\n$$\n\nIn particular, $\\zeta $ belongs to the interval $\\left[ 0,1\\right] $ which\nis a convenient domain for defining neural network. In turn, $h$ is always positive and is normalized \nto be around one: we will parameterize it with neural network in the way \nthat ensures that it is nonnegative."
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Here is the  Fischer-Burmeister (FB) in TensorFlow:\n",
                "min_FB = lambda a,b: a+b-tf.sqrt(a**2+b**2)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## Parameterizing decision functions with neural network\n\nThere are many different decision functions that we can approximate for characterizing the solution, including\nconsumption, next-period income, etc. We chose to approximate the two functions that\nwe defined earlier: the share of consumption, $\\zeta \\equiv \\frac{c}{w}$, and the\nnormalized Lagrange multiplier $h$.\nSince the model is stationary, \nwe look for a decision rule \n\n$$\n\\begin{equation*}\n\\left( \n\\begin{matrix}\n\\zeta  \\\\ \nh%\n\\end{matrix}%\n\\right) =\\varphi (s;\\theta ),\n\\end{equation*}\n$$\n\nwhere $s=(r, \\delta, q, p, w)$ is the 5-dimensional state space, and $\\varphi$ is a function to be determined.\n\nA common approach in computational economics is to approximate an unknown function $\\varphi$ using some flexible function family $\\varphi(...;\\theta)$ parameterized by a vector of coefficients $\\theta$, e.g., a polynomial family. Neural networks are just a special family of approximating functions. A distinctive feature of neural networks is that they have a nonlinear dependence of the approximation function on the coefficients $\\theta$. TensorFlow contains a submodule keras, which makes it easy to build such a network. Below, we build the multilayer perceptrion: a 2 hidden layers 32x32x32x2 network with relu activation functions and linear outputs"
        },
        {
            "cell_type": "code",
            "execution_count": 6,
            "metadata": {},
            "outputs": [],
            "source": [
                "# construction of neural network\n",
                "layers = [\n",
                "    tf.keras.layers.Dense(32, activation='relu', input_dim=5, bias_initializer='he_uniform'),\n",
                "    tf.keras.layers.Dense(32, activation='relu'),\n",
                "    tf.keras.layers.Dense(32, activation='relu'),\n",
                "    tf.keras.layers.Dense(2, activation=tf.keras.activations.linear)\n",
                "]\n",
                "perceptron = tf.keras.Sequential(layers)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 7,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAIECAIAAAChZcOwAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdezxUef8A8A9DRFFR6EHybG2r7MS6ltfGYwnbRUVKSR6yq3ZTbbvsZnsSLVb7pDylki5baF2iZGslSi26PVK61/arxJZpmcZMGJzfH999zs66jGOYC/t5/zXzPWe+5/t1xmfO+Z7vRYmiKEAIoZ4oy7sACKGBAYMFQogRDBYIIUYwWCCEGFERfVNWVvbvf/9bXkVBCCmUdevW2dvb02//dGXx7NmzrKwsmRcJDWbl5eXl5eXyLoVUVFdXD+L/l6ysrGfPnommqHTeKTMzU1blQYOft7c3DNIvVUZGho+Pz6CsGgAoKSl1SME2C4QQIxgsEEKMYLBACDGCwQIhxAgGC4QGj8ePH/P5fClljsECKZzc3FwjI6M7d+7IuyD9IzU1VUlJydfXNy4u7syZM6KbCgoK8vLyyOu9e/euXr3az8/P2dm5pKSEYeY8Hm/EiBFK/zN//nxNTc2KioqEhATRMaJVVVVxcXGrVq1SUlJau3atZBXp4tEpQvKlqak5ZswYdXV16R2itrbWwMBAevl3lpiYqKOjI5qSlJQEACEhIQCQnp7e0tKyY8cOAIiPj3d0dDx16tTMmTN7zDYlJWXBggWmpqbkraurKwBYWFg0NDSEhYV9++23JH3KlClTpkwBgJMnT0pcBQwWSOG4uLi4uLhIL//6+vqlS5eePXtWeofoTEXlT/9r+fn5RUVFdB+NnJwcXV1d8jowMPCLL75ITU3tMVi0tbUdP378zJkzHTIHACcnp4sXL+7cuXPVqlWi6RoaGhJXAW9D0F+LQCBYtGjRL7/8Iscy8Hi8wMDAyMhIOsXS0vL+/fui+3TuE9VZdnZ2ZWVlYGDgkSNHXr9+3WHrunXrNm/e3I81xWCBFEt9fX1KSoqLi0tubi4AXL9+/fPPPzc1NeXz+UFBQbq6ujY2NuQf4Pbt2xs2bDAzM6upqfH09Bw1apSNjQ3pWp6enq6lpWVkZAQAXC43KiqKxWKRYQ45OTl37tzhcDgrVqzYunUrAPz8889GRkanTp2SWR2Tk5PV1NTMzMzolPDw8MLCQvK6qqoKANzd3XvMp7i4mM/nf//9935+fmZmZgUFBaJbNTU1raystmzZ0l/FxmCBFMuvv/5669atwsLCtrY2ANDX179+/frjx4/DwsLWrVtXWFh448aNDRs2AMDhw4d37dp1//79rVu3rlmzZt++fY8fP3Z2dq6trV28eDE9AkpbW/vrr7+ePHkyebtkyRI2m62rq5ucnLx+/XoA4HK5r169qq+vl1kds7KybG1tu9t67NgxGxubhQsX9phPUlISn8+/evXq8uXLa2tr586d26FV2N7ePjs7m/wl+w6DBVIs77zzzty5c+m3+vr61tbWABAZGWlmZjZ16lRra+tr164BQExMjIeHh7KyclxcnKOj4/z585OSkgQCwe7du6HTzbmmpmZ3R/Tw8ODxeL6+vtKq0p+1t7dfvXq1Q2Mnramp6fTp05mZmcrKjP43VVRU3nvvvQMHDmRmZjY3N5MwStPT0+Nyubdv3+6HcmOwQAqoQ3Mdi8USTTQ0NOTxeOS1hoYGi8VSVVUlbz09PdXU1G7evNnbI5JDyEZ9fb1QKBw5cmSXW4uLi8PDw42NjXub7fz58728vK5fvy6aOGLECAB48eKFZEXtAIMFGjxUVFTGjh3b2toq74KIQwJTd7cGNTU1/v7+kuU8Y8aMpqYm0RRyedLe3i5Zhh1gsECDikAgmDRpkrxLIY62tra6unpDQ0OXW01MTJg8B+lOh7r/9ttvAKCvry9xhqIwWKDBo7a2tq6uzsvLCwBUVFQaGxvpH/DGxkb6B1ZZWbmxsVH0g/3128uEkpLStGnTampqutzq7Owscc7nz58PCAgQTeFwOFpaWnTjbh9hsEAKp7a2FgDq6urIWy6XCwD0zcXLly8FAgG9c3Nzc2VlJXkdHR3t7+9vY2MDAObm5g0NDTExMffv34+Ojm5ubr53715FRQUAjB07lsPhXLt27dy5cwKBoLCwcOTIkbKc88rX17e0tLTzkj0nTpwwMTGhqwMAK1eudHBwePjwYedMLly4YGdnl5KS0tzcDAC5ublDhw718/MT3ae0tHTBggX91SKDwQIplqKiosTERABISUk5c+bM2bNnyeiJjRs31tXVHT58+PLly69fv46MjCRXDaqqqocOHVq4cGFQUJCBgUFKSgrJZ82aNbNnz46Li/P393dzc5s+ffrs2bOrq6sBICQkxNDQ0NfXl8PhkCZSTU1NupVUBpYtW6ajo9N5tkGBQNDc3NzS0kKnPH36tKysbN++fZ0zMTY2Hj58+BdffDFt2rSIiAgAOHTokOgOb968KS0tDQsL67dyUyJ++OGHDikI9ZGXl5eXl5eUMg8KClJXV5dS5j1i+P9y5MgRAGhoaBBNvHLlypw5c5gcpaSkJDY2VoLiRURExMfHd0icNGnSmjVrmHwcAH744QfRFLyyQEgW3rx5I/rWysrK19d327Zt4j/F4/Hy8vLIYLNeOXXqlFAoJL3ORAmFwt5mReufgWQ8Hm/48OH9khVCzDU2NgqFQoqi+vIEQTZCQkIcHBymTp1KN2H6+PgUFBScOHFizpw53X3qxo0bmzdv7u0A3MrKSi6XGxsbS6fcunXr9OnTdXV1fRoqInqZIcFtyO7du99///2//e1vvfqUNOTk5BgaGt6+fVveBaHOnDkTGBhI/ryurq5HjhyR9hEzMjLo7sOrV6+uqKiQ9hGZk95tyK5du0g/yKCgoAsXLkjjEOIN7tt26HQb0tdg0dra6uDgoK+v3w+l65uCggJLS8tffvlFeoeoqalhvvPo0aMBoLq6WjblKSsrA4CpU6dK73CSkWqbhXz91YJFX9ssWCyWoaFhHzPpFy4uLteuXRs/fryU8ieTIDDfX0tLCwC0tbVlUx7SsVd6h0MIGzgZkWASBHIXLaV76c7lkerhEAKJg8Xx48eDg4PDwsI+/fRT0oWGoChq9+7dISEhtra2rq6uDx48ALFTEpCtAQEBcXFxc+fOpedH6jIfMRR8EgQZlKdHL168WLFiRVRU1IoVK+bNm/fq1StyHocPH66kpJSQkEAe75eVlRkYGHzzzTfdnYXnz5/HxsZOmTLlt99+mzlz5rhx40hWaPATvSdheA+Wmppqa2v75s0biqLq6up0dXXpNouYmJiDBw9SFNXa2mpmZqavr8/n82traz/44AMAWLVq1a1btyoqKtTU1BYtWkQ+MnHixIsXL1IUJRAIHBwcxOQjpki3b98m05BmZWVRFCXmiOHh4SNGjGCxWGvXri0uLs7OztbV1dXQ0CD3/66uroaGhnS25ubmdnZ25PWsWbNMTEzoTfn5+UOHDk1NTe2uSG+99RYANDY2yqY8d+/eBQBHR8fuyuPo6Ojj40Nes9nspUuXktfh4eEAcOXKFfK2ubnZ1tZWzFk4derUpEmTWCzWv/71r71799rY2Dx//lzMqcE2iwEK+t7AyefzDQwM0tLS6JR58+aRYPH8+XM9Pb22tjaSvnHjRgA4evQoRVFffvklAHA4HLLJwcFhwoQJFEW1tLQoKSlt376dpOfk5IjPR4xz587RwULMESmK8vX1VVVVbWlpIW/JPIgbN26kKMrT01P0n9POzq67f06KolpbW8WURzRYyKA8PQYLJyenb775hrxesmTJu+++S14/e/ZMRUUlKCiIvD158mRUVBQl9iyQZz0PHjwQU30aBosBqnOw6HU/iwsXLtTW1pqbm9Mpampq5EVpaalQKPzoo4/oTUFBQUOHDoWupiQg3d1VVVVnzpy5Zs2aqqqq2NhYT09P8fmI0eMkCHQHe7lMgiCD8ohXVFQEAE1NTampqZcvX6b+NzDB0NDQ29v7yJEjMTExurq6GRkZ//rXv0DsWVBVVVVRUSHRkImsrKxB3JgyiKvWQa+DBfkFGzJkSOdNd+7c0dTUTE5O7lWG2dnZK1asSE5OzsnJycjIcHJykiwfiSnaJAhSKk9bW9u333579erV1atX29raig5MWLt2bXp6+t69e9evX8/hcMi88v14Fuzs7CReq0KRlZWVJSQkkOuLwcfHx6dDSq+DBQkTT548mThxYodNGhoa1dXV1dXVog9T6+rqSI+D7qioqKSmpn744YefffaZm5vb9evXJcunLxRtEoT+Lc+DBw/+9re/zZs3b8yYMdnZ2QDQYWCStbX19OnTd+7cOWnSpNmzZ5PEfjwLhoaGTGaUHIgSEhIGa9U6B4tePw159913AUA0mra3t5Pxf+bm5hRFiY5ye/To0a5du8Tk1tzcvHfvXgDw9fUtLy+nKKq4uFiCfPpCSpMgkOt8qtMwZCmVp7sDURT18ccfV1RUFBQUODo6kkTSRVp0t88++6ympuazzz7z9vYmKTI+C0jx9frKYvr06U5OTgcPHnzvvff8/f1v3bp18eLFurq69PT0OXPmWFtbp6WlNTU1zZs37/Xr18eOHTt69CiInZJg//79ISEhLBZr7Nix2tralpaWtra23eUjhgSTILDZbOg0CUJWVlZMTMzChQszMjKam5ufPXtWUVFhYWFBT4LA4/FsbGzITAEpKSnkv7ozso4Dl8sdNmyYDMpD8u8w/xKXy129evXIkSNJW8mhQ4dsbGyuXLly69atFy9e3LhxQ09PT09PDwDmzJljbGzMZrPpiWRdXFy6OwskfjU0NJCeYOivQrS1k2HrLpfLDQgI0NPTMzY23rRpU3BwcEBAAJm7/dWrV0uWLBkzZszo0aOXLVtGHqoVFhaamJgAwMqVK1++fPn999+T/59Nmzbx+Xxra+uZM2fGxsYGBwcnJyeTQ3SZjxhnz559//33AcDKyqqgoEDMEVtbW4OCgoYMGbJ27Vpvb+/AwMCoqKj29na6arNnzx42bJidnd2VK1eWL1++dOnSEydOUBRVWVlpaGg4ceLEzMxMiqKKiooMDAxyc3M7F6a4uHjlypXkz+vu7n706FFplyc3N9fBwYEckc1mu7q6uri4TJo0idwz7tmzh1xfDB8+3M7OrrCw8Mcff9TV1fXy8qIf1lAU9dFHH5Gq0bo8C3v37iV3In5+fv/97397/Lbg05ABCvp9bMhAJN9JEDpThPK0t7dbWVmRvjP9C4PFANU5WAyktU7FNK3t37+fbplDEjh79uw//vEPqa5FjGTg8ePHY8aMEbNISl8MpGBBt0f0kaJNgiDH8ly8ePGjjz6aPHlyVVVVSUmJjI/+F5Gamrp06dLFixez2WxLS0vRNZ8LCgqam5vJ79zevXurqqrq6+tramr+9a9/kdvqHvF4PCMjI9JiBQBTp06tqKioqKg4f/58aGgo/Y2qqqrKz89/+vTprl271qxZ0+OkO10aSMGiXyQlJZ05c6atrS04ONjf35++1f9rlkdHR6epqem///3vgQMH6FW8B5ba2loDAwNFyES8xMTEDquQJSUlAQCZBSs9Pb2lpWXHjh0AEB8f7+joeOrUqR5XUQeAlJSUBQsWkK4xAODq6goAFhYWDQ0NYWFh3377LUmfMmXKlClTAODkyZOS10H0nmRw34MhuZBqm8Vvv/32j3/8Q16Z9GUOzpMnT4r+Wby9vUNCQshrMjDPz8+vx5xbW1sdHR3JZWlnmzdv/s9//tMhEefgRH9FEswbIKVMeovH4wUGBkZGRtIplpaW9+/fF92HyT1pdnZ2ZWVlYGDgkSNHyKN6UevWrdu8eXM/Vg2DBVIU2dnZn3zyyfr1693d3SMiIshyGMzH6ctl8gHJJCcnq6mpmZmZ0Snh4eGFhYXkdVVVFQC4u7v3mE9xcTGfz//+++/9/PzMzMwKCgpEt2pqalpZWW3ZsqXfyi16mYG3IajfMbwN2bZt27Rp08jQWw6HM2HChBkzZpD+JgzH6cts8gGaxLch9vb23t7e3e0fGhpqY2NDj/cVTygUXr16dfny5crKyurq6h3moI2KitLW1hYdHo23IWhge/nyZURExMcff0yG3uro6Hz11Vfnz59PTU0FAA0NDdGdu3suGBMT4+HhoaysHBcX5+joOH/+/KSkJIFAsHv3buaZAICHhwePx/P19e17vbrU3t5+9erVDo2dtKamptOnT2dmZpI1jXukoqLy3nvvHThwIDMzs7m5ecOGDaJb9fT0uFzu7du3+6HceBuCFEF5eTmfzzc2NqZTZs2aBQDFxcW9ykcukw/0Vn19vVAoHDlyZJdbi4uLw8PDRf8UDM2fP9/Ly+v69euiiaQ//osXLyQragcYLJD8PXnyBP635DdB30H0JVtFm3yAIJGIHhzYQU1Njb+/v2Q5z5gxo6mpSTSFXJ7017LPGCyQ/JE52Tu32/d9nL6iTT4AANra2urq6h2G/NFMTEz60jevQ2VJ/NXX15c4Q1EYLJD82dvba2lpkcmWierqaoFAQJbq6tW8AaKkNPlAHykpKU2bNq27iyZ6vTIJnD9/PiAgQDSFw+FoaWlNnjxZ4jxFYbBA8qejoxMXF/fzzz+fPXuWpOzYscPf39/JyQkAzM3NGxoaYmJi7t+/Hx0d3dzcfO/evYqKCgCgx+mfO3eODPkng/1JJh0G+zPMpLCwcOTIkVlZWdKrr6+vb2lpKdVpCpITJ06YmJjQ5QeAlStXOjg40DMwirpw4YKdnV1KSgp5xpybmzt06FA/Pz/RfchECv3VBIPBAimEjz/+OCcn59tvv/300083btyor69/4MABsmnNmjWzZ8+Oi4vz9/d3c3ObPn367Nmzq6urASAkJMTQ0NDX15fD4ZDnHaqqqocOHVq4cGFQUJCBgUFKSkpvM2GxWJqamnQrqTQsW7ZMR0dHdGZDQiAQNDc3kzUZiKdPn5aVlXWY2YwwNjYePnz4F198MW3atIiICAA4dOiQ6A5v3rwpLS0Vnb6or0Sfo2I/C9TvZDlEXcaD/fvS3fvKlStz5sxhcpSSkpLY2FgJihcREREfH98hEftZIKTo3rx5I/rWysrK19e3x9GfPB4vLy+PDDbrlVOnTgmFwvXr13dIFwqFvc2K9pcbdYoGMUWbfEBUSEiIg4PD1KlT6SZMHx+fgoKCEydOkHbcLt24cWPz5s29nWeksrKSy+XGxsbSKbdu3Tp9+nRdXV1fhopgsECDhKJNPkBbsmTJkiVLutxERpSLMX36dAmOyGazyXyutMmTJ5NnIqIRpLcwWKBBIiQkRILLdcQctlkghBjBYIEQYgSDBUKIEQwWCCFGumjgzMjIkH050GBFekkOyi9VWVkZDNKqdU20h9ZgXQ8aISSBDj04lajer9yLBjeyLPhf6AcTMYNtFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRjBYIIQYwWCBEGIEgwVCiBEMFgghRpQoipJ3GZCcpaampqSktLe3k7ePHz8GgPHjx5O3ysrKgYGBS5YskVv5kGLAYIHgxo0bbDZbzA6VlZXvvvuuzMqDFBMGCwQAMGnSpHv37nW56a233nrw4IGMy4MUELZZIAAAPz8/VVXVzumqqqoBAQGyLw9SQHhlgQAAfvnll7feeqvLL8ODBw/eeust2RcJKRq8skAAAKamppaWlkpKSqKJSkpKVlZWGCkQgcEC/W7ZsmUsFks0hcViLVu2TF7lQYoGb0PQ716+fGlgYEA/QAUAZWXlmpoaPT09OZYKKQ68skC/GzNmzIwZM+iLCxaL5ejoiJEC0TBYoD/4+fmJXmn6+fnJsTBI0eBtCPrD69evR48e3dLSAgCqqqovX74cMWKEvAuFFAVeWaA/aGlpubm5qaioqKioeHh4YKRAojBYoD9ZunRpW1tbW1sbDgZBHeBtCPqTpqYmXV1diqI4HM7QoUPlXRykQBQoWHToEYQQAgDF+Q9VkXcB/mTNmjX29vbyLsWgVVZWlpCQ8MMPP4jf7fr160pKSuLHoSogHx+fQfb9IedL3qUQQSkMAPjhhx/kXYrBjISJHncTCoVCoVAG5elfg+/7w/B8yYxiXVkgRaCigt8K1AV8GoIQYgSDBUKIEQwWCCFGMFgg1P8eP37M5/PlXYp+hsEC9SA3N9fIyOjOnTvyLogUFRQU5OXlkdd79+5dvXq1n5+fs7NzSUkJwxx4PN6IESOU/mf+/PmampoVFRUJCQmUwnSU6CNs90Y90NTUHDNmjLq6uvQOUVtba2BgIL38xUtKSgKAkJAQAEhPT29padmxYwcAxMfHOzo6njp1aubMmT1mkpKSsmDBAlNTU/LW1dUVACwsLBoaGsLCwr799lspVkBWMFigHri4uLi4uEgv//r6+qVLl549e1Z6hxAjPz+/qKgoMzOTvM3JydHV1SWvAwMDv/jii9TU1B6DRVtb2/Hjx8+cOdP5qbOTk9PFixd37ty5atWqfi+8jOFtCJIngUCwaNGiX375RS5H5/F4gYGBkZGRdIqlpeX9+/dF92EyCiE7O7uysjIwMPDIkSOvX7/usHXdunWbN2+WVx37EQYLJE59fX1KSoqLi0tubi4AXL9+/fPPPzc1NeXz+UFBQbq6ujY2NuTf4Pbt2xs2bDAzM6upqfH09Bw1apSNjU15eTkApKena2lpGRkZAQCXy42KimKxWKRfdk5Ozp07dzgczooVK7Zu3QoAP//8s5GR0alTp2RQu+TkZDU1NTMzMzolPDy8sLCQvK6qqgIAd3f3HvMpLi7m8/nff/+9n5+fmZlZQUGB6FZNTU0rK6stW7b0a9nlQd5dSP8Ag667rqKRoPvw7du3165dCwBZWVkURdXW1n7wwQcAsGrVqlu3blVUVKipqS1atIiiqPDw8BEjRrBYrLVr1xYXF2dnZ+vq6mpoaNTU1FAU5erqamhoSGdrbm5uZ2dHXs+aNcvExITelJ+fP3To0NTU1N7WToLvj729vbe3d3dbQ0NDbWxs2tramGQlFAqvXr26fPlyZWVldXX127dvi26NiorS1tZubW3tVfEUrbs3Xlkgcd555525c+fSb/X19a2trQEgMjLSzMxs6tSp1tbW165dA4CYmBgPDw9lZeW4uDhHR8f58+cnJSUJBILdu3cDgIaGhmi2mpqa3R3Rw8ODx+P5+vpKq0r/097efvXqVR0dnS63NjU1nT59OjMzU1mZ0f+IiorKe++9d+DAgczMzObm5g0bNohu1dPT43K5t2/f7odyyw8GC9SDDo12ZEZfOtHQ0JDH45HXGhoaLBaLXtnM09NTTU3t5s2bvT1ihxUJpKS+vl4oFI4cObLLrcXFxeHh4cbGxr3Ndv78+V5eXtevXxdNJHOOvXjxQrKiKggMFkhaVFRUxo4d29raKu+CdI2EpLa2ti631tTU+Pv7S5bzjBkzmpqaRFPI5YnoMgsDEQYLJEUCgWDSpEnyLkXXtLW11dXVGxoautxqYmLSl9mYOtT6t99+AwB9fX2JM1QEGCyQtNTW1tbV1Xl5eQGAiopKY2Mj/TPe2NhI/8wqKys3NjaKflA2v8BKSkrTpk2rqanpcquzs7PEOZ8/f77DatIcDkdLS2vy5MkS56kIMFigHtTW1gJAXV0decvlcgGAvrl4+fKlQCCgd25ubq6srCSvo6Oj/f39bWxsAMDc3LyhoSEmJub+/fvR0dHNzc337t2rqKgAgLFjx3I4nGvXrp07d04gEBQWFo4cOTIrK0sGVfP19S0tLaU6dcc+ceKEiYkJXREAWLlypYODw8OHDztncuHCBTs7u5SUlObmZgDIzc0dOnRohyVXSktLFyxYIJu2GOnBYIHEKSoqSkxMBICUlJQzZ86cPXuWjKHYuHFjXV3d4cOHL1++/Pr168jISHLVoKqqeujQoYULFwYFBRkYGKSkpJB81qxZM3v27Li4OH9/fzc3t+nTp8+ePbu6uhoAQkJCDA0NfX19ORwOaSLV1NSkW0mlatmyZTo6OqQziCiBQNDc3EzWTyGePn1aVla2b9++zpkYGxsPHz78iy++mDZtWkREBAAcOnRIdIc3b96UlpaGhYVJoQayJe9nt38A7GchZdJ+bh8UFKSuri69/MWT7Ptz5cqVOXPmMNmzpKQkNja29+WiIiIi4uPjJfgg9rPoZ/RzO4QkYGVl5evru23bNvG78Xi8vLw8MtisV06dOiUUCtevXy9pARXIAA4We/bsmTFjxjvvvCPvgsCxY8ecnJzI2ORp06Y5ODhYWFjY2dmFhYU9evRI3qWTncbGRjLTr7wL0js+Pj6TJ08+ceKEmH1u3LixefNmLS2tXuVcWVnJ5XJjY2P7VkCFIe9Lmz9ALy8jW1tbHRwc9PX1pVck5sjt97hx4+iUy5cvu7m5sVisr776imGXYWmT6mXtrl27SG/IoKCgCxcuSOkoYvT2+6P4FO02ZAAPUWexWIaGhl02UMse6b8suoSXtbV1fn6+n5/fN998M2zYsC+//FJ+pZOFkJAQCa7S0QAygG9DFEqXHXiUlZV37do1ZsyY6Ojop0+fyr5UCPWjgRcsjh8/HhwcHBYW9umnn5IuAARFUbt37w4JCbG1tXV1dX3w4AGIHVJNtgYEBMTFxc2dO5ee36XLfEDSodPa2toLFy4UCAQZGRkyKCRCUiTfuyBRwOCeMzU11dbW9s2bNxRF1dXV6erq0m0WMTExBw8epCiqtbXVzMxMX1+fz+eLGVJNUdTEiRMvXrxIUZRAIHBwcBCTD9XT0GnSa3jSpEmdNx05cgQAAgICZFBI8RTtHrh/Mfn+DCyKdr4UqSg9nWw+n29gYJCWlkanzJs3jwSL58+f6+np0e2IGzduBICjR49SFEUaCzgcDuTW6ioAACAASURBVNnk4OAwYcIEiqJaWlqUlJS2b99O0nNycsTnQ1GUmPkIxASLn376CQCcnZ1lU0gxFO3L178wWEjbQGrgvHDhQm1trbm5OZ2ipqZGXpSWlgqFwo8++ojeFBQURJobOw+pJm2iqqqqM2fOXLNmTVVVVWxsrKenp/h8QNKh06R/9MSJE2VTyB6RG6JBqaysTN5F6E+KVp2BFCzu3r0LAEOGDOm86c6dO5qamsnJyb3KMDs7e8WKFcnJyTk5ORkZGU5OTpLlIx6ZRJ/NZitIIX18fCT7oOJLSEhQrGXHB5eB1MBJwsSTJ086b9LQ0KiuriadHWj02KfuqKiopKampqamqqiouLm53blzR7J8xKAoKisrS1VV1c3NTUEKKe+LWWmBQXobojgGUrB49913AUD0L9je3k7GL5mbm1MUJTpW59GjR7t27RKTW3Nz8969ewHA19e3vLycoqji4mLx+YgZOk11023xu+++u3nzZlhY2Lhx42RTSISkRY6BswNg8Mvg5OTEYrF27drF5/MvX748duxYAEhLS2tsbCRzQ86fP//w4cM7d+50dnauq6ujKOrTTz8FkbbDf/zjH1paWhRFNTU1WVhYkDbLlpYWXV3dsrKy9vb27vI5c+aMlpZWZmZmlwV79uwZABgbG9Mp//d///fpp58qKSmFhoaSxkgxmfdXIcVTtAaz/sXk+zOwKNr5UqSiMDjZXC43ICBAT0/P2Nh406ZNwcHBAQEBhYWFbW1tr169WrJkyZgxY0aPHr1s2bLnz59TFFVYWGhiYgIAK1eufPny5ffffz9s2DAA2LRpE5/Pt7a2njlzZmxsbHBwcHJyMjlEl/lQFFVUVGRgYJCbm9u5VLm5uU5OTiT4Ojg4ODs7e3h4uLu7r1u3rrKyUnRPaRdSPEX78vUvDBbSpkhFGXQnW9Eo2pevfw2+74+ina+B1GaBEJIjDBYIIUYwWKC/ooKCAjI/IADs3bt39erVfn5+zs7OJSUlDHPIzMy0tLQcNmwYm80+fvy46Ka0tDQrKystLS1bW9sff/yRJFZUVCQkJFADbbIPURgsUL8RHdcn30zES0pKevTo0ezZswEgPT29paVlx44dhw8fdnNzc3R0JN3zxTt48GBJScmBAwfy8vJYLJa3tzc9lm/btm1Hjhzx8/P75z//WVVVNWvWLLJ4qoWFBZvNHtAzcWKwQP2jvr5+6dKlipCJePn5+UVFRfTUGzk5OfSqgoGBgRRFpaamis9BKBQ+fPgwMTGRzWY7OTnt27dPKBReunQJABobG0+ePJmfnx8aGpqQkFBYWKikpBQfH08+6OTkNHz48J07d0qtctI1kLp7I4UlEAgWLVpED6uXYybi8Xi8wMDAoqIiOsXS0pJeNp3ocW0hZWXlTZs20W/J/GCk58ulS5diY2PpHOzt7S0sLETnZ1q3bp2pqam7u7upqWlfKyNzeGWBupCdnf3JJ5+sX7/e3d09IiKCrIiRnp6upaVlZGQEAFwuNyoqisVi2dvbA0BOTs6dO3c4HM6KFSu2bt16+/btDRs2mJmZ1dTUeHp6jho1ysbGhsy4zzwTkHQOETGSk5PV1NTMzMzolPDwcDpYVFVVAYC7u7v4TFgsluj6r2lpaYmJiW+//TYAODs7k6hB09bWJp1oCE1NTSsrqy1btvS1JnIh72e3f4BB95xc0TB8br9t27Zp06a1tLRQFMXhcCZMmDBjxoz29naKolxdXQ0NDek9zc3N7ezsyOtZs2aZmJiQ1+Hh4SNGjGCxWGvXri0uLs7OztbV1dXQ0KipqWGeCdXTHCIdMPn+2Nvbe3t7d7c1NDTUxsaG+YSpPB4vMjJST0/vp59+6nKH1tbW0aNH79+/XzQxKipKW1tbzHQHNOxngRTay5cvIyIiPv74Y7LMj46OzldffXX+/HlyJ6+hoSG6M5l5tLOYmBgPDw9lZeW4uDhHR8f58+cnJSUJBILdu3czzwQAPDw8eDyer69v3+sFAO3t7VevXiV3DZ01NTWdPn06MzOTLGLcIz6fHxkZeenSpfr6+pkzZ+7fv7/zPsePH586dery5ctFE/X09LhcLt1QMoBgsEB/Ul5ezufzjY2N6ZRZs2YBQHFxca/yIWuL0QuLeXp6qqmp3bx5s7fl6ccl/+rr64VC4ciRI7vcWlxcHB4eLlpx8TQ1NePj4/Pz869duzZq1KjOdxb19fXR0dGHDx/u0AgyYsQIAHjx4kXvayBnGCzQn5AZAMiq3wR9B9GXbFVUVMaOHUuvkCoXJO7QizN3UFNT4+/vL0G2U6ZMCQ0Nffz4sVAoFE1fu3ZtQkKCnp5eh/3JlYtsFn/uXxgs0J+MHz8eADo/kpg0aVIfcxYIBH3PpC+0tbXV1dXJBIidmZiY9PgcpDtTpkwxNDQUXZ91586dnp6e77//fuedSSDW19eX7FhyhMEC/Ym9vb2WllZubi6dUl1dLRAI5syZAwAqKiqNjY30j3NjYyP9C6msrNzY2NhdtrW1tXV1dV5eXr3NpB9/gcl6cd1dIjk7O0uc8927d8nfh0hLSxs6dCiZA5EQfTrL4XC0tLQmT54s8eHkBYMF+hMdHZ24uLiff/757NmzJGXHjh3+/v5kDL65uXlDQ0NMTMz9+/ejo6Obm5vv3btXUVEBAGPHjuVwONeuXTt37pxAIACA5ubmyspKkkl0dLS/v7+NjU2vMiksLBw5cmRWVlZ/1c7X17e0tJTq1Of6xIkTJiYmdGkBYOXKlQ4ODl0uYdXQ0BAQEHDs2DGSz8OHD8+fPx8XF0e2/vjjj4mJiUKhcM+ePXv27Nm9e/fKlSvJjJBEaWnpggUL+rEtRnbk/TjmD4CPTqWM+aO43NxcV1fXTz755Ouvv/7uu+/Ic1OKorhc7uzZs4cNG2ZnZ3flypXly5cvXbr0xIkTFEVVVlYaGhpOnDiRzA8UFBQ0ZMiQtWvXent7BwYGRkVFSZCJmDlEOmPy/WlpaZkwYQKJF6LS09P19fUvX75Mp3z44YfKysphYWGdM+HxeLNmzdLR0Xn//fejoqKOHDlCVnilKOry5cudZ05WU1N79eoV2UEgEIwaNeru3btMaqRoj04VqSgYLKRMll++oKAgdXV12RyLYPj9uXLlypw5c5hkWFJSEhsb2+dy/UlERER8fDzDnRUtWOBtCPprsbKy8vX13bZtm/jdeDxeXl5e/67eeurUKaFQuH79+n7MU5YwWCCpaGxsJBfn8i5IF3x8fCZPnnzixAkx+9y4cWPz5s1aWlr9ddDKykoulxsbG9tfGcoeDiRD/S8pKenMmTNtbW3BwcH+/v4ODg7yLlFHrq6u4neYPn16/x6RzWaz2ez+zVPGMFig/hcSEtK/F/BIEeBtCEKIEQwWCCFGMFgghBjBYIEQYkRJcR5uKSkp2dnZGRoayrsgg1Z1dXV5eTkZoDH4ZGVlDbLvDzlfCvQfqjhF8fb2lncREAAAGaZhYWEh74IgAIDMzEx5F+F3ChQskIJYuHAhAGRkZMi7IEixYJsFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiBIMFQogRDBYIIUYwWCCEGMFggRBiREXeBUDyJxAImpub6bctLS0AUF9fT6eoqalpaGjIoWRIkShRFCXvMiA527Vr16pVq8TssHPnzpUrV8qsPEgxYbBAUFdXZ2Bg0NbW1uVWFotVW1s7evRoGZcKKRpss0AwevRoZ2dnFovVeROLxfrggw8wUiDAYIGIpUuXdnmNSVHU0qVLZV8epIDwNgQBAPB4vNGjR4s2cxJDhgypq6vT0tKSS6mQQsErCwQAMHz48NmzZ6uqqoomqqiozJ07FyMFIjBYoN8tWbKktbVVNKWtrW3JkiXyKg9SNHgbgn7X0tKiq6vL4/HolGHDhnE4HDU1NTmWCikOvLJAvxsyZIi3t/eQIUPIW1VVVR8fH4wUiIbBAv3B19eXdN8EAKFQ6OvrK9/yIIWCtyHoD+3t7fr6+nV1dQCgq6v766+/dtn5Av014ZUF+oOysrKvr++QIUNUVVWXLFmCkQKJwmCB/mTx4sUtLS14D4I6k/Wo07KysmfPnsn4oIg5iqJ0dHQA4PHjx//3f/8n7+KgbhkZGdnb28v0kJRseXl5ybR6CA1SXl5eMv7nlcN8Fl5eXpmZmbI/LhLD29sbAMh5uX37NgCYmZnJuUz9JCMjw8fHhxpcDfnkfMkYTn6DOho0YQL1L2zgRAgxgsECIcQIBguEECMYLBBCjGCwQAgxgsECSSg3N9fIyOjOnTvyLkg/KygoyMvLI6/37t27evVqPz8/Z2fnkpIShjlkZmZaWloOGzaMzWYfP35cdFNaWpqVlZWWlpatre2PP/5IEisqKhISEhT/4S4GCyQhTU3NMWPGqKurS+8QtbW10su8S0lJSY8ePZo9ezYApKent7S07Nix4/Dhw25ubo6Ojj/99FOPORw8eLCkpOTAgQN5eXksFsvb2/vBgwdk07Zt244cOeLn5/fPf/6zqqpq1qxZhYWFAGBhYcFms8PCwqRatb7DYIEk5OLicu3atfHjx0sp//r6ehnPFZyfn19UVBQSEkLe5uTkkP5pABAYGEhRVGpqqvgchELhw4cPExMT2Wy2k5PTvn37hELhpUuXAKCxsfHkyZP5+fmhoaEJCQmFhYVKSkrx8fHkg05OTsOHD9+5c6fUKtcPsFMWUkQCgWDRokW//PKLzI7I4/ECAwOLioroFEtLS/LLT1NSUhKfibKy8qZNm+i3ZJSNtbU1AFy6dCk2NpbOwd7e3sLC4uHDh/TO69atMzU1dXd3NzU17WtlpAOvLJAk6uvrU1JSXFxccnNzAeD69euff/65qakpn88PCgrS1dW1sbEh/+q3b9/esGGDmZlZTU2Np6fnqFGjbGxsysvLASA9PV1LS8vIyAgAuFxuVFQUi8Uig6NycnLu3LnD4XBWrFixdetWAPj555+NjIxOnTolpRolJyerqamJ9l4NDw+ng0VVVRUAuLu7i8+ExWKpqPzxA5yWlpaYmPj2228DgLOzM4kaNG1tbRMTE/qtpqamlZXVli1b+loTqcFggSTx66+/3rp1q7CwkKxjpq+vf/369cePH4eFha1bt66wsPDGjRsbNmwAgMOHD+/atev+/ftbt25ds2bNvn37Hj9+7OzsXFtbu3jxYnrcpLa29tdffz158mTydsmSJWw2W1dXNzk5ef369QDA5XJfvXolugJr/8rKyrK1te1u67Fjx2xsbBYuXMgwt8bGxs2bN2/fvn3ixIld7tDW1nbz5s0O91n29vbZ2dndLQ0ndxgskCTeeeeduXPn0m/19fXJz2ZkZKSZmdnUqVOtra2vXbsGADExMR4eHsrKynFxcY6OjvPnz09KShIIBLt37waADusta2pqdndEDw8PHo8npVk22tvbr169Su4aOmtqajp9+nRmZqayMqP/Fz6fHxkZeenSpfr6+pkzZ+7fv7/zPsePH586dery5ctFE/X09LhcLt1QomgwWCAJiV5vAwCZVotONDQ0pCcK19DQYLFY9KIknp6eampqN2/e7O0RpTdzV319vVAoHDlyZJdbi4uLw8PDjY2NGeamqakZHx+fn59/7dq1UaNGdb6zqK+vj46OPnz4cIdGkBEjRgDAixcvel8DWcBggWRNRUVl7NixHdYokS8Shrq7/q+pqfH395cg2ylTpoSGhj5+/FgoFIqmr127NiEhQU9Pr8P+5Mqlvb1dgmPJAAYLJAcCgWDSpEnyLsUftLW11dXVGxoautxqYmLS43OQ7kyZMsXQ0FB0qbedO3d6enq+//77nXf+7bffAEBfX1+yY0kbBgska7W1tXV1dWTONBUVlcbGRvonvbGxkf5dVVZWbmxsFP2g9H5ylZSUpk2bVlNT0+VWZ2dniXO+e/funDlz6LdpaWlDhw719PSkU0SfznI4HC0tLbqVV9FgsEASIt0ryboBAMDlcgGAvrl4+fKlQCCgd25ubq6srCSvo6Oj/f39bWxsAMDc3LyhoSEmJub+/fvR0dHNzc337t2rqKgAgLFjx3I4nGvXrp07d04gEBQWFo4cOTIrK0tK1fH19S0tLe3c5/rEiRMmJiZ04QFg5cqVDg4Ool0kaA0NDQEBAceOHSP5PHz48Pz583FxcWTrjz/+mJiYKBQK9+zZs2fPnt27d69cufLu3bv0x0tLSxcsWKCwk6pjpywkiaKiosTERABISUn5+9//rqysTMZTbNy4cdOmTadPn758+XJjY2NkZGRERAQAqKqqHjp0qLq6WktLy8TEhDxVBYA1a9ZcvXo1Li4uPz8/MTHx0aNHra2t1dXVFhYWISEhJ0+e9PX13bJlC2ki1dTU7LB0cz9atmxZXFxceXl5h1lwBQJBc3MzvfYSADx9+rSsrGzfvn2xsbEdMlFRUeFwOMHBwdu3b3dxcRk/fnx+fj5p9L1y5YqXl9ebN29IHxNCTU2Nvpx58+ZNaWlpaWmplCrYD2Q856eXl5fsJxpFPZLqeQkKClJXV5dS5j364YcfGH7Pr1y5MmfOHCZ7lpSUxMbG9q1cHUVERMTHxzPcWS7/RwPmNkR0wV6EpMHKysrX13fbtm3id+PxeHl5efQQkn5x6tQpoVBIup8prAEQLPbs2TNjxox33nlH3gX5XUNDQ0RExJdffslw/2PHjjk5OSkpKZFWNAcHBwsLCzs7u7CwsEePHkm1qAqisbFRKBRSCj8EGwB8fHwmT5584sQJMfvcuHFj8+bNWlpa/XXQyspKLpfb+aZG0QyAYBEUFNTe3q4gfWDz8vI++uijLVu2dGioF2P+/PlHjhwBgHHjxpWWll68eLGioiIxMfHGjRtvv/32hg0bFPa5er9ISko6c+ZMW1tbcHDwxYsX5V2cnrm6uoo+v+hs+vTp/Tswn81mL1q0qB8zlJIBECxYLJahoaG8S/G72bNnJycn9/ZTpBfz0KFD6RRra+v8/HwfH59vvvmGbi0flEJCQjgcDkVRycnJDg4O8i4OktwACBaKRk1Nrbcf6bJLj7Ky8q5du8aMGRMdHf306dP+KBpCUqS4weL48ePBwcFhYWGffvqp6IxJFEXt3r07JCTE1tbW1dWVTEMkZog02RoQEBAXFzd37lwXFxcx+UhMsgHU2traCxcuFAgEGRkZCls1hH4n46cvDB/5pKam2travnnzhqKouro6XV1dfX19sikmJubgwYMURbW2tpqZmenr6/P5/Nra2g8++AAAVq1adevWrYqKCjU1tUWLFpGPTJw48eLFixRFCQQCBwcHMfkwqUJTUxMAfPLJJ6KJ+fn5Q4cOTU1N7fIjpB/xpEmTOm8izRkBAQHyrdogfqTN/NHpACKX86WIwYLP5xsYGKSlpdEp8+bNI8Hi+fPnenp6bW1tJH3jxo0AcPToUYqiyOMJcntMUZSDg8OECRMoimppaVFSUtq+fTtJz8nJEZ9Pj7oMFhRFtba2dvcRMcGCTOvo7Ows36phsBhY5HK+FLEH54ULF2pra83NzekUupmgtLRUKBR+9NFH9KagoCDScNh5iDTpkKuqqjpz5sw1a9ZUVVXFxsaSbvli8pGYZL10SS/piRMnyr1q5eXlclluV9qqq6tBTisJS095ebmdnZ2MD6qIwYL0lh8yZEjnTXfu3NHU1Ozt84js7OwVK1YkJyfn5ORkZGQ4OTlJlo80kKn02Wz24KsaGmQUMViQMPHkyZPOU5JpaGhUV1dXV1eLPkytq6sbPXq0mAxVVFRSU1M//PDDzz77zM3N7fr165Ll0+8oisrKylJVVXVzc8vKypJv1ezs7DIzM/tYIwWUkZHh4+MzyKomlwslRXwa8u677wIAudUk6E5Z5ubmFEWJrrDw6NGjXbt2icmtubl57969AODr61teXk5RVHFxsQT59EhM3yqqm86L33333c2bN8PCwsaNG6fIVUMIFPPKYvr06U5OTgcPHnzvvff8/f1v3bp18eLFurq69PT0OXPmWFtbp6WlNTU1zZs37/Xr18eOHTt69CiIHSK9f//+kJAQFos1duxYbW1tS0tLW1vb7vLpEZ/Ph06zKhUWFi5YsCAlJYVM09AB6e4pOmT7yZMn33333X/+85/Q0NDIyEgAcHFxkXvVEBJHxg2qDFtxuVxuQECAnp6esbHxpk2bgoODAwICyFzSr169WrJkyZgxY0aPHr1s2bLnz59TFFVYWEhmVV+5cuXLly+///77YcOGAcCmTZv4fL61tfXMmTNjY2ODg4OTk5PJIbrMp0cFBQVkRmZTU9M9e/bU1NSQ9KKiIgMDg9zc3M4fyc3NdXJyIn9tBwcHZ2dnDw8Pd3f3devWVVZWiu4px6rh05CBRS7nS4mS7fAecq81yG4gB4FBfF5Im4WMv+fSJpfzpYi3IXIkphVw//79ZAlMhP6aFLGBU47quoeR4i8CV1HvDgYLJAv9sh66DBZVx1XUxcBggaSuX9ZDl8Gi6riKunjYZoGkq1/WQ5fBouq4inqP8MoC9U52dvYnn3yyfv16d3f3iIiI5uZm6M166Aq7qDquot4zGT+qHcTP8wc0hudl27Zt06ZNa2lpoSiKw+FMmDBhxowZ7e3tFEW5uroaGhrSe5qbm9vZ2ZHXs2bNMjExIa/Dw8NHjBjBYrHWrl1bXFycnZ2tq6uroaFBeqwwzITqaU4AUQz7Wdjb23t7e3e3NTQ01MbGhh7L2yMejxcZGamnp/fTTz91uUNra+vo0aP3798vmhgVFaWtrS1m+DINZ/dGCu3ly5cREREff/wxWbxDR0fnq6++On/+PLmTZ7geumIuqo6rqDOBwQIxVV5ezufzRRcTnzVrFgAUFxf3Kh8FXFQdV1FnAoMFYurJkyfwv8V7CfoOoi/ZKsKi6riKOhMYLBBT48ePB4DOjyT6vh663BdVx1XUmcBggZiyt7fX0tLKzc2lU6qrqwUCAVllo1froYtShEXVcRV1JjBYIKZ0dHTi4uJ+/vnns2fPkpQdO3b4+/uTMbXM10MHhVxUHVdR7xF2ykK98PHHHxsYGHz77be5ubkjR47U19en/xOYr4cOCrmoOq6i3jMZP6rFfhaKSZbnRcaLquMq6v0Fb0MQ+h2uoi4eBgska4q8qDquoi4GtlkgmRJdVN3f318Bl0p2dXUVv8P06dP794hsNpvNZvdvntKAwQLJVEhISP9ewCOZwdsQhBAjGCwQQoxgsEAIMYLBAiHECAYLhBAjcngakpWVJfEYPiRVg/i8DL6qdblQplTJekWysrKyZ8+eyfKIqLdIF8a1a9fKuyBIHCMjow7DWKRN1sECKb6FCxcCQEZGhrwLghQLtlkghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGMFgghBjBYIEQYgSDBUKIEQwWCCFGVORdACR/ly5dqqyspN/+8ssvALB37146hc1m29rayqFkSJEoURQl7zIgOTt58uTs2bNZLJaysjIAkK+EkpISALS3t7e1teXl5c2aNUvOpUTyhsECgVAo1NXVff36dZdbtbS06urqhgwZIuNSIUWDbRYIVFVVFy9e3GU4ELMJ/dVgsEAAAIsXL25paemcLhQKfX19ZV8epIDwNgQBALS3t48dO/bFixcd0kePHv3rr7+Stgz0F4dfAgQAoKys7Ofn1+F2Y8iQIcuXL8dIgQj8HqDfdb4TaWlpWbx4sbzKgxQN3oagP0yYMOHhw4f0W1NT00ePHsmxPEih4JUF+sPSpUtVVVXJ6yFDhvj7+8u3PEih4JUF+sPDhw8nTJhAv713797EiRPlWB6kUPDKAv3hrbfeYrPZSkpKSkpKbDYbIwUShcEC/cmyZctYLBaLxVq2bJm8y4IUC96GoD+pqakxMjKiKOrZs2d/+9vf5F0cpEBkHSz+/e9/l5WVyfKIqLfOnTsHAI6OjnIuBxLL3t5+3bp1sjyirG9DysrKysvLZXxQ1KPy8nL6vBgbG48bN06+5elH1dXVWVlZ8i5FPysvL5f9j64c5rOws7PLzMyU/XGRGN7e3gBAzstvv/0GAKNGjZJzmfpJRkaGj4/PIPvKkfMlYzj5Depo0IQJ1L/waQhCiBEMFgghRjBYIIQYwWCBEGIEgwWSUG5urpGR0Z07d+RdkH5WUFCQl5dHXu/du3f16tV+fn7Ozs4lJSUMc8jMzLS0tBw2bBibzT5+/LjoprS0NCsrKy0tLVtb2x9//JEkVlRUJCQkKH73SAwWSEKamppjxoxRV1eX3iFqa2ull3mXkpKSHj16NHv2bABIT09vaWnZsWPH4cOH3dzcHB0df/rppx5zOHjwYElJyYEDB/Ly8lgslre394MHD8imbdu2HTlyxM/P75///GdVVdWsWbMKCwsBwMLCgs1mh4WFSbVqfYfBAknIxcXl2rVr48ePl1L+9fX1S5culVLmXcrPzy8qKgoJCSFvc3Jybt++TV4HBgZSFJWamio+B6FQ+PDhw8TERDab7eTktG/fPqFQeOnSJQBobGw8efJkfn5+aGhoQkJCYWGhkpJSfHw8+aCTk9Pw4cN37twptcr1A+xngRSRQCBYtGgRWe5INng8XmBgYFFREZ1iaWlJfvlpZC0VMZSVlTdt2kS/1dHRAQBra2sAuHTpUmxsXoMnWQAAIABJREFULJ2Dvb29hYWF6FRD69atMzU1dXd3NzU17WtlpAOvLJAk6uvrU1JSXFxccnNzAeD69euff/65qakpn88PCgrS1dW1sbEh/+q3b9/esGGDmZlZTU2Np6fnqFGjbGxsSNfy9PR0LS0tIyMjAOByuVFRUSwWy97eHgBycnLu3LnD4XBWrFixdetWAPj555+NjIxOnTolpRolJyerqamZmZnRKeHh4XSwqKqqAgB3d3fxmbBYLBWVP36A09LSEhMT3377bQBwdnYmUYOmra1tYmJCv9XU1LSystqyZUtfayI1GCyQJH799ddbt24VFha2tbUBgL6+/vXr1x8/fhwWFrZu3brCwsIbN25s2LABAA4fPrxr16779+9v3bp1zZo1+/bte/z4sbOzc21t7eLFi0loAABtbe2vv/568uTJ5O2SJUvYbLaurm5ycvL69esBgMvlvnr1qr6+Xko1ysrKErNE47Fjx2xsbBYuXMgwt8bGxs2bN2/fvr27OUHa2tpu3rzZ4T7L3t4+Ozub/EkVEAYLJIl33nln7ty59Ft9fX3ysxkZGWlmZjZ16lRra+tr164BQExMjIeHh7KyclxcnKOj4/z585OSkgQCwe7duwFAQ0NDNFtNTc3ujujh4cHj8aS0iEl7e/vVq1fJXUNnTU1Np0+fzszMZDjROZ/Pj4yMvHTpUn19/cyZM/fv3995n+PHj0+dOnX58uWiiXp6elwul24oUTQYLJCERK+3AYDFYokmGhoa8ng88lpDQ4PFYtGze3p6eqqpqd28ebO3RySHkIb6+nqhUDhy5MgutxYXF4eHhxsbGzPMTVNTMz4+Pj8//9q1a6NGjep8Z1FfXx8dHX348OEOjSAjRowAgM6rtygIDBZI1lRUVMaOHdva2irvgvyBhKHurv9ramokm7t4ypQpoaGhjx8/FgqFoulr165NSEjQ09PrsD+5cmlvb5fgWDKAwQLJgUAgmDRpkrxL8QdtbW11dfWGhoYut5qYmPT4HKQ7U6ZMMTQ0pK+qAGDnzp2enp7vv/9+553J5AD6+vqSHUvaMFggWautra2rq/Py8gIAFRWVxsZG+ie9sbGR/l1VVlZubGwU/aD0fnKVlJSmTZtWU1PT5VZnZ2eJc7579+6cOXPot2lpaUOHDvX09KRTRJ/OcjgcLS0tupVX0WCwQBIi3Svr6urIWy6XCwD0zcXLly8FAgG9c3Nzc2VlJXkdHR3t7+9vY2MDAObm5g0NDTExMffv34+Ojm5ubr53715FRQUAjB07lsPhXLt27dy5cwKBoLCwcOTIkdKb88rX17e0tLRzn+sTJ06YmJjQhQeAlStXOjg4iHaRoDU0NAQEBBw7dozk8/Dhw/Pnz8fFxZGtP/74Y2JiolAo3LNnz549e3bv3r1y5cq7d+/SHy8tLV2wYIH0mmb6CDtlIUkUFRUlJiYCQEpKyt///ndlZWUynmLjxo2bNm06ffr05cuXGxsbIyMjIyIiAEBVVfXQoUPV1dVaWlomJibkqSoArFmz5urVq3Fxcfn5+YmJiY8ePWptba2urrawsAgJCTl58qSvr++WLVtIE6mmpqbo9Xz/WrZsWVxcXHl5Of00lxAIBM3NzaILOz59+rSsrGzfvn2xsbEdMlFRUeFwOMHBwdu3b3dxcRk/fnx+fj5p9L1y5YqXl9ebN29Ep5VUU1OjL2fevHlTWlpaWloqpQr2A0q2vLy8vLy8ZHxQ1COpnpegoCB1dXUpZd6jH374geH3/MqVK3PmzGGyZ0lJSWxsbN/K1VFERER8fDzDneXyf4S3IQj9zsrKytfXd9u2beJ34/F4eXl59BCSfnHq1CmhUEi6nymsARMs6If2aMBpbGwUCoWUwg/BBgAfH5/JkyefOHFCzD43btzYvHmzlpZWfx20srKSy+V2vqlRNAMgWOzZs2fGjBnvvPOOvAsC0M18BOIdO3bMycmJrAk4bdo0BwcHCwsLOzu7sLCwv8Ia5UlJSWfOnGlrawsODr548aK8i9MzV1dX0ecXnU2fPr1/B+az2exFixb1Y4ZSMgCCRVBQUHt7uyJ0mO9uPgLx5s+ff+TIEQAYN25caWnpxYsXKyoqEhMTb9y48fbbb2/YsEFhO+H0i5CQEA6HQ1FUcnKyg4ODvIuDJDcAnoawWCxDQ8Mun1TJEpmPgExDAAA+Pj4ODg7x8fEffPBBj58lQx6GDh1Kp1hbW+fn5/v5+X3zzTfDhg378ssvpVdyhPrFALiyUBA9zkcgRpf9/5SVlXft2jVmzJjo6OinT5/2Z1kRkgLFDRbHjx8PDg4OCwv79NNPRadXoyhq9+7dISEhtra2rq6uZM4yMfMpkK0BAQFxcXFz5851cXERk48Y4ucjkGy2BW1t7YULFwoEgoyMDDlWDSFGZPyoluHz4dTUVFtb2zdv3lAUVVdXp6urq6+vTzbFxMQcPHiQoqjW1lYzMzN9fX0+n19bW0tuB1atWnXr1q2Kigo1NbVFixaRj0ycOPHixYsURQkEAgcHBzH5MK9Ia2vr6NGj9+/fT97m5+cPHTo0NTW1y53JoINJkyZ13kSaMwICAuRbtUHc/4V5P4sBRC7nSxGDBZ/PNzAwSEtLo1PmzZtHgsXz58/19PTa2tpI+saNGwHg6NGjFEWR237SlkZRlIODw4QJEyiKamlpUVJS2r59O0nPyckRnw9D2dnZLi4u7e3tdEpra2t3O4sJFmQOWGdnZ/lWDYPFwCKX86WIDZwXLlyora01NzenU9TU1MiL0tJSoVD40Ucf0ZuCgoJIw2Hn+RRIg4KqqurMmTPXrFlTVVUVGxtLxvCIyYcJMh/BqVOnRBsjJOvST4ZUTJw4Ue5Vy8rKknhspeIbfFUjI/FkSRGDBRlaM2TIkM6b7ty5o6mpmZyc3KsMs7OzV6xYkZycnJOTk5GR4eTkJFk+tO7mI5AAWXeDzWbLvWp2dnZr167t7acUX1lZWUJCArm+GDR67GYqDYoYLEiYePLkSef5CzU0NKqrq6urqw0NDenEurq60aNHi8lQRUUlNTX1ww8//Oyzz9zc3K5fvy5ZPoSY+Qh6i6KorKwsVVVVNze3rKws+VbN0NCQ+RyTA0tCQsIgq1pmZqbsD6qIT0PeffddABD9KaA7ZZmbm1MUJbocy6NHj3bt2iUmt+bm5r179wKAr69veXk5RVHFxcUS5EOImY9ATN8qqpuezt99993NmzfDwsLGjRsn96ohJJ4iXllMnz7dycnp4MGD7733nr+//61bty5evFhXV5eenj5nzhxra+u0tLSmpqZ58+a9fv362LFjR48eBbHzKezfvz8kJITFYo0dO1ZbW9vS0tLW1ra7fMQg8xEsX758z549AEBR1I0bN8zMzD744IPCwsIFCxakpKR0eSdJJnERnd/hyZMn33333X/+85/Q0NDIyEgAcHFxkWPVEOqZjBtUGbbicrncgIAAPT09Y2PjTZs2BQcHBwQEkInnX716tWTJkjFjxowePXrZsmXPnz+nKKqwsJB0eVi5cuXLly+///77YcOGAcCmTZv4fL61tfXMmTNjY2ODg4OTk5PJIbrMR4zLly93biZUU1N79eoVRVFFRUUGBga5ubmdP5ibm+vk5ET2d3BwcHZ29vDwcHd3X7duXWVlpeie8qoa8/MyEOHTkP6iRMl2LKC3tzfI6Y4LiTGIz0tGRoaPj4+Mv+fSJpfzpYi3IXIkphVw//79ZL1chP6aFLGBU47quoeR4i+ioKCATBEIAHv37l29erWfn5+zs3NJSQnDHDIzMy0tLYcNG8Zms48fPy66qcspDioqKhISEhT/2geDBZIF0dE98s1EvKSkpEePHpEfhvT09JaWlh07dhw+fNjNzc3R0ZF0txXv4MGDJSUlBw4cyMvLY7FY3t7e9Nic7qY4sLCwYLPZog+wFBMGCyR19fX1HRb1lFcm4uXn5xcVFdHz5eXk5NArCQYGBlIUlZqaKj4HoVD48OHDxMRENpvt5OS0b98+oVB46dIl+N8UB/n5+aGhoQkJCWSug/j4ePJBJyen4cOH79y5U2qV6wfYZoGkSyAQLFq0iB4mK8dMxOPxeIGBgUVFRXSKpaVlh8mNeuwzrqysvGnTJvotWTyVDFbucYqDdevWmZqauru7m5qa9rUy0oFXFqh3srOzP/nkk/Xr17u7u0dERDQ3NwNAenq6lpaWkZERAHC53KioKBaLRebUz8nJuXPnDofDWbFixdatW2/fvr1hwwYzM7OamhpPT89Ro0bZ2NiQ2fGZZwKSzgkgRnJyspqampmZGZ0SHh5OB4uqqioAcHd3F58Ji8USXQI2LS0tMTHx7bffhp6mOAAATU1NKyurzgujKhAZP6odxM/zBzSG52Xbtm3Tpk1raWmhKIrD4UyYMGHGjBlk6K2rq6uhoSG9p7m5uZ2dHXk9a9YsExMT8jo8PHzEiBEsFmvt2rXFxcXZ2dm6uroaGho1NTXMM6F6mhNAFMN+Fvb29t7e3t1tDQ0NtbGxocfy9ojH40VGRurp6f30009d7tBhigMiKipKW1tbzPBlGi4FgBTay5cvIyIiPv74Y7LSj46OzldffXX+/HlyJ6+hoSG6M5lJsLOYmBgPDw9lZeW4uDhHR8f58+cnJSUJBILdu3czzwQAPDw8eDyer69v3+sFAO3t7VevXiV3DZ01NTWdPn06MzOTLFzcIz6fHxkZeenSpfr6+pkzZ+7fv7/zPsePH586dery5ctFE/X09LhcLt1QomgwWCCmysvL+Xy+sbExnTJr1iwAKC4u7lU+ZHkxem0xT09PNTW1mzdv9rY8/bjMX319vVAoHDlyZJdbi4uLw8PDRSsunqamZnx8fH5+/rVr/9/e3Uc1kd0PA78koSC4QQTl5aAbqLosyEYs8rJyVJqCwKEuRVxsXEAUqdC6oIcebaF7fIFCFntAOQgI2rUItgKCQNRdWVCqARY5CLLiG0vXZqFCtkAJwRBgnj9un9mUQBxC3uD3/fxF7kzu3MvAd2bu3JfW5cuXKz5Z4CkOioqKpjWCLFu2DCH06tWruddAGyBYAKq+/fZb9P9X+sbIJ4j5ZMtgMGxtbcmRLzqB485sM8j39vZGRkaqkO369evj4+N7enpkMpl8+mxTHOA7F72d7R2CBaDK3t4eIaT4SsLR0XGeOUskkvlnMh9mZmbGxsZ4QjNFLBZL5blz1q9fb2dnJ79Eq5IpDnAgtra2Vu1YmgbBAlDl5eXFZDIrKyvJFKFQKJFI8JI8DAZDLBaTF2exWExeIWk0Gh53O6O+vr6BgQE8WndOmajxCozXf5rtFonD4aic85MnT+SXLFIyxQFCSCQSMZlMZ2dnlQ+nURAsAFUWFhY8Hu/+/ftffvklTjl79mxkZCQeU+vi4jI0NJSWlvbs2bOUlBSpVPr06dO2tjaEkK2trUgkam1tvXPnDh5cL5VK29vbcSYpKSmRkZHu7u5zyqS2ttbc3LysrExdteNyuQKBgFDoc11VVcViscjSIoTi4uK8vb1nXAViaGgoKirq2rVrOJ8XL17cvXuXx+PhrXiKA5lMlp+fn5+fn5eXFxcXh+eFwwQCwc6dO9XYFqNe0CkLzMHBgwdtbGw+/fTTyspKc3Nza2tr8j8hISHhwYMHPB6Pz+dnZ2d3d3dPTEwIhUJXV9fY2Niamhoul5uamorfdxgaGl66dEkoFDKZTBaLlZSUNNdM6HS6qamp/O39PEVERPB4vKamJtyzgySRSKRS6fj4OJny8uXLxsbGwsJCxdVJGQyGSCSKiYk5c+aMr6+vvb09n8/HPS9aWlpCQ0PHxsZwpxLMyMiIvJ0ZGxsTCAQCgUBdNVI/Lb+qhX4W+kmb5yU6OtrY2Fg7xyLmMp9FS0vLjh07qOzZ0NCQnp4+v3JNl5ycnJGRQXFn6GcBgC65ublxudw3zoU7MjJSXV1NDiFRi5s3b8pkssTERDXmqXYQLIC2icVimUxG6OWI7LCwMGdn56qqKiX7dHR0nDx5kslkquug7e3tw8PDig81+gbaLIBW5ebm3r59e3JyMiYmJjIyUg/XVffz81O+w+bNm9V7RDabzWaz1ZunJkCwAFoVGxur3ht4oDXwGAIAoASCBQCAEggWAABKIFgAACjRQQOnUCi8evWq9o8LlBAKhQihRXleGhsb0aKr2rS1bLVEy53AtL9OPACL0uJfkQzoP7zg+CK7FIP5gzYLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlECwAAJRAsAAAUALBAgBACQQLAAAlBgRB6LoMQMeKi4svXLgwNTWFP/b09CCE7O3t8UcajbZ///49e/borHxAP0CwAKijo4PNZivZob29/b333tNaeYB+gmABEELI0dHx6dOnM25as2bN8+fPtVweoIegzQIghFB4eLihoaFiuqGhYVRUlPbLA/QQ3FkAhBD65ptv1qxZM+Mfw/Pnz9esWaP9IgF9A3cWACGEHBwcNm7caGBgIJ9oYGDg5uYGkQJgECzAf0VERNDpdPkUOp0eERGhq/IAfQOPIeC/+vv7bWxsyBeoCCEajdbb22tlZaXDUgH9AXcW4L9Wrly5detW8uaCTqdv27YNIgUgQbAAPwgPD5e/0wwPD9dhYYC+gccQ8IP//Oc/K1asGB8fRwgZGhr29/cvW7ZM14UC+gLuLMAPmEymv78/g8FgMBiBgYEQKYA8CBbgf3z00UeTk5OTk5MwGARMA48h4H+8fv3a0tKSIAiRSLRkyRJdFwfoE0K7QkNDdV1jABaD0NBQLf/zMrRfSU9Pz8OHD2v/uECJzMxMhBA+Lw8fPjQwMFA+DnUBaWxszMrK+tvf/qbrgqgTPl9apoNgYWdn9+GHH2r/uECJ0tJShBA+LyEhIQghBkMHfxsakpWVtcj+5PD50rLF8wcB1GUxhQmgRvA2BABACQQLAAAlECwAAJRAsAAAUALBAqiosrJy1apVXV1dui6Imn3xxRfV1dX45/Pnz3/88cfh4eEcDqehoYFiDqWlpRs3bly6dCmbzb5+/br8ppKSEjc3NyaT6eHhcePGDZzY1taWlZVF6H33SAgWQEWmpqYrV640NjbW3CH6+vo0l/mMcnNzu7u7f/7znyOErly5Mj4+fvbs2aKiIn9//23btn3++edvzOGzzz5raGj485//XF1dTafTd+3aRU53nJmZefny5fDw8H379nV2dgYFBdXW1iKEXF1d2Wz20aNHNVq1+YNgAVTk6+vb2tpKLi+idoODgx999JGGMp8Rn8+vq6uLjY3FHysqKh4/fox/3r9/P0EQxcXFynOQyWQvXrzIzs5ms9k+Pj6FhYUymay5uRkhJBaLa2pq+Hx+fHx8VlZWbW2tgYFBRkYG/qKPj89bb72Vk5OjscqpAbxRB/pIIpHs3r37m2++0doRR0ZG9u/fX1dXR6Zs3LgRX/lJ0+YoVUSj0Y4fP05+tLCwQAht2rQJIdTc3Jyenk7m4OXl5erq+uLFC3LnI0eOODg4BAQEODg4zLcymgF3FkAVg4ODFy5c8PX1raysRAg9fPjwt7/9rYODw+joaHR0tKWlpbu7O/5Xf/z4cVJSkpOTU29vb3Bw8PLly93d3ZuamhBCV65cYTKZq1atQggNDw+fOnWKTqd7eXkhhCoqKrq6ukQi0YEDB06fPo0Qun///qpVq27evKmhGhUUFBgZGTk5OZEpx44dI4NFZ2cnQiggIEB5JnQ6Xb5LW0lJSXZ29jvvvIMQ4nA4OGqQzMzMWCwW+dHU1NTNzS01NXW+NdEYCBZAFf/617++/vrr2trayclJhJC1tfXDhw97enqOHj165MiR2trajo6OpKQkhFBRUdG5c+eePXt2+vTphISEwsLCnp4eDofT19f3y1/+EocGhJCZmdkf/vAHZ2dn/HHPnj1sNtvS0rKgoCAxMREhNDw8/P333w8ODmqoRmVlZR4eHrNtvXbtmru7O/U+42Kx+OTJk2fOnFm3bt2MO0xOTj569Gjac5aXl1d5eTn+leohCBZAFe++++4HH3xAfrS2tsaXzRMnTjg5OW3YsGHTpk2tra0IobS0tMDAQBqNxuPxtm3bFhISkpubK5FI8vLyEEImJiby2Zqams52xMDAwJGRES6Xq4nqTE1NPXjwAD81KHr9+vWtW7dKS0tpNEr/L6OjoydOnGhubh4cHNy+ffvFixcV97l+/fqGDRv27t0rn2hlZTU8PEw2lOgbCBZARdOGkOCZfslEOzu7kZER/LOJiQmdTidXPAsODjYyMnr06NFcjzhtpQI1GhwclMlk5ubmM26tr68/duzY6tWrKeZmamqakZHB5/NbW1uXL1+u+GQxODiYkpJSVFQ0rREET0326tWruddAGyBYAG1jMBi2trYTExO6LsgPcBia7f6/t7c3MjJShWzXr18fHx/f09Mjk8nk0w8fPpyVlaU4czq+c5FfjUGvQLAAOiCRSBwdHXVdih+YmZkZGxsPDQ3NuJXFYr3xPchs1q9fb2dnJ7+ObE5OTnBw8JYtWxR3/ve//40Qsra2Vu1YmgbBAmhbX1/fwMAAnjONwWCIxWLyki4Wi8nrKo1GE4vF8l/U3CXXwMDg/fff7+3tnXErh8NROecnT57s2LGD/FhSUrJkyZLg4GAyRf7trEgkYjKZZCuvvoFgAVSEu1cODAzgj8PDwwgh8uGiv79fIpGQO0ul0vb2dvxzSkpKZGSku7s7QsjFxWVoaCgtLe3Zs2cpKSlSqfTp06dtbW0IIVtbW5FI1NraeufOHYlEUltba25uXlZWpqHqcLlcgUCg2Oe6qqqKxWKRhUcIxcXFeXt7y3eRIA0NDUVFRV27dg3n8+LFi7t37/J4PLz1xo0b2dnZMpksPz8/Pz8/Ly8vLi7uyZMn5NcFAsHOnTs11zQzT9ApC6iirq4uOzsbIXThwoUf//jHNBoNj6f45JNPjh8/fuvWra+++kosFp84cSI5ORkhZGhoeOnSJaFQyGQyWSwWfquKEEpISHjw4AGPx+Pz+dnZ2d3d3RMTE0Kh0NXVNTY2tqamhsvlpqam4iZSU1NT+ft59YqIiODxeE1NTeTbXEwikUilUryWCvby5cvGxsbCwsL09PRpmTAYDJFIFBMTc+bMGV9fX3t7ez6fjxt9W1paQkNDx8bGcB8TzMjIiLydGRsbEwgEAoFAQxVUAy3P+RkaGqr9iUbBG2n0vERHRxsbG2so8zfCs29S2bOlpWXHjh1U9mxoaEhPT59fuaZLTk7OyMiguLNO/o8WzGMI+R4OAA1xc3PjcrlvnAt3ZGSkurqaHEKiFjdv3pTJZLj7md5aAMEiPz9/69at7777rq4LgpDS0cezuXbtmo+Pj4GBAW5F8/b2dnV19fT0PHr0aHd3t6YLrA/EYrFMJiP0fgg2QigsLMzZ2bmqqkrJPh0dHSdPnmQymeo6aHt7+/DwsOJDjb5ZAMEiOjp6ampKH/rAKhl9rERISMjly5cRQm+//bZAILh3715bW1t2dnZHR8c777yTlJSkt+/V1SI3N/f27duTk5MxMTH37t3TdXHezM/PT/79haLNmzerd2A+m83evXu3GjPUkAXQwEmn0+3s7GZsfNYmcvQx/lhYWPiTn/ykubl57dq1b/wu7sUsv8DXpk2b+Hx+eHj4H//4x6VLl/7ud7/TULF1LjY2Vr137EBXFsCdhZ5QMvr4jWbs0kOj0c6dO7dy5cqUlJSXL1+qqZgAaIr+Bovr16/HxMQcPXr00KFD8jMmEQSRl5cXGxvr4eHh5+eHHwSUDJHGW6Oiong83gcffODr66skHyWUjD5Gqg6gNjMz+/DDDyUSydWrV3VYNQAo0fLbF4qvfIqLiz08PMbGxgiCGBgYsLS0tLa2xpvS0tI+++wzgiAmJiacnJysra1HR0f7+vp+9rOfIYR+/etff/31121tbUZGRrt378ZfWbdu3b179wiCkEgk3t7eSvKhUoWRkZETJ05YWVl9/vnnZCKfz1+yZElxcfGMX8H9iB0dHRU34eaMqKgo3VZtEb/Spv7qdAHRyfnSx2AxOjpqY2NTUlJCpvziF7/AweK7776zsrKanJzE6Z988glC6K9//StBEPixXyQS4U3e3t5r164lCGJ8fNzAwODMmTM4vaKiQnk+yonF4sTExMDAwB/96EcIoQsXLpCbJiYmZvuWkmCBp3XkcDi6rRoEi4VFJ+dLHxs4//73v/f19bm4uJApRkZG+AeBQCCTyX71q1+Rm6Kjo3HDoeIQadwmamhouH379oSEhM7OzvT0dNwtX0k+yuHRxwihzs7OrVu3pqam7tu3D29SrZcu7iW9bt06nVdNKBTip6FFprGxESG0yKomFArt7Oy0fFB9DBa4tzy+dE/T1dVlampaUFAwpwzLy8sPHDhQUFBQUVFx9epVHx8f1fKRh0cfHz9+XCaTzacPMp5Kn81m67xqTU1NYWFhc/3WQrH4qoZH4mmTPjZw4jDx7bffKm4yMTERCoVCoVA+kRzLNBsGg1FcXFxcXMxgMPz9/bu6ulTLZxrF0cdzRRBEWVmZoaGhv7+/zqsGjyELiPYjBdLPYPHee+8hhPA5xshOWS4uLgRByK+w0N3dfe7cOSW5SaXS8+fPI4S4XG5TUxNBEPX19Srko2ja6GMlfauIWTov/ulPf3r06NHRo0fffvttvaoaADPQfkSkcgXz8fGh0+nnzp0bHR396quvbG1tEUIlJSVisRh3bQgJCSkqKsrJyeFwOAMDAwRBHDp0CMm1Av70pz9lMpkEQbx+/dqktrjJAAAHx0lEQVTV1RW3Po6Pj1taWjY2Nk5NTc2Wz2wGBwf37t1bXl4+NTVFEMTz58/9/PzEYjHeevv2bSaTWVpaOuN3//nPfyKEVq9eTab84x//OHTokIGBQXx8PG6MVFIkTVeNgAbOhQbehvxgeHg4KirKyspq9erVx48fj4mJiYqKwnNJf//993v27Fm5cuWKFSsiIiK+++47giBqa2vxrOpxcXH9/f1/+ctfli5dihA6fvz46Ojopk2btm/fnp6eHhMTU1BQgA8xYz5KjIyMBAUFWVhYbNmy5dSpU5cvX8bjHbC6ujobG5vKykrFL1ZWVvr4+ODQ7O3tzeFwAgMDAwICjhw50t7eLr+nrqpGQLBYaHRyvgwI7Q7v2bVrF0KotLRUmwcFb7SIz8vVq1fDwsK0/HeuaTo5X/r4NkSHVqxYMdumixcv4iUwAfi/CYLF/5jrCxEA/u/Qx7chAOiPL774As8YiBA6f/78xx9/HB4ezuFwGhoaKOZQUlLi5ubGZDI9PDxu3LiBE9va2rKyshbWwxEEC6AN8kMBdZvJnOTm5nZ3d+PHzytXroyPj589e7aoqMjf33/btm24q75ymZmZly9fDg8P37dvX2dnZ1BQEJ7O29XVlc1my7/h1n8QLIDGDQ4OTlvUU1eZzAmfz6+rqyMn46ioqCAXFty/fz9BEMXFxcpzEIvFNTU1fD4/Pj4+KyurtrbWwMAADxdACPn4+Lz11ls5OTmaq4J6QZsF0CyJRLJ7925yTL0OM5mTkZGR/fv319XVkSkbN26UX+MDzTJNibzm5ub09HRyNy8vL1dXV/lpnI4cOeLg4BAQEODg4KC+smsK3FmAuSkvL//Nb36TmJgYEBCQnJwslUoRQleuXGEymatWrUIIDQ8Pnzp1ik6n4zn1Kyoqurq6RCLRgQMHTp8+/fjx46SkJCcnp97e3uDg4OXLl7u7u+PZ8alnglSdQIS6goICIyMjJycnMuXYsWNksOjs7EQIBQQEKM+Ew+FMmx7JzMwM95rBTE1N3dzcFBdD1VNa7texiDv/LGgUz0tmZub7778/Pj5OEIRIJFq7du3WrVtxl1Y/Pz87OztyTxcXF09PT/xzUFAQi8XCPx87dmzZsmV0Ov3w4cP19fXl5eWWlpYmJia9vb3UMyHeNIGIPNU6ZXl5ee3atWu2rfHx8e7u7uQ8ABRNTEysWLHi4sWL8omnTp0yMzNTMr/BjGApAKDX+vv7k5OTDx48iMfOWVhY/P73v7979y5+dDcxMZHfGU87qigtLS0wMJBGo/F4vG3btoWEhOTm5kokkry8POqZIIQCAwNHRka4XO7866VoamrqwYMHeOZERa9fv75161ZpaSlex5i669evb9iwYe/evfKJVlZWw8PDZGuIPoNgAahqamoaHR1dvXo1mRIUFIQQqq+vn1M+eHkxcrRucHCwkZHRo0eP5loezS3zNzg4KJPJzM3NZ9xaX19/7Ngx+d8DxTxTUlKKioqmtXQsW7YMIfTq1SuVS6s1ECwAVXjSALzSN0Y+QcwnWwaDYWtrSy6Sqg9wGJpt9Yne3t7IyMi55nn48OGsrCwrK6tp6fj2ZEEsBwHBAlBlb2+PEFJ8JeHo6DjPnCUSyfwzUSMzMzNjY2M8GaIiFov1xvcg0+Tk5AQHB2/ZskVxEw6+1tbWKpRTyyBYAKq8vLyYTGZlZSWZIhQKJRIJntSDwWCIxWLyaiwWi8mrJY1GE4vFs2Xb19c3MDCAZ3OZUyaauxrjteNmu2PicDhzyq2kpGTJkiV40kNM/hWsSCRiMpnOzs6qFVWbIFgAqiwsLHg83v3797/88kuccvbs2cjISDwA38XFZWhoKC0t7dmzZykpKVKp9OnTp21tbQghW1tbkUjU2tp6584diUSCEJJKpe3t7TiTlJSUyMhId3f3OWVSW1trbm5eVlamocpyuVyBQEAodMeuqqpisVhk4RFCcXFx3t7esy2CdePGjezsbJlMlp+fn5+fn5eXFxcXhyeOxAQCwc6dOzXX/qJG0CkLzMHBgwdtbGw+/fTTyspKc3Nza2trHo+HNyUkJDx48IDH4/H5/Ozs7O7u7omJCaFQ6OrqGhsbW1NTw+VyU1NT8fsOQ0PDS5cuCYVCJpPJYrGSkpLmmgmdTjc1NZ3PnIbKRURE8Hi8pqYm3NGDJJFIpFLp+Pg4mfLy5cvGxsbCwkLFxUpbWlpCQ0PHxsZwRxLMyMiIvGcZGxsTCAQCgUBDtVAzLb+qhX4W+kmb5yU6OtrY2Fg7xyLmMflNS0vLjh07qOzZ0NCQnp6uwiGSk5MzMjJU+CL0swBAj7i5uXG53MzMTOW7jYyMVFdXq7Ce682bN2UyWWJioqoF1DYIFkDbxGIxnpFQ1wV5s7CwMGdn56qqKiX7dHR0nDx5kslkzinn9vb24eFhxScXfQZtFkCrcnNzb9++PTk5GRMTExkZ6e3tresSvYGfn5/yHTZv3qxCtmw2m81mq1QinYFgAbQqNjZWhTt2oA/gMQQAQAkECwAAJRAsAACUQLAAAFCigwbOpqYmvEQK0B+4i+GiPC94jehFVrWmpiZPT08tH1TbwWJa51mgJ7T/l6c1dnZ2OllzXKM8PT21/6+k7eULAQALFLRZAAAogWABAKAEggUAgBIIFgAASv4fY4QFJ9GblAEAAAAASUVORK5CYII=\n",
                        "text/plain": [
                            "<IPython.core.display.Image object>"
                        ]
                    },
                    "execution_count": 7,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "# this cell requires graphviz (replace by UML plot)\n",
                "tf.keras.utils.plot_model(perceptron, to_file='model.png', show_shapes=True)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "Next, we create the decision rule which takes as input 5 vectors of the same size $n$ for the states $r$, $\\delta$, $q$, $p$, $w$ and returns two vectors of size $n$ for $\\zeta$ and $h$, respectively. We use different nonlinear transformations for the two decision functions, namely, sigmoid/logistic and exponential ones:\n\n\n$$\n\\begin{equation*}\n\\varphi(s;\\theta)=\\left(\\begin{matrix}\\frac{1}{1+e^{-\\mathcal{N}_1(s;\\theta)}}\\\\ \n\\exp(\\mathcal{N}_2(s;\\theta))\\end{matrix}\\right)\n\\end{equation*}\n$$\n\n\nwhere $\\mathcal{N}_1$ and $\\mathcal{N}_2$ denote first and second component of neural network output. This transformation ensures $\\zeta\\in[0,1]$ and $h>0$"
        },
        {
            "cell_type": "code",
            "execution_count": 8,
            "metadata": {},
            "outputs": [],
            "source": [
                "def dr(r: Vector, δ: Vector, q: Vector, p: Vector, w: Vector)-> Tuple[Vector, Vector]:\n",
                "\n",
                "    # we normalize exogenous state variables by their 2 standard deviations \n",
                "    # so that they are typically between -1 and 1 \n",
                "    r = r/σ_e_r/2\n",
                "    δ = δ/σ_e_δ/2\n",
                "    q = q/σ_e_q/2\n",
                "    p = p/σ_e_p/2\n",
                "    \n",
                "    # we normalze income to be between -1 and 1\n",
                "    w = (w-wmin)/(wmax-wmin)*2.0-1.0\n",
                "\n",
                "    # we prepare input to the perceptron\n",
                "    s = tf.concat([_e[:,None] for _e in [r,δ,q,p,w]], axis=1) # equivalent to np.column_stack\n",
                "\n",
                "    x = perceptron(s) # n x 2 matrix \n",
                "\n",
                "    # consumption share is always in [0,1]\n",
                "    ζ = tf.sigmoid( x[:,0] )\n",
                "    \n",
                "    # expectation of marginal consumption is always positive\n",
                "    h = tf.exp( x[:,1] )\n",
                "    \n",
                "    return (ζ, h)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "Finally, as an illustration, we plot the initial guess of decision rules against $w$. Note that the coefficients of the perceptron are initialized with random values, so that each run will provide a different plot. Here, we are using TensorFlow in an eager mode, i.e., calculations are returned immediately, so that the library essentially behaves in the same way as numpy, and is in fact mostly compatible with it."
        },
        {
            "cell_type": "code",
            "execution_count": 9,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "WARNING:tensorflow:Layer dense is casting an input tensor from dtype float64 to the layer's dtype of float32, which is new behavior in TensorFlow 2.  The layer has dtype float32 because it's dtype defaults to floatx.\n",
                        "\n",
                        "If you intended to run this layer in float32, you can safely ignore this warning. If in doubt, this warning is likely only an issue if you are porting a TensorFlow 1.X model to TensorFlow 2.\n",
                        "\n",
                        "To change all layers to have dtype float64 by default, call `tf.keras.backend.set_floatx('float64')`. To change just this layer, pass dtype='float64' to the layer constructor. If you are the author of this layer, you can disable autocasting by passing autocast=False to the base Layer constructor.\n",
                        "\n"
                    ]
                }
            ],
            "source": [
                "wvec = np.linspace(wmin, wmax, 100)\n",
                "# r,p,q,δ are zero-mean\n",
                "ζvec, hvec = dr(wvec*0, wvec*0, wvec*0, wvec*0, wvec)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 10,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEYCAYAAAC+xZqSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5dnH8e+dENawZwHDJrIGBCVhU6AEUFARbcUCFhQV0LpAXVqVV8GlKmgtdasbtUorjYqKiCD7IliQHUnCKkjYsgdIWLLM/f4xI40hgRAycybJ/bmuuZyZ88ycXw7O3HPOec7ziKpijDGmcgtwOoAxxhjnWTEwxhhjxcAYY4wVA2OMMVgxMMYYgxUDY4wxWDEwlYyIzBeRO86x/G0ReaqE77VcRMaUXTpjnGPFwJR7IrJPRAaUpK2qXqeqH3peN1pEVhVafq+qPldGuVqLSKyIpIjIMRHZJSKvi0iTsnh/Y8qSFQNjvEBEWgFrgUPAlapaB7ga2AP0cjKbMUWxYmAqlJ9/7YvIX0QkQ0T2ish1BZYvF5ExItIeeBvoKSJZIpLpWf6BiPzZc7++iMz1/LLP8Nwv6a/6p4HVqvqwqh4AUNVkVf2bqsYWzFoov3oKCSJSzfN37BeRJM8hrBqeZSGePJkiki4i34pIgGfZYyJyUESOi8gOEel/EZvUVBJWDExF1B3YAYQALwH/EBEp2EBVE4B7gf+qarCq1ivifQKAfwLNgWbASeCNEmYYAHxWuvhnTAXaAFcArYAIYJJn2SPAASAUCAcmAioibYEHgK6qWhsYCOy7yBymErBiYCqin1T1PVXNBz4EGuP+wrwgqpqmqp+p6glVPQ48D/yqhC8PAY78/EBEHvD8is8SkffO92JP8RoLPKSq6Z71vwAM9zTJxf13NVfVXFX9Vt0DjeUD1YBIEQlS1X2quqekf7OpvKwYmIrozJewqp7w3A2+0DcRkZoi8o6I/CQix4CVQD0RCSzBy9Nwf1n/nOMNz97H34CgErw+FKgJbPAUkUzgG8/zAC8Du4GFIvKjiDzuWc9u4A+4D1Mle05gX1KSv9dUblYMTGV2viF7HwHaAt09J4D7eJ6X4l9yxhLgN+dpk437C9/9piKNCixLxX1YqoOq1vPc6qpqMICqHlfVR1S1JXAj8PDP5wZUdaaq9sJ9eEtxH24y5pysGJjKLAloIiJVi1leG/cXcqaINAAmX8B7Pw30FpG/ikgEuE/6Au0LtNkCdBCRK0Skuuc1AKiqC3gPmCYiYZ7XR4jIQM/9wSLSynM46Rjuw0P5ItJWRPqJSDXglCd//gXkNpWUFQNTmS0F4oAjIpJaxPK/ATVw/0pfg/swTYmo6k6gB9AE2CIix4HVuLuaPlWgzbPAYmAXsKrQ2zyG+1DQGs9hqsW491QAWnseZwH/Bf6uqstxny+Y4sl8BAjDfXLZmHMSm9zGGGOM7RkYY4yxYmCMMcaKgTHGGKwYGGOMwYqBMcYYoIrTAUorJCREW7RoUeSy7OxsatWq5dtAJWTZSseylZ4/57NspXMx2TZs2JCqqqFnLVDVcnmLiorS4ixbtqzYZU6zbKVj2UrPn/NZttK5mGzAei3iO9UOExljjLFiYIwxxoqBMcYYrBgYY4zBh8VARAJFZJOIzC1iWTUR+VhEdovIWhFp4atcxhhjfLtnMAFIKGbZ3UCGqrYCpmHjrxtjjE/5pBh4JhG/AZheTJObcE9PCDAL6F94zlpjjKnsXC4XJ0+e9Mp7+2QIaxGZBbyIe7KQR1V1cKHl24BBqnrA83gP7tmlUgu1GweMAwgPD4+KjY0tcn1ZWVkEB1/wLIc+YdlKx7KVnj/ns2wll5iYyCuvvELjxo25//77S50tJiZmg6pGn7WgqIsPyvIGDMY98QZAX2BuEW3igCYFHu8BGp7rfe2is7Jn2UrHn7Op+nc+y3Z+LpdLX3jhBa1WrZrWq1dPp0+fXm4vOrsaGCIi+4BYoJ+I/LtQmwNAUwARqQLUBdJ9kM0YY/yaiLB582YGDx5MfHw8d999t1fW4/VioKpPqGoTVW0BDAeWqurIQs3mAHd47g/1tLEp2IwxlVJ2djZ//OMfSUhw97mZMWMGs2bNonHjxl5bp2MD1YnIs7h3V+YA/wD+JSK7ce8RDHcqlzHGOGnRokXcc8897N27l4iICNq3b0+1atW8vl6fFgN1T9i93HN/UoHnTwG3+jKLMcb4k7S0NB555BE+/PBD2rRpw4oVK+jTp4/P1m9XIBtjjB947bXX+Oijj/i///s/tmzZ4tNCAOV4PgNjjCnvEhMTSUlJoUuXLjz22GMMHTqUyy+/3JEstmdgjDE+5nK5ePPNN4mMjOTOO+9EValZs6ZjhQCsGBhjjE/Fx8fTu3dvHnjgAa666ipmz56NPwy4YIeJjDHGR9atW0evXr0IDg5mxowZjBw50i8KAdiegTHGeF1mZiYAUVFRTJw4kYSEBEaNGuU3hQCsGBhjjNccP36c8ePH07p1a5KSkggICGDy5MmEhYU5He0sVgyMMcYL5s+fT8eOHXnjjTcYMWIENWvWdDrSOdk5A2OMKUO5ubmMHj2amTNnEhkZyerVq+nZs6fTsc7L9gyMMaYMBQUFERgYyNNPP83GjRvLRSEAKwbGGHPR9u7dy4033kh8fDwAH374IZMnT/bJmEJlxYqBMcaUUn5+PtOmTaNjx44sX76c7du3A/hVL6GSsmJgjDGlsHXrVnr27MnDDz9M3759iYuL4ze/+Y3TsUrNTiAbY0wpzJw5k3379vGf//yHYcOGlcu9gYJsz8AYY0ro22+/ZdWqVQBMmjSJhIQEhg8fXu4LAVgxMMaY8zp69Ci///3v6dOnD08//TQANWvWpGHDhs4GK0NeLwYiUl1EvheRLSISJyLPFNFmtIikiMhmz22Mt3MZY0xJfPnll0RGRvLuu+/y0EMP8eWXXzodySt8cc7gNNBPVbNEJAhYJSLzVXVNoXYfq+oDPshjjDElMn/+fG6++WY6derE7Nmz6dq1q9ORvMbrxcAzsX2W52GQ52aT3Rtj/JKqsmfPHgAGDhzI+++/z8iRIwkKCnI4mXeJ+7vayysRCQQ2AK2AN1X1sULLRwMvAinATuAhVU0s4n3GAeMAwsPDo2JjY4tcX1ZWFsHBwWX5J5QZy1Y6lq30/Dmfv2U7ePAgr7zyCrt37+att94iIiLC6UhFupjtFhMTs0FVo89aoKo+uwH1gGVAx0LPNwSqee7fCyw933tFRUVpcZYtW1bsMqdZttKxbKXnz/n8JVtubq5OnTpVq1evrnXq1NF33nlHlyxZ4nSsYl3MdgPWaxHfqT7tTaSqmcByYFCh59NU9bTn4XtAlC9zGWMqr6ysLLp168Zjjz3GwIEDiY+PZ9y4cQQEVK7Olr7oTRQqIvU892sAA4Dthdo0LvBwCJDg7VzGmMrN5XIBEBwcTO/evZk1axZffPGF3x4a8jZflL7GwDIR2QqsAxap6lwReVZEhnjajPd0O90CjAdG+yCXMaaSWrp0KZdffvmZgeVeffVVbrnllgpx8Vhp+aI30VbgyiKen1Tg/hPAE97OYoyp3DIyMnj00Ud5//33adWqFVlZWed/USVRuQ6KGWMqrc8//5zIyEg+/PBDHnvsMbZu3Uq3bt2cjuU3bKA6Y0ylsHbtWi655BLmzZvHlVeedbCi0rM9A2NMheRyuXj77bdZunQpAM8++yxr1661QlAMKwbGmApnx44d9O3bl9///vfMnDkTgGrVqlGlih0MKY4VA2NMhZGTk8Pzzz9P586d2bZtG//85z957733nI5VLlgxMMZUGJ988glPPvkkN910EwkJCYwePbpSdxe9ELbPZIwp17Kzs/nhhx/o0aMHt912GxEREcTExDgdq9yxPQNjTLm1cOFCOnbsyA033EBWVhYBAQFWCErJioExptxJTU3l9ttvZ+DAgVSvXp0vv/zSr0Y/LY/sMJExplxJSkqiY8eOZGZm8tRTTzFx4kSqV6/udKxyz4qBMaZcyM7OplatWoSHhzN+/HhuvvlmLr/8cqdjVRh2mMgY49fy8/N5/fXXadasGXFxcQA89dRTVgjKmBUDY4zfiouLo1evXowfP56uXbtSu3ZtpyNVWFYMjDF+6c9//jNXXnklu3bt4t///jfz58+nWbNmTseqsKwYGGP8UnZ2Nr/97W9JSEjgd7/7nV085mV2AtkY4xeOHz/OxIkTuemmmxgwYADPP/98pZt60km+mPayuoh8LyJbPLOZPVNEm2oi8rGI7BaRtSLSwtu5jDH+4+uvvyYyMpI333yT9evXA1gh8DFfbO3TQD9V7QxcAQwSkR6F2twNZKhqK2AaMNUHuYwxDsvIyGDEiBEMHjyYunXr8t133/H44487HatS8noxULef55YL8ty0ULObgA8992cB/cUOEBpT4a1atYrPP/+cZ555ho0bN9KjR+HficZXRLXw97IXViISCGwAWgFvqupjhZZvAwap6gHP4z1Ad1VNLdRuHDAOIDw8PCo2NrbI9WVlZfntpemWrXQsW+n5W77Dhw+TmJhIt27dOHbsGMePHyciIsLpWGfxt+1W0MVki4mJ2aCq0WctUFWf3YB6wDKgY6Hn44AmBR7vARqe672ioqK0OMuWLSt2mdMsW+lYttLzl3y5ubn6yiuvaM2aNbVZs2aak5PjN9mKUlGzAeu1iO9Un56hUdVMYDkwqNCiA0BTABGpAtQF0n2ZzRjjPVu2bKFnz5488sgj9OvXj1WrVhEUFOR0LFOAL3oThYpIPc/9GsAAYHuhZnOAOzz3hwJLPRXMGFPO7dq1i+joaPbv38/HH3/MnDlzaNq0qdOxTCG+uM6gMfCh57xBAPCJqs4VkWdx767MAf4B/EtEduPeIxjug1zGGC9KTEykadOmtG7dmjfffJNbbrmFhg0bOh3LFMPrxUBVtwJXFvH8pAL3TwG3ejuLMcb7MjMzeeyxx/jggw/YuHEjHTp0YNy4cU7HMudhVyAbY8rMF198wf33309SUhIPPfQQLVq0cDqSKSErBsaYi6aqjBgxgo8//pjOnTszZ84coqPP7r1o/Jdd722MKbWf+3mICO3bt+fFF19k3bp1VgjKISsGxphS2bVrF/369WPRokUATJ48mccff9y6jJZTVgyMMRckNzeXqVOn0qlTJzZt2kRGRobTkUwZsHMGxpgS27BhA2PGjGHz5s38+te/5o033uCSSy5xOpYpA1YMjDEltm7dOpKSkvjss8/4zW9+43QcU4asGBhjzmnJkiVkZGQwdOhQxo0bx4gRI6hbt67TsUwZs3MGxpgipaenc9dddzFgwABeeukl92BmAQFWCCooKwbGmF9QVT755BPat2/PjBkzeOKJJ1ixYoXNQVzB2WEiY8wvrF+/nmHDhhEVFcWCBQu44oornI5kfMD2DIwxuFwu1q1bB0DXrl2ZO3cua9assUJQiVgxMKaS2759O3369OGqq65i9+7dANxwww1UqWIHDioTKwbGVFI5OTn8+c9/pnPnzsTHxzN9+nQuu+wyp2MZh1jpN6YSys3NpUePHmzatIlhw4bx6quvEh4e7nQs4yDbMzCmEjl9+jQAQUFBjBw5kjlz5hAbG2uFwPhk2sumIrJMRBJEJE5EJhTRpq+IHBWRzZ7bpKLeyxhTet988w1t27Zl4cKFADz88MPceOONDqcy/sIXh4nygEdUdaOI1AY2iMgiVY0v1O5bVR3sgzzGVCpHjx5l5MiRfPTRR7Rr1446deo4Hcn4IV9Me3kYOOy5f1xEEoAIoHAxMMaUsU8//ZSxY8dy4sQJJk2axMSJE6lWrZrTsYwf8ukJZBFpgXs+5LVFLO4pIluAQ8Cjqhrnw2jGVEgZGRlERETw8ccf07FjR6fjGD8mP89U5PUViQQDK4DnVfXzQsvqAC5VzRKR64FXVbV1Ee8xDhgHEB4eHhUbG1vkurKysggODi7rP6FMWLbSsWwlk5+fz+zZswkODmbgwIG4XC6OHz/ut+MJ+dO2K6yiZouJidmgqmdPRaeqXr8BQcAC4OEStt8HhJyrTVRUlBZn2bJlxS5zmmUrHct2fj/88IN2795dAR0xYsSZ5/0lX1EsW+lcTDZgvRbxneqL3kQC/ANIUNW/FtOmkacdItINdy+nNG9nM6YiOH36NJMmTaJLly7s2bOHmTNn8tFHHzkdy5QzvjhncDUwCvhBRDZ7npsINANQ1beBocDvRSQPOAkM91QwY8x5rF69mueee46RI0cybdo0QkJCnI5kyiFf9CZaBZxz7FtVfQN4w9tZjKkojh07xsqVKxk8eDD9+vVjy5YtdOrUyelYphyzK5CNKWe++uorIiMjGTp0KElJSQBWCMxFs2JgTDmRlJTEsGHDGDJkCPXr12fFihU2jIQpMzZQnTHlQFZWFp06dSIzM5PnnnuOP/3pT1StWtXpWKYCsWJgjB9LSUkhNDSU4OBgpkyZQs+ePWnXrp3TsUwFZIeJjPFDeXl5/OUvf6F58+YsWLAAgDvvvNMKgfEa2zMwxs9s3ryZu+++m40bNzJkyBAbRsL4hO0ZGONHXnjhBaKjozl48CCffvops2fPJiIiwulYphKwYmCMHwkNDeX2228nPj6eoUOH4rkw3xivs2JgjIMyMzMZO3Ys7733HgBjx47l/fffp0GDBg4nM5WNnTMwxiGff/45999/P8nJyTRr1szpOKaSs2JgjI8dOnSIBx54gC+++IIrrriCr7/+mi5dujgdy1RydpjIGB/74YcfmD9/PlOmTOH777+3QmD8gu0ZGOMDO3fu5LvvvmP06NEMHDiQffv22VASxq/YnoExXpSbm8uLL75Ip06d+NOf/sTx48cBrBAYv2PFwBgvWb9+PdHR0UycOJHBgwezZcsWateu7XQsY4pkh4mM8YLk5GR69+5NgwYN+OKLL7j55pudjmTMOV3QnoGIPFzgftsSvqapiCwTkQQRiRORCUW0ERF5TUR2i8hWEbEzaqZc2rZtGwBhYWF8/PHHxMfHWyEw5UKJioGI1BORfwJDReQ+EekFPF7CdeQBj6hqe6AHcL+IRBZqcx3Q2nMbB7xVwvc2xi+kpaUxevRoLr/8cpYsWQLAkCFDqFu3rsPJjCmZ8x4mEpFo3JPZ3ykiA4FUoBPweUlWoKqHgcOe+8dFJAGIAOILNLsJmOGZ93iNp/g09rzWGL+lqixdupRhw4aRnp7OxIkTufrqq52OZcwFK8k5gxlAFwBVXSAiIUBjVZ17oSsTkRbAlcDaQosigMQCjw94nrNiYPzaHXfcwb/+9S+6du3KokWLbPpJU26J+8f4ORqIbFTVLud77rwrEgkGVgDPq+rnhZZ9Dbyoqqs8j5cAf1LVDYXajcN9GInw8PCo2NjYIteVlZVFcHDwhcTzGctWOv6UzeVyARAQEMCCBQtISUlhxIgRBAYGOpysaP607QqzbKVzMdliYmI2qGr0WQtU9Zw3YBZwXaHntp3vdYXaBwELgIeLWf4OMKLA4x249z6Kfc+oqCgtzrJly4pd5jTLVjr+ki0uLk6vuuoqfeedd8485y/ZiuPP+Sxb6VxMNmC9FvGdWpITyA8CL4rITBEZLyLvAHtKWoXEPQbvP3Cfd/hrMc3mALd7ehX1AI6qnS8wfuT06dM888wzXHHFFWzfvp06deo4HcmYMnXecwaqelhEooCbgSuArcDD537VL1wNjAJ+EJHNnucmAs087/82MA+4HtgNnADuvID3N8ar1q1bx5133klcXBy33XYb06ZNIywszOlYxpSpEl10pqr5wGee2wVR93mAc87Q4dl1uf9C39sYX8jIyODYsWPMnTuXG264wek4xniFXYFsTBHmzZvHzp07+cMf/sC1117Lrl27qFatmtOxjPEaG5vImAJSUlK47bbbuOGGG/jggw/Izc0FsEJgKjwrBsbg7lX3r3/9i/bt2zNr1iyefvpp1q5dS1BQkNPRjPmFI0dPcTrv3JcElIYdJjIG+PHHH7nrrrvo2rUr06dPJzKy8Igpxvhebr6L+EPH2PBTBhv3Z7BpfyYHM0/y4JXVGFjG67JiYCqt/Px8FixYwPXXX89ll13GqlWr6Nq1KwEBtsNsnJGencPGnzLYsD+DDT9lsPVAJqdy3Rc5XlK3Ol2a1+fuXpdS+9jeMl+3FQNTKW3dupUxY8awbt061qxZQ/fu3enevbvTsUwloqrsScliw08ZrN/nLgA/pmQDUCVA6HBJHW7r1pyo5vXp0rwejevWOPPa5ct/KvM8VgxMpXLq1Cn+/Oc/M3XqVOrXr8/MmTPp1q2b07FMJXAqN5+tB46y/qd0Nni+/DNPuDso1K8ZRFTz+tzSpQnRzevTqUk9alT17fAmVgxMpaGq9O3bl7Vr1zJq1CimTZtGw4YNnY5lKqj07BzW70tn/U8ZrN+Xzg8Hj5Kb7z7x2zK0FtdGhhPdvAFdmtfnstBauAdrcI4VA1PhHTt2jODgYAICAnjkkUeoU6cOAweW9ek3U5mpKvvTT/D93nQ2/JTBun3p7PEc8qkaGMDlTepy19WXEt2iAVHN69OgVlWHE5/NioGp0L788kvuu+8+Jk2axD333MOtt97qdCRTAeS7lITDx1i3L531+9xf/snHTwNQt0YQ0c3rMzSqKdEt6nN5RF2qB/nniLYFWTEwFdKRI0cYP348n376KZ06daJLF5tJ1ZTeyZx8NiVmnPni37Q/k6zTeQBE1KvBVZc1pOulDejaogGtQoMJCHD2kE9pWDEwFc6sWbMYO3YsJ0+e5Pnnn+ePf/yjXTxmLsiJnDy2puTx33kJrN2bzraDR8lzKSLQNrw2N195CV1bNCC6RQMi6tU4/xuWA1YMTIVTp04dOnXqxLvvvkvbtm2djmPKgdx8F1sSM1m1O5XvdqexKTGD3HylauA+OjWpy9g+Lenaoj5RzRpQt2bF/GFhxcCUe3l5eUybNo2TJ08yadIkrr32Wq655hrHe2cY/+VyKfGHj/HdnlS+25PGur3pZOfkIwIdL6nL3b1aEpx9gLuHxPi8i6dTrBiYcm3Tpk2MGTOGjRs3cuutt6KqiIgVAnOWAxknWLUrlW93p/Ld7lQyPH38W4UF85suTbi6VUN6tGxIvZrunj7Llx+pNIUArBiYcurEiRM888wzvPLKK4SGhvLpp59yyy23WBEwZxw/lcuaH9NZtSuFb3el8mOqu6tneJ1q9GsXTu/WIVx1WUPC6lR3OKl/8HoxEJH3gcFAsqp2LGJ5X+BL4OfBNj5X1We9ncuUb3v37uVvf/sbd955Jy+99BL169d3OpJxWF6+i60Hj7JqVyqrdqWycX8GeS6lRlAgPVo2YGSP5vRuHUKrsGD70VAEX+wZfAC8Acw4R5tvVXWwD7KYciwjI4PPPvuMMWPG0KFDB3bt2kWzZs2cjmUctD/tBN/uTuHbnal8tyeVY6fyzhz3H9enJb1bh9KleT2qVak8h3tKy+vFQFVXikgLb6/HVFyqyooVKxg+fDipqan06dOHNm3aWCGohI6dyuW/e9L41nPo56e0E4B7RM9BHRvRu3UoV7cK8csrfP2dv5wz6CkiW4BDwKOqGud0IOMfDh48yP3338+XX35Jly5dmD9/Pm3atHE6lvGRfJey9UAmK3em8u2uFDYlZpLvUmpWDaRny4aMvqoFvVuH+sXYPuWduOei9/JK3HsGc4s5Z1AHcKlqlohcD7yqqq2LeZ9xwDiA8PDwqNjY2CLXl5WVRXBwcBmlL1uWreTy8/MZOXIkGRkZ3Hbbbfzud78jMND/dvf9bbsV5s/5isqWdtLFD6n5bEvNJz4tnxN5IECLOgF0CAmkQ8NAWtcPoIqXr/Itb9utpGJiYjaoanTh5x0vBkW03QdEq2rqudpFR0fr+vXri1y2fPly+vbte8E5fcGynd+PP/5IixYtCAgIOLMnkJiY6BfZiuIv2604/pxv+fLl9OzVm+/3prN8RwordqawOzkLgEZ1qtOnTYhjh378fbuVNpuIFFkMHD9MJCKNgCRVVRHphnte5jSHYxkH5Obm8tJLL/Hcc8/x2muvMW7cOK677joAEhMTHU5nylJi+gmW70hm1oZT/H7JIk7m5lM1MIDuLRswvGtT+rQJpbX1+vEpX3Qt/Q/QFwgRkQPAZCAIQFXfBoYCvxeRPOAkMFx9sbti/Mq6deu4++67+eGHH7j11lsZMmSI05FMGTqdl8/3e9NZtj2F5TuTz8zoFVpDuDW6GX3bhtKjZUNqVnX892ml5YveRCPOs/wN3F1PTSX10ksv8cQTT9CoUSNmz57NTTfd5HQkUwYS00+wfGcKK3Yks3p3mvvXf5UAul/agN91b07ftqHs37aOmJjzHj02PmBl2Djm56EjunTpwrhx45gyZQp169Z1OpYppZw8F+v2pbN8RzLLdvzv2H/TBjUYGtWEmHZn//pPtMNAfsOKgfG5tLQ0HnroISIiInjxxRcZMGAAAwYMcDqWKYUjR0+xfEcyS7cns3p3Ktk5vzz2H9MujJYh1u2zPLBiYHxGVYmNjWXChAlkZGTw1FNPOR3JXKB8l7JpfwbLdiSzdHsKCYePAe6Lvm66MoKYtmFcdVlDalWzr5byxv7FjE8kJiZy7733Mm/ePLp37857773H5Zdf7nQsUwIZ2Tms2JnC0u3JrNyVQuaJXAIDhOjm9Xn8unb0bRtK2/Da9uu/nLNiYHzi2LFj/Pe//2XatGk8+OCDfnnxmHFTdY/1v2y7+/DP5sRMXAoNa1Wlf7twYtqF0rt1KHVrVMxJXiorKwbGa+Li4vjiiy948skn6dChA/v37/fbKzoruxM5eazencbS7cks257MkWOnAOjUpC4P9mtNTLswOkXULZdz+5qSsWJgytzp06d58cUXeeGFF6hTpw5jxoyhUaNGVgj8zE9p2e5f/ztSWPNjGjl5LoKrVaF36xBi2oXRt20oYbVtrP/KwoqBKVPfffcdY8aMISEhgd/97ndMmzaN0NBQp2MZ/tf1c+n2ZJbt+N+FXy1DajGye3P6tw+ja4sGVK0S4HBS4wQrBqbMZGdnM2TIEGrVqsW8efPODCVhnJN8/BQrD+QS+9xD4i8AABdQSURBVK8NfLsrxd31s0oAPVo2ZFSP5sS0DaNFSC2nYxo/YMXAXLSVK1fSq1cvatWqxVdffUXHjh2pXbu207EqJVUl7tAxlm5PZklCElsOHAWgUZ3MM10/r25lwz6Ys9n/EabUkpOTmTBhArGxscyYMYNRo0bRs2dPp2NVOqdy8/luTyqLE5JZmuA++SsCnZvU49Fr21Anaz+jbuxnXT/NOVkxMBdMVZkxYwYPP/wwx48f55lnnmHYsGFOx6pUjhw9debX/+o9qZzKdVGraiC9W4fSr30YMW3DCK1dDYDlyw9aITDnZcXAXLB7772Xd999l6uuuor33nuPyMhIpyNVeC6XsvXgUZYmJLFkezJxh9xX/japX4PhXZsR0y6MHi0b2Fy/ptSsGJgSycvLIz8/n2rVqjFs2DA6d+7MvffeS0CA9TzxlqzTeazalcKSBHfvn9SsHAIEops34PHr2tGvXZiN+W/KjBUDc15btmxhzJgx9O/fnylTptCvXz/69evndKwKKTH9BEs8v/7X/phOTr6L2tWr0LdtGP3bhfGrNqHUt8nejRdYMTDFOnXqFM899xwvvfQS9evXp0uXLk5HqnDyXcrG/RksSXAf/9/lGfa5ZWgtRl/dgn7twohqXp+gQNsDM97li5nO3gcGA8lFzYEs7n3cV4HrgRPAaFXd6O1c5tzWrVvHyJEj2blzJ6NHj+aVV16hQYMGTseqEI6dymXlzv8d/sk8kUuVAKHbpQ0Y3q0Z/dtZ33/je77YM/gA90xmM4pZfh3Q2nPrDrzl+a9xUM2aNQkICGDRokU210AZ2J92gsUJSSzdnszavWnk5iv1awYR0zaM/u3D6NMmlDrVbeA34xxfTHu5UkRanKPJTcAMz7zHa0Sknog0VtXD3s5mfmnVqlXMnj2bv/3tb3To0IG4uDg7QVxKP4/7v7jQ4Z9WYcHcdfWlDIgMp0uz+gTawG/GT/jDOYMIILHA4wOe56wY+Mjhw4d58MEH+eyzz+jcuTPHjx+ndu3aVgguUNbpPFbuTGFxQhLLd6SQnp3zi8M/A9qH0byhHf4x/kncP8i9vBL3nsHcYs4ZfA28qKqrPI+XAH9S1Q1FtB0HjAMIDw+Pio2NLXJ9WVlZfjtCpj9lU1XmzZvHW2+9RU5ODiNGjGDUqFFUqeIPvxF+yZ+2W0GpJ12sTTxBwtEqJKTnk69QKwg6hQRyRVgVOoYEUivI2V///rrtwLKV1sVki4mJ2aCq0YWf94dP/QGgaYHHTYBDRTVU1XeBdwGio6O1b9++Rb7h8uXLKW6Z0/wpW2pqKrfccgtRUVG89957HDp0yG+yFeYv283lUrYcyGRJQjKLE5LYfiQbEFqGVOeuXmH0bx9OdPP6VPGj3j/+su2KYtlKxxvZ/KEYzAEeEJFY3CeOj9r5Au/Jzc3lo48+4vbbbyckJIQ1a9bQqlUrAgICOHSoyBpc6Z3MyWfV7lQWx7v7/6dmnT4z7eP/Xd+e2sf3MfyGvk7HNOai+KJr6X+AvkCIiBwAJgNBAKr6NjAPd7fS3bi7lt7p7UyV1YYNGxgzZgybN28mLCyM66+/njZt2jgdyy+lZ+ewJCGJhfFJfLsrhVO5LmpXq0KftqEM8Iz9U6+m++Kv5cv3O5zWmIvni95EI86zXIH7vZ2jMjtx4gSTJ0/mr3/9K+Hh4Xz++edcf/31TsfyO/vTTrAw/ggL45NYvy8dl0LjutX5bXRTrokMp/ulDW3iF1Nh+cNhIuNlv/71r1m4cCFjx47lpZdeol69ek5H8gs/j/2/MM5dALYfOQ5Au0a1eSCmFddENqJjRB0b+8dUClYMKqj09HRq1KhBjRo1eOqpp5g4cSK/+tWvnI7luLx8F9/vTWdhfBIL445w6Ogp9+BvLRrw5A3tuTayEc0a1nQ6pjE+Z8WgglFVZs2axYMPPsgdd9zB1KlT6dWrl9OxHHUqN5+VO1NYEJfEku1JZJ7IpVqVAPq0CeWha9rQv304DWzwN1PJWTGoQA4cOMD999/PnDlziIqK4rbbbnM6kmOOnsxl2fZkFsQdYfmOFE7m5lOnehX6tw9nYIdG9GkTYlM/GlOAfRoqiNmzZ3P77beTl5fHX/7yFyZMmOCXF49505Gjp1jkOQH83z1p5LmUsNrVuCUqgoEdGtGjZUMb/dOYYlSub4sKSFUREdq0aUPv3r15/fXXadmypdOxfEJV2ZF0nMXxSSyK/9/k7y1DajGmd0uu7RDOFU3qEWDj/xhzXlYMyqmcnBymTp3K7t27+fDDD4mMjOTrr792OpbXncrNZ82PaSzbnsyS7ckcyDgJQOem9fjjwLYM7BDOZaE2+5cxF8qKQTm0du1axowZw7Zt2xg2bBg5OTlUrVpxT4AezDzJ0v25/PvDdazencbJ3HyqBwVw9WUh3B/Tiv7twgirU93pmMaUa1YMypGsrCyefPJJXnvtNSIiIpgzZw433nij07HKXF6+iw0/ZbB0RzLLt6ewI8nd/79pg+PcGt2EmHZh9GzZkOpBNvm7MWXFikE5kp2dzb///W/uu+8+XnjhBerUqeN0pDKTlnWaFTtTWLo9mZU7Uzh2Ko8qAULXFg34v+vbU+vYPkbcEGOHf4zxEisGfi41NZW///3vPPnkk4SHh7Nz584KMf2ky6VsO3SUZdtTWLYjmS0HMlGFkOBqDOzQiH7twujVOoTantm/li/fb4XAGC+yYuCnVJWZM2fyhz/8gaNHjzJo0CC6detWrgtBRnYOK3elsGJHCit3pZCalYMIdG5Sjwn9W9OvXRgdL6lrvX+McYAVAz/0008/ce+99/LNN9/QvXt3pk+fTseOZ80L5PdcLmXrwaMs257Mip0pZ37916sZxK/ahNK3bSh9WofSMLia01GNqfSsGPgZVeXmm29m165dvPbaa9x3330EBpafE6VHT+Sycpf70M+KHSmkZbt//XdqUo/x/Vrzq7ahdG5Sz+b+NcbPWDHwE3FxcVx66aXUrFmT6dOnExISQvPmzZ2OdV4/X/i1dLu758+G/Rnku5T6Z379h9GnTaiN/WOMn7Ni4LDTp0/z/PPPM2XKFB577DGee+45oqKinI51Tqfz8lnzYzqL45NYuj2Zg5nuC786XFKH+/peRky7MPv1b0w545NiICKDgFeBQGC6qk4ptHw08DJw0PPUG6o63RfZnLR69WrGjBnD9u3bGTVqFBMmTHA6UrEysnNYtsM97++KHSlk5+RTIyiQXq1DeLBfK2LahRFuF34ZU275YtrLQOBN4BrgALBOROaoanyhph+r6gPezuMvXn/9dcaPH0/z5s355ptvGDhwoNORzrInJYslCUksTkhmw0/uwz9htasx5IpLuCYynKsuC7ELv4ypIHyxZ9AN2K2qPwJ4Jr6/CShcDCqF3NxcAAYMGMBDDz3Es88+S3BwsMOp3FyqfL83ncUJSSyOT+LH1GwA2jd2H/4Z0D6cyyOs66cxFZEvikEEkFjg8QGgexHtbhGRPsBO4CFVTSyiTbmVlJTEhAkTSEpK4pprrqF9+/b89a9/dToW2afz+HZXCovik1nwwwmycv9LUKDQo2VD7riqBf3bh9Gkvs38ZUxFJ+756L24ApFbgYGqOsbzeBTQTVUfLNCmIZClqqdF5F7gt6rar4j3GgeMAwgPD4+KjY0tcp1ZWVl+82tbVVmwYAF///vfOXXqFMOGDePOO+8kIMC5cfUzT7nYnJLPxuR84tPyyXNBzSrQob7S9ZLqXB4aSI0q/vXr35/+TQvz52zg3/ksW+lcTLaYmJgNqhp91gJV9eoN6AksKPD4CeCJc7QPBI6e732joqK0OMuWLSt2mS8lJibqgAEDFNBevXppQkKCI9lcLpduP3xM31i6S4e8sUqbPzZXmz82V6+eskSfnrNNV+9O0Zy8fL/ZbkWxbKXnz/ksW+lcTDZgvRbxneqLw0TrgNYicinu3kLDgV/MxygijVX1sOfhECDBB7m8rmrVquzZs4e///3v3HPPPQQEBHDkyBGfrDsv38W6fRksik9icUIS+9NPAO5x/x+9tg0DIsNpG17bxvsxxgA+OGegqnki8gCwAPev/vdVNU5EnsVdoeYA40VkCJAHpAOjvZ3LWzZv3szbb7/Nm2++SVhYGDt27CAoKMgn6846ncfKnSks8vT/P3oyl6pVArj6sobc86uWDGgfbt0/jTFF8sl1Bqo6D5hX6LlJBe4/gfvwUbl18uRJnn32WV5++WUaNmzII488QuvWrb1eCJKPn2JxfDKL4o+wencaOfku6tUMon/7MK5pH06fNqHUqmbXFhpjzs2+JcrAihUrGDt2LLt27eKuu+7i5Zdf9uroontTs1kYd4QFcUfYlOge/K1pgxqM6tmcayLDiW5enyo28bsx5gJYMbgIqkp+fj5jx44lPz+fxYsX079/f6+sJ/7wMRZsO8KCuKQzM391jKjDQwPacE1kOO0a2fF/Y0zpWTEopa+++op+/fpRq1YtvvrqK5o2bUrNmmXXHz/fpWzcn8GCbUf4Ju4IBzJOEiDQtUUDJg2O5NoO4db/3xhTZqwYXKBDhw7xwAMP8MUXX/Dyyy/z6KOP0rZt2zJ575w8F9/tSWVBXBKL4pNIzTpN1cAArm7VkAf7tWJA+3Ab+98Y4xVWDErI5XIxffp0/vjHP5KTk8OUKVPKZGC57NN5rNiZwoK4IyzdnszxU3nUrBpITLswBnZoREzb0DNTPxpjjLdYMSihRx99lGnTptG3b1/effddWrduXer3ysjOYXFCEgvikvh2Vwqn81zUrxnEdR0bMbBDI65uZQPAGWN8y4rBOeTm5nLixAnq1q3LuHHjiIyM5O677y7VidojR0+xMP4IsetOsmPhYvJdyiV1qzOiWzMGdmhE1xbWA8gY4xwrBsVYv349Y8aMoXXr1nz66ae0a9eOdu3aXdB77EvNZkGc+wTwpv2ZADSqJdzTpyWDOjbi8oi61gPIGOMXrBgUkp2dzeTJk5k2bRrh4eFMnjy5xK9VVbYfOc4329zXAGw/8r8uoI9e24ZBHRtxIH4DffteWFExxhhvs2JQwKZNm7jlllvYu3cv48aNY+rUqdSrV++cr3G5lM0HMs90Af0p7QQiEN28Pk8NjmRgoS6gByrlLA7GGH9nxaCAiIgIGjduzAcffECfPn2KbZeX72Lt3nS+2XaEhfFHSDp2mqBA4arLQrj3V+5JYEJrWxdQY0z5UamLgaryySef8J///IfPPvuMsLAwVq9eXWTbU7n5rN6dyjfbjrAoIYnME7lUDwqgb5swBnVsREy7MOrWsC6gxpjyqdIWg8TERO677z7mzp1LdHQ0aWlphIWF/aJN9uk8lu9I4Zu4IyxNSCI7J5/a1aswoH04Azs04ldtQqlR1bqAGmPKv0pXDFwuF2+99RaPP/44LpeLV155hQkTJhAY6P5SP3oylyUJSczfdoSVO93XADSoVZUbO1/CoI6NuOqyEKpWsS6gxpiKpdIVg9OnT/Pqq6/Ss2dP3nnnHS699FLSsk6zKP4g87cd4bs9qeTmK4091wAM6tiIri0aEGiTwBtjKrBKVwxq1KjBypUroUZdFsYnMXHRGtbuTcOl0KxBTe66+lIGdWxE5yb1CLACYIypJHxSDERkEPAq7pnOpqvqlELLqwEzgCggDRimqvu8keXrrYf55+q9bNifgSq0Cgvm/phWDOrYiMjGdewiMGNMpeT1YiAigcCbwDXAAWCdiMxR1YI97u8GMlS1lYgMB6YCw7yR51DmSbJz8nloQBuu69iI1uG1vbEaY4wpV3yxZ9AN2K2qPwKISCxwE1CwGNwEPO25Pwt4Q0REVbWsw9zd61LG9mlZ1m9rjDHlmi+6xUQAiQUeH/A8V2QbVc0DjgINvRHGzgMYY8zZxAs/vn+5ApFbgYGqOsbzeBTQTVUfLNAmztPmgOfxHk+btELvNQ4YBxAeHh4VGxtb5DqzsrIIDg72xp9z0Sxb6Vi20vPnfJatdC4mW0xMzAZVjT5rgap69Qb0BBYUePwE8EShNguAnp77VYBUPIWquFtUVJQWZ9myZcUuc5plKx3LVnr+nM+ylc7FZAPWaxHfqb44TLQOaC0il4pIVWA4MKdQmznAHZ77Q4GlntDGGGN8wOsnkFU1T0QewP3rPxB4X1XjRORZ3BVqDvAP4F8ishtIx10wjDHG+IhPrjNQ1XnAvELPTSpw/xRwqy+yGGOMOZsNsmOMMcaKgTHGGB90LfUWEUkBfipmcQjuHkn+yLKVjmUrPX/OZ9lK52KyNVfV0MJPltticC4isl6L6kfrByxb6Vi20vPnfJatdLyRzQ4TGWOMsWJgjDGm4haDd50OcA6WrXQsW+n5cz7LVjplnq1CnjMwxhhzYSrqnoExxpgLYMXAGGNM+S4GIjJIRHaIyG4RebyI5dVE5GPP8rUi0sKPso0WkRQR2ey5jfFhtvdFJFlEthWzXETkNU/2rSLSxY+y9RWRowW226Si2nkhV1MRWSYiCSISJyITimjjyHYrYTantlt1EfleRLZ4sj1TRBtHPqclzObY59Sz/kAR2SQic4tYVrbbraihTMvDDfegd3uAlkBVYAsQWajNfcDbnvvDgY/9KNto4A2Htl0foAuwrZjl1wPzAQF6AGv9KFtfYK4D26wx0MVzvzaws4h/U0e2WwmzObXdBAj23A8C1gI9CrVx6nNakmyOfU49638YmFnUv11Zb7fyvGdwZjpNVc0Bfp5Os6CbgA8992cB/cU3M96XJJtjVHUl7tFhi3MTMEPd1gD1RKSxn2RzhKoeVtWNnvvHgQTOnrHPke1WwmyO8GyLLM/DIM+tcK8VRz6nJczmGBFpAtwATC+mSZlut/JcDPxqOs1SZAO4xXM4YZaINPVBrpIqaX6n9PTs2s8XkQ6+Xrlnd/xK3L8kC3J8u50jGzi03TyHOjYDycAiVS12u/n4c1qSbODc5/RvwJ8AVzHLy3S7lediUFQFLFzVS9LGG0qy3q+AFqraCVjM/yq8P3Bqu5XERtxjq3QGXgdm+3LlIhIMfAb8QVWPFV5cxEt8tt3Ok82x7aaq+ap6BdAE6CYiHQs1cWy7lSCbI59TERkMJKvqhnM1K+K5Um+38lwMDgAFq3QT4FBxbUSkClAX3xyCOG82VU1T1dOeh+8BUT7IVVIl2baOUNVjP+/aq3uejCARCfHFukUkCPeX7Ueq+nkRTRzbbufL5uR2K5AhE1gODCq0yKnP6XmzOfg5vRoYIiL7cB9m7ici/y7Upky3W3kuBv48neZ5sxU6ljwE93FefzEHuN3TO6YHcFRVDzsdCkBEGv18XFREuuH+fzjNB+sV3DPyJajqX4tp5sh2K0k2B7dbqIjU89yvAQwAthdq5sjntCTZnPqcquoTqtpEVVvg/v5YqqojCzUr0+3mk5nOvEH9eDrNEmYbLyJDgDxPttG+yAYgIv/B3bskREQOAJNxnzxDVd/GPSvd9cBu4ARwpx9lGwr8XkTygJPAcB8V+KuBUcAPnmPMABOBZgWyObXdSpLNqe3WGPhQRAJxF6BPVHWuP3xOS5jNsc9pUby53Ww4CmOMMeX6MJExxpgyYsXAGGOMFQNjjDFWDIwxxmDFwBhjDFYMjDHGYMXAmIviubBwhed+FxFREWnoGfNmm4jUdDqjMSVhxcCYi5OJe9hogAeBNUB93KNNLlLVEyJS36lwxpSUFQNjLs5RoKaINMR9Retq3MVgHPCWp800h7IZU2JWDIy5CKr68/DCY3EPD3Ac6AQEqupOERkEtBORR53KaExJWDEw5uK5cA9i9gVwDHgUeNuzLBX4t6r+xaFsxpSIFQNjLl4OMN8zwcgxoBbw85y1nXBPe2qMX7OB6ozxIs+Il7cAU1TVn4YpN+YXrBgYY4yxw0TGGGOsGBhjjMGKgTHGGKwYGGOMwYqBMcYYrBgYY4zBioExxhisGBhjjMGKgTHGGOD/AdPuY24JPxdJAAAAAElFTkSuQmCC\n",
                        "text/plain": [
                            "<Figure size 432x288 with 1 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "plt.plot(wvec, wvec, linestyle='--', color='black')\n",
                "plt.plot(wvec, wvec*ζvec)\n",
                "plt.xlabel(\"$w_t$\")\n",
                "plt.ylabel(\"$c_t$\")\n",
                "plt.title(\"Initial Guess\")\n",
                "plt.grid()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## Residuals in the model's equations \n\nTo identify the unknown decision functions for $\\zeta$ and $h$, we use two modelp's equations,\nnamely, the definition of normalized Lagrange multiplier and the FB function representing the KT conditions, respectively:\n\n$$\n\\begin{gather*}\nh=\\beta \\overline{r}E_{\\epsilon }\\left[ \\frac{u^{\\prime }\\left(\nc^{\\prime }\\right) }{u^{\\prime }(c)}\\exp \\left(   \\delta ^{\\prime }-\\delta\n+r^{\\prime }\\right)  \\right] , \\\\\n\\Psi_{FB}\\left( 1-\\zeta ,1-h\\right) =0\n\\end{gather*}\n$$\n\nwhere $\\epsilon'=(\\epsilon'_r,\\epsilon'_\\delta,\\epsilon'_q,\\epsilon'_p)$.\n\nWe do not need to include the definition $\\zeta = \\frac{c}{w}$ because we will impose it to hold exactly in the solution by setting $c=w\\zeta $ and $c^{\\prime }=w^{\\prime }\\zeta ^{\\prime }$.\n\nWe next construct the residuals in the above two equations which we will minimize. For given vectors of current state $s=(r,\\delta ,q,p,w)$ we use next-period shocks $\\epsilon'$ to compute realized residuals:\n\n$$\n\\begin{equation*}\n\\begin{matrix}\nR_1(s,\\epsilon')=\\beta \\overline{r} \\frac{u^{\\prime }\\left(\nc^{\\prime }\\right) }{u^{\\prime }(c)}\\exp \\left( \\delta ^{\\prime }-\\delta\n+r^{\\prime }\\right) -h,\n\\\\ \nR_2(s)=\\Psi_{FB}\\left( 1-\\zeta ,1-h\\right), \n\\end{matrix}\n\\end{equation*}\n$$\n\nwhere the transition equation is $w^{\\prime }=\\left( w-c\\right) \\overline{r%\n}\\exp (r)+\\exp (y)$."
        },
        {
            "cell_type": "code",
            "execution_count": 11,
            "metadata": {
                "lines_to_next_cell": 0
            },
            "outputs": [],
            "source": [
                "def Residuals(e_r: Vector, e_δ: Vector, e_q: Vector, e_p: Vector, r: Vector, δ: Vector, q: Vector, p: Vector, w: Vector):\n",
                "\n",
                "    # all inputs are expected to have the same size n\n",
                "    n = tf.size(r)\n",
                "\n",
                "    # arguments correspond to the values of the states today\n",
                "    ζ, h = dr(r, δ, q, p, w)\n",
                "    c = ζ*w\n",
                "\n",
                "    # transitions of the exogenous processes\n",
                "    rnext = r*ρ_r + e_r\n",
                "    δnext = δ*ρ_δ + e_δ\n",
                "    pnext = p*ρ_p + e_p\n",
                "    qnext = q*ρ_q + e_q\n",
                "    # (epsilon = (rnext, δnext, pnext, qnext))\n",
                "    \n",
                "    # transition of endogenous states (next denotes variables at t+1)\n",
                "    wnext = tf.exp(pnext)*tf.exp(qnext) + (w-c)*rbar*tf.exp(rnext)\n",
                "\n",
                "    ζnext, hnext = dr(rnext, δnext, qnext, pnext, wnext)\n",
                "    cnext = ζnext*wnext\n",
                "\n",
                "\n",
                "    R1 = β*tf.exp(δnext-δ)*(cnext/c)**(-γ)*rbar*tf.exp(rnext) - h\n",
                "    R2 = min_FB(1-h,1-ζ)\n",
                "\n",
                "    return (R1, R2)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## The expected squared sum of residuals\n\nWe construct the objective function for minimization as the squared sum of\ntwo residuals in the two model's equations on a given 5-dimensional domain $s=(r,\\delta ,q,p,w)$: \n\n$$\n\\begin{equation*}\n\\Xi (\\theta )=E_{s}\\left[ \\left( E_{\\epsilon }\\left[ R_1(s,\\epsilon')\\ \\right] \\right) ^{2}+v\\left( R_2(s)\\right) ^{2}\\right] ,\n\\end{equation*}\n$$\n\nwhere $v$ is the exogenous relative weights of the two residuals in the\nobjective function. We placed the first residual $R_1(s,\\epsilon')$ under the expectation\noperator $E_{\\epsilon }$ across next-period shocks $\\epsilon' =\\left(\n\\epsilon' _{r},\\epsilon' _{\\delta },\\epsilon' _{q},\\epsilon' _{p}\\right) $ as is\nrequired by the definition of $h$; the second residual $R_2(s)$ does not include random variables  \nand requires no expectation operator. The value of the objective function $\\Xi (\\theta )$\ndepends on the coefficients $\\theta $ because these coefficients determine the choices via  \n\n$$\n\\begin{equation*}\n\\left( \n\\begin{matrix}\n\\zeta  \\\\ \nh%\n\\end{matrix}%\n\\right) =\\varphi (s;\\theta ).\n\\end{equation*}\n$$\n\nA shortcoming of the constructed objective function is that it requires\na potentially costly evaluation of two nested expectation operators: for\neach random grid point $s=(r,\\delta ,q,p,w)$, we need to construct a\nseparate approximation of the expectation function $E_{\\epsilon }\\left[\n R_1(s,\\epsilon')\\right] $ by considering a potentially large\nnumber of next period shocks $\\epsilon' =\\left( \\epsilon' _{r},\\epsilon'\n_{\\delta },\\epsilon' _{q},\\epsilon' _{p}\\right) $. In particular, if there are \n$n$ grid points and $J$ next-period shocks, we have $n\\times J$ function\nevaluations."
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "## All-in-one expectation function\n\nIn the paper, we introduce a technique which we call an all-in-one expectation operator\nthat makes it possible to merge the two expectation operators into a single\none. This technique relies on a simple result from probability theory that\nsays that for two random variables $a$ and $b$, which are independent and\nfollow the same distribution, we have $E[a]^{2}=E[a]E[b]=E[ab]$. \n\nTherefore, we replace $\\left( E_{\\epsilon }\\left[  R_1(s,\\epsilon')\\right] \\right) ^{2}$ by the product of two residuals constructed\nby using two uncorrelated random draws $\\epsilon' _{1}$ and $\\epsilon' _{2}$,\nand therefore, we can pull the expectation out of squares \n\n\n$$\n\\begin{equation*}\nE_{\\epsilon _{1}}\\left[ R_1(s,\\epsilon'_1)\\right]\nE_{\\epsilon _{2}}\\left[  R_1(s,\\epsilon'_2)\\right]\n=E_{\\epsilon _{1},\\epsilon _{2}}\\left[ R_1(s,\\epsilon'_1)R_1(s,\\epsilon'_2) \\right]. \n\\end{equation*}\n$$\n\nWith that result, we can re-write the objective function as just one\nexpectation operator: \n\n$$\n\\begin{equation*}\n\\Xi (\\theta )=E_{s,\\epsilon _{1},\\epsilon _{2}}\\left[ \\underset{\\xi (\\omega\n;\\theta )}{\\underbrace{R_{1}\\left( s,\\epsilon'\n_{1}\\right) R_{1}\\left(\ns,\\epsilon' _{2}\\right) +v\\left(\nR_{2}\\left( s\\right) \\right) ^{2}}}\\right] \\equiv E_{\\omega }\\left[ \\xi\n(\\omega ;\\theta )\\right],\n\\end{equation*}\n$$\n\nwhere $\\omega =(s,\\epsilon' _{1},\\epsilon' _{2})$. Therefore, we wrote the\nobjective function of the deep learning method as a single expectation\noperator $E_{\\omega }\\left[ \\xi (\\omega ;\\theta )\\right] $ of a function $%\n\\xi (\\omega ;\\theta )$ that depends on a vector-valued random variable $%\n\\omega $. We approximate $\\Xi (\\theta )$ by using Monte Carlo simulation:\n\n$$\n\\begin{equation*}\n\\Xi (\\theta )\\approx \\Xi ^{n}(\\theta )=\\frac{1}{n}\\sum_{i=1}^{n}\\xi (\\omega\n_{i};\\theta ),\n\\end{equation*}\n$$\n\ni.e., we draw $n$ random draws of $\\omega =(s,\\epsilon' _{1},\\epsilon' _{2})$\nand compute the average of the objective function"
        },
        {
            "cell_type": "code",
            "execution_count": 12,
            "metadata": {},
            "outputs": [],
            "source": [
                "def Ξ(n): # objective function for DL training\n",
                "\n",
                "    # randomly drawing current states\n",
                "    r = tf.random.normal(shape=(n,), stddev=σ_e_r)\n",
                "    δ = tf.random.normal(shape=(n,), stddev=σ_e_δ)\n",
                "    p = tf.random.normal(shape=(n,), stddev=σ_e_p)\n",
                "    q = tf.random.normal(shape=(n,), stddev=σ_e_q)\n",
                "    w = tf.random.uniform(shape=(n,), minval=wmin, maxval=wmax)\n",
                " \n",
                "    # randomly drawing 1st realization for shocks\n",
                "    e1_r = tf.random.normal(shape=(n,), stddev=σ_r)\n",
                "    e1_δ = tf.random.normal(shape=(n,), stddev=σ_δ)\n",
                "    e1_p = tf.random.normal(shape=(n,), stddev=σ_p)\n",
                "    e1_q = tf.random.normal(shape=(n,), stddev=σ_q)\n",
                "\n",
                "    # randomly drawing 2nd realization for shocks\n",
                "    e2_r = tf.random.normal(shape=(n,), stddev=σ_r)\n",
                "    e2_δ = tf.random.normal(shape=(n,), stddev=σ_δ)\n",
                "    e2_p = tf.random.normal(shape=(n,), stddev=σ_p)\n",
                "    e2_q = tf.random.normal(shape=(n,), stddev=σ_q)\n",
                "    \n",
                "    # residuals for n random grid points under 2 realizations of shocks\n",
                "    R1_e1, R2_e1 = Residuals(e1_r, e1_δ, e1_p, e1_q, r, δ, q, p, w)\n",
                "    R1_e2, R2_e2 = Residuals(e2_r, e2_δ, e2_p, e2_q, r, δ, q, p, w)\n",
                "\n",
                "    # construct all-in-one expectation operator\n",
                "    R_squared = R1_e1*R1_e2 + R2_e1*R2_e2 \n",
                "    \n",
                "    # compute average across n random draws\n",
                "    return tf.reduce_mean(R_squared)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "So far, we have been using TensorFlow in the eager execution mode as if it was numpy: result of each operation is computed immediately."
        },
        {
            "cell_type": "code",
            "execution_count": 13,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "<tf.Tensor: id=676, shape=(), dtype=float32, numpy=0.21182436>"
                        ]
                    },
                    "execution_count": 13,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "n = 128\n",
                "v = Ξ(n)\n",
                "v"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "Note that the intermediate results are still stored as special tensorflow objects (tensors) and they can be converted to a regular value easily."
        },
        {
            "cell_type": "code",
            "execution_count": 14,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "0.21182436"
                        ]
                    },
                    "execution_count": 14,
                    "metadata": {},
                    "output_type": "execute_result"
                }
            ],
            "source": [
                "v.numpy()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "# Model training\n\nWe are now ready to perform minimization of the objective $\\Xi_n$, hence to solve (or to train) the model using stochastic optimization - the stochastic gradient descent method, and in particular, its version called Adam.\n\nSome readers might have noticed that we have never defined in this notebook a vector of coefficients $\\theta$ in the code and that the objective function $\\Xi_n$ doesn't seem to depend on it, but on $n$ only.\n\nThis is not quite true. The coefficients $\\theta$ have actually been created, and initialized with random values. They are attached to the perceptron object and stored in the computational graph created by TensorFlow. We can easily recover these variables and print these values."
        },
        {
            "cell_type": "code",
            "execution_count": 15,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "[<tf.Variable 'dense/kernel:0' shape=(5, 32) dtype=float32, numpy=\n",
                        "array([[ 0.27676374,  0.39069772,  0.26160532,  0.22172564, -0.14331418,\n",
                        "         0.05529916, -0.2997829 ,  0.12813216,  0.2777946 , -0.14456508,\n",
                        "        -0.14910114,  0.18854314,  0.2899565 ,  0.11027104, -0.00051203,\n",
                        "         0.37444663,  0.01240587,  0.04889569, -0.02504221, -0.31196538,\n",
                        "        -0.05168563, -0.25079435,  0.34215814, -0.0137023 ,  0.32700056,\n",
                        "         0.2783717 , -0.25808793,  0.292719  ,  0.25731295,  0.21283293,\n",
                        "        -0.26468152, -0.2413572 ],\n",
                        "       [-0.18395486, -0.26174784, -0.03679186, -0.27664554,  0.19304872,\n",
                        "        -0.29115328,  0.0895026 , -0.05567339,  0.06700936,  0.3385952 ,\n",
                        "         0.20031953, -0.16950646, -0.0635725 ,  0.01131013, -0.25897276,\n",
                        "        -0.26392013,  0.13591498,  0.25545698, -0.18999511,  0.1419881 ,\n",
                        "         0.25043184, -0.40154755, -0.09462759, -0.28182727, -0.16937858,\n",
                        "        -0.33178714, -0.39249578,  0.01656675,  0.38928336, -0.0028846 ,\n",
                        "        -0.12537593,  \n"
                    ]
                }
            ],
            "source": [
                "θ = perceptron.trainable_variables\n",
                "print( str(θ)[:1000] ) # we truncate output"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "Now it is time to choose an optimizer. In TensorFlow, the optimizer object is in charge of performing the optimization steps, given the computed gradient. For the stochastic gradient descent, the updating rule would be:\n$$\\theta \\leftarrow \\theta(1-\\lambda) - \\lambda\\nabla_{\\theta} \\Xi_n(\\theta)$$\nwhere $\\lambda$ is a learning rate.\nFor Adam, the learning rate evolves over time and can be specific to each coefficient."
        },
        {
            "cell_type": "code",
            "execution_count": 16,
            "metadata": {
                "lines_to_next_cell": 2
            },
            "outputs": [],
            "source": [
                "from tensorflow.keras.optimizers import Adam, SGD"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 17,
            "metadata": {},
            "outputs": [],
            "source": [
                "variables = perceptron.trainable_variables\n",
                "optimizer = Adam()\n",
                "# optimizer = SGD(λ=0.1) # SGD can be used in place of Adam"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "We use the GradientTape functionality from TensorFlow to compute the gradient from the objective, and supply it to the optimizer. Also, to speed up computation, we specifically instruct TensorFlow to run each training step in a \"graph\" mode. This is done by using a`tf.function` decorator. As a result, all the computations will be pre-optimized in the graph and will be ran efficiently in a scalable way."
        },
        {
            "cell_type": "code",
            "execution_count": 18,
            "metadata": {},
            "outputs": [],
            "source": [
                "@tf.function\n",
                "def training_step():\n",
                "\n",
                "    with tf.GradientTape() as tape:\n",
                "        xx = Ξ(n)\n",
                "\n",
                "    grads = tape.gradient(xx, θ)\n",
                "    optimizer.apply_gradients(zip(grads,θ))\n",
                "\n",
                "    return xx"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "We repeat the training_step K times with the following function (each training step is called an epoch)."
        },
        {
            "cell_type": "code",
            "execution_count": 19,
            "metadata": {},
            "outputs": [],
            "source": [
                "def train_me(K):\n",
                "\n",
                "    vals = []\n",
                "    for k in tqdm(tf.range(K)):\n",
                "        val = training_step()\n",
                "        vals.append(val.numpy())\n",
                "        \n",
                "    return vals"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 20,
            "metadata": {},
            "outputs": [
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "100%|██████████| 50000/50000 [01:11<00:00, 703.45it/s]\n"
                    ]
                }
            ],
            "source": [
                "# with writer.as_default():\n",
                "results = train_me(50000)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "On a modern CPU, optimization should be done within a few minutes. It would be dramatically faster on hardware adapted to deep-learning. To see how the training has performed, we can plot the empirical errors against the number of epochs."
        },
        {
            "cell_type": "code",
            "execution_count": 21,
            "metadata": {},
            "outputs": [
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "/opt/pablo/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt\n",
                        "  \"\"\"Entry point for launching an IPython kernel.\n"
                    ]
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXgUVdo28Pvpzp6QQBISCCGEEBYhgEDYRCSgyCbquMwAOjOoI6KjfvPq6ID7KCgz7o6++uKouCPiBoLsRER2kFX2PSAkkAU6+3K+P7KQrZPeq7rr/l0Xl11Vp6qfOnaerj516hxRSoGIiHyfSesAiIjIM5jwiYgMggmfiMggmPCJiAyCCZ+IyCCY8ImIDMJP6wCaEh0drRITEx3aNz8/H6Ghoa4NyMuwDlgHRj9/wJh1sHXr1nNKqdb11+s64ScmJmLLli0O7Zueno60tDTXBuRlWAesA6OfP2DMOhCR442t12WTjoiMF5HZeXl5WodCROQzdJnwlVILlVJTIiIitA6FiMhn6DLhExGR6zHhExEZBBM+EZFB6DLh86YtEZHr6bJbplJqIYCFqampdzt7rPIKhbKKCpSVK5RVKJSVV6C8QqG0QqG8vGpbharaXu91uUJ5hYKCQvtWIegQFYoAP11+RxIRNUuXCb/a7lN56PrEDw7tW1ZegYqli+DK4f79TIKEqBB0jglDcvW/1i3QKSYUIQG6rkoiIn0n/OiwQEwekujQvhknTqJTxw4wm0zwMwv8TAI/s6nqv1XLNdtMMJuk1ra6+1QohePn83Eo04JDmRYczLRgxd5MlFdc+jZp1zK45ksgvlUwosICERUagKiwAESGBiAyJAB+Zv46ICLt6Drht4kIwvQxlzm0b3r6WaSldXVZLH0TWtVZLimrqPMlcCjLgoNnLdh49DyKSisaPUbLEH9EhgYgvlUIUuLCkdIuAilxEWgfGQwRaVBeKYWS8goEmE2NbicisoeuE76eBfiZ0Dm2BTrHtqizvqJCIbewFNn5xThnKUF2fgnO55fgvKW48rWlBEfP5WP2miMoq/qFEB7khx5xEUiIDEFOQQmyLMU4ZylG1sViFJVWJvyWIf5oGeKPViEBGN4tBrcP6oCwwMb/9ymlMG/LSZzJKsMwpfhlQUQAdJrwRWQ8gPHJyclah2I3k0kqm3BCA5AcY71ccVk5DpyxYPfpPOw+lYfdpy9g5b5MRIUGoHWLQHRICEHrFoGICPbHxeIy5OaXIqegBGcuFGHWD/vwdvph3DmkIyZfkYiIEP+a45ZXKDz53W58tvEEAGDJ6bW4L60TxqS0YZMSkcHpMuG7speOXgX6mdEzPgI94+0fPmL7yVy8ueogXl1xAO/+dARje7bB9b3bITWxFf7x1U58t/007hmWhPLsDKw6U44HPv8FbcKDMGFAe4zq0QYdo0MR5G92w1kRkZ7pMuFT0y5v3xL//XN//Hr6At5bexSLd53BvC0ZCPI3oai0Ao+M6oq/Dk9GevpZTJ84DKv2ZeLjDcfx2oqDeG3FQYhU3mROah2GhMhgRAT7o0NUKK7vHccvAiIfxoTvxbrHhePl3/fGzNIUrNqXicW7fsOQ5GhMHJBQU8ZsEozsHouR3WNxMrsA20/m4nCWBUey8nHknAU7M3JxobAUFQp4cel+DOvSGhUVCpMGJiA1MVLDsyMiV2PC9wFB/maM7dkWY3u2bbJc+8gQtI8MabBeKYX1R87j7fTDWHvwHApLy7F49294/8/9cUVyNDIvFmHLsRyMSWnDG8BEXowJnyAiuKJTNK7oFA0AOGcpxqR3N+DODzfj7dv74ZVlB7DrVB6+uvcK9OvQqsH+ZeUVuOWd9Zh8RSJu7NPO0+ETkY3YbYMaiA4LxOd3D0JiVCju+GAzdp3KQ4CfCR+tP4bCknLM23wSb646WPPg2c5Tedh+Mhcf/HwUQOUXxujX1mDDkfMangUR1ccrfGpUVFggPrt7EO79ZCsGdIyEpbgMn2w4jk1Hs/FbXhEAYNepPBSUlKNlSAAAYEdGHo5kWfDVtgzsO3MRczedQNfYFsjIKazpjfTtL6dw7Hw+/nZNF83OjciodJnwvbkfvi+JDA3AF/cMBgAcybLg4/XHERroh8/uHoiVezPx3tqjEAGUquz1czqvECNe/rFm/2+3n8bSPWdRWFqOVQ8PQ1LrMPzti+0AgEkDEhATHqTJeREZlS4TvhH64XubpNZhWP33NMSEByLQz4zBSVGYOCAB247n4NGvdmJcr7bo074lDpy1oLS8ArERQXjy290oLC0HALy56lCdq/oFO07j2u5t8I+vduLFW3shvlXDm8lE5Fq6TPikT7V7+IgIkmPCkBQdCktxGUantEFcy2CM6Vm5XSmFTq1D0SMuAq8s248P1x/H17+cqtn/xaX7MWPRXgDAswt/xTu394PJxB5ARO7Em7bkFJNJcOeVHRHXMrjO+uqePxHB/rg1tX2dbX0SWqK47NIAc8t+PYvHvtmFvMJSnLMUY2dGrkdiJzIaXuGT2/WIC8f1veNwWdtwxIYHIju/BL+cyMW4nm3x0q298bcvfsHczSex/WQuMnIKYSkuw5pHhqNFkB/CgvzgZxL2/ydyASZ8cjsRwRsT+9Qs5+SXYN+Zi3hs7GUIDjDjuRtSsHTPWew7c7GmzCPzd2Dj0WyM7tEGS/acwTu398XolLZ4ZsEenMguwPuT+2txKkRejU065HGtQgPw0q29ERla2Z0zJjwIT4yrO+/BxqPZAIAle84AAKZ+sg1HsiyYs+4YVu3LxCNf7kBxWblnAyfycrzCJ124sU87bDyajVv6xeNkdkHNDd3aanf5/HJrBoZ3i8Gmo9mIDgvA/SM6ezJcIq/EhE+6EB0WiHf/lFqzPGlgAno9s6xmkpjG3PfptprXTPhEzWOTDulSSIAfwoP966wb0NH66J0jXkrHwbOV9wDyi8uQebHIrfEReSMmfNKtKVclAQAC/So/pp/9ZSAWPXglxvVqi6/vu6JO2SPn8jHy1TU4lVuIHk8vxYCZK/HGyoMej5lIz3TZpMOhFQgApg7rhHuuSsKZC0U4bymBn9mEHnEReGtSXwDA1d1isHJfZp19hsxaVfP6leUHAAD7hpTXdO008+EuMjBdXuErpRYqpaZERNg//R/5FhFB24hgpLRr+FkY3q2JSYNr+XDdMQx7Mb3OlwGREenyCp/IFrcNTEC7VsEYnBSFQ5kWXPeftQCAkAAzCkouddn8+fB5nMotBABk55fUdAclMhpdXuET2UJEMLxrDIL8zWhfNfjajBtTsO3JkXXKrTmQVfP68W92eTRGIj3hFT75hIgQfxybNa5muX1kME5mFzYodyjTgsKSchSVlqMVr/TJYHiFTz7pH6O71byeOqxTzeuDmRZc9tQS9HluObo+8QMSpy3CjpMcrI2MgQmffNJ1veJw9IWx+GBUCKaN6YYWgQ1/zFaP2PnUd7s9HR6RJpjwyWeJXBpls6S8wmq5nIJSHMmyYOvxbFwoKvVUeEQex4RPhvDhnQNww+VxjW47kV2AES//iJvfXo8pH23xcGREnsOET4YwKCkKr0+4NETz3mdHN1puw5FsZOQUeCosIo9iwidDeX9yKr6cOhjBAWarZa7812oUlJR5MCoiz2DCJ0MZ0S0W/RMrB2H77q9DrJbbmZHnqZCIPIYJnwyrd/uWmHFjSqPbJszeAKWsD81M5I2Y8MnQbhuYgBUPXYX9M0Y36Lp555zNKG9iPH4ib8OET4YmIkiOaYFAPzM+/svAOttW78/C/5v7C0qb6NJJ5E08lvBFJElE3hOR+Z56TyJ79GoXgVv6xddZ9/3O39D58R9wKNOCxGmL8NPBLCt7E+mfTQlfRN4XkUwR2V1v/WgR2S8ih0RkWlPHUEodUUrd5UywRO5kMkmdYRhqu+aVyvl0//jeJk+GRORStl7hzwFQp+OyiJgBvAVgDIDuACaKSHcR6Ski39f7Z9vA5UQaCzA3/yexdM8ZD0RC5Hpia08EEUkE8L1SKqVqeTCAZ5RSo6qWpwOAUuqFZo4zXyl1SxPbpwCYAgCxsbH95s6da1N89VksFoSFhTm0r69gHThWB+tOl6FLKxMyLlbgtW3FjZaZMzrUFeG5HT8DxqyD4cOHb1VKpdZf78zwyO0AnKy1nAFgoJWyEJEoADMB9BGR6da+GJRSswHMBoDU1FSVlpbmUHDp6elwdF9fwTpwrA5ql+6RchZ3NzLcQnz3VOQVlqJP+5Yw6XjaRH4GWAe1OZPwG/uUW/25oJQ6D2CqE+9H5HHDu7ZudH11m/6AxEh8cc+gmkHaiPTMmV46GQDa11qOB3DauXAqich4EZmdl8enHUlbfmaT1UHXAGDTsWx8tP64ByMicpwzCX8zgM4i0lFEAgBMALDAFUFxEnPSkxdu6tnk9qcX7PFQJETOsbVb5ucA1gPoKiIZInKXUqoMwP0AlgLYC2CeUoqffPI5IQF+mD91MP46vPEumwCQOG0RisvKrW4n0gObEr5SaqJSqq1Syl8pFa+Ueq9q/WKlVBelVCel1ExXBcUmHdKb1MRIPDKqG1b/Pc1qma5PLGGXTdI1XQ6twCYd0quO0U13x7zn463YejwHZy8UeSgiItvpMuET6dnBmWPQvW241e03v70OA59f6cGIiGzDhE9kJ3+zCf++pVez5U7lFnogGiLb6TLhsw2f9C6lXQSOzRqHZ2/oYbXMkFmrMOj5lVh36JwHIyOyTpcJn2345C3+OKgDDj8/FgvvvxIp7Ro285y5UIRJ/93IyVRIF3SZ8Im8hYjAbBL0jI/AO7f3s1qu4/TFOJR50YORETXEhE/kIvGtQjC2Zxur2695ZQ2yLjY+GBuRJ+gy4bMNn7zV/97WD4+N7WZ1e/+ZK7DnND/XpA1dJny24ZM3m3JVJxx9YazV7ePeWOvBaIgu0WXCJ/J2IoKDM8dgRLfG5/5JnLYIidMW4ccDWTUTpReWcGgGci9nhkcmoib4m0148ZZe6DdjhdUyf36/7pSJr/y+N27qG2+lNJFzeIVP5EZRYYE2PaRV7aF5O9wYDRmdLhM+b9qSL/l9anvMuDHF7v3m/HwUf3qfk6aT6+gy4fOmLfma2wd1wP/90Xo//dpy8ksAAM8s/BVrDmS5MywyGF0mfCJfNKhjFCKC/Zst1+e55Th49tJDWnxKl1yFCZ/IQyJC/LHj6Wux4qFhzZYd+eqamtcHzlrcGRYZCBM+kYeVVVTUvH54ZJdmy496bQ0Spy3CnXM2uzMsMgAmfCIP69Q6DGNS2mDxg0PxwNWdbd5v1b5MPDxvB05mF7gxOvJlukz47KVDvszfbMLbt/dD9zjrk6hY89W2DNz76VaczC5A3+eWY/epyr8RpRTb+qlZukz47KVDRrLnn6Owf8Zom8vvPnUBQ/+9Gtn5JbjuP2txKrcQHacvxr2fbHNjlOQL+KQtkcZCA537M0zfnwkAWMIJ1KkZurzCJzKit2/rCwCY+bsUBPnb/qf5+De73RUS+Rhe4RPpxJiebXFs1jgAwG0DOyBx2iK7j/FbXiGOZuVDAeib0MrFEZK3Y8In0qlNj12Nez7Zil9O5Nq8z+AXVtW8Ht87Dje3dUdk5K3YpEOkUzHhQfjmviEO779q71kUlikOu0w1mPCJdG7Rg1fiQTv661fLLynHvSsKcNlTS9wQFXkjXSZ89sMnuqRHXAQeGtkFG6Zf7fAxOJcuATpN+OyHT9RQbHggAKBnO/v/LvrPXMEHs4g3bYm8hYhgw/Sr0SrUH4F+Zlz3n5+w+9QFm/fvOH0x7hzSEb3bR+BCYSluH9QBIuLGiElvmPCJvEibiKCa198/MBRDZq3CqdxCm/d//+ejNa9DAvxwcz9Op2gkumzSISLbrHl0uMP7Pvxl5XSKxWXldn1pkPfiFT6RFzObBE9e1x3Pff+rQ/vXfrjryPNjYTKxiceX8QqfyMvddWVHHJs1DsdmjcMdQxJxq4PNNM8s3IOKCt7Y9WVM+EQ+5OnxPfDirb0d2vej9cc5AJuPY8In8kH/urmnQ/vd9+m2OvPpkm9hwifyQX/on4AXb+mFYfH236Yb+eoabD6W7YaoSGu8aUvko25NbY/WlsP4MaPM/n3fWY9WIf7IKSjFgRljEOBnwl8+3IwNR7Kx+5+j3BAteYIur/A5tAKR66T/PQ0/Txth9345BaUAgBeX7gMArNibCUux/V8epB+6TPgcWoHIdRKjQ9GuZbDD+7/701EMmbWq+YKke7pM+ETkejdcHode8RGYclWS3fvWfjCruOzScMv5xWU4b+HAbN6CbfhEBvH6hD41r0+cL3C4C2bXJ5Zg33OjEeRvxqjX1iAjp7Bmpi7SN17hExnQm5P64MERyQ7v/+ryAygqLUdGDodk8CZM+EQG5Gc24aFru2LFQ8Mc2v//1hxBtyfrTqxSUFKGzItFrgiP3IQJn8jAkmPCXHasG9/6GQNmrnTZ8cj1mPCJDK56HB5nHThrcUE05E5M+EQEAHjp1t5oEeRYP47V+zJdHA25AxM+EQEAbukXj13PjEKA2f60cMeczTWvE6ctwm95DW/mVlQoFJaUN1hPnsOET0R1HJg5BvOnDnbqGINfWFUzcfqJ8wUor1B4adl+XPbUEj6tqyEmfCJqIDUxEsdmjcOmx6/GhP7tHTpG/5krsP/MRVz14mq8snw/vt52CgBwobDUlaGSHZjwiciqmBZBmHVzL2x67GqH9h/12hoAwLrD56FQObkK503XDp+0JaJmBQeYndr/UKYFF4sqm3IEzPha4RU+ETWrRZA/Nj12Ne4c0tGh/auTPcArfC15LOGLyI0i8q6IfCci13rqfYnINWLCg/DU+O548OrOTh2n+mYueZ5NCV9E3heRTBHZXW/9aBHZLyKHRGRaU8dQSn2rlLobwGQAf3A4YiLS1EMjuzi1/3X/WYui0nIs2X0GczedaLLsyezKHj7kGra24c8B8CaAj6pXiIgZwFsARgLIALBZRBYAMAN4od7+dyqlqp/MeKJqPyIyqNrj8EwYkACgcthlkwj8q54DyMgpwNB/r8bUYZ0wbUw3TeL0NaKUbd+eIpII4HulVErV8mAAzyilRlUtTwcApVT9ZF+9vwCYBWC5UmpFE+8zBcAUAIiNje03d+5cW8+lDovFgrAw140T4o1YB6wDd53/L5llOJZXgasT/PHg6gKnjvXBqBCICCYvyUdcqOD5oSEAgMO55XhuQxGSIkx4arDjE7gY8TMwfPjwrUqp1Prrneml0w7AyVrLGQAGNlH+AQDXAIgQkWSl1DuNFVJKzQYwGwBSU1NVWlqaQ8Glp6fD0X19BeuAdeCu8699xO59LuKaV9Y4fKw7lhbgnquSABzB6XxVE2/EiRxgwzq0CA9HWtoQh4/vzs9AXmEplu05g1tTHXtWwdOcuWnb2L12qz8XlFJvKKX6KaWmWkv2ROR9kmNa4Nu/Op6Qgcrhlr3Ro/N34JH5O7HntGPzb3+68TgSpy1CUalnhpxwJuFnAKj9tRYP4LRz4VTiJOZE3uXy9i3x8q293XNwpZCTX+KeYzsps6rHUVFphUP7v77iIIDKXwqe4EzC3wygs4h0FJEAABMALHBFUJzEnMj73NwvHu/+qUGzsdN2ZOShz3PLcSjzosuP7Swbb4Hqhq3dMj8HsB5AVxHJEJG7lFJlAO4HsBTAXgDzlFJ73BcqEeldl1jnb45WX/XWd+zcpZvDWReLcbHI+8fkqf6+qN0+nldQigtuOjebbtoqpSZaWb8YwGKXRoTKJh0A45OTHZ9zk4g8LyEyBFOuSsL1veNwyzvrHGrqeHXFAfTt0BI7TubWWW+qdXnaf+YKtAzxx/antHuGc9+ZC05/6ahGMn7vZ5cBgFsmhtfl0Aps0iHyTiKCx8ZehpR2Edj33Bj87RrHnsr943ub8NKyAw2OXVtugfVkO/KVH/HpxuMOvbetRr/2Ew5n5VfF5uhRqgaU89D4QrpM+ETkG+5LS8YLN/V0ybHMdmTVg5kWPP7N7uYLaqz6Cn/1/kzY+kyUM5jwichtAvxMmDggAcdmjcP9w51rojU1kfCtdWtcd+gcKmxMpB/8fBQ/7PqtZjk7vwQbj5y3L0g7VUf26PydmLv5ZJNlXUGXCZ/dMol8z99HdcX+GaMd3t9kJd+vOZCFbk8uwZZj2Q22TfrvRnx/pPl29oycAvxz4a+499Ntl/Z9dwP+MHuDw/HaovZV/fSvd+Gxb3a59f10mfDZhk/kmwL9zNjyxDUO7Xs8u/EhHNYeOgcA+OngOZSVN7xJfDa/MqmesxRb/SVwKqfhHLz7ztjeDdRVLfCfbWx6MDln6TLhE5Hvigj2d2i/6V/vQuK0RViw49LznUWl5dhe1Zvn9ZUHcf9nv1jtr586YwUmvdv4FXtTjT7Otq2vO3QOidMWYfmvZ+16X3dgwicij3L2avjBz3+peX0o04JNRy815SzZc6bBuD61m/63najb1bOaO++XTvrvRgDA84v3evR9G6PLhM82fCLf1dTNV3u9/ePhZsvY8m6qiWttW5JybiNDIyilkL4/s5HSdct4ki4TPtvwiXyXK6c4XLTzt+YL1VNaXoHV+zPrDnjWSN6tjvPxb3c3O5bPr6cv4OVl++sk8GW/nsXkDzbXLB89l99gP0836XAScyLyqNoPUH0xZZDbe8KcLahA4rRFNcudH/+hzvZJAxNwbffYBvsJKhPy51WzcjX1PMGLS/cDqBzyYdbNvSrf90KRk5G7ni6v8InIGAYmReHAjDFufY8DOU0P7/DZxhNYc+BcnXX5xWWoPbPi+sOV289eKML0r3ehtJHeQADq9KVvqrUmI6cAv+UVevwSn1f4RKQJc1XH+gA/E2bd1BPTvnZvH3R7/JZXt5vmsfMFeGv1IWw8mo01B7KQ1rW1Q8fNyClAVGggrvzXagBAWKBnU7AuEz4HTyPybU+MuwxXdo6uWZ4wIAFZF4vx8vIDTezlGcVljffVr262ac7HG47jcKYF87Y0fHK2OtFX401b8KYtka/7y9AkdGsTXmfd7/q20yiaur10uj6xpImSlZqaoerJb3djzrpjKChpfhYr9sMnIkOKbxWi2Xt/8POxOstz1h1rtFy1939uerteMeETkW6M7x2ndQgAgE82ND3EQf2x+h3FB6+IyLC6tWmhdQgeVdhE05C1ewnOYMInIt24d1gnzJ86GPcMS9I6FM1ttzIMhDN0mfA5tAKRMZlMgtTESEwakKB1KD5JlwmfvXSIjK1DVCi2PTlS6zB8ji4TPhFRZGgAkqJDtQ5DM/Xn8HUFJnwi0i/PzO2tS27I90z4RER65I7vOiZ8ItKtNyf2xdBaQzCQc5jwiUi3useF490/pQIAftdHu6EXfIUuB08jIqoW5G/G9qdGokWQPx4a2QVD/726+Z18QFmF6x/D1eUVPvvhE1FtLUMCYDYJ2kdqN96OpxmmDZ/98ImoOX0SWmodgtfRZcInIrJmRLcYAMBHdw6Av9l3+226Y1w1tuETkVd578+pUKpyGIbtT12LHk8v1Tokr8ErfCLyKiICU9X0iKGBfvhgcn+NI/IeTPhE5NWGd4vBsv+5SuswXM4wN22JiOzRJdb3xtF3Rxs+Ez4R+YSfHh2OxCjjdNt0BBM+EfmE9pEh+MtQTpzSFCZ8IvIZ1c0gkwZyApXGMOETkc8Y2DESADA2pS2+nDqYA6/Vo8uEz6EViMgRXWJb4NiscbiyczT6J0bi47sGah2Srugy4XNoBSIi19NlwicicqU/De6gdQh2Yz98IiI7zRkdimdvSNE6DF1gwiciQ+gVzyZiJnwiMoQF91+Jt2/ri5Yh/lqHohkmfCIyjDE922Jsz7Zah6EZJnwiMpR/Xt9D6xA0w4RPRIbib/aOtPdbXpHLj+kdZ05EZDBFpeUuPyYTPhEZVnJMmNYhWCVu6IjPhE9EhrXioWFah2CVuOHRKyZ8IjKcoZ2jMaBqoLUF9w/ROJrGJUaHuvyYnMSciAyn9qBqveJbYssT1yB1xgoNI2rIHfeWeYVPRIYXHRaIYH+z1mHUIW5oxPdYwheRy0TkHRGZLyL3eup9iYi8kUmrhC8i74tIpojsrrd+tIjsF5FDIjKtqWMopfYqpaYC+D2AVMdDJiLyfVqOljkHwOjaK0TEDOAtAGMAdAcwUUS6i0hPEfm+3r+Yqn2uB7AWwEqXnQERkQ9yR7dMUUo1XwqAiCQC+F4plVK1PBjAM0qpUVXL0wFAKfWCDcdapJQaZ2XbFABTACA2Nrbf3LlzbYqvPovFgrAw/fax9QTWAevA6OcP2F4HU5bno6TqWac2IYIzBdZzY8cIE47mVbgqxEY9MzgIiRGO3VcYPnz4VqVUg5YUZ3rptANwstZyBgCr84mJSBqAmwAEAlhsrZxSajaA2QCQmpqq0tLSHAouPT0dju7rK1gHrAOjnz9gex283Oo0Hvj8FwDAX6/tgSe/3W21bHh4OJCX66oQG5WamoqUdq4d0tmZhN/YDw6rX4lKqXQA6U68HxGR24zvHYd1h8/j800nmm0/d0f7eoP30NmTthkA2tdajgdw2rlwKnEScyLSQnUTt7UeMsv+5yoAqHloy5309qTtZgCdRaSjiAQAmABggSuC4iTmRKSFipqEDyyvSu61dYltgVUPD8Ojo7q6PRbNrvBF5HMA6wF0FZEMEblLKVUG4H4ASwHsBTBPKbXH9SESEXlGdR8WEaBzbItGk25S6zD4mU0Y0S3GrbG4I+Hb1IavlJpoZf1iNHED1lEiMh7A+OTkZFcfmojIqqFdWuPLrRnoEVfZuiCwfmMyNjzQrbHEtAhy+TF1OZaOUmohgIWpqal3ax0LERnH9b3jkNa1NcKDKue9FZFLl/31PD6uOzpEheJIlgXztmS4PJbI0ACXH1OXCZ+ISCvVyR5oujdOWKAfpg7rBABuSfjuwIRPRGRFdTt6UutQ9O/QfM+cdi2DcSq30M1ROU6Xo2WyWyYR6UF118jFDw7Fv27pZbXcuJ5tAQC/69POI3E5SpcJn90yiUgXbOwp8/qEy7HzmWvxt2s6o1e8bXnryuRoJwJzjC4TPhGRHtjaM9LPbEJ4kD/8zCZ8OU65HbwAAAWTSURBVHVwk2XH9ar8NeDqYRNswYRPRGRFdRu+jWNMVu7TzNdE97bhTkTkHF0mfLbhE5EeVCdvZX2YMMeP7YkBeerRZcJnGz4R6YFDV/h2JPLBSVH2BeQkXSZ8IiI96BgdCsB10w3WTvACIK1ra5cc11bsh09EZMVHdw7AjoxcBAe4boJzWyedcgde4RMRWREVFogR3WLt2qf2b4FWIZee2k3r2hrP3tDjUjkBRna379jO0mXC501bIvIFEcH+6BFX2Stnzh0D0Dm2RZ3tSa0vTb0466aebo9Hl006HDyNiLyV1Gvv//q+K1BW3nwzjicmVdFlwici8gX3DOuEQD8zAmtl2vpN+K9PuBzxrYI9Eg8TPhGRC9W+vp84IKGJcpUlb7i8cvydrIvF7gwLABM+EZEutG4RiJUPD0NIgBlHsvLd8h5M+ERELtRcl/2mWvM7Vd3EbRvhniYe9tIhItIAh1aowqEViIhcT5cJn4jIW9XvlqknTPhERB6k4cgKTPhERFrQ4ncAEz4RkUEw4RMRGYQuEz67ZRKRrxqUVDlmzqBOnp38BNBpwme3TCLyZgM7RuJ/b+vb+LakKOyfMRpXdIr2cFR80paIyOW+uGdwk9sD/Vw3oYo9dHmFT0RErseET0RkEEz4REQGwYRPRGQQTPhERAbBhE9EZBBM+EREBsGET0RkEKK0HKuzGSKSBSAXQO0xFiKaWK79OhrAOReGU/99nSlrbXtj621Z54k6sOf8bSnvyjpoqj60qgNbyjpTB778d2Btm5HqwNb11pY7KKVaN9hbKaXrfwBm27pc7/UWd8bhTFlr2xtbb8s6T9SBPefv6Tpopj40qQNbyjpTB778d8A6sH19c8v1/3lDk85CO5brb3NnHM6Utba9sfW2rPNEHdh7XE/WQXOfEVdx5WegqTL2/j+vv+ztfwfWthmpDmxdb9dnX9dNOs4QkS1KqVSt49AS64B1YPTzB1gHtXnDFb6jZmsdgA6wDlgHRj9/gHVQw2ev8ImIqC5fvsInIqJamPCJiAyCCZ+IyCAMk/BFJFREPhSRd0XkNq3j8TQRSRKR90RkvtaxaEVEbqz6//+diFyrdTxaEJHLROQdEZkvIvdqHY9WqvLBVhG5TutYPMmrE76IvC8imSKyu9760SKyX0QOici0qtU3AZivlLobwPUeD9YN7Dl/pdQRpdRd2kTqPnbWwbdV//8nA/iDBuG6hZ11sFcpNRXA7wH4TFdFO3MBAPwDwDzPRqk9r074AOYAGF17hYiYAbwFYAyA7gAmikh3APEATlYVK/dgjO40B7afv6+aA/vr4Imq7b5iDuyoAxG5HsBaACs9G6ZbzYGNdSAi1wD4FcBZTwepNa9O+EqpNQCy660eAOBQ1RVtCYC5AG4AkIHKpA94+XlXs/P8fZI9dSCV/gXgB6XUNk/H6i72fg6UUguUUlcA8JmmTTvrYDiAQQAmAbhbRHwiH9jCT+sA3KAdLl3JA5WJfiCANwC8KSLj4N7HrrXW6PmLSBSAmQD6iMh0pdQLmkTnGdY+Aw8AuAZAhIgkK6Xe0SI4D7H2OUhDZfNmIIDFGsTlSY3WgVLqfgAQkckAzimlKjSITRO+mPClkXVKKZUP4A5PB6MBa+d/HsBUTwejEWt18AYqv/iNwFodpANI92wommm0DmpeKDXHc6Hogy/+lMkA0L7WcjyA0xrFogWjnz/AOgBYBwDroAFfTPibAXQWkY4iEgBgAoAFGsfkSUY/f4B1ALAOANZBA16d8EXkcwDrAXQVkQwRuUspVQbgfgBLAewFME8ptUfLON3F6OcPsA4A1gHAOrAVB08jIjIIr77CJyIi2zHhExEZBBM+EZFBMOETERkEEz4RkUEw4RMRGQQTPhGRQTDhExEZBBM+EZFB/H85/o8e6aygCAAAAABJRU5ErkJggg==\n",
                        "text/plain": [
                            "<Figure size 432x288 with 1 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "plt.plot(np.sqrt( results) )\n",
                "plt.xscale('log')\n",
                "plt.yscale('log')\n",
                "plt.grid()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "In the training graph, one can see that errors are decreasing until they reach $2*10^{-5}$ on average (the errors are volatile both because they depend on a specific random draw and because only two such random draws are used to approximate the expectation function). The numbers in the graph represent the mean of the squared residuals. If we take the square root, we get $5*10^{-3}$ for the $L_2$ norm, wich is still quite accurate approximation for a model with 5 states variables and a kinked solution. We show the constructed decision rule below."
        },
        {
            "cell_type": "code",
            "execution_count": 22,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEYCAYAAAC+xZqSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUdfb/8dchPYRQQggICO6iaMBFCKJIEZRFcBVUUNRVbHxRFyw/O+66llVXV10suLIoFrAERQRFLKyABRUEpCkWEOkdQno/vz/uDQ7DJIQkM3eSnOfjMY/MzL2Z+547d+6593PvfK6oKsYYY+q3Bl4HMMYY4z0rBsYYY6wYGGOMsWJgjDEGKwbGGGOwYmCMMQYrBgGJiIpIhwqGfyci/coZ1k9ENldmXC+IyJ9F5OMQTOdlEXkw2NOpTUQkW0R+53WOyhKRPiLyowfTTRaRH0Uktpzh94nIq9WcRqXem4jcLSIvVGdaNUVEFojIqEqMFyMiP4hIiyN5/TpVDETkVxEpFJHmfs8vd1fw7avwmoes1FS1k6ouqMz/H8m4oaCqr6nqQK9ziEgrEZksIttEJMtdeO8XkYZeZ6sJgb64qpqgqr8EYVrRIvKEiGx2C856ERlf3ddV1c9VtWNNZDxCdwEvqWp+sCZQ2femqg+r6mFXwOFEVQuAF4E7j+T/6lQxcK0HLil7ICInAnHexTH+RKQZ8BXO59JTVRsBfwSaAL/3MlstNQ7oDvQAGgH9gW89TVRFIhIDXAFUa8u/thGRyBp+ydeBK9z5WTmqWmduwK/A34BvfJ57HPgroEB797kFwCifca4EvvB5rEAHYDRQBBQC2cB7PtMZ4N6PA14G9gHfA7cDm/0ylY17H/AWzoKeBawCjsP5Mu8ENgEDff63MTAZ2AZsAR4EInwzu+9vH04RHOz3nn5xp7Me+HM57/U04Btgv/v3NJ9hC4B/AAvd1/kYaO4z/C1gu/u/nwGdfIa9DDxYzuf0oPveG1TwWVYpFxDrzt89QIb7vyn+n4XP5/Gqe7+9+7lf5X4O+4DrgJOBle5rTfCbvwuBZ9yMPwBnusMeAkqAfJzlZoLvcuXz2U4BdgEbcJbbBpX5bAPMq9nAzRUMvwtY586r74Hz3edj3PfV2WfcZCAPaAH049Bl+TZ3fuwHpgGxPsPvwFlWtwKj/N7v2e60s3CW5dvKydoXWOv33DHAp+7/zgUmlH1u7vBTgS/d97IC6OczrBnwkptpHzDTfd7/vd3p5soCfvT5LO/zm9YQ4Dt3WguAEyo7f/zeU9nyMx7Yi/Od8J9We3ceRpaz3roaWOO+r4+Adn7T+Bk4vbLrz7q4Z/A1kCgiJ4hIBDCCKm5lqOok4DXgX+rs4p8bYLR7cbZmfw+chbNVU5FzgalAU5ytt49w9tBaAw8A//UZ9xWgGKcwdQUG4nzJypyCs+A2B/4FTBZHQ+BpnBVII5wV63L/IO4W+vvuuEnAv4H3RSTJZ7RLcVaQLYBonIW9zAfAse6wZTjzqjIGADNUtTTQwGrmugJnRdvW/d/rcFZulXUKznsaATyJsyExAOgEXCQip/uN+wvO/L8XmCEizVT1r8DnwFh3uRkbYDrPuDl/B5wOjHTfj+9rH/LZlpP5a+AWEfmLiJwYYLx1QB93evcDr4pIK3WaE2bgsycNXAR8qqo7y5nWRcAgnBX0H3BWaojIIOAWnHnVwX1PviYD17rLY2dgXjmvf6L7vn29DizFmRf/wOc7JiKtcZaVB3FW/LcBb4tIsjvKVCAe5/NrgbPyPYiIdATGAie7+c7CWbH7j3cc8AZwM07RnAO8JyLRPqMFnD/lKFt+WuBsQFSaiJwH3A1c4Gb53M3maw3QpbKvWReLATgLwEicpocfcCp+sFwEPKSqe1V1E84KrCKfq+pHqlqMs2WdDDyiqkVAOtBeRJqISAowGGeLL8f9co4HLvZ5rQ2q+ryqluAUjlZAijusFOgsInGquk1VvwuQ5U/Az6o6VVWLVfUNnPnlW/ReUtWfVDUPeBM4qWyAqr6oqlnuSuU+oIuIND7cDMNZSW+rYHh1chW5r99BVUtUdamqZlYiU5l/qGq+qn4M5ABvqOpOVd2C84Xr6jPuTuBJVS1S1Wk4K7E/HW4CPhsp49z59yvwBHC5z2gVfbb+/gk8CvwZWAJsEZEDK0xVfUtVt6pqqZvzZ5wmJXBWtL7F4FL3ufI87b7WXuA9fpvvF+F8Jt+pai5O0fFVBKSKSKKq7lPVZeW8fhOcrXMARORonL2ze1S1QFU/c6db5jJgjqrOcd/fXHcenC0irXC+Q9e50yxS1U8DTLMEZy8pVUSiVPVXVV0XYLwRwPuqOtf9vj6O0zJwWiXmTyBbVfUZdxk/kg0WgGuBf6rqGndd8jBwkoi08xknC2d+VkpdLgaX4lTlKUGe1lE4zQplNhxm/B0+9/OA3e4XvuwxQALQDogCtolIhohk4Ow1+J4hsL3sjvsFBEhQ1RycBfc69//fF5Hjy8nun3cDzl7KIdMAct1siEiEiDwiIutEJJPftqQOOnhfjj04K7fyVDkXzmf/EZAuIltF5F8iElWJTGX8Px//xwk+j7eouz/uk/GoSkyjOc7ejO97LPf9+X627lkw2e7tO3d4iao+q6q9cL78DwEvisgJACIy0j2Jomw56sxvn9M8IE5ETnFXJCcB71SQvbz57v898L0PMAynqWiDiHwqIj3Lef19OMc9yhwF7HOX6TK+860dcGHZe3PfX2+c5astsFdV91XwflDVtThb+/cBO0UkXUQCfY4HLZfunu0mKrdcBuI/j45EO+Apn/e8FxC/LI1wmrMqpU4WA1XdgNPOejbObrC/HJxdxzItK3q5w0xuG85CV+boymSshE1AAU5beBP3lqiqnSrzz+7exx9xvhQ/AM8HGG0rzkLl62gqtyd1KTAUp1mgMU77JjgL5OH8DzhfRMpb/qqcy936u19VU3G22M7B2UuEI/vcK6O1X5PM0TjZoeLlZjfOlrLve6zs+/vcbXpKCLQsqGqeqj6Ls1JNdVfwz+M0gySpahNgNe7n5K7Q3sTZO7gUmK2qWf6vWwnbgDY+j32/E6jqN6o6FGdjZqY7zUBW4hxH833dpn5nmfl+xzYBU32+I01UtaGqPuIOayYih906VtXXVbU3zmeiOHta/g5aLt3Pvi1Vb3nwX0aOZPnchNPs5vu+41T1S59xTsA5hlIpdbIYuK4BzvDboiizHLhAROLF+T3BNRW8zg6cdt3yvAmME5GmItIGuKHKiX2o6jacA6NPiEiiiDQQkd/7tVkHJCIpIjLE/QIV4BzELAkw6hzgOBG5VEQiRWQEkIpzQPJwGrmvvQdnAX64cu8McI4BJAKvlO3WikhrEfm3iPyhOrlEpL/bbh4BZOKsdMve+3LgYhGJEpHuwPAjyBxIC+BG9/UuxPnyzXGHlbvcuHuCbwIPiUgjdx7cQhWPbYnIzeL8viXOnV9X4Hw+3wINcVY6u9xxr8LZM/D1Os6e5J+puImoIm8CV7nH6uKBv/vkixbn9y2N3eaVTAIvjwCLgSbusYCyDbslwP3u6/Tm4ObCV4FzReQsd2811p0Xbdzv0AfAf9zvZ5SI9PWfoIh0FJEzxDnzJh9nDzBQvjeBP4nIme7e5q0434EvA4xbFcuBviJytNvcOq6CcSfirHc6ue+hsbsM4j5ujXMM5evKTrzOFgNVXaeqS8oZPB7nDKEdOO2xFR34nIyzhZUhIjMDDL8fZ9dxPc7Ke2rVUx9iJE5zwvc4W3rTqbh5pUwDnAV1K87u4+nAX/xHUtU9OFvOt+Ks1O8AzlHV3ZWYxhSc973FzVfphc5tTz0NZ0W9SESygE9wzsBYW81cLXHmUybOAbRP+W0lew/Ogf59OJ9bVVd8ZRbhHGzejdM0M9zNDvAUMFxE9olIoONIN+BsCf6Cc+bQ6zjnhldFHs4xh+1uljHAMFX9RVW/d4d9hbO8n4hzFssBqrrIzXIUzsrziKnqBzjHy+YDa93pgbOyBOd4yK9uk+J1OG39gV6nEOdMNN/hl+IcbN2Lc6B+is/4m3D2UO/GKXibcM7oK1u3XY6znP2Ac4zn5gCTjQEewZl323GK/N0Bsv3o5nrGHfdc4Fw3c7W5xzum4ewdLaWCjR9VfQdn7yXdnaercY6PlLkUeMU9nlcpcnCTpzGmMkTkSpzT/Hp7nSUcuccrVgMx7gHOI/nfsrNjulbhwGq95+7hrAD6avlnhR2izu4ZGGNCS0TOd5tymuJstb53pIUAQFV3qerxVgiqxj3r6vgjKQRgxcAYU3OuxWmqWYfT5n69t3HMkbBmImOMMbZnYIwxxoqBMcYYoKZ7yguZ5s2ba/v27QMOy8nJoWHD8OwJ2bJVjWWrunDOZ9mqpjrZli5dultVkw8ZoEfYM2i43NLS0rQ88+fPL3eY1yxb1Vi2qgvnfJataqqTDVii9aTXUmOMMUfIioExxhgrBsYYY6wYGGOMIYTFwO1R8FsROaTzJRGJEZFpIrJWRBZJFS5cb4wxpupCuWdwE04vkoFcg3MBiw44PYoG6kvcGGNMkISkGLj9/P8JeKGcUYbidCUNTvfDZ/pdNMQYY+q90tJS8vKC039fSPomEpHpONdpbQTcpqrn+A1fDQxS1c3u43XAKerXf72IjAZGA6SkpKSlp6cHnF52djYJCRVdbc47lq1qLFvVhXM+y1Z5W7Zs4fHHH6dly5aMGTOmytn69++/VFW7HzIg0I8PavKGc5GS/7j3++FcVs9/nO+ANj6P1+Fcos9+dBZClq1qwjmbanjns2yHV1RUpI899pjGxsZqYmKiPv/887X2R2e9gCEi8iuQDpwhIv6X99uMe81UEYnEuabu3hBkM8aYsJWfn0+vXr24/fbbGThwIN9//z2jRo0KyrSCXgxUdZyqtlHV9sDFwDxV9b/k3bvAFe794e441re2MaZeKi0tBSA2Npb+/fuTnp7OzJkzad26ddCm6dnvDETkAREZ4j6cDCSJyFqcC4Pf5VUuY4zx0tdff03Xrl1ZssS5hPsjjzzCiBEjCPY5NSHttVRVFwAL3Pt/93k+H7gwlFmMMSac5OTkcM899/Dkk0/Spk0bcnJyQjp9+wWyMcZ47JNPPuHEE09k/PjxXH/99axevZrTTz89pBlq7fUMjDGmrvjqq6+IjIzk008/pW/fvp5ksD0DY4zxwKxZs/jwww8BuPPOO1mxYoVnhQCsGBhjTEjt2LGDESNGcN555/H0008DEBUVRVxcnKe5rBgYY0wIqCqvvvoqqampzJw5kwcffJBZs2Z5HesAO2ZgjDEhMHfuXC6//HJ69uzJ5MmTOeGEE7yOdBDbMzDGmCApLS1l1apVAPzxj39k+vTpfP7552FXCMCKgTHGBMVPP/1Ev3796NmzJ1u3bkVEGDZsGBEREV5HC8iKgTHG1KDi4mIeffRR/vCHP7Bq1SomTJhAq1atvI51WHbMwBhjakh+fj69e/dm6dKlnH/++Tz77LO1ohCA7RkYY0y1+XYsd9ZZZ/HWW28xY8aMWlMIwIqBMcZUy5dffkmXLl345ptvAHjooYcYPny4x6mOnBUDY4ypguzsbG666SZ69+5NZmZm0C5HGSpWDIwx5gjNnTuXzp0788wzzzBmzBhWr17taVcSNcEOIBtjzBFavHgxsbGxfPbZZ/Tu3dvrODXC9gyMMaYSZsyYwZw5cwC44447WL58eZ0pBBCCYiAisSKyWERWiMh3InJ/gHGuFJFdIrLcvQXnIp/GGHOEtm/fzvDhwxk2bBjPPvss4HQsFxsb63GymhWKPYMC4AxV7QKcBAwSkVMDjDdNVU9yby+EIJcxxpRLVXnllVdITU1l9uzZ/POf/2TmzJlexwqaoBcDdWS7D6Pcm13s3hgT1pYsWcKVV15Jp06dWLFiBXfddRdRUVFexwqakBwzEJEIEVkO7ATmquqiAKMNE5GVIjJdRNqGIpcxxvgqLS1l5cqVAHTv3p133nmHTz/9lI4dO3qcLPhENXQb6SLSBHgHuEFVV/s8nwRkq2qBiFwHXKSqZwT4/9HAaICUlJS09PT0gNPJzs4mISEhGG+h2ixb1Vi2qgvnfOGUbePGjTz22GOsXbuWqVOnEhsbGzbZ/FVnvvXv33+pqnY/ZICqhvQG3AvcVsHwCGD/4V4nLS1NyzN//vxyh3nNslWNZau6cM4XDtkKCwv14Ycf1piYGG3atKm+8sorWlpaGhbZylOdbMASDbBODfrvDEQkGShS1QwRiQMGAI/6jdNKVbe5D4cAa4Kdyxhj8vPz6dWrF8uWLWP48OFMmDCBlJQUr2N5IhQ/OmsFvCIiETjHKN5U1dki8gBOhXoXuFFEhgDFwF7gyhDkMsbUU6WlpTRo0IDY2FgGDx7MX//6Vy644AKvY3kq6MVAVVcCXQM8/3ef++OAccHOYowxX3zxBaNHj+bll1/m5JNP5sEHH/Q6UliwXyAbY+qFrKwsxo4dS58+fcjLy6OgoAAR8TpW2LBiYIyp8z766CM6derEf/7zH2688UZWrVpFnz59vI4VVqyjOmNMnbd06VIaNmzIF198wWmnneZ1nLBkewbGmDpHVZk+ffqBjuVuv/12li9fboWgAlYMjDF1yrZt2xg2bBgXXnghzz33HOB0LBcTE+NxsvBmxcAYUyeoKi+99BKpqal88MEHPProo7zzzjtex6o17JiBMaZOmDt3LldffTV9+vThhRde4LjjjvM6Uq1iewbGmFqrpKSE5cuXA/DHP/6RWbNmsWDBAisEVWDFwBhTK61Zs4Y+ffrQq1cvtm7diogwZMgQGjSw1VpV2FwzxtQqRUVFPPjgg5x00kn8+OOPTJo0iVatWnkdq9azYwbGmFojPz+fnj17snz5ckaMGMHTTz9NixYtvI5VJ9iegTEm7JWUlAAQGxvLkCFDmDlzJunp6VYIapAVA2NMWPvss8/o3LkzixY5F0i8//77GTp0qMep6h4rBsaYsJSZmcn111/P6aefTkFBAcXFxV5HqtOsGBhjws6cOXPo1KkTkyZN4pZbbmHVqlX06tXL61h1mh1ANsaEnZUrV5KYmMj06dM55ZRTvI5TL9iegTHGc6rKtGnTmD17NgC33XYby5Yts0IQQkEvBiISKyKLRWSFiHwnIvcHGCdGRKaJyFoRWSQi7YOdyxgTHrZs2cJ5553HxRdfzKRJkwCIjIy0juVCLBR7BgXAGaraBTgJGCQip/qNcw2wT1U7AOOBR0OQyxjjIVVl9uzZpKamMnfuXB5//HHrWM5DQS8G6sh2H0a5N/UbbSjwint/OnCm2PXojKnTPv74Y5544gm6devGypUrufXWW4mIiPA6Vr0lqv7r5SBMRCQCWAp0AJ5V1Tv9hq8GBqnqZvfxOuAUVd3tN95oYDRASkpKWnp6esDpZWdnk5CQUOPvoyZYtqqxbFUXTvlKSkpYt24dxx13HKrKggUL6NevX1heizic5pu/6mTr37//UlXtfsgAVQ3ZDWgCzAc6+z3/HdDG5/E6IKmi10pLS9PyzJ8/v9xhXrNsVWPZqi5c8q1atUp79Oih8fHxumXLFlUNn2yB1NVswBINsE4N6dlEqpoBLAAG+Q3aDLQFEJFIoDGwN5TZjDHBUVhYyP3330+3bt345ZdfeOGFF6xjuTAU9N8ZiEgyUKSqGSISBwzg0APE7wJXAF8Bw4F5bgUzxtRieXl5nHrqqaxcuZJLL72UJ598kuTkZK9jmQBC8aOzVsAr7nGDBsCbqjpbRB7A2V15F5gMTBWRtTh7BBeHIJcxJkhKSkqIiIggLi6O888/n4ceeohzzjnH61imAkEvBqq6Euga4Pm/+9zPBy4MdhZjTPDNnz+fa6+9lilTpnDqqady3333eR3JVIL9AtkYUyP279/PtddeyxlnnOF7coipJawYGGOq7f333yc1NZUXXniB22+/nRUrVtCzZ0+vY5kjYB3VGWOqbfXq1SQlJTFr1iy6dz/0FHYT/mzPwBhzxFSV119/nffeew+AW2+9lSVLllghqMWsGBhjjsjmzZs599xz+fOf/8zkyZMBp2O56Ohoj5OZ6rBiYIyplNLSUv773/+SmprK/PnzGT9+PG+//bbXsUwNsWMGxphKmTt3Ltdddx1nnnkmkyZN4ne/+53XkUwNsj0DY0y5iouLWbJkCQADBw5kzpw5zJ071wpBHWTFwBgT0KpVq+jZsyd9+/Zl69atiAiDBw8Oyx5GTfVZMTDGHKSgoIC///3vdOvWjQ0bNvDyyy9bx3L1gB0zMMYckJeXR48ePVi9ejWXX34548ePJykpyetYJgRsz8AYQ3FxMQBxcXFcdNFFzJkzhylTplghqEesGBhTz33yySccf/zxfPXVVwDcc889DB482ONUJtSsGBhTT2VkZDBq1CgGDBhAgwYNaNDAVgf1mX36xtRD7733Hqmpqbz88svcddddrFixglNOOcXrWMZDdgDZmHpozZo1tGjRgvfee4+0tDSv45gwEPQ9AxFpKyLzRWSNiHwnIjcFGKefiOwXkeXu7e+BXssYUzWqytSpU5k1axYAt9xyC998840VAnNAKPYMioFbVXWZiDQClorIXFX93m+8z1XVrotnTA3bsWMHf/rTn/jggw8477zzGDp0KJGR1ihgDhaKy15uA7a597NEZA3QGvAvBsaYGlRaWsrEiRO57bbbaNCgAU8//TRjxozxOpYJUxLKS9OJSHvgM6Czqmb6PN8PeBvYDGwFblPV7wL8/2hgNEBKSkpaenp6wOlkZ2eTkJBQw+lrhmWrGst25BYvXsydd97JSSedxJ133knLli29jnSIcJ13UHez9e/ff6mqHnrhibJrlQb7BiQAS4ELAgxLBBLc+2cDPx/u9dLS0rQ88+fPL3eY1yxb1Vi2yikqKtLFixerqmppaal+9NFHOm/ePI9TlS+c5p2/upoNWKIB1qkhObVURKJwtvxfU9UZAQpSpqpmu/fnAFEi0jwU2YypK8pODz399NPZtm0bIsLAgQOtYzlTKaE4m0iAycAaVf13OeO0dMdDRHq4ufYEO5sxdUF+fj5/+9vf6N69O1u2bGHq1KnWsZw5YqE4paAXcDmwSkSWu8/dDRwNoKoTgeHA9SJSDOQBF7u7M8aYCuTl5dG9e3e+//57Ro4cyfjx42nWrJnXsUwtFIqzib4AKtxPVdUJwIRgZzGmriguLiYyMpK4uDguvfRSunfvzllnneV1LFOLWXcUxtQyH3/8MR07djzQsdxf//pXKwSm2qwYGFNL7Nu3j6uuuoqzzjqLqKgoIiIivI5k6hArBsbUAjNnziQ1NZWpU6dy9913s3z5cnr06OF1LFOH2G/SjakF1q5dS6tWrfjggw846aSTvI5j6iDbMzAmDKkqr7zyCjNnzgTg5ptvZtGiRVYITNBYMTAmzGzYsIFBgwZx5ZVX8uqrrwIQGRlJVFSUx8lMXWbFwJgwUVpayjPPPEOnTp348ssvmTBhAm+++abXsUw9YccMjAkTc+fO5cYbb2TQoEFMnDiRdu3aeR3J1CO2Z2CMh4qKili0aBEAAwcOZO7cucyZM8cKgQk5KwbGeGTZsmX06NGDfv36HehYbsCAAdaxnPGEFQNjQiwvL49x48bRo0cPtm/fzmuvvWYdyxnP2TEDY0IoLy+Pbt268cMPP3D11Vfz+OOP07RpU69jGWPFwJhQKCoqIioqiri4OEaOHMnJJ5/MgAEDvI5lzAHWTGRMkH344Yccd9xxLFy4EIBx48ZZITBhx4qBMUGyZ88eRo4cyeDBg4mLiyM6OtrrSMaU64iKgYjc4nO/Y83HMaZumDFjBqmpqbzxxhvcc889fPvtt5x88slexzKmXJU6ZiAiTYDxQEcRyQdWAtcAV1Xif9sCU4CWQCkwSVWf8htHgKeAs4Fc4EpVXXYE78OYsPLLL7/Qtm1bPv74Y7p06eJ1HGMO67B7BiLSHShS1auA+4FFwLHAIRe2L0cxcKuqngCcCowRkVS/cQa7r3ksMBp4rpKvbUxYUFVefPFF3nnnHVSVm2++ma+//toKgak1KtNMNAUoAVDVj4ANwC5Vfa8yE1DVbWVb+aqaBawBWvuNNhSYoo6vgSYiYidem1ph69atDBw4kGuuuYbXX38dcDqWi4y0k/VM7VGZYpCvqvllD1R1N/BAVSYmIu2Brjh7F75aA5t8Hm/m0IJhTFgpKSnhySef5JprrmHRokVMnDiRadOm2S+ITa0kqlrxCCLTgcmq+oHPc6tVtfMRTUgkAfgUeEhVZ/gNex/4p6p+4T7+BLhDVZf6jTcapxmJlJSUtPT09IDTys7OJiEh4UjihYxlq5pwzLZ48WLuvPNOunfvzh133EFycrLXkQIKx3lXxrJVTXWy9e/ff6mqdj9kgKpWeANaAcuB14Ebgf8Csw73f36vEQV8BNxSzvD/Apf4PP4RaFXRa6alpWl55s+fX+4wr1m2qgmXbAUFBfrll1+qqmppaanOmzdP582b53GqioXLvAvEslVNdbIBSzTAOvWwzUSqug1IA94GknHOJLq0slXIPVNoMrBGVf9dzmjvAiPFcSqw352uMWFjyZIlnHzyyZxxxhls2+Ysnv3797dmIVMnVOoIl6qW4BSDt6swjV7A5cAqEVnuPnc3cLT72hOBOTinla7FObX0sKesGhMqeXl53HvvvTzxxBOkpKSQnp5uHcuZOifopzuocxygwk0nd9dlTLCzGHOk8vLy6Nq1Kz/++COjRo3iscceo0mTJl7HMqbG2blvxgRQWFhIdHQ0cXFxXHXVVQeah4ypq6xvImP8vP/++xx77LEsXLgQVeXOO++0QmDqPCsGxrh2797NZZddxjnnnENCQgKxsbF2cNjUG1YMjAGmT5/OCSecwLRp07j33ntZtmwZaWlpXscyJmTsmIExwMaNGznmmGOYPHkyJ554otdxjAk52zMw9ZKq8vzzz/P2287Z0jfddBNffvmlFQJTb1kxMPXOunXrOPPMMxk9ejTTpk0DICIiwjqWM/WaFQNTb5SUlPDvf/+bE088kaVLlzJp0qQDxcCY+s42hUy98b///Y9bb72Vc889l+eee0ilpmMAABmtSURBVI7Wra1jXGPK2J6BqdMKCwsPXIh+4MCBzJ8/n1mzZlkhMMaPFQNTZy1evJi0tDQGDBhwoGO5fv362W8HjAnAioGpc3Jzc7n11lvp2bMn+/bt46233qJVq1ZWBIypgB0zMHVKbm4uJ510Ej///DPXXnstjz76KI0bN/Y6ljFhz4qBqRMKCgqIiYkhPj6eUaNG0aNHD/r16+d1LGNqDWsmMrXee++9R4cOHfj8889RVW6//XYrBMYcISsGptbatWsXl1xyCUOGDKFJkyY0bNgQEbFjA8ZUgRUDUyu9+eabnHDCCbz99tvcf//9LF26lG7dunkdy5haK+jHDETkReAcYKeqdg4wvB8wC1jvPjVDVR8Idi5Tu23evJkOHTrw4osvkpqa6nUcY4KqpFTJyC1kb04hu7MLySzQGp9GKA4gvwxMAKZUMM7nqnpOCLKYWqq0tJSJEyeSnJzMsGHDuOmmm7jpppuIiIjwOpoxR6SopJTMvCIy84vZn1dEZl4RWfnFZOYXsT+viN1ZBezOLmB3duGBv3tzCij1Wf//pUsMQ2o4VyiugfyZiLQP9nRM3fXzzz9zyy23sGLFCkaMGMEFF1xgRcB4rqiklH25hezLKWJfbqG75V5ERl4h+3OdFXtGrvM4w328P6+I3MKSCl83NqoByY1iaJ4QQ9tm8XQ9uglJDWNonhBNUkIMSQnR7Fm3qsbfjzjXog8utxjMrqCZ6G1gM7AVuE1VvyvndUYDowFSUlLS0tPTA04vOzubhISEmohe4yxb5ZWUlPDWW2/x0ksvERUVxfXXX8/ZZ58ddgeIw22++QvnfOGUrVSVnCLILFD2Fyo7M/MokBiyCpXMQiXL55ZZqOQVl/9akQ0gIUqIj3L+Niy7RUJ8lBDvPo6PhLhI53G8Oyw2gsMu49WZb/3791+qqt39nw+HYpAIlKpqtoicDTylqsce7jW7d++uS5YsCThswYIFYXtqoWWrvA8//JDBgwczdOhQLrvsMoYPH+51pIDCbb75C+d8wchWUFzibJHnlm2xO00xGXmFB7bO9+cVu80zTnNNZl4Re3MKKS49dH0Y0UBIahhNs4bRJCVE06xhDEkNo2kaH02zhlE0iXeGNY2PpmnDKJrERRMXHdw91+rMNxEJWAw8/9GZqmb63J8jIv8RkeaqutvLXMYbBQUFLF68mD59+jBo0CA+//xzevXqxaeffup1NOORopJS9rjt57vc9vQ9OYXsyS5gT3Yhe90Dq3tzCtmXU0hOBc0wEQ2ExnFRJMZGkhgXRWJsFCmJsSTGRpGUEH2geaZ5Qgzr16xgcP/eNI6LokGD8NobDQbPi4GItAR2qKqKSA+c0133eBzLeODrr7/m6quvZv369axfv56WLVvSu3dvr2OZGlZUUnpgC33tvhL0h51k5Dlt77uzC9iRWcDOrHx2ZhawK7uAvTmFAV8nJrIBzRNiaOZutf8+OeGgrfWm8dE0iY+icVwUTeKd5xpGR1S6mbFgUwOaNoyuybce1kJxaukbQD+guYhsBu4FogBUdSIwHLheRIqBPOBiDUXblQkbOTk5/O1vf+Opp56iTZs2zJgxg5YtW3ody1RSbmExu7MK2Z1TwO4sZ+XtbLkXOgdYc50t9r25hWTkFJFV4NfYvuibA3ejIoTkhBiSE2NplxRP9/ZNadEoluaNommeEONsuTd0DqLGH8GK3RxeKM4muuQwwyfgnHpq6qGyjuXWrl3L9ddfzyOPPEJiYqLXsQyQX1TCrqwCtmfms21/Pjv257MjM58dWQXsyMxnV1YBOzPzy22WiY+OOLDV3iQ+mmOaNzxoi71JfBQbfl5D31O70yTO2YKvL00y4cjzZiJTP/l2LHfttdfSo0cP+vbt63WsOquopJS9OYXsOrDl7rS3+x5kLWu62ee2wQc6BTI2qgEtE2Np0SiW1KMS6dcx+cAWe7Lb1u4cZI0mNurwB1EXZPzMSW2bBOMtmyNkxcCE3MyZMxk7dixvvPEGvXv35tZbb7Xd/SpSVTJyi9iSkcf2/flsy8xnW0beQe3uW/fmkPXhBwH/v4FA0/hoGrtt680TounQIuHAFn3zhGhaNo6jVeNY90BrpH1WdZQVAxMyO3bs4IYbbuCtt96iS5cuNGrUyFYsARQWlx7Yct/jniHjtL0XkZFb+NuZNdkFbMvIJ6/o4C34yAZCcqMYWiTGcnRSPEdF59Gl4zEHzpJp7m65JzWMoVFspDXLGMCKgQmR9PR0xowZQ3Z2Ng899BC33347UVFRXscKucz8Irbsy2Pb/jy2ZDhb8dsz3bb4TOe0yYzconL/v3Gccwpk84QYjm/ZiP4dW3BUkzhaN4mllbsFn5QQQ4TPCt45J/24ULw9U4tZMTAhsX37djp27MjkyZM54YQTvI5To1SV/blF7M4pIMPtniAjr4idWfls35/P1ox8tmbksXlfLpn5B59JE9lAaNEohpTGsXRITuC03yeRnBBzoNuBsuaapvHRNI6LOmglb0xNsmJggqK0tJTnnnuOFi1aMHz4cG644QZuuOGGWtenkKqyN6eQjXtz2ZKRx5Z9eWzNyHPOrskqYMf+fHZl5VPy0ccB/z8xNpJWjeM4qkks3ds3pU3TOHdL3vnb3G8r3hivWDEwNe7HH39k1KhRfPHFF1xyySUMHz48LIuAqpJVUMzurAL3DBrnR08b9uSycW8Ov+7OZePeXLL9zotvFBt54IBqh+Tm5O/bQdfUDgd+ANXkwMHYGBrG2FfM1A62pJoaU1RUxBNPPMF9991HXFwcL730EldccYVnB4lLSpWtGXls3Ous1Dfvy2VrRj5bMpw2+11ZBeQXlR7yf1ERQttm8bRrFk+PY5pxdLN4jm4WT5tmztZ8YuzBxzoWLFhAvz6/C9XbMiYorBiYGjNv3jzGjRvHsGHDmDBhQkh+RZxfVMLOzAK2ZOSxcW8OG/fm8uvuXNbuzGb97hwKS35b2Uc2EFISY2ndJI6ubZuSkhhDi0axJDeKoWnDaJrFR9MsIZqWibHWdGPqHSsGplry8/NZvHgxffv25ayzzmLhwoWcdtppNfb6JaXKln15rNxVzPqF69m4N5dNe52Dsdv257M/7+AzbyIbCG2axtGhRQL9OiZzTPOGHJ3kbNm3ahxnK3ljymHFwFTZl19+yTXXXMP69ev59ddfadmyZZUKQWmpsnW/e1A20zkD56cdWazZlsVPO7IoKHa37pd+T3x0BG2bxtO2WRwnt2/mbN0nxnJU4zjaJcXTqnEskRF2aW9jjpQVA3PEsrOzufvuu5kwYQJt27Zl1qxZlWoSKiopZe3ObH7emc3aHVn8vDObX3blsH5PDoXFB7fdJzWM5oRWiVx+ajuOS2lExuafOH9Ab5onRNsP1YwJAisG5ojk5ubSpUsX1q9fz5gxY3j44Ydp1KjRQeOUlKrbdp/Dr3tyWLcrm9VbMvl+W+aBlX4DgXZJDfl9ckNO75hM+6SGtG4aR8vEWFomxpIYd3C3Bwty1pHcKCak79WY+sSKgamUvLw84uLiiI+PZ8yYMZxyyin06tULgK0ZeSzdsI/lmzJYtXk/q7fuP6iTs0YxkaQelcjIU9txYpvGdGzZiGOaNyQmMvxONzWmvrJiYA7r7bff5oYbbuCNN97g1F696X3+lSzftI/XXlvKsg0ZbM/MB5weLTsd1ZiLurel01GJ/C65Ie2SGpLU0Jp2jAl3VgxMubZt28bYsWOZ+f5HHH/mhfx3dQl/+WTugf7r2zaL4+RjmtG9XVO6Hd2U41s1IsoO3hpTK4XiSmcvAucAO1W1c4DhAjwFnA3kAleq6rJg5zIVe2HK69z9nzeJaN+ddjdfRQ7Cuv2lnNe1Nf06tqDr0U1onmBt+MbUFaHYM3gZ50pmU8oZPhg41r2dAjzn/jUeyC8q4aWFvzJ+TQLx/f6PVo0iGdqtHYM6t+QPrRtbd8fG1FGhuOzlZyLSvoJRhgJT3Osefy0iTUSklapuC3Y285vi4mKenvUVN3+SRUZRA848viU3nHEsXdo2sfZ+Y+oBCcW1591iMLucZqLZwCOq+oX7+BPgTlVdEmDc0cBogJSUlLT09PSA08vOziYhIaHG8tekcMy27OdNPPv1HkqSjyUmZzs392vPCUnhdaZPOM63MuGcDcI7n2Wrmupk69+//1JV7e7/fDgcQA602RmwQqnqJGASQPfu3bVfv34BX9C5mEfgYV4Lp2y5+QVc8fAUFucmQeNYesg60p8cS0QYHgQOp/nmL5yzQXjns2xVE4xs4VAMNgNtfR63AbZ6lKXeWL1lP9e/9AWbCo+iWcEGptx8Lns2rQvLQmCMCb5w+Oa/C4wUx6nAfjteEDz7MrMZ89+PGfrsQvIlhlt6JPDtM3/hxN+3Pfw/G2PqrFCcWvoG0A9oLiKbgXuBKABVnQjMwTmtdC3OqaVXBTtTffXCrPk89PGvaKMWnNupOQ8O60bj+Pp3HWJjzKFCcTbRJYcZrsCYYOeoz7bv3seIh99gQ3Q7kAhu6RrFjSPs7F1jzG/C4ZiBCaIPV2ziuhcWoPFt6VC8kfR//JnkpolexzLGhBkrBnXU1j2ZPP6/X5jx7RaaNW7Izac1YeTZ53odyxgTpqwY1DGlpaXc8/w7vPp9IZHxidx4RgfGnNHBegg1xlTIikEdsuKnDYx8ajb7G7UnqmgP/x58HEP6dvQ6ljGmFrBiUAeoKv/vmTd5Z71AXCtOi9/BS/dfSWxMtNfRjDG1hBWDWm7jnlzGvbOShVsTiM3fyMRr+tC/+yG9fhhjTIWsGNRShUXFjPrXq3yV05yYqCgeHNqJi08eRKQdGzDGVIEVg1rog4XfcuOriylq3IZm+Vt5//7LaNU4zutYxphazIpBLZKdm8/If05laX4yxDRmRNtcHrn+Gho0CIdeRYwxtZkVg1pixaYM/vLyQrYUHUVy/q+8ftv5HNfuKK9jGWPqCCsGYW7P/izufn0hczeWkNwohjt6xvKXodZ7hzGmZlkxCGP/nTGPf87bCAnJDO3cnH8M70ZirHUsZ4ypeVYMwtDWXXsZ8fAbbIppDwi3p0Uz5kLrWM4YEzxWDMLMB8s3cv2Ln6GxbTmudCPpD11GUuNGXscyxtRxVgzCxOZdGTz2yXpmLd9KUmI8t/RO4s+DhngdyxhTT1gx8FhpaSl3T3ybN34sITK+ETcPOJa/9OtAdKSdLmqMCZ2QrHFEZJCI/Cgia0XkrgDDrxSRXSKy3L2NCkUur337w3q6jH2O9I3xRBVl8vSfjuLmAcdZITDGhFwoLnsZATwL/BHYDHwjIu+q6vd+o05T1bHBzhMOVJUbn5rGuxsiIK4VfRJ28uI/riY6ynbUjDHeCMXapwewVlV/ARCRdGAo4F8M6oUdOaVc+vwivtreiLiCjTz/f/3o0/UEr2MZY+q5UBSD1sAmn8ebgUDnSQ4Tkb7AT8D/U9VNAcaptQqLirn60al8lZ1EfEwR/zy/MyNOHmxdSRhjwoI416MP4gRELgTOUtVR7uPLgR6qeoPPOElAtqoWiMh1wEWqekaA1xoNjAZISUlJS09PDzjN7OxsEhISav7NVNGiHzYy8dsctOnRxO9bx0Pnd6ZpbPgVgXCbb74sW9WFcz7LVjXVyda/f/+lqtr9kAGqGtQb0BP4yOfxOGBcBeNHAPsP97ppaWlanvnz55c7LJSycvJ06F8n6dG3vaNH3/i6jntuun7yySdexypXuMy3QCxb1YVzPstWNdXJBizRAOvUUGyefgMcKyLHiEg0cDHwru8IItLK5+EQYE0IcgXV8k0ZDHziE5YXH0WL/C18cms/Hr5umDULGWPCUtCPGahqsYiMBT7C2ep/UVW/E5EHcCrUu8CNIjIEKAb2AlcGO1ew7NmfxV2vLeSTTSWkJMZyd+9YRp9jHcsZY8JbSM5lVNU5wBy/5/7uc38cTvNRrfaf6f/jXwu2QEJzzj8xmQeGdaWRdSxnjKkF7MT2GrB5xx5G/DOdLbHtAWVcjziuvaCH17GMMabSrBhU0/vfbuQvL30OsW05Xjfx+kN/to7ljDG1jhWDKtq4cx+PffIr763YSnLjeG7rk8zFA61jOWNM7WTF4AiVlpZyx3+m89baUiLjGnHLH4/jutN/b/0JGWNqNSsGR2DJ9+u4csJHZCe2I7pwM89cdDyDTjvW61jGGFNtVgwqobRUGfvkNN7fHAFxKfRL3MXz/7jGOpYzxtQZtjY7jF92ZXPXjFUs3tmI+PwNTL72TE7r0tHrWMYYU6OsGJQjv6CQqx6Zyjf5KcTFRPHIBZ25qLt1LGeMqZusGAQwc8E33PbmtxQntia5YCuzx11OSmKs17GMMSZorBj4yMrJ49KHprCyqCVEJXD57wp4YNQ1tjdgjKnzrBi4lm7Yx9hXvmRbaRtaFvzK63cM4/dtWnodyxhjQqLeF4Od+/Yz7rUvmbellFaJsfytbzyjzv6T17GMMSak6nUxeHrax/z7822Q0JzhXZK574JuJMTU61lijKmn6uWab8O2XVz8yJtsi2uPaCl/7RnPqKHWsZwxpv6qd8XgvWUbGPvKQohpSyc28cYjl9O4UUOvYxljjKfqXTHYnVdKSmIst/dL4cIzrWM5Y4yBelgMRvZsz2WntiMqwk4XNcaYMiFZI4rIIBH5UUTWishdAYbHiMg0d/giEWkfrCwRDcQKgTHG+An6WlFEIoBngcFAKnCJiKT6jXYNsE9VOwDjgUeDncsYY8xvQrGJ3ANYq6q/qGohkA4M9RtnKPCKe386cKaISAiyGWOMAURVgzsBkeHAIFUd5T6+HDhFVcf6jLPaHWez+3idO85uv9caDYwGSElJSUtPTw84zezsbBISEoLxdqrNslWNZau6cM5n2aqmOtn69++/VFW7HzJAVYN6Ay4EXvB5fDnwjN843wFtfB6vA5Iqet20tDQtz/z588sd5jXLVjWWrerCOZ9lq5rqZAOWaIB1aiiaiTYDbX0etwG2ljeOiEQCjYG9IchmjDGG0Bwz+AY4VkSOEZFo4GLgXb9x3gWucO8PB+a5FcwYY0wIBP13BqpaLCJjgY+ACOBFVf1ORB7A2V15F5gMTBWRtTh7BBcHO5cxxpjfBP0AcrCIyC5gQzmDmwO7yxnmNctWNZat6sI5n2Wrmupka6eqyf5P1tpiUBERWaKBjpaHActWNZat6sI5n2WrmmBks5/iGmOMsWJgjDGm7haDSV4HqIBlqxrLVnXhnM+yVU2NZ6uTxwyMMcYcmbq6Z2CMMeYIWDEwxhhTu4tBOF0noQrZrhSRXSKy3L2NCmG2F0Vkp9tBYKDhIiJPu9lXiki3MMrWT0T2+8y3v4coV1sRmS8ia0TkOxG5KcA4nsy3Smbzar7FishiEVnhZrs/wDiefE8rmc2z76k7/QgR+VZEZgcYVrPzLVCHRbXhhvNr5nXA74BoYAWQ6jfOX4CJ7v2LgWlhlO1KYIJH864v0A1YXc7ws4EPAAFOBRaFUbZ+wGwP5lkroJt7vxHwU4DP1JP5VslsXs03ARLc+1HAIuBUv3G8+p5WJptn31N3+rcArwf67Gp6vtXmPYNwvk5CZbJ5RlU/o+KOAIcCU9TxNdBERFqFSTZPqOo2VV3m3s8C1gCt/UbzZL5VMpsn3HmR7T6Mcm/+Z6148j2tZDbPiEgb4E/AC+WMUqPzrTYXg9bAJp/Hmzn0C3BgHFUtBvYDSWGSDWCY25wwXUTaBhjulcrm90pPd9f+AxHpFOqJu7vjXXG2JH15Pt8qyAYezTe3qWM5sBOYq6rlzrcQf08rkw28+54+CdwBlJYzvEbnW20uBoEqoH9Vr8w4wVCZ6b4HtFfVPwD/47cKHw68mm+VsQynb5UuwDPAzFBOXEQSgLeBm1U1039wgH8J2Xw7TDbP5puqlqjqSTjd1/cQkc5+o3g23yqRzZPvqYicA+xU1aUVjRbguSrPt9pcDML5OgmHzaaqe1S1wH34PJAWglyVVZl56wlVzSzbtVfVOUCUiDQPxbRFJApnZfuaqs4IMIpn8+1w2bycbz4ZMoAFwCC/QZ5fz6S8bB5+T3sBQ0TkV5xm5jNE5FW/cWp0vtXmYhDO10k4bDa/tuQhOO284eJdYKR7dsypwH5V3eZ1KAARaVnWLioiPXCW4T0hmK7gdLW+RlX/Xc5onsy3ymTzcL4li0gT934cMAD4wW80T76nlcnm1fdUVcepahtVbY+z/pinqpf5jVaj8y3o1zMIFg3j6yRUMtuNIjIEKHazXRmKbAAi8gbO2SXNRWQzcC/OwTNUdSIwB+fMmLVALnBVGGUbDlwvIsVAHnBxiAp8L5xLtq5y25gB7gaO9snm1XyrTDav5lsr4BURicApQG+q6uxw+J5WMptn39NAgjnfrDsKY4wxtbqZyBhjTA2xYmCMMcaKgTHGGCsGxhhjsGJgjDEGKwbGGGOwYmBMtbg/LPzUvd9NRFREktw+b1aLSLzXGY2pDCsGxlRPBk630QA3AF8DTXF6m5yrqrki0tSrcMZUlhUDY6pnPxAvIkk4v2hdiFMMRgPPueOM9yibMZVmxcCYalDVsu6F/w+ne4As4A9AhKr+JCKDgONF5DavMhpTGVYMjKm+UpxOzN4BMoHbgInusN3Aq6r6uEfZjKkUKwbGVF8h8IF7gZFMoCFQds3aP+Bc9tSYsGYd1RkTRG6Pl8OAR1Q1nLopN+YgVgyMMcZYM5ExxhgrBsYYY7BiYIwxBisGxhhjsGJgjDEGKwbGGGOwYmCMMQYrBsYYY7BiYIwxBvj/nsnNzK0fkT8AAAAASUVORK5CYII=\n",
                        "text/plain": [
                            "<Figure size 432x288 with 1 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "wvec = np.linspace(wmin, wmax, 100)\n",
                "ζvec, hvec = dr(wvec*0, wvec*0, wvec*0, wvec*0, wvec)\n",
                "\n",
                "plt.title(\"Multidimensional Consumption-Savings (decision rule)\")\n",
                "plt.plot(wvec, wvec, linestyle='--', color='black')\n",
                "plt.plot(wvec, wvec*ζvec)\n",
                "plt.xlabel(\"$w_t$\")\n",
                "plt.ylabel(\"$c_t$\")\n",
                "plt.grid()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "```\n  ___ ___                        _____           ._.\n /   |   \\_____ ___  __ ____   _/ ____\\_ __  ____| |\n/    ~    \\__  \\\\  \\/ // __ \\  \\   __\\  |  \\/    \\ |\n\\    Y    // __ \\\\   /\\  ___/   |  | |  |  /   |  \\|\n \\___|_  /(____  /\\_/  \\___  >  |__| |____/|___|  /_\n       \\/      \\/          \\/                   \\/\\/\n```"
        }
    ],
    "metadata": {
        "jupytext": {
            "cell_metadata_filter": "-all",
            "encoding": "# -*- coding: utf-8 -*-",
            "main_language": "python",
            "notebook_metadata_filter": "-all"
        },
        "kernelspec": {
            "display_name": "Python 3",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.7.4"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 4
}