expansion_suite.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5e7b1d0",
   "metadata": {
    "id": "f5e7b1d0"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/lphansen/RiskUncertaintyValue/blob/main/expansion_suite.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7735eea6",
   "metadata": {},
   "source": [
    "# Numerical Solution: Expansion Suite"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef624dcb",
   "metadata": {
    "id": "ef624dcb"
   },
   "source": [
    "We previously covered the small-noise expansion method and the recursive utility framework. This section covers the usage of `Expansion Suite`, an open source Python toolbox which uses these methods to solve discrete nonlinear DSGE models.\n",
    "\n",
    "\n",
    "\n",
    "<!-- For this notebook \n",
    "\n",
    "- Section 1 provides a five-minute guide to solving the DSGE model using the expansion suite, as well as how to \n",
    "    - compute shock elasticities and IRF\n",
    "    - approximate and simulate variables based on model solutions. \n",
    "\n",
    "\n",
    "- Section 2 provides examples of `LinQuadVar` computations, including\n",
    "    - constructions, operations, and evaluations\n",
    "    - adjusting time periods\n",
    "    - computing expectations -->\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f590ad3-c8d4-4fba-93f4-bf14f4f3364b",
   "metadata": {},
   "source": [
    "<!-- **Table of Contents**\n",
    "* [1. Five-minute guide to solving the DSGE model using the Expansion Suite](#s1)\n",
    "    * [1.1 Equilibrium Condition Function](#s1-1)\n",
    "    * [1.2 Model Steady State Function](#s1-2)\n",
    "    * [1.3 Model Variable Dimension Array](#s1-3)\n",
    "    * [1.4 Model Parameter Calibration Tuple](#s1-4)\n",
    "    * [1.5 Log Consumption Growth Function](#s1-5)\n",
    "    * [1.6 Solve the BY model using the Expansion Suite](#s1-6)\n",
    "    * [1.7 Model Solutions](#s1-7)\n",
    "    * [1.8 Approximate Variables](#s1-8)\n",
    "    * [1.9 Model Simulations](#s1-9)\n",
    "    * [1.10 Shock Elasticities and IRF](#s1-10)\n",
    "* [2. Examples to use `LinQuadVar` in Computation](#s2)\n",
    "    * [2.1 Define a `LinQuadVar` ](#s2-1)\n",
    "    * [2.2 `LinQuadVar` Operations](#s2-2)\n",
    "    * [2.3 `LinQuadVar` Split and Concat](#s2-3)\n",
    "    * [2.4 Evaluate a `LinQuadVar`](#s2-4)\n",
    "    * [2.5 Adjust Conditional Information Time](#s2-5)\n",
    "    * [2.6 Compute the Expectation of time $t+1$](#s2-6) -->\n",
    "\n",
    "\n",
    "```{contents}\n",
    ":depth: 1\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e247ca61",
   "metadata": {},
   "source": [
    "Consider the framework previously introduced:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}[{N_{t+1}^*Q_{t+1}^* \\psi_1(X_t, J_t, X_{t+1}, J_{t+1})}|{\\mathfrak{A}_t}] - \\psi_2(X_t, J_t)  &= 0 \\tag{1}\\\\\n",
    "  \\phi\\left(X_t, X_{t+1}, J_t, {\\sf q}W_{t+1}, {\\sf q} \\right) &= 0 \n",
    "\\tag{2} \\end{align}\n",
    "$$\n",
    "\n",
    "The `Expansion Suite` uses the function `uncertain_expansion` to approximate a solution to the above system locally. It takes the following inputs:\n",
    "\n",
    "```{list-table}\n",
    ":header-rows: 1\n",
    "\n",
    "* - Input\n",
    "  - Type\n",
    "  - Description\n",
    "* - `eq`\n",
    "  - *callable*\n",
    "  - Function which takes state and jump variables and outputs error of equilibrium conditions, i.e. RHS of equation (1)\n",
    "* - `ss`\n",
    "  - *callable*\n",
    "  - Function which takes model parameters and outputs steady-state values of state and jump variables\n",
    "* - `var_shape`\n",
    "  - *tuple of int*\n",
    "  - Triple of integers $(n_J,n_X,n_W)$ where $n_J$ is the number of jump variables, $n_X$ is the number of state variables, and $n_W$ is the number of shocks.\n",
    "* - `args`\n",
    "  - *tuple of float or ndarray*\n",
    "  - Array of preference and model parameters\n",
    "* - `gc`\n",
    "  - *callable*\n",
    "  - Function which takes in variables and outputs log-growth of consumption\n",
    "* - `args`\n",
    "  - tuple of float/ndarray\n",
    "  - Preference and model parameters\n",
    "* - `gc`\n",
    "  - *callable*\n",
    "  - Function which takes in variables and outputs log-growth of consumption\n",
    "* - `approach`\n",
    "  - *string*\n",
    "  - If equal to '1', then solve the model using approach 1 (default, see Appendix C); If equal to '2', then use approach 2 (see Appendix D). The function raises an exception if other values are provided\n",
    "* - `init_util`\n",
    "  - *dict, or None*\n",
    "  - Dictionary for initialization of the expansion coefficients $\\mu^0, \\Upsilon_0^2, \\Upsilon_1^2,$ and $\\Upsilon_2^2$. The keys are: `mu_0`, `Upsilon_0`, `Upsilon_1`, `Upsilon_2`. If *None*, zero matrices are used\n",
    "* - `iter_tol`\n",
    "  - *float*\n",
    "  - Iteration stops when the maximum absolute difference between the approximated variables in the current and previous iteration is below `iter_tol`\n",
    "* - `max_iter`\n",
    "  - *int*\n",
    "  - Maximum number of iterations\n",
    "```\n",
    "\n",
    "The output is of class `ModelSolution`, which stores the coefficients for the linear-quadratic approximation for the jump and state variables; continuation values; consumption growth; and log change of measure, as well as the steady-state values of each variables. \n",
    "\n",
    "We will now walk through the computation using a simple example.\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f6c2319",
   "metadata": {},
   "source": [
    "# 1 AK Model with Recursive Utility"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1141ee3e",
   "metadata": {},
   "source": [
    "## 1.1 Model Setup <a class=\"anchor\" id=\"s1-1\"></a>\n",
    "\n",
    "Consider the following recursive utility maximization problem with AK technology:\n",
    "\n",
    "$$\n",
    "\\max_{C_t, I_t} V_t = \\left[ (1 - \\beta) \\left(C_t\\right)^{1-\\rho}  + \\beta  \\left( R_t \\right)^{1-\\rho} \\right]^{\\frac 1 {1-\\rho}}\\tag{3}\n",
    "$$\n",
    "such that\n",
    "\n",
    "$$\n",
    "C_t + I_t = \\alpha K_t\n",
    "$$ \n",
    "\n",
    "$$\n",
    "\\log K_{t+1} - \\log K_t = \\frac 1 \\psi \\log \\left[ 1 + \\psi \\frac {I_t}{K_t} \\right] + \\beta_k Z_t - \\eta_k + \\sigma_k W_{t+1} \n",
    "$$ \n",
    "\n",
    "$$\n",
    "Z_{t+1} = \\exp( - \\beta_z) Z_t + \\sigma_z W_{t+1} \n",
    "$$ \n",
    "\n",
    "where $\\theta_1$ and $\\theta_2$ are adjustment cost parameters which determine the efficiency of converting investment to capital, and $\\alpha_k$ is depreciation.\n",
    "$Z_t$ governs the (scalar) conditional mean of (stochastic) technology growth while $\\sigma_k$ governs the stochastic volatility of technology growth. \n",
    "$W_{t+1}$ is a two-dimensional Wiener process which provides the exogenous shocks to this model.\n",
    "\n",
    "For the purposes of computation, we divide the first two equations by $K_t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 First-Order Conditions <a class=\"anchor\" id=\"s1-2\"></a> \n",
    "\n",
    "Let ${\\widehat G}_t = {\\widehat K}_t,$ $X_t = Z_t,$ and $D_t = \\frac {I_t }{K_t}.$ We can rewrite the previous equations as: \n",
    "\n",
    "$$\n",
    "{\\widehat C}_t = \\log \\left( \\alpha - D_t  \\right) + {\\widehat G}_t \n",
    "$$ \n",
    "\n",
    "$$\n",
    "{\\widehat G}_{t + 1} - {\\widehat G}_t =  \\frac 1 \\psi \\log \\left[ 1 + \\psi D_t \\right] + \\beta_k X_t - \\eta_k + \\sigma_k W_{t+1}\n",
    "$$ \n",
    "\n",
    "Apply the general form of FOC and costate equations and notice in this special case $\\psi^x_d = 0$, we have \n",
    "\n",
    "The FOC is  \n",
    "$$\n",
    "    \\mathbb{E} \\left[N^*_{t + 1}Q^*_{t + 1} \\beta  \\left(\\frac{C_{t + 1}}{G_{t + 1}}\\right)^{-\\rho} \\left(\\frac{G_{t + 1}}{G_{t}}\\right)^{1-\\rho} \\left(\\frac{C_{t}}{G_{t}}\\right)^{\\rho} \\frac{1}{1 + \\psi D_t} \\widetilde{MG}_{t+1} \\right] - 1 = 0\n",
    "$$ \n",
    "\n",
    "The law of motion of costate $\\widetilde{MG}$ is \n",
    "\n",
    "$$\n",
    "    \\mathbb{E} \\left[N^*_{t + 1}Q^*_{t + 1} \\beta \\left(\\frac{C_{t + 1}}{G_{t + 1}}\\right)^{-\\rho} \\left(\\frac{G_{t + 1}}{G_{t}}\\right)^{1-\\rho} \\left(\\frac{C_{t}}{G_{t}}\\right)^{\\rho} \\widetilde{MG}_{t+1} \\right ] - \\widetilde{MG}_t + \\frac{C_t}{G_t} = 0\n",
    "$$\n",
    "\n",
    "Finally, we can express our system of equations as\n",
    "\n",
    "$$\n",
    "\\psi_1 = \n",
    "\\begin{bmatrix}\n",
    "0\\\\ \n",
    "\\beta  \\left(\\frac{C_{t + 1}}{G_{t + 1}}\\right)^{-\\rho} \\left(\\frac{G_{t + 1}}{G_{t}}\\right)^{1-\\rho} \\left(\\frac{C_{t}}{G_{t}}\\right)^{\\rho} \\frac{1}{1 + \\psi D_t} \\widetilde{MG}_{t+1} \\\\\n",
    "\\beta \\left(\\frac{C_{t + 1}}{G_{t + 1}}\\right)^{-\\rho} \\left(\\frac{G_{t + 1}}{G_{t}}\\right)^{1-\\rho} \\left(\\frac{C_{t}}{G_{t}}\\right)^{\\rho} \\widetilde{MG}_{t+1}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\psi_2 = \n",
    "\\begin{bmatrix}\n",
    "\\frac{C_t}{G_t} + D_t - {\\mathbf a}\\\\\n",
    "1\\\\\n",
    "\\widetilde{MG}_t - \\frac{C_t}{G_t}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\\phi=\\begin{bmatrix} \n",
    "&\\frac{G_{t+1}} {G_t}  -  \\left( 1 + \\psi D_t \\right)^{\\frac{1}{\\psi}} + exp( \\beta_k X_t - \\eta_k + \\sigma_k W_{t+1} ) \\\\\n",
    "&X_{t+1} - \\exp( - \\beta_z) X_t - \\sigma_z W_{t+1} \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "\n",
    "The model has thre\n",
    "e jump variables (In the expansion suite, $J_t$ is used to represent the array of jump variables):\n",
    "- $\\log C_t - \\log G_t$ : Log consumption over capital\n",
    "- $\\log D_t$ : Log investment over capital\n",
    "- $\\log MG_t$ : Log co-state variable\n",
    "\n",
    "and two state variables (In the expansion suite, $X_t$ is used to represent state variables):\n",
    "- $\\log G_{t} - \\log G_{t-1}$ : Log capital growth process\n",
    "- $X_{t}$ : Exogenous technology"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36420e45",
   "metadata": {},
   "source": [
    "## 1.2 First-Order Conditions <a class=\"anchor\" id=\"s1-2\"></a>\n",
    "\n",
    "In general, write a triangular system with stochastic growth as:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "X_{t+1}(\\mathbf{q}) & =\\psi^x\\left[D_t(\\mathbf{q}), X_t(\\mathbf{q}), \\mathbf{q} W_{t+1}\\right] \\\\\n",
    "\\log G_{t+1}(\\mathbf{q})-\\log G_t(\\mathbf{q}) & =\\psi^g\\left[D_t(\\mathbf{q}), X_t(\\mathbf{q}), \\mathrm{q} W_{t+1}\\right]\n",
    "\\end{aligned}\n",
    "$$\n",
    "where $D_t$ is a date $t$ decision vector for the planner, $G_t$ is a growth variable, $X_t$ is the state vector. In addition, we impose\n",
    "$$\n",
    "\\widehat{C}_t(\\mathrm{q})=\\kappa\\left[D_t(\\mathrm{q}), X_t(\\mathrm{q})\\right]+\\widehat{G}_t(\\mathrm{q}) \\text {. }\n",
    "$$\n",
    "\n",
    "In what follows for notational convenience, we will leave the q implicit. We use homogeneity to rewrite the utility recursion:\n",
    "$$\n",
    "\\widehat{V}_t-\\widehat{G}_t=\\frac{1}{1-\\rho} \\log \\left[(1-\\beta) \\exp \\left[(1-\\rho)\\left(\\widehat{C}_t-\\widehat{G}_t\\right)\\right]+\\beta \\exp \\left[(1-\\rho)\\left(\\widehat{R}_t-\\widehat{G}_t\\right)\\right]\\right]\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\widehat{R}_t-\\widehat{G}_t=\\frac{1}{1-\\gamma} \\log \\left(\\mathbb{E}\\left[\\exp \\left((1-\\gamma)\\left[\\left(\\widehat{V}_{t+1}-\\widehat{G}_{t+1}\\right)+\\left(\\widehat{G}_{t+1}-\\widehat{G}_t\\right)\\right]\\right) \\mid \\mathfrak{A}_t\\right]\\right)\n",
    "$$\n",
    "\n",
    "Differetiate the utility recursion with respect to $D_t$, we obtain the first-order condtions for $D_t$:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& (1-\\beta) \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{G}_t\\right)\\right]\\exp \\left[(1-\\rho)\\left(\\widehat{C}_t-\\widehat{G}_t\\right)\\right] \\kappa_d\\left(D_t, X_t\\right) \\\\\n",
    "& +\\beta \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{G}_t\\right)\\right] \\exp \\left[(1-\\rho)\\left(\\widehat{R}_t-\\widehat{G}_t\\right)\\right]\\mathbb{E}\\left(\\exp \\left[(1-\\gamma)\\left(\\widehat{V}_{t+1}-\\widehat{R}_t\\right)\\right] \\psi_{d^{\\prime}}^x\\left(D_t, X_t, W_{t+1}\\right) M X_{t+1} \\mid \\mathfrak{A}_t\\right) \\\\\n",
    "& \\beta \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{G}_t\\right)\\right] \\exp \\left[(1-\\rho)\\left(\\widehat{R}_t-\\widehat{G}_t\\right)\\right]\\mathbb{E}\\left(\\exp \\left[(1-\\gamma)\\left(\\widehat{V}_{t+1}-\\widehat{R}_t\\right)\\right] \\psi_{d^{\\prime}}^g\\left(D_t, X_t, W_{t+1}\\right) \\mid \\mathfrak{A}_t\\right) \\\\\n",
    "& =0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the co-state $MX_{t+1} = \\frac{\\partial (\\widehat{V}_{t+1}-\\widehat{G}_{t+1})}{\\partial X_{t+1}}$.\n",
    "\n",
    "Using the envelop theorem, the corresponding co-state equation is given by \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "MX_t = & (1-\\beta) \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{C}_t\\right)\\right] \\kappa_x\\left(D_t, X_t\\right) \\\\\n",
    "& +\\beta \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{R}_t\\right)\\right] \\mathbb{E}\\left(\\exp \\left[(1-\\gamma)\\left(\\widehat{V}_{t+1}-\\widehat{R}_t\\right)\\right] \\psi_{x^{\\prime}}^x\\left(D_t, X_t, W_{t+1}\\right) M X_{t+1} \\mid \\mathfrak{A}_t\\right) \\\\\n",
    "& \\beta \\exp \\left[(\\rho-1)\\left(\\widehat{V}_t-\\widehat{R}_t\\right)\\right]\\mathbb{E}\\left(\\exp \\left[(1-\\gamma)\\left(\\widehat{V}_{t+1}-\\widehat{R}_t\\right)\\right] \\psi_{x^{\\prime}}^g\\left(D_t, X_t, W_{t+1}\\right) \\mid \\mathfrak{A}_t\\right) \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "In the AK model with production, let $MK_{t+1}$ be the co-state variable for the capital law of motion (4) and  $MC_{t+1}$ be the co-state variable for the consumption-investment allocation (6). You can access a step-by-step derivation {download}`here<../documents/AKderive.pdf>`. \n",
    "\n",
    "The FOCs are:\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\\left[\\left(\\frac{S_{t+1}}{S_t}\\right) \\frac {MK_{t+1}}{MC_{t+1}}  \\left(\\theta_1 \\theta_2\\left[1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)\\right]^{\\theta_1-1} \\exp(B_{t+1} - B_t)\\right) \\mid \\mathfrak{A}_t\\right] - 1  = 0,\\tag{7}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}\\left[\\left(\\frac{S_{t+1}}{S_t}\\right) \\left(\\frac {MK_{t+1}}{MC_{t+1}} \\right) \\left[\\left( \\frac{K_{t+1}}{K_t} \\right) - \\left(\\theta_1 \\theta_2\\left[1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)\\right]^{\\theta_1-1}  \\exp(B_{t+1} - B_t)\\right)\\left(\\frac {I_t} {K_t} \\right)  \\right] \\mid {\\mathfrak A}_t \\right] \\\\ \n",
    " - \\frac{MK_t}{MC_t} + {\\mathbf a}  = 0, \\tag{8}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "<!-- $$\n",
    "\\begin{align*}\n",
    "& MC_t \\mathbb{E}\\left[\\left(\\frac{S_{t+1}}{S_t}\\right) \\left(\\frac {MK_{t+1}}{MC_{t+1}} \\right) \\left[\\left( \\frac{K_{t+1}}{K_t} \\right) - \\left(\\theta_1 \\theta_2\\left[1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)\\right]^{\\theta_1-1} \\exp(B_{t+1} - B_t)\\right)\\left(\\frac {I_t} {K_t} \\right)  \\right] \\mid {\\mathfrak A}_t \\right]\\cr &  - MK_t + {\\mathbf a} MC_t = 0\n",
    "\\end{align*}\n",
    "$$ -->\n",
    "\n",
    "<!-- Equations 4 and 5 are the law of motion equations, while equations 6, 7, and 8 are the forward-looking equations. --> \n",
    "\n",
    "<br>\n",
    "\n",
    "## 1.3 Adapting the Model <a class=\"anchor\" id=\"s1-3\"></a>\n",
    "We want to write our FOCs and law of motion equations in the form of equations 1 and 2.\n",
    "\n",
    "First note that the stochastic discount factor can be expressed as:\n",
    "\n",
    "$$\n",
    "\\frac{S_{t+1}}{S_t} = \\beta \\left(\\frac{V_{t+1}}{R_t}\\right)^{1-\\gamma} \\left(\\frac{V_{t+1}}{R_t}\\right)^{\\rho-1} \\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho}\n",
    ":= N_{t+1}^*Q_{t+1}^* \\beta\\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho}\n",
    "$$\n",
    "\n",
    "Hence, we can rewrite equations 7 and 8 as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbb{E}\\left[ N_{t+1}^*Q_{t+1}^* \\beta\\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho}\\frac {MK_{t+1}}{MC_{t+1}} \\frac{K_{t+1}}{K_t} \\frac{1}{1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)} \\mid \\mathfrak{A}_t\\right] - 1  &= 0\\\\\n",
    "\\mathbb{E}\\left[N_{t+1}^*Q_{t+1}^* \\beta\\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho} \\frac {MK_{t+1}}{MC_{t+1}} \\frac{K_{t+1}}{K_t} \\frac{1\n",
    "+(\\theta_2-1)\\frac{I_t}{K_t}}{1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)} \\mid \\mathfrak{A}_t\\right] - \\frac{MK_t}{MC_t}+{\\mathbf a}&= 0\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Finally, we can write express our five equations as:\n",
    "\n",
    "$$\n",
    "\\psi_1 = \n",
    "\\begin{bmatrix}\n",
    "0\\\\ \n",
    "\\beta\\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho}\\frac {MK_{t+1}}{MC_{t+1}} \\frac{K_{t+1}}{K_t} \\frac{1}{1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)}\\\\\n",
    "\\beta\\left(\\frac{C_{t+1}}{C_t}\\right)^{-\\rho} \\frac {MK_{t+1}}{MC_{t+1}} \\frac{K_{t+1}}{K_t} \\frac{1\n",
    "+(\\theta_2-1)\\frac{I_t}{K_t}}{1+\\theta_2\\left(\\frac{I_t}{K_t}\\right)}\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\psi_2 = \n",
    "\\begin{bmatrix}\n",
    "\\frac{C_t}{K_t} + \\frac{I_t}{K_t} - {\\mathbf a}\\\\\n",
    "1\\\\\n",
    "\\frac{MK_t}{MC_t}-{\\mathbf a}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\\phi=\\begin{bmatrix} \n",
    "&\\frac{K_{t+1}}{K_t} - \\left[1 + \\theta_2 \\left({\\frac {I_t} {K_t}}\\right) \\right]^{\\theta_1} \\exp( - \\alpha_k + Z_t - {\\frac 1 2} \\mid \\sigma_k \\mid^2  + \\sigma_k\\cdot W_{t+1})\\\\\n",
    "&Z_{t+1} - {\\mathbb A} Z_t + \\mathbb{B} W_{t+1}\n",
    "\\end{bmatrix}$$\n",
    "\n",
    " \n",
    "The model has three jump variables (In the expansion suite, $J_t$ is used to represent the array of jump variables):\n",
    "- $\\log C_t - \\log K_t$ : Log consumption over capital\n",
    "- $\\log I_t - \\log K_t$ : Log investment over capital\n",
    "- $\\log MK_t - \\log MC_t$ : Log co-state variable\n",
    "\n",
    "and two state variables (In the expansion suite, $X_t$ is used to represent state variables):\n",
    "- $\\log K_{t} - \\log K_{t-1}$ : Log capital growth process\n",
    "- $Z_{t}$ : Exogenous technology\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5e59370",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "# 2 Inputs <a class=\"anchor\" id=\"s2\"></a>\n",
    "\n",
    "## 2.1 Libraries <a class=\"anchor\" id=\"s2-3\"></a>\n",
    "\n",
    "Begin by installing the following libraries and downloading `RiskUncertaintyValue`, which contains the functions required to solve the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0d30383",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fatal: destination path 'RiskUncertaintyValue' already exists and is not an empty directory.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>.container { width:97% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "workdir = os.getcwd()\n",
    "!git clone https://github.com/lphansen/RiskUncertaintyValue \n",
    "workdir = os.getcwd() + '/RiskUncertaintyValue'             \n",
    "sys.path.insert(0, workdir+'/src')                        \n",
    "import numpy as np\n",
    "import autograd.numpy as anp\n",
    "from scipy import optimize\n",
    "np.set_printoptions(suppress=True)\n",
    "np.set_printoptions(linewidth=200)\n",
    "from IPython.display import display, HTML\n",
    "from BY_example_sol import disp\n",
    "display(HTML(\"<style>.container { width:97% !important; }</style>\"))\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "from lin_quad import LinQuadVar\n",
    "from lin_quad_util import E, concat, next_period, cal_E_ww, lq_sum, N_tilde_measure, E_N_tp1\n",
    "from uncertain_expansion import uncertain_expansion\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06a5ca90-4eec-48d3-862a-4a00d766b58e",
   "metadata": {},
   "source": [
    "## 2.2 Equilibrium Condition Function <a class=\"anchor\" id=\"s2-2\"></a>\n",
    "\n",
    "The first input for `uncertain_expansion` requires us to define the equilibrium conditions for the model. Therefore we define a function `eq` which takes in a list of jump and state variables as inputs, and outputs the equilibrium equations in the form of equations 1 and 2. The inputs required:\n",
    "\n",
    "```{list-table}\n",
    ":header-rows: 1\n",
    "\n",
    "* - Variable\n",
    "  - Type\n",
    "  - Description\n",
    "  - Corresponds to\n",
    "* - `Var_t`\n",
    "  - array-like\n",
    "  - Vector of jump and state variables in the current period\n",
    "  - $(X_t, J_t)$\n",
    "* - `Var_tp1`\n",
    "  - array-like\n",
    "  - Vector of jump and state variables in the next period\n",
    "  - $(X_{t+1}, J_{t+1})$\n",
    "* - `W_tp1`\n",
    "  - array-like\n",
    "  - Vector of shock variables in the next period\n",
    "  - $(W_{t+1})$\n",
    "* - `q`\n",
    "  - float\n",
    "  - Perturbation parameter\n",
    "  - $q$\n",
    "* - `mode`\n",
    "  - string\n",
    "  - By default, this function returns the LHS of equation (1) and (2). Set mode to 'psi1', 'psi2', or 'phi' to return the corresponding function in the equilibrium conditions\n",
    "  -  \n",
    "* - `args`\n",
    "  - tuple of float/ndarray\n",
    "  - Preference and model parameters\n",
    "  -  \n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba75562f",
   "metadata": {
    "id": "ba75562f"
   },
   "source": [
    "\n",
    "While writing out `eq`, ensure that: \n",
    "\n",
    "- The first components of `Var_t` and `Var_tp1` are fixed as `q_t` or `q_tp1`.\n",
    "- State variables follow jump variables.\n",
    "- The number of jump variables (except `q_t` and `q_tp1`) equals the number of forward-looking equilibrium conditions. \n",
    "- The number of state variables equals to the number of state evolution equations.\n",
    "- All values are expressed as type `float`. Otherwise, this may cause type errors.\n",
    "\n",
    "Note that the log change of measure is skipped in the specification of equilibrium conditions.\n",
    "\n",
    "For the example model, we write:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6ae5d2d3",
   "metadata": {
    "id": "6ae5d2d3",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def eq_onecap_2d(Var_t, Var_tp1, W_tp1, q, mode, *args):\n",
    "    \n",
    "    # Unpack model parameters\n",
    "    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k, sigma_z1, sigma_z2 = args\n",
    "    w1_tp1, w2_tp1 = W_tp1.ravel()\n",
    "\n",
    "    # Unpack model variables\n",
    "    q_t, cmk_t, imk_t, mk_t, Z_t, gk_t = Var_t.ravel()\n",
    "    q_tp1, cmk_tp1, imk_tp1, mk_tp1, Z_tp1, gk_tp1 = Var_tp1.ravel()\n",
    "\n",
    "    # Intermedite varibles that facilitates computation\n",
    "    sdf_ex = anp.log(β) - ρ*(cmk_tp1+gk_tp1-cmk_t)\n",
    "    g_dep = -α_k + beta_k*Z_t + sigma_k*w1_tp1 - 0.5 * (sigma_k**2)\n",
    "\n",
    "    ## Forward-looking conditions\n",
    "    psi1_1 = 0.\n",
    "    psi1_2 = anp.exp(sdf_ex + gk_tp1) * 1/(1 + ϕ_2*anp.exp(imk_t)) * mk_tp1\n",
    "    psi1_3 = anp.exp(sdf_ex + gk_tp1) * mk_tp1\n",
    "\n",
    "    psi2_1 = -a + anp.exp(cmk_t) + anp.exp(imk_t)\n",
    "    psi2_2 = 1.\n",
    "    psi2_3 = mk_t - anp.exp(cmk_t)\n",
    "\n",
    "    # State evoluion processes\n",
    "    phi_1 = Z_tp1 - beta_z*Z_t - sigma_z1*w1_tp1 - sigma_z2*w2_tp1\n",
    "    phi_2 = gk_tp1 - ϕ_1 * anp.log(1.+ϕ_2*anp.exp(imk_t)) - g_dep\n",
    "    \n",
    "    if mode == 'psi1':\n",
    "        return np.array([psi1_1, psi1_2, psi1_3])\n",
    "    \n",
    "    return anp.array([\n",
    "        psi1_1 * anp.exp(q_tp1) - psi2_1,\n",
    "        psi1_2 * anp.exp(q_tp1) - psi2_2,\n",
    "        psi1_3 * anp.exp(q_tp1) - psi2_3,\n",
    "        phi_1, phi_2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d554abc5",
   "metadata": {
    "id": "d554abc5"
   },
   "source": [
    "## 2.3 Steady State Function <a class=\"anchor\" id=\"s2-3\"></a>\n",
    "\n",
    "This function takes model parameters as input and outputs the steady states for each variable following the variable order specified in `Var_t` and `Var_tp1`. To find the steady state, we remove the time subscripts from each variable and solve the system of equations. \n",
    "\n",
    "Note that the first entry of the output array will always be 0, since the steady state of `q_t` is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e0d871c4",
   "metadata": {
    "id": "e0d871c4",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def ss_onecap_2d(*args):\n",
    "\n",
    "    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k, sigma_z1, sigma_z2 = args\n",
    "\n",
    "    def f(imk):\n",
    "        g_k = ϕ_1 * np.log(1.+ ϕ_2 * np.exp(imk)) - α_k - 0.5 * (sigma_k**2)\n",
    "        sdf_ex = np.log(β) - ρ*g_k\n",
    "        cmk = np.log(a - np.exp(imk))\n",
    "        mk = 1./(anp.exp(sdf_ex  + g_k)* 1/(1 + ϕ_2*anp.exp(imk)))\n",
    "        return mk - anp.exp(cmk) - anp.exp(sdf_ex + g_k) * mk\n",
    "\n",
    "    imk = optimize.bisect(f,-10,-2.5)\n",
    "    cmk = np.log(a - np.exp(imk))\n",
    "    g_k = ϕ_1 * np.log(1. + ϕ_2 * np.exp(imk)) - α_k - 0.5 * (sigma_k**2)\n",
    "    sdf_ex = np.log(β) - ρ*g_k\n",
    "    mk = 1./(anp.exp(sdf_ex + g_k)* 1/(1 + ϕ_2*anp.exp(imk)))\n",
    "\n",
    "    return np.array([0., cmk, imk, mk, 0., g_k])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aadc21b2",
   "metadata": {
    "id": "aadc21b2"
   },
   "source": [
    "## 2.4 Specifying Number of Variables <a class=\"anchor\" id=\"s2-4\"></a>\n",
    "\n",
    "We need to specify the number of jump variables, state and shock variables as an array `(n_J, n_X, n_W)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6f4f50b9",
   "metadata": {
    "id": "6f4f50b9",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "var_shape = (3, 2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "715b17d5",
   "metadata": {
    "id": "715b17d5"
   },
   "source": [
    "## 2.5 Model Parameters <a class=\"anchor\" id=\"s2-5\"></a>\n",
    "\n",
    "Next, we need to specify all the model parameters using a Python tuple. \n",
    "<!-- Here, we use the specification described in <insert reference here> -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "43a4ac46",
   "metadata": {
    "id": "43a4ac46",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "rho_option = 1.5\n",
    "σ_k = 0.477*2\n",
    "sigmaz1 = 0.011*2\n",
    "sigmaz2 = 0.025*2\n",
    "\n",
    "delta = 0.01\n",
    "rho = rho_option\n",
    "alpha = 0.0922\n",
    "\n",
    "a = alpha\n",
    "ϕ_1 = 1/8\n",
    "ϕ_2 = 8\n",
    "beta_k = 0.01\n",
    "beta_z1 = np.exp(-0.056)\n",
    "beta_z2 = np.exp(-0.145)\n",
    "beta_z3 = np.exp(-0.1)\n",
    "\n",
    "sigma_k = σ_k*0.01\n",
    "sigma_z1 = sigmaz1\n",
    "sigma_z2 = sigmaz2\n",
    "\n",
    "α_k = 0.04\n",
    "β = np.exp(-delta)\n",
    "γ = 12.0\n",
    "ρ = rho \n",
    "\n",
    "args = (γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z1, sigma_k, sigma_z1, sigma_z2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "238e1c2a",
   "metadata": {
    "id": "238e1c2a"
   },
   "source": [
    "## 2.6 Log Growth Consumption <a class=\"anchor\" id=\"s2-6\"></a>\n",
    "\n",
    "Finally, we need to specify how to compute the log-growth of consumption $\\log{C_{t+1}/C_t}$ using the defined variables `Var_t` and `Var_tp1` in the equilibrium conditions. \n",
    "\n",
    "- In habit formation models, $C_t$ is replaced with $U_t$, in which case we would specify the log growth process for $U_t$ correspondingly.\n",
    "- For some models with preference shocks, $C_t = \\lambda_t \\tilde{C}_t$, we would need to specify the log growth process for $\\lambda_t \\tilde{C}_t$. \n",
    "\n",
    "In this example, we simply use the decomposition:\n",
    "\n",
    "$$\n",
    "\\log{C_{t+1}/C_t} = \\log{C_{t+1}/K_{t+1}} + \\log{K_{t+1}/K_t} - \\log{C_t/K_t}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2e1398d3",
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def gc_onecap_2d(Var_t, Var_tp1, W_tp1, q, *args):\n",
    "    \n",
    "    # Unpack model parameters\n",
    "    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k, sigma_z1, sigma_z2 = args\n",
    "\n",
    "    # Unpack model variables\n",
    "    q_t, cmk_t, imk_t, mk_t, Z_t, gk_t = Var_t.ravel()\n",
    "    q_tp1, cmk_tp1, imk_tp1, mk_tp1, Z_tp1, gk_tp1 = Var_tp1.ravel()\n",
    "\n",
    "    # Compute log consumption growth\n",
    "    gc_tp1 = cmk_tp1 + gk_tp1 - cmk_t\n",
    "    \n",
    "    return gc_tp1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "020b982b",
   "metadata": {
    "id": "020b982b"
   },
   "source": [
    "## 2.7 Computing ModelSol <a class=\"anchor\" id=\"s2-7\"></a>\n",
    "\n",
    "Now we are able to use steps 2.1 to 2.5 as inputs to the `uncertain_expansion` function. The solution is stored in a class named `ModelSolution`(Please refer to uncertainexpansion.ipynb for details.) under the Linear Quadratic Framework. \n",
    "\n",
    "See [src/uncertain_expansion.py](https://github.com/lphansen/RiskUncertaintyValue/blob/main/src/uncertain_expansion.py) for the source code of the expansion suite.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "217c9448",
   "metadata": {
    "id": "217c9448",
    "outputId": "aab57b59-6dce-4452-dfde-a6ee003cdbd7",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1: error = 0.0211794947\n",
      "Iteration 2: error = 3.46944695e-17\n"
     ]
    }
   ],
   "source": [
    "# Solve the One-Capital Adjustment Model\n",
    "eq = eq_onecap_2d\n",
    "ss = ss_onecap_2d\n",
    "gc = gc_onecap_2d\n",
    "var_shape = (3, 2, 2)\n",
    "args = (γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z1, sigma_k, sigma_z1, sigma_z2)\n",
    "approach = '1' # See Exploring Recursive Utility Appendix for difference about approach '1' and '2'\n",
    "ModelSol = uncertain_expansion(eq, ss, var_shape, args, gc, approach = '1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "961bdc55",
   "metadata": {
    "id": "961bdc55"
   },
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "# 3 Outputs\n",
    "\n",
    "## 3.1 List of Outputs <a class=\"anchor\" id=\"s3-1\"></a>\n",
    "\n",
    "We now examine the contents of `ModelSol`, which contains the attributes listed below. Each approximation is stored in a class `LinQuadVar`, which contains the coefficients for $X_{1,t}, X_{2,t}, X_{1,t}'X_{1,t}, W_{t+1}, W_{t+1}'W_{t+1}, X_{1,t}'W_{t+1}$ and the constant.\n",
    "\n",
    "\n",
    "```{list-table}\n",
    ":header-rows: 1\n",
    "\n",
    "* - Input\n",
    "  - Type\n",
    "  - Description\n",
    "* - `JXn_t`\n",
    "  - *LinQuadVar*\n",
    "  - Approximation of jump and state variables at time $t$. Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation. The variables are indexed in the order specified in Section 2. \n",
    "* - `Jn_t`\n",
    "  - *LinQuadVar*\n",
    "  - Same as `JXn_t` but limited to jump variables.\n",
    "* - `Xn_tp1`\n",
    "  - *LinQuadVar*\n",
    "  - Same as `JXn_t` but limited to state variables.\n",
    "* - `JXn_t_tilde`\n",
    "  - *LinQuadVar*\n",
    "  - Same as `JXn_t` but using distorted measure. This variation is also available for `Jn_t` and `Xn_tp1`.\n",
    "* - `util_sol`\n",
    "  - *dict*\n",
    "  - Dictionary containing solutions of the continuation values, including $\\mu^0, \\Upsilon_0^2, \\Upsilon_1^2,$ and $\\Upsilon_2^2$ etc.\n",
    "* - `vmrn_tp1`\n",
    "  - *LinQuadVar*\n",
    "  - Approximation of continuation values $\\widehat{V}^1_{t+1}-\\widehat{R}^1_t$ . Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation. \n",
    "* - `gcn_tp1`\n",
    "  - *LinQuadVar*\n",
    "  - Approximation of consumption growth $\\widehat{C}_{t+1}-\\widehat{C}_t$ . Replace `n` with `0,1,2` to get the zeroth, first and second-order contribution. Omit `n` to get the full approximation. \n",
    "* - `ss`\n",
    "  - *dict*\n",
    "  - Steady states for state and jump variables\n",
    "* - `log_N_tilde`\n",
    "  - *LinQuadVar*\n",
    "  - Approximation for the log change of measure\n",
    "```\n",
    "\n",
    "For example, we can obtain the coefficients for the first-order contribution of $\\log{C_t/K_t}$ in the following way, noting that `cmk_t` was listed as the first jump variable when we specified the equilibrum conditions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3fa57ffa",
   "metadata": {
    "id": "3fa57ffa",
    "outputId": "1cdb816b-4ed7-4c19-e009-60fdcf049004",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': array([[-0.02117949]]), 'x': array([[ 0.04344641, -0.        ]])}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Log consumption over capital approximation results, shown in the LinQuadVar coefficients form\n",
    "ModelSol['JX1_t'][0].coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4f132c",
   "metadata": {},
   "source": [
    "We can also display the full second-order approximation of $\\log{C_t/K_t}$ using the `disp` function, which renders a `LinQuadVar` object in LATEX form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a21c4b7d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\log\\frac{C_t}{K_t}=-3.683+\\begin{bmatrix}0.04379&-1.477e-19\\end{bmatrix}X_t^1+\\begin{bmatrix}0.02172&-7.372e-20\\end{bmatrix}X_t^2+X^{1T}_{t}\\begin{bmatrix}-0.0002216&5.549e-21\\\\-4.919e-21&-1.892e-38\\end{bmatrix}X^1_{t}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Log consumption over capital approximation results, shown in the Latex analytical form\n",
    "disp(ModelSol['JX_t'][0],'\\\\log\\\\frac{C_t}{K_t}') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c2d4266",
   "metadata": {},
   "source": [
    "Another example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d5b18ec2",
   "metadata": {
    "id": "d5b18ec2",
    "outputId": "1b2ebb60-d294-4281-e573-a76dcbc55b4d"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\log\\frac{K_{t+1}^2}{K_t^2}=\\begin{bmatrix}0.9455&-0\\end{bmatrix}X_t^2$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Log capital growth process second order approximation results\n",
    "disp(ModelSol['X2_tp1'][0],'\\\\log\\\\frac{K_{t+1}^2}{K_t^2}') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0190ee8e",
   "metadata": {
    "id": "0190ee8e"
   },
   "source": [
    "\n",
    "## 3.2 Approximate Related Variables <a class=\"anchor\" id=\"s3-2\"></a>\n",
    "\n",
    "We can use the method `ModelSol.approximate` to approximate variables which can be expressed as functions of the state and jump variables. For example, suppose we want to approximate log investment growth $\\log{I_{t+1}/I_t}$. We can express this as:\n",
    "\n",
    "$$\n",
    "\\log{I_{t+1}/I_t} = \\log{I_{t+1}/K_{t+1}} + \\log{K_{t+1}/K_t} - \\log{I_t/K_t}\n",
    "$$\n",
    "\n",
    "We use this expression to write the function `gi_tp1_approx`, which is used as an input in `ModelSol.approximate` below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "145a47f6",
   "metadata": {
    "id": "145a47f6",
    "outputId": "ae3e6f0d-39c8-47a6-b2b3-934eaa6cffaa",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': array([[0.00035532]]),\n",
       " 'x': array([[ 0.01018563, -0.        ]]),\n",
       " 'w': array([[ 0.00917057, -0.00083961]])}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def gi_tp1_approx(Var_t, Var_tp1, W_tp1, q, *args):\n",
    "\n",
    "    # Parameters for the model\n",
    "    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z1, sigma_k, sigma_z1, sigma_z2 = args\n",
    "\n",
    "    # Variables: (q_t, q_tp1 is excluded when using the method in `ModelSolution`)\n",
    "    cmk_t, imk_t, mkmc_t, Z_t, gk_t = Var_t.ravel()\n",
    "    cmk_tp1, imk_tp1, mkmc_tp1, Z_tp1, gk_tp1 = Var_tp1.ravel()\n",
    "    \n",
    "    gi_tp1 = imk_tp1 + gk_tp1 -imk_t\n",
    "    \n",
    "    return gi_tp1\n",
    "\n",
    "n_J, n_X, n_W = var_shape\n",
    "gi_tp1_list = ModelSol.approximate(gi_tp1_approx, args)\n",
    "gi_tp1 = {'gi_tp1':gi_tp1_list[0],\\\n",
    "        'gi0_tp1':gi_tp1_list[1],\\\n",
    "        'gi1_tp1':gi_tp1_list[2],\\\n",
    "        'gi2_tp1':gi_tp1_list[3]}\n",
    "gi_tp1['gi1_tp1'].coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd28ede",
   "metadata": {
    "id": "bbd28ede"
   },
   "source": [
    "<br>\n",
    "\n",
    "## 3.3 Simulate Variables <a class=\"anchor\" id=\"s3-3\"></a>\n",
    "Given a series of shock processes, we can simulate the path of our state and jump variables using the `ModelSolution.simulate` method. Here, we simulate 400 periods of i.i.d standard multivariate normal shocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "58417c59",
   "metadata": {
    "id": "58417c59"
   },
   "outputs": [],
   "source": [
    "Ws = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 400)\n",
    "JX_sim = ModelSol.simulate(Ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "738d5ef7",
   "metadata": {
    "id": "738d5ef7"
   },
   "source": [
    "<br>\n",
    "\n",
    "## 3.4 Shock Elasticities and IRF <a class=\"anchor\" id=\"s3-4\"></a>\n",
    "\n",
    "`ModelSolution` contains a method `elasticities` used to compute the shock elasticities of jump and state variables. We do this by defining a `log_SDF_ex` function, which expresses the log stochastic discount factor exclusive of the change of measure ${N}^*_{t+1}$ and variable $Q_{t+1}^*$. Shock elasticities are great tools for nonlinear impulsive response analyses, and in linear models, the shock elasticities correspond to the traditional impulse response function (IRF). \n",
    "\n",
    "`ModelSolution` also has a method `IRF` to compute the traditional IRF of all jump and state variables. Below we calculate the IRF of all variables with respect to the first shock for 400 periods.\n",
    "\n",
    "The figure plots elasticities and IRF of investment rate with regard to shock 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "c2fb03dd",
   "metadata": {
    "id": "c2fb03dd",
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def log_SDF_ex(Var_t, Var_tp1, W_tp1, q, *args):\n",
    "\n",
    "    # Log stochastic discount factor exclusive of the change of measure N and Q.\n",
    "\n",
    "    # Parameters for the model \n",
    "    γ, β, ρ, a, ϕ_1, ϕ_2, α_k, beta_k, beta_z, sigma_k1, sigma_k2, sigma_z = args\n",
    "\n",
    "    # Variables: (q_t, q_tp1 is excluded when using the method in `ModelSolution`) \n",
    "    cmk_t, imk_t, mk_t, Z_t, gk_t = Var_t.ravel()\n",
    "    cmk_tp1, imk_tp1, mk_tp1, Z_tp1, gk_tp1 = Var_tp1.ravel()\n",
    "\n",
    "    sdf_ex = anp.log(β) - ρ*(cmk_tp1+gk_tp1-cmk_t)\n",
    "    \n",
    "    return sdf_ex\n",
    "\n",
    "expo_elas, price_elas = ModelSol.elasticities(log_SDF_ex, T=40, shock=0)\n",
    "states, controls = ModelSol.IRF(T = 400, shock = 0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABBYAAAFPCAYAAAD5m00yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAC/9klEQVR4nOzdeXhMZ/vA8e9kExIhiYglhIQQYitiCbHvSy0tqrai1pZaquiq1dpq6WuntdXy2pUS1BavpUpr34pshKQkskwi28z5/ZFfpiaZxCDJZLk/1+Vq55znnHNPRo6Ze57nvlWKoigIIYQQQgghhBBCvAIzUwcghBBCCCGEEEKI/EsSC0IIIYQQQgghhHhlklgQQgghhBBCCCHEK5PEghBCCCGEEEIIIV6ZJBaEEEIIIYQQQgjxyiSxIIQQQgghhBBCiFdmYeoACrupU6eye/fuLMe0adOGZcuW5VJE+c/ixYtZsmRJlmOqV6/OL7/8AsDAgQMJDQ3l2LFjORJPUlIST58+xdnZGYBdu3Yxbdo0NmzYQKNGjYw6x7lz5xg0aBCzZs2iV69euu3379+nQoUK2Rbr48eP+emnnzhx4gQPHz7E1taWKlWq0KdPHzp16oS5uXm2XetVpP0cjHH79m2qVatGz549mT17dg5HJsS/DN3HzczMKFq0KO7u7vTv35+ePXu+8Dw5fW9KT+6dr07unUKIV5V2b3v+PpXZ54EiRYrg6OhIkyZNmDhxIqVKldLtM+YzxPnz57Gzs8tyjFqtZu3atfz222+EhIRQpEgRKleuTI8ePejZsydFihR5hWeZfR48eECbNm2MGnv06FEGDRpE+fLl+fnnn3M4MpGeJBbyiGnTpmFvb29wX9myZXM5mvxp1KhRuLm5GdxXsmTJXIkhNDSUoUOHMnLkSN0/Fg0bNmTu3Lm4u7sbfR53d3fmzp3LG2+8odv2xRdfEBgYmG03yvPnzzNmzBhSUlJ48803qVGjBjExMZw4cYJJkybxyy+/sHDhQmxtbbPleq8i7efwvClTpuDm5saoUaMyjJ87dy4VK1bMrfCE0PP8fVxRFNRqNXv37mXq1Kk8ffqUoUOHZnn8qFGjePbsWW6EmuG6cu80ntw7hRA5Jf3nAbVazdmzZ9m5cyfXrl1jx44dWFlZZXnM84oWLZrl9e7cucPw4cOJiIigS5cu9O/fn4SEBH7//Xe++uorduzYwbJlyyhduvTrP7lX5ODgkOF+NmvWLCD1uacfO3369Bc+b5EzJLGQR7Rt2xYXFxdTh5GvNW3a1OhvtXLKgwcPCAoK0ttWoUKFl/6mrFSpUrz55pt6206dOkX58uVfN0Qg9du7UaNGUbp0aX788Ue98w4fPpytW7fy5Zdf8vnnn7Nw4cJsuearMPRzmDJlisHtgMFtQuQWQ/fxt956i86dO7N06VIGDBiQ4Q3h83x8fHI6RIPk3mk8uXcKIXKSoX9H3n33Xb766iu2bNnCkSNH6Ny58wuPMUZMTAwjRoxAo9GwY8cOqlevrts3ZMgQjh07xkcffcS4cePYvHkzZmamWUFfrFixDPeoH374ATB872rbtm2uxCUykhoLQhRCixYtIj4+nsWLFxt8w923b1969+7NgQMH+OOPP0wQoRAFg7W1Na1bt0atVnPnzh1ThyNek9w7hRCmkLac7vLly9l2zp9++omHDx8ya9YsvaRCmtatWzNq1CguXryoWxInRFYksZCPnDx5kmrVqjF+/Hi97Z9//jnVqlXj5MmTQOqN4NNPP2X79u20adOGunXr0q9fP37//fcM57xw4QJDhgyhXr161KtXj0GDBnH+/Hm9MdHR0UydOpWWLVvi5eVF27ZtmT9/PomJiboxAwcOpHXr1hnOn377wIEDGTZsGAsXLqRevXo0adKE27dvA3D37l3Gjh1LgwYNqFOnDv369eN///vfq//AXsLBgwcZMGAA9evXx8vLi9atWzN37lySkpJ0Y5KSkvj2229p06YNXl5etGjRghkzZhAdHQ2krplLW9M6bdo0qlWrptterVo1zp07p3euxYsX0759e2rXrk2HDh1YtWoVGo0GSF0fW61aNXbt2gVAtWrVCA0N5Y8//tBt79u3L82aNUOr1eo9l3v37lGtWjU2bdpk8LkmJCRw6NAhGjduTJUqVTL9maRN2967dy8AnTp1olu3bhnGpcW6Z88e3bZdu3bRo0cPatWqRePGjZk6dSr//POPbv+DBw+oVq0a69at45133sHLy4shQ4ZkGouxqlWrxtSpU/Uer169mlWrVtGyZUvq1KnDwIEDCQ4OJjAwkGHDhlG3bl1at27Nhg0bMpzvRc9DCGOoVCoA3e9369at+eyzz5g+fTq1atXC19eXyMhIg/fRe/fuMX78eBo1akT9+vUZOHAgFy5c0Bsj9065d74uuXcKkbelTe1XFCXbzrlnzx5cXV1p3rx5pmMGDRqEpaWl7n42fPhwGjVqREpKit64tHvT83V7jh8/Tr9+/ahTpw4NGzbkww8/JDAwUO+4atWqsWjRIkaNGoWXlxedO3fOcO6X1bp1awYOHKj3+Ouvv2b79u106NCB2rVr07t3b65cucLjx48ZP3489erVo3nz5ixcuDDDvw3GPA+RShILeURMTAyRkZEG/6S9YfL19aVnz54cPHhQ96bx9OnTbNu2jX79+uHr66s735kzZ/j666/p0KED48ePJzIykuHDh+t9g3L06FEGDhzIo0ePGD16NKNHj+bRo0cMGTKEo0eP6sZ99NFHHD9+nLfffpsvv/wSb29vVq1axcyZM1/puf7111/s37+fjz/+mJ49e1KlShVu375N3759uXv3LiNHjmTChAmkpKQwYsQIDhw4YNR5Y2NjM/0ZPv8mN73t27czfvx4ihcvzuTJk5kyZQrly5fnp59+YtWqVbpxaTelLl268OWXX9KhQwe2bdvGhAkTgNT1wGlrV/v27ZthPdjzxo4dy5IlS6hduzbTpk3jjTfeYP78+ZkeM3fuXOzt7XFzc2Pu3Lk0bNiQrl278vjx4wyJoP3792NhYUGnTp0MnuvatWskJydTt27dTOOD1DW6pUuX1n2I6datG3///Tf37t3TG3fgwAGsra11U8+WLFnCtGnTqFixItOmTaNv37789ttv9OvXj8jISL1jf/jhB5ydnZk+fbrBN97Z4eeff2bnzp0MHTqUIUOG8Ndff/Hhhx8yePBgypcvz9SpU7G3t+fbb7/V+/14mechRGa0Wi1//PEHVlZWerUC9u/fz61bt/j000/p06cPDg4OGY4NCgqiT58+/P777wwYMICJEycSFRXFe++9x5UrVwDk3in3Trl3ClEIpL3vr1GjRoZ9mX2GSP8B+XlhYWGEhYW98H5ma2uLl5eX3v0sKiqKM2fO6I1L+/cm7X60a9cuRo8eTdGiRfn4448ZMmQIFy9epE+fPhk+lK9fv56EhAQ+++wz+vTpg4VF9q/UP3r0KD/88ANvvfUWH3zwAQEBAXz44Ye89957mJmZMXXqVDw8PFixYoXe7IyXeR4CUIRJffLJJ4qHh0eWf27cuKEbHxUVpfj4+Cjt2rVTIiMjlZYtWypt27ZV4uLidGNatWqleHh4KL/99ptuW0REhNKgQQOlT58+iqIoSnJysuLr66u0aNFCiY2N1Y2Ljo5WmjdvrjRv3lxJSkpSnjx5onh4eCg//vijXtxTp05VBg8erHs8YMAApVWrVhmeX/rtAwYMUDw8PJTff/89w7j0zyM5OVnp37+/0rRpUyUxMTHTn+F//vOfF/4Mn/9ZpI+pY8eOSt++fRWtVqt3bV9fX6Vr1666bbVr11ZmzJihd+2FCxcqvXr1UtRqtaIoivL7778rHh4eys6dO3Vjdu7cqfecT5w4oXh4eCjLly/XO9ekSZOUmjVrKlFRUQbP06pVK2XAgAG6x0+ePFFq1KihfPnll3rn6dChgzJs2LBMf14HDhxQPDw8lE2bNmU6Jk3Pnj2VevXqKYqiKCEhIYqHh4eyZMkS3f6UlBSlcePGykcffaQbU716deX777/XO8/t27eVmjVrKt9++62iKIpy//59xcPDQ2nXrp2SlJT0wjie5+HhofdzSL/vk08+0Xtcp04d5fHjx7pt48aNUzw8PJR58+bptgUFBSkeHh7KggULXup5CKEo/97Hr1+/rkRERCgRERHKP//8o1y8eFEZP3684uHhoXz33Xe68a1atVKqV6+uBAcH650n/b1p/PjxSu3atZWgoCDdtsjISKV+/frKuHHjdMfIvVPuncaQe6cQeUPave35+5Shf0ciIiKU4OBgZePGjUrdunWVTp066f3ev+gzxP379zON4fLly4qHh4cyd+7cF8b74YcfKh4eHsrTp08VtVqt1KlTR5k2bZremB49eihvv/22oiiKEhsbq7zxxhvKhAkT9Mb8888/SsOGDZUxY8botnl4eCj169dXoqOjXxjH81q1amXwc0favufvda1atVKqVaum3Lp1S7dtzpw5ioeHh+4erCiKEhcXp9SsWVOZOHHiSz8PkUqKN+YR8+bN02sh87znKzWXKFGCGTNmMGbMGN5++23CwsLYtGkTxYoV0zvGzc1Nr3iJg4MDb775Jhs3biQiIoLQ0FDCwsKYPHmyXuVqOzs7BgwYwPz587l27Ro1a9akWLFibN68GRcXF5o3b06xYsV01VhfhbW1NQ0bNtQ9fvr0KX/88QcDBw4kISGBhIQE3b527doxa9Ysrl69Sv369bM87yeffGJwjRiQ6XZIna767Nkz3XRlgIiICOzs7IiPj9dtK1OmDAcOHNAtB7Gzs+Ojjz7io48+etFT1nPixAnMzMwYMGBAhvhHjx5tdCVxR0dHGjduzG+//cbnn3+Oubk5N27cIDAw0GDV7zTK/0+jMyYjbGlpqRtfoUIF6tWrx8GDBxk7diwAZ8+eJTIyUpeh/u2339BqtbRu3Vrvm6lSpUrh6enJiRMnmD59um5748aNsbS0NOr5vqp69erp/W5VqlQJSP27lSat6FHaVN2XfR5CAAZbSlpZWTFw4EAmTZqkt71ixYpZVuHXarX4+/vTokULXF1dddvt7e3ZvHkz9vb2cu98Ln65d2Y/uXcKkfsM/TtStGhRWrduzeeff27w9z6zzxBOTk4vvJ4x97O0MYqiYGNjQ5s2bTh69CjJyclYWloSGBjIjRs3+Oyzz4DU2dRqtZq2bdvq3QfMzc1p3Lgx/v7+pKSk6M5bp06dF7bEfF0VK1bULbMDqFy5MqB/PytWrBiOjo48fvz4lZ6HkK4QecYbb7xhdEXXNm3a0L59ew4fPsw777yj11YrjaH1n66uriiKQmhoKA8ePAD+/cV6XlrbsYcPH1KvXj2+/vprPv/8c8aNG4eVlRXe3t60b9+eHj16vFJv25IlS+pVlr1//z6QOu0ys3Zgjx49euF5a9as+UqVzS0tLTl//jy//vorAQEBhISEEBERAaBXnOurr77io48+Ytq0aXz++efUrVuXdu3a0bt3b4oXL2709UJDQ3F0dMzwJtjJycmofwSe161bN06dOsX58+dp3Lgx+/fvp0iRIllWxE1rGfTkyZMXnv+ff/7RazHUtWtXvvnmG+7du4e7uzsHDhygZMmSuvV5ISEhAPTr18/g+dL/g2ho+nd2c3R01Huc9g/A89dO6zmf9kHgZZ+HEKD/5s7MzAw7Ozvc3d0N3ifT/71MLyoqivj4eL2kQhoPDw8A3XIIuXfKvTMnyL1TiNyX9u9IcnIy//vf/9i0aROdOnXiq6++yvQ998t8hkjzsvczKysrXUvLrl278uuvv3L27Fl8fX05cOAA5ubmum4VafeBtOVuhkRGRupiMMX9LO3elf7a5ubmGe5nxj4PIYmFfCk+Pp4bN24AqWuu4uPjM8xYMPSPd1qthud/aQxJ25d2jm7dutG8eXOOHDmCv78/Z86c4dSpU2zevJnt27dn2T4t7ZrPS/tlTj/m3XffzfRNXVaFsl7X/PnzWbVqFTVq1KBu3bq8+eab1KtXj2+++UbvTXmTJk04fvy47s/p06eZNWsW69atY9euXUbfGDUajd43fK+jbdu2WFtb4+fnR+PGjfHz86Nly5ZZfnOXNgslfQG49B4+fMjDhw91PeUBOnfuzKxZs/Dz82PkyJEcOXKEjh076v6upK3nW758OdbW1i+MP/3fhZyQWSY5q9fgZZ+HEPByb+5e9Hc/7b6YVXsvuXe+Orl3vpjcO4XIfc//O5I2Y23mzJlERUWxbNmybLsHlilTBhcXlxfez549e8b169epV6+ebluzZs2wt7fHz88PX19f/Pz8aNKkie7De9p94Jtvvsn038QSJUro/j+v38+MfR5CEgv50oIFCwgNDWXKlCnMmzePBQsW6KYfpUnLsj0vODgYc3NzXFxcSE5OBiAgICDDuLRiJGXKlCEuLo6bN29StWpV3nrrLd566y2SkpKYN28eGzZs4NSpU7Ru3RozMzODRb6MyYSmfbNlbm5O06ZN9fbdvXuXBw8e6KrhZrfQ0FBWrVrFm2++maH41/OxJyUlcfPmTcqUKUOXLl3o0qULWq2WtWvXMnfuXPbv369XgTYr5cqV48yZM8TFxWFjY6Pbfv36ddasWcPo0aONjt/W1pZWrVpx/Phx3nnnHUJDQ5k2bVqWxxQtWpR27dqxb98+bt68iaenp8Fxa9euBaB79+66bQ4ODjRt2pSjR49Sq1YtoqOj9fanvZZly5bNcF5/f3+jpyqbWkF5HiL/sre3x9ramuDg4Az7fvrpJ548ecJ7770HyL1T7p15R0F5HkLkFQMHDuTs2bMcPXqU9evXZ0sXmDTdu3dn2bJlHD16lDZt2hgcs3XrVuLj4/XuV5aWlnTs2JGDBw/y999/c+fOHYYPH67bn3YfSLvvPe/cuXNotdosv5TMKwrK88hN0hUin/nzzz/ZtGkTffr0YdiwYfTu3ZuNGzdmyDhevXqVS5cu6R4/efKEvXv30rhxY0qUKEHNmjVxcnJiy5YtqNVq3Ti1Ws3mzZtxcnLCy8uLO3fu8O6777Jjxw7dGCsrK11V2rQsY6lSpYiIiCA8PFw37tq1awbfFKdXunRpvLy82L17t97xycnJTJ8+nXHjxr1265nMpLU7S/+tnr+/P0FBQbrrRkVF0bdvX1auXKkbY2ZmRq1atXT/D//+PLKqxNuiRQu0Wi3bt2/X275lyxb8/PwyrbVhZmZm8LzdunUjPDyclStXUrx4cVq0aJHlcwZ0tTXGjx+vWxbzvD179rBx40Y6duxIkyZNMlzvxo0bbNq0iXLlyumt327VqhUAK1eu1JsVc/PmTUaPHs369etfGFteUFCeh8i/LCws8PHxwd/fX+/b/+joaH766SdCQkLk3vn/5N6ZdxSU5yFEXvL1119TokQJFi1apFs+nB1GjBiBi4sLn332mW4m9PNOnz7NwoULqVu3bobaD926dePp06csWLAAa2trvVoFTZs2pUiRIvz444+6LzIBwsPDGTNmDN9//322zbzISQXleeQmmbGQRxw5ckS3dsmQN998k8TERD799FMcHByYPHkykPom58iRI3z66af88ssvuqmHVlZWvP/++wwePBhra2s2b96MVqtlypQpQGq28fPPP+ejjz6id+/evPXWWwDs2LGDf/75h//85z+YmZlRp04dGjRowMKFC3n06BHVqlXj0aNHbNy4ETc3N90bp7T1Vu+//z7vvPMOERER/Pzzz1SqVEnvlzEzn332GYMHD6Z379688847lCxZkv3793P58mUmTZqU5c8mzZkzZwgLC8t0f7t27TIsGalSpQrlypVjxYoVJCYmUqZMGa5cucLu3bspUqQIcXFxQGryo1u3bmzevJlnz55Rr149oqKi2LhxI6VKldK1J0uLc+/evSiKYrAIT+vWrfHx8WH27NncuXOHWrVqcfHiRfbs2cPYsWMpWbKkwfgdHBy4desWmzdvxtvbW/eG3tfXl5IlS3LgwAF69eplVPa0dOnS/Pjjj4wePZpu3brRs2dPPD09iY+P5+TJk5w6dYpmzZoZbCnatm1bihUrhr+/PyNGjNC7qXp4eDBw4EB+/vlnoqKiaNu2re7nZGNjw/jx418YW15QUJ6HyN8mTZrE22+/zdtvv827776Lra0t27ZtIz4+Xlf4UO6dcu/MSwrK8xAiLylVqhSTJ0/m888/58svv2TNmjXZct6iRYvy008/MWLECPr06UO3bt2oW7cuKSkpnDt3jsOHD1OjRg3+85//ZFiu8MYbb1C+fHmOHz9Oly5d9GaROTg4MHHiRGbNmkXfvn3p3r07KSkpbN68mcTERD755JNsiT+nFZTnkZsksZBHvKjLwptvvsnixYsJDAxk3rx5uuqp9vb2fPzxx3z66af88MMPur/kdevWpUuXLixbtozY2FgaNGjApEmT9Cp8d+jQgTVr1rBs2TKWLl2KhYUFderU4dtvv6VBgwZA6tqjpUuXsmTJEo4fP87WrVspUaIE7du3Z/z48bo3Yq1ateKLL75gw4YNfPvtt1SuXJmvvvqK8+fPc+LEiRc+/3r16rFlyxYWL17M2rVrSUlJoXLlysyePdvgG0xDVqxYkeX+o0ePZnhzbGVlxapVq5g9ezYbNmxAURQqVqzI9OnTSUlJ4dtvv+XatWt4eXnxzTffUKFCBfbv38/+/fspWrQoTZo0YcKECbo1wu7u7gwcOJBdu3Zx9epVgwXRzMzMWLZsGcuWLWPfvn3s3buXihUr8sUXX/DOO+9kGv+HH37Il19+yXfffcfYsWN1b44tLS3p0KEDW7dupWvXrkb9rCC1Cu++ffvYuHEjR44cYffu3VhbW1O1alXmzZtH165dDa7vLlasGK1bt+bXX381eL1PP/0UNzc3/vvf/zJnzhyKFy9OgwYNGD9+PO7u7kbHZ2oF5XmI/Mvd3Z2tW7eyYMECfvzxR8zMzKhduzZz5syhatWqgNw75d6Z9xSU5yFEXvL222+zZ88eTp8+zZ49e+jRo0e2nLdSpUrs2rWLrVu3sn//fg4fPoy5uTlubm588cUX9O7d22DRSJVKRbdu3VixYoXB+9mQIUNwdnZm7dq1LFy4EGtra2rWrMm8efNe2KkoLykozyO3qJSsqviJfKl169aUL18+0yrhouD58ssvOXr0KP7+/rlSBEcIIQoCuXcKIYQQ2UNqLAiRz0VFRXHw4EF69Oghb4yFEMJIcu8UQgghso8shRAin7p+/To//vgjV65cITExkXfffdfUIQkhRJ4n904hhBAi+8mMBSHyqeLFi3P27FlSUlL4/vvvKVu2rKlDEkKIPE/unUIIIUT2kxoLQgghhBBCCCGEeGUyY0EIIYQQQgghhBCvTBILQgghhBBCCCGEeGVSvDEPefo0Dq1Wf2WKo6MtERFqE0UkXpa8XvmPMa+ZmZkKe3ubXIpImEr6e7D8Puc/8prlL8a+XnIPLhwMvQ9OT37H8wd5nfKH7L4HS2IhD9FqFYM31BfdZEXeIq9X/iOvmQDD92D5u5H/yGuWv8jrJdJk9j7Y0DiR98nrlD9k5+skSyGEEEIIIYQQQgjxyiSxIIQQQgghhBBCiFcmiQUhhBBCCCGEEEK8MkksCCGEEEIIkUf9+uuvdOnShdq1a9OpUyf27NmT5fi4uDhmzJiBj48P9erV4/333ycoKEhvTEpKCosWLaJFixbUqVOH/v37c+XKlQznWr9+Pe3ataN27dr07NkTf3//145PCFEwSWJBCCGEEEKIPMjPz4/Jkyfj4+PD0qVL8fb25pNPPuHgwYOZHjNhwgQOHjzI5MmTmTNnDuHh4QwaNIjY2FjdmG+//ZZ169bx/vvvs3DhQszNzRkyZAj379/Xjfnxxx+ZM2cOPXv2ZPHixVSoUIExY8Zw8eLF14pPCFEwqRRFkZKdeUREhDpDZU4np+I8fhybyREir5HXK/84ez2MXf73iIxJxMGuCL1auNOkZhmDY83MVDg62uZyhCK3pb8Hy+9z/iOvWd7n93swlcra4elqr3u9bgY/JehRDJ0auxo8pjDfg9u1a4eXlxcLFy7Ubfvoo4+4ffs2fn5+GcZfuHCBd999l9WrV+Pr6wtAZGQkbdq0YfTo0YwYMYIHDx7Qvn17Pv/8c9555x0AkpKS6NChA76+vsyYMYP4+Hh8fX3p168fkydPBkBRFPr160fx4sX58ccfXym+rBh6H5ye/I5DSHgsszZe0PtZqVQqHOysiYxJILOPdtkxxthzODkU43FkfI5fJ7eeT0G9jpmZGdPefYMKzsUNjk1j7D1Y2k0KIQqds9fDWO93i6QULQARMYms97sFkGlyQQghxOurVNaO5XuuMbqHF05OxbkZ/FT3WOi7f/8+ISEhTJw4UW97hw4d8PPz4/79+1SoUEFv3+nTp7GxscHHx0e3zcHBgYYNG3Ly5ElGjBjB77//jkajoUOHDroxVlZWtGzZkhMnTgBw+fJlYmNjad++vW6MSqWiXbt2LFy4kKSkJMLDw186PgFR6kQWbb9EWETqB++X/TCYlGJojEJY5LMXXDk7xhh3jtB/4nLlOrn1fArudTSs3HuDme83esExxpHEghCiQEubmRARk4jj/89M2OV/T5dUSJOUomWX/z1JLAghRDZ5fnbC8+pVLcWSXVfxv/yQG4GRjO7hlWGMgICAAAAqV66st93VNXVmR2BgYIYP7gEBAbi6umJubq63vWLFiroZBAEBAZQoUQIHB4cM53348CEJCQm6a7u5uWUYk5KSwv3793nw4MFLx1eYRKkTWbLrColJGh5HPdMlDTRa0J+YkV0fOoV4eQ8j4rgfHvvCWQvGkMSCEKLAMjQz4af9NzOdahkRk5ib4QkhRIH2/OwEt3J2HP4jhH1ngrAwNyMhScMfN8Lp1KiiJBUykVYTwdZWfwqyjY0NAGq1OsMxarU6w/i0Y9LGZzUGUos/po1N22ZozKvElxVjl7s4Ob3+B6CcEhmTwIwffyf0nxiSUxResLJDiDzhpwM3WTqlzWufRxILQogCI/3shMRkTYaZCVmt33S0K5LTIQohRIFkaHZCUrIGFycbFm67hFab+iGriKUZbuXsCHgYQxefyhz6PRgvN0dJLhiQ9g23SqUyuN3MLGMN9qxKp6WNz2zM89dLm6L/ojEvG19W8nONhbS6B5JMEPlRSLiav649zHTWgtRYEEIUKoZmJ2TFysJML+lgZWFGrxbuORqjEEIUVGmzE4Z29iQ+MZnjf4Vy72EMkJpMSNEoNK7pTFOvMqzae4MPe9fGt0FF3MoU181qkOSCvuLFU9/kp//mPy4uTm//82xtbXVLFNIfkzazwNbWVncOQ+e1tbWlePHiKIqid1z6a79KfAXJ80sdQp/EmzocIV5LdtRakMSCECLfST8zoaevG9uO3c0wOyEzz9daMKYrhBBCiFSGZiZcufeEs9fDKeNQlP/svAKASgX1PZyoVrEke08H0b5hRY5fDMXKwlwvieDpas/oHl4EPYqRxEI6abULQkJCqFatmm57cHCw3v70x5w9ezbDjIPg4GDdeDc3N6KiooiOjqZEiRJ6Y1xcXLCystK7do0aNfTGWFlZUa5cOd35Xya+giIkPJZv1p9HY9zbDiHyvMdRr1/DQxILQoh8xdDMhB9/vWn08WkzE5rULEOTmmXy7LRKIYTIi9JmJozoXgONRuHwhfvcDHoKQAlbKyqXtSPwUQxdmrji6eqgNxuhuqs9y/dco1ENZ71zerraS1LBAFdXV1xcXDh48CDt2rXTbT98+DCVKlWiXLlyGY5p1qwZK1as4MyZM7rOEJGRkVy4cIGRI0cC0LRpUwAOHTpEnz59gNR2k/7+/jRr1gyAevXqUaxYMQ4dOqRLLCiKwm+//UbDhg2xsrJ6pfjyu7SODiHhL+p6kDsszc1Y+XFLU4eRgby3yh+y+3WSxIIQIk9LPzvhWVLGugmQ+u2YoWWjNtbmWFtZ6HWFkJkJQgjxYulnJyiKQlhkHHY2VizYelk3rra7Ix28K6JoFVbsvU63ppU4fjGUmLhkmZ3wmsaOHcu0adMoUaIELVu25NixY/j5+bFw4UIgNWkQEhJClSpVsLW1pWHDhnh7ezNx4kQmT55MyZIlWbx4McWLF+edd94BoHz58vTs2ZOZM2cSHx+Pq6sra9euJTo6muHDhwNQtGhRhg4dyrJlyzA3N6dOnTrs3LmT69evs2HDBqPjK0ii1Il8ueYPYuOTs/W8FUrbMmOod7aeUwhTkMSCECLPepm6CYpiuG5C/3bVJJEghBCvIG12wsD2HjyJTuD4X6E8iUnAzAzKOhbjUUQ8XZq40ruFOzeDn7Ji73WZnZDNevXqRVJSEmvWrGH79u1UqFCBOXPm0LlzZwBOnDjBtGnT2LBhA40apa6PXrJkCbNnz2bu3LlotVrq16/PokWL9JY9fP3119jZ2bFq1Sri4+OpWbMma9eu1bWKBPjggw8wNzdn27Zt/Pjjj1SpUoVly5ZRv359o+MrKKLUiXz+0zninqW88jkkgSAKOpWSVflYkasMVcOVqUT5i7xeryf97IT4xBSeJWqMOvb5ugkvMzvBmNfM2Gq4In9Lfw+W3+f8R16zV2OobsK1wAh+vx5OeGS8rgijmZmKtm+4UM21JGsP3KJVvfIcvxiqm4WQ/hw3g58S9CiGTo1dM1wTjH+95B5cOOTVrhCvk1QorMkEuRfnD9l9D5YZC0KIPOFluzo8L33dBCGEEMZLm5kwuocXDsWLsOd/gfxxMxyF1LoJVcqX4G5oNF0au+pmIqSfmWCoq4PMThD53askFfJq3QMhcpokFoQQuS79zIReLdxfqquD1E0QQohXY2h2gkarpaKzLQu2XkLz/98Yu5ezo0uTSlhaqFi594aubkJ0XJLUTRCFxubf/jY6qVBYZycIkUYSC0KIXGWwq8O+Gxi7JkvqJgghxKt7fnZCGYdi7Dxxl7PXU2cnFLUy51mShnYNXHinrQc3g59KVwdRaIWEx3Lh9uMXjpMZCkKkksSCECJHpZ+dkJicsauDAqj+/7/pyewEIYR4eYZmJtwMfkrgo2i6Na3Ewm2XSNGk3nUrl7WjbtVS/Hb+Pt2aVuD4xVDqVnUi6FGMzE4QhVJIeCxfrT3/wnEyS0GIf0liQQiRY16qqwPS1SE/+/XXX1m+fDn379+nfPnyjBw5kh49ehh17Jw5c7h58ybr1q3T237hwgXefffdDONbtmzJypUrAXTVyPfs2UNYWBhlypShe/fuDB8+HCsrq9d9WkLkW8/PTPB0tedaQATL9lyjhI0V4U+fYWmuAqD1G+WpX6201E0Q4v9FqROZuV6SCkK8LEksCCFyzC7/e0bXTXjVrg7C9Pz8/Jg8eTKDBg2iefPmHDlyhE8++QRra2s6duyY5bEbN25kzZo1NGnSJMO+27dvU6xYMdauXau33c7OTvf/M2fOZO/evYwZM4ZatWpx9epVli5dSlhYGF9//XX2PEEh8qG02QXLdl/FxcmWv+9HoQD2xVW0b1CBM9fD6Pj/XR1SNIrMTBDi/63zu8mL3rpIUkGIjCSxIITIFumXPHRtWsnozg7S1SF/W7BgAZ06dWL69OkANG/enOjoaH744YdMEwvh4eHMnTuXAwcOULx4cYNjbt26RdWqValbt67B/dHR0Wzbto3JkyczfPhwAF2CYv78+UyaNEmvb7sQBZWhZQ9nr4dx7M8HPEtM4fb9KEqXLMqADh6Yq1Qs/+W61E0QwoCQ8Fiu3IvMcowkFYQwzMzUAQgh8r+0JQ9piYSImETWH7yd6Xgba3Mc7YoAqTMVBneqLgmFfOr+/fuEhITQvn17ve0dOnQgICCA+/fvGzxu4cKF3Lhxg7Vr1+Lp6WlwzM2bN6lWrVqm11ar1fTr14/WrVvrbXdzc9PFJkRhkLbs4WbwU4LDYpm7+S9W77tBwMMYzMzMaFmvPPGJKZibmREUFpvp7AQhCrtF2y5nuV+SCkJkTmYsCCFemjEFGSG1wrhGq0jdhAIsICAAgMqVK+ttd3V1BSAwMJAKFSpkOG748OG4ublhZmbG0qVLM+zXarXcuXMHe3t7evbsyZ07dyhVqhSDBg3ivffeQ6VSUb58eb766qsMxx49ehRLS0tdDEIUdJ6u9nRrWkmvXWS1CiV58FjNmJ618HS1p2H10lI3QYgsXA+MICouKdP9klQQImuSWBBCvJSXKcj4LEnD+91qSN2EAiw2NhYAW1tbve02NjZA6qwCQ6pUqZLleQMDA0lISCAwMJCJEydib2/P0aNHmTt3Lmq1mnHjxhk87rfffmP37t0MGjQo0yUWmXF0tM2wzcnp5c4hTK+gv2Y7j92hasWS1K7ihKIoXLnzhB/3XSPoYQxFLM3RaDX0alkFOxsr3ThI/bmULFmUOyFR+Oahn1FBf71E/rFk19VM91mam0lSQYgXkMSCEOKlvGxBRqmbULApSuq3oyqVyuB2M7NXW3Hn7OzM6tWr8fT0xMkp9YNRkyZNSEhIYPXq1QwdOjRDMuPw4cNMmjSJ+vXrM2nSpJe+ZkSEGq3236anTk7Fefw49pXiF6ZRGF4zJ7sizFp3ng7eFbl09zH3QmNQAfWqlOJOaDTtG1bg8LlgRvfwomwJa72fR9kS1pStVSbP/IyMfb3MzFQGE39CZJfrgREkJht+b2Nb1IL/jPfN5YiEyH8ksSCEyFT6JQ+9Wri/dEFGUbClzQpIPzMhLi5Ob//LsrW1xdc34xu5li1bsn37dgIDA6lVq5Zu+7p165gzZw7e3t4sXbqUIkWKvNJ1hcgrDBVkvBn8lHM3wrArZslO/3sUtTKniKU53XxcOXjuvlHtIoUQGWU1W8G+uHUuRiJE/iXFG4UQBhkqyPjTrzczHS8FGQuntNoKISEhetuDg4P19r+s27dvs3nzZpKTk/W2JyQkAGBv/++Hpe+++45Zs2bRuXNnVq9enWEmgxD50fMFGQGOXLjPgq2XOHn5Ec+SNNR2d+RZkob2DSugQiUFGYV4RSHhsVnOVpAlEEIYR2YsCCEMMrTkQasoWJirUKlUJEtBRkFqkUYXFxcOHjxIu3btdNsPHz5MpUqVKFeu3CudNzg4mBkzZuDs7EybNm102w8cOICLiwvly5cHYNGiRaxfv5733nuPTz75JMOSDCHyq7TkwNJdV7EpasHjqASKWVvQo1llyjgUZdW+m3RrWonjF0OlIKMQr+E/O69kuk9mKwhhPEksCCGAjMseMlvykKJReL+bpxRkFDpjx45l2rRplChRgpYtW3Ls2DH8/PxYuHAhAJGRkYSEhFClShWjZxO0bNkSLy8vPv/8cyIjIylTpgz79u3j2LFjLF68GJVKxe3bt1m5ciVeXl507NiRy5f124S9zPWEMCVDyx5+vx7GwXMhxCemEJ+YgqerPePeqk3Awxi9ZQ6y7EGIVxelTiQyk/c7MltBiJcjiQUhxEt1epCCjCK9Xr16kZSUxJo1a9i+fTsVKlRgzpw5dO7cGYATJ04wbdo0NmzYQKNGjYw6p5WVFatXr2bRokUsWbKEyMhIqlatypIlS2jbti2Q2gFCq9Vy7do1+vbtm+EcmzZtokGDBtn3RIXIIWnLHkb38KJCaVs2HLzFhduPMVOlzghrXd+FU1ceEfAwhqBHMZkue5DEghAvZ73frUz3yWwFIV6OSkkr3S1MLn1FcigcFa4Lkvz6en287LRRRRmtLMwKXO0EY14zqUheOEhXiPwvP79m1wIiWLL7KlqtQopGoVLZ4jx++owxPWvh6WrPzeCnBW5mgnSFEM8z9D44vez+HR82+xiGriidIF5Pfr4XFybZfQ+WGQtCFELGLnsAdPtlyYMQQrweg50egiI5fS2Muw+iSfr/AnIt6pajdMmiemNlZoIQ2SskPNZgUgFktoIQr0ISC0IUMi+77GHeGJ/cCk0IIQq055c8eLrac/LyQ34+dBuNVsHRzpqiRcxpW7+CFGQUIhcszaTFpI9XGYZ1rZHL0QiR/0liQYhCxlC3B0OsLMzo1cI9FyISQojCIW3WwbI9V3EqUZSgsFiKWJrRsm45zt38hw961ZaCjELkgih1Io+jEwzuu3zvSS5HI0TBYGbqAIQQuetFMxTS/lvQaikIIYSpKYpCRHQCSclagsJiqVSmOPPG+OBgZ51pQUYhRPbbceJepvtkGYQQr0ZmLAhRgD1fS8G+eBFKlyya6VhZ9iCEENknfT2FB4/VrPzlOqFP4jAzU+Fbpxx//f2Y+/+o6dTYNcPxsuxBiJxz6a7hWQnlSxWTFpNCvCJJLAhRQKWvpfA0NpGnsYm4lLIhPOoZyc8th5BlD0IIkb3S6ikM6+LJ7ZAoDv0RggJYWpgx/q3a1KjkQKMazrLkQYhcFqVOJD4hJcN2czMVHhXk91CIVyWJBSEKqMxqKTxLSmFIp+p6XSGk24MQQmQvT1d7OnhX4D87r6AoYGluRkPP0vjUKiudHoQwoX2nAw1u12gV7oZG53I0QhQcklgQooB4ftmDg10RIjOppRARk0iTmmUkkSCEENnAUAvJv27/w94zQYSEq7Etaon6WTIdG1Wkp69bhuNlyYMQuetaYKTB7bIMQojXI8UbhSgA0pY9pBVmzCypAP8WaBRCCPH60pY83Ax+iqIo7Dp5jyW7r/HgHzVNazoD0K1pJY5fDOVm8FMTRyuEiHuWnGGbubksgxDidcmMBSEKAGkhKYQQpqFrIbn7KjZFLfnn6TPKOhajo3dFtp+4p6ufIC0khTC9kPBY4hM1GbZrNLIMQojXJTMWhCgApIWkEEKYhqIoRMUmkpis4Z+nz6hZyYFvhjVC/SxZWkgKkces2ns9wzZzcxWt6pWXZRBCvCaZsSBEPpO+lkLlMsUzHSstJIUQIvukr6cQE5/Esl1X+ftBNGZmKlq/UZ4/bv7D7ftR0kJSiDwmSp3Iw4j4DNtltoIQ2UMSC0LkI+lbSEbGJBIZk0jpktY8VSdJC0khhMhBafUURvfwIiExhZ/23yA+UYOl+f+3kKzsQP1qpWXJgxB5kKFuECoVtKxbnoEdqpkgIiEKljy/FCIuLo4ZM2bg4+NDvXr1eP/99wkKCnrhcVevXmXgwIHUq1ePZs2asWDBApKT9Yu1BAUFMWrUKBo0aECjRo348ssvUavVemOePHnCpEmTaNSoEfXr12fixIk8fvw40+vOmTOHIUOGvMpTFeKFMquloNEqDOlUXZY9CCFEDvJ0tWdYF08WbrvE4l1XSUjSUL+aEx/1qUONyg66MbLkQYi852ZwVIZtigJ/35eiqkJkhzw/Y2HChAlcvXqVKVOmYGNjw5IlSxg0aBD79++neHHDU8CDg4MZMmQI9erVY9GiRdy7d4+FCxeiVqv54osvAIiOjmbw4ME4OTkxZ84cIiIimDdvHmFhYaxcuRKAlJQUhg0bRnx8PF999RUpKSnMnz+f4cOHs3PnTiws9H98GzduZM2aNTRp0iRnfyii0MqsloK0kBRCiJwX+CiGLUfvkKJRAOjU2JXeBmaGyZIHIfIeM1XGbdINQojsk6cTCxcuXMDf35/Vq1fj6+sLQIMGDWjTpg1btmxhxIgRBo9btWoVxYsXZ9myZVhZWdGiRQusra2ZOXMmI0eOxNnZmU2bNhETE8OePXuwt0+9oTg7OzNixAguX75MnTp12L9/P7du3eLAgQO4u6e+cfD09KRr164cPnyYzp07AxAeHs7cuXM5cOBApskOIV7F8/UUihezzHSctJAUQojs9Xw9Ba2icOiPEHaeuIeFuRlFi1jQtr4Lxy+GUqOSgyQRhMjjpL6CEDnP6KUQW7Zs4enT3J0qdPr0aWxsbPDx+bf4nIODAw0bNuTkyZNZHteqVSusrKx02zp27IhGo+HUqVO6MQ0bNtQlFQCaNWuGjY0N/v7+ujFVqlTRJRUA3eO0MQALFy7kxo0brF27Fk9Pz9d/4kLwbz2FtFkKsfGpS3nM06XcpZaCEEJkv7R6CudvhbNg6yW2H78HKjAzU/FBr1r09HVjdA8vlu+5xs1gmUotRF6WWX0F6QYhRPYxOrEwY8YMmjdvzogRI9i7dy/x8RmzftktICAAV1dXzM3N9bZXrFiRwMCMNwiAZ8+e8ejRIypXrqy33cHBAVtbW91xAQEBGcaYm5vj4uKS5RhD1x8+fDj79++ncePGL/8khchEZvUUrK3MpJaCEELkME9Xezo2qsjyPde5FfyUIpbmNKtVlg9715YWkkLkM1JfQYicZ/RSiAMHDvDrr7/i5+fHlClTsLa2plWrVnTt2hVfX18sLTOfpm1ISkoK+/fvz3R/qVKlUKvV2NraZthnY2OTochimtjYWIAXHhcbG2vUmCpVqhgcExwcrHtsaMyrcHTMGA+Ak5Msr8hPsuv1yqyeQnyChv9+2zVbriFSye+YEOJ5Wq3C3tOB7DsdRPFilsTGJ9O+YQV6+rplGCv1FAqXlJQUrl69yqNHj/D29sba2hqNRkOJEiVMHZrIgtRXECLnGZ1YcHNzY9y4cYwbN44bN26wf/9+Dh48iJ+fH3Z2drRv354uXbrQuHFjVCoDv73pJCYmMmXKlEz3e3t7Z5msMDMzPNlCUVILKhmKQVEUveOya0x2iYhQo9UqetucnIrz+HFstl9L5IxXfb2er6XgYFcEFyebTMc62BWRvxPZyJjXzMxMlWniTwiRvz1fSwEgJi6JhdsuERyupmYle4LD1XRrWonjF0OpLkmEQs3Pz49vv/2WiIgIANasWUNSUhLjx4/ngw8+YPjw4SaOUBgi9RWEyB2vVLyxRo0a1KhRg48//phbt26xatUqdu7cyc6dO3FycqJnz570798fZ2fnTM9hY2PD7du3s7zOuHHjePDgQYbtcXFxBmcbwL8zFQzNaIiPj9cVV7S1tTU4Ji4ujvLly79wTGbXF+JlpdVSSFv2EBmTSGRMIuUcivIkJlFvOYTUUxBCiOyVVkthdA8vLMxV/GfHFeISUmhQ3YlbwVGM7uGFp6s91V3tdeMkuVD4nDp1ikmTJvHGG28wfPhwZs+eDYCLiwseHh7Mnz8fJycn3nzzTRNHKtLbdzoQczMVmue+vDM3U+FbpxwDO1QzYWRCFCyv9LW7oij88ccfzJgxg1GjRnHgwAHs7Ox4++23ad26Nf/973/p2LEjx48ff63gKleuzP3793WzENIEBwcbrH0AqQkLZ2dnvaUKABEREajVat1xlStXzjBGo9Hw4MGDLMcAhISEZHp9IV5WZrUUElO0DO5UXeopCCFEDvJ0tWf0mzX5z44rzNr4F/GJKQzuWI3KZez0kghST6FwW7p0KV5eXmzYsEEveeDu7s7mzZupV68e69evN2GEIjN3Q2P0kgoAGq3MVhAiu71UYuH8+fN8/fXXNG/enMGDB7N7927q1avHsmXLOH36NF9//TVfffUVhw4dwtHRkW+++ea1gmvWrBkxMTGcOXNGty0yMpILFy7QtGnTTI/z8fHh+PHjJCUl6bYdOnQIc3NzvL29dWPOnTtHVFSUbsypU6eIj4/XnbtZs2bcuXOHgIAA3Zi7d+8SEBCQ5fWFeBmZ1VKIiEmkSc0yzBvjw5qprZk3xkeSCkIIkc2SkjWcuhpGYrIGgA7eFWlRtzydGrtmmJng6WpPp8aupghTmNjNmzfp0qWLwaWwFhYWdO3aNdPC4sK0hnXJ2LHN0sKMiX3qmCAaIQouo5dCNGvWjIiICMzNzWncuDFdu3alXbt22NhkXAtesmRJatWqxe+///5awTVs2BBvb28mTpzI5MmTKVmyJIsXL6Z48eK88847unF3794lKSmJGjVqAP92aRgxYgSDBw8mKCiIBQsW0KdPH8qVKwdA//792bhxI0OGDGHs2LFERUUxb948fH19eeONNwDo3LkzK1asYPjw4UyaNAlFUZg/fz5Vq1alU6dOr/XcROGVvp6ChbmKFI2SYVzaTAUhhBDZI309hciYBOZtuUj402dYWZjR3rsCJy4+pJaboyx3EHosLS1JSUnJdH9UVNRLFzIXuWPV3usZtmkVhb2ng2QphBDZyOjEgouLC6NGjaJz5844ODi8cPz777/PxIkTXys4gCVLljB79mzmzp2LVqulfv36LFq0SK/67owZMwgNDeXYsWNA6rS0NWvWMHfuXMaNG4e9vT3vvfceH374oe4YBwcHNmzYwHfffcfkyZOxsbGhY8eOegUlraysWLt2Ld9++y2fffYZVlZW+Pj4MHXqVCwsXqk8hSjkDNVTgNReys+v+JFaCkIIkf3S11P4YccV4hNSsLIwY/zbdf6/w4OD1FIQGXh7e7Njxw4GDBiQYd8///zD5s2bqV+/vgkiE1mRwo1C5B6Vkr6AQSb27NlDgwYNcHFxMbj/3r17HDlyhJEjR2ZrgIWJdIXI/170en287LTBpQ821uZYW1kQEZOIo10RerVwl2UPuUS6Qog06e/Bcv/Nf4x5zW4GP+U/O66QlKwBFdT3cKLVGy56SYSbwU8JehQjyx5ymLG/Y3nhHnzv3j369u2Lo6Mjvr6+bNy4kXfffRdzc3N2795NUlISW7ZswdMz47R7YRxD74PTe9n78s+HbnHy8iMp3JjL5N/P/CG778FGf+0+bdo05s2bl2li4Y8//mDp0qWSWBAiC5nVU4hL0LD4oxa5HI0QQhQuWq3CpTtPdPUUOnpXoE+rqhnGeUpbSZGOu7s7mzZtYubMmfz8888AbNy4EQAvLy8+++wzSSrkQVK4UYjck2li4f79+3z22We6jgyKorB8+XK2bduWYayiKNy6dQsnJ6eci1SIfE5RFIoWseBZYsY1mlJPQQghclZCUgqr9t7g0t0nWFqY0aFhBU5cekgtt1KSRBBGqVatGj///DNRUVGEhISg1WopX768vP/Nwyb0qcOkpaf1lptK4UYhckamiYUKFSrg6urKqVOnAFCpVERGRvLs2bMMY83NzXFzc9OrYSCE0C/UWMTSjMRkLWYq0Eo9BSGEyDHpizQ+jU1kzuY/+edpAkUszRn3Vu3UWQmVpJ6CeHklS5akZMmSJCcnc/r0aczNzWnSpInU38qDdpy4R/pF31K4UYickeUd8Ouvv9b9f/Xq1Zk+fTrdunXL8aCEKAjSF2pMTNZibqaiee2yXA2IkHoKQgiRQ54v0mhjbcH3/72E+lkyXpUd9NpIerraM7qHF0GPYiSxILKUlJTEzJkzefDgAWvWrCEpKYm+ffty69YtIHWpxPr163F0dDRxpOJ5V+49ybBNCjcKkTOMTq2m3TiFEMbZ5X9Pl1RIo9EqXA2IYN4YHxNFJYQQBV9awmDJziskpWjRKgpDOlXHt045g2MlqSBeZMmSJWzbto3evXsDqUXNb968yaBBg/D09GT27Nn88MMPel/KCdOKUieSmKz/PszSwoy5o5pQwlaWoAqR3TJNLKTvArFnzx6jTtijR4/siEuIfC+zQo2ZbRdCCJF9IqITeJaUWqSxXUMXg0kFIYzl5+fHW2+9xcyZMwE4dOgQxYsXZ8qUKVhYWHD//n22b99u4ijF8/adDszQZUKrlWUQQuSUTBMLU6dO1esCMXXqVFQqFVl1p1SpVJJYEILUVmUqwNBvixRqFEKInKMoCgd+D2anfwDmZiraN6zA/648om4VJ5mZIF5ZWFgYdevWBeDZs2ecP3+eli1b6uoqlC1blpiYGBNGKNKTjhBC5K5MEwsbNmzA3d1d77EQwrC0Io2RMYnYFLUkLiGZErZWxCWkkPzccggp1Cjygi1bttCxY0fs7eVDlsj/ni/UqNUqbDlyhyN/PsBMBePfro1XZUe83BylSKN4LaVKleLJk9T1+v/73/9ISkqiZcuWuv23b9+mdOnSJopOGDKhTx0+XnZGL7kgHSGEyDmZJha8vb0NPo6NjcXW1haVSgXAvXv3cHJyws7OLgfDFCLvSl+kUf0sGZUKujerTBFLc11XCCnUKPKKGTNm8O2339K0aVO6du1K27ZtKVasmKnDEuKVpBVqHNGtBucP3uZ/l0IxM1PxVgs3vCqnFtKTIo3idTVq1Ij169dTpEgRNm3aRNGiRWnbti0xMTHs3LmTbdu20a9fv2y/blxcHN9//z2HDx8mPj6eBg0a8Omnn1KpUqUsj7t69Spz587l2rVr2NjY0KtXLz788EMsLS11Y4KCgpg9ezYXLlzA3Nycjh078vHHH2Nra6sb8+TJE2bNmsWpU6dISUmhRYsWTJs2LdMWm3PmzOHmzZusW7cuO57+a9l3OjDDjAXpCCFEzjG6eKNWq2Xu3Lls3ryZPXv24ObmBsDy5cs5ePAgo0aN4oMPPsixQIXIqwwVaVQU2H8miHljfCSRIPKcAwcO8Ouvv+Ln58eUKVOwtramVatWdO3aFV9fX703nkLkdZ6u9gzr4skPO66g0SpYWZoxvndtPCs5ZBgnSQXxqqZPn054eDhz5syhWLFizJw5Ezs7O/7880/mzJlDw4YNc+R98IQJE7h69SpTpkzBxsaGJUuWMGjQIPbv30/x4sUNHhMcHMyQIUOoV68eixYt4t69eyxcuBC1Ws0XX3wBQHR0NIMHD8bJyYk5c+YQERHBvHnzCAsLY+XKlQCkpKQwbNgw4uPj+eqrr0hJSWH+/PkMHz6cnTt3ZmivuXHjRtasWUOTJk2y/efwKm7fj8qwTTpCCJFzjE4s/PTTT6xbt45u3bpRokQJ3fahQ4dibW3N0qVLKVWqVI5ka4XIy6RIo8hv3NzcGDduHOPGjePGjRvs37+fgwcP4ufnh52dHe3bt6dLly40btxYNztNiLwqPiGF/WeDdd9MdmhYMUNSQYjXZWdnx9q1a4mMjMTW1hYrKysAPD092bp1K3XqZP/0+gsXLuDv78/q1avx9fUFoEGDBrRp04YtW7YwYsQIg8etWrWK4sWLs2zZMqysrGjRogXW1tbMnDmTkSNH4uzszKZNm4iJiWHPnj26ZXHOzs6MGDGCy5cvU6dOHfbv38+tW7c4cOCAbnm0p6cnXbt25fDhw3Tu3BmA8PBw5s6dy4EDBzJNdphCtQolefgkHgBzMxW+dcrJTAUhcpCZsQN37NhBr169mDdvnl6P3ho1ajBz5ky6d+/Opk2bciRIIfIqraJQxNLc4D4p0ijygxo1avDxxx9z9OhR9uzZQ7Nmzdi5cydDhw6lRYsWLFy4kPDwcFOHKYRB6mfJzPvvRQIeRmNtZU7fdh4cvxjKzeCnpg5NFFAODg66pAJAsWLFciSpAHD69GlsbGzw8fm3RbWDgwMNGzbk5MmTWR7XqlUrvTg7duyIRqPh1KlTujENGzbUq7XTrFkzbGxs8Pf3142pUqWKXs21tMdpYwAWLlzIjRs3WLt2LZ6enq//xLNBlDqRU1fDdI81WoVTVx8RrZYvfYTIKUbPWAgLC8vyxvnGG29w8ODBbAlKiLwsrVBjREwiVpZmJCVrMTNT6bU0kiKNIr9QFIXz58/j5+fH8ePHCQsLo0SJEnTo0AEzMzP++9//smHDBhYsWECrVq1MHa4oxJ4v0ggQrU7k25//JCI6AWsrcz7oXRvfBhWpWMpGCjWKHHHy5En27dvHkydP0Gg0GfarVCrWr1+fbdcLCAjA1dUVc3P9LzAqVqyIn5+fwWOePXvGo0ePqFy5st52BwcHbG1tCQwM1J27e/fuemPMzc1xcXHRG5P+PGnXTxsDMHz4cNzc3DAzM2Pp0qUv/0RzgLSaFCL3GZ1YKFOmDH/99Rd9+/Y1uP/q1at6MxmEKIjSF2pMStb+//S6sly5F0FkTCIOUqRR5ANpyYTDhw8TERFBkSJFaNWqFd26dcPX11e3dvajjz7irbfe4ptvvpHEgjCptCKNo3t44WxflG83XOCpOonabg50aOSqSyJIoUaREzZt2sTMmTMBcHR01JsN8CpSUlLYv39/pvtLlSqFWq3WK6SYxsbGBrVabfC42NhYgBcel1aM/UVjqlSpYnBMcHCw7rGhMa/C0TFjPIY4Ob14uUVQuNpgq8mg8FijjhevT37O+UN2vk5GJxa6du3KsmXLcHd3591338XGxgZIzYxu27aNXbt28f7772dbYELkRYYKNWq0ClfuRTBvjA9OTsV5/DjWRNEJYZxmzZoRERGBubk5jRs3pmvXrrRr1053X39eyZIlqVWrFr///rsJIhXiX2kJg2W7r6LVKjxL0tC/rQdtG7gYHCtJBZGdNmzYQPXq1Vm9ejWlSpV67fMlJiYyZcqUTPd7e3tnWUjXzMzwamZFSf0wbag+jqIoesdl15jsEhGhzjDLID1j32d92KuWwVaT43rVkvdpuUDeD+cPxr5OZmYqoxJ/RicWRo0axZUrV1iwYAE//PADDg4OmJmZ6aaD+fj4MHbsWGNPJ0S+JIUaRUHg4uLCqFGj6Ny5Mw4OLy5y9/777zNx4sRciEyIrDmVsNYlFZrXKWswqSBETnj06BHTp0/PlqQCpH7rf/v27SzHjBs3jgcPHmTYHhcXZ3C2Afw7U8HQjIb4+HhdcUVbW1uDY+Li4ihfvvwLx2R2/bxCWk0KkfuMTixYWlqyevVq/P39OXHiBA8fPkSj0dCiRQt8fX1p06aNVA8XBVpyihZLCzOS081YACnUKPKXfv360aBBg0yTCvfu3ePIkSOMHDkSSC3wKISpPYl+xswNF3iWpMG3Tjn++vsxjWs8lZkJIldUrFiRJ0+e5Oo1K1euzNmzZ1EURe89dnBwsMHaB5CasHB2dtZbqgAQERGBWq3WHVe5cuUMYzQaDQ8ePKBDhw66MX///XeGa4SEhORYwcrsIq0mhch9Lz2PqUWLFnz55ZesXLmSH3/8ka+//pq2bdtKUkEUaCkaLcv3XCM5JbWmwvOkUKPIb6ZNm8alS5cy3f/HH3/kmQJcQgBERCcwc8MFYuKTGdjBgyGdqjO6hxfL91yTDhAiV4wYMYKff/6ZO3fu5No1mzVrRkxMDGfOnNFti4yM5MKFCzRt2jTT43x8fDh+/DhJSUm6bYcOHcLc3Bxvb2/dmHPnzhEVFaUbc+rUKeLj43XnbtasGXfu3CEgIEA35u7duwQEBGR5/bygWoWSuv83N1fRql551kxtzYyh3qYLSogCLtMZC+fPn8fd3V33jdb58+eNOmHDhg2zJzIhTOz57g9pMxUGtPegaBEL3XZHKdQo8oH79+/z2Wef6dbeKorC8uXL2bZtW4axiqJw69YtnJyccjtMIXSe7wARGZPA3C1/EfcshXoepWhVL3X5gxRpFLnpzz//xMbGhjfffJPKlSvj4OCQ4Uu17O4K0bBhQ7y9vZk4cSKTJ0+mZMmSLF68mOLFi/POO+/oxt29e5ekpCTd7LLhw4ezf/9+RowYweDBgwkKCmLBggX06dOHcuXKAdC/f382btzIkCFDGDt2LFFRUcybNw9fX1/eeOMNADp37syKFSsYPnw4kyZNQlEU5s+fT9WqVenUqVO2Pc/slqHVpCa11WR3n0qUsJUZpkLklEwTCwMHDmTevHl069ZN9zirWQlp07Ru3ryZ/VEKkcvSd39Im6lQtIgFTWqWkUSCyFcqVKiAq6urrn+5SqUiMjKSZ8+eZRhrbm6Om5sbH374YW6HKYROWgeIge092OkfQHRcEkUszWlbv4LeOCnSKHLL//73PyC1S9qzZ88IDQ3NlesuWbKE2bNnM3fuXLRaLfXr12fRokWUKFFCN2bGjBmEhoZy7NgxANzd3VmzZg1z585l3Lhx2Nvb89577+nd1x0cHNiwYQPfffcdkydPxsbGho4dO+oVlLSysmLt2rV8++23fPbZZ1hZWeHj48PUqVN1nYPyon2nA9EqUl9BiNymUhTFYPnV3bt306BBAypUSP1HfNeuXUYtd+jZs2f2RliIGKqGK1VVTePjZacNFmR0tCvCvDE+mR4nr1f+Y8xrZmw13PyievXqeoljkSr9PVh+n03rz9v/sGzPNcxUKopYmjO2V60XJhHkNctfsrsiucjfsqsrxJdr/uD+PxmLTlYobStLIXKJ3Ivzh1zrCpE+QdCrV68sT6TRaHj48OELLyhEfiDdH0RBduvWLVOHIESW4hOS2Xc6CBWpLX3b1HeRmQkiT9BoNFy7do3Q0FCsrKwoW7YsNWvWNHVY4jkzhnqzYOslrgVGAqk1Fnxrl5PZCkLkMKPnMXl6ejJv3jy6du1qcP/u3buZNWsWf/75Z7YFJ4SpWFuZk5CkybBduj+I/GjPnj00aNAAFxcX3WNj9OjRI+eCEiITickaFu24woPHaopYpS5/OH4xlOqy7EGY2PHjx5kxYwbh4eG6mjUqlYrSpUvz5Zdf0rp1axNHKCC1xsKN54q6So0FIXJHpomF8PBwzp49q3usKArnz58nJSUlw1itVsu+ffukM4TIt54v1Fi0iAUJSRrMVCq9NXrS/UHkV1OnTmXevHm6xMLUqVNRqVRkshIOSH2z/DKJhV9//ZXly5dz//59ypcvz8iRI40+fs6cOdy8eZN169bpbb9w4QLvvvtuhvEtW7Zk5cqVGbafOHGCsWPHcv36daPjFqb3fKHGFI2WpbuucvdBNBbmKj7oVRtPV3uqu9qzfM81RvfwkuSCMIkLFy7w4Ycf4ujoyIQJE3B3d0dRFAICAti8eTPjxo1jw4YNusKHwnT2nQ7MsKRCaiwIkfMyTSw4ODiwYsUKgoKCgNQ3mVu3bmXr1q2ZnmzgwIHZHqAQOS19ocZniSmYqcC3TlmuBkRI9weR723YsAF3d3e9x9nJz8+PyZMnM2jQIJo3b86RI0f45JNPsLa2pmPHjlkeu3HjRtasWUOTJk0y7Lt9+zbFihVj7dq1etvt7OwyjP3rr790VctF/pJWqHFU95qcuPyQa4GRmJup6N3CXZdEkA4QwtQWL15M+fLl2bFjB8WLF9fb179/f3r37s3y5ctZvXq1iSIUaW7fj8qwTaNRuBsanfvBCFGIZJpYsLS0ZM2aNTx48ABFURg8eDAjR47Exydj4TozMzMcHBxwc3PL0WCFyAm7/O/pkgpptApcDYjIslCjEPlFWt/y9I9jY2OxtbXVzTa7d+8eTk5OBj+4Z2XBggV06tSJ6dOnA9C8eXOio6P54YcfMk0shIeHM3fuXA4cOJDhTXqaW7duUbVqVerWrZvptRMSElizZg3Lli3DxsbmpeIWeYOnqz2j36zJou1XSNZoKWJpxri36mRIIEgHCGFKV65cYezYsQbvV7a2trz11luSVMgjqlUoycMn8YDUVxAiN5lltbNcuXJ4e3vTqFEjZs2axVtvvYW3t3eGPw0aNJCkgsi3pFCjKGy0Wi2zZ8/Gx8eHwMBA3fbly5fTtGlTlixZYvS57t+/T0hICO3bt9fb3qFDBwICArh//77B4xYuXMiNGzdYu3Ytnp6eBsfcvHmTatWyfjN44MAB1q9fz7Rp0xgwYIDRcYu85WbIU5I1qQne9g0rSgJB5DsqlYrk5GRTh1HoRakT+d+VMN3jtPoK0Wp5TydETssysfC8nj17UqFCBb2+50+fPmXTpk1s2bKFqKionIhPiBxXwsbK4HYp1CgKqp9++ol169bRoUMHvV7oQ4cOpUePHixdupT//ve/Rp0rICAAgMqVK+ttd3V1BdBLXDxv+PDh7N+/n8aNGxvcr9VquXPnDmFhYfTs2RMvLy9atmzJmjVr9JY7NGrUiKNHjxqsxSDyhxMXQ/n1TDCW5mZ0a+rK8Yuh3Hyu8JoQeUGdOnXYsWMH8fHxGfap1Wq2b99OrVq1TBCZeN6+04FoNOlnoabWVxBC5Cyju0LExMQwYcIEYmJi2L59O2q1mt69e/Po0SMURWHZsmVs3ryZChUq5GS8QmSriOgEklIydn+QQo2iINuxYwe9evXiu+++09teo0YNZs6cSXJyMps2baJfv34vPFdsbGr/Y1tb/f7GacsS1OqMvcQBqlSpkuV5AwMDSUhIIDAwkIkTJ2Jvb8/Ro0eZO3cuarWacePGAVC+fPkXxmgsQz2anZwML9MQ2eOP62H8fOg2FuYqvhjWiLrVStO4TnnmbLjAJ4MaULuK00ufU16z/CW/vF4ffPABgwYNomvXrgwYMIBKlSoB6Io3hoeHM2PGDNMGKbh9P4r0lXakvoIQucPoxMKiRYs4d+4cI0aMAFLfmD58+JApU6bg5eXFxx9/zKJFi5g/f36OBSvE63q++4N98SL//82nit4t3DhxMVQKNYpCISwsjDp16mS6/4033uDgwYNGnev5lmuGtpuZGT0xTo+zszOrV6/G09MTJ6fUD5dNmjQhISGB1atXM3To0AzJjNcVEaHWqyTu5FScx49js/Ua4l8BD2OYu/kvStpaMahjdco7FOXx41jKlrBmZPeaXLoZTtkS1i91TnnN8hdjXy8zM5XBxF9uatCgAYsXL+brr79m7ty5unueoig4OTmxcOHCTGdgidxTrUJJwiLj0WrB3EyFbx2pryBEbjE6sXDs2DEGDBig+5boyJEjODo6MnToUADefffdDJW7hchL0nd/eBqbut6uaxNXujSpRJcmlUwYnRC5p0yZMvz111/07dvX4P6rV6/i6Oho1LnSCpmln5kQFxent/9l2dra4uvrm2F7y5Yt2b59O4GBgTLtOJ95vq1k+NN4Fm2/TNEi5vjULkedKqX0xkqhRpEXtWnThpYtW3L9+nUePHgApM6aqlmzJhYWRr+lFjkkSp3IqathaP9/JYRGm1pfobtPJUrYyvJWIXKa0V8lRUREULVqVSB16uulS5f0OkTY29vr1V8QIq8x1P0BUhMOQhQmXbt2Ze/evaxatUqXAAB49uwZ69evZ9euXXTr1s2oc6XVVggJCdHbHhwcrLf/Zd2+fZvNmzdnKIaWkJAApP6bI/KXtLaSf97+h4VbL5Oi0ZKsUSSBIPIVc3NzXFxcqFixIpUrV8bNzU2SCnnEvtOBaNO1HJb6CkLkHqPvhM7Ozrrq3keOHEGj0dCyZUvd/r/++ouyZctme4BCZBfp/iBEqlGjRnHlyhUWLFjADz/8gIODA2ZmZjx58gSNRoOPjw9jx4416lyurq64uLhw8OBB2rVrp9t++PBhKlWqRLly5V4pxuDgYGbMmIGzszNt2rTRbT9w4AAuLi7ZWltB5A5PV3tGdPNk0Y6roChYW1kwtmctSSyIfOPChQt8//33XLlyRbfcy9zcnMaNGzNlyhQ8PDxMHGHhdjc0Bo1GP7Eg9RWEyD1GJxZatWrF+vXrUavV7N+/nxIlStC6dWvCw8NZvXo1v/zyC2PGjMnJWIV4LY52RQwmEaT7gyhsLC0tWb16Nf7+/pw4cYKHDx+i0Who0aIFvr6+tGnTJkPNhKyMHTuWadOmUaJECVq2bMmxY8fw8/Nj4cKFAERGRhISEkKVKlWMrovQsmVLvLy8+Pzzz4mMjKRMmTLs27ePY8eOsXjx4peKT+QNiqJw9nq4ro5Fm/ouklQQ+ca5c+cYNmwYxYoVo3///lSqVAmNRkNQUBD79u3jnXfeYcuWLZJcMKEZQ72JUicyaclpFMDSwoy5o5rIMgghconRiYWPP/6YZ8+esWPHDpydnfnqq6+wtrbm77//ZtOmTXTv3l1X2FGIvKipVxn2nQnW2ybdH0Rh1qJFC1q0aPHa5+nVqxdJSUmsWbOG7du3U6FCBebMmUPnzp0BOHHiBNOmTWPDhg00atTIqHNaWVmxevVqFi1axJIlS4iMjKRq1aosWbKEtm3bvnbMIvftPxvM2evhWFma0aFhRY5fDKW61FIQ+cSiRYsoX748W7ZswcHBQW/f2LFj6dOnDwsWLGDFihUmilBA6nKItDkLacsgpHijELlDpShK+q4sLyUpKYno6Ghd1W7x6tJXJAepcJ1dnkQ94+v1FzBTgbm5GU9jc6b7g7xe+Y8xr1leqEj+Os6fP4+7u7vuzfD58+eNOq5hw4Y5GVaeI10hcs6ft/9h6e5rWJirmPB2HTwrOXAz+CnL91xjdA+vbEsuyGuWv+SnrhD16tVj/PjxDBkyxOD+H3/8keXLl/Pnn3/mbmAFiKH3well9XcmSp3IJyvOkvxcPS2ZtWAaci/OH7L7Hvza1WasrKx0SYUbN25Qo0aN1z2lENni+daS5mYqzFQwY1gjyjgUM3VoQuSqgQMHMm/ePF1BxoEDB2a5lEBRFFQqFTdv3sytEEUB83wHiOCwWFb/eoPixSypW8URz0qpCS5PV3tG9/Ai6FGMzFoQeZ6joyMRERGZ7k9MTMz2Frji5ew7HZghMSGzFoTIPUYnFpKTk1m1ahWHDx8mPj4erfbfbKBGoyEuLg61Wi1vREWekL61pEaroDJXEfgoRhILotCZNWsWdevW1T3+7rvvpEaByFFpHSAGdvDgv0fvUsTSHK1WoXFN/SLP0lZS5BejRo3im2++oV69erRu3Vpv3+XLl9mwYQOTJk0yUXQC/r94o1aKNwphKkYnFhYtWsRPP/1EmTJlsLOz4++//6ZBgwY8fvyY0NBQrK2t+fTTT3MyViGMZqi1ZIpGYZf/vWxd+iBEftCzZ0+9x7169cpyvEaj4eHDhzkZkijgPF3teb+rJz/svIoKKGJpzthe0gFC5F+XLl3C0dGRsWPH4ubmhru7O5aWlty/f5+rV69iZWXFr7/+yq+//qo7RqVSsX79ehNGXbhM6FOHj5edQaNVsDBXMW90U1kCIUQuMjN24MGDB/H29ubYsWOsXr0agC+++IJDhw6xcuVKUlJSsLS0zLFAhXgZ0lpSiMx5enrqvflNb/fu3fTo0SP3AhIFjqIo/HHrH7RaBY1WkQ4QIt87c+YMAGXLluXZs2dcu3aNixcv8uTJE8qWLYujoyMPHjzQ+5PWpl3kjn2nA3VtQBUF9p4OMm1AQhQyRs9YCA8P57333sPMzAxnZ2ccHR25ePEiHh4etGjRgp49e7Jt2zb69u2bk/EKYZQSNlZExyVl2C6tJUVhFB4eztmzZ3WPFUXh/PnzpKSkZBir1WrZt2+fLJUQr+XExVBOXw3DysKMDt7SAULkf8eOHTN1CCILUepETl0NI20lhEarcOrqI7r7VJJZC0LkEqMTC9bW1nozEipWrMjff/+te1y7dm0OHTqUvdEJ8QriEpL1aoCkkdaSorBycHBgxYoVBAUFAanTc7du3crWrVszPWbgwIG5FJ0oaO4+iGbTb39jbqZiXO/a1KjsQHVX+2zvACFEXnDnzh3MzMxwd5f3F6a073QgWkUKNwphSkYnFjw9PTl58qRuRoKbmxsXL17U7Q8PD5dvuITJaRWFVXtvEJ+oobtPJU5ffURETM60lhQiv7C0tGTNmjU8ePAARVEYPHgwI0eOxMfHJ8NYMzMzHBwccHNzM0GkIr+LUieydM9VihaxYGgXT2pUlg4QomBQFIXVq1cTGBjIrFmz0Gq1jBo1iv/9738ANGnShMWLF2NjY2PiSAunu6ExaDRSuFEIUzI6sfDuu+/y0Ucf0b9/f1atWkWXLl3YuXMn06ZNw83NjXXr1ulVHRciNz3fWhLAp1YZejR3o0dz+XAkBEC5cuUoV64ckNolokGDBlSoUMHEUYn87vm2kikaLcv2XCPuWTKt6rlQr6qT3ljpACHys59++okFCxbQvHlzAPz8/Dh58iQdOnSgatWqrF69mqVLlzJlyhQTR1o4zRjqTZQ6kYlLTgNgaWHG3FFNZBmEELnI6OKNHTt25JtvviEqKoqiRYvStGlT3n//fXbv3s38+fOxs7Nj2rRpORmrEAaltZZ8vjDj+Zv/cPZ6mAmjEiLv6tmzJxUqVODZs2e6bU+fPmXTpk1s2bKFqKgo0wUn8pW0tpI3g5+y9ehd7j6IxsLcjLpVS5k6NCGy1e7du2nXrp2ugPmBAwcoWrQoc+bM4YMPPqB///4cPHjQxFEWbvtOB+r+P20ZhBAi9xg9YwHg7bff5u2339Y9njRpEu+88w7R0dG4u7tjZWWV7QEK8SKGWksmpWiltaQQmYiJiWHChAnExMSwfft21Go1vXv35tGjRyiKwrJly9i8ebPMaBAvlLbEYfHOKyQkabC0MOPD3rVlZoIocO7fv8+QIUMASE5O5uzZs3h7e2NtbQ2Au7s7T548MWGEhVta8cY0Go0UbxQitxk9YyEz5cqVw9PTU5IKwmSktaQQL2fRokWcO3dON6V3x44dPHz4kI8//pgNGzZgZmbGokWLTBukyDdK2FiRlKIBoEPDCpJUEAWSnZ0darUagHPnzhEfH4+vr69uf0hICKVKyUwdU9l3OhCt1nDxRiFE7sh0xkKbNm1e+mQqlYojR468VkBCvKwiluYkJmsybJfWkkIYduzYMQYMGMC4ceMAOHLkCI6OjgwdOhRIramzdu1aU4Yo8onEZA2Ltl9G0UK7Bi6cuPQQz0oOklwQBU69evXYuHEj5cuXZ8WKFVhYWNC+fXuSk5M5fvw4W7ZsoW3btqYOs9C6GxqDRivFG4UwpUwTC2lFvoTIy/64GU5isgYzM5VeplpaSwqRuYiICKpWrQpAbGwsly5donPnzrr99vb2evUXhMjMst1XeRKdwFst3enc2JW6VZ2kraQokKZPn86wYcMYN24cKpWKKVOm4OTkxLlz5xg3bhxubm6MHz/e1GEWWhP61OHjZWfQaBUszFXMG91UlkAIkcsyTSz8/PPPuRmHEC8t/Gk86/xu4V7Ojpb1yrPnfwHSWlIIIzg7O3P//n0gdbaCRqOhZcuWuv1//fUXZcuWNVF0Ir84ey2MqwGRNKrhTOfGroC0lRQFV9myZdm7dy83btzA2dkZZ2dnAKpXr86CBQto1aoVRYsWNXGUhde+04EoSuoXTIoCe08HMbBDNRNHJUTh8lLFG7OSkJBAZGSkzHQQOer5tpLmZirMzWDkmzUpVaIoPrXkg5AQxmjVqhXr169HrVazf/9+SpQoQevWrQkPD2f16tX88ssvjBkzxtRhijzm+daSjyLi2HDoNi6lbHBxstEbJ20lRUFlYWFB7dq10Wq1PHnyBDs7O0qUKKE340vkvrTCjWkTVzVaKdwohClkWryxZ8+enDx5Um9bcnIyR48e5enTpxnGHz58+JXqMghhrPRtJTVaBa0Cdx7I+jkhXsbHH39Mly5d2LFjB3Z2dixcuBBra2vCw8PZtGkT3bp1Y8SIEaYOU+Qxaa0lr957wvI911Cp4Kk6EbdyJUwdmhC5Ijg4mA8//JD69evj6+vLn3/+ydmzZ3n77be5cOGCqcMrtPadDkSrSOFGIUwt08TCzZs3iY7W/8CmVqv54IMPuHXrVo4HJkR6htpKpmgUdvnfM1FEQuRPVlZWzJw5k0uXLnHo0CGaNGkCpE7pPXnyJHPmzJFOPyKDtGUOS3Zd48HjOFQqGNOzlsxOEIVCUFAQb7/9Nn/88QfNmzfXTbs3NzcnICCAoUOHcunSJdMGWUjdDY1Bo5HCjUKY2ku3m1TSZQSFyC3SVlKInGVlZYWTkxMAN27cMHE0Ii+KT0ghWZOa4G1bX1pLisJjwYIFWFtbc+DAAb766ivd+2Fvb28OHDhAqVKlWLJkiYmjLJxmDPVmzdTW2NmkJsRb1SvPmqmtmTHU28SRCVG4ZFuNBSFymoNdESINJBGkraQQLyc5OZlVq1Zx+PBh4uPj0Wr/nQmk0WiIi4tDrVZz8+ZNE0Yp8pqnsYn89Ot1zFQqOjWuyPGLoVSXegqikPj999957733cHR0zLAk2NnZmf79+/Pjjz+aKDoRpU4kJi4JQOorCGEiLz1jQQhTqVo+4zpeaSspxMtbtGgRixcvJjo6mqJFixIaGkrZsmWxsLAgLCyM5ORkPv30U1OHKfIQraKweMdlEpK1vNe5Gr1buDO6hxfL91zjZnDGuktCFDRJSUnY2dllut/S0pLERJlBaSr7Tgfq/l/qKwhhGpJYEPlCSHgsF24/xrVMcd0MBUe7IgzuVF3aSgrxkg4ePIi3tzfHjh1j9erVAHzxxRccOnSIlStXkpKSgqWlpYmjFHnJkQsPCApX075BBXxqpXZ/er61pBAFXfXq1Tl27JjBfSkpKezdu5dq1aS9oSlEqRM5dSVM91ijSe0KEa2WRI8QuUmWQog8LylZw8q917EtZsmkvnWxLSofeIR4HeHh4bz33nuYmZnh7OyMo6MjFy9exMPDgxYtWtCzZ0+2bdtG3759TR2qyAMe/KNmx4l71K1Sir5tqujtk9aSorAYOXIkY8aMYfLkybouaKGhoRw9epSffvqJGzdusGjRItMGWUhl1RViYAdJ9giRW7JMLAQEBHD+/Hnd49jYWABu376NhYX+offuSWV+kb3OXg9jl/89XXHGzo0rSlJBiGxgbW2tNyOhYsWK/P3337rHtWvX5tChQ6YITeQBfr8HU6msHZ6u9iSnaFi57zpWlmZUcLZFpVKZOjwhTKJVq1Z8++23fPfdd+zfvx+Azz//HEVRKFKkCJ988gkdOnQwcZSF093QGDRa6QohhKllmVhYsWIFK1asyLB9zpw5GbYpipIjbzji4uL4/vvvdUXGGjRowKeffkqlSpWyPO7q1avMnTuXa9euYWNjQ69evfjwww/13kwHBQUxe/ZsLly4gLm5OR07duTjjz/G1tZWN+bJkyfMmjWLU6dOkZKSQosWLZg2bZqucjrA48eP+eGHHzh9+jRRUVFUrlyZ999/n06dOmX7z6OwOHs9jPV+t/TaSx658IDyTray9EGI1+Tp6cnJkyd1MxLc3Ny4ePGibn94eLh8gCzEKpW1Y/mea4zu4cWlO08IfRxHUStzqleUmQmicOvVqxft27fnzJkzhISEoNVqKV++PE2bNsXeXn4/TGXGUG8WbrvE1YBI6nuUYmyv2qYOSYhCKdPEwgcffJCbcWRqwoQJXL16lSlTpmBjY8OSJUsYNGgQ+/fvp3jx4gaPCQ4OZsiQIdSrV49FixZx7949Fi5ciFqt5osvvgAgOjqawYMH4+TkxJw5c4iIiGDevHmEhYWxcuVKIHXN3LBhw4iPj+err74iJSWF+fPnM3z4cHbu3ImFhQVJSUkMHz6c2NhYxo0bR+nSpTl06BAfffQRGo2Grl275trPqiDZ5X9PL6kAkJSiZZf/PUksCPGa3n33XT766CP69+/PqlWr6NKlCzt37mTatGm4ubmxbt066tata+owhYmk1U5YsusKzxI1WFqY8UHv2rLkQQjA1taW9u3bG9z3yy+/8Oabb+ZyRCJKnciNoNQispfvRRKtTpSOEEKYQJ5OLFy4cAF/f39Wr16Nr68vAA0aNKBNmzZs2bKFESNGGDxu1apVFC9enGXLlmFlZUWLFi2wtrZm5syZjBw5EmdnZzZt2kRMTAx79uzRZZmdnZ0ZMWIEly9fpk6dOuzfv59bt25x4MAB3N1TOw94enrStWtXDh8+TOfOnTl58iS3bt1i+/bt1K6dmiH18fHh4cOHrF69WhILryjCQFvJrLYLIYzXsWNHvvnmG9auXUvRokVp2rQp77//vq6QY7ly5Zg2bZqJoxSmVKlMcdKWLLdr4CJJBVEopaSkcOTIES5fvoyiKNSoUYMuXbpgbm6uNy40NJQvvviCM2fOSGLBBPadDkT5/xuWgtRWEMJU8nRXiNOnT2NjY4OPj49um4ODAw0bNuTkyZNZHteqVSusrKx02zp27IhGo+HUqVO6MQ0bNtSbutasWTNsbGzw9/fXjalSpYouqQDoHqeNsbGxoW/fvtSqVUsvBjc3N0JCQl7j2RdumdVSSOsIIYR4PW+//TYHDhzQvUGeNGkSx44dY/fu3Rw8eFDvvicKn5V7r5OQpKFZrbKcvPxIWkqKQiciIoK33nqLCRMmsHbtWtatW8cnn3xCjx49iI7+d+3++vXr6dq1K6dPn+aNN94wYcSFU5Q6kVNXw0grsSAdIYQwnTydWAgICMDV1TVDZrhixYoEBgYaPObZs2c8evSIypUr6213cHDA1tZWd1xAQECGMebm5ri4uGQ5Jv31mzRpwtdff623Hjk5ORl/f3+qVq36ks9YAMTEJ5GcoiH9Cm8rCzN6tZAPO0LklHLlyuHp6amXlBWFz4Hfg7lyLwLv6qUZ2sWT0T28WL7nmiQXRKGyYMECbt26Rb9+/di2bRv79u1j8uTJPHjwgG+++YakpCTGjh3L7NmzsbKy4ptvvmHTpk2mDrvQyaojhBAid5ms3WRKSoquqq4hpUqVQq1W6xVSTGNjY4NarTZ4XFrnihcdFxsba9SYKlWqGBwTHBycaezff/89QUFBLF26NNMxhjg6ZowHwMnJcC2JgmrNhvNotAoDO3vidzaIJ0+fUcq+KIM6edKyfgVTh/dChe31KggK+muW1hrtZahUKo4cOZID0Yi8LD4hhf1ngnCwK8Kwrp7AvzUXgh7FyJIIUWicPXuW9u3b8+WXX+q2Va1aFWtra77//nu++uorjh49Sps2bfj6669xdHQ0YbSF193QGDQa6QghRF5gssRCYmIiU6ZMyXS/t7e3XgeH9MzMDE+2SFtjZaiiuaIoesdl15jnt8+bN49169YxbNgw2rZtm2n8hkREqNGma5fj5FScx49jX+o8+dmFW/9w6vJDevq60bJ2WVrWLqu3P6//LArb61UQGPOamZmpMk385QflypUzdQgin9h67A4JyRom9auHpcW/swU9Xe0lqSAKlYiICJo0aZJhu6+vL9988w2//PILn332GQMGDDBBdCLNjKHeAIxbdBJ1Qgqt6pWX+gpCmIjJEgs2Njbcvn07yzHjxo3jwYMHGbbHxcUZnG0A/85UMDSjIT4+XtdJwtbW1uCYuLg4ypcv/8Ix6a+flJTE1KlT2b9/P8OGDcsyaSL0nb0exi7/e0TEJKJSpdZR6Ny4oqnDEqLA+Pnnn00dgsij/H4PplJZOzxd7bkaEMH/rjzC27M0t0Oe4lbOztThCWEyiYmJBt9rpm3r37+/JBXyiCh1IuqEFABOXX1Ed59K0hVCCBN46RoLKSkpXLx4kQMHDvDkyRPUarVeEZvsVLlyZe7fv6+bhZAmODjYYO0DSE1YODs7Z1iqEBERgVqt1h1XuXLlDGM0Gg0PHjzIcgxASEiI3vXVajXvvfcefn5+TJ8+XZIKL+Hs9TDW+93SdXtQFIiJS+aPm/+YODIhCq+EhAQePnxo6jBELqhU1o7le65x6c5j1vndwtHOmhtBT6lUVpIKQmQlrVuZML29p/6tuyb1FYQwnZdKLPj5+dGyZUv69+/PpEmTuHPnDn/++SctWrTgxx9/zPbgmjVrRkxMDGfOnNFti4yM5MKFCzRt2jTT43x8fDh+/DhJSUm6bYcOHcLc3Bxvb2/dmHPnzhEVFaUbc+rUKeLj43XnbtasGXfu3CEgIEA35u7duwQEBOjGaDQaRo8ezeXLl1mwYAGDBw/OludeWOzyv0dSilZvW7JGyy7/eyaKSIiCp2fPnhk66SQnJ3P06FGePs1YkO/w4cOvVJdB5D9p9ROW7bnO09hEniWmMLqHlyx7EOIFLCxMNulXPCdKncjpq2G6x9IVQgjTMTqxcOrUKSZNmkSlSpX45JNPdLMIXFxc8PDwYP78+fzyyy/ZGlzDhg3x9vZm4sSJbN++nd9++40hQ4ZQvHhx3nnnHd24u3fvcuPGDd3j4cOH8/jxY0aMGMHx48dZu3Yts2bNok+fPrq1xv3798fKyoohQ4bw22+/sX37dj7++GN8fX117YI6d+6Mq6srw4cPZ//+/fz666+8//77VK1alU6dOgHw3//+lz/++INevXpRtmxZLl26pPtz+fLlbP15FERpMxWM3S6EeHk3b97MMLNMrVbzwQcfcOvWLRNFJfIMRSFFk5rgbVPfRZIKQvy/qKgoHj58qPcnLCz1Q2xkZGSGfTLTK/dJVwgh8g6j061Lly7Fy8uLDRs2EB0dzezZswFwd3dn8+bNDBo0iPXr1/Pmm29ma4BLlixh9uzZzJ07F61WS/369Vm0aBElSpTQjZkxYwahoaEcO3ZMF9OaNWuYO3cu48aNw97envfee48PP/xQd4yDgwMbNmzgu+++Y/LkydjY2NCxY0e9ZQxWVlasXbuWb7/9ls8++wwrKyt8fHyYOnWqLlN96NAhALZu3crWrVv1Yjc3N9dLeIiMHO2KGEwiONrJ2jghclr6ZWai8ElK1rD61xuoVNCpUUWOXwyluhRqFAKA7777ju+++87gvsmTJ2fYplKp5H1fLrsbGoNGK10hhMgLjE4s3Lx5kwkTJhjshmBhYUHXrl2ZN29etgYHUKJECWbNmsWsWbMyHWOoMFmDBg3Ytm1bluf28PBg3bp1WY4pW7YsS5YsyXT/hg0bsjxeZK2pVxn2ndGvY2FlYUavFu4mikgIIQqPdQdvEaVO4u1W7nRq5ErNyo4s33NNlkOIQq9nz56mDkEYYcZQb36/HsaqfTcY17s2dauWMnVIQhRaRicWLC0tSUlJyXR/VFRUlu0hhUgvOUXLhduPsS1qgZWFOZGxiTjaFaFXC3ea1Cxj6vCEEKJAe/BYzbkb4dSsZE+nRq7AvzUXgh7FSGJBFGpZfaGVm+Li4vj+++85fPgw8fHxNGjQgE8//ZRKlSpledzVq1eZO3cu165dw8bGhl69evHhhx/qvVcPCgpi9uzZXLhwAXNzczp27MjHH3+s1w3jyZMnzJo1i1OnTpGSkkKLFi2YNm0aTk5OujGPHz/mhx9+4PTp00RFRVG5cmXef/993bLhnPa/K48AOHcjTBILQpiQ0YkFb29vduzYYbC1zj///MPmzZupX79+tgYnCrYDvwfzKCKeCX3qUMvN0dThCCFEoaFVFNYfvIWNtSUjutfU2+cpSyGEyDMmTJjA1atXmTJlCjY2NixZsoRBgwaxf/9+XQv19IKDgxkyZAj16tVj0aJF3Lt3j4ULF6JWq/niiy8AiI6OZvDgwTg5OTFnzhwiIiKYN28eYWFhrFy5EkjtBDds2DDi4+P56quvSElJYf78+QwfPpydO3diYWFBUlISw4cPJzY2lnHjxlG6dGkOHTrERx99hEajoWvXrjn684lSJ/L3/SgA/vr7CdHqRGk1KYSJGJ1YmDBhAv369aN79+74+vqiUqk4evQoJ06cYPfu3SQlJTFu3LicjFUUAGevh7HL/56uroJ7OTtJKgghRC7zvxjKvdAYhnf1pHgxK1OHI4Qw4MKFC/j7+7N69Wpde8sGDRrQpk0btmzZwogRIwwet2rVKooXL86yZcuwsrKiRYsWWFtbM3PmTEaOHImzszObNm0iJiaGPXv2YG+fmkh0dnZmxIgRXL58mTp16rB//35u3brFgQMHcHdPXaLq6elJ165dOXz4MJ07d+bkyZPcunWL7du3U7t2bSC189rDhw9ZvXp1jicWni/eqCW1aOPADtVy9JpCCMOM7gpRpUoVNm3aROnSpfn5559RFIWNGzeyfv16KlasyLp16/D09MzJWEU+d/Z6GOv9bukVawz5R83Z62FZHCWEyA4BAQGcP39e9+fixYsA3L59W2/7+fPnuXdP2r0WRH6/B3Mz+ClPYxPZ4X+PGpXssbctgt/vwS8+WAiR606fPo2NjQ0+Pj66bQ4ODjRs2DBDC+H0x7Vq1Qorq3+Thh07dkSj0XDq1CndmIYNG+qSCpDaZt3GxgZ/f3/dmCpVquiSCoDucdoYGxsb+vbtS61atfRicHNzIyQk5DWe/YtFqRM5dTWMtDrE0mpSCNMyesbCnTt38PDw4OeffyYqKoqQkBC0Wi3ly5fXW2clRGZ2+d8jKUWrty05Rcsu/3tSU0GIHLZixQpWrFiRYfucOXMybFMUBZVKlRthiVxUqawdy/dco5xjMVI0Cj5eZVj+y3VG9/AydWhCCAMCAgJwdXXF3Nxcb3vFihXx8/MzeMyzZ8949OgRlStX1tvu4OCAra0tgYGBunN3795db4y5uTkuLi56Y9KfJ+36aWOaNGlCkyZN9PYnJyfj7+9P1apVX+LZvrysWk3KrAUhcp/RiYUhQ4bQs2dPJk+eTMmSJSlZsmQOhiUKIkNtJbPaLoTIHh988IGpQxB5gKerPe0aVmD3yQCqVyzJlqN3pfuDECaSkpLC/v37M91fqlQp1Gq1XiHFNDY2NqjVaoPHxcbGArzwuNjYWKPGVKlSxeCY4ODMZzp9//33BAUFsXTp0kzHGOLomDEeQ5ycUmtLBIWr0WgytpoMCo/VjRGmI69B/pCdr5PRiYX4+HhcXFyy7cKi8HG0K2IwieBoJ0V2hMhJklgQAMkpGk5deYhtUQtuhUTRrWklSSoI8QrCw8MJCwvDzc2NIkWKYGFhYbAde1YSExOZMmVKpvu9vb2z7LaW2fWU//8G39CsM0VR9I7LrjHPb583bx7r1q1j2LBhtG3bNtP4DYmIUKPVKlmOcXIqzuPHqcmTzwc1AOCLn84RGZvIt8Mb6Qo3po0RpvH86yTyLmNfJzMzlVGJP6PvgoMHD2bNmjVcuHDB2EOE0NOsVtkM26wszOjVwt3AaCGEENnJ7/cQHkcloNEqdGtaieMXQ7kZ/NTUYQmRb/z555/06tWLli1b0q9fP65du8Yff/xBy5YtOXDgwEudy8bGhtu3b2f65+eff8bW1pa4uLgMx8bFxRmcbQD/zlQwNKMhPj5e10nC1tbW4Jjnz23MmDRJSUlMmjSJn376iWHDhmWZNMluj6MTiE9IYe/poFy7phAiI6NnLFy7do3Hjx8zcOBArK2tKVmyZIZspUql4siRI9kepMj/tFqFS3cjKFrEHGsrC57GJuJoV4ReLdylvoIQQuSwx1HP2HcmCAtzFR/0qo2nqz3VXe1ZvueaLIcQwghXrlzhvffeo2zZsgwePJj169cDUKJECSwsLJg8eTI2Nja0aNEi265ZuXJlzp49m6HuTXBwsMHaB5CasHB2ds6wVCEiIgK1Wq07rnLlyhnGaDQaHjx4QIcOHXRj/v777wzXCAkJoU6dOrrHarWakSNH8tdffzF9+nQGDx78ak/4FUSpE0lM0gBw6uojuvtUknaTQpiI0TMWEhMT8fLyokGDBnh5eeHi4kK5cuX0/pQtm/EbaSEA/C8/JDg8loEdqjF/rA9rprZm3hgfSSoIIUQu2HLkDgDDu9bQJRE8Xe0Z3cOLoEcxpgxNiHzhhx9+wMXFhV9++YURI0bolhzUqlWLvXv34u7uzsqVK7P1ms2aNSMmJoYzZ87otkVGRnLhwgWaNm2a6XE+Pj4cP36cpKQk3bZDhw5hbm6Ot7e3bsy5c+eIiorSjTl16hTx8fG6czdr1ow7d+4QEBCgG3P37l0CAgJ0YzQaDaNHj+by5cssWLAgV5MKAL+c+je2tMKNQgjTMHrGws8//5yTcYgCLDY+iV3+96hWoSSNPJ1NHY4QQhQqV+494dLdJ7zd0h3vdPdgT1d7ma0ghBEuXrzImDFjsLa25tmzZ3r7bG1t6dOnD//5z3+y9ZoNGzbE29ubiRMn6oqnL168mOLFi/POO+/oxt29e5ekpCRq1KgBwPDhw9m/fz8jRoxg8ODBBAUFsWDBAvr06UO5cuUA6N+/Pxs3bmTIkCGMHTuWqKgo5s2bh6+vL2+88QYAnTt3ZsWKFQwfPpxJkyahKArz58+natWqdOrUCYD//ve//PHHH/Tt25eyZcty6dIlXVwqlUpvZkN2i1IncvpquO5xWrtJmbUghGkYnVgQ4lXt9A/gWaKGAe09pIWdEELkouQUDZt/u0NZx2K0a1jB1OEIka9ZWVllui8xMRGtVpvp/le1ZMkSZs+ezdy5c9FqtdSvX59FixZRokQJ3ZgZM2YQGhrKsWPHAHB3d2fNmjXMnTuXcePGYW9vz3vvvceHH36oO8bBwYENGzbw3Xff6ZZxdOzYUa82gpWVFWvXruXbb7/ls88+w8rKCh8fH6ZOnYqFRepHiEOHDgGwdetWtm7dqhe7ubk5N27cyPafSZp9pwN1M0fSSLtJIUzH6MRCmzZtjBp39OjRVw5GFBxnr4exy/+erguEV2UHyjsZ10ZICJE7UlJSuHr1Ko8ePcLb2xtra2s0Go3eG1aR//j9HkylsnZ4utpz8FwI/0Q94+1W7vx2/j6dGruaOjwh8qU6derw66+/MmjQoAz74uPj2b59O7Vq1cr265YoUYJZs2Yxa9asTMcYmlXcoEEDtm3bluW5PTw8WLduXZZjypYty5IlSzLdv2HDhiyPz0l3Q2PQaDO2m7wbGm2iiIQo3IyusZC+nkK5cuUoU6YMFhYWhIaGYmlpScuWLXMwVJFfnL0exnq/W3qtJW/fj+Ls9TATRiWEeJ6fnx8tW7akf//+TJo0iTt37vDnn3/SokULfvzxx5c+36+//kqXLl2oXbs2nTp1Ys+ePUYfO2fOHIYMGZJh+4ULF6hWrVqGPyNHjtQbt379etq1a0ft2rXp2bMn/v7+Lx1/QVKprB3L91zj9xth7D8bjEeFEvj9HkKlsnamDk2IfGvcuHHcuHGDAQMGsGfPHlQqFVeuXGHDhg28+eabPHjwgFGjRpk6zEJlxlBvPuyVmsz5YkgD1kxtzZqprZkx1NvEkQlROGVLjYVr164xfPhwXUEYUbjt8r9HUor+dMDkFC27/O9JsUYh8oBTp04xadIk3njjDYYPH87s2bMBcHFxwcPDg/nz5+Pk5MSbb75p1Pn8/PyYPHkygwYNonnz5hw5coRPPvkEa2trOnbsmOWxGzduZM2aNTRp0iTDvtu3b1OsWDHWrl2rt93O7t8PyD/++CMLFizggw8+oGbNmuzcuZMxY8awceNG6tWrZ1T8BU1aUcaF2y6hVRRCH8cxpmctqaUgxGuoV68eK1eu5Msvv2TOnDkALFy4EAAnJycWLlxI48aNTRliofTP09R6F0q6mQtCiNyXLTUWvLy8GDBgAEuXLtW1qBGF1/MzFYzZLoTIXUuXLsXLy4sNGzYQHR2tSyy4u7uzefNmBg0axPr1641OLCxYsIBOnToxffp0AJo3b050dDQ//PBDpomF8PBw5s6dy4EDB3R91dO7desWVatWpW7dugb3x8fHs2LFCoYOHcqYMWMA8PX1pV+/fixduvSVZl4UFGYqSNGkvtFu/YaLJBWEyAY+Pj789ttv3Lhxg5CQELRaLeXLl8fLy0tXc0Dkrgu3/wHgxKWHVC4ny/iEMCWjl0K8SKlSpQgKCsqu04l8zNHOcCXezLYLIXLXzZs36dKlC2ZmGf8JsLCwoGvXrgQGBhp1rvv37xMSEkL79u31tnfo0IGAgADu379v8LiFCxdy48YN1q5di6enZ6ZxVquWeQGuy5cvExsbq3dtlUpFu3btOHv2rF6rtcJEqyis87uFSgWdG1fk+MVQbgY/NXVYQhQIKpWKmjVr0qlTJ7p06YKNjQ3BwcGmDqtQilInEvgoFoDfb4QTrZYvsIQwpWxJLDx+/JgtW7boWtiIwq19w4oZtllZmNGrhbsJohFCpGdpaUlKSkqm+6OiorC0tDTqXGn9zStXrqy33dU1tUhgZgmKtHZomU0d1mq13Llzh7CwMHr27ImXlxctW7ZkzZo1uirgadd2c3PLcO2UlJRMkxoF3c4T9wh/+ozOjVx5q2UVRvfwYvmea5JcEOI1KIrCqlWrmDZtGpB6jxoxYgTdu3ena9euDB06lLi4OBNHWbg83xUirRuEEMJ0XrsrRFJSEpGRkWg0Gr788stsC0zkX4GPYjBTgZ1NEaLUiTjaFaFXC3epryBEHuHt7c2OHTsYMGBAhn3//PMPmzdvpn79+kadKzY29dsiW1v9ri82NjYAqNVqg8dVqVIly/MGBgaSkJBAYGAgEydOxN7enqNHjzJ37lzUajXjxo3TnTvtWumvXRjf5CcmaTh+MZQyDkXp2SI14ZJWcyHoUYwsiRDiFf30008sWLCA5s2bA6m1ZU6ePEmHDh2oWrUqq1evZunSpXrtGkXOiVIncupqGGmVFTQahVNXH9HdpxIlbGWGrBCmYHRiIbPZCObm5jRq1IiuXbtKVwjBvYfR/H4jnK5NXenlKzMUhMiLJk6cSN++fenevTu+vr6oVCqOHj3KiRMn2L17N0lJSYwbN86oc6V9W6RSqQxuN7TcwhjOzs6sXr0aT09PnJycAGjSpAkJCQmsXr2aoUOHoihKhutmFdOLODpmbInr5GS4/kNetfnQLRKSNEzoXx/n0v8WuXRyKo6vCePKTfntNSvs8svrtXv3btq1a8fixYsBOHDgAEWLFmXOnDlYW1sTFxfHwYMHJbGQS/adDkSr6BdsTJu1MLBD5kvohBA5J1u6QggBqW/m/3v0DiVsrOgsvdKFyLPc3d3ZtGkTM2fO1N3bN27cCKQW4/3ss88yrXuQXlrhxfQzE9JmC2RWmPFFbG1t8fXN+FG4ZcuWbN++ncDAQIoXL46iKMTFxenNmHjVa0dEqNE+V1ncyak4jx/HvlL8phAZk8DOY3doUL00TrZW+Sr27JLfXrPCztjXy8xMZTDxl5vu37+va4ubnJzM2bNn8fb2xtraGki9rz558sSEERYud0Nj0Gj0EwsajcLd0GgTRSSEeK0StsnJyZw+fRpzc3OaNGkiFXELqbPXw9jlf0/X9aF57bJYW8nfBSHysmrVqvHzzz8TFRWlV908bXaAsdJqK4SEhOgVWkwrZpa+9oKxbt++zZ9//snbb7+tV+8hISEBAHt7e71r16hRQ+/aVlZWhaLuj9/vwVQqa4enqz27TgagVRTeqFoKv9+D6SQJXiGyjZ2dnS6Beu7cOeLj4/WSnyEhIZQqVcpU4RU6M4amtrifvup3KpS2ZXQPLxNHJIQweo5qUlISX3zxBUOHDtU97tOnD6NHj2bEiBH06NGDiIiIHAtU5E1nr4ex3u+WXivJczfCOXs9zIRRCSFe5Nq1a0yYMAGNRkPt2rWpW7cua9asYdy4cdy7d8/o87i6uuLi4sLBgwf1th8+fJhKlSq98of74OBgZsyYwcmTJ/W2HzhwABcXF8qXL0+9evUoVqwYhw4d0u1XFIXffvuNhg0bYmVl9UrXzk8qlbVj+Z5rHP/rAWeuhfGGhxObj9yhUlm7Fx8shDBavXr12LhxI4cPH2bBggVYWFjQvn17kpOTOXz4MFu2bKFRo0amDrPQUT9L5u/7UdIRQog8wOjEwpIlS9i2bRtly5YFYM+ePdy8eZOBAwfy3Xff8fjxY3744YccC1TkTbv875GUotXblpSiZZe/8R9MhBC568KFC/Tv35/Tp0/z9Om/nQKcnJz4888/eeutt7h165bR5xs7diy//vorX3/9NSdPnuSrr77Cz8+P8ePHAxAZGcmlS5cyLeRoSMuWLfHy8uLzzz9n+/bt/O9//2PKlCkcO3aMqVOnolKpKFq0KEOHDmXVqlX85z//wd/fn/Hjx3P9+nXGjh1r/A8kH/N0tWf0mzXZfOQOVhZm3Ah6yugeXlKkUYhsNn36dIoUKcK4ceO4efMmkyZNwsnJib/++otx48bh5OSku+eJ3BP3LJnouCTpCCFEHmD0fHU/Pz/eeustZs6cCcChQ4coXrw4U6ZMwcLCgvv377N9+/YcC1TkTc/PVDBmuxDC9H744QcqV67M+vXrKVmypG770KFD6d27NwMHDmT+/PmsXr3aqPP16tWLpKQk1qxZw/bt26lQoQJz5syhc+fOAJw4cYJp06axYcMGo7/Rs7KyYvXq1SxatIglS5YQGRlJ1apVWbJkCW3bttWN++CDDzA3N2fbtm38+OOPVKlShWXLlhnd1aIgSEjSoNEqaLQKHbwrSlJBiBxQtmxZ9u7dy40bN3B2dsbZ2RmA6tWrs2DBAlq1akXRokVNHGXh8iQqXtcVQjpCCGF6RicWwsLCqFu3LgDPnj3j/PnztGzZUldXoWzZssTExORIkCLvcrQrYjCJ4GgnN3Yh8qqbN28yceJEvaRCmhIlStCnTx9d5XNj9evXj379+hnc16tXL3r16pXpsZkVB3ZwcODrr7/O8roqlYoxY8YwZswY44MtQLRahc1H7qBSQZfGrhy/GEp1V3tJLgiRAywsLKhdu7bethIlSuiSqCJ3/XIqSPf/0hFCCNMzOrFQqlQpXbXb//3vfyQlJem1l7x9+zalS5fO9gBF3ta+YUW2HL2jt83KwoxeLaTVpBB5lYWFhd4SiPTUajVarTbT/SLv2OF/l4iYBN70qcSbzd3wrOTA8j3XZDmEEK9p2rRpL32MSqXiu+++y4FoRHpR6kTO3QzXPdZoFJm1IISJGZ1YaNSoEevXr6dIkSJs2rSJokWL0rZtW2JiYti5cyfbtm3L9NsqUXCF/BOLSgUlbIoQpU7E0a4IvVq406RmGVOHJoTIRKNGjdi4cSPdu3enQoUKevvCw8PZuHEj3t7eJopOGCs5RYP/pYeUcShK92apHTI8Xe0Z3cOLoEcxklgQ4jXs3r37pY+RxELu2Xc6EEW/26TMWhDCxIxOLEyfPp3w8HDmzJlDsWLF+Oabb7Czs+PPP/9kzpw5NGzYkA8++CAnYxV5TOhjNWeuhdGuQQX6talq6nCEEEYaP348b7/9Nt27d8fX15dKlSqhUqkICQnB398flUrFxIkTTR2meIFjf4XyLFHDBz2roVKpdNs9ZSmEEK/tZQrYitx3NzQGjVY/s6DRKNwNjTZRREIIoxMLdnZ2rF27lsjISGxtbXVtvDw9Pdm6dSt16tTJsSBF3rTTPwBrK3O6Nq1k6lCEEC/Bzc2NXbt2sXDhQk6ePKlr12htbY2Pjw8TJ07E3V2WM+Vl8Qkp/HomiJqVHfCs5GDqcIQQIlfNGOrNuRvhrNx7nZnDG1GulI2pQxKi0DM6sZCmRIkSXLt2jdDQUKysrChTpowkFQqhOw+iuHT3CT193bAtamnqcIQQL8nV1ZVFixahKApPnz5Fq9Vib2+Pubm5qUMTRjj4RzBxCSm8JfVshMgVxtZcmDVrVg5HItLEJ6YAUMz6pT/OCCFywEv9Jh4/fpwZM2YQHh6Ooii6qZelS5fmyy+/pHXr1jkSpMg7zl4PY5f/PSJiElGpoKSNlalDEkK8BpVKhYODfOOd1/n9HkylsnZ4utoTpU7k8Pn7VK9YkhtBkbiWKW7q8IQo8F5Uc8HR0VHupbksIvoZAEkpUmxYiLzA6MTChQsX+PDDD3F0dGTChAm4u7ujKAoBAQFs3ryZcePGsWHDBt54442cjFeY0NnrYaz3u6W7gSsKbPrtbywszKRYoxB5WJs2bZg+fTpt2rTRPX4RlUrFkSNHcjo0YaRKZe103R4u3PqHlBQt9/9R082nsqlDE6JQMFRzQavV8vjxYw4cOMDKlSv5/vvvTRBZ4XUtIBKAQ+eCGdihuomjEUIYnVhYvHgx5cuXZ8eOHRQvrv/tSP/+/enduzfLly9n9erV2R6kyBt2+d/LkBVOStGyy/+eJBaEyMPKlStHsWLF9B6L/CWt28Oy3VeJT0jBwsKMMT1rSZFGIUzIzMwMZ2dn/q+9Ow+Lqmz/AP4dlgFlEQQVRATFDcWFXBLBfUWzUF+zzBTLLFLrLcnU/JVLuWCKZWa9pWalZZqSGyolUuBSuAKKxr4ILuzrDDNzfn8Qk8Mi28AZ4Pu5Lq+Xec45M/fwNPcLN8+5n/nz5yMpKQkbN27Enj17xA6rRcjOlyHlQT4AIDQiHU+7d+E2k0Qiq3Fh4caNG1i0aFGFogIAmJqa4j//+Q+LCs1cRq6sVuNEpBu+++47jceffvopLC35C2lT4+xgCQszIxQUKzDatSOLCkQ6pHfv3ggICBA7jBbj0e0muc0kkW7Q09YTSSQSlJSUaOvpSAe1Na+8EmxVxTgR6aZp06bh888/FzsMqqXQG3eR+qAATnbmOB95D7cSs8QOiYj+ERISAhMT7kzQGLLzZQiNSEfZZpNKpYDQiDTk5PMPXURiqvGKhf79++PQoUOYPXu2xpJaAMjPz8fBgwfRt29frQdIuqN/VysEX7urMSY10MN0diUnalIyMzNhbW0tdhhUC7cSs/Dt6dsw1NfDkhn9kPqgQN1zgSsXiBpeVbtCyOVy3L59G7GxsZg7d24jR9UyHQuLh6psucI/uGqBSHw1LiwsXrwYc+fOxVNPPYU5c+bA0dERANTNG+/du4c1a9Y0VJwkMqVKhZuJWbA0M4JEAmTmymBlboTpI53YX4GoiZk6dSoOHDiAYcOGoVOnTmKHQzVwPeYhFEoBk4d2hnlrKcwdpPDxckFCWi4LC0SNoKpdIfT09GBtbQ1vb2/897//bdygWqiY1FwolZqFBaVSQExqjkgRERFQi8LCoEGDsH37dqxduxZ+fn7qrSYFQUC7du3g7++PoUOHNligJK6LUfdwL6sIi6f3xRM92okdDhHVg56eHuLi4jBx4kR07twZVlZW0NPTvDNOIpFg7969IkVI5WXmyWAk1cekJzurx5wdLFlUIGokle0KQeJY89IQAMBb20PRv5s1vD25IwSRLqhxYQEo3aJs1KhRiIqKQkpKCgDAzs4Offr0gYFBrZ6KmhClSoVjYQno3MEUrt25fJqoqQsLC1M3b5TJZLh79241V5CYUu7nIzz6Pp4a5gDTVoZih0NEALKysqCvrw9zc3OxQ2mxCmUKtDbm7x9EuqLWn0Z9fX106tQJenp66q9ZVGjezkem4352Ed6Y0U+9UoWImq6zZ8+KHQLVwi+h8WhlpI8JgztXfzIRNZiEhARs374dISEhKCgoAABYWFhgwoQJWLJkCXvXNKIShRIlChX+ir6PiYPtudUkkQ6oVUUgPDwcH3/8MW7cuAHhn6Yp+vr6GDp0KJYtW4YePXo0SJAkHoWydLWCo40Z+nezEjscIqqHkpISxMTEQKFQoFu3bmjVqpXYIVE1ku7l4fKdB3ja3ZGrFYhEFB0djTlz5qC4uBjDhw+Ho6MjlEolEhMTcejQIQQHB+PAgQOwtbUVO9QWoVCmBABk5BSzaSORjqhxYeHSpUt4+eWX0bp1a8yePVudUBMSEnDs2DE8//zz+OGHH1hcaCYuRKXjcEgsMnJLt+4Z1LMdVysQNWHffPMNduzYgfz8fACAVCrF7NmzsXTpUq460zGBFxPhaGsOZwfLf1YrGMDRxgyBFxPhOdRB7PCIWqTNmzfD2NgYP/30E7p27apx7M6dO5g7dy4+/vhjbNmyRaQIW5b0zAL116ERaXja3ZGrFohEVuOfJrdt2wY7Ozv88MMPaNu2rcaxRYsW4dlnn8XWrVvxxRdfaD1IalwXotKxNzAacoVKPXb2SirsO5hxBwiiJiggIAAbN26EnZ0dnnnmGejp6eHSpUv45ptvoFQqsXLlSrFDpEc42ppjZ0Akpo/oiqt/P4S7iw12n4yGj5eL2KERtVhXr17Fa6+9VqGoAAA9evTA3Llz2fC2Ef0anqL+mltNEukGvepPKRUdHY3nn3++QlEBAKytrTF79mz89ddfWg2OxHE4JFajqAAAcoUKh0NiRYqIiOpj//79GDBgAE6ePIlVq1Zh5cqVCAgIwMSJE3HgwAHI5XKxQ6RHODtYwsfLBfuC7sBQX4LrsRnw8XLhDhBEIjIxMYFCoajyuFQqrbC7DjWM7HwZrsc8VD9WKgWERqQhJ18mYlREVOMMaGVlhYyMjCqPy2QymJqaaiUoElfZ7Q81HSci3RYbG4upU6fCyOjfZaISiQTe3t6Qy+WIi4sTMTqqjImxAZQqASVKAaNd7VhUIBLZnDlzsHfv3kq3nbx79y727duHF198UYTIWp5jYfFQ/dPrrUzZqgUiEk+Nb4V47bXXsG7dOri6umLMmDEax65fv45vv/0WS5cu1XqA1PiszI0qLSJYmfPeNaKmqKioCGZmZhXGO3XqBEEQkJubK0JU9Dg//Po3AGDSk/YIvpqKXg6WLC4QiUilUsHExAQzZsyAh4cHunbtCqlUiqSkJJw9exYSiQTJyclYsWKF+hqJRIL169eLGHXzFJOaC5XmwloolQJiUnPECYiIANSisHDt2jVYWVlh0aJF6Nq1K5ycnGBoaIjk5GRERERAKpXi+PHjOH78uPoaiUTC+82aoGkjuuLr47c0xqQGepg+0kmkiIioPlQqVaXNV/X19QEASqWysUOixwiLTMPt5Gw86dwBz47ujr5drbEzIJK3QxCJ6JNPPlF/HRISgpCQkArnHDlyROMxCwsNY81LQ/Db5RTsC7oD/yUeaGMiFTskIkItCgvnz58HANja2qKoqAiRkZHqY2Vb66SkpFR6LTUtRoalv2yYtjJEflEJrMyNMH2kExs3EhE1gt/CU6CvJ8Hz47oD+LfnQkJaLgsLRCKp7BYIEk+xvLTfRSupvsiREFGZGhcWzp4925BxkI4QBAHHLySivWUrrH9lKPT0uMUkUXOQnZ2Nu3fvaozl5JQuG83MzKxwDAA6duzYKLHRvzJyipF8Px+jXO1g/shf4Zx5KwQRkVr2P40aC4tLIDVkcYFIF9S4sBAaGgoPD48qj8vlcnzyySd45513tBIYiSMqPhOJ6Xnw9uzFogJRM7J+/foql+T6+vpWGJNIJLh582ZDh0XlnLqUBACYNKSzyJEQUXkBAQEICwvDgwcPoCp/kz94C3Bjik7MBgAcO5/IbSaJdESNCwsLFizA7NmzsWzZMhgbG2scu3LlClauXInExEQWFpq44xcSYWlmhGEuvO2BqLmYNm2a2CFQDeQUyPH7jbtwc7GBVRvj6i8gokbj7++PL7/8EoaGhrCysuLWkiLKzpchLaMAABAakYan3R3RxpQNxonEVuPCwsyZM7F//36EhYXBz88P/fv3R3FxMbZs2YJ9+/bB3NwcH374YUPGSg3sTnI27iRn4/mx3WGgz//DJGouNmzYIHYIVANn/kqCQqnC5KEOYodCROUcOXIEHh4e2L59O1q1aiV2OC3asbB4lO02WbbNJFctEImvxoWFdevWYdKkSVi1ahVmz56NWbNm4Y8//kBKSgpmzJgBX19fWFhYNGCo1FAuRKXjcEgsMnJlkEgAYyPeq0ZE1NACLybC0dYczg6WKCguQfCVVPS0t8DVOw/gyeICkU7Jz8/HxIkTWVQQWWZuMUIj0vFPXQFKpcBVC0Q6olZ/lnZ3d8eRI0fQsWNH7N+/HykpKVi2bBk+/PBDFhWaqAtR6dgbGI2M3NImOIIA7DtzBxei0kWOjIioeXO0NcfOgEjcSszCb5dTUCxXIulePhxtzcUOjYjKGT58OC5evCh2GC3ej0G3oSpbrvCPslULRCSuGq9YAICwsDCsXbsWycnJmDBhAm7evInNmzcjKSkJS5cuhampaUPFSQ3kcEgs5ArNBkRyhQqHQ2K5vSQRUQMq20ZyZ0AkZHIF9PUkWDS9L3d/INJB//d//4f58+dj6dKlGDduHKysrCCRVGxyPXjwYBGiazmiEzKhVGoWFpRKATGpOSJFRERlalxYeOutt3Dq1ClYW1vj888/x5gxY1BUVIQtW7Zg//79CAoKwooVKzBlypSGjJe0rGylQk3HiYhIe5wdLOHQwQxRCZlw72vLogKRjrp79y7y8vJw4sQJnDx5ssJxQRAgkUhw69YtEaJrOT5dOhoPHuRh6Y4w9HFsi5emOIsdEhH9o8aFhcDAQEybNg0rVqyAuXnpMs1WrVph1apV8PT0xHvvvQdfX18WFpoYK3OjSosIVua8T42IqKFFxWfgZmImLM2McD3mIW4lZrG4QKSD1q5di9zcXLz88stwdHSEgUGtFv2SlhXLFewJRqRjapwVv/rqKwwfPrzSYwMHDsTRo0exfft2rQVGjWOymwO+O31HY0xqoIfpI51EioiIqGW4lZiFHUciIQjAC+N7oJWRAXYGRMLHy4XFBSId8/fff2Px4sV45ZVXxA6lxRMEAcUyJYylLO4Q6ZIaN2+sqqhQRiqVYunSpfUOqLyCggKsWbMG7u7ucHV1xSuvvIKEhIRqr4uIiMCLL74IV1dXeHh4YOvWrSgpKdE4JyEhAa+99hoGDRqEJ598Eh988AHy8/M1znn48CGWLl2KJ598EgMHDsTbb7+NBw8eaJxz//59+Pr6ws3NDU888QRef/11JCYm1vu9N4a8wtLvSRsTKYDSlQrzPHuxvwIRUQOLT8uBhakU7S1bYUA3a3XPhYS0XLFDI6JybGxsoKfHrbh1gaxECQHAxah05OTz1l0iXVFlhly8eDHCw8M1xlQqFaKjo1FUVFTh/KNHj6J3795aD7Cst4Ovry82bdqEe/fuYe7cucjLy6vymsTERHh7e8PIyAjbtm3DSy+9hD179mjs5Z6Tk4N58+bh4cOH2LRpE5YuXYqTJ09qFEcUCgVefvll3LhxA6tXr8bq1atx5coVLFiwAAqFAgAgk8mwYMECRERE4P3338eWLVtw//59zJkzB7m5uv3DYYlCid8up6CfkxX8l3hg9/Ix2Py6O4sKRESNwKljG6RnFmHiYHvo6ZU2gXN2sORWk0Q6aMGCBdi7dy9iYmLEDqXFK5YrAQAPc4q5GwSRDqlyDdGvv/6KiRMnaozl5ORg2rRp2L17N9zc3CpcI5Tb/qW+wsPDERISgq+++gojRowAAAwaNAhjx47FDz/8gIULF1Z63f/+9z+YmZnh888/h1QqxciRI2FsbIwPP/wQr776Kjp06IB9+/YhNzcXAQEBsLQsXXLaoUMHLFy4ENevX0f//v1x4sQJREdH4+TJk3ByKr01wNnZGU899RTOnDmDyZMnIzg4GLdv38bPP/8MFxcXAED37t0xduxYnD59GjNnztTq90SbzkemI6+wBBMH24sdChFRi3P6z2SYtjLEsL62YodCRNWIjo6Gnp4enn76adjb28Pa2hr6+pr3+EskEuzdu1ekCFuOe1mF6q9DI9LwtLsj2piyNxiR2Gq9pkvbxYPHCQsLg4mJCdzd3dVjbdu2xeDBg/H7778/9rrRo0dDKpWqxyZNmgSlUonQ0FD1OYMHD1YXFQDAw8MDJiYmCAkJUZ/TrVs3dVEBgPpx2TkeHh7Yv3+/uqgAAIaGhgAAuVxen7ffoFSCgDN/JaNzB1P04r28RESNKi2jANdiHmLME3YwMmQDMiJdFxwcDD09PdjY2KCkpARpaWlISUnR+JecnCx2mC3Cr+Ep6q9VgsBVC0Q6Qqe7nsTFxcHBwaFCRbhz584IDAys9JqioiKkpaWhS5cuGuNt27aFqakp4uPj1c/99NNPa5yjr6+PTp06aZxT/nnKXr/sHFNTUwwcOBAAUFJSgtjYWGzatAkWFhYYP358Hd5147gRm4G0jEIsnNq70n2YiYio4Zz5KxkG+noY80QnsUMhoho4e/as2CEQgMzcYlyPeah+rFQKXLVApCNEKywoFAqcOHGiyuPW1tbIz8+HqalphWMmJiYVmiyWKeu9UN11eXl5NTqnW7dulZ5TWXPGJUuWqCvaH330Edq3b1/l+6uMlVXFeACgXTuzWj1PTZw9eB3WFq3gOdwJBvpsRqRNDTFf1LA4Z9SYcgvkCItIxzAXG5ibSKu/gIharIKCAnz88cc4c+YMCgsLMWjQILz33ntwdHR87HURERHw8/NDZGQkTExMMH36dCxZskS9qhYobWK+ceNGhIeHQ19fH5MmTcI777yj8fPxw4cPsWHDBoSGhkKhUGDkyJFYsWIF2rVrpz7n/v378PPzQ1hYGGQyGYYOHYp3330XDg7a7xfzY9BtqMqtni5btfDixJ5afz0iqjnRCgsymQzLli2r8viQIUM0kl95VXXmLbtVo7K/wguCoHGdts4p88orr2DevHk4evQoVqxYAQCYPn16le+hvIyMfKhUmsmyXTszPHhQdaPK2rgQlY7DIbHIyC3toPukc3tkZRZo5bmplDbnixpHTeZMT09SZeGPqCYCLybC0dYczg6WOHslBQqlCt07mSPwYiKbNRLpoICAgDpd5+XlpdU43nrrLURERGDZsmUwMTHBZ599hrlz5+LEiRMwM6u8KF7WxNzV1RXbtm1DbGws/P39kZ+fj/fffx/Av03M27Vrh02bNiEjIwObN29Geno6vvzySwD/NjEvLCzE6tWroVAosGXLFixYsAA///wzDAwM1E3MZTIZ3n//fRgbG2PHjh2YM2cOTpw4AXNzc61+P6ITMqFSaY4plQJiUnO0+jpEVHuiFRZMTExw+/btx57zxhtvICUlpcJ4QUFBpasNgH9XKlS2oqGwsFCdhE1NTSs9p6CgAHZ2dtWeU9nrl90S4ebmhtTUVHz55Ze1Kiw0pAtR6dgbGA254t9sfOXvh7gQlc5dIIiIGpijrTl2BkRiwVPOOHslFU4dzXHgbCx8vFyqv5iIGt3y5ctrdauoIAiQSCRaLSywiXlFny4djQOnb+H7M3fgv8RDvV06EYnvsYWF7Oxs3L17V/04J6e0GpiZmakxDgBZWVlaD65Lly64cOGCOlmXSUxMrLT3AVBasOjQoUOFWxUyMjKQn5+vvq5Lly4VzlEqlUhJSVHvhtGlSxfcuXOnwmskJSWhf//+AICbN28iPj4eU6ZM0TinT58+iIiIqOU7bjiHQ2I1igoAUKJQ4XBILAsLREQNzNnBEj5eLvj00A3ISpRQqgQsnt4XzmyeS6STHt2iXCzVNTGvqrBQVRPzNWvWIDQ0FDNmzKi2iXn//v2rbWI+efJkUZqYF8lKt3xvJWXjWyJd8tjCwvr167F+/foK476+vg0W0KM8PDzwxRdf4Pz58+qkmpmZifDwcLz66qtVXufu7o7g4GAsW7ZMnVRPnz4NfX19DBkyRH3O7t27kZ2dDQsLCwBAaGgoCgsLMWzYMPXrnzhxAnFxcejatSsAICYmBnFxcfDx8QEAXLx4EX5+fujbty86d+4MoLRAcfHiRfTo0UP735Q6Krv9oabjRESkXT07W8BAXwJZCTBuoB2LCkQ6bNq0aWKHwCbmVSiWK6EnkcDQgD3CiHRJlYUFXUiogwcPxpAhQ/D222/D19cXFhYW2L59O8zMzPD888+rz4uJiYFcLkfv3r0BAAsWLMCJEyewcOFCzJs3DwkJCdi6dSueffZZdOzYEQAwe/ZsfP/99/D29saiRYuQnZ2NzZs3Y8SIEXjiiScAAJMnT8YXX3yBBQsWYOnSpRAEAVu2bEH37t3h6ekJoLSHwnfffQcfHx8sWbIExsbG2LdvH+7cuYPdu3c38nesalbmRpUWEazM2UGXiKgxnLqUhIJiBVy7WyP46l30cmjL4gJRC8Um5hXVtJeRRF8PrY0N0L69dvs3kHaxMXbToM15qrKwoAtLwADgs88+w8aNG+Hn5weVSoWBAwdi27ZtaNOmjfqcNWvWIDU1Vb0VkJOTE3bv3g0/Pz+88cYbsLS0xPz587FkyRL1NW3btsW3336L9evXw9fXFyYmJpg0aZJGQ0mpVIo9e/bgo48+wqpVqyCVSuHu7o7ly5fDwKD0W2dhYYHvv/8eH3/8MdauXYuCggL069cPe/fuxaBBgxrpu1S96SOdsPvELSgfaQ4pNdDD9JFOj7mKiIi04VZiFo78HofWxgZ47RkXxKTmYGdAJHy8XFhcIGqB2MS8osqamJfXrp0Z7j3IR7FcgZj4h9xiUkexmXnTUNN5qmkTc9GaN9ZUmzZtsGHDhscWOr777rsKY4MGDcJPP/302Ofu0aMHvvnmm8eeY2tri88+++yx59jZ2cHf3/+x54htaO8OOBgcg7zCEihVAqzMjTB9pBP7KxARNYKIuAwoVQLGDewEQwM9dc+FhLRcFhaIWiA2Ma+7+PQ8KJTcYpJI1+h8YYG0IzoxC9n5csyf3AvD+3UUOxwiohZFqRSgryfBKFc79ZizgyWLCkRUJTYxrygztxgZucUAgNCINDzt7shVC0Q6gl1PWohfL6fAtJUhhvbuIHYoREQtSrFcgdCIuxjUqz0s+AMwEdWQh4cHcnNzcf78efVYWRPzskbjlSlrYv7orgyVNTG/dOkSsrOz1edU1sT877//RlxcnPqcsibmZedcvHgRS5cuRVJSkvqchmxi/mPQbfxztwdUQumqBSLSDSwstAAPs4twLeYhRg7oCEMDbs1DRNSYzkemo0imxLiBncQOhYiakEebmB88eBBBQUHw9vautIn5zZs31Y8XLFiABw8eYOHChQgODsaePXuwYcOGCk3MpVIpvL29ERQUhIMHD+Kdd96p0MTcwcFB3RT9+PHjeOWVVyo0Mbe1tYWPjw9OnTqFc+fO4bXXXsOdO3fw9ttva/X7kZ0vw29/PlrAEBAakYacfO5wRqQLWFhoAc5eTYUEEox+ZAkuERE1PEEQ8NvlFDjamKFrR3YwJ6La+eyzzzBmzBj4+flh+fLlsLGxwTfffFOhifnixYvVj8uamBcWFuKNN97Anj17MH/+fLz33nvqc8qamFtYWMDX1xf+/v6YNGmSRs+wsibmvXv3xqpVq7Bu3Tq4urpi165dFZqY9+jRA2vXrsWbb76J4uJi7N27F08++aRWvxfHwuKhEjSbO3LVApHukAiC8Pj2q9RoKuuGW9+uqrISJXx3hMHZwRKvT+tb3xCpGuyC2/TUZM5q2g23JTt+/Dh27tyJ5ORk2NnZ4dVXX4WXl1eNrt20aRNu3br12Ga6+fn5mDp1KoYNG4aPPvpIPZ6dnQ0/Pz+cPXsWMpkM/fr1w1tvvYUBAwbU+j2Uz8Ha+DxHxWdiy4FreHmKM9z72tbruah6zMFNi7Y7klPTVt2uEB/s/hPJ9ys2k7Rvb4o1Lw1pyNColpiLm4YWtysE1c/FqHQUFCswlktwiaiBBAYGwtfXF3PnzsXw4cPx66+/4t1334WxsTEmTZr02Gu///577N69G25ubo89b8OGDbh7967GWElJCby9vZGcnIzFixejR48eCA4Oxty5c/Htt9/Wqbigbb+GJ8O8tSGGOLO/DRFRfax5aQjatTPDf5YfxyjXjpg1prvYIRHRI1hYaKYuRKXjcEgsMnJl0NeTqDvoEhFp29atW+Hp6YmVK1cCAIYPH46cnBx88sknVRYW7t27Bz8/P5w8eVK9/VlVQkJCEBgYWOG84OBg3Lp1C9u2bVPf7+vu7o6HDx9i48aN+PHHH7Xw7mov8GIiHG3NYWVuhBuxGZgyzBExqTlISMuF51AHUWIiImoOVCoBshIljAzZM4xI17DHQjN0ISodewOjkZFb2sxGqRLw7anbuBCVLnJkRNTcJCcnIykpCRMmTNAYnzhxIuLi4pCcnFzpdf7+/rh58yb27NkDZ2fnKp8/JycHq1atwjvvvANzc80eBfHx8QCAUaNGaYwPHjwYV69eRU5OTh3eUf052ppjZ0AkDp2LhZ6eBPbtTLAzIBKOtuyxQERUH/ISJQDASMrCApGuYWGhGTocEgu5QqUxJleocDgkVqSIiKi5KtuGrPye6g4OpX+ZL/vlv7yyLuNDhw597POvW7cOTk5OeO655yocs7Ut7VlQ/haJsmJGSkpKDd6B9jk7WOLlKc4Iv/0AHdq2wndn7sDHywXODpaixENE1FwUyRUAwBULRDqIt0I0Q2UrFWo6TkRUV3l5pU1/TE01m/qYmJgAKG26WJlu3bpV+9xBQUH47bffcOzYMUgkkgrHR4wYASsrKyxbtgxr166Fvb09/vjjDxw+fBgAUFRUVKv3UlljonbtHn+bRlWUsRkAgLsPCzFrfA+MGNS5Ts9DtVfXOSNxcL6oNu5lFgIo3WqSiHQLCwvNkJW5UaVFBCtzIxGiIaLmrGxjofK/+JeN6+nVbWFcZmYmPvjgAyxbtgydOlXefNbCwgI7d+7Eu+++i+nTpwMAevXqhTfffBNr166FsbFxrV5Tm7tCHPz1DvQkwBQ3B5wIjUdnaxOuWGgE7ETetHBXCKqtX/5ZfXsjNgPjB9uLHA0RPYq3QjRD00Z0rTAmNdDD9JFOIkRDRM1ZWUPF8isTCgoKNI7X1urVq+Hk5IT//Oc/UCgUUChKl78KgqD+GgD69++PU6dO4dy5cwgKCsIvv/yiXi3x6D7vjSn4SgruZRVh9BOdMG2EE3y8XLAzIBK3ErNEiYeIqDnIzpfhYmRpv7DbyVnIyedKXCJdwsJCM9Ta2BAAYNqq9H+tzI0wz7MX3PrYiBkWETVDZb0VkpKSNMYTExM1jtfW6dOn8eeff8LFxQV9+vRBnz59kJqaip9//hl9+vRBSkoKsrKycPjwYWRnZ8PW1hadO5febhAVFYU2bdpUudKhoYVGpMFAX4Jpw0uLvM4OlvDxckFCWq4o8RARNQfHwuLVq+EEATgaliBuQESkgbdCNEPnrqaijYkUm18fBgN91o6IqOE4ODigU6dOOHXqFMaPH68eP3PmDBwdHdGxY8c6Pe+hQ4cqjPn4+KBfv37w8fFB+/btkZeXhxUrVmD9+vWYMWMGgNJbKI4dO4axY8dW2pehoRUWK5D6sABufWzQ2vjf/4t1drDkrRBERHWUnS9DaEQ6lP/crqZUCQiNSMPT7o5oY8pbfYl0AQsLzczDnCJE/LNvOosKRNQYFi1ahBUrVqBNmzYYNWoUzp49i8DAQPj7+wMo/WU/KSkJ3bp1q9DksSp9+/atMCaVSmFpaak+ZmVlBU9PT2zduhVSqRSmpqb49NNP1TGJ4eLNdMhLVBjlaifK6xMRNUfHwuKhEjQbNqoEAUfDEvDixJ4iRUVEj+Jvns3M79fTAAAj+tuKHAkRtRTTp0/HmjVrEBoaikWLFuHPP//Epk2bMHnyZADAuXPnMGvWLERFRWn9tT/88EOMHTsWGzZswLJly2BjY4N9+/aJchuEIAg4dzUVDh3M4GjDTvdERNoSk5pbYScIpVJATGqOSBERUXkSQRC4X4uOKN+RHKhdh2uFUoV3dp6HQwcz/Hdm/4YIkarBjuRNT03mjB3JW4b67goRk5qD9d9dxtxJPTFqAFcsiIE5uGnhrhD0qMp+Di4v7OY97Doahc/+O0LjdjPSLczFTYO2czBXLDQj12MeIidfzh9oiYhEcO5qKoyl+njSuYPYoRARNUvFciUAwEjKX2GIdA0/lc3IuWt30dbcCP2crMQOhYioRckvKsFf0fcxtI8NWhnxr2hERA2hWKaAgb4e9PX4KwyRruGnspm4n1WIqPhMjOjXEXp6jd8JnYioJTsfmY4ShQqjBtRtFwwiIqpesVwJI0P++kKki/hnlSbuQlQ6DofEIiNXBgBo3YpTSkTUGAIvJsLR1hy9Olsg5FoqnDqao6BYgcCLifAc6iB2eEREzU6xXAEjqb7YYRBRJVjya8IuRKVjb2C0uqgAAIeCY3EhKl3EqIiIWgZHW3PsDIhEUHgy0jIK0cPeAjsDIuFoay52aEREzVJuvhx5hSXIyZdVfzIRNSoWFpqwwyGxkCtUGmNyhQqHQ2JFioiIqOVwdrCEj5cLDp2LhYG+BH/cSIOPlwucHSzFDo2IqFmKT8tBiUKFo2EJYodCROWwsNCEPbpSoSbjRESkXQ4dTKFSCVAoBYx2tWNRgYiogWTny5CRUwwACI1I46oFIh3DwkITZmVuVKtxIiLSrl9CE6ASgOH9bRF8NRW3ErPEDomIqFk6FhYPQSj9WiUIXLVApGNYWGjCpo90gn65HSCkBnqYPtJJpIiIiFqOW4lZ+O1KCqzbGMN7Ui/4eLlgZ0AkiwtERFqWnS9DaMS/PcSUSoGrFoh0DAsLTdjQ3h3Q2sgABvqlxQUrcyPM8+wFtz42IkdGRNT8Xfv7IVQqAeMH2UMikah7LiSk5YodGhFRs3IsLB6qsuUK/+CqBSLdwr0Jm7A7ydnIKyrBK0/1hpsLiwlERI1JIgH09SQY2qeDeszZwZJ9FoiItCwmNRdKpWZhQakUEJOaI1JERFQeCwtNWOiNNLQy0scTPduJHQoRUYuiUKpwPjIdrt2tYdZaKnY4RETN2pqXhgAAFm4OxoTBnfGfUbztl0jX8FaIJqpIpsBft+/jSecOMDLUFzscIqIW5XrMQ+QXlcCjX0exQyEiahEUShUUSgFGUv7cS6SLWFhoov6Kvg95iQru/WzFDoWIqMX540YaLM2M4NKlrdihEBG1CLISJQDwD2pEOoqFhSbqjxt30dHaBF1tzcUOhYioRcnKkyEiLgPDXGygV25nHiIiahgyeVlhgb++EOkifjKboLsPCxCbmguPvraQSPhDLRFRY7oQlQ5BADz6csUYEVFjUa9Y4K0QRDqJhYUmKCwiDfp6Eu4EQUTUyARBwB830tCjUxt0aNta7HCIiFoM3gpBpNu4K0QTciEqHT+HxCIzVwZDAz3cTMiEWx8WF4iIGlrgxUQ42prDQF+Ce5mFmDLUGbcSs5CQlgvPoQ5ih0dE1Oxl5BQDABQKlciREFFluGKhibgQlY69gdHIzJUBAEoUKuwNjMaFqHSRIyMiav4cbc2xMyASx84nwEiqD7PWhtgZEAlH9rkhImoUv1+/CwD4M/qeyJEQUWVYWGgiDofEQl6uQitXqHA4JFakiIiIWg5nB0v4eLkgMi4T7S1aYdeJW/DxcoGzg6XYoRERNXvZ+TLcTMgCAFyPyUBOvkzkiIioPBYWmoiM3MoTaFXjRESkXc4Olhjcqz2S7+djtKsdiwpERI3kWFg88E+/cgHA0bAEMcMhokqwsNBEWJkb1WqciIi061ZiFm4lZmHqMEcEX03FrcQssUMiImr2svNlCI1Ih1IpAACUSgGhEWlctUCkY1hYaCKmj3SC1EBzuqQGepg+0kmkiIiIWo5biVnYGRAJHy8XTBvRFT5eLtgZEMniAhFRAzsWFg+VIGiMqQSBqxaIdAwLC02EWx8bzPPspV6hYGVuhHmevbgrBBFRI0hIy9XoqVDWcyEhLVfkyIiImreY1Fz1aoUySqWAmNQckSIiospwu8kmxK2PDQsJREQiqGxLSWcHS/ZZICJqYGteGqL+ul07Mzx4kCdiNERUFa5YICIiIiIiIqI6Y2GBiIiIiIiIiOqMhQUiIiIiIiIiqjMWFoiIiIiIiIiozlhYICIiIiIiIqI6464QOkRPT1KrcdJNnK+mp7o545y2DJXNM+e+6eGcNS01mS/OactQ03nmfw9NA+epadBmDpYIgiBUfxoRERERERERUUW8FYKIiIiIiIiI6oyFBSIiIiIiIiKqMxYWiIiIiIiIiKjOWFggIiIiIiIiojpjYYGIiIiIiIiI6oyFBSIiIiIiIiKqMxYWiIiIiIiIiKjOWFggIiIiIiIiojpjYYGIiIiIiIiI6oyFBR11/PhxTJkyBf369YOnpycCAgLEDon+oVKp8MMPP2Dq1KlwdXXFuHHjsGHDBuTn56vPCQ0NxYwZM9C/f3+MGTMGu3fvFjFietTixYsxfvx4jTHOFz2K+Vd3Mf82fczBVBfMy+LSVu6NiIjAiy++CFdXV3h4eGDr1q0oKSlpzLfSotQ139Z1nlhY0EGBgYHw9fWFu7s7duzYgSFDhuDdd9/FqVOnxA6NAHz99ddYt24dRo0ahR07dmD+/PkICAjAm2++CQC4cuUKXnvtNXTt2hXbt2/H1KlT4efnh127dokcOf3yyy8ICgrSGON80aOYf3Ub82/TxhxMdcG8LD5t5N7ExER4e3vDyMgI27Ztw0svvYQ9e/Zgw4YNYr2tZq2u+bZe8ySQzhk3bpzw3//+V2PszTffFCZNmiRSRFRGpVIJgwcPFlavXq0xfuLECaFHjx7CzZs3hXnz5gkzZ87UOO7n5ycMGjRIkMlkjRkuPSI9PV0YPHiwMGLECGHcuHHqcc4XPYr5V3cx/zZtzMFUV8zL4tJW7l25cqUwcuRIjc/1vn37BGdnZyE9Pb3h30gLUp98W5954ooFHZOcnIykpCRMmDBBY3zixImIi4tDcnKySJERABQUFODpp5/GU089pTHetWtXAMDff/+N8PDwSucvNzcXV65cabRYSdOqVavg7u4ONzc39ZhMJuN8kRrzr25j/m3amIOpLpiXxaet3BsWFobRo0dDKpWqz5k0aRKUSiVCQ0Mb+F20LPXJt/WZJxYWdExcXBwAoEuXLhrjDg4OAID4+PhGj4n+ZWpqilWrVmHgwIEa47/++isAoHfv3igpKeH86ZiDBw8iKioK//d//6cxnpyczPkiNeZf3cb823QxB1NdMS+LTxu5t6ioCGlpaRXOadu2LUxNTTmPWlSffFvfeTLQQvykRXl5eQBKP8SPMjExAQCNJimkG65fv47//e9/GDduHOdPB6WmpmLDhg3YsGED2rZtq3GM80WP4n8PTQ/zr+5jDqb64H8juqm2ubeqc8rO4zxqR33zbX3niSsWdIwgCAAAiURS6bieHqdMl1y+fBkLFixAp06d8OGHH1Y5f2U4f41LEASsXLkSI0eOxMSJEys9DnC+qBTzb9PC/Kv7mIOpvpiXdU9dcu/jzhEEgfOoBdrIt/WdJ65Y0DFmZmYAKlZgCwoKNI6T+E6ePInly5fD0dERX3/9NSwtLfHw4UMAFeev7DHnr3Ht27cPt2/fxrFjx6BQKAD8m1gVCkWVnzfOV8vE/Nt0MP82DczBVF/My7qlrrm37C/glf3Fu7CwkPOoBdrIt/WdJxYWdEzZPS1JSUno2bOnejwxMVHjOIlrz5492LRpE4YMGYIdO3aoP2idO3eGvr4+kpKSNM4ve8z5a1ynT59GVlYWPDw8Khzr06cPVq9ezfkiNebfpoH5t+lgDqb6Yl7WHfXJvSYmJujQoYN63spkZGQgPz+f86gF2si39Z0nrjvRMQ4ODujUqVOFvXnPnDkDR0dHdOzYUaTIqMzBgwexceNGeHp64uuvv9ao3hkZGWHQoEE4c+aMukoIlH7YzczM4OLiIkbILdaaNWtw6NAhjX+jR4+GjY0NDh06hEmTJnG+SI35V/cx/zYtzMFUX8zLukEbudfd3R3BwcGQy+Ua5+jr62PIkCGN92aaKW3l2/rME1cs6KBFixZhxYoVaNOmDUaNGoWzZ88iMDAQ/v7+YofW4mVkZOCjjz6CnZ0dXnjhBdy8eVPjeOfOneHj44P58+fjrbfewrRp03D16lXs2rULS5cuRatWrUSKvGUq2wrpURYWFpBKpejbty8AcL5IA/Ov7mL+bXqYg0kbmJfFpa3cu2DBApw4cQILFy7EvHnzkJCQgK1bt+LZZ59lgUgLtJVv6zNPEuHRkgXpjB9//BG7d+9GWloa7O3tsXDhQnh5eYkdVosXEBCAd999t8rjfn5+eOaZZxAUFIRPP/0U8fHx6NChA1544QW89NJLjRgpVWX58uW4fPkygoKC1GOcL3oU869uYv5tHpiDqS6Yl8WjzdwbHh4OPz8/3Lp1C5aWlvDy8sKSJUtgaGjY0G+jRaprvq3rPLGwQERERERERER1xh4LRERERERERFRnLCwQERERERERUZ2xsEBEREREREREdcbCAhERERERERHVGQsLRERERERERFRnLCwQERERERERUZ2xsEA6a/ny5ejZs2e1/5YvX47Dhw+jZ8+euHTpkqgxX7p0qdIYXVxcMHbsWHz00UfIycnR2utp+33ryveRiMTF/Fs95l8iaijMwdVjDtY9BmIHQFSVWbNmwc3NTf348uXLOHDgAGbNmoWBAweqxzt37gxra2v4+fnByclJjFArGD9+PMaPH69+LJfLERkZie+//x7h4eE4ePAgDAzq//EbPHiwTr1vImoemH+rx/xLRA2FObh6zMG6h4UF0lmurq5wdXVVP1YqlThw4AAGDBiAZ555psL59vb2jRneY/Xs2bNCjDNnzoSpqSm+/vprnD59GlOmTKn369jb2+vU+yai5oH5t3rMv0TUUJiDq8ccrHt4KwRRI5o8eTIA4OrVqyJHQkTUsjD/EhGJhzm4+WNhgZqF8vdFlT2Ojo7GG2+8AVdXVwwdOhSbNm2CUqnEkSNHMHHiRAwYMADPPfccoqOjNZ4vJycH69atw/Dhw+Hi4gJPT0/s3bsXgiDUK049vdKPnEKhqNVrbd++HX379kVQUBDc3d3h6uqKgwcPVno/WFFREbZs2YIxY8bAxcUFY8aMwccff4yioiKNWDIyMrBixQoMHToUAwcOxPvvvw+5XF6v90dELQ/zL/MvEYmHOZg5WFfwVghq1hYuXIiBAwdi+fLlOHPmDHbv3o07d+7g9u3bmDdvHgRBwM6dO/HGG2/g5MmTMDAwQGFhIebMmYO0tDTMnj0bNjY2uHjxItavX4+EhAR88MEHdY7nwoULAIA+ffoAQK1eS6FQYNWqVXj55Zchl8sxcOBAXLt2TeP55XI55s+fj2vXrmH69OlwcXHBjRs38NVXX+Hy5cv49ttvYWhoCJlMhjlz5iAlJQVz585Fu3btcOTIEZw8ebLO742I6FHMv8y/RCQe5mDm4MbGwgI1awMGDIC/vz+A0iVYbm5uOH/+PI4ePYru3bsDAAoKCvDFF18gJSUFjo6O2LVrF+Lj4/Hzzz+jZ8+eAIDZs2dj69at+PLLLzFr1iz06tXrsa9bVFSEzMxM9ePMzEyEhYVh+/btsLW1VS8Hq81rqVQqzJkzBwsXLlQ/b/mk+vPPP+Pq1atYsWIFvL291c/XrVs3bN68GQcPHsTs2bNx8OBBxMXFYceOHRg3bhwA4Nlnn8XMmTORl5dXl281EZEG5l/mXyISD3Mwc3Bj460Q1KyVJQwAMDMzQ9u2beHo6KhOqADQqVMnAMCDBw8AAGfOnEGPHj3Qrl07ZGZmqv+VPVdwcHC1r7tr1y64ubmp/02ZMgXr169Hz549sWvXLpiYmNTptTw8PB77umfPnoWpqSleeOEFjfG5c+fC1NQUv/32GwDg999/h7W1tcb3p3Xr1pg5c2a1742IqCaYf0sx/xKRGJiDSzEHNx6uWKBmzdraWuOxgYEBrKysNMb09fUBlFZDASApKQnFxcUa2/w8Ki0trdrXfeaZZ+Dl5QUAkEgkMDY2hr29fYV4avta5WMvLyUlBfb29jA0NNQYl0qlsLe3R2pqKgAgNTW10k66Xbp0eezzExHVFPNvKeZfIhIDc3Ap5uDGw8ICNWtlCfNREonksdcolUoMHDgQixcvrvR4+/btq31de3t7DBs2rNrzavtaZY1vqvK4xjoqlUqdbCUSCWQyWa2uJyKqDebffzH/ElFjYw7+F3Nw42BhgagcOzs7FBQUVEiKOTk5uHDhAhwcHHT2tezs7HDt2jWUlJRoVGzlcjlSUlIwaNAgAKVL38LDw6FQKGBg8G8aSE5Orse7ISKqH+ZfIiLxMAdTfbDHAlE5Y8aMQXR0NM6dO6cxvnPnTrz55pv4+++/dfa1xowZg/z8fOzbt09jfP/+/SgoKMCoUaMAABMmTEBeXh4OHjyoPqekpAQ//fRTnd4HEZE2MP8SEYmHOZjqgysWiMp59dVXcebMGSxevBjPPfccunfvjsuXL+OXX37BiBEjMGLECJ19rZkzZ+LIkSPYuHEj7ty5AxcXF0RGRuLw4cPo37+/ujHNM888g59++gnr1q1DbGwsHB0dcfToUXXzHiIiMTD/EhGJhzmY6oOFBaJyLCwscODAAXz66ac4deoUDhw4gI4dO+L111/HwoULq73HS8zXkkql+Oabb7Bjxw4EBgbi6NGjsLGxwauvvgofHx/10jB9fX18/fXX8Pf3R2BgIAoLCzFixAh4e3vjrbfe0tr7IyKqDeZf5l8iEg9zMHNwfUgEdqogIiIiIiIiojpijwUiIiIiIiIiqjMWFoiIiIiIiIiozlhYICIiIiIiIqI6Y2GBiIiIiIiIiOqMhQUiIiIiIiIiqjMWFoiIiIiIiIiozlhYICIiIiIiIqI6Y2GBiIiIiIiIiOqMhQUiIiIiIiIiqjMWFoiIiIiIiIiozv4fpJai0sXT/jMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "time_periods = list(range(1, 41))  # 40 time periods starting from 1\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))  # Adjust the figure size as needed\n",
    "\n",
    "# Plotting Exposure Elasticity in the first subfigure\n",
    "axs[0].plot(time_periods, expo_elas[:, 1], label='Exposure Elasticity', marker='o')\n",
    "axs[0].set_xlabel('Time Period')\n",
    "axs[0].set_ylabel('Exposure Elasticity')\n",
    "axs[0].set_title('Exposure Elasticity Over Time')\n",
    "axs[0].grid(True)\n",
    "\n",
    "# Plotting Price Elasticity in the second subfigure\n",
    "axs[1].plot(time_periods, price_elas[:, 1], label='Price Elasticity', marker='x')\n",
    "axs[1].set_xlabel('Time Period')\n",
    "axs[1].set_ylabel('Price Elasticity')\n",
    "axs[1].set_title('Price Elasticity Over Time')\n",
    "axs[1].grid(True)\n",
    "\n",
    "time_periods = list(range(1, 401))  # 400 time periods starting from 1\n",
    "\n",
    "# Plotting IRF in the third subfigure\n",
    "axs[2].plot(time_periods, controls[:, 1], label='IRF', marker='^')\n",
    "axs[2].set_xlabel('Time Period')\n",
    "axs[2].set_ylabel('Impulse Response')\n",
    "axs[2].set_title('IRF Over Time')\n",
    "axs[2].grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab36438",
   "metadata": {
    "id": "eab36438"
   },
   "source": [
    "<br>\n",
    "<br>\n",
    "\n",
    "# 4 Using `LinQuadVar` in Computation <a class=\"anchor\" id=\"s2\"></a>\n",
    "\n",
    "In the previous section, we saw how to use `uncertain_expansion` to approximate variables and store their coefficients as `LinQuadVar` objects. In this section, we explore how to manipulate `LinQuadVar` objects for different uses.\n",
    "\n",
    "To aid our examples, we first extract the steady states for the state evolution processes from the previous model solution:\n",
    "\n",
    "See [src/lin_quad.py](https://github.com/lphansen/RiskUncertaintyValue/blob/main/src/lin_quad.py) for source code of `LinQuadVar` definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b5d7ee6b",
   "metadata": {
    "id": "b5d7ee6b",
    "outputId": "882669aa-3be7-4c20-dd54-1273a7c1826d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': array([[-3.93129532],\n",
       "        [ 1.00634613]])}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_J, n_X, n_W = ModelSol['var_shape']\n",
    "X0_tp1 = LinQuadVar({'c':np.array([[ModelSol['ss'][1]],[ModelSol['ss'][2]]])}, shape = (2, n_X, n_W))\n",
    "X0_tp1.coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5270666",
   "metadata": {
    "id": "a5270666"
   },
   "source": [
    "## 4.1 `LinQuadVar` Operations <a class=\"anchor\" id=\"s4-1\"></a>\n",
    "We can sum multiple LinQuads together in two different ways. Here we demonstrate this with an example by summing the zeroth, first and second order contributions of our approximation for capital growth. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "205dd148",
   "metadata": {
    "id": "205dd148",
    "outputId": "75bfa88f-3b99-4535-cf0d-121591eae233"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\log\\frac{K_{t+1}}{K_t}=-3.931+\\begin{bmatrix}0.9861&0\\end{bmatrix}X_t^1+\\begin{bmatrix}0.00011&0.00025\\end{bmatrix}W_{t+1}+\\begin{bmatrix}0.493&0\\end{bmatrix}X_t^2$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gk_tp1 = X0_tp1[0] + ModelSol['X1_tp1'][0]  + 0.5 * ModelSol['X2_tp1'][0] \n",
    "disp(gk_tp1,'\\\\log\\\\frac{K_{t+1}}{K_t}') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f901b7a3",
   "metadata": {},
   "source": [
    "In the next example, we sum together the contributions for both capital growth and technology:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b5747b24",
   "metadata": {
    "id": "b5747b24",
    "outputId": "b23404ea-0275-4a57-812e-a675e5f3d37d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'c': array([[-3.93129532],\n",
       "        [ 1.00667219]]),\n",
       " 'w': array([[0.00011, 0.00025],\n",
       "        [0.00477, 0.     ]]),\n",
       " 'x': array([[0.98609754, 0.        ],\n",
       "        [0.9628084 , 0.        ]]),\n",
       " 'xx': array([[ 0.        ,  0.        ,  0.        ,  0.        ],\n",
       "        [-0.21924476, -0.        ,  0.        ,  0.        ]]),\n",
       " 'x2': array([[0.49304877, 0.        ],\n",
       "        [0.48000401, 0.        ]])}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lq_sum([X0_tp1, ModelSol['X1_tp1'], 0.5 * ModelSol['X2_tp1']]).coeffs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1952e3",
   "metadata": {
    "id": "5c1952e3"
   },
   "source": [
    "## 4.2 `LinQuadVar` Split and Concat <a class=\"anchor\" id=\"s4-2\"></a>\n",
    "`split` breaks multiple dimensional LinQuad into one-dimensional LinQuads, while `concat` does the inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d085149f",
   "metadata": {
    "id": "d085149f",
    "outputId": "a6bdf64b-17a6-48c5-b159-5715a34e0174"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<lin_quad.LinQuadVar at 0x7f9ad3e40b80>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gk_tp1, Z_tp1 = ModelSol['X1_tp1'].split()\n",
    "concat([gk_tp1, Z_tp1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e67b04c1",
   "metadata": {
    "id": "e67b04c1"
   },
   "source": [
    "## 4.3 Evaluate a `LinQuadVar` <a class=\"anchor\" id=\"s4-3\"></a>\n",
    "We can evaluate a LinQuad at specific state $(X_{t},W_{t+1})$ in time. As an example, we evaluate all 5 variables under steady state with a multivariate random normal shock."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5782676a",
   "metadata": {
    "id": "5782676a",
    "outputId": "0e8e4471-839b-4679-c7a2-fd8d6a8e6879"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.79592804],\n",
       "       [-3.88503366],\n",
       "       [ 1.03543888],\n",
       "       [-0.00025322],\n",
       "       [-0.00712572]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1 = np.zeros([n_X ,1])\n",
    "x2 = np.zeros([n_X ,1])\n",
    "w = np.random.multivariate_normal(np.zeros(n_W),np.eye(n_W),size = 1).T\n",
    "ModelSol['JX_tp1'](*(x1,x2,w))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5051f96c",
   "metadata": {
    "id": "5051f96c"
   },
   "source": [
    "## 4.4 Next period expression for `LinQuadVar` <a class=\"anchor\" id=\"s4-4\"></a>\n",
    "`ModelSol` allows us to express a jump variable $J_t$ as a function of $t$ state and shock variables. Suppose we would like to compute its next period expression $J_{t+1}$ with shocks. The function `next_period` expresses $J_{t+1}$ in terms of $t$ state variables and $t+1$ shock variables. For example, we can express the $t+1$ expression for the first-order contribution to consumption over capital as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "7e088a47",
   "metadata": {
    "id": "7e088a47",
    "outputId": "116af9b7-9791-4391-b7a4-8d9ee133baff"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\log\\frac{C_{t+1}^1}{K_{t+1}^1}=-0.09143+\\begin{bmatrix}11.71&-1.63e-40\\end{bmatrix}X_t^1+\\begin{bmatrix}0.001306&0.002969\\end{bmatrix}W_{t+1}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmk1_tp1 = next_period(ModelSol['J1_t'][0], ModelSol['X1_tp1'])\n",
    "disp(cmk1_tp1, '\\\\log\\\\frac{C_{t+1}^1}{K_{t+1}^1}') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "fa8ab664",
   "metadata": {
    "id": "fa8ab664",
    "outputId": "f4244fa2-2a8a-4a42-8a33-2ab7407c7c8d"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\log\\frac{C_{t+1}^2}{K_{t+1}^2}=-0.02162+\\begin{bmatrix}1.136&2.455e-41\\end{bmatrix}X_t^1+\\begin{bmatrix}0.0001267&0.000288\\end{bmatrix}W_{t+1}+\\begin{bmatrix}11.71&-1.63e-40\\end{bmatrix}X_t^2+X^{1T}_{t}\\begin{bmatrix}-51.15&4.666e-39\\\\-7.867e-39&-2.048e-59\\end{bmatrix}X^1_{t}+X^{1T}_{t}\\begin{bmatrix}-0.01141&-0.02594\\\\1.007e-43&2.289e-43\\end{bmatrix}W_{t+1}+W_{t+1}^{T}\\begin{bmatrix}-6.365e-07&-1.447e-06\\\\-1.447e-06&-3.288e-06\\end{bmatrix}W_{t+1}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmk2_tp1 = next_period(ModelSol['J2_t'][0], ModelSol['X1_tp1'], ModelSol['X2_tp1'])\n",
    "disp(cmk2_tp1, '\\\\log\\\\frac{C_{t+1}^2}{K_{t+1}^2}') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daa0c4c1",
   "metadata": {
    "id": "daa0c4c1"
   },
   "source": [
    "## 4.6 Compute the Expectation of time $t+1$ `LinQuadVar` <a class=\"anchor\" id=\"s4-6\"></a>\n",
    "\n",
    "Suppose the distribution of shocks has a constant mean and variance (not state dependent). Then, we can use the `E` function to compute the expectation of a time $t+1$ `LinQuadVar` as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3140272c",
   "metadata": {
    "id": "3140272c",
    "outputId": "418a75e2-ff56-4bd7-8471-67634f26654c"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\mathbb{E}[\\log\\frac{C_{t+1}^2}{K_{t+1}^2}|\\mathfrak{F_t}]=-0.02167+\\begin{bmatrix}1.14&2.451e-41\\end{bmatrix}X_t^1+\\begin{bmatrix}11.71&-1.63e-40\\end{bmatrix}X_t^2+X^{1T}_{t}\\begin{bmatrix}-51.15&4.666e-39\\\\-7.867e-39&-2.048e-59\\end{bmatrix}X^1_{t}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E_w = ModelSol['util_sol']['μ_0']\n",
    "cov_w = np.eye(n_W)\n",
    "E_ww = cal_E_ww(E_w, cov_w)\n",
    "E_cmk2_tp1 = E(cmk2_tp1, E_w, E_ww)\n",
    "disp(E_cmk2_tp1, '\\mathbb{E}[\\\\log\\\\frac{C_{t+1}^2}{K_{t+1}^2}|\\mathfrak{F_t}]')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18a1686c",
   "metadata": {
    "id": "18a1686c"
   },
   "source": [
    "Suppose the distribution of shock has a state-dependent mean and variance (implied by $\\tilde{N}_{t+1}$ shown in the notes), we can use `E_N_tp1` and `N_tilde_measure` to compute the expectation of time $t+1$ `LinQuadVar`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a58482f7",
   "metadata": {
    "id": "a58482f7",
    "outputId": "b8f2fc85-7c37-4d01-a5ed-c6b4d6d940a0"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\mathbb{\\tilde{E}}[\\log\\frac{C_{t+1}^2}{K_{t+1}^2}|\\mathfrak{F_t}]=-0.02167+\\begin{bmatrix}1.141&2.451e-41\\end{bmatrix}X_t^1+\\begin{bmatrix}11.71&-1.63e-40\\end{bmatrix}X_t^2+X^{1T}_{t}\\begin{bmatrix}-51.16&4.666e-39\\\\-7.867e-39&-2.048e-59\\end{bmatrix}X^1_{t}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_cm = N_tilde_measure(ModelSol['util_sol']['log_N_tilde'],(1,n_X,n_W))\n",
    "E_cmk2_tp1_tilde = E_N_tp1(cmk2_tp1, N_cm)\n",
    "disp(E_cmk2_tp1_tilde, '\\mathbb{\\\\tilde{E}}[\\\\log\\\\frac{C_{t+1}^2}{K_{t+1}^2}|\\mathfrak{F_t}]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef066880",
   "metadata": {
    "id": "ef066880"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a2dbdf-d34b-4a19-b80e-5627814a7866",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "interpreter": {
   "hash": "1c67a9c81f4b1f676ee52b818d3e176be5c9dd44f82a35d6054954eb94c37bc5"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}