opt_savings.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ce97681c",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Interpolations, Optim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ba34209c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lin_interp"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Linear interpolation with constant extrapolation outside the bounds.\n",
    "\"\"\"\n",
    "function lin_interp(x, fx)\n",
    "    return linear_interpolation(x, fx, \n",
    "                extrapolation_bc = Interpolations.Flat())\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "d50af698",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "create_opt_savings_model (generic function with 1 method)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function create_opt_savings_model(; β=0.9, \n",
    "                                    γ=2.0,\n",
    "                                    R=1.2,\n",
    "                                    w_size=100,\n",
    "                                    w_max=10)\n",
    "    w_grid = LinRange(1e-8, w_max, w_size)\n",
    "    u(c) = c^(1 - γ) / (1 - γ)\n",
    "    v_init = u.(w_grid)\n",
    "    return (; β, u, R, w_grid, v_init)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f469f076",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "T"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "The Bellman operator\n",
    "\n",
    "    (Tv)(w) = min_{0 ≤ c ≤ w} { u(c) + β v(R(w - c))}\n",
    "\n",
    "\"\"\"\n",
    "function T(v, model)\n",
    "    (; β, u, R, w_grid) = model\n",
    "\n",
    "    v_new = similar(v)\n",
    "    σ = similar(v) \n",
    "    v = lin_interp(w_grid, v)\n",
    "\n",
    "    for (i, w) in enumerate(w_grid)\n",
    "        result = maximize(c -> u(c) + β * v(R * (w - c)), 0.0, w) \n",
    "        v_new[i] = Optim.maximum(result)\n",
    "        σ[i] = Optim.maximizer(result)\n",
    "    end\n",
    "\n",
    "    return (; σ, v_new)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "aa559ed3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vfi (generic function with 1 method)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function vfi(model;\n",
    "               tolerance=1e-6,    \n",
    "               max_iter=10_000,  \n",
    "               print_step=100)      \n",
    "\n",
    "    k = 0\n",
    "    error = tolerance + 1\n",
    "    v = model.v_init\n",
    "\n",
    "    while (error > tolerance) & (k <= max_iter)\n",
    "        (; v_new) = T(v, model)\n",
    "        error = maximum(abs.(v_new - v))\n",
    "        if k % print_step == 0\n",
    "            println(\"Completed iteration $k with error $error.\")\n",
    "        end\n",
    "        v = v_new\n",
    "        k += 1\n",
    "    end\n",
    "\n",
    "    if error <= tolerance\n",
    "        println(\"Terminated successfully in $k iterations.\")\n",
    "    else\n",
    "        println(\"Warning: hit iteration bound.\")\n",
    "    end\n",
    "\n",
    "    σ, v_new = T(v, model)\n",
    "    return σ, v_new\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "75032b6e",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_greedy"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\" Get a v-greedy policy \"\n",
    "function get_greedy(v, model)\n",
    "    (; β, u, R, w_grid) = model\n",
    "    σ, v_new = T(v, model)\n",
    "    return σ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "cda2be6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "T_σ"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "The policy operator\n",
    "\n",
    "    (T_σ v)(w) =  u(σ(w)) + β v[R(w - σ(w))]\n",
    "\n",
    "\"\"\"\n",
    "function T_σ(v, σ, model)\n",
    "    (; β, u, R, w_grid) = model\n",
    "\n",
    "    v_new = similar(v)\n",
    "    v = lin_interp(w_grid, v)\n",
    "    for (i, w) in enumerate(w_grid)\n",
    "        c = σ[i]\n",
    "        v_new[i] = u(c) + β * v(R * (w - c))\n",
    "    end\n",
    "    return v_new\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "0e12f36c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_value"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"Approximate lifetime value of policy σ.\"\n",
    "function get_value(v_init, σ, m, model)\n",
    "    v = v_init\n",
    "    for i in 1:m\n",
    "        v = T_σ(v, σ, model)\n",
    "    end\n",
    "    return v\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "4b859670",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "policy_iteration"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"Howard policy iteration routine.\"\n",
    "function policy_iteration(model; \n",
    "                          tolerance=1e-6, \n",
    "                          max_iter=1_000,\n",
    "                          m=10,\n",
    "                          print_step=10)\n",
    "    v = model.v_init\n",
    "    error = tolerance + 1\n",
    "    k = 1\n",
    "    while error > tolerance && k < max_iter\n",
    "        last_v = v\n",
    "        σ = get_greedy(v, model)\n",
    "        v = get_value(v, σ, m, model)\n",
    "        error = maximum(abs.(v - last_v))\n",
    "        if k % print_step == 0\n",
    "            println(\"Completed iteration $k with error $error.\")\n",
    "        end\n",
    "        k += 1\n",
    "    end\n",
    "    return get_greedy(v, model), v\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "8b6196b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "using LaTeXStrings\n",
    "fontsize=12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "999fd2b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = create_opt_savings_model()\n",
    "(; β, u, R, w_grid, v_init) = model;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "8c685069",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving via VFI. \n",
      "\n",
      "\n",
      "Completed iteration 0 with error 9.000000634421515e7.\n",
      "Completed iteration 100 with error 2390.5260684490204.\n",
      "Completed iteration 200 with error 0.06349575519561768.\n",
      "Completed iteration 300 with error 1.6689300537109375e-6.\n",
      "Terminated successfully in 306 iterations.\n",
      "  0.171449 seconds (139.28 k allocations: 11.239 MiB, 73.90% compilation time)\n",
      "elapsed time (ns):  171449305\n",
      "gc time (ns):       0\n",
      "bytes allocated:    11784958\n",
      "pool allocs:        139169\n",
      "non-pool GC allocs: 113\n",
      "minor collections:  0\n",
      "full collections:   0\n"
     ]
    }
   ],
   "source": [
    "println(\"Solving via VFI. \\n\\n\")\n",
    "@timev σ_star, v_star = vfi(model);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "fec4fee2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Solving via HPI.\n",
      "\n",
      "\n",
      "Completed iteration 10 with error 44654.35739028454.\n",
      "Completed iteration 20 with error 1.1860822439193726.\n",
      "Completed iteration 30 with error 3.159046173095703e-5.\n",
      "  0.098984 seconds (38.95 k allocations: 3.262 MiB, 93.19% compilation time)\n",
      "elapsed time (ns):  98983597\n",
      "gc time (ns):       0\n",
      "bytes allocated:    3420697\n",
      "pool allocs:        38889\n",
      "non-pool GC allocs: 58\n",
      "minor collections:  0\n",
      "full collections:   0\n"
     ]
    }
   ],
   "source": [
    "println(\"\\n\\nSolving via HPI.\\n\\n\")\n",
    "@timev σ_star_hpi, v_star_hpi = policy_iteration(model);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "51e39b40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApsAAAG8CAYAAACR0253AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABc0UlEQVR4nO3dd3RUdcLG8e+dSTLpk94gNAUNHQkiVVBAXEGxYC+rgCCIlFWRVVexADZElB5WURd11VXZV11lraAgiIAgCqKUACE9M6mTMvf9A80S6ZDkpjyfc+acd+7c8szk3cPj797fvYZpmiYiIiIiIjXAZnUAEREREWm4VDZFREREpMaobIqIiIhIjVHZFBEREZEao7IpIiIiIjVGZVNEREREaozKpoiIiIjUGJVNEREREakxPlYH+COv18v+/fsJCQnBMAyr44iIiIjIH5imSX5+PgkJCdhsxx67rHNlc//+/SQmJlodQ0RERESOIzU1laZNmx5znTpXNkNCQoCD4UNDQy1OIyIiIiJ/5Ha7SUxMrOxtx1Lnyubvp85DQ0NVNkVERETqsBO55FEThERERESkxqhsioiIiEiNUdkUERERkRqjsikiIiIiNabOTRA6WRUVFZSVlVkdo0Hw9fXFbrdbHUNEREQakHpbNk3T5MCBA7hcLkzTtDpOg2AYBk6nk7i4ON1QX0RERKpFvS2bLpeLvLw8oqOjCQoKUjk6TaZpUlhYSGZmJgEBAYSFhVkdSURERBqAelk2TdMkIyOD0NBQoqKirI7TYAQEBODxeMjIyMDpdKrAi4iIyGmrlxOEKioqqKio0E3fa0BoaGjl7ysiIiJyuupl2SwvLwfAx6deDszWab//pr//xiIiIiKno16Wzd/pNG/1028qIiIi1emky+aXX37J0KFDSUhIwDAM3n333Sqfm6bJww8/TEJCAgEBAfTr148ffvihuvKKiIiISD1y0mWzsLCQTp068cILLxzx8yeffJJZs2bxwgsvsG7dOuLi4hg4cCD5+fmnHVZERERE6peTLpsXX3wxjz32GFdcccVhn5mmyezZs7n//vu54ooraN++PUuXLqWoqIhly5ZVS+DGKCUlBcMwCA4OrrL8888/xzCMI77WrFljUVoRERGR/6nWGTY7d+7kwIEDDBo0qHKZw+Hg/PPP5+uvv2b06NGHbePxePB4PJXv3W53dUaq9/bt28fdd99NQkICLpfriOtMnz6d/v37V1nWvn372ognIiIiFsvIyAAgJibG4iRHVq0ThA4cOABAbGxsleWxsbGVn/3RjBkzcDqdla/ExMTqjFTvjRkzhr59+zJw4MCjrtO6dWvOO++8Kq8/joKKiIhIw2KaJq+88gpJSUmMHTvW6jhHVSOz0f84o9k0zaPOcp46dSoul6vylZqaWhOR6qVXX32VL774gnnz5lkdRUREROqQ3bt3c/HFF3PzzTeTk5PDjh07jnoG1GrVWjbj4uIADhvFzMjIOGy083cOh4PQ0NAqr9NRWFh41FdJSckJr1tcXHzK6xYVFZ3Wd4CDv9nEiROZOXMmTZs2Pea648aNw8fHh9DQUC666CJWrVp12scXERGRuqeiooI5c+bQrl07PvroIxwOB9OnT2fdunU4nU6r4x1RtZbNli1bEhcXx4oVKyqXlZaW8sUXX9CzZ8/qPNRRBQcHH/V15ZVXVlk3JibmqOtefPHFVdZt0aLFUdft27dvlXXbtm172t9j7NixnHXWWdxxxx1HXcfpdDJhwgQWLlzIZ599xnPPPUdqair9+vXjo48+Ou0MIiIiUnekpqbSp08fJkyYQGFhIb1792bTpk1MnToVX19fq+Md1UlPECooKGDHjh2V73fu3MnGjRuJiIigWbNmTJw4kenTp9O6dWtat27N9OnTCQwM5Prrr6/W4A3Z22+/zb///W82bNhwzJusd+nShS5dulS+79OnD5dffjkdOnTg3nvv5aKLLqqNuCIiIlILwsPD2bdvHyEhITzxxBOMHj0am63uP5/npMvmt99+W2Xm8+TJkwG45ZZbeOmll7j33nspLi5m7Nix5Obm0r17dz7++GNCQkKqL/UxFBQUHPUzu91e5f3vs7eO5I9/vF27dp3wulu3bj1GwmMrKChg3LhxjB8/noSEBPLy8oCDI8QAeXl5+Pr6EhQUdMTtw8LCGDJkCAsWLKC4uJiAgIBTziIiIiLW2rx5M+3atcNmsxEcHMybb75JfHx8vZpQbZimaVod4lButxun04nL5Trq9ZslJSXs3LmTli1b4u/vX8sJa9auXbto2bLlMde57LLLDnty06HGjBnDwoULKS4uPunfpyH/tiIiIvVFYWEhDz74IM899xzz5s074u0jrXQife131XqfTTl9cXFxfPbZZ4ctnzlzJl988QUffvghUVFRR90+NzeX//u//6Nz584qiyIiIvXQJ598wqhRo9i5cycAmzZtsjjR6VHZrGP8/f3p16/fYctfeukl7HZ7lc+uv/56mjVrRnJyMlFRUfz8888888wzpKen89JLL9VaZhERETl9ubm53H333fz9738HIDExkQULFvCnP/3J4mSnR2WzHuvYsSNvvPEGCxYsoKCggIiICHr37s0rr7xCt27drI4nIiIiJ2jFihXcfPPNHDhwAMMwGDduHNOnT6+1OS81SddsShX6bUVERGrft99+S/fu3WnTpg0pKSn06tXL6kjHdDLXbNb9+fIiIiIiDYxpmlWuxUxOTub9999nw4YNdb5oniyVTREREZFa9OuvvzJw4EDOPfdcfvzxx8rlgwcPbpBnFVU2RURERGpBRUUFs2bNon379nzyySfY7Xa+//57q2PVOE0QEhEREalh33//PSNHjmTdunUA9O/fn8WLF3PGGWdYnKzmaWRTREREpAY9/vjjdO3alXXr1uF0OklJSeGTTz5pFEUTNLIpIiIiUqN8fX0pLy9n2LBhzJ07l4SEBKsj1SqVTREREZFqlJ+fT1paGm3atAFg8uTJdOzYkcGDB1uczBo6jS4iIiJSTT788EPatWvHZZddRklJCQA+Pj6NtmiCyqaIiIjIacvKyuLGG2/kT3/6E6mpqZSUlLB7926rY9UJKpsiIiIip8g0TV5//XWSkpL4xz/+gc1mY9KkSWzZsoWzzjrL6nh1gspmHXL55ZcTEBBAXl7eUde54YYb8PX1JT09HcMwjviKioqqXP/hhx/GMAyysrJq4RuIiIg0HgUFBVx66aVcd911ZGVl0b59e1avXs2sWbMICgqyOl6dobJZh4wYMYKSkhKWLVt2xM9dLhfvvPMOQ4YMITY2FoCrrrqK1atXV3l99NFHtRlbRESkUQoKCsLj8eDr68u0adNYv3495557rtWx6hzNRq9DLr74YhISEvj73//O2LFjD/v8tddeo7i4mBEjRlQui42N5bzzzqvNmCIiIo3W9u3biY2Nxel0YhgGixcvpqCggHbt2lkdrc7SyGYdYrfbueWWW1i/fj2bN28+7PMXX3yR+Ph4Lr74YgvSiYiINF5lZWXMnDmTjh07MmXKlMrlzZs3V9E8DpXNOua2227DMAz+/ve/V1m+detW1q5dyy233ILdbq9cbpom5eXlVV6madZ2bBERkQbru+++49xzz2Xq1Kl4PB527dpFWVmZ1bHqjQZ1Gn3o86vIzPdYHaNSdIiDf4/vfVLbnHnmmfTt25dXX32VJ598El9fX4DK8nnbbbdVWX/evHnMmzevyrLFixczcuTI00guIiIixcXFPPzwwzzzzDNUVFQQERHBs88+y0033YRhGFbHqzcaVNnMzPdwwF1idYzTNmLECG6++WaWL1/OlVdeSXl5Oa+++ip9+vShdevWVda9+uqrueeee6osa9GiRS2mFRERaXg2b97MlVdeyc8//wzANddcw3PPPVc5QVdOXIMqm9EhDqsjVHGqea666irGjx/Piy++yJVXXskHH3xAeno6TzzxxOHHiI4mOTn5dKOKiIjIIeLi4sjJySEhIYH58+dz6aWXWh2p3mpQZfNkT1nXVQEBAVx33XUsXryYtLQ0/v73vxMSEsLw4cOtjiYiItJgrVmzhu7du2MYBtHR0fzf//0fSUlJOJ1Oq6PVa5ogVEeNGDGCiooKnnrqKT744AOuvfZaAgMDrY4lIiLS4KSnp3PNNdfQo0cP3n777crl5513nopmNWhQI5sNSXJyMh07dmT27NmYplnl3poiIiJy+kzT5JVXXmHSpEnk5ORgt9v55ZdfrI7V4Ghksw4bMWIEpmnStm1bunfvbnUcERGRBmPXrl1cfPHF3HLLLeTk5NC5c2fWrl1b5R6aUj0Ms47dlNHtduN0OnG5XISGhh5xnZKSEnbu3EnLli3x9/ev5YQNm35bERFp6F555RXuuOMOCgsLcTgcPPTQQ9x9992VtxuU4zuRvvY7nUYXERGRRiU2NpbCwkJ69+5NSkoKZ511ltWRGjSdRhcREZEGrbS0lG+//bby/aBBg/jvf//LF198oaJZC1Q2RUREpMFau3YtXbt2pX///qSmplYuv/DCC7HZVINqg35lERERaXAKCwuZPHkyPXr0YMuWLfj7+2umuUVUNkVERKRB+eSTT+jQoQPPPvssXq+Xm266iR9//JF+/fpZHa1R0gQhERERaRBM02TMmDEsWrQIgGbNmrFw4UIGDx5scbLGTSObIiIi0iAYhoHT6cQwDO688062bNmiolkHaGRTRERE6q20tDQKCws588wzAXj44Ye54oorOO+88yxOJr/TyKaIiIjUO6ZpsmTJEpKSkrjxxhupqKgAIDAwUEWzjlHZFBERkXrl119/ZeDAgYwcORKXy0VFRQVZWVlWx5KjUNkUERGReqG8vJxnnnmG9u3b88knnxAQEMDTTz/N6tWriY2NtTqeHIXKZh300ksvYRhGlacdHGrIkCG0aNGi8n2LFi0wDKPyFRwcTPfu3Xn55ZerbNevXz/at29fk9FFRERqxIEDB+jRowd33303xcXF9O/fn82bN/OXv/wFHx9NQanLVDYbiF69erF69WpWr15dWVZvueUW5s+fb3U0ERGR0xYVFYVpmjidTlJSUvjkk08444wzrI4lJ0D/KdBAhIWFVbkgesCAATRv3pxZs2Zxxx13WJhMRETk1HzzzTd06tQJf39/fHx8WLZsGcHBwSQkJFgdTU6CRjYbqLCwMM466yx2795tdRQREZGTkp+fz/jx4+nRowePPvpo5fI2bdqoaNZDGtmswyoqKigvLz9suWmax922rKyM3bt3Ex0dXRPRREREasSHH37ImDFj2LNnDwBZWVmYpolhGBYnk1PV8Mrm1y/A6rnHXy++E1z/etVly66FtE3H37bHOOh55//ee/LhhXOPvc4pONZ9wpo3b17lvWmalcV07969PPzww2RkZHDPPfecVgYREZHakJWVxaRJk3j11VcBaNmyJYsWLWLAgAEWJ5PT1fDKpicf8vcffz1nk8OXFWWd2Lae/KrvTfPw7f64zil4+eWXSUpKOmz5pEmTSE1NrbLsgw8+wNfXt/J9QEAA48eP57HHHjvtHCIiIjXps88+4+qrryYrKwubzcbEiRN55JFHCAoKsjqaVIOGVzYdIRByAtdzBEYdedmJbOsIqfreMA7f7o/rnIKkpCSSk5MPW+50Og8rm7179+bZZ5/FMAwCAwM544wz8PPzO+0MIiIiNa1ly5YUFxfTvn17lixZwrnnnnv8jaTeaHhls+edp376+o+n1U+UIwT+8uOpbVtNnE7nEYupiIhIXeP1evnyyy/p168fcPB+0Z9++imdO3fWQEkDpNnoIiIiUmu2b99O//796d+/P5999lnl8nPPPVdFs4FS2RQREZEaV1ZWxsyZM+nYsSNffvklQUFB7Nu3z+pYUgsa3ml0ERERqVO+++47RowYwcaNGwEYNGgQCxcurPLoZWm4DPNEbtpYi9xuN06nE5fLRWho6BHXKSkpYefOnbRs2RJ/f/9aTtiw6bcVEZHqNHPmTB544AEqKioIDw9n9uzZ3HTTTbpvZj13In3tdxrZFBERkRrTvHlzKioquPrqq5kzZw6xsbFWR5JaprIpIiIi1cblcrF9+3a6desGwLXXXkvz5s3p2bOnxcnEKpogJCIiItVi+fLltG3bliFDhpCTkwOAYRgqmo2cyqaIiIiclvT0dK655houu+wy9u/fT2hoqGaaSyWVTRERETklpmlWPlr5n//8J3a7nXvvvZfvv/+eDh06WB1P6oh6fc1mHZtI3yDoNxURkRNRWlrKpZdeykcffQRA586dWbJkCeecc47FyaSuqZcjm76+vgAUFRVZnKTh+f03/f03FhERORI/Pz+aNGmCw+FgxowZrF27VkVTjqhe3mcTIC0tjby8PGJiYggMDNT9uk6TaZoUFRWRkZFBWFgY8fHxVkcSEZE6ZuvWrYSEhJCYmAhAbm4uGRkZnHXWWRYnk9rWKO6zGRcXB0BGRobFSRqWsLCwyt9WREQEDp4ynzlzJo899hgDBgzg/fffxzAMwsPDCQ8Ptzqe1HH1tmwahkF8fDwxMTGUlZVZHadB8PX1xW63Wx1DRETqkLVr1zJixAi2bNkCgN1up6ioiKCgIIuTSX1Rb8vm7+x2uwqSiIhINSssLOTBBx/kueeew+v1Eh0dzZw5c7jmmmt06ZqclGqfIFReXs4DDzxAy5YtCQgIoFWrVjzyyCN4vd7qPpSIiIjUgG3bttGhQweeffZZvF4vN954I1u3buXaa69V0ZSTVu0jm0888QQLFixg6dKltGvXjm+//ZZbb70Vp9PJhAkTqvtwIiIiUs2aNWuGj48PiYmJLFy4kIsvvtjqSFKPVXvZXL16NZdddhmXXHIJAC1atOC1117j22+/re5DiYiISDX5+OOPufDCC7Hb7QQEBLB8+XKaNGlCSEiI1dGknqv20+i9e/fmk08+Yfv27QBs2rSJVatW8ac//emI63s8Htxud5WXiIiI1I60tDSuuOIKLrroIubMmVO5/Oyzz1bRlGpR7SObU6ZMweVycfbZZ2O326moqODxxx/nuuuuO+L6M2bMYNq0adUdQ0RERI7BNE2WLFnC3XffjcvlwsfHh+LiYqtjSQNU7SObb7zxBq+++irLli3ju+++Y+nSpTz99NMsXbr0iOtPnToVl8tV+UpNTa3uSCIiInKIX375hQEDBjBq1ChcLhfJycmsX7+ev/71r1ZHkwao2p8glJiYyH333ce4ceMqlz322GO8+uqr/PTTT8fd/mTuSC8iIiIn5/XXX+e2226juLiYgIAAHn30USZMmICPT72/G6LUIkufIFRUVITNVnXA1G6369ZHIiIidUC7du0oKyvjggsuYNGiRZxxxhlWR5IGrtrL5tChQ3n88cdp1qwZ7dq1Y8OGDcyaNYvbbrutug8lIiIix1FSUsKXX37JoEGDAOjQoQPffPMNXbp00T0zpVZU+2n0/Px8HnzwQd555x0yMjJISEjguuuu429/+xt+fn7H3V6n0UVERKrHV199xciRI/n5559Zt24dXbp0sTqSNBCWnkYPCQlh9uzZzJ49u7p3LSIiIicgPz+fqVOnMm/ePEzTJC4ujuzsbKtjSSNV7bPRRURExDoffvgh7dq1Y+7cuZimyW233cbWrVsZMGCA1dGkkdLUMxERkQbijjvuYMGCBQC0bNmSRYsWqWSK5TSyKSIi0kC0a9cOm83G5MmT2bx5s4qm1AnVPkHodGmCkIiIyInZu3cv6enpdO3aFQCv18vmzZvp1KmTxcmkoTuZvqaRTRERkXrG6/Uyf/582rZty9VXX01hYSEANptNRVPqHJVNERGRemT79u3069ePsWPHkp+fT0xMDDk5OVbHEjkqlU0REZF6oKysjBkzZtCxY0dWrlxJUFAQzz33HKtWrSIxMdHqeCJHpdnoIiIidVxeXh79+/dn48aNAFx00UUsXLiQ5s2bWxtM5ARoZFNERKSOczqdNGvWjIiICF5++WU+/PBDFU2pNzSyKSIiUgd98cUXtGvXjqioKAzDYOHChdhsNmJiYqyOJnJSNLIpIiJSh7hcLkaPHk2/fv2YNGlS5fK4uDgVTamXNLIpIiJSRyxfvpw77riD/fv3AxAUFERFRQV2u93iZCKnTmVTRETEYunp6dx1113885//BODMM88kJSWF888/3+JkIqdPZVNERMRCX331FZdeeik5OTnY7XbuvvtuHnroIQICAqyOJlItVDZFREQs1LZtW/z8/OjcuTNLlizhnHPOsTqSSLXSBCEREZFaVFFRwb/+9S9M0wQgPDyczz77jLVr16poSoOksikiIlJLtm7dSu/evbnyyit59dVXK5efffbZ+Pr6WphMpOaobIqIiNSw0tJSHnnkEbp06cKaNWsICQnB6/VaHUukVuiaTRERkRq0du1aRowYwZYtWwAYMmQI8+bN0/PMpdHQyKaIiEgNmTVrFj169GDLli1ERUWxbNkyli9frqIpjYrKpoiISA3p1q0bpmly00038eOPP3LddddhGIbVsURqlU6ji4iIVJPc3FzWrl3LRRddBECfPn3YsmULbdu2tTiZiHU0sikiIlIN3n77bZKSkrj88sv55ZdfKperaEpjp7IpIiJyGvbv388VV1zBVVddRXp6Os2bNyc/P9/qWCJ1hsqmiIjIKTBNk5SUFNq2bcs777yDj48PDz74IBs3bqRz585WxxOpM3TNpoiIyEnyer1ccskl/Oc//wEOTgRKSUmhY8eOFicTqXs0sikiInKSbDYb3bt3JyAggGeeeYbVq1eraIochWH+/nDWOsLtduN0OnG5XISGhlodR0REBIDvv/8ewzDo0KEDAB6Ph3379tGqVSuLk4nUvpPpaxrZFBEROYaSkhIeeOABunbtys0330xZWRkADodDRVPkBOiaTRERkaP46quvGDlyJD/99BMArVq1orCwkLCwMGuDidQjGtkUERH5g/z8fO6880569+7NTz/9RGxsLG+99RZvv/22iqbISdLIpoiIyCF2795Nnz59SE1NBeC2227jqaeeIiIiwuJkIvWTyqaIiMghEhMTad68OT4+PixatIgBAwZYHUmkXlPZFBGRRs00Td5++20GDx5McHAwNpuN119/nbCwMIKCgqyOJ1Lv6ZpNERFptFJTUxk6dCjDhw/nwQcfrFzepEkTFU2RaqKyKSIijY7X62X+/Pm0a9eO999/Hz8/P6KioqyOJdIg6TS6iIg0Ktu2bWPUqFGsXLkSgB49elQ+41xEqp9GNkVEpFEwTZN//etfdOrUiZUrVxIUFMScOXNYuXKliqZIDdLIpoiINBrnnXce/v7+9OvXj4ULF9K8eXOrI4k0eBrZFBGRBqu4uJjXX38dAMMwSEhIYP369Xz44YcqmiK1RGVTREQapM8//5yOHTty3XXX8eGHH1YuP+OMMzAMw8JkIo2LyqaIiDQoLpeL0aNH079/f3bs2EGTJk2w2+1WxxJptFQ2RUSkwXjvvfdo27YtixYtAmD06NH88MMPDBo0yOJkIo2XJgiJiEiDMGHCBObMmQNA69atWbx4Meeff77FqUREI5siItIgXHjhhdjtdqZMmcKmTZtUNEXqCI1siohIvbRr1y62b99eeYr80ksvZfv27bRq1criZCJyKI1siohIvVJRUcFzzz1H+/btufbaa0lPT6/8TEVTpO7RyKaIiNQbW7duZcSIEaxZswaArl27UlxcbHEqETkWjWyKiEidV1payrRp0+jcuTNr1qwhJCSE+fPn89lnn9GiRQur44nIMWhkU0RE6rSSkhK6d+/O999/D8CQIUOYP38+TZs2tTiZiJwIjWyKiEid5u/vT58+fYiOjua1115j+fLlKpoi9YhhmqZpdYhDud1unE4nLpeL0NBQq+OIiIgFVqxYQcuWLTnzzDMByM/Px+PxEBUVZXEyEYGT62sa2RQRkTojJyeHW2+9lUGDBjFq1Ch+Hw8JCQlR0RSpp1Q2RUTEcqZp8uabb9K2bVteeuklDMOgffv2lJaWWh1NRE6TJgiJiIil9u/fz7hx43j33XcBOPvss1myZAk9e/a0NpiIVAuVTRERscx3333HBRdcgMvlwsfHh6lTp3L//ffjcDisjiYi1URlU0RELNO+fXsSExNp06YNKSkpdOzY0epIIlLNdM2miIjUmvLyclJSUigrKwPAz8+P//znP6xevVpFU6SBUtkUEZFasWnTJs477zxGjRrFk08+Wbm8SZMm2O12C5OJSE1S2RQRkRpVUlLC/fffT3JyMuvXr8fpdOqm7CKNiK7ZFBGRGrNq1SpGjRrFTz/9BMDll1/O3LlziY+PtziZiNSWGhnZ3LdvHzfeeCORkZEEBgbSuXNn1q9fXxOHEhGROmru3Ln07duXn376idjYWN566y3+9a9/qWiKNDLVPrKZm5tLr1696N+/Px9++CExMTH88ssvhIWFVfehRESkDhs4cCAOh4Prr7+ep59+mvDwcKsjiYgFqr1sPvHEEyQmJvLiiy9WLmvRokV1H0ZEROqYrKwsPvroI2644QYA2rRpw/bt20lMTLQ4mYhYqdpPoy9fvpzk5GSGDx9OTEwMXbp0YfHixdV9GBERqSNM0+S1114jKSmJm266idWrV1d+pqIpItVeNn/99Vfmz59P69at+eijjxgzZgx33XUXL7/88hHX93g8uN3uKi8REakfUlNTGTp0KNdffz1ZWVm0a9dOT/8RkSoM0zTN6tyhn58fycnJfP3115XL7rrrLtatW1flv3Z/9/DDDzNt2rTDlrtcLkJDQ6szmoiIVBOv18uCBQu47777yM/Px8/PjwceeIApU6bg5+dndTwRqWFutxun03lCfa3aRzbj4+Np27ZtlWVJSUns2bPniOtPnToVl8tV+UpNTa3uSCIiUs2GDRvGuHHjyM/Pp0ePHmzYsIEHH3xQRVNEDlPtZbNXr15s27atyrLt27fTvHnzI67vcDgIDQ2t8hIRkbpt2LBhBAUFMWfOHFauXHnYIIOIyO+qfTb6pEmT6NmzJ9OnT+fqq69m7dq1LFq0iEWLFlX3oUREpJasX7+e/Px8+vXrB8Ctt97K4MGDSUhIsDaYiNR51T6y2a1bN9555x1ee+012rdvz6OPPsrs2bMrb4UhIiL1R3FxMVOmTOHcc8/lxhtvrJzEaRiGiqaInJAaeVzlkCFDGDJkSE3sWkREasnnn3/OqFGj2LFjBwC9e/emrKzM4lQiUt/UyOMqRUSk/nK5XIwePZr+/fuzY8cOmjRpwvLly3n99deJjIy0Op6I1DM1MrIpIiL1U2ZmJp07d2b//v0AjBkzhpkzZ+J0Oi1OJiL1lcqmiIhUio6Opk+fPnz33XcsXryY888/3+pIIlLPqWyKiDRipmny6quvMmDAAOLj4wGYP38+/v7+BAQEWJxORBoCXbMpItJI7dq1i8GDB3PzzTczfvz4yuXh4eEqmiJSbVQ2RUQamYqKCp577jnat2/Pxx9/jMPhIDk5Ga/Xa3U0EWmAdBpdRKQR+eGHHxg5ciRr1qwBoG/fvixevJg2bdpYnExEGiqVTRGRRuI///kPl156KWVlZYSGhvLkk08yatQobDad5BKRmqOyKSLSSPTu3Zv4+Hg6derEvHnzaNq0qdWRRKQR0H/Oiog0UIWFhcyZM6fyWszg4GDWrl3Le++9p6IpIrVGI5siIg3QihUruP3229m1axcOh4PRo0cDEBsba3EyEWlsNLIpItKA5OTkcOuttzJo0CB27dpFs2bNaNmypdWxRKQRU9kUEWkATNPkrbfeom3btrz00ksYhsH48ePZsmULgwYNsjqeiDRiOo0uIlLPmabJfffdx5NPPglAUlISS5YsoUePHhYnExHRyKaISL1nGAbDhw/H39+fBx98kA0bNqhoikidoZFNEZF66JdffuGbb77h+uuvByA5OZndu3cTExNjcTIRkao0sikiUo+Ul5fz9NNP06FDB2699VZ+/PHHys9UNEWkLtLIpohIPbFp0yZGjhzJt99+C8AFF1yAv7+/xalERI5NI5siInVcSUkJDzzwAMnJyXz77bc4nU5SUlL473//q9saiUidp5FNEZE6zOv10rNnTzZs2ADAFVdcwQsvvEB8fLzFyURETozKpohIHWaz2bj22mtJS0tj7ty5XHHFFVZHEhE5KYZpmqbVIQ7ldrtxOp24XC5CQ0OtjiMiUus++OADIiMj6d69O3BwUlB+fj7h4eEWJxMROehk+pqu2RQRqSMyMzO58cYbueSSS/jzn/+Mx+MBwMfHR0VTROotnUYXEbGYaZq89tprTJgwgaysLGw2G5dccgler9fqaCIip01lU0TEQqmpqdxxxx28//77AHTo0IElS5bQrVs3i5OJiFQPlU0REYts376d5ORk8vPz8fPz44EHHmDKlCn4+flZHU1EpNqobIqIWKR169acd955FBQUkJKSQtu2ba2OJCJS7TRBSESklpSVlfHss8/icrkAMAyDN954g5UrV6poikiDpbIpIlIL1q9fT7du3Zg8eTJTpkypXB4eHo7dbrcwmYhIzVLZFBGpQUVFRdx7772ce+65bNq0iYiICHr37m11LBGRWqNrNkVEasjnn3/OqFGj2LFjBwDXXnstzz33HDExMRYnExGpPSqbIiI14MUXX+S2224DoEmTJsyfP5+hQ4danEpEpPapbIqI1IAhQ4YQHR3NFVdcwRNPPIHT6bQ6koiIJVQ2RUSqQXp6Ov/4xz+YPHkyANHR0fz0009ERERYnExExFqaICQichpM0+Sll14iKSmJv/zlL7z99tuVn6loiohoZFNE5JTt2rWL22+/nRUrVgDQpUsXzjjjDItTiYjULRrZFBE5SRUVFTz33HO0a9eOFStW4O/vz8yZM1m7di2dO3e2Op6ISJ2ikU0RkZNgmiY33HADb7zxBgB9+/Zl8eLFtGnTxuJkIiJ1k0Y2RUROkGmaANx2222EhoayYMECPvvsMxVNEZFj0MimiMhxfPPNN+zZs4fhw4cDMGjQIHbv3k1YWJi1wURE6gGNbIqIHEVhYSGTJk2iR48e3HbbbaSmplZ+pqIpInJiNLIpInIEH3/8MaNHj2bXrl0AXH755QQGBlobSkSkHlLZFBE5RE5ODpMnT2bp0qUANGvWjIULFzJ48GCLk4mI1E8qmyIivyksLKRDhw7s378fwzC48847mT59OsHBwVZHExGpt3TNpojIb4KCgrjhhhtISkriq6++Ys6cOSqaIiKnyTB/v5dHHeF2u3E6nbhcLkJDQ62OIyINmGmaLFmyhB49etCuXTsAiouLsdlsOBwOi9OJiNRdJ9PXNLIpIo3Sjh07uPDCCxk1ahQjRoygoqICgICAABVNEZFqpLIpIo1KeXk5Tz/9NB07duSzzz4jICCAa665xupYIiINliYIiUijsWnTJkaMGMH69esBuPDCC1m0aBGtWrWyOJmISMOlsikijcKqVavo378/5eXlhIWFMWvWLP785z9jGIbV0UREGjSVTRFpFM477zw6depEixYteP7554mPj7c6kohIo6BrNkWkQcrPz2fatGmUlJQA4OPjw6effspbb72loikiUos0sikiDc4HH3zAmDFjSE1NpbS0lMcffxxAt1MTEbGARjZFpMHIysrixhtv5JJLLiE1NZVWrVpxwQUXWB1LRKRRU9kUkXrPNE2WLVtGUlIS//jHP7DZbEyePJnNmzdz4YUXWh1PRKRR02l0Ean3Hn30UR566CEAOnTowJIlS+jWrZvFqUREBDSyKSINwM0330xERASPPvoo3377rYqmiEgdopFNEal3tm3bxgcffMCkSZMAaNGiBbt27SIkJMTiZCIi8kca2RSReqOsrIzp06fTqVMnJk+ezGeffVb5mYqmiEjdpJFNEakX1q9fz4gRI9i0aRMAgwcP1mMmRUTqAY1sikidVlRUxJQpU+jevTubNm0iIiKCV155hQ8++IDmzZtbHU9ERI5DI5siUmeZpsnAgQP5+uuvAbj22mt57rnniImJsTiZiIicqBof2ZwxYwaGYTBx4sSaPpSINDCGYXDnnXfSpEkTli9fzmuvvaaiKSJSz9ToyOa6detYtGgRHTt2rMnDiEgD8u6772K32xk6dChwcDRz6NChBAcHW5xMRERORY2NbBYUFHDDDTewePFiwsPDa+owItJApKenc/XVV3P55ZczcuRIcnJygIOjmyqaIiL1V42VzXHjxnHJJZcwYMCAY67n8Xhwu91VXiLSeJimyUsvvURSUhJvvvkmdrud2267jYCAAKujiYhINaiR0+ivv/463333HevWrTvuujNmzGDatGk1EUNE6rhdu3Zx++23s2LFCgC6dOnCkiVL6NKli8XJRESkulT7yGZqaioTJkzg1Vdfxd/f/7jrT506FZfLVflKTU2t7kgiUgelpaXRvn17VqxYgb+/PzNnzmTt2rUqmiIiDYxhmqZZnTt89913ufzyy7Hb7ZXLKioqMAwDm82Gx+Op8tkfud1unE4nLpeL0NDQ6owmInXMiBEj+OWXX1i8eDGtW7e2Oo6IiJygk+lr1V428/Pz2b17d5Vlt956K2effTZTpkyhffv2x9xeZVOkYfJ4PDz55JPccsstNGvWDIDi4mIcDgc2m54vISJSn5xMX6v2azZDQkIOK5RBQUFERkYet2iKSMO0Zs0aRowYwdatW1m9ejXvv/8+hmFoEpCISCOg4QQRqTEFBQVMnDiRnj17snXrVqKjo7nlllusjiUiIrWoVh5X+fnnn9fGYUSkDvn444+5/fbbKy+rufnmm5k1axaRkZEWJxMRkdqkZ6OLSLV78803ufrqqwFo1qwZixYt4qKLLrI4lYiIWEFlU0Sq3ZAhQzjrrLMYNGgQjz/+OCEhIVZHEhERi6hsishp27dvH88//zyPP/44drudgIAAvvvuOwIDA62OJiIiFlPZFJFT5vV6SUlJ4Z577sHtdhMbG8ukSZMAVDRFRARQ2RSRU7Rjxw5GjRpVOQGwW7duDBgwwNpQIiJS5+jWRyJyUsrLy3nqqafo0KEDn3/+OQEBATzzzDOsXr2aDh06WB1PRETqGI1sishJueOOO0hJSQHgwgsvZNGiRbRq1criVCIiUldpZFNETspdd91FTEwMS5YsYcWKFSqaIiJyTBrZFJFjWrVqFRs2bGD8+PEAdOjQgd27d+Pv729xMhERqQ80sikiR+R2uxk3bhx9+vRh0qRJbNiwofIzFU0RETlRGtkUkcO8//77jBkzhr179wLw5z//mRYtWlgbSkRE6iWVTRGplJmZycSJE1m2bBkArVq1YtGiRVx44YUWJxMRkfpKZVNEACgrK6N79+7s3LkTm83G5MmTmTZtmm7OLiIip0XXbIoIAL6+vkyaNIkOHTqwZs0annrqKRVNERE5bYZpmqbVIQ7ldrtxOp24XC5CQ0OtjiPSYHm9XubPn09SUhIXXHBB5bLy8nL8/PwsTiciInXZyfQ1nUYXaYR++uknRo0axapVq2jVqhWbN28mMDAQm82moikiItVKp9FFGpGysjKmT59Op06dWLVqFcHBwUyaNEm3MhIRkRqjkU2RRmL9+vWMGDGCTZs2ATB48GAWLFhA8+bNLU4mIiINmcqmSCOwefNmzj33XLxeLxEREcyePZsbb7wRwzCsjiYiIg2cyqZII9C+fXuGDBlCUFAQs2fPJiYmxupIIiLSSOiaTZEGKC8vj0mTJpGVlQWAYRj885//ZNmyZSqaIiJSqzSyKdLAvPvuu4wdO5a0tDSys7N5+eWXAXA4HBYnExGRxkgjmyINRHp6OldffTWXX345aWlptGnThpEjR1odS0REGjmVTZF6zjRNli5dSlJSEm+++SZ2u52pU6eyadMm+vbta3U8ERFp5HQaXaSemz17NpMnTwbgnHPOISUlhS5dulicSkRE5CCNbIrUc3/+859p1aoVTzzxBN98842KpoiI1Cka2RSpZ7Zs2cLSpUt58sknMQyD8PBwfvzxRz1mUkRE6iSNbIrUEx6Ph4cffphzzjmHp59+mldeeaXyMxVNERGpqzSyKVIPrFmzhhEjRrB161YALr30Ui688EKLU4mIiByfRjZF6rCCggImTpxIz5492bp1K9HR0bz++uu8++67NGnSxOp4IiIix6WRTZE67KqrruKjjz4C4Oabb2bWrFlERkZanEpEROTEqWyK1GH3338/27ZtY8GCBVx00UVWxxERETlpKpsidYRpmrz11lvk5uZy++23A9CnTx+2b9+Or6+vxelEREROjcqmSB2wf/9+xo4dy3vvvUdAQAADBgygVatWACqaIiJSr2mCkIiFvF4vixYtIikpiffeew8fHx/uueceTf4REZEGQyObIhbZsWMHo0aN4vPPPwfg3HPPJSUlhQ4dOlgbTEREpBqpbIpYwOVykZycjMvlIjAwkMcee4y77roLu91udTQREZFqpbIpYgGn08nEiRNZtWoVixYtqrw+U0REpKExTNM0rQ5xKLfbjdPpxOVyERoaanUckWpRUlLCo48+yhVXXEHXrl0BKC8vx263YxiGxelEREROzsn0NY1sitSwlStXMmrUKLZt28b777/P+vXrsdvt+Pjof34iItLwaTa6SA1xu92MGzeOvn37sm3bNuLi4njooYd0XaaIiDQqGloRqQHvv/8+Y8aMYe/evQCMHDmSp556irCwMGuDiYiI1DKVTZFq9vHHHzNkyBAAWrVqxeLFi7ngggssTiUiImINlU2RajZgwADOP/98unXrxrRp0wgMDLQ6koiIiGV0zabIaUpNTWXUqFEUFBQAYLPZ+O9//8tTTz2loikiIo2eRjZFTpHX62XBggVMmTKFgoICQkNDeeaZZwA001xEROQ3+hdR5BT89NNPjBo1ilWrVgHQs2dPRo4caXEqERGRuken0UVOQllZGY8//jidOnVi1apVBAcH8/zzz7Ny5UqSkpKsjiciIlLnaGRT5CTcd999zJo1C4CLL76YBQsW0KxZM4tTiYiI1F0a2RQ5CX/5y18444wzeOWVV3j//fdVNEVERI5DI5six/Dpp5/y8ccfM3PmTAASEhLYtm2bngIkIiJyglQ2RY4gLy+Pe+65h5SUFAD69evH4MGDAVQ0RUREToLKpsgfvPvuu4wdO5a0tDQAxo4dS8+ePS1OJSIiUj+pbIr85sCBA4wfP5633noLgDZt2pCSkkKfPn0sTiYiIlJ/qWyKAKZpMnDgQLZs2YLdbufee+/lb3/7G/7+/lZHExERqdc0G10EMAyDxx9/nHPOOYd169Yxffp0FU0REZFqYJimaVod4lButxun04nL5SI0NNTqONJAVVRUMGfOHGJjY7n++uurLNcEIBERkWM7mb6m0+jS6GzZsoWRI0fyzTffEB4ezsCBA4mOjgY001xERKS66TS6NBoej4eHHnqIc845h2+++YbQ0FBmzpxJZGSk1dFEREQaLI1sSqOwZs0aRowYwdatWwG49NJLmTdvHk2aNLE4mYiISMOmsikN3u7du+nduzcVFRXExMTw/PPPM3z4cAzDsDqaiIhIg6eyKQ1e8+bNGTNmDPn5+cyaNUunzUVERGpRtV+zOWPGDLp160ZISAgxMTEMGzaMbdu2VfdhRI4qOzubkSNH8vPPP1cumzNnDkuXLlXRFBERqWXVXja/+OILxo0bx5o1a1ixYgXl5eUMGjSIwsLC6j6USBWmafLPf/6Ttm3bsmTJEsaMGVP5mc2muXAiIiJWqPbT6P/5z3+qvH/xxReJiYlh/fr19O3bt7oPJwLAvn37GDduHO+99x4Abdu25bHHHrM4lYiIiNT4NZsulwuAiIiII37u8XjweDyV791ud01HkgbE6/WSkpLCPffcg9vtxtfXl7/+9a9MnToVh8NhdTwREZFGr0bLpmmaTJ48md69e9O+ffsjrjNjxgymTZtWkzGkAVu6dCmjR48G4NxzzyUlJYUOHTpYnEpERER+V6OPqxw3bhzvv/8+q1atomnTpkdc50gjm4mJiXpcpZyQ0tJS+vbtyzXXXMNdd92lJwCJiIjUgjrxuMrx48ezfPlyvvzyy6MWTQCHw6HTnXLCNm7cyKxZs1iyZAm+vr74+fnx9ddfawKQiIhIHVXt/0Kbpsmdd97Jv/71Lz799FNatmxZ3YeQRqikpIS//vWvJCcn88orrzBr1qzKz1Q0RURE6q5qH9kcN24cy5Yt47333iMkJIQDBw4A4HQ6CQgIqO7DSSOwcuVKRo4cyfbt2wG46qqruOWWWyxOJSIiIiei2q/ZPNojAF988UX+/Oc/H3f7k7kGQBo2t9vN1KlTmTdvHgDx8fHMnTuXyy+/3OJkIiIijZul12zW4HwjaWRGjx7N66+/DsDIkSN56qmnCAsLszaUiIiInBQ9G13qrGnTprF582bmzJnDBRdcYHUcEREROQUqm1InmKbJsmXL2L59e+V9V9u0acPmzZuPemmGiIiI1H0qm2K51NRUxowZwwcffIBhGAwdOpTk5GTg6NcAi4iISP2ge8aIZbxeL3PnzqVt27Z88MEH+Pn58eijj9KpUyero4mIiNQrpmny2baMOjl3RiObYomffvqJkSNH8tVXXwHQq1cvFi9eTFJSksXJRERE6pd9v/7Imnfn8peMwcy+pgvDujSxOlIVKptS6zweD/379+fAgQMEBwczc+ZM7rjjDt2cXURE5CRUlJez7o3H6bh9LlcaHj63OZn2bz8GtI0l2FF3Kl7dSSKNhsPh4PHHH+ett95iwYIFNGvWzOpIIiIi9crOH76h7J07Oa98O/w2vWGc40OuuPquOlU0oQZu6n66dFP3hqeoqIiHHnqI888/nyFDhgD/ux+rJgCJiIgcn9drkuPKI2/PD2St+xddU1/C16g4+JlpsC76Ctrd/AzBoeG1ksfSm7qLHOrTTz9l1KhR/Prrr7z22mv8/PPPBAQEqGSKiIgchcdTzI+r3uOtrOZsyzVJc5WQ4fbQwfsjbzumcSZUjmbutjWlePBsup870MrIx6SyKTUiLy+Pe+65h5SUFACaNm3K/PnzCQgIsDiZiIhI3VPqKeHHr/+NZ+PbnO36ks4UklI6nnXeHpXrpBsRlf93mWnn28Rb6HLDY/gHBFkR+YSpbEq1e+eddxg3bhxpaWkA3HHHHcycOVOXRYiIiBwiJ2Mfv37zf5g7PqGNaxWdKKzy+SX2Nfyftwdhgb7EhfrTJDSMtQVDqAhpSkKPa+hx9jkWJT85KptSrdatW8cVV1wBHHwCUEpKCn369LE4lYiIiPVKyir4fq+LnFV/54zdr3FG2S8kG4dPnSkwA/jJ2Zu4Tleztc9FBPodWtd61V7gaqKyKdWqW7du3HTTTSQmJvLggw/i7+9vdSQREZFaV1JUwO6ta8n7ZR3LjfP5Lq2M7en5VHhN7rD/zGDfHZXXXQIUmv78GNoLe4crSOpzOcl1/NT4yVDZlNOyc+dO7rnnHp5//nni4+MBWLp0qSYAiYhIg1dW6iHnwB6yUrdTsP9HzKyfCXD9QpRnD3HeDM76bdTyKY8PP5pnV273pbcjU3idX+wtyYjpRUj7wbTueiHJ/oFWfZUapbIpp6SiooI5c+bwwAMPUFRUhK+vL6+99hqg2xmJiEj9V1yYT1baLlzpuyjJTiWtPIR19i4ccJVwwF3CAVcJK0pvItYoIvZIOzjkn8IOtp1sMJNoHRNM+yZOzmvZgayEIZyR0JwzausLWUhlU07ali1bGDlyJN988w0A/fr145FHHrE4lYiIyMkxTZOMtFTSf15P0b4fMLK2EZr/C3Flewgnn0Qg8bd1V1R05eWyiCrbZ/qF4TSKDttvoenPft9E8oLPxBvXiavaDeLeszoT4Gev+S9VB6lsygnzeDxMnz6dGTNmUFZWRmhoKE8//TQjRozQoyZFRKTOKnDnkvbrFlx7f2KlX19+zS5iZ1YBOzMLuce7hD/7fHzcfcQaOZX/t2FATIiDn2znkG9vhScoASOqNUFNkoht2YGouGa01r+LlVQ25YQ9/fTTlSOYl112GXPnzqVJkyYWpxIRETl4n8r01J/JSf2J4rSfMLJ3EJy/k+jSVGLIofVv691RMo9Mwiq322E//N+xdCLJ9kugyBFDeVAchCZgjzqTt8/qSUKYP9HBDnzsNmBArXy3+k5lU07YhAkTWL58OXfffTdXXXWVrs0UEZFa5c7LJmP3j+wpcvCTJ5w92UXsySnCk7mLf3rGkGiYlae9j6alkUamGYbNgMSIQIzQ81htlmKPTcLZrD0JZ3Yi1hlx5Osw5ZSobMpRffzxx7zyyissXboUm81GcHAwa9asUckUEZEa4SkpImPvL+Tt+5mijJ14c3fhl7+XkOK9RJenEU4+ocCH5cN4pvzqyu18CMB0GMDh96zMJZQDvk3JD2qBN+JM7jy7HwktzqJZRBB+Pr+f6r62Vr5fY6WyKYfJyclh8uTJLF26FIALLriAW2+9FdBMcxEROXWVZXL/r7jyslnj6MXe3CL25hazN7eYmSXT6GfbdNzRyeZGepX3QQEBfOt7Lj6OAMpCErFHtya0SRJxrToQHhVHeM19JTkBKptSyTRN3nzzTcaPH09GRgaGYXDXXXcxfPhwq6OJiEgdZ5om7uJyDqTtpXjnGjzZuzFz9+BbsI/gkjQiy9OJIq9yhrfLDORGT0qVfez1iYIjzKvxmgaZRgTZfgkUBjYjMLYbL7TrQvOIIJpFBOIM9AUG1cbXlFOgsikA7Nu3j7Fjx7J8+XIA2rZtS0pKCj169LA4mYiI1AUlRQVk7t+J68BOirP2UJ6bij1/H+/6DWVdURxprhIKPOVcYPuOv/s9fdz9OY0iQinEzcEn5UQF+5Hh35n1VFAa0gxbRHMColsS3uRMopueSWxAkK6jrKdUNgWA66+/ni+//BJfX1/++te/MnXqVBwOh9WxRESkFnjKyjjgLiXNVUKaq5g0Vwmdt8wkoGgfIZ50IisyCcdd5b6Tv3u1tDk/e4Mr3+83o454jCzCyPaJpSAggdLgJtjCmvFCu2QSYqJpEhb42z0oB9bYdxTrqGwKAM888wwTJkxg4cKFtG/f3uo4IiJSTQrcOWTv34U7I5Xi7D1UuPZhy0/DUXSA4NJMIioy+bKiPRPL7qyy3VeOT2liZB93/wlGFg4fGwlhATQJC6BlSASrC2/HHp5IYFRznPGtiG7SiqiAII5cQ6WhU9lshMrLy5k1axYA9957LwDJycmsWrVKE4BEROqJivJycjL2kntgN4VZqXhy9/Nl6FAO5HvIcHs44C7hJtdCbjHeJ/g4+4o/5Iblv0szI2liZFNu2sgyIsn1jabIP5ayoARwNsUR2YyQ2BaMTjybKeGRf/j3Q5dgyf+obDYyGzduZMSIEXz33Xf4+vpy1VVX0apVK0AzzUVE6gLTNHEVl5Hu9pDuLsGVsYfoHW9jFKThV5ROUGkmYeVZRJq5RBsm0Ydse3tJi8prIAHS7cHge/RjlZl2so0IKgKiuOTseOJD/Ylz+pMQFkCAsYiM8DAiYxOJ8/Ehrua+sjRwKpuNRElJCY888ghPPvkkFRUVhIeH8+yzz9KyZUuro4mINAqm14vblUNu+h4KMlMpztlHhSsNCg7gV5ROgCeT0PJsni6/hnfKzqvcrrWxlxWOuYfv8AjjAzFGLm7zYNkM8LXjCjqDH7ydKPaPoSwoHiM0HkdkM4KiEomIa0F4TJPKItnzsL3FV9dXl0ZOZbMRWLlyJSNHjmT79u0ADB8+nDlz5hAXp/9OFRE5XabXS152OnkZeyjI2kdJ7n6yyvz4xtGLjPwSMtwe0vNLWFI4njbGXpzH2V+sN6PK+3QzrMp7r2mQYzjJtUdR6IjGExCDNzgeu7MJ01r3JzomnthQf0L9fTCMwcCUav2+IidLZbOBy87O5qKLLqK4uJj4+HjmzZvHsGHDrI4lIlLnlZeVklVYQWZB6cHSmO/Bf/cXJBz4L37FmQSWZuEszyHCzCHcqKhy4/D13ta8VFr1mdsuv8AjjkYeKo9gYoN96RMTRWyoP7GhDmJDHGwomEtgRAJhcc2JjE0kytdPk22k3lDZbOAiIyOZNm0a27dv56mnniIsLMzqSCIiljFNk4KiYvL2/4w78+AoZJnrAGZBOj6FGTg8WQSXZuH05hJuurnYM59cQiu3v93+FZf7vlt1p0c6nU1elffOAF92+ybhbwRR4h9NeWAcRkgsvmHxBEUmEhrTjIi4RMICgrgVuPWwPeqSJ6m/VDYbmMzMTCZOnMjYsWPp1asXAHfffbcm/4hIg1bqKSEnYy/uzL0UZe/Hk5eGNz8dW2E6X/v1YmVFOzILPGTme0go38unjruPv1MDYow8cs3/lc0Ms+qDD3MIJc8WQYFfFB7/aMoDYzBC4vCNaMbbZ/cgJsSf6BAH/r529IQbaaxUNhsI0zRZtmwZEyZMIDs7m02bNvH9999js9lUNEWkXiovKyUvM428zL3kFXnY6duarIJSMvM9ZBZ4uHLP47T0/ESYNwcnhcTBEWdMf1nmx3cV/zulnUHYMY9bavqQY4Tj8omgV5NwukY1IybEQUyIPwm+rdjOYJwxTYmIaUqEn4OIav3WIg2PymYDsGfPHsaMGcOHH34IQMeOHVmyZAk22xEeMCsiYqEKr0l2oYes/FKyCjyU791A4L5VUJCBT3EmDk82wWU5OL25hJn5RBkmUcBGbyvuKX2syr5u8NtLc1vqcY8ZbeQBEB7oS3SIg+jgCL7JvwgzIAJC4vAJjcM/PJ7gqCaExyQSGh5NnM1GHPC3I+4x6TR/BZHGRWWzHvN6vcybN4+pU6dSUFCAn58ff/vb37j33nvx9T3GjdVERKpReVkZuZn7cGXuoygnDY/rABXudCjMPKRA5hLqzeViz/Qqp6JH2D/gQd9Xj7zjQ07KRBuuwz7ONJ0UmQ5ybWHk2yMockRRFhCNNygGe0gsfmHxBEbEc3F8K66Pboqfz6H/Aa6bjovUFpXNeuzf//4348ePB6BXr16kpKRw9tlnW5xKRBoCT0kRuRn7cGftpyg3jTJXOuX56RiFmXzrdy5feTuQVeAhq8BDcPFeVjomVbm5+NFE4yLjkHnbmebhNwLymL7kGGHk+4RT6BdFqX8U5UHxPNK2HdHBDqJCHEQHO4gOvIDAwAACq/F7i0j1U9msxy699FKuvPJKLrjgAsaMGaPT5iJyVKbXS2GBC1fmPvKz95NXVMYO//ZkFXjILjh4SvuavdM5w7MVp5lHKEVHvQZyTXk5q8sTKt+XHPfOkQcLZK4RRvsYX+LCY4gKdhAd4qCZTwTrPWcQEJ5AcGQ8YTGJhISGE2+zHXZL8V6n9QuIiFVUNuuRb7/9lgcffJDXX38dp9OJYRi89dZbVscSEYuVV3jJyD/4LOzCPZvw3bMS3Gn4FqYRUJJBWHkm4d5cgo3Symdkb/G24IHS6VX2c4NvGon2/cc9XhQHT2n7+9qICnYQFRzGuqLzMR2hVARGYwuOwSc0joDwuCoFMs5m44kj7lGntEUaMpXNeqCoqIi//e1vPPvss3i9Xh5++GGeffZZq2OJSC0oLswnK20X7vQ9FGfvoSxvPzb3PvyKDhBUmskE7mVbYSBe8+D6t9o/5CHfVw7f0R9uShF1hGsgs3CSbwaQZwunwCecYr8Iyvwj8QbFYAuOwdcZR2B4LOfHnMGW+OYE+dkPudvF8ur94iLSYKhs1nGffvopo0aN4tdffwXg+uuv569//avFqUTkdJUUFZCTsY/8rH1keENINWNJd5eQ7i6hMCeNyWn3EOHNwkkhicfYj68nDa95RuX7A+bhN+JxEUSuLYJCnwhKHJGUBUThDY7nyaSORIX4ERl08DrIyICB+Dv8CKmB7ysijZfKZh2Vl5fH3XffzZIlSwBo2rQpCxYs4JJLLrE4mYgcien14i6pILPAQ3aBh+zCUoJ2/Bu/nG3YirLwK8nCvzSH4PI8wry5hBjF/H7V40flw5hVfnXlvhyU8rz/7uMe02satAosxgh3EhvqT7zTn1aOML4tjScgshmhsc2Jim+OMyjkBK6qFBGpGSqbddT9999fWTTHjh3LjBkzCA0NPc5WIlKdSsoqyMnOoDDtZ4pyD+BxpePNz4DCTOwl2Tg8OQSV5RBSkUeWGcKfPDOqbP+S76v0sG86fMd/OKUdTdVT2h78SDMjMG2+uHyiKPKPpSwoHkLj8QtvSmBUImGxzYmMa8ZsP8cRkmsqjYjUHSqbddTDDz/Mpk2bmDFjBn369LE6jkiDUF5WSl52OvnZ+ynMOUCJO4NydwZmQSb24ix8S7KZFzCG7SWhZBeUUuAp53r7J0z3XXLcfdsoP2xZ1hHGE/PNAFy2MPJ9Iih2RFLmH0V4VDIzWnUgOthBnPPgCGVE4C8YNhsJh+1BRKR+UdmsA0zT5KWXXmLVqlWVo5nR0dGsWrXK4mQidZvp9eJ25eDO2kdBzgGK8zIoc6ez26cVP9jakFVYSnaBhzJ3FikFY3Ee8kSao8nOHcRus83/3pvHPqNQZtrJNZzk+UTQJzGS6BB/IoP9iAx24CwbyybbTQSGxxESlUBYVDwhAUG6JlJEGhWVTYv9+uuvjB49mv/+978ADB8+nMGDB1ucSsQ6RQUu8rIO4Mp3s9+n2cF7QBZ6yCkoJXnXIhIKvieoLJeQijzCTBdOo+Kw8cPvyoeytPy6yvd2TML93Yedvj6SSMON09+XyCA/ooIdxPt14pvCYXgDo7EFR+MbGnPw0YYRcTijmhAaFkmMzUYMcPgc8DMOP4CISCOjsmmRiooK5syZwwMPPEBRURH+/v488sgjDBgwwOpoItWqpKyCnMJSsgtKyS70ULZ/CyH7vsQszMSnOBs/Tw4BZbmElOfiNN0EGh4CgWJvPCNKn6myr/N8N9LRvuF/C45SHiNxV3nv6+vLNqMVht2HYt9wPI5IKgKiMIKjsAdH4x8WR2B4HKFRCcyNSsDPz+8Pe7z89H8IEZFGSmXTAlu2bGHEiBGsXbsWgH79+rF48WLOPPNMi5OJHF9puZfc/EIK9v5AUW4anrx0ygsyf7vuMRtfTw4BpbkEV+QS5nVxWekj/GI2qdz+GvtnPOG7+Mg7P6Q8Rhj5h32c/dujDctMO3lGKPk2J4W+EXgcEZQHRGEGRuETEkOL6CTead6NqGAHkcF+BPr5ABdX588gIiInSGWzllVUVHDllVeyfft2QkNDefrppxkxYoQeNSmWKSv14Mo6gDvnAIU5B/C40ynPP1gebcXZ+JVk84VPL97z9iK7wIO7pJxw3GzwH3P8nRsQQT6/HLLoSNdAek2DPCMEty2MQp8wSn67mfiks1oTGeIgMsiPiCA/Ynw64AoOIDQskmib7YSexS0iItZS2axldrudF154gblz5zJ37lyaNGly/I1ETsLBGdcHyM9OozA3nRJXOuXuTNxlNj4PuYTsAs/B09qFpTya/zd6sYkoOOakmXXlYews71z5Po9gKkwDu2EedRsXQbiMMNrGBxERHktksIOoID+a+ESzvjgBhzOGoIh4QiPjCYuMI8LHhz/ejvy8w/YaeVK/hYiIWE9ls4YVFBRw//33k5SUxJgxB0eCBg4cyMCBAy1OJvXFH8vjrsD2ZBUbZP9WGBPTP6Nf9usEVeQR6nURRsERy+Mubyy3l3aosqzA1xfsx88QabgI8fchMujgLOuIID9Wu4djdwRiBEXjExKNwxlDYFgszqgmOKPicPo5cALTjrjHc0/ptxARkfpHZbMGffzxx9x+++3s3r2bkJAQrrnmGsLDw62OJRbzVlSQW1BEVjFkFXjIKvCQn5fNGTtexF6Uha8nB//S3KOWx0meJ/nZbFq5v6vse7nD94fjHjfCqDppJtjhwwHflmyjgGK/cEr9wqkIjMIIOnjdo19ozMFb9kTGMSQqniv9A/+wx+TT+yFERKRRUNmsAdnZ2UyePJmXX34ZgBYtWrBw4UIVzQbONE3cednkpO0kP2MPJdmpVLj2Yy/Yj6M4neDSTJwVuYSZbpaUD2dexWWV24ZSyPf+fz+h40Qabn4+5Ox11m/XQBaYAbhsoRTYww8rj/bgaByhMfxf615E/DYy6e9rBy6qzp9ARETkMCqb1cg0Td58803Gjx9PRkYGhmFw11138dhjjxEcHGx1PDkNFeXl5GbsIzd9NwWZeyjN3UtuiZeP/f9EmquEdHcJaa4SXjUeoKvt52PvzIAoo+rjCd0EUmra8TMqKpcVmAHk2ZwU2sP+Vx4DoriiRTJXR7cm4rf7QEb496YkcBLBAUHo/8tERKSuUdmsRjt27OC6667D6/XStm1bUlJS6NGjh9Wx5ChKPSXkZaWRWRFMdglkF3rILijFsf8bzt7/Lg5PNoFlOYRW5BJuuogyKqpcB7nHG82Y0k5V9nnA9+ij178/acZtDyc4rAlXNWlKVLCDqGA/IoP9+LHgJYJCIwiJjMcZGXfU8qj/jxIRkfpEZbMatW7dmqlTp+Lj48PUqVNxOBxWR2pUKrwmuUUHbx6e43Ljt+MDytyZUJiFUZyFX0k2/qU5BJfn4TRdhFJIDHCLZwY/ms0r9zPM9j03+f2n6s6PcPPwWCMXMAGDID87cU5/0u3nsc6MoDw4HpszAUdEIsHRiYTHtSA8Kp4Yu50Y4Ezg6sP22PSwJSIiIvWdyuZp+Pnnnxk3bhzPPvss7dq1A+Cxxx6zOFXD4fWauHMzcWXupTDnACWudMrcGXgLs7AVZeFbko3jt/L4ttmf50sGY/52LWMgJWz1v/uEjhNpuA92xt9kHfLww3LTVjkame8XgycwFm9wPPawpgRENmVF697EhQcR4u/72xb9qufLi4iINBAqm6egvLycWbNm8dBDD1FSUsKECRMqn20uR2d6vbjzsnFl7Tt4/8e8A5S6M/EWZGIrysJV7stLgX/+7bGGpeQWlfKy/TF62Y8/0zqyPK2yaAIU4aDY9CPAKD1s3XwzAJctjAJ7GMV+EfQ8szmdo88kMvjgjcNjHe3ZxYCDt/CJiCHabtfNw0VERE6RyuZJ2rhxIyNGjOC7774DDt4zc+HChRansobp9ZLvzsWdlUZBThrFeQdHHjeG9iPN4zh4OruwlPY5HzGi6O+EmW6cRsUh44ZVpZkRfJ02tMqybPvhT5v5oyLTQai/nU7hYUQF+f1WGh1867oX/4AgHKExBEbEExIZR1hUPCH+gYQcsn2XI+61+RGXioiIyMlR2TxBxcXFPPLIIzz11FNUVFQQHh7Os88+y80334xhHOGCvnrI9HopyM/DnZ1GlseHDDOMnEIPWQWlFLpyuODXJ/ArzSWoLJeQijzCTDehRjl/rIOPeHzZYraqfB9u8xDtl3vE6x4PFYEbMHH42A/Osg7yI8fbjXXeUCoCoiAwEltIDA5nNAFhv5fHBAKDQrgcuPywPZ592r+JiIiInB6VzRP08ssvM3PmTACGDx/OnDlziIuLszjVsZmmSWFpBTkFpZUzrf1/+Q/2nO0YRdnYS7JxeHIIPKQ8hhhlhABvlV/Os+XDK/floJR7/f9wqcBRymOkkV/lGsgcWwQZROC2h1PkG0apXwTlAZGYgQfv/+jnjMU/LIaQiAS2JJxJkMPnkALfu1p/ExEREaldKpsnaMSIEXzwwQfceuutDBs2zLIcRaXl5GalU3jgZ4pyD1DqzqAiPxOjMBN7SQ4OTw4Bv5XHdNPJZZ5Hq2z/su9iets3H3nnh5THKKo+bcaDH/lmACFGMWWmnTwjFLc9jCKfMDy/lUcCo7EFR3Fb837cFdOCiCAHkcF+hDj+hGH8hZjq/jFERESkzlPZPIr333+f2bNn8+9//xt/f398fHx47733qv04RQUu8rIOUJCTVlkevfmZUJSNvTgLP08OswPuZHtxKNmFHkrKvNxoX8Fjvi8ed992yg5blnWEKybLTRt5Rij5tjAKfQ+Wx+Dwc/lL8zZEBPsRGXTwXpC55Z/gjYwl1BlBtM2mSTMiIiJyXDVWNufNm8dTTz1FWloa7dq1Y/bs2fTp06emDldtMjMzmTBhAq+99hoAzz33HFOmTDnh7YsKXORlplGQe4DivAN4XAfL48+OtmwwksgpPDhpxpufwT9LxhBoePjjE6f/KNe1m33mmZXvs81jT5qpMA1yDScuewQ9WkYQEeIgMuhgaQwsu53vjKvxD4slKDwOZ2QcoeHRRNntVW5Y3vWIe444wV9BRERE5KAaKZtvvPEGEydOZN68efTq1YuFCxdy8cUXs3XrVpo1a1YThzxtpmnyj3/8g4kTJ5KdnY3NZmPy5MmMGnkrabu3k5+dhruohF/92/42y9pDdmEpF6U+R6vizQSX5xFmugg0So9YHteXX8q/yv9XEn2xE+jvOaFsUTY3UQEHb8sTGeQgwbcz3xReiTcwCltwFL6hMfg7YwkKjyU0Mh5nRExleXztsL21PsVfSEREROTkGaZ56N0Jq0f37t0555xzmD9/fuWypKQkhg0bxowZM465rdvtxul04nK5CA09/m1vTtePaW52bPiSn96fj70og4QgL82CvTTxLyHSzCXYKK5cd7u3CYNKn6qy/Yu+T9Dfvum4x3mtvD9Ty0cBYDMgPNCPV7gf0+7A4xdOmX/kwfIYdLA8OpyxBEfEVZZHm91evV9cRERE5BSdTF+r9pHN0tJS1q9fz3333Vdl+aBBg/j6668PW9/j8eDx/G+Ez+12H7ZOTfrLPzdxTsZ/eSx61eEf/mG2daRxeLZsnFSYBi4jBLctjEKfsMPKo09IDG1j27CiaUcigvwIC/TDbjOAgTXzpURERETqiGovm1lZWVRUVBAbG1tleWxsLAcOHDhs/RkzZjBt2rTqjnHC4p3+HEg//FrEItNBti2SfJ8ISvwiKA2IxBsYw6NJ7YgM/u0ayGA/Ivx6YwQHEeHjoysaRURERP6gxiYI/fFG56ZpHvHm51OnTmXy5MmV791uN4mJiTUV6zBDOyWQHTOAdSVO/CMTCYlpRkR8S0JCw0m02Q5bv2etJRMRERGp/6q9bEZFRWG32w8bxczIyDhstBPA4XDgcDiqO8YJG9alCXRpgmqkiIiISPU7fOjuNPn5+dG1a1dWrFhRZfmKFSvo2VOFTkRERKQxqZHT6JMnT+amm24iOTmZHj16sGjRIvbs2cOYMWNq4nAiIiIiUkfVSNm85ppryM7O5pFHHiEtLY327dvzwQcf0Lx585o4nIiIiIjUUTVyn83TUdv32RQRERGRk3Myfa3ar9kUEREREfmdyqaIiIiI1BiVTRERERGpMSqbIiIiIlJjVDZFREREpMaobIqIiIhIjVHZFBEREZEao7IpIiIiIjVGZVNEREREaozKpoiIiIjUmBp5Nvrp+P3pmW632+IkIiIiInIkv/e0E3nqeZ0rm/n5+QAkJiZanEREREREjiU/Px+n03nMdQzzRCppLfJ6vezfv5+QkBAMw6iVY7rdbhITE0lNTT3uw+Sl7tHfr/7T37D+09+w/tPfsH6r7b+faZrk5+eTkJCAzXbsqzLr3MimzWajadOmlhw7NDRU/wOrx/T3q//0N6z/9Des//Q3rN9q8+93vBHN32mCkIiIiIjUGJVNEREREakxKpuAw+HgoYcewuFwWB1FToH+fvWf/ob1n/6G9Z/+hvVbXf771bkJQiIiIiLScGhkU0RERERqjMqmiIiIiNQYlU0RERERqTEqmyIiIiJSYxp92Zw3bx4tW7bE39+frl27snLlSqsjyQmaMWMG3bp1IyQkhJiYGIYNG8a2bdusjiWnaMaMGRiGwcSJE62OIidh37593HjjjURGRhIYGEjnzp1Zv3691bHkBJWXl/PAAw/QsmVLAgICaNWqFY888gher9fqaHIUX375JUOHDiUhIQHDMHj33XerfG6aJg8//DAJCQkEBATQr18/fvjhB2vC/qZRl8033niDiRMncv/997Nhwwb69OnDxRdfzJ49e6yOJifgiy++YNy4caxZs4YVK1ZQXl7OoEGDKCwstDqanKR169axaNEiOnbsaHUUOQm5ubn06tULX19fPvzwQ7Zu3cozzzxDWFiY1dHkBD3xxBMsWLCAF154gR9//JEnn3ySp556iueff97qaHIUhYWFdOrUiRdeeOGInz/55JPMmjWLF154gXXr1hEXF8fAgQPJz8+v5aT/06hvfdS9e3fOOecc5s+fX7ksKSmJYcOGMWPGDAuTyanIzMwkJiaGL774gr59+1odR05QQUEB55xzDvPmzeOxxx6jc+fOzJ492+pYcgLuu+8+vvrqK50RqseGDBlCbGwsS5YsqVx25ZVXEhgYyCuvvGJhMjkRhmHwzjvvMGzYMODgqGZCQgITJ05kypQpAHg8HmJjY3niiScYPXq0JTkb7chmaWkp69evZ9CgQVWWDxo0iK+//tqiVHI6XC4XABERERYnkZMxbtw4LrnkEgYMGGB1FDlJy5cvJzk5meHDhxMTE0OXLl1YvHix1bHkJPTu3ZtPPvmE7du3A7Bp0yZWrVrFn/70J4uTyanYuXMnBw4cqNJtHA4H559/vqXdxseyI1ssKyuLiooKYmNjqyyPjY3lwIEDFqWSU2WaJpMnT6Z37960b9/e6jhygl5//XW+++471q1bZ3UUOQW//vor8+fPZ/Lkyfz1r39l7dq13HXXXTgcDm6++War48kJmDJlCi6Xi7PPPhu73U5FRQWPP/441113ndXR5BT83l+O1G12795tRSSgEZfN3xmGUeW9aZqHLZO678477+T7779n1apVVkeRE5SamsqECRP4+OOP8ff3tzqOnAKv10tycjLTp08HoEuXLvzwww/Mnz9fZbOeeOONN3j11VdZtmwZ7dq1Y+PGjUycOJGEhARuueUWq+PJKapr3abRls2oqCjsdvtho5gZGRmH/ReB1G3jx49n+fLlfPnllzRt2tTqOHKC1q9fT0ZGBl27dq1cVlFRwZdffskLL7yAx+PBbrdbmFCOJz4+nrZt21ZZlpSUxNtvv21RIjlZ99xzD/fddx/XXnstAB06dGD37t3MmDFDZbMeiouLAw6OcMbHx1cut7rbNNprNv38/OjatSsrVqyosnzFihX07NnTolRyMkzT5M477+Rf//oXn376KS1btrQ6kpyECy+8kM2bN7Nx48bKV3JyMjfccAMbN25U0awHevXqddjtxrZv307z5s0tSiQnq6ioCJutahWw2+269VE91bJlS+Li4qp0m9LSUr744gtLu02jHdkEmDx5MjfddBPJycn06NGDRYsWsWfPHsaMGWN1NDkB48aNY9myZbz33nuEhIRUjlI7nU4CAgIsTifHExISctj1tUFBQURGRuq623pi0qRJ9OzZk+nTp3P11Vezdu1aFi1axKJFi6yOJido6NChPP744zRr1ox27dqxYcMGZs2axW233WZ1NDmKgoICduzYUfl+586dbNy4kYiICJo1a8bEiROZPn06rVu3pnXr1kyfPp3AwECuv/5660KbjdzcuXPN5s2bm35+fuY555xjfvHFF1ZHkhMEHPH14osvWh1NTtH5559vTpgwweoYchL+/e9/m+3btzcdDod59tlnm4sWLbI6kpwEt9ttTpgwwWzWrJnp7+9vtmrVyrz//vtNj8djdTQ5is8+++yI//bdcsstpmmaptfrNR966CEzLi7OdDgcZt++fc3NmzdbmrlR32dTRERERGpWo71mU0RERERqnsqmiIiIiNQYlU0RERERqTEqmyIiIiJSY1Q2RURERKTGqGyKiIiISI1R2RQRERGRGqOyKSIiIiI1RmVTRERERGqMyqaIiIiI1BiVTRERERGpMSqbIiIiIlJj/h/+5P1WPLkmTwAAAABJRU5ErkJggg==",
      "text/plain": [
       "Figure(PyObject <Figure size 800x520 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5.2))\n",
    "ax.plot(w_grid, w_grid, \"k--\", label=L\"45\")\n",
    "ax.plot(w_grid, σ_star, lw=2, label=\"VFI\")\n",
    "ax.plot(w_grid, σ_star_hpi, \"--\", lw=2, label=\"HPI\")\n",
    "ax.legend(fontsize=fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eadb05c1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "julia",
   "notebook_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Julia 1.9.1",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}