{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"# Asset Replacement Model\n",
"\n",
"**Randall Romero Aguilar, PhD**\n",
"\n",
"This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.\n",
"\n",
"Original (Matlab) CompEcon file: **demdp02.m**\n",
"\n",
"Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
"\n",
" !pip install compecon --upgrade\n",
"\n",
"Last updated: 2022-Oct-22\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## About\n",
"\n",
"Profit-maximizing entrepreneur must decide when to replace an aging asset."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from compecon import BasisSpline, DPmodel, qnwnorm, NLP"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Model Parameters\n",
"\n",
"The maximum age of the asset is $A=6$, production function coefficients are $\\alpha =[50.0, -2.5, -2.5]$, net replacement cost is $\\kappa = 40$, long-run mean unit profit is $\\bar{p} = 1$, the unit profit autoregression coefficient is $\\gamma = 0.5$, the standard deviation of unit profit shock is $\\sigma = 0.15$, and the discount factor is $\\delta = 0.9$.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"A = 6\n",
"ɑ = np.array([50, -2.5, -2.5])\n",
"κ = 40\n",
"pbar= 1.0\n",
"ɣ = 0.5\n",
"σ = 0.15\n",
"δ = 0.9"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### State space\n",
"The state variables are the current unit profit $p$ (continuous) and the age of asset $a\\in\\{1,2,\\dots,A\\}$. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Discrete states are : ['a=1', 'a=2', 'a=3', 'a=4', 'a=5', 'a=6']\n"
]
}
],
"source": [
"dstates = [f'a={age+1:d}' for age in range(A)]\n",
"print('Discrete states are :', dstates)"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"Here, we approximate the unit profit with a cubic spline basis of $n=200$ nodes, over the interval $p\\in[0,2]$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"n = 200\n",
"pmin, pmax = 0.0, 2.0 \n",
"basis = BasisSpline(n, pmin, pmax, labels=['unit profit'])"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Action space\n",
"There is only one choice variable $j$: whether to keep(0) or replace(1) the asset."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"options = ['keep', 'replace']"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Reward Function\n",
"The instant profit depends on the asset age and whether it is replaced. An asset of age $a=A$ must be replaced. \n",
"If the asset is replaced, the profit is $50p$ minus the cost of replacement $\\kappa$; otherwise the profit depends on the age of the asset, $(\\alpha_0 + \\alpha_1 a + \\alpha_2 a^2)p$.\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"def profit(p, x, i, j):\n",
" a = i + 1\n",
" if j or a == A:\n",
" return p * 50 - κ\n",
" else:\n",
" return p * (ɑ[0] + ɑ[1]*a + ɑ[2]*a**2 )"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### State Transition Function\n",
"\n",
"The unit profit $p$ follows a Markov Process \n",
"\\begin{equation}p' = \\bar{p} + \\gamma(p-\\bar{p}) + \\epsilon \\end{equation}\n",
"where $\\epsilon \\sim N(0,\\sigma^2)$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"def transition(p, x, i, j, in_, e):\n",
" return pbar + ɣ * (p - pbar) + e"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"The continuous shock must be discretized. Here we use Gauss-Legendre quadrature to obtain nodes and weights defining a discrete distribution that matches the first 10 moments of the Normal distribution (this is achieved with $m=5$ nodes and weights)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"m = 5\n",
"e, w = qnwnorm(m,0,σ**2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"On the other hand, the age of the asset is deterministic: if it is replaced now it will be new ($a=0$) next period; otherwise its age will be $a+1$ next period if current age is $a$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"h = np.zeros((2, A),int)\n",
"h[0, :-1] = np.arange(1, A)"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"## Model Structure"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"model = DPmodel(basis, profit, transition,\n",
" i=dstates,\n",
" j=options,\n",
" discount=δ, e=e, w=w, h=h)"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### SOLUTION\n",
"The ```solve``` method returns a pandas ```DataFrame```."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving infinite-horizon model collocation equation by Newton's method\n",
"iter change time \n",
"------------------------------\n",
" 0 3.7e+02 0.1566\n",
" 1 4.3e+01 0.3285\n",
" 2 1.2e+01 0.4877\n",
" 3 5.7e-01 0.6446\n",
" 4 4.5e-13 0.8045\n",
"Elapsed Time = 0.80 Seconds\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" unit profit | \n",
" i | \n",
" value | \n",
" resid | \n",
" j* | \n",
" value[keep] | \n",
" value[replace] | \n",
"
\n",
" \n",
" | \n",
" unit profit | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | a=1 | \n",
" 0.000000 | \n",
" 0.000000 | \n",
" 0 | \n",
" 234.982583 | \n",
" -5.684342e-14 | \n",
" keep | \n",
" 234.982583 | \n",
" 206.124123 | \n",
"
\n",
" \n",
" | 0.001001 | \n",
" 0.001001 | \n",
" 0 | \n",
" 235.058576 | \n",
" -6.159590e-07 | \n",
" keep | \n",
" 235.058576 | \n",
" 206.209036 | \n",
"
\n",
" \n",
" | 0.002001 | \n",
" 0.002001 | \n",
" 0 | \n",
" 235.134563 | \n",
" -6.366804e-07 | \n",
" keep | \n",
" 235.134563 | \n",
" 206.293949 | \n",
"
\n",
" \n",
" | 0.003002 | \n",
" 0.003002 | \n",
" 0 | \n",
" 235.210544 | \n",
" -2.317537e-07 | \n",
" keep | \n",
" 235.210544 | \n",
" 206.378862 | \n",
"
\n",
" \n",
" | 0.004002 | \n",
" 0.004002 | \n",
" 0 | \n",
" 235.286520 | \n",
" 4.292312e-07 | \n",
" keep | \n",
" 235.286520 | \n",
" 206.463774 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig1, ax = plt.subplots(figsize=[10,5])\n",
"\n",
"ax.set(title='Action-Contingent Value Functions', xlabel='Net Unit Profit',ylabel='Value')\n",
"\n",
"cc = np.linspace(0.3,0.9,model.dims.ni)\n",
"\n",
"\n",
"\n",
"for a, i in enumerate(dstates[:-1]):\n",
" ax.plot(S.loc[i,'value[keep]'] ,color=plt.cm.Blues(cc[a]), label=f'Keep {i}')\n",
"\n",
"ax.plot(S.loc['a=1','value[replace]'], color=plt.cm.Oranges(0.5),label='Replace')\n",
"\n",
"for a, i in enumerate(dstates[:-1]):\n",
" if pmin < pcrit[a] < pmax:\n",
" ax.plot(pcrit[a], vcrit[a], color='C2', marker='.', markersize=9)\n",
" ax.annotate(f'$p^*_{a+1}$', (pcrit[a], vcrit[a]), xytext=(pcrit[a]-0.2, vcrit[a]+50), arrowprops=dict(width=1))\n",
" print(f'Age {a+1:d} Profit {pcrit[a]:5.2f}')\n",
"\n",
"ax.legend(fontsize=\"xx-small\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Plot Residual"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvYAAAEsCAYAAACorPf0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABi4UlEQVR4nO3deVwU5R8H8M/ucoOgKJeCqSSogIh5a2lK3iaSmuZ9ZFr+ShPv7NK8yvAoyytvTUkzFTPvNPPAI03xxINDRRTkWmBhd35/ECvLLrgsu+yyfN6vl78f+8w8M98dhum7z37nGZEgCAKIiIiIiKhCExs7ACIiIiIiKjsm9kREREREZoCJPRERERGRGWBiT0RERERkBpjYExERERGZASb2RERksjhxGxGR9iyMHQARVV47d+7E9OnT1dolEgmqVq2KRo0aYeDAgejUqVOZ9nPmzBkMHToUHTp0wIoVK1T2PWjQIHz66adl2r6piI+PL9WxioqKgqOjowEj0l1qairCw8PRtm1bvPHGG8p2X19f2NnZ4eLFi0aMLl/Hjh2RkJCgcZm1tTWqVKkCX19fhISE4M033yy3uIYMGYKzZ8/il19+QUBAwAvX79mzJ27duoXDhw/D09OzHCI0rd8jkTlhYk9ERufl5YUmTZooX+fk5CApKQmnT5/GiRMnMHnyZIwePdp4AVYwdnZ2WiX4lpaW5RCNbubPn4+dO3eiVatWxg7lhYKDg2Fra6vSlpGRgZiYGJw8eRInT55ETEwMJk6caKQIiaiyYGJPREbXrFkzzJ8/X6391q1bGDBgAMLDw9GzZ0+4u7sbIbqKp1q1avjmm2+MHUaZFFeCs2/fPojFplVFOn36dI0j3YIgYPPmzZg9ezZWrVqFfv36lcuI+IIFC5CVlVVuo+9EZDpM6+pIRFRI/fr10bNnT+Tl5eHYsWPGDodMgLe3N+rWrWvsMLQiEokwePBgNG7cGHK5HCdPniyX/dasWRPe3t6wtrYul/0RkelgYk9EJs3NzQ0AIJVK1ZZdv34dH374IVq3bg1/f3906dIFS5Ys0biuNpYtWwZfX1+cOXMGv/76K9588000btwY7du3R3h4OORyOTIzMzFnzhy0a9cOTZs2xYABA3Du3Dm1bWVmZuLHH3/EW2+9hWbNmsHf3x+vvvoqJk2ahFu3bhW733379uGtt95CYGAgWrZsiUmTJiEuLk6n91Nat27dwkcffYTWrVujadOmGDNmDO7cuYNp06bB19cX8fHxynV9fX0RFBSkcTs9e/ZUWx8Azp8/jwkTJqB9+/bw9/dH06ZNERoainXr1kEul6ts+9dffwUAfPTRR8pjU9J+o6OjlbH7+/ujU6dO+Oqrr/D06VOV9c6cOQNfX18sW7YMV69exbvvvotmzZohKCgIQ4cOxdmzZ3U7eCXw8PAAkH/fQGEKhQIRERHo168fgoKC0LRpUwwePBiHDh3SuJ1Dhw5h2LBhaNu2LRo3bowuXbpg3rx5SE5OVllvyJAh8PX1xb///qvSfvPmTXz44Ydo06YNgoKCMHbsWNy5c0fjvgz5+yUiw2EpDhGZtILkpHHjxirthw4dwoQJE5Cbmws/Pz80a9YMV69exfLly3Hs2DFs2LABVapU0WmfK1euxF9//YWgoCC0atUKZ86cwY8//oisrCxcuHAB9+/fR9OmTZGUlISLFy9i+PDh2LVrF15++WUA+Un9gAEDcPPmTdSqVQstW7ZETk4Orly5gr179+Lo0aPYu3cvatasqbLf9evX4/Dhw/D19cVrr72GixcvYu/evYiKisLBgwcNOgJ74cIFjB49GpmZmQgICICHhwfOnj2LgQMHon79+mXe/s6dOzFjxgxYWFggKCgIgYGBSEhIwJUrV3D16lXEx8fjk08+AQD06tUL//zzD+Li4tCsWTN4eHigRo0axW77999/x+TJk5Gbm4vAwEC4u7sjOjoaGzZswB9//IH169erjfL/888/WLlyJapXr46WLVsiNjYWZ86cwYULF7Bt2zb4+fmV+T0D+R9ICz74FZwfQH6ZTlhYGCIjI+Hg4ICmTZtCJBIhKioKH3zwAT788EN88MEHyvUPHDiADz/8EFZWVmjWrBns7e1x5coVrFu3DseOHcOuXbvU6vwLu3DhAkaNGgWpVKo8RlFRURg4cCAkEkmZ32dpfr9EZEACEZGR7NixQ/Dx8RGmTp2q0i6TyYQHDx4I3377reDj4yO89957KssfP34sNG3aVAgICBCOHz+ubM/NzRVmzZol+Pj4CDNnzlS2nz59WvDx8RHGjBmjtu8vvvhC2bZ06VLBx8dH8PHxEX7//Xdl+7Fjx5Tt3bp1Ex4/fqxcNmnSJMHHx0f4+uuvlW0//PCD4OPjI0yfPl2Qy+XK9oyMDGHw4MGCj4+PsHz5crX9+vr6Cr/88ovK+r179xZ8fHyEX3/99YXHMy4uTvDx8RFef/31F65bWG5urtClSxfBx8dH2LRpk7L92bNnQr9+/ZTvPS4uTrnMx8dHaNKkicbt9ejRQ2V9qVQqvPLKK0KTJk2EGzduqKx7+PBhwcfHRwgMDBRkMpmyferUqWq/B037ffTokRAQECD4+fkJhw8fVrbn5eUJ8+fPF3x8fISQkBBBoVAIgvD8XPDx8RFmzZol5OTkKPvMnDlT+XvTxuuvv652XARBEBQKhZCeni6cOXNGGDhwoODj4yO8+eabQl5ennKdzZs3Cz4+PkL//v1Vzqe4uDghODhY8PX1FaKiopTtnTp1Eho1aiTExMQo22QymTBy5EjBx8dHiIiIULYXnGOXL19WHovu3bsLPj4+wvbt25XrpaenK9ct799vSdsnIt2xFIeIjO7XX3+Fr6+v8p+/vz86dOiAH3/8EQEBAQgPD1dZ/5dffkFGRgZGjRqFV199VdluYWGBmTNnokaNGti1a5da6YO2WrRoga5duypft2/fHnZ2dgCA8ePHw8XFRbmsc+fOAIDY2Fhlm62tLV577TV89NFHKjd62tvbK6c9fPjwodp+W7Vqhbfeektl/ZCQEABQK6soSUJCgsrx1PRv2bJlyvWjoqJw9+5dtGjRAoMGDVK2Ozk5YeHChRCJRFrvW5MnT56gffv2eO+99+Dj46OyrGPHjnB3d0dWVhZSUlJKve1t27YhJycHgwcPRseOHZXtEokEU6ZMQcOGDREdHY3Tp0+r9LOzs8OMGTNgZWWlbCt476U51gDQqVMnlWPboEEDvPLKKxgyZAjOnz+PVq1aYeXKlSoj4+vWrYNIJMLChQtVzidPT09MnToVgiBgw4YNyvakpCRYWFjA2dlZ2WZpaYmpU6di9uzZxZbNAPklMrdv30br1q3Rr18/ZbuDgwPmzp1r0r9fIiodluIQkdEVne5SLpfj2bNnuHLlCv79918MGjQIK1euVJZjREVFAQBatmypti1ra2s0b94cv//+O/755x+0b9++1PEULfsBgKpVq0IqlaJBgwYq7QXzwOfk5Cjbhg0bhmHDhqmsl5aWhuvXrytjz83NVdtHYGCgWltB0lea+wa0me7S19dX+fP58+cBAK+99praenXq1EGjRo1w9epVrfdflJeXFxYtWqTSplAocP/+fVy6dEl57DQdkxcpKHMp+IBVmEgkQteuXXHt2jWcO3cOrVu3Vi7z9fWFjY2NyvoF51dp79EoPN1lSkoKTp06Bblcjtdffx2TJk1SK2VKTEzE/fv3UbNmTbz00ktq22vTpg3EYrHyXAHyZ47666+/EBoain79+qFDhw5o2LAhfHx81JLpogqOUbt27dSWeXl5wcfHBzdu3CjVey66DUP9fomodJjYE5HRFTfdpVQqxdSpU3HgwAFMnToVa9asAQA8evQIANSS56ISExN1isfJyUmtrWBUs2rVqhrbNe1706ZNiIqKwr1795SjlQXrCxqmc9R0T0DBKK+m9YtT2ukuk5KSAKDY6US9vLzKlNgD+fEfPnwYv/32G27evImEhARlolfSMXmRgthr1aqlcXlBe8F6BTQ9mEuXYw2oT3cZExODESNG4OjRo2jUqJHGxB4AHjx4oPIBq6jk5GTk5ubC0tISs2fPxgcffIDo6GgsXrwYixcvhqurK4KDgzFo0CCV+v2iCt57wY3oRdWqVatMiT1guN8vEZUOE3siMll2dnb4/PPPcfjwYfz1119ITEyEm5ubcoaNzp07l3hDacFsJKVV1psJT506hbFjxyI7Oxuenp5o2bIl6tevj4CAACQnJ2PatGka+5W1JEJXL9qvhUXp/lNRdAYUuVyOsWPH4vjx47C2tkZAQADatGmDBg0aoHnz5vjf//6H27dvlzpu4MXJokKhAACVkhvAsMfa29sbS5cuxTvvvIPvv/8etWvXVpZUAc+Pj4uLywsfwJWXlwdLS0vUrFkTO3fuxJkzZ3Do0CGcOHEC9+7dw5YtW7B9+3Z8++236NKli07xmvLvl4hKh4k9EZm06tWro1q1anjy5AkePnwINzc3uLq64t69exgzZgwCAgKMHaIKQRAwa9YsZGdn45tvvkGvXr1Ulm/fvt1IkRWv4ANQQkKCxuWPHz9WaxOJRMVOYZienq7yevfu3Th+/DiCgoLwww8/oFq1airL09LSdAkbAODq6oq7d+8iPj5e44h0wZSM1atX13kfumjSpAlGjhyJVatWYc6cOWjdurUyvoLyqqpVq5bqmxWRSIRWrVopPwzEx8dj5cqV2LZtW4mJfcF+Hzx4oHF50W8zCvZlCr9fIiod3jxLRCYtNTUVz549A/C8VKRp06YAgOPHj2vsM3LkSAwcOBA3b94slxgLS05ORlxcHDw9PdWSegDKhxQVjCSbgoJ7FTTNn/7s2TNcunRJrd3Ozg45OTnIyMhQaU9ISFBLFAv6DxgwQC3pu3XrlvKDQ+HRd21H1F955RUAwMGDB9WWCYKgbG/evLlW29On8ePHw9PTE+np6ViwYIGy3dPTE+7u7rh3757GZxRcu3YNXbt2xZQpUwDkH9PevXvjvffeU1nP09MTs2bNglgsVpanaVLwQeDw4cNqy54+fYorV66otRv690tEhsHEnohMVk5ODr766ivk5eXhlVdeUSb2/fv3h7W1NVauXKmS3AuCgO+++w4nT55EYmIivL29yz3mqlWrwsbGBo8ePcL169eV7QqFAmvXrsX+/fsBqN5sa2xNmjRBUFAQ/v33X3z//ffKdplMhlmzZmmMteCGzY0bNyrbsrKy8MUXX6itW/B7O3bsmEpyFxcXh48//lj5uvB+Ckpnio4OF9W/f3/Y2Nhg06ZNOHLkiLJdoVBg0aJFuHbtGnx9fZUfBsuTjY0NZsyYAQCIjIzEqVOnlMsGDx6M3NxcTJ48WeUbkZSUFMycORN3795VPuegZs2aePbsGY4fP67yHgFg3759UCgUJX5zFRgYiMaNG+PSpUv48ccfle05OTmYOXOmxptaDf37JSLDYCkOERnduXPnEBYWptKWlZWFqKgopKamomrVqioJRa1atTBnzhxMmzYN7777Lho1agRPT0/cvHkT9+7dg52dHRYvXqyXB++UlkQiwaBBg7BmzRr069cPLVq0gI2NDa5evYqHDx/i5Zdfxu3bt/HkyRODxZCSkqJ2PDXp3bu3crrQhQsXYtCgQVi6dCn279+PunXr4vLly8rjX/CtSYERI0bg4sWLWLx4MY4ePQoXFxecP38eFhYWaN26tUoSGxISgp9++gm///67MtFOTk7GhQsXIBaLUbt2bcTGxiIpKUl5E2idOnUAAIsWLcKJEycwYsQIjVM6enh4YO7cuZgyZQrGjRuHwMBAeHh4IDo6GrGxsXB3d0d4eLjKtKPlqVOnTmjfvj3+/PNPzJ49G7/99hssLS0xcuRInD9/HkePHkXXrl3RuHFjWFtbIyoqCpmZmWjevDnGjRsHIP/biy+//BJjx47FuHHj0LhxY3h4eODhw4e4fPkybG1tMXXq1GJjEIlEmD9/PoYOHYrw8HDs27cPderUwT///IP09HTUr19f7WnIhv79EpFhcMSeiIwuLi4Oe/bsUf7bu3cv/v77b7i5uWHEiBHYvXu32swib775Jn7++Wd07twZjx49wtGjR6FQKBAaGopdu3ZpnLKyvHz88ceYNm0aXnrpJZw/fx7//PMP3NzcMH36dOzatQs1atTAP//8g+TkZIPsXyqVqhzP4v7FxMQo+9SuXRs7duxA37598fTpU/z555+oXbs2NmzYoDLPeoEuXbrg+++/R5MmTXDjxg2cO3cObdu2xS+//KI2u467uzs2b96MTp06ITMzE8ePH0dSUhK6deuGX375BSNGjAAAHD16VNlnwIAB6NmzJ2QyGU6cOFFiWVWPHj2U50JsbCyOHDkCsViM0aNH49dffzXKNzeFzZw5E1ZWVoiJicG6desA5H8A/P777/HZZ5/B29sbly5dwrlz5/DSSy9h+vTpWLNmjcqN4e3bt8eaNWvQrl07xMbG4vDhw3j06BHefPNN7Ny584X3mnh7e2P79u0IDQ1FcnKy8ve7fv161K5dW219Q/9+icgwRAKL3oiIqAQ9e/bErVu3cPjwYZVpHYmIyLRwxJ6IiIiIyAwwsSciIiIiMgNM7ImIiIiIzABr7ImIiIiIzACnu9QThUIBuZyfkYiIiIjIcCwti5/KmYm9nsjlAp49kxo7DCIiIiIyYy4uVYpdxhp7IiIiIiIzwMSeiIiIiMgMMLEnIiIiIjIDTOyJiIiIiMwAE3siIiIiIjPAxJ6IiIiIyAwwsadKg89iIyIiInPGeezJ7CnkAo6tvY7cbDk6vtsQltbFP9iBiIiIqKLiiD2ZvbgryXhyPwOpiVm4efKRscMhIiIiMggm9mT2cqR5z3/OyithTSIiIqKKi4k9mT1B8by2XiQSGTESIiIiIsNhYk9mr/BNsyKe8URERGRECQnx6NKlvUG2zTSHzJ4gf57Yi8UcsSciIiLjOHPmFMaPH4PMzEyDbJ+z4pDZKzzLpYiJPRERUYVy9WEaVp+OhVQmN+h+7KwkGN2qNvw8HF+47v79kdi6dRMePnwAe3t7DBs2EnXreiMs7EO1dTt37obJk2cgKuoMFi78CkOGjMC33y4wxFtgYk/mT1G4xp6JPRERUYWy9UIC/rqTXC77sreSYE6PkhP7hIR4LF78NZYvX4N69bxx8uQJzJo1FXv3HsTBgyeK7efr2wA///wrnjxJ0nfYSkzsyezx5lkiIqKKa2DTWsiUyctlxH7gK54vXM/V1Q3r1v0Md3d3pKSkwMLCAjKZDGlpabCzsy+2n6Ojkz7D1YiJPZk9lcSed5UQERFVKH4ejgjv42/sMJQkEgl27NiG/fsj4eRUFY0a+QEAEhMTMXz4QLX1g4O7IixsWrnExsSezJ6geP4zb54lIiKisjh8+ABOnz6JzZsj4OjohPT0dOzbtweurq7Yv/+YUWPj+CWZPdXpLpnYExERke7S09NhaWkFCwtLSKVSLF++FACQl2f8h2AysSezx5tniYiISF+6desJV1dXhIR0w8CBobCxsUGdOvUQG3vf2KGxFIfMH2vsiYiISF9sbW0xf/63Ovf38KiJv/46p8eInmOaQ2avcI29CByxJyIiIvPExJ4qFQHCi1ciIiIiqoAqTGIfFxeH8ePHo0WLFmjRogWmTJmC5OQXP6xA234nTpzAO++8g8DAQAQFBWH48OH4559/DPBOqNwVHqRnXk9ERERmqkLU2KekpGDYsGGQyWQYPXo05HI51qxZgxs3biAiIgJWVlZl6nf27Fm8++67qF+/PiZOnIi8vDxs2bIFgwcPxpYtW9C4cePyfLtkQIVnyCEiIiIyJxUisV+3bh0ePXqEPXv2wNvbGwAQGBiIESNGYNeuXejfv3+Z+s2dOxceHh7Yvn07bG1tAQAhISHo3r07wsPDsXbt2nJ4l2QohR82W7jenoiIiMicVIhSnMjISLRo0UKZnANAmzZtULduXURGRpapX2pqKq5fv46uXbsqk3oAqFGjBpo3b46LFy8a4B1R+Xqe2XPEnoiIiMyVySf2qampiIuLg5+fn9oyPz8/XLlypUz9HBwcsH//fgwfPlxtvZSUFEgkkrK9ATItzOuJiIjITJl8Yp+YmAgAcHNzU1vm4uKCjIwMpKen69xPIpGgTp06autdv34dFy5cQFBQkD7eBhmRSikOE3siIiIyUyaf2GdmZgKASplMAWtrawCAVCrVW7+CvlOnTgUAjBkzRoeoyVSxFIeIiIiM4fHjREydOhHdunVEnz7dsW7dar3vw+QTe4XixXc7isXqb0PXfllZWRg3bhyuX7+OMWPGoEWLFtoFSqaLN88SERGRkc2b9yU8PGph9+4/sHz5auzfH4kDB/brdR8mPyuOvb09ACAnJ0dtWUFbwTpl7ZeWlob33nsPFy5cwFtvvYWJEyeWLXgyCSrPmuWIPRERUYVikXgRdueWQCTLMOh+BCsHSJt9hDy3F5dh798fia1bN+Hhwwewt7fHsGEjUbeuN8LCPlRbt3Pnbpg0aRqsra0xdOgIWFpawsOjJtq2fQ1Xr15G585d9fYeTD6xr1mzJgAgKSlJbdnjx4/h6OgIOzu7Mvd7+vQpRo0ahWvXruHtt9/GF198AZFIpNaXKjbm9URERBWL7aXVsL53qFz2JVg6IL3zdyWuk5AQj8WLv8by5WtQr543Tp48gVmzpmLv3oM4ePBEsf3mz/9W+XNeXh7Onj2Ffv0G6i12oAIk9o6OjvD09MTVq1fVlkVHR8Pf37/M/TIyMpRJ/fDhwzF9+nT9vQEyusK5PGvsiYiIKpaswNEQ5WaWy4h9VuDoF67n6uqGdet+hru7O1JSUmBhYQGZTIa0tDTY2alXkRQll8sxd+4XsLGxQdeuPfQRupLJJ/YA0LlzZ2zYsAExMTHKOen//vtv3L17F6NGjSpzvy+//BLXrl3D0KFDmdSbI6GYn4mIiMjk5bkFIa3HOmOHoSSRSLBjxzbs3x8JJ6eqaNQof2r1xMREDB+uPgIfHNwVYWHTAOTfy/npp9OQkpKCb79dBisrK73GJhIqwBBmcnIyevbsCYlEgpEjRyInJwerV69G7dq18fPPP8PKygpxcXG4cOECmjZtCi8vL637xcTEoHv37qhSpQpmzJihcd763r17vzDG3Fw5nj3TPMsOGdf5PfcRc/YxAKDBq+5o3NnLyBERERFRRXXw4H5s2PATvv9+FRwdnZCeno5u3V5HRMRueHjULLZfWloaJk78ANWrV8cXX8zTOHOjNlxcqhS7rEKM2Ds7O2PTpk2YN28eli5dChsbGwQHB2PKlCnKTzpRUVGYPn065s2bp0zstel39uxZAEB6enqxo/XaJPZkyp5/djX9j7FERERkytLT02FpaQULC0tIpVIsX74UQH7dfEk++2w63NzcMXv2fIM9ALVCjNhXBByxN13nd99DTFT+TdQ+bd3QpGttI0dEREREFVVWVha++GImLlw4D1tbW3Ts+Aaios5g3Lj/oW3bVzX2uXPnNoYOHQArK2tIJM+nW+/cuRsmT55Rqv2XNGLPxF5PmNibrnO/3cOdc/8l9m3c0KQbE3siIiKqmEpK7E3+AVVE+sSPsURERGSumNiT2Sv8pZSgYGZPRERE5omJPZm/Qrk8R+yJiIjIXBU7K06DBg1K/eRVkUiE6OjoMgdFpE8quTwzeyIiIjJTxSb2ISEhpU7siUwSR+yJiIioEig2sZ8/f355xkFkMCo19szsiYiIyEzppcZeoVAgPT0dBw8e1MfmiPSLI/ZERERUCWj95Nn09HTMmzcPBw4cgFQqLXbk89q1a3oLjkjvmNkTERGRmdJ6xP7rr7/Gzp078fLLL6N169YQBAE9e/ZE69atYWFhAWtrayxbtsyQsRLpRHW6SyMGQkRERJXW7du3MG7cSLzxxmvo1+9N7N27S+/70HrE/tixY3jjjTewbNkypKSkoHXr1hgyZAgaN26Ma9euYdCgQbhz547eAyQqM5VSHNMbsRcEAbJjhyGyd4BVi1bGDoeIiIgMYNasqXjrrf5YvnwN7tyJwZgxwxAY2BReXrX1tg+tE/vk5GS0bdsWAFCtWjW4ubnh8uXLaNy4MRo2bIi+fftiz549GDNmjN6CI9KHwqm8Ceb1yD1zCumfzgAAVNvyCyR6/AMnIiIi/du/PxJbt27Cw4cPYG9vj2HDRqJuXW+EhX2otm7nzt0wefIMrF27BdbW1gCA1NRnsLS0hK2trV7j0jqxt7e3h0LxvI6hdu3auHnzpvK1j48PduzYodfgiPTCxDP7nMMHlD/n/nuZiT0REVEh155FY+PttcjKkxp0P7YWdhjy8gg0rNqoxPUSEuKxePHXWL58DerV88bJkycwa9ZU7N17EAcPnii2n42NDQCgZ8838OxZCkaMeBc1arjo9T1ondg3btwYv//+O95++21IJBK8/PLLOHPmDARBgEgkwt27d2FlZaXX4Ij0QXW6SyMGUqxCz4swzQCJiIiMZsfdbTj9+GS57Mvewh4zm3xe4jqurm5Yt+5nuLu7IyUlBRYWFpDJZEhLS4Odnf0L97Fjx17cvRuDiRPHo1Ejf7Ru3VZP0ZcisR85ciRGjRqFLl26YMeOHejTpw+2bt2KkSNH4qWXXsKOHTvQqVMnvQVGZAgmmTcXfhAc7+4lIiJS8VbdtyGVS8tlxP6tOv1fuJ5EIsGOHduwf38knJyqolEjPwBAYmIihg8fqLZ+cHBXhIVNU762trZGgwaN0LlzV5w+fdI4iX3r1q2xcuVKrF+/Ho6OjmjcuDE+//xzLFiwAKdOnUJgYCCmT5+ut8CI9MbES3Eg5og9ERFRcRpWbYS5zb42dhhKhw8fwOnTJ7F5cwQcHZ2Qnp6Offv2wNXVFfv3H9PYRyrNxLBh72DNmg1wdHQCAMhkuXB2dtZrbFon9gDQrl07tGvXTvl6wIABCA0NRXZ2NhwdHfUaGJG+mP50l4UTe+NFQURERC+Wnp4OS0srWFhYQiqVYvnypQCAvLy8YvvY2dnDzc0Na9aswP/+9zGio6/g2LHD+P77lXqNrcxPnrWysmJSTxWGKU53qVKKw8yeiIjIpHXr1hOurq4ICemGgQNDYWNjgzp16iE29n6J/T79dDYePEhAz57B+PrruZg160vUq/eyXmPTesS+QYMGEKkkIJrxybNkagRB888mo9DflUl+8CAiIiIlW1tbzJ//ban7ubq64euvlxggoue0TuxDQkLUEnu5XI4nT57g4sWLcHd3R79+/fQeIJFemWDiLCpcY68wyVohIiIiqgC0Tuznz59f7LJHjx5h4MCBsLOz00tQRPpUoaa7JCIiItJRmWvsAcDd3R3vvPMO1q5dq4/NEemXSimOCWb2Is6KQ0RERGWnl8QeyH+a1sOHD/W1OSKDMMm8mYk9ERER6UGpprvUJDc3F9euXcO6detQt25dfcREpFcqubIp5s1M7ImIiEgP9DorzpIlhr3Tl0gnKjX2TJyJiIjIPJVpVhwAEIvFqFGjBrp37w5fX1+9BkekD6b+4FmO2BMREZE+6GVWHCKTVvjmWYUJJs5M7ImIiEgP9HbzLJGpYo09ERERmQpBEDB+/BisWbNC79sudsRe2yfNFsUnz5LpqUA19qYeHxEREZVJRMRWXL78D4KCXtH7totN7DXV1B86dAg5OTlo164d6tWrB4VCgbi4OPz5559wcHAw6JNn4+LisGDBApw9exYA0KFDB0ybNg3Ozs567/fJJ5/g/v372Lhxo/7eABmPyjz2xgujWOLnX5yZ/AcPIiIiwv79kdi6dRMePnwAe3t7DBs2EnXreiMs7EO1dTt37obJk2cAAGJj72P37l/x6qvtDRJXsYl90Zr6jRs34ujRo/jtt9/UprWMj4/HO++8o9MIvzZSUlIwbNgwyGQyjB49GnK5HGvWrMGNGzcQEREBKysrvfWLiIhAREQEWrRoYZD3QuVPMPG7Z0UsxSEiIipWbvRVSNevgSCVGnQ/Ijs72A0bBctGfiWul5AQj8WLv8by5WtQr543Tp48gVmzpmLv3oM4ePBEsf3kcjnmzfsCEyZMxh9/7NN3+ABKcfPs6tWrMXz4cI1z1Xt6emLw4MFYv349JkyYoM/4AADr1q3Do0ePsGfPHnh7ewMAAgMDMWLECOzatQv9+/cvcz+5XI4ffvgB3333nd7jJ9MhKIwdgSZM7ImIiIqTHbEVuX//VT77sreH5aezS1zH1dUN69b9DHd3d6SkpMDCwgIymQxpaWmws7Mvtt/WrRtRr97LaNashfET+/T09GJHxgFAoVBAJpPpJaiiIiMj0aJFC2VyDgBt2rRB3bp1ERkZWWxir22/nJwc9OvXDzdu3EBISAhOnTplkPdBxiGY+jz2Kl90mWB8RERERmTTbyAUUmm5jNjb9Bv4wvUkEgl27NiG/fsj4eRUFY3+G+FPTEzE8OHq/YODuyI0tB8iI3djzRrDlnlrndg3adIEGzduRI8ePeDm5qay7Pbt21i3bp1ByldSU1MRFxeHLl26qC3z8/PDsWPHytwvJycHGRkZCA8PR/fu3dGxY0d9hU+moALV2MMUp+MkIiIyIstGfnBa8K2xw1A6fPgATp8+ic2bI+Do6IT09HTs27cHrq6u2L//mMY+69evwZMnSQgN7QEAyM7OhlgsRkpKCsLCpuktNq0T+48//hhDhgxB9+7d0b59e3h5eUEmk+Hu3bv466+/UKVKFUyZMkVvgRVITEwEALUPEwDg4uKCjIwMpKeno0qVKjr3c3BwwIEDB2BhofXhoApEKPaFqWApDhERUUWRnp4OS0srWFhYQiqVYvnypQCAvLy8YvsMGzYKw4aNUr7+6qvP4e7ugVGj3tNrbFpnsv7+/oiIiMDSpUtx7NgxSP/7OsTBwQG9evXCRx99BHd3d70GBwCZmZkAAFtbW7Vl1tbWAACpVKqW2Jemn1gshljMKf3NlsmX4hSuxTHB+IiIiEipW7eeOHv2FEJCusHW1hYdO76BOnXqITb2Pry8ahs1tlINUb/88stYunQpBEFASkoKRCIRqlWrZqjYAOTX7r+IpqRc135kfgRTL8VRyetNMUAiIiIqYGtri/nzy1YaNHPm5/oJpgidak9EItEL54/XF3v7/LuLc3Jy1JYVtBWso49+ZN5MfsReiw+kRERERJoUm9g3bNgQCxcuRK9evQBo9yRakUiE6OhovQZYs2ZNAEBSUpLassePH8PR0RF2dnZ660dmSCjmZxMkyJnYExERkW5KfPJs7dq1VV4b6gFUJXF0dISnpyeuXr2qtiw6Ohr+/v567UfmR4CJ19gXnglHITdeHERERFShFZvYz5s3T+V10SfRlqfOnTtjw4YNiImJUc5J//fff+Pu3bsYNWqU3vuRmTH1GvvCT81iKQ4RERHpqEzzO+bm5uLkyZMQi8Vo06aNwaaLfPfdd/Hbb79h+PDhGDlyJHJycrB69Wr4+fmhd+/eAIC4uDhcuHABTZs2hZeXl9b9yPyp3DxrivPEF07mWYpDREREOtJ6WhiZTIZPP/0UI0eOVL7u378/xo0bh/feew8hISF4+vSpQYJ0dnbGpk2b0KBBAyxduhTr169HcHAwVq9erXwablRUFKZMmYKoqKhS9aNKwNQT+8LTcbIUh4iIiHSk9RD7d999h+3bt+Ott94CAOzatQvXrl3D0KFD0bBhQ8yfPx9LlizBl19+aZBA69Wrh1WrVhW7PDQ0FKGhoaXup8mRI0dKHR+ZLkElcTZiIMUQVEbsmdgTERGRbrRO7H///Xf07dsXc+bMAQD88ccfyqfNWlhYIC4uDhEREQYLlEgfFCY+Ys8aeyIiItKV1qU4jx49QpMmTQAAWVlZiIqKQuvWrZV19R4eHkhLSzNIkERlYuqlOIWTeZbiEBERkY60HrGvUaMGnjx5AgA4ceIEZDIZOnTooFx+48YNuLq66j1AorJSLcUxwcS+cHy8eZaIiMgs3blzG8OHvwNra2tl28yZn6NDh05624fWiX3Lli2xfv16WFtbY/PmzbC1tUVwcDDS0tKwY8cObN++HQMGDNBbYESGIAj5yb1IXP7PZCgWS3GIiIjM3u3bt9C6dVssWBBusH1ondjPmDEDiYmJWLBgAezs7DB79mw4Ojri/PnzWLBgAZo3b47x48cbLFAiXRWdu14QBIhgQok9S3GIiIiK9TQ+A9HHHiAvx7CDXxbWYjTqUBPVPR1euO7+/ZHYunUTHj58AHt7ewwbNhJ163ojLOxDtXU7d+6GyZNn4Pbtm6hX72VDhK6kdWLv6OiItWvXIjk5GQ4ODsrpIhs2bIht27YhMDDQYEES6ZNCAYglxo6ikMKfPFiKQ0REpOLW34l4eCO1XPZlaS1B9X4lJ/YJCfFYvPhrLF++BvXqeePkyROYNWsq9u49iIMHTxTb7/btW5DJZOjbtxcAoHfvUAwZMkKv8Zf6iVLOzs5ITEzEo0ePUK9ePVhbWyMgIECvQRHpk1BkyN7U6uwLT3cpsBSHiIhIRf02bsiVyctlxN6ntdsL13N1dcO6dT/D3d0dKSkpsLCwgEwmQ1paGuzs7Ivt5+joBF/fhggJeQuPHj3ElCkTUb16DXTv3kt/76E0K58/fx5fffUVrl27BgD46aefIJfLMWPGDEybNg3du3fXW2BEelO0FMfEEnuVyfVZikNERKSiuqcDXh3sY+wwlCQSCXbs2Ib9+yPh5FQVjRr5AQASExMxfPhAtfWDg7siLGwaPv/8K2Vb3br1EBraD3/9ddw4if3ly5cxYsQIeHh4YNiwYVi/fj0AwMnJCRYWFggLC4O9vT3at2+vt+CIDMHkEnsFb54lIiKqKA4fPoDTp09i8+YIODo6IT09Hfv27YGrqyv27z+msU9OTjZWr16BYcNGwcEhv9QnLy9XWdquL1rPY79kyRJ4enrit99+w5gxY5TlDQEBAdi9eze8vb2xYsUKvQZHVFp5uQqc330PUbvuIk+WP/pd9OZZk8udBT55loiIqKJIT0+HpaUVLCwsIZVKsXz5UgBAXl5esX2srW1w9uwprF27Cnl5ebhz5zZ27oxAly76rXbROrG/ePEiQkNDYWNjA5FIdUYRBwcH9O/fH7du3dJrcESlFffvU8REJeHu+Se4eyH/uQtFM3tTHrFnjT0REZFp69atJ1xdXRES0g0DB+bnxnXq1ENs7P0S+82ZsxAxMbfQo0cnhIV9hMGDh6N167Z6ja1UNfYlfV2Qk5MDBZMSMrL0J9nKnzOScwColdibXmKvUmPPvyEiIiJTZmtri/nzvy11Py+v2li8eLkBInpO6xH7wMBA7N27V+MyqVSKiIgIzo5DRleQzANAdnpu/g9qpTimldirfNBgKQ4RkX4UrcMkqgS0Tuw//PBDREdHY/Dgwdi1axdEIhEuX76MDRs2oHfv3oiPj8fYsWMNGSvRC8mkz+vbsjPzE3tTn+6SNfZERPpldfcgqq8JgO3574wdClG50jqxDwoKwooVK/Do0SMsWLAAgiAgPDwcc+fORXZ2NsLDw9GqVStDxkr0QrKs54l9Xo7mJNnkEnvOY09EpFdO+0ZAnPMMDqfnGzsUonJVqhr7tm3b4uDBg4iOjkZsbCwUCgVq1aoFf39/WFiU+llXRHony3qezOfJ/kuSTX5WnMLTXXLEnoiIiHRT6mxcJBLBz88Pfn5+Ku1yuRwrVqzA+++/r7fgiEpLlq0+Yl+0zFKQm1hmr+DNs0RERFR2JSb2aWlpiIiIwKVLlyAIAho1aoRBgwbB0dFRZb1Lly7hk08+we3bt5nYk9Eo5ILK46ZzZQWj36qZvTFvnhUEAZmLFkDx9Ans/zcRkpq1VD95mNqHDiKiik4QgCLTdBOZq2IT+9jYWAwZMgSPHz9W3nx48OBBbN68Gdu3b0fNmjUhk8nwzTffYPPmzZDL5ejVS3+PxCUqrdxs1QdD5MkUEARB/QFVcuMl9nlX/kX2bzvzX1hZwfGLuSzFISIyJEUuINHv0z2JTFWxN89+++23ePz4MSZNmoSTJ0/iwoULWLx4MRQKBb788kukpaVh8ODB2LBhA2rWrIk1a9bg66+/Ls/YiVTIsoskxQIgz1Wo19jnGTGxv/rv85+vXAagesOsUMJT64iISAtFRnNEeVlGCoSo/BU7Yn/+/Hn06dMHo0ePVrZ17doVOTk5mDlzJiZPnox///0XQ4cOxaRJk2BtbV0uARMVp/BUlwXkeYLaA6qMOWIvf/RQ+bOQI/vvh0LxMLEnIiqbIom8KC8LgrWTkYIhek6hUGDt2lXYt28PsrOz0Llzd/zvfxMhFms9SeULFZvYp6SkICgoSK29WbNmyMvLw5kzZ/D999+jY8eOeguGqCzURuwBKOQKtdEbeZ7x6tgVT54ofxbS0/JH6zliT0SkN+KcZ6oNuRyxJ9OwefMGnDp1EqtXbwAAfPTROOzfH4nu3fVXyl5sYp+XlwcbGxu1djs7OwDAqFGjmNSTScnNKmbE3oRq7BVPkwq9UEDIzFD94JGbW/5BERGZEVFOmuprluJUeE/u3cbl/b8iN9uwv0tLG1s07toHNeq8/MJ19++PxNatm/Dw4QPY29tj2LCRqFvXG2FhH6qt27lzN0yePAN79+7C9Omfolo1ZwDAwoVLYGEh0et70Hny+aZNm+ozDqIyKzyHfQFFrvrovMKIM88UHrEHACEtrciIPRN7Q9oasxE3Uq9hlM9YeDnUNnY4RGQATOzNT/TR3xF/5UK57MvSxhavDh9f4joJCfFYvPhrLF++BvXqeePkyROYNWsq9u49iIMHT2jsI5VKkZAQjwcPErBw4VfIyspCz569MXLkGL3Gr3Nir896ICJ9kGkcsddw86yRRuwFhQKKJ0kqbYrUVEAo9EGDpTg6O/bwCBZfWYgA50B83nQuJCLVUZDbaTex6sYPAIDMvEx83WKJMcIkIgMTZ6leZ0V52UaKhPSl0evdkJeTXS4j9g1f7/bC9Vxd3bBu3c9wd3dHSkoKLCwsIJPJkJaWBjs7e419MjLSAQBHjx7Cjz/+hMzMTHz88Xh4eNQsn1IcALhz5w6ioqJU2tLT8wO7ceOGxqfNNm/eXG/BEZWG5sReQNHbZ+VGmhVHSE1VS9wFaSaEwvPq5zKx19XiK18jLTcNJxNP4GrKv2js3ERl+eXkf5Q/n38ShTxFHizEfGI2kbkRZz5Wec3EvuKrUedldBw72dhhKEkkEuzYsQ3790fCyakqGjXKf2hrYmIihg8fqLZ+cHBXjBr1HgBg8ODhcHR0gqOjE3r3DsXJk8fLL7H/8ccf8eOPP2pctmDBAo3t165dK3tURDrQePOsxhF745TiqNTX/0fIzlYZsWcpjm6Sc5KRlpuqfB397KpaYp+QGa/yOj4zDnWq1C2P8IioHEkyE1UbWIpDenb48AGcPn0SmzdHwNHRCenp6di3bw9cXV2xf/+xYvtVqeKIjIwM5Wu5XK52H2BZFZvYjx9fcn0RkanRfPOsQv3mWSON2BetrwcAZGcBCtXpLgWZDLIzpyB5qQ4sar9UfgFWYHfTY1ReP5Q+UFvncbbqKN6d9NtM7InMkFiqmtizxp70LT09HZaWVrCwsIRUKsXy5UsB5E88U5LOnbti3brVCAgIhFSaid9+24nRo8fqNbYKk9jHxcVhwYIFOHv2LACgQ4cOmDZtGpydnfXST9ftk+kouHnW0lqC3Jz8nzUl8caqsVckPVZrE6RZqjX2goCM8K+Rs/c3wMYG1dZtgaSWZzlGWTHdKZLYJ2Y9VFvnSbbqNyZ30mPQEW8YNC4iKn/qpThM7Em/unXribNnTyEkpBtsbW3RseMbqFOnHmJj78PLq/iJGT74YAKWL1+KIUP6Qy7PQ58+/dC584tr+kujQhSYpqSkYNiwYZDJZBg9ejTkcjnWrFmDGzduICIiAlZWmh8VrW0/XbdPpqWgxt6miqUysc8fsVdN5PNy1Ut2yoOmEXuh6Ig9kJ/UA0B2NrJ/2wn799WnziJVRUfsnxWdxxpAeq7qTBmxGfcNGRIRGYk485HKa1Gu1EiRkLmytbXF/PnflrqflZUVJkwIw4QJYQaIKl+FSOzXrVuHR48eYc+ePfD29gYABAYGYsSIEdi1axf69+9fpn66bp9MizKxd7BE+pP8m6UUCkGtxv7JvQxciLwPj/pO8PCpWm7xKZ7mJ/YiO3sI0kwAgJAlVR2xLyLnxJ+wG/c/iESicomxolJL7GUpausUTezjMpnYE5kdQQFJmurftkiWrvxZnBYHhW11wNKuvCMjKhcVYs7KyMhItGjRQpl0A0CbNm1Qt25dREZGlrmfrtsn48lMycE/v8fi+l8PlU+SLSjFsbZ//nlVU9lNckImbp9+jBMbbyHxTpra8lKTyzQ2C3l5Kk+SLUjsxe4egCR/KkYhK1vtybiFKeLjoIiLLXuMZkwuyHEv/a5KW6rsmdo6GbkZKm0JmfHIU3AWIiJzIk6PV5sFR5yTf2O9w7HpqL6xNapvaAXJk2hjhEdkcCaf2KempiIuLg5+fn5qy/z8/HDlypUy9dN1+1Q8QRDUyl+KkuZl4ve4vdh4ey3uZ9wrdhuCIEChUB3RzkjOxpFV13Dz70Rc/iMel/bHQZ6rgPy/h1FZ21s+346i5Fgu7L6v/GAgCAJk0lzIMnPyR/rzGwGF5tIdUeZjVN0Ziho/esMxcjhEWcn/7VMB6dZNSO4ejOTuwcj+7ReIpElQJOXf0CWuUR2i/8q7BFmO6nSXGsjOnMpfNycHuZcuIvfqv1Ck6+EDSQWUp8jD6ccn8dPNlfjr0Z8QBAEPMhOQo8gBALjbegAAchQ5kMlzlP0yczOV057Wq5L/AT5PyMND6QMoBAUOxP+O76LDsS9uDzJzMzXvXJ4Li4dRsIw9BlG2+jcCRFTOFHKI0xMgSbkN/Pch3fKR+kOMRDmpkCTfgu3VjQAAcXYynH5/FyhcopOXDcgyi73eE1UUJl+Kk5iYnwy5ubmpLXNxcUFGRgbS09NRpUoVnfrpun1T8POEj5ArTwdQvmUaxaehusXxJ05qtVeRIIdIkEMQW0Io9Jn0yh/A9ch05Fnk/46u7c+GQmwDADi9NgtyiW2xW81+Bmz9SAZAAUFspbJdEQSU9G7zVQXwBnAdwP4wCBD91w+Ab738/z+0D8LhfRDZWgJN/CHIMiHyzU8uhVtXILKzAJr4F7+L4wcgHD9Qzr9l02YJ4C5O4y7WAwBGoEGhpU4AgG2/j1PpU3id9v/9fOz3z1V6PcW/2IlfDREyERmcCPnX7M6qzdczgd3z1NvPmtYkIVTxWEqqYMBi03rYodaJ/dChQzFu3Di0bt1a4/IjR45g0aJFei9dyczMHz2ztVVPzqytrQHkP6a3aOKtbT9dt29sSQ9iIct7AsA4c7IbgwBAEAEQ1Msn8iRiQMj/XSpEUP4sFz//uTiKgoxZUJ05Qde5czT2K1LqX9zPpd4uERERGYUsT4YnD+NRw8N0Zq8rNrHPyspCSsrzr5vPnj2LN954Ay+9pD6vtkKhwPHjxxEfH6+2rKyKlmFoIharVxRp20/X7RubS83acHSvg+ykjBeuKyopJSx20YvSSOG/bZdMDgGK/9aVABAJIggiASIBkEAMMYBcKKAQASIhf3sKUX7RhFgAxP/tQRBZQi6uAkFsBwh5sMxLgGVeIqQ2rwAiCURCNgRR/ii9Vc4NyKx9AQCWubHItVSfesoyNxYSeSqyrRsAovzSHbEiExJ5cv7NrCIJBEgAkQUA4fnovQAIIpGyLl4QiSFADDEUEAkKSOT5Hz7yLIEcGwF2mSJYFHnmVIa9BLZZckgUQI6VCBZyAZIi3/7mWoogtRPDKVV1gdRWDLkEkMgBi//uH1B+4DFzAvLPDwlEEEOEHMhVzlIxRLCACLL/PuxKIFKuYQ0JciBX/pwLhfK8xH/btRYk+b87KPD8/Bep/K+gPOMFtb8rkbLYRwTNfxn8aFY6leCkNiEV/ewUFxnkksMCYsjV/k4FiKCABBKoDhApVL6tfb42UUlsXBxMKqkHXpDYh4SEID09/25ykUiEuXPnYu7cuRrXFwQBbdu21XuA9vb2AICcnBy1ZQVtBevo0k/X7ZuCkE++MnYIWslT5EEkEkEikpR9W7kKSCQiiMT5l97jG27i0a1UlXWCJ0/AweVXAQD+wbVw5VACAKBKdRukP82/qarr//zh6GoLWXYensZlws7REo6utsXOPiNOi4dYmog8F39AYv3iQHOlsEy6jMzL8Uj7aqHKIqfvVyFj4VzI79+FVYeOkN+7C/k91Zs/rbv2QJWZn0G6eQOka1ZAXM0Z9hMnw7rday/edyWxJ3YXwq88P7ajfceirqMPpkZ9rLbupIBp+OHf+QCARS2X4Z+nF7Dx9lrl8nEN/oe+9dQfA05EFYMoKxl255dCkvEQ0sajkFezBRyOTILttW0q62W0+xxZgaPzZ895dgeChR0UVWoaKWoi/Ss2sXd2dsbXX3+Nf//9F4Ig4Pvvv8cbb7wBX19ftXXFYjGcnZ3Ro0cPvQdYs2b+H1xSUpLassePH8PR0RF2durTVmnbT9ftk/YsxPq7lcPCUvXbE9e6VdQSe/uqz587UHhWHE//ahBLRHBwtoGja37plZWNBTzqO71wvwpHTygcS/Gp3NIOuTVbwdIlF6LvVkJIfZbfLpHAor4PYJX/LYEgk2mcFUdcwwUAYDdoKGz7DQAsLTnlZRE9vXrjTtptRMbtRtMazdDnpX6ISb+tcd303OfT3dlIbNGrdgj+iN+Hx9mJaF6jJULrckpboopMsHVGZrvPVdusq6qtl1e9Yf4PIjHk1V42fGBE5azEjKt9+/Zo3749AODBgwcYMGAAAgMDyyWwAo6OjvD09MTVq1fVlkVHR8PfX/NNh9r203X7ZBpc6qre+2BlK4GVrYXyHqrCib1YIobf6+U7MiOytIR1l27I3r4VAGDZtBlEtrYQWf734SM3F9BQDiauXv35NviANI1EIhE+8g/D+EYTIPnvw6OdRPOH8AyVxN4GNWxcsPa1zXic/Rie9l56+TaJiEyLYK0+aKNM7InMlNbF4/PmzSv3pL5A586dcerUKcTEPH8Izd9//427d++ie/fuZe6n6/bJ+Kp52MPC6vlpXDASL/6vVEchf540G2vA227YSFi2bA2LgEA4TJqa32j5fMRe03Sc4uo1yjPECk1S6BshOwvNiX2a7Pm3OjaS/HsxbC3s8JJDHSb1RGZKYeOs8lpu7wbB1rmYtYnMQ6lqJI4fP449e/bgyZMnkMvV53oViURYv3693oIr8O677+K3337D8OHDMXLkSOTk5GD16tXw8/ND7969AQBxcXG4cOECmjZtCi8vL637lWY9Mj1iiQgudavg4Y38xK1GbQdlu0IuaHxAVXkTOzrB6RvV6bCUo/BajNiT9myLSexTCj2JtiCxJyLzJq9SS+V1notxBieJypPWif3mzZsxZ84cAED16tVhVY7lAc7Ozti0aRPmzZuHpUuXwsbGBsHBwZgyZYoyjqioKEyfPh3z5s1TJvba9CvNemSaAjp54tkDKSysJXi5Zf7zCApurlXkFZr5xIRK1JUPqMqV5c/CUwRH7HVT3Ij9M5XEvvjnGhCR+VA4qs6IlufWxDiBEJUjrRP7DRs2oEGDBli1ahVq1Cj/pKNevXpYtWpVsctDQ0MRGhpa6n6lXY9MT1UPO/Sa0kSlTSwpKMUpNGJvSpn9f6U4kOUCGp48y8ReN5ZiS1iKLZGrUJ1jNDXnmfJna21mNSKiCk/uVEflda77K8YJhKgcaV1j//DhQ7z99ttGSeqJSktTjb0pKbh5VtOIvdWr7SGyYbmIrmw1jMin5aYByC/DEYtM77kURGQAYgmkTfOfLpvr0Ry5NVsZOSAiw9N6xL527dp48uSJIWMh0htlKY7cNEtxCqa7hEymDMzq9U6w6RUCyyCOKpWFrYWdMpEvkJGX/yA31tcTVS6ZradBGjgKglUVQMwb5cn8aT10NWbMGGzcuBG3bt0yZDxEeqG5FMdIwWigMmL/382zIocqsGreEiIL/c37XxlpmvJSIeTf7M/6eqLKR7BzASz4oZ4qB60ziPPnz8Pe3h69e/dG3bp14ezsrPbAHEPNikNUWppH7E0osy9UYy9I8keRCmKmsrG1KD5554g9ERGZM60T+xMnTgAA3N3dkZWVhYSEBIMFRVRWz0fsTbTG3ir/Bk4hVwaR9X83c7L2Wy+Km/IS4Ig9ERGZN60T+yNHjhgyDiK9Uib2JjrdpbLGPi8PKHgmhEkFWHHZFvP0WQCw4dfxRERkxnQaIkxMTMSlS5eQnp4OmUwGhYYH7BAZk6ZSHFOssQcAIScn/wcJb+zSh+Lmsgc4Yk9EROatVIn9+fPnERoaig4dOmDAgAG4cuUKzp49iw4dOmDfvn2GipGo1DRPd2k6mb2ooMYeAGT/JfZiluLog6bpLguwxp6IiMyZ1pnE5cuXMWLECGRmZmLYsGHKdicnJ1hYWCAsLAx//vmnQYIkKi2RhllxTKrSRcMTjUWcik0vePMsERFVVlon9kuWLIGnpyd+++03jBkzBoKQnzAFBARg9+7d8Pb2xooVKwwWKFFpiE2+FMdSvVHCEXt94M2zRERUWWmdSVy8eBGhoaGwsbFRmzbQwcEB/fv35xz3ZDI0zWMvMqXM3lJ9xJ419vpR4s2zHLEnIiIzVqohQisN5QMFcnJyeBMtmQyRphp7E8rrRVYaRuxZiqMXJd08W1L9PRERUUWndWIfGBiIvXv3alwmlUoRERGBgIAAvQVGVBYFI/ZyU53uUsOIvYg3z+pFiTfPllB/T0REVNFpnUl8+OGHiI6OxuDBg7Fr1y6IRCJcvnwZGzZsQO/evREfH4+xY8caMlYirWmc7tKEiDR9+8Uae70oqcaeI/ZERGTOtH5AVVBQEFasWIHPPvsMCxYsAACEh4cDAFxcXBAeHo5WrVoZJkqiUiq4eVZQFL551oSG7DXdPMtSHL0oaVackpYRERFVdFon9gDQtm1bHDx4ENHR0YiNjYVCoUCtWrXg7+8PC4tSbYrIoDTePGtCeb2IN88aTEk3z5a0jIiIqKIr1Xf/Dx48wKJFi+Dp6Ylu3bqhR48eOHv2LBYtWoSnT58aKkaiUiuYx95Uabp5ljX2+lHyk2c5Kw4REZkvrTOJmzdvok+fPli7di0ePnyobE9LS8PmzZsREhKCuLg4gwRJVFoFpTiFFZ2m1ag4Ym8wJdXRl1R/T0REVNFpndgvWrQI9vb2iIyMRIMGDZTtYWFhiIyMhKWlJb755huDBElUWmINI/amlNdrfEAVR+z1gtNdEhFRZaV1JvHPP/9g2LBhqFOnjtoyLy8vDB48GFFRUfqMjUhnIk0j9hrajEbjk2c5Yq8P1iWU2/DmWSIiMmdaJ/aCICAnJ6fE5dnZ2XoJiqisNJfiGCGQYmgasWeNvX6IRcUfRxuO2BMRkRkr1QOqtm3bhrS0NLVlmZmZiIiIQGBgoF6DI9KVxtF5E0rsoWkee053aXAsxSEiInOm9RyV48ePx+DBg9GzZ0/06tULL730EkQiEWJjYxEZGYmkpCTMmzfPkLESaU1zjb3pZPaap7vkiL2+uNi4Iin7sUqbRCSBpVhDCRQREZGZ0DqxDwwMxNq1a7FgwQL89NNPEITn84M3aNAA8+bNQ1BQkEGCJCotjTX2ppPX8wFVBta3ztv44foylTZbiZ1JfbgjIiLSN60T+5SUFDRr1gwRERFITk5GQkICFAoFPDw84OrqasgYiUrN5EfsJZL8m2XlctU20ou36r6Nlq5tcCDhd2yJ2QAAsJZYGzkqIiIiw9I6se/Tpw/69euHDz74AM7OznB2djZkXERlovH+SdPJ6/NZWqok9pzuUn/EIjFqO7yEqlbVlG0WYj4dm4iIzJvWmURycjJcXFwMGQuR3og1JMkmNd0lNNTZs8Ze7xytHJU/i0v3oG0iIqIKR+v/0vXq1Qvbtm1DfHy8IeMh0gtNI/YmVImTr2idPWvs9c7JsqryZ1MqxSIiIjIErb+bFovFuHPnDrp06YLatWujevXqaqOiIpEI69ev13uQRKVl6jX2QP5c9kLh16yx1ztPey/lz33rDDBiJERERIandWJ/8uRJVKuWX6+ak5ODBw8eGCwoTbZs2YINGzbgwYMHeOmllzB27Fj06NHDIH1jY2PRs2dPrFq1Ci1bttTXW6BypOkBVSZZY18Ya+z1rpa9J2a/Mh/Z8my87hFs7HCIiIgMSuvE/siRI4aMo0Rr1qzBwoUL0bVrVwwfPhwHDx7Exx9/DJFIhO7du+u1b1paGj744IMSn7JLpk/zdJemldmLij6kiiP2BtHW7TVjh0BERFQudBoiTExMxKVLl5Ceng6ZTAaFQqHvuJTS0tLw3XffoWfPnliyZAkGDBiAVatWoVmzZli4cCHkhWcVKWPfmJgYvP3227h586bB3g+VD82lOEYIpASiIiP2RV8TERERlUapEvvz588jNDQUHTp0wIABA3DlyhWcPXsWHTp0wL59+wwS4JEjRyCVSjFw4EBlm1gsxjvvvIOHDx/i4sWLeun766+/onfv3nj27Bn69etnkPdC5UfjiL2JzYqDorPiWHA6RiIiItKd1on95cuXMWLECGRmZmLYsGHKdicnJ1hYWCAsLAx//vmn3gO8cuUKAMDPz0+lvVGjRirLy9r3xo0b6NGjB/bs2YOmTZuWPXAyKk0j9qZWYy+yKjJiz8SeiIiIykDrTGLJkiXw9PTEzp07IZVKsW7dOgBAQEAAdu/ejYEDB2LFihVo3769XgN8/PgxnJycYGtrq9JeMKd+STfxlqbvxx9/DKuiNc9UYWmusTdCICVRG7FnKQ4RERHpTuvE/uLFi3j//fdhY2ODrKwslWUODg7o378/li5dqvWOk5KSSlxuZ2cHe3t7ZGZmwsbGRm15QVvRWAorTV8m9eZF06w4JnfzbNGaeo7YExERURmUKpMoKfnNyckp1U207dq1K3H52LFjMXHiRCgUCo0JWUFbSclaWfpSxaaxnt7Uft1F/p548ywRERGVhdaJfWBgIPbu3YuhQ4eqLZNKpYiIiEBAQIDWO54zZ06Jyxs2bAgAsLe3R3Z2ttrygtF2e3v7YrdRlr5UsVWIB1QVLb2RcMSeiIiIdKd1JvHhhx9iyJAhGDx4MDp16gSRSITLly/j1q1b2LhxIx48eIAvvvhC6x1rO/OMh4cHUlNTIZPJVL4xePz4MQDAzc3NIH2pYqsYNfZFb57lPPZERESkO61nxQkKCsKKFSvw6NEjLFiwAIIgIDw8HHPnzkV2djbCw8PRqlUrvQfo5+cHQRBw7do1lfaC1yV9S1CWvlSxVYga+6KlbSzFISIiojIo1Xf/bdu2xcGDBxEdHY3Y2FgoFArUqlUL/v7+sDDQjX/t27eHtbU1Nm7ciMDAQAD5tfNbtmxBrVq10KRJE4P0pYqtQtTY8+ZZIiIi0qMXZhK5ubm4ffs28vLy8PLLL8PW1hZ+fn5qc8MbSrVq1TBmzBgsW7YMgiCgVatW+OOPP3Du3DksXrwYEsnz8oVDhw4BAIKDg0vdl8yLxhp7E3tAlajIdJci1tgTERFRGZSYSaxbtw7ff/89MjIyAOTPivPOO+9g0qRJBhuh1+SDDz6Ara0tNm/ejIMHD6JOnTpYsmQJunbtqrLe3LlzATxP7EvTl8yLSEORmYlV4gBWHLEnIiIi/REJgiBoWrBr1y5MmzYNtWrVwuuvvw6xWIwzZ87g5s2bGDJkCGbMmFHesZq03Fw5nj2TGjsM+o80NQd7v7ms0tbtowBUqaH+XANjkW5cB+nK5crXNU6cNWI0REREVBG4uFQpdlmxQ4RbtmxBkyZNsH79elhbWwMABEHAxIkTsW3bNoSFhfGhTmSyKkKNPeetJyIiIn0qdlacmJgY9OrVS5nUA/mzigwfPhwymQx37twplwCJdKF5VhwjBFICkR2fo0BERET6U2xin5WVhSpV1If6PT09IQgC0tLSDBoYUVlonMfe1G6e5QPSiIiISI+KTewVCoXGeb8LZpKRy+WGi4qojDQ/edYIgZRAZGdn7BCIiIjIjGj9gCqiikRzjb1pZfYcsSciIiJ9KnF+vWfPnuHBgwcqbampqQCA5ORktWUAULNmTT2GR6SbilFjzxF7IiIi0p8SE/u5c+cq54YvKiwsTK1NJBIhOjpaP5ERlYFILIJIBBSezFVTaZkxiV1clT9btXvNiJEQERGROSg2se/Tp095xkGkdyKxCIK8cGZvvFg0ETtVhf3HU5D96y+w6T/Q2OEQERFRBVfsA6qodPiAKtOzc/Z55MkUytchM4JgZcunuxIREVHFVdIDqnjzLJmtojfQmtp0l0RERET6xMSezFbRG2iZ2BMREZE5Y2JPZkttxJ55PREREZkxJvZktoo+pErTFJhERERE5oKJPZktlRF7EUtxiIiIyLwxsSezVXiEnqP1REREZO6Y2JPZEokL/8zEnoiIiMwbE3syW4Vr7HnjLBEREZk7JvZktgqP0nPEnoiIiMwdE3syW4VH7FljT0REROaOiT2ZLY7YExERUWXCxJ7MllglsTdiIERERETlgOkOmS2RhCP2REREVHkwsSezxXnsiYiIqDJhYk9mizX2REREVJkwsSezpVJjz7yeiIiIzBwTezJbrLEnIiKiyoSJPZktlRp7CRN7IiIiMm9M7MlsSSyfn94csSciIiJzx8SezJbEgqU4REREVHlUmMR+y5Yt6Nq1Kxo3boxevXohMjJSr30zMjIwZ84cvPbaa/D390fHjh0RHh4OmUymz7dB5Uhs8fz0llhUmFOdiIiISCcVIttZs2YNvvjiC/j6+mLGjBlwdXXFxx9/jH379umlryAIGD9+PLZs2YKOHTti5syZaN68OVasWIGJEyca8q2RARUuxSk8ek9ERERkjkSCIAjGDqIkaWlpaN++PTp27IhFixYBABQKBYYMGYKEhAQcPnwYEomkTH0PHTqEDz74AJ9++ikGDRqk7B8eHo4ff/wRW7ZswSuvvFJinLm5cjx7JtXTuyZ9uHb8If49GA8AqOlbFe0G1zdyRERERERl4+JSpdhlJj9if+TIEUilUgwcOFDZJhaL8c477+Dhw4e4ePFimfuePXsWABASEqLSv1u3bgBQ4j7IdEksC43Sc8CeiIiIzJzJJ/ZXrlwBAPj5+am0N2rUSGV5Wfq+//772LVrF+zt7VXWS0lJAQBYWFjoGj4ZEevqiYiIqDIx+Yz18ePHcHJygq2trUq7i4sLAODBgwdl7lu1alVUrVpVrf/WrVsBAEFBQTrHT8bDxJ6IiIgqE6Ml9klJSSUut7Ozg729PTIzM2FjY6O2vKAtKyur2G2Upe+uXbvwxx9/oGXLlggMDCwxVjJNFlbPE3t5rsKIkRAREREZntES+3bt2pW4fOzYsZg4cSIUCgVEIvUC6YI2TcsK6Nr38OHD+OSTT+Di4oL58+eXGCeZLkeX59/UCAqTvkeciIiIqMyMltjPmTOnxOUNGzYEANjb2yM7O1ttecFoe9G6+MJ06bt3715MmzYN9vb2WLVqFWrWrFlinGS6qtSwQXUve6Q8kMK/Uy1jh0NERERkUEZL7Pv166fVeh4eHkhNTYVMJoOVlZWy/fHjxwAANzc3vfX9+eef8cUXX8DJyQnr1q1DgwYNtH4/ZHpEYhFeH90QeTlyWNma/O0kRERERGVi8ncX+vn5QRAEXLt2TaW94HVAQIBe+u7atQufffYZatSogc2bNzOpNxNisYhJPREREVUKJp/Yt2/fHtbW1ti4caOyTaFQYMuWLahVqxaaNGlS5r4xMTGYNWsWnJ2dsXHjRnh7exvq7RARERERGYTJD2VWq1YNY8aMwbJlyyAIAlq1aoU//vgD586dw+LFi1WeOnvo0CEAQHBwcKn6Llu2DDKZDK+++iouXbqES5cuqcTg6+vLEXwiIiIiMmkiQRBMfroQQRDw008/YfPmzXjy5Anq1KmD999/H127dlVZr2PHjgDynzhbmr5t27bFkydPit3/+PHj8b///a/EGHNz5Xj2TKrL2yMiIiIi0oqLS5Vil1WIxL4iYGJPRERERIZWUmJv8jX2RERERET0YkzsiYiIiIjMABN7IiIiIiIzwBp7IiIiIiIzwBF7IiIiIiIzwMSeiIiIiMgMMLEnIiIiIjIDTOyJiIiIiMwAE3siIiIiIjPAxJ6IiIiIyAwwsSciIiIiMgMWxg6AgLi4OCxYsABnz54FAHTo0AHTpk2Ds7OzXvrpun1Tpev7OXHiBH744QdcvXoVYrEYgYGBmDBhApo0aaKyXt++ffHvv/+q9e/SpQuWLl2qt/dRnnQ9Ztoei8p+jsXHx6NTp04lbnPDhg1o2bIlAPM8xwp88sknuH//PjZu3PjCdSvrNayw0hyvynwNK1Ca41VZr1+FaXO8Kvv1S9u/q6JM9frFxN7IUlJSMGzYMMhkMowePRpyuRxr1qzBjRs3EBERASsrqzL103X7pkrX93P27Fm8++67qF+/PiZOnIi8vDxs2bIFgwcPxpYtW9C4cWMAgCAIiImJQXBwMDp37qyyjVq1ahn8/RmCrsdM22PBcwxwdnbGwoUL1dpzcnIwe/ZsVK9eHQ0aNABgnudYgYiICERERKBFixYvXLeyXsMKK83xqszXsAKlOV6V9fpVmLbHqzJfv7T9uyrKpK9fAhnVt99+KzRs2FC4ffu2su3kyZOCj4+PsG3btjL303X7pkrX99O7d2+hQ4cOglQqVbYlJSUJzZs3F4YPH65si42NFXx8fIQdO3YY5g0Yga7HTNtjwXOseHPmzBEaNGggREVFKdvM8RzLy8sTli1bJvj6+go+Pj7C4MGDX9insl7DBEG341WZr2G6HK/Kev0SBN2OlyaV4fql7d9VUaZ8/WKNvZFFRkaiRYsW8Pb2Vra1adMGdevWRWRkZJn76bp9U6XL+0lNTcX169fRtWtX2NraKttr1KiB5s2b4+LFi8q227dvA4DK9is6Xc8BbY8FzzHNbty4gU2bNqFPnz5o1qyZst3czrGcnBz06dMHy5YtQ+/eveHm5qZVv8p6DdPleFXma5iu51dlvX7peryKqgzXr9L8XRVlytcvJvZGlJqairi4OPj5+akt8/Pzw5UrV8rUT9ftmypd34+DgwP279+P4cOHqy1LSUmBRCJRvr516xaA5xctqVSqh8iNpyzngDbHgudY8cLDw2FjY4MJEyaotJvbOZaTk4OMjAyEh4djwYIFsLB4cYVnZb2GAbodr8p8DdPleAGV8/oF6H68iqoM16/S/F0VZurXLyb2RpSYmAgAGj9Ru7i4ICMjA+np6Tr303X7pkrX9yORSFCnTh21ftevX8eFCxcQFBSkbLt16xbs7e0xb948BAUFISgoCMHBwRVy5AbQ/ZgB2h0LnmOaXb9+HUePHsWAAQPg6uqqsszczjEHBwccOHAA3bt317pPZb2GAbodr8p8DdPleAGV8/oF6H68Cqss16/S/F0VZurXL948a0SZmZkAoPIVUAFra2sA+Z+Gq1SpolM/XbdvqvT5fjIzMzF16lQAwJgxY5Ttt2/fRmZmJtLT07Fw4UKkpaVhw4YN+Pjjj5Gbm4uQkBA9vJPyU5Zjps2x4Dmm2datWyGRSDB48GC1ZeZ2jonFYojFpRsjqqzXMEC346VJZbmG6Xq8KuP1C9DP+VWZrl9FFfd3VXQdwHSvX0zsjUihULxwHU1/oNr203X7pkpf7ycrKwvjxo3D9evX8d5776nMGNC/f38oFAoMGjRI2dajRw/07NkTX3/9NXr16lXs13OmqCzHTJtjwXNMXXZ2Nnbv3o2OHTtqnCXC3M4xXVTWa5i+VKZrmK4q4/VLHyrz9aukv6vCTP36VbnOWBNjb28PIL8mrqiCtoJ1dOmn6/ZNlT7eT1paGkaOHIkzZ87grbfewsSJE1WWDxw4UOWCBQA2Njbo3bs3njx5orxxqKIoyzHT5ljwHFN35swZSKVSdO3aVeNyczvHdFFZr2H6UNmuYbqqjNcvfais168X/V0VZurXLyb2RlSzZk0AQFJSktqyx48fw9HREXZ2djr303X7pqqs7+fp06cYOnQoLly4gLfffhtfffUVRCKRVvsueJBERbtRyBDnQOFjwXNM3Z9//glLS0t06NChVPuuqOeYLirrNaysKuM1TN/M+fqlD5Xx+lXavytTv34xsTciR0dHeHp64urVq2rLoqOj4e/vX6Z+um7fVJXl/WRkZGDUqFG4du0ahg8fji+//FLtDzcxMRE9evTAd999p9b/7t27AABPT88yvovypesx0/ZY8BxTd+HCBQQEBMDBwUFtmTmeY7qorNewsqis1zBdVNbrlz5UtuuXNn9XRZn69YuJvZF17twZp06dQkxMjLLt77//xt27d0u8q13bfrpu31Tp+n6+/PJLXLt2DUOHDsX06dM1ruPm5ob09HREREQgIyND2f7w4UPs3LkTLVu2hIuLi/7eTDnR5ZiV5ljwHHsuNzcXt2/fRqNGjTQuN9dzTBeV9Rqmq8p8DSutynz9KovKeP3S5u9KE1O+fokEQRAMsmXSSnJyMnr27AmJRIKRI0ciJycHq1evRu3atfHzzz/DysoKcXFxuHDhApo2bQovLy+t+5VmvYpCl+MVExOD7t27o0qVKpgxY4bGG3t69+4NADh06BA++OAD1K9fH/369UNmZiY2b96M3NxcbN26tUI+lEPXc0zbY8FzzEvZNzY2Fm+88QYmT56M0aNHa9y+OZ5jhRXcdLdx40ZlG69hxdPmeFX2a1hh2p5flfX6VZS2xwuofNcvbf+uKtz1yyDPs6VSiYmJEUaPHi00adJEaNWqlTB16lTh6dOnyuU7duzQ+AjnF/Ur7XoVRWmP15YtWwQfH58S/xV28OBBoW/fvoK/v7/QrFkz4f3331d5HHRFpOs5pu2xqOznWIFLly4JPj4+wtatW0vcvjmeYwVef/11tUfY8xpWPG2OF69hz5Xm/Kqs16/CSnO8Ktv1S9u/q4p2/eKIPRERERGRGWCNPRERERGRGWBiT0RERERkBpjYExERERGZASb2RERERERmgIk9EREREZEZYGJPRERERGQGmNgTEREREZkBJvZERCZs2rRp8PX1xebNmzUuj4+Ph6+vL5YtW6bT9p8+fQqpVFriOsuWLYOvry/OnDlTYgzTpk0r9f4L3l9hMpkMiYmJJfY7c+YMfH191f75+/ujU6dO+Oqrr5CamlrqeEry999/o1u3bvD398c777yDnTt3ajwucXFxet0vEZG2LIwdABERvVh4eDi6dOmCGjVq6G2bf/75J8LCwvDrr7/Czs5Ob9stjbfffhutW7dWvk5ISMDIkSPx3nvvITQ09IX933jjDbzxxhvK1zKZDFeuXMGmTZtw7tw5REREwMKi7P+pUygUmDRpEiQSCaZPnw53d3f4+Phg4cKF8Pb2Vq43atQouLi4YP78+WXeJxFRaTGxJyKqANLT0zFv3jwsWrRIb9u8fPky0tLS9LY9XQQFBSEoKEj5Oj4+Hvfu3dO6v6+vL3r37q3S1q9fPzg4OGD16tX4448/0KNHjzLHmZSUhOTkZIwYMQKDBg1Stnt5eams99dff6FPnz5l3h8RkS5YikNEVAF07NgRe/fuxalTp4wdSoXQvXt3AMDFixf1sr3c3FwAgL29vV62R0RkCEzsiYgqgE8++QS2trb4/PPPIZPJXrj+xYsXMWLECOWI+MiRI3H58mXl8mnTpuG7774DAHTq1AlDhgzRa7xDhgzBqFGjcPz4cYSGhiIgIAAdOnTAsmXLoFAoVOIoqLHfuXMnhg4dCgCYPn26Wu19aYjF+f95y8vLA5B/n0BAQAAOHjyItm3bIigoCBEREQCAlJQUfP7553j11Vfh7++PLl26YOXKlZDL5cq+nTp1AgB89913yrr6wjX2BfcZAMCvv/5a4j0JRESGwsSeiKgCqFWrFt5//33cu3cPK1euLHHdkydPYsiQIUhPT8dHH32EcePG4cGDBxg0aBDOnTsHIL+2vaA2ffr06Rg7dqzeY7558yYmTJiAli1b4pNPPoGXlxe+++47bN26VeP6zZs3V8bx9ttvY+HChTrvu+CbDT8/P2VbXl4ePvnkEwwbNgyjRo3CK6+8gtTUVAwYMAC//PILunTpgunTp8Pb2xuLFi3CpEmTAOTX8U+fPl35c9G6egBwdnZWxtusWTON6xARGRoTeyKiCmLEiBGoX78+Vq5cifv372tcR6FQ4LPPPkNAQAC2bduG4cOHY8yYMdi5cyc8PDwwZ84cAPm17QUjzMHBwWjbtq3e4338+DG++eYbTJ06FW+//TZWr14NJycn7NmzR+P6Xl5eaNOmDQCgSZMmarXzmmRlZSE5OVn57/bt21i/fj2WLVsGDw8PZUkOkH9sBg8ejDFjxmD8+PGoV68eVq1ahXv37mHx4sX45JNPMGjQICxfvhzvvPMOfv/9d/z5559o0KABgoODATyv6S96E7OdnZ0yXi8vL43rEBEZGhN7IqIKwtLSUlmK8+WXX2pcJzo6GnFxcQgODkZqaqoy4c3Ozsbrr7+Oa9eu4dGjRwaJTyQSqby2tbVFhw4dlK+tra1Rt25dPHnyRG/7XLNmDVq3bq3816NHD8ydOxe+vr5Ys2aNWk18u3btVF4fOXIE3t7eysS9wPvvvw8AOHz4sN5iJSIyNM6KQ0RUgTRr1gx9+vTBzp07ERkZicDAQJXlsbGxAICFCxcWW8ry8OFDuLu7a71Pa2trAFDWnBdV0G5lZaXSXrVqVWWtewErKyuVGvuy6t27N0JCQgDkf7CwsbGBl5dXsaPl1atXV3kdHx+PV199VW09FxcXODo6IiEhQW+xEhEZGhN7IqIKZvLkyThy5AjmzZuH1atXqywrSJo/+ugjNGnSRGP/evXqlWp/jo6OAIDMzEyNywseBFWwXoGiSb0hFC7f0UbRmARBKHZdhUIBS0tLnWMjIipvLMUhIqpgnJ2dERYWhqSkJCxevFhlWa1atQDk13y3adNG5Z+DgwPkcjlsbGxKtb+XX34ZAHDr1i2Ny2/evAkAqF+/finfifHVqlULd+7cUWtPSkpCRkYGPDw8jBAVEZFumNgTEVVAffv2RdOmTXH06FGVdn9/f7i4uGDjxo0qI+wZGRmYMGECpk+fDolEAuD56HVJo9YA0LhxY7i4uCAiIgIpKSkqy2QyGbZu3Qo7Ozu1+nVdFMSmz3Kdkrz++uu4c+cODh06pNJeMPNQ4XsEtCEWi8stdiKioliKQ0RUAYlEInz++ecIDQ1VztUO5N9gO2vWLEyYMAGhoaHo27cvrK2tERERgQcPHuCbb76BhUX+pd/Z2RkAsHr1arz22mvKudqLsrKywueff44JEyagT58+6Nu3Lzw8PPD06VPs2bMHt2/fxuzZs5XbK4tq1aoBAHbv3g1BENCnTx9lvIbw3nvv4cCBA5gwYQIGDhyIOnXq4PTp0zhw4AA6d+6M9u3bl2p7zs7OOHv2LLZv34527dqhZs2aBoqciEgdR+yJiCooX19f5QOdCuvSpQt++uknuLm5Yfny5ViyZAns7e3xww8/oGfPnsr1evTogTZt2mDnzp345ptvStxXcHAwNm/ejCZNmmD79u347LPPsGHDBtSqVQvr169H37599fKevL29MWTIEFy5cgVz587FgwcP9LLd4lStWhXbtm1DSEgI9u3bh/nz5yMmJgZTpkxRK3PSRlhYGPLy8jB79mycPXtW/wETEZVAJLzoO1giIiIiIjJ5HLEnIiIiIjIDTOyJiIiIiMwAE3siIiIiIjPAxJ6IiIiIyAwwsSciIiIiMgNM7ImIiIiIzAATeyIiIiIiM8DEnoiIiIjIDDCxJyIiIiIyA/8HLulY23+DjcgAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,1, figsize=[12,4])\n",
"S['resid2'] = 100*S['resid'] / S['value']\n",
"S['resid2'].unstack(level=0).plot(ax=ax)\n",
"ax.set(title='Bellman Equation Residual',\n",
" xlabel='Net Unit Profit', \n",
" ylabel='Percent Residual')\n",
"\n",
"ax.legend(fontsize=\"x-small\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### SIMULATION"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"Collapsed": "false"
},
"outputs": [],
"source": [
"T = 50\n",
"nrep = 10000\n",
"sinit = np.full(nrep, pbar)\n",
"iinit = 0\n",
"data = model.simulate(T,sinit,iinit, seed=945)\n",
"data['age'] = data.i.cat.codes + 1"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Print Ergodic Moments"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Ergodic Means\n",
"\tPrice = 1.00\n",
"\tAge = 2.01\n",
"\n",
"Ergodic Standard Deviations\n",
"\tPrice = 0.17\n",
"\tAge = 0.83\n"
]
}
],
"source": [
"frm = '\\t{:<10s} = {:5.2f}'\n",
"\n",
"print('\\nErgodic Means')\n",
"print(frm.format('Price', data['unit profit'].mean()))\n",
"print(frm.format('Age', data['age'].mean()))\n",
"print('\\nErgodic Standard Deviations')\n",
"print(frm.format('Price', data['unit profit'].std()))\n",
"print(frm.format('Age', data['age'].std()))"
]
},
{
"cell_type": "markdown",
"metadata": {
"Collapsed": "false"
},
"source": [
"### Plot Simulated and Expected Continuous State Path"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"Collapsed": "false"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" time | \n",
" _rep | \n",
" i | \n",
" unit profit | \n",
" j* | \n",
" age | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0 | \n",
" 0 | \n",
" a=1 | \n",
" 1.0 | \n",
" keep | \n",
" 1 | \n",
"
\n",
" \n",
" | 1 | \n",
" 0 | \n",
" 1 | \n",
" a=1 | \n",
" 1.0 | \n",
" keep | \n",
" 1 | \n",
"
\n",
" \n",
" | 2 | \n",
" 0 | \n",
" 2 | \n",
" a=1 | \n",
" 1.0 | \n",
" keep | \n",
" 1 | \n",
"
\n",
" \n",
" | 3 | \n",
" 0 | \n",
" 3 | \n",
" a=1 | \n",
" 1.0 | \n",
" keep | \n",
" 1 | \n",
"
\n",
" \n",
" | 4 | \n",
" 0 | \n",
" 4 | \n",
" a=1 | \n",
" 1.0 | \n",
" keep | \n",
" 1 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" i | \n",
" unit profit | \n",
" j* | \n",
" age | \n",
"
\n",
" \n",
" | _rep | \n",
" 0 | \n",
" 1 | \n",
" 2 | \n",
" 0 | \n",
" 1 | \n",
" 2 | \n",
" 0 | \n",
" 1 | \n",
" 2 | \n",
" 0 | \n",
" 1 | \n",
" 2 | \n",
"
\n",
" \n",
" | time | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" a=1 | \n",
" a=1 | \n",
" a=1 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" keep | \n",
" keep | \n",
" keep | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" | 1 | \n",
" a=2 | \n",
" a=2 | \n",
" a=2 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" keep | \n",
" keep | \n",
" keep | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
"
\n",
" \n",
" | 2 | \n",
" a=3 | \n",
" a=3 | \n",
" a=3 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" 1.000000 | \n",
" replace | \n",
" replace | \n",
" replace | \n",
" 3 | \n",
" 3 | \n",
" 3 | \n",
"
\n",
" \n",
" | 3 | \n",
" a=1 | \n",
" a=1 | \n",
" a=1 | \n",
" 1.203344 | \n",
" 0.796656 | \n",
" 1.203344 | \n",
" keep | \n",
" keep | \n",
" keep | \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
"
\n",
" \n",
" | 4 | \n",
" a=2 | \n",
" a=2 | \n",
" a=2 | \n",
" 1.101672 | \n",
" 0.898328 | \n",
" 1.101672 | \n",
" keep | \n",
" keep | \n",
" keep | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
"
\n",
" \n",
"
\n",
"