{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"#This code replicates Model SIM (from Godley and Lavoie 2007)\n",
"#Created by Marco Veronese Passarella\n",
"#April 24th, 2020"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%html\n",
""
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#Import packages\n",
"%matplotlib inline\n",
"from ipywidgets import interactive \n",
"import matplotlib.pyplot as plt \n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#Set number of periods\n",
"nperiods=100\n",
"\n",
"#Define variables as matrices of zeros (or ones)\n",
"y = np.zeros([6,nperiods],dtype = float)\n",
"cons = np.zeros([6,nperiods],dtype = float)\n",
"gov = np.zeros([6,nperiods],dtype = float)\n",
"tax = np.zeros([6,nperiods],dtype = float)\n",
"yd = np.zeros([6,nperiods],dtype = float)\n",
"h_h = np.zeros([6,nperiods],dtype = float)\n",
"h_s = np.zeros([6,nperiods],dtype = float)\n",
"n = np.zeros([6,nperiods],dtype = float)\n",
"w = np.ones([6,nperiods],dtype = float)\n",
"sav = np.zeros([6,nperiods],dtype = float)\n",
"\n",
"#Identify coefficients\n",
"alpha1=0.6*np.ones([6,nperiods],dtype = float)\n",
"alpha2=0.4\n",
"theta=0.2"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#Create model using three loops (time, iterations and scenarios)\n",
"i = 1 #Time index\n",
"z = 1 #Iteration index\n",
"j = 0 #Scenario index\n",
"\n",
"for j in range(0,6): \n",
" for i in range(1,nperiods):\n",
" for z in range (1,20):\n",
" \n",
" if i > 70 and j == 1: alpha1[j,i]=0.70\n",
" if i > 70 and j == 2: alpha1[j,i]=0.75\n",
" if i > 70 and j == 3: alpha1[j,i]=0.80\n",
" if i > 70 and j == 4: alpha1[j,i]=0.85 \n",
" if i > 70 and j == 5: alpha1[j,i]=0.90 \n",
" \n",
" y[j,i]=cons[j,i]+gov[j,i] #Output (income) \n",
" cons[j,i]=alpha1[j,i]*yd[j,i] + alpha2*h_h[j,i-1] #Consumption\n",
" n[j,i]=y[j,i]/w[j,i] #Employment\n",
" tax[j,i]=theta*w[j,i]*n[j,i] #Taxes\n",
" yd[j,i]=w[j,i]*n[j,i]-tax[j,i] #Disposable income\n",
" h_s[j,i]=h_s[j,i-1]+gov[j,i]-tax[j,i] #Supply of money\n",
" h_h[j,i]=h_h[j,i-1]+yd[j,i]-cons[j,i] #Demand for money\n",
" sav[j,i]=h_h[j,i]-h_h[j,i-1] #Saving (not included in original model)\n",
" if i > 5: gov[j,i] = 20 #Government spending (trigger)\n",
" z += 1\n",
" i += 1\n",
" j += 1"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"#Create auxiliary variable (time)\n",
"time = np.linspace(0, nperiods, nperiods)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[0, 50, 0, 110]"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xV9f348dc7iywgEECBAGHJEGSqKFhnRVBx1EGL4KyjrmqX2n7r+KrVn7bqt6221t3WMhQHaFEcuNkEkSWbhJVFCNnJve/fH+ckuYTMe29yk9z38/E4j3vP+pz3Pbn3nc/5nHM+R1QVY4wx4SEi1AEYY4xpOZb0jTEmjFjSN8aYMGJJ3xhjwoglfWOMCSOW9I0xJoyEfdIXkb4iUiAikaGOpb0RkWtE5MtGLvuKiDwcaDnBIiIPiMi/3Pf2HTHtRtgkfRHZKSLF7o+3cuilqrtVNVFVPUHYxggR+UBEskUk4BsgROR4EflQRA6KSJ6IrBKRqe68M0Qkw2fZJSKiIjKqRhlvu9PPCDSecBXM70hzcv/Og1pree1JKCoiwRI2Sd91ofvjrRz2Brn8cmAucH2QylsALAaOAXoAdwD59Sz/PTCrckREkoEJQFaQ4jHGtHHhlvSPIiKpbo0myh3vLyKfi8hhEflIRP5aeZjfEFXdrKovAuuDEFc3oD/wD1Utc4evVLW+2sW/gSt9miF+DLwFlNWznVdE5FkR+a979POViBwrIk+7RxibRGSMz/LD3KOKPBFZLyLTfOYli8i7IpIvIsuBgTW2NVREFotIrohsFpErmrZL5M8icsiN6WyfGdeKyEb3b7ZdRG7ymddNRBa68eaKyBciEuHO6yUib4pIlojsEJE76thwze/IEhH5X3dfHXaPxrr5LD9BRL52t7m2KUdZDezfJSJyg894VW1TRD53J691/45XVh4Nish97tHnThGZ4W95dcT7U599v0FExjbiczT1O7dTRO51yz8oIi+LSGyNGLa6f993RaSXzzwVkZtFZIu77l9FRHzmX+fGf1Cco/R+Da0rIsOAvwGnuPHnuctPdWM8LCJ7ROSXDf29Q0JVw2IAdgLn1DI9FVAgyh3/BngSiAEm4dSs/9XEbQ1ydm1A8QqwBVgIXAwcU2P+GUCGz/gS4AbgQ2CKO205cAqQAZxRx3ZeAbKBcUAs8AmwA+eIIRJ4GPjUXTYa2Arc5+6fs4DDwBB3/mycI50EYASwB/jSnZcApAPXAlHAWHe7x/vE8XAdMV4DVAB3uTFcCRwCurrzz8f5ByPA6UARMNad9wecH2i0O5zmLhcBrAJ+736WAcB2YLK73gOVf/daviNLgG3AcUCcO/6YO683kANMdbfxQ3e8uzv/HmBhHZ+zof27BLihxn750mdcgUE1viMVwJ+ADu6+KfS3vFrivdz9G5/o7tNBQL9GfI5XaOR3zue3+x3QB+gKfIX7XXHLzsb5PnUA/gx8XuMzLASSgL44R73nufMuduMchvOd/B3wdSPXPWJfudP2Aae577vgfgdb2xBuNf233ZpHnoi8XXOmiPTF+QL/Xp2a9ZfAuy0eJe5/DDgT5wv/R2CfOEcggxtY9TVglogMAZJU9ZtGbO4tVV2lqiU4RwYlqvqaOm3Yc4DKWtcEIBEnwZWp6ic4P4ofi3N08SOcfVeoqt8Br/ps4wJgp6q+rKoVqroaeBO4rBHxAWQCT6tquarOATbjJHtU9T1V3aaOz3D+8Z3mrlcO9AT6uet+4e7bE3ES8UPuZ9kO/AOY3sh4XlbV71W1GOcf3Wh3+lXA+6r6vqp6VXUxsBLnnwCq+piqXlBHmXXu30bGVJf/UdVSd9+8BzTlCKs+NwD/T1VXuPt+q6ruonGfo7HfuUp/UdV0Vc0FHvEpawbwkqquVtVS4F6cGniqz7qPqWqequ4GPqX6b3UT8AdV3aiqFcCjwGjf2n4969amHBguIp1U9aD7HW91wi3pX6yqSe5wcS3zewG5qlrkMy29OQJxD7krTyj/rbZlVDVDVW9T1YE4NahCnKRen/k4tZ/bgX82MpwDPu+LaxlPdN/3AtJV1eszfxdO7bY7Tm0pvca8Sv2Ak33+6ebh/GCPbWSMe9xk7Vt2LwARmSIiS93D+zycBFvZ3PIETm3uQ3Gafu7xiadXjXjuwzl/0hj7fd4XUb2P+gGX1yh3Es4/nobUt3/9dVBVC2uU16uuhZuoD84RT02N+RyN/c5Vqvm9qvwMvfD5nqlqAc6Rle+26vtbPePzd8rFOWJpzLq1+RHOd2+XiHwmIqfUs2zIRIU6gFZmH9BVROJ9En+f5tiQqj6KU7No7PLpIvJX4D8NLFckIv8FbqFGm3oQ7AX6iEiEzw+6L84J5CycpoQ+wCafeZXSgc9U9Yd+bru3iIhP4u8LvCsiHXCOGGYB76hquXsUJwCqehj4BfALETke+FREVrjx7FDVho6cmiod+Keq/tSPdevbv+D804/3Wb4x/zC7iEiCT+Lvi9NU4m95vtKp/TvW0Ofwh+/vsK+7jcpt+bbDJwDJOM1ODUkHHlHVf/sRz1FX56nqCuAiEYkGbsM5AmyW/BGIcKvp18s9NF0JPCAiMe5/6gsbu757kicWpx0TEYl1k1KTiUgXEXlQRAaJSIQ4JwqvA5Y2YvX7gNNVdac/267HMpxE8WsRiRbnBOWFwGz3sHw+zr6LF5HhwNU+6y4EjhORme660SJyontSrDF6AHe4612O0w77Ps6+7oD7T0dEpgDnVq4kIhe4+1Bwzs943GE5kC8ivxGROBGJFOeS2xP93TmufwEXishkt8xYcU6opjRi3Tr3rzs/DbjU3b+DOPoqsQM45yZqetD9Pp+G08w2L8DyKr0A/FJExrnf/UFu00hDn8Mft4pIioh0xfl+z3Gnvw5cKyKj3d/ao8CyRn73/wbc61YGEJHO7nerMQ4AKSJS+VuPEZEZItJZVcup/q61Opb0jzYD5+RnDs4JpTlAaeVMtznmtDrW7YdzaFp59U4xTtuzP8pwTiB+hPMF+s6N45qGVlTVvVr/VT5+UdUyYBowBefk2bPALFWtrNnfhnP4ux/nZN3LPusexknG03FqZ/uBx3ESdmMsAwa7230EuExVc9xy78CpVR0EfsKR52EG4+zDApyT9M+q6hL3n9SFOG20O9xyXwA6NzKeWqlqOnARTmLKwqlN/gr3t+Y26/23jnUb2r9P4XwvDuCcL6lZQ30AeNVtrqhst9+Ps1/2usvfHGB5vvHOw/lbvI5zovZtnJPrDX0Of7yOc65muzs87MbwMfA/OEd7+3COPBp1XkZV38L5Ds4Wkcrf2JRGxvMJzu98v4hku9NmAjvdsm7GOb/T6siRzaSmJhGZA2xS1ftDHYsxTeHWsP+lqo05ymi1RGQnzlVGH4U6lvbAavo1uE0OA90mlfNwam1HXeljjDFtkZ3IPdqxOG3TyTjXt9+iqmtCG5IxxgSHNe8YY0wYseYdY4wJI62ieadbt26ampoa6jCMMaZNWbVqVbaqdm/KOq0i6aemprJy5cpQh2GMMW2KiOxqeKkjWfOOMcaEEUv6xhgTRizpG2NMGGkVbfq1KS8vJyMjg5KSklCH0iixsbGkpKQQHR0d6lCMMaZOrTbpZ2Rk0LFjR1JTU5HqB920SqpKTk4OGRkZ9O/fP9ThGGNMnVpt805JSQnJycmtPuEDiAjJyclt5qjEGBO+Wm3SB9pEwq/UlmI1xoSvVp30jTHGBJclfWOMCSOW9I0xJoxY0q/Hxo0b6d+/P16v85hPr9fLueeey2uvNfRscmOMaZ0s6ddj2LBhDB06lIULFwJw3333MWTIEGbNmhXiyIwxxj+t9jr91uKuu+7iqaeeory8nK+++opPPvkk1CEZY4zfrKbfgHPPPZeMjAzuvfde5s6da3fcGmPatDZT00+9572gl7nzsfMbtdypp57KmDFj6NmzZ9BjMMaYltRg0heRl4ALgExVHeFO6wrMAVKBncAVqnpQnDuUngGmAkXANaq6OhiBNjZBN4cNGzZw7bXXhmz7xhgTLI1p3nkFOK/GtHuAj1V1MPCxOw4wBRjsDjcCzwUnzNBav349I0aMCHUYxhgTsAaTvqp+DuTWmHwR8Kr7/lXgYp/pr6ljKZAkIm26TSQ9PZ2kpCQSExNDHYoxxgTM3zb9Y1R1H4Cq7hORHu703kC6z3IZ7rR9NQsQkRtxjgbo27evn2E0vz59+rB9+/ZQh2FMm+f1KhVexeNVKrxe91WrXn3ne9V5rRpUUVU8XqrmVy6jylHLVM5XpdZXb9W4om5sldOA6uV93jvTneVUQXHe47Ocou5r9ThV40fOA3e+z3jlNpx5PmVUveeo5Zsq2Cdya+t1rNbwVPV54HmA8ePHB/ARjDE1eb1KSYWH4jIPRWUeSis8lJR7KSn3UFxe/b60wktphYfSci+lFV7KKsfd9+Ue57XM4zPu8VJeoZR7nfEKj1LmvlZ4vJR7ndcKr1LhcRJzudeLKkRHCpERQlREhPvqjB8xiBDhzosQZ1pEhBApECGV74WICKrnS+VA9XiEMy44y4k7X9xlRKRqnjNdEAGhejnfaSK45Un1PEAipHoa1etR+b7mPHe8UtWyVeO+7+WIcdxtVi7nD3+T/gER6enW8nsCme70DKCPz3IpwF4/t2FMWPF6lYKyCg4VlZNfUs7hkgp3KK9+La2gqNRDYWkFBaUVFJZVUOiOF5U5Cb24zENJhYcOURHERUcSFx1JbEwksVGRxEZHEFs5LTqSDlERdKh8dYf4mEi6xMcQExVBTFQE0ZHOa0ykVI07gxAdGUFURAQxUU4ij6qaJkRVvVYneRNcV/mxjr9J/13gauAx9/Udn+m3ichs4GTgUGUzkDHhRFU5VFxOdkEZuYVl5BaWklNYRm5BGTmFZRwsKiOvqJy84nLyi8vJKyojv6SCuOhIOsVG0Skumo6xUXSMrXx13neKjaZX5zgSOkSR2CGS+JgoEjpEkdAhkoSYKOJiIqsSfYQlWVOLxlyy+R/gDKCbiGQA9+Mk+7kicj2wG7jcXfx9nMs1t+JcsmnXOZp2p6isgr15Jew7VMzevGIO5JeSebiEzPxSMg+XkuUOHaIj6JbYga4JMXRNiCHZfU3pEscJKZ3pEh9D5/hoOsdFkxQXTae4aKIj7X5J07waTPqq+uM6Zp1dy7IK3BpoUMaEUkm5h/TcInblFLE71xkyDhaxx030xWUeenaOpVdSHD07x3Fs5w4M7tGRiQO70aNTB7onxtKjUwdioyND/VGMOUqbuSPXmGDyepW9h4rZmlnA1swCtmUVsC2rkF05hRwsKiclKY6+yfH07eoMpwxMpndSHD07x9I1IcaelGbaLEv6pt07WFjGxn35rN+bz4Z9+WzJPMz2rEI6xkYxqEcig7onMrxnJy48oRf9uiVwbKdYO+lo2i1L+qZdySkoZfXuPNbtOcSGvfls2HuI/JIKhvfsxPBenZgwoCuzTunHwB6JdIq1zvNM+LGk3wjz5s3jySefpLi4mKKiImbOnMn9998f6rDCnserbN5/mFW7D7Jm10FW7T5IbmEZo/skMSoliUvH9uZ/LhhGny7xdiWLMS5L+g149dVX+fOf/8zbb79NSkoKBQUFPPdcu+hSqM3xepWN+/P5emsOX27NZtWug/To1IGxfbtwYv+u3HT6QAb3SLQEb0w9LOnXIz8/n7vvvpsVK1aQkpICQGJiIr/61a9CHFn4SM8t4qut2Xy5NZuvt+XQOS6aiYOS+fFJfXj6ytF0SYgJdYjGtCmW9Ovx1ltvcfLJJzNgwIBQhxI2VJV1ew7x4foDfLhhP7mFZUwc1I0fHNede6cOo3dSXKhDNKZNaztJ/4HOzVDmoXpnr1+/ntGjRwd/u+YI5R4vy7bn8uGG/SzecIDY6EjOPf4Y/nDpCYzpk2TNNcYEURtK+vUn6OaQkJBAcXFxi283HKgqa9LzeHNVBu+t20e/5ATOHX4M/7z+JAZ2T7Tr4I1pJm0n6YfA1KlTmT59OnfddRfHHHMMpaWlvPbaa/z0pz8NdWht1p68Yt5ancH81XtQ4NIxvVl4+yRSusSHOjRjwoIl/XqceOKJPPDAA0yePBmPx0NFRQVXXeVPv3bhrazCy/vr9jF3ZTob9uUzdWRPnrh8FGP7JlmN3pgWZkm/ATNnzmTmzJmhDqNNyiko5fVlu/nn0l0MPiaRGSf34+xhPaxPGmNCyJK+CbpN+/N5+cud/Pe7fUwZ0ZPXrj+Jocd2CnVYxhgs6ZsgUVW+3JrN3z7bxpYDBcyc0I9Pf3kGyYkdQh2aMcaHJX0TsBU7c3nig81kHy7l1jMHceGoXsREWb/wxrRGlvSN377NyOPJD79ne1YBd549mEvG9CbKHgJiTKtmSd802ab9+fzpw+/5NuMQt541iCvHj7eavTFthCV902h5RWU8vmgzizcc4ObTB/B/Px5jV+IY08ZY0jcNUlXeWJXB44s2c/7IY/n4F6fTOc76ojemLbKkb+q1ef9hfvf2OkorvLx8zYmMTGmGPpCMMS3Gkn4jfPDBBzz88MMUFxdTWlrKuHHjePLJJ+nWrVuoQ2s2haUV/N/HW5i3KoO7fngcPzmprz1C0Jh2wJJ+A+bNm8fDDz/MW2+9xYABA/B4PDz22GOUlJSEOrRms3R7Dr+Yu5aT+3flg5//gO4d7Vp7Y9oLS/r1KCws5Pbbb+fjjz+u6lM/MjKS3/72tyGOrHlUeLw88/EWZq9I5//96ATOHNoj1CEZY4LMkn493n//fUaNGsXxxx8f6lCaXXpuEXfOXkNChyjeu2MSPTrGhjokY0wzaDNJ/9m0Z3lubfWzaWdfMBuA6QunV027ZdQt/Gz0zzhr7llkFWcBMKzrMOZeOJcHvn6AN7e8WbXsx5d/TI/4+muy69evZ8SIEVXjd9xxB5988gmJiYksXbo0KJ+rNXh37V4efHc9N58+kOsn9beHlhjTjomqhjoGxo8frytXrjxi2saNGxk2bFiIInI8/vjjZGZm8sc//rFq2qJFi3jhhRd44403jlq+NcTcFIWlFdz/7npW7TrI/00fY1fmGNPGiMgqVR3flHXsNsp6TJ48mfnz57N3717AuV598eLFjB07NsSRBW53ThHT/vIlAiy8fZIlfGPCRJtp3gmF0aNH8/DDD3PeeecRGRlJdHQ048ePb/P966/alcvN/1rNHWcNYuYpqaEOxxjTgizpN2DGjBnMmDEj1GEEzTtpe3hwwQb+eMUozhxiV+cYE24s6YcJVeUvn2xl9op0/n3DyQzraQ81MSYcWdIPA6UVHu6dv44tBwp462en0qOTXY5pTLgK6ESuiNwlIutF5DsR+Y+IxIpIfxFZJiJbRGSOiMT4W35ruLKosVprrHlFZcx8cTmFpRXMuWmCJXxjwpzfSV9EegN3AONVdQQQCUwHHgeeUtXBwEHgen/Kj42NJScnp9UmU1+qSk5ODrGxrSuh5haWMf35pYzs3ZlnZ4wjPsYO7IwJd4FmgSggTkTKgXhgH3AW8BN3/qvAA8Bzta5dj5SUFDIyMsjKygowxJYRGxtLSkpKqMOocrCwjBkvLOOMIT34zXlDELEbroypogrqrX6lcry2aeoM1HitWqaWeUe9cvR032lHzaeeZXzm+cHvpK+qe0TkSWA3UAx8CKwC8lS1wl0sA+hd2/oiciNwI0Dfvn2Pmh8dHU3//v39DS+s5RWVcdWLyzhtcDdL+KZhXi94SqGiBCrK3NdSZ5qnzJnm8RkqSsFTDt5yd1qF8+otd957y33mV4DXneatAK/Hfa0cvEeOq7d6GfVUL6Oe6nnqcV+1xvTKpF05z1v3AICARIC4r1XjNae5A76vEbVMq/lKA/PFJ46a793xI6ZT/f6IeU3jd9IXkS7ARUB/IA+YB0ypZdFa/yWp6vPA8+DcketvHOZIh4rKmfnick4ZkMy9U4Zawm8vKsqgNN8ZygqhtADKCqD0sDNeVuC8lhdBWRGUF0J5sc/7EqgorvHqDp4yiOwAUbEQ1aF6iOwAUTEQ6TNEdYDIaOd9RLT7vnI8qvp9THz1/Iio6qFqPBIk0mdeZI1p7qtE+EyPdJNtJERE+LyvnO6OVybsI6b7DIizfntwc9N/34E075wD7FDVLAARmQ+cCiSJSJRb208B9gawDdMEh4rLmfnSMsanduG35w+zhN/aeD1QlAtF2VCUA8UH3SGv+n1JnjNecshJ8CVuoveUQ2wn6NARYjpCh0SISYSYBHdagjNEJ0B8MkTHOe9j4t338U5Sj4478rUy0dt3JWwEkvR3AxNEJB6needsYCXwKXAZMBu4Gngn0CBNw/JLypn10nLG9u3C7y8Ybgm/pXi9ThI/vB8KDriv+6Eg03lfmO3ML8x2EnlsZ0jo5iTmuK4Q1wXikpyhc29nPLYzxCZBh05uou/kJGn7m5ogCKRNf5mIvAGsBiqANTjNNe8Bs0XkYXfai8EI1NStoLSCq19azgm9O3P/hZbwg6qsCA7uhLzdcCgd8vfAoQw4tAfyMyB/n5OYE4+Fjsc4r4k9oOsA6HsKJHR3k3w3iO/qNDkYE0KttpdN0zger3LDqyvo0TGWx3400hK+P0oLIPt7yNkKuTvg4I7q15JDkNQXkvo5NfHOKdAppfp9x14Q3bou1TXhw59eNu3C7TbuD+9vpLTCy8OXjLCE35CSfDjwHWRtgqzvIXuz81qUDcmDoNtgp4bebyKMngFd+ztJvb2c9DMGS/pt2uzlu/l4UyZv/exUoiMtMVVRddrT93/rDPu+hf3rnDb37kOhx3DofhwMOMN5TepnzS4mbFjSb6O+2ZbDkx9uZs5Np5AU73dPF+1DeQnsS4P05ZC+DDJWOFe79DwBjj0Bhk2DM3/r1OYj7Stvwpv9AtqgndmF3P6fNTx95RgGdk8MdTgtrzgPdn4Ju75yEn3mBuh2HPQ5CY6/BCY/6rTDW3OXMUexpN/GHCou5/pXV/DzcwYzaXC3UIfTMspLIGM5bF/iDFmbnQSfOgl++CD0GuNco26MaZAl/TakwuPlttdXc9rg7lw1oV+ow2leeemw+X34fpFTm+8+1GmDP+dBJ+FHdQh1hMa0SZb025CH39uIiPC789vOw9cbTRUyN8Km92DTAifpH3cejLsWLnvZuXnJGBMwS/ptxMJv9/LZ91m8c9tEotrTlTr7v4NvZ8PGhU43BUPPh3MfcW5sspOuxgSd/aragL15xdz/znpeuuZEOsVGhzqcwBVkwrp5kPYfp7+ZUVfCFa/BsSPt5KsxzcySfivn9Sq/mLuW6yb1Z1SfNtzEUVEKm/8La/8Du76BoVNh8iOQeprd/GRMC7Kk38r944vteLzKzacPDHUo/jm8H1a8AKtegR7DYNRP4EcvOr1EGmNanCX9Vuy7PYd4/vPtvHPbRCIj2lizx941sPQ5+P4DGHk5XPtfp5sDY0xIWdJvpYrLPNw5ew2/v3A4KV3iQx1O43g9sGmhk+zz0uHkG2HK4053wcaYVsGSfiv16PsbGdG7MxeNrvVpk62L1wvr58OSPzgJ/pRbYeiFdvWNMa2Q/SpboY83HuCTTZm8f+dpoQ6lfqrODVSfPOJ0Lzz1CRhwpl2BY0wrZkm/lck6XMo989fx15+MpXNcK708UxW2fwqfPOxclXPW72DIFEv2xrQBlvRbEVXl3vnfcsX4FE7q3zXU4dRubxp8cJ9zrf2Z98LwS+ySS2PaEEv6rcjiDQfYkV3IszPGhTqUo5Uccmr26992avajZ1ibvTFtkFXRWomScg8PLdzAQxeNICaqFf1ZVOHbefCXk5ymnFuXwbirLeEb00bZL7eVeHbJNkb1SWLioFbUXXL2FnjvF1CUC1f+C/qcGOqIjDEBsqTfCuzKKeSf3+xsPVfrVJTB5084d9L+4Fdw0o1WszemnbBfcivw4IIN3PiDgfTsHBfqUCBnG7xxHXTsCbd8BZ16hToiY0wQtaLG4/D00YYD7Mwp5PpJ/UMbiCqs+Te8+EMYcxX8+D+W8I1ph6ymH0Il5R4eXLieRy8ZGdqTtyWHYOFdcGADXL0Ajjk+dLEYY5qV1fRD6Lkl2xjZuzOnDe4euiDSV8DfToPYJLjxU0v4xrRzVtMPkd05Rbz2zU7euyNEJ29V4aun4Zu/wgVPw7ALQhOHMaZFWU0/RB5csJ4bThtAr6QQnLytKIW3b3FutLrxM0v4Ddi7dy+dOnViz549VdM+//xzBg4ceMS0plqzZg3Tpk1j5MiRDBkyhEmTJvHRRx8FFOumTZuYPn0648ePrxoee+yxgMo07Ysl/RD4eKNz5+0Np4Xg5G1RLvzzEig9DNe+D53bQC+eIdarVy9mzJjBM888A8DmzZuZNWsW8+fPp3dv//bfl19+ybRp07j99ttZt24dmzdv5rnnnuPgwYN1rpOamlpvmcuXL+ecc85h1qxZrFy5kpUrV7JgwQIyMzP9itG0U6oa8mHcuHEaLjwer/7wT0v0ow37W37j2VtVnxmj+sFvVT2elt9+G7Zjxw5NTk7WrVu36nHHHaeLFi3yu6yysjJNTU3VefPmNWm9fv361TmvoqJChwwZoi+//LLfcZm2B1ipTcy31qbfwt5bt4/4mCjOGtqjZTe862uYe7XTSdr461p22+1Aamoq559/PuPGjeOJJ55g8uTJR8y/7LLL2Lp1a63rfvPNN8TFVTfjLVq0CFXlRz/6UdDiW7JkCfn5+cycOTNoZZr2KaCkLyJJwAvACECB64DNwBwgFdgJXKGqdR+zhhGPV3n6o+/5/YXHIy3ZDfHaOU7PmJc+D4PObrnttjNnn302u3fv5qc//elR8954441Gl5OWlsa4ceMa9R2YNm0au3fvBpxzC6NHjwYgKiqKlStXHlHmqFGjiIyMBGDp0qXcfPPNVf8IHnzwwUbHZ9q3QGv6zwCLVPUyEYkB4oH7gI9V9TERuQe4B/hNgNtpFxZ+u5ek+Bh+MLgF+9f58ilY8ZJ7/f3wlttuO7R27VpOPLH2/oeaUtPv2LEjzpF5w959992q96mpqaSlpdW6XHz8kY/UnDBhAmlpaUyZMoWhQ4c2alsmTDS1PahyADoBOwCpMX0z0NN930Cel20AABloSURBVBPY3FBZ4dCmX17h0TOf+FS/+D6r5Tb6xZ+cNvxDe1tum+3YmWeeqa+//nrA5WzevFm7dOmiK1asqJqWlpam77//fr3r1demv337dk1KStJvvvmmalpmZqYmJSXptm3bAo7ZtE60cJv+ACALeFlERgGrgDuBY1R1n/sPZZ+I1Np4LSI3AjcC9O3bN4Aw2oZ30vbSLbEDEwclt8wGv/o/WP0aXPMedOrZMtts59LS0hgzZkzA5Rx33HHMmTOH22+/ncLCQkpKSujXr19Al1b279+fN954g5///OcUFBSQmJhIXFwcDz30EAMGDAg4ZtN+iDbyMPOoFUXGA0uBiaq6TESeAfKB21U1yWe5g6rapb6yxo8fr77tk+1NhcfL2X/6jMcuPYFTBrZA0v/mr7D8ebjGLsk0pj0TkVWqOr4p6wRynX4GkKGqy9zxN4CxwAER6ekG1BMI+4uE56/ZQ8/OsS2T8Jc+5yT8qxdawjfGHMXvpK+q+4F0ERniTjob2AC8C1ztTrsaeCegCNu4co+XP3+yhbvOOa75N7bs77D0WeekbVKf5t+eMabNCfTqnduBf7tX7mwHrsX5RzJXRK4HdgOXB7iNNu3NVRn065rAyQOauZa//B/w9V/gmoWQ1P7PkRhj/BNQ0lfVNKC29iS7GBwoq/Dy50+28n8/Ht28G1o7B756xqnhd+nXvNsyxrRpdkduM5q7Mp2BPRIZ169r821k9zLnxqurF0DXED+IxRjT6lmHa82ktMLDs59u5a5zBjffRvJ2w9yZcPFzduOVMaZRLOk3k/mr93DcsR0Z07feq1X9V3oYXp8OE++E485tnm0YY9odS/rNQFV55aud3DCpmW6K8XrgzRsgZRxM+FnzbMMY0y5Zm34z+GZ7Dh7V5rv79qMHoLQArvgntGTHbcaYNs+SfjN49eudXH1qavP0pLnmX7BxAfz0E4iKCX75xph2zZp3giw9t4hlO3K5dEwz3A2762tYfD/8ZC7EN+MVQaZWzfF4w2Bp6JGOn3/+OR07dmT06NGMHj2aMWPG8MQTT+DxeEIYdWA+//zzo3oQXbZsWVj05RUIS/pB9q+lu7hsbAoJHYJ8EJW/F+Zd4/SJ370F7u41R/Dn8YYtqaFHOq5evZpLLrmEtLQ00tLS+PDDD3n33Xd55JFHghZDQ49zDLbVq1cf1dX18uXLGTt2bIvG0dZY0g+i4jIPc1emM+uU1OAW7PXCWzfD+OvtISghUF5ezsyZM3nqqaf44Q9/WDV95MiRXH65/zec79u3jyuvvJIxY8YwevRovvvuOyZOnIjX6/WrvN/85je89NJLbNu2jWnTpvH3v/+dUaNGAU6CrHwAC0D37t359a9/zbx58/yOv6mC/XlXr17NSSeddMS0FStWWNJvgCX9IHo7bQ/j+nWhb3J8wws3xdJnoaIETvtFcMs1jdLYxxtedtllVc0nNYfi4uIjllVVLr74YqZNm8aaNWs46aSTmDZtGr/73e+IiKj/Z/n222+zf//+o6b7PtLxl7/85RGPdFy9evVR3UInJCSQk5PT0McPikA+L9T+mVevXs2jjz5Kampq1TB79mxL+g1pagf8zTG0h4eoeL1enfzUZ8F/SMq+b1Uf76+asz245ZpGe+ihh/TSSy8Naplff/21nnDCCVXj999/v55++umqqpqXl6fXXnutpqSk1Lpuv379dPHixbXOe/XVV/WMM844YlpRUZFGRUVpbm7uEdOfe+45Peuss+qNs6FYLrzwQh01apSOGjVKo6Ojq97X/E3X93lVVb/44gu97rrrdMaMGXrrrbc2+JmLioo0JiZGCwoKjpgWGRmpe/bsqZp25ZVXanl5eb2fsSEN7YNQwh6MHjpLt+dS4Q3yZZrlxfDmT+HcR6yLhRBq7OMNm/LIxOXLlzNhwoSq8bS0NK67znlgfefOnXnppZc455xzai1r586ddcZQ2yMd165dS69evejSpfpGQY/Hw7PPPsvdd99d72dqKJbGPs6xvs8LMGnSJCZNmgTARRddVPUgmEo1P/PatWtJTU0lISGhatq6devo1q0bvXr1ApwK7Y4dO3jiiSdYtGgRzz33HMOHN/3O9Yb2QVtjST9IXvl6R/Av0/zoAegxFEZND16ZpsmmTp3KQw89xMqVKxk/3ulfcO3atezdu5cpU6ZULdeUh6NHR0dXPfD8nXfeYdGiRVxwwQUBx7pmzZqjHtxesz0/KyuLu+++m759+zJr1iy2bdvGTTfddMQ6kydP5le/+lXA8VRq7Od97733GDZs2BEJvzarV6+uOl9Rae3atUc07WzcuJEJEyZw77330q1bN3bs2FFr0m+Jz9+aWNIPgoyDzmWaf7oiiL1pbv0INi6EW760G7BCrDkeb3jFFVcwZ84chg4dSp8+ffj3v//NnXfeyYknnnhUMmuK2h7puGbNGr7++mvGjh1LREQEsbGxXHnlldxyyy1EREQwcODAZr/0tDGf95VXXmHnzp2N2q9r1qxpMOl/9dVXVf+U16xZw/nnn19rWS3x+VuVprYHNcfQ1tv0H31/gz60YH3wCizIUn1yiOr2z4JXpmlTfvazn2nv3r31pptu0u3bQ3s+pyViWbBgQdU2brrpJs3MzAy4zGuvvVazs7NVVfWCCy5QVdX9+/fr008/3eSyWtPfwxd+tOn7/YzcYGrLz8gtLvMw8fFPeOtnp9IvOaHhFRqiCrNnQPJAOPd/Ay/PGFNl4cKFREZGHtEs15b584xca94J0DtpexjbNyk4CR9g9atwaDdc/nJwyjPGVAnGeZO2zq7TD4Cq8orbz05QHN4PHz0Il74AUR2CU6YxxviwpB+ANel5lFZ4mTSoW3AK/PB/YNzVzhU7xhjTDKx5JwBvrd7DpWN6B+cyzZ1fOR2q3bY88LKMMaYOVtP3U1mFl/fW7ePiYPSm6SmH938Jkx+BmCCdGzDGmFpY0vfTZ99nMbB7An26BqGfneX/gMQeMPyiwMsyxph6WPOOn95ak8ElY1ICL+jwfvjiSbh2kd2EZYxpdlbT98Oh4nK++D6b80f2DLywxb+HMVdZH/lh7qOPPqrqkXPq1KmhDse0Y1bT98P76/YxaXA3OsdHB1bQrq9h55dwq528DXfnnHNOnZ2VGRNMVtP3w1tr9gR+AtdTAe/9Es59GDrU37mUMcYEiyX9JkrPLWLLgcOcOaRHYAWteAESkuH4S4ITmGl2mzZtYvr06YwfP75qCKTTteYu15jaWNJvonfS9nD+CT2JiQpg1x0+AJ//P5j6pJ28bSOWL1/OOeecw6xZs1i5ciUrV65kwYIFZGZmtspyjamLdbjWBKrK2X/6jCcuG8W4fl0aXqEu794OHTo51+WbVs/j8XD88cdzzz33cM0117T6ck34sA7Xmtm6PYfweJWxfZP8LyR3O2xcALevDl5gplktWbKE/Px8Zs6c2eCyTXl6VlPKNSZYLOk3wfzVe7h4dIDdLnz+JJx0I8R3DV5gplmlpaUxatQoIiMjAVi6dCk333xzVcJ+8MEHq5ZtytOzmlKuMcEScNIXkUhgJbBHVS8Qkf7AbKArsBqYqaplgW4n1Mo9XhZ+u5c3bj7V/0JytsHm/8Ida4IXmGl28fFH3nU9YcIE0tLSmDJlCkOHHtk5XlNq+k0p15igaepTV2oOwN3A68BCd3wuMN19/zfglobKaAtPzvp443695K9fBlbImzeqfvpYcAIyLWb79u2alJSk33zzTdW0zMxMTUpK0m3btrW6ck34wI8nZwVU0xeRFOB84BHgbnHaPc4CfuIu8irwAPBcINtpDeav3sMlYwPodiF7C2xdbLX8Nqh///688cYb/PznP6egoIDExETi4uJ46KGHGDBgQKsr15j6BNq88zTwa6CjO54M5KlqhTueAdR6F5OI3AjcCNC3b98Aw2heh0vK+WxzFv970Qj/C/nscZhwC8R2Dl5gpsWcffbZnH322W2mXGPq4vfF5iJyAZCpqqt8J9eyaK3XhKrq86o6XlXHd+/e3d8wWsR/v9vPhIHJdEmI8a+AzE2w7VM4+ebgBmaMMU0USE1/IjBNRKYCsUAnnJp/kohEubX9FGBv4GGG1jtpe7jq5H7+F/DZ43DKrdChY8PLGmNMM/K7pq+q96pqiqqmAtOBT1R1BvApcJm72NXAOwFHGWLfHyhgrL83Yx3YADu/cC7TNMaYEGuObhh+g3NSdytOG/+LzbCNFuP1KgcLy+jqb9POZ4/BqXdYp2rGmFYhKDdnqeoSYIn7fjtwUjDKbQ0OFpWRGBtFdKQf/x/3r4PdS+HiNn/xkjGmnbAO1xqQXVBGt8QO/q285DGYeKc999YY02pY0m9ATkEp3RL9aNrZ9y1krITx1wU/KGOM8ZMl/QZkFZT6V9Nf+ixMuBmi4xpe1hhjWogl/Qb41bxTkAWb34exVzdPUMYY4ydL+g3I9qd5Z9UrMPwi60nTGNPqWNJvQPbhJjbveMph5Ytw0k3NF5QxxvjJkn4Dcgqb2Lyz8V3oOhCODaCfHmOMaSaW9BuQXVBKt45NSPrL/g4nWy3fGNM6WdJvgNO808g2/b1rIH8vDJnavEEZY4yfLOnXQ1XJbkrzzrLn4cTrIdKeQmmMaZ0s6dcjv6SCmMgIYqMjG164IAs2v2eXaRpjWjVL+vVo0t24q1+BYdPsMk1jTKtmSb8ejb4xy1MOK16yE7jGmFbPkn49shvbBcPGBdC1Pxw7svmDMsaYAFjSr4dzuWYjmneWP2+1fGNMm2BJvx7Zh0tJTmigpr83DfLSYcj5LROUMcYEwJJ+PbILyxq+MWu5XaZpjGk7LOnXI/twKd3ru3qnKBc2LbTLNI0xbYYl/Xo0eCJ3/XwYeDYkJLdcUMYYEwBL+vVo8JLNtXNg1PSWC8gYYwJkSb8e2QWlJNfVvJO7HQ7ugIFntWxQxhgTAEv6dSgqq8DjVRI71HGC9tu5cPylEBndsoEZY0wALOnXIcdt2hGRo2eqwtrZMOrKlg/MGGMCYEm/Dln19aOfsQIioqDX2JYNyhhjAmRJvw71Xq757Rynll/bUYAxxrRilvTrkF1QVvvduBVlsP4tGHlFywdljDEBsqRfh5y6+t3Zuhi6DYEu/Vo+KGOMCZAl/TrUeWNWZdOOMca0QZb061DrjVnFebDtUxh+UWiCMsaYAFnSr0NWbTX9DW/DgDMgrksoQjLGmID5nfRFpI+IfCoiG0VkvYjc6U7vKiKLRWSL+9omM2R2bY9KtG4XjDFtXCA1/QrgF6o6DJgA3Coiw4F7gI9VdTDwsTve5uTUbN45uAuyNsGgH4YuKGOMCZDfncCr6j5gn/v+sIhsBHoDFwFnuIu9CiwBfhNQlC2srMJLUVkFneOcLhbKveWsW/FnOO4MyPmOrrFdSe2cyubczRSWFwIQGRHJqO6j2F+4n70Fe6vKOq7LcURIBJtyN1VNOybhGHon9ua77O8o85QBEBsVy/Dk4aTnp5NVnFW17PHdjqekooRteduqpqV0TKFHfA/WZK5BVQHoGNORwV0Gsz1vO3mleVXLju4xmtySXHbn766a1r9zfxJjElmXta5qmn0m+0z2mdrmZ2qqoDz5Q0RSgTHAMuAY9x8CqrpPRHrUsc6NwI0Affv2DUYYQZNTWErXhBgiIpybr0rLS3g6/UNIHgirn2Zir4ncNOom5m6ey5a8LQDER8fzt3P+xor9K5j3/byqsn578m+Jj47n6dVPV007v//5XDn0Sl767iWyi7MB6J3Ymz+c9geWZCxh8a7FVcs+efqTZBZlHrH+jGEzmJw6mb+s+Qvl3nIAjk8+nt+c9BsWbl/IygMrq5Z94dwX2JizkX+s+0fVtNtG38bw5OFHlGmfyT6Tfaa2+ZmaSir/s/ldgEgi8BnwiKrOF5E8VU3ymX9QVett1x8/fryuXLmyvkVa1LqMQ9wz/1veu+M0CsoKWJz2Apd88wrcscbuwjXGtBoiskpVxzdlnYCu3hGRaOBN4N+qOt+dfEBEerrzewKZgWwjFJwulZ32/O8Pfs8bW96EE6zbBWNM2xfI1TsCvAhsVNU/+cx6F6h8fuDVwDv+hxcavlfubM3byqCCPBhxaYijMsaYwAXSpj8RmAmsE5E0d9p9wGPAXBG5HtgNXB5YiC0vu6CM7m5Nf0vWdwwqKYLkwSGOyhhjAhfI1TtfAnW1d5ztb7mtQXZBKcd2igVgVtfRxCV+DRF2H5sxpu2zTFaLbLezNVUl/8C3JB9zQqhDMsaYoLCkX4vsglKSEzqQU5LDzenvwLEjQx2SMcYEhSX9WmQfdu7G3Zq3lYEVihxrNX1jTPtgSb8WOYVO887WnE0MKiqAHsNCHZIxxgRFUO7IbU88XiWvqJyu8TGMiOzEiIjOEBMf6rCMMSYoLOnXkFtYRue4aKIiIxhZUkxkj1GhDskYY4LGmndqqHxilqpy+ro/crj7caEOyRhjgsaSfg1OFwwx7C/cT4x66dirSd1aGGNMq2ZJv4bKfvS3HNzCoNJSOHZEqEMyxpigsaRfQ2XzTlTJIc4q80LiMaEOyRhjgsZO5NaQ5d6Ne6rGcGqnYdazpjGmXbGafg2VN2bd++2zpHcbEOpwjDEmqCzp15BdUErXhEg+KtxB155jQh2OMcYElSX9GnIKS/FG5pDsVRJ625U7xpj2xZJ+DdmHy/B6DjCypMT60DfGtDt2IteHqpJbWMZZcUmcF9ELIm33GGPaF6vp+zhUXE5sdAQLvp/Dju52EtcY0/5Y0vdReY3+vzKXUpo8MNThGGNM0FnS95FdUEZyYiTpFYX07zMp1OEYY0zQWdL3kV1QSkJCFj0rPHToZZdrGmPaHztT6SP7cCmjoqP4SW4kxHYOdTjGGBN0VtP3kV1QRqJ3Dd4eQ0MdijHGNAtL+j6yC0r5uvwzvu3ULdShGGNMs7Ck7yO7oIz93oMM7nliqEMxxphmYUnfR2bBYbKknD79fhDqUIwxplnYiVwfFYczeay0iOiudo2+MaZ9spq+S1VJLt7CxM6DrQ99Y0y7ZUnfVVjmoazLV/yzU0KoQzHGmGZjSd+VU1CKNzaHwd3tmbjGmPbLkr4ru6CU7OhiBqVMDHUoxhjTbCzpuzIP5nNuYSG9+54W6lCMMabZNEvSF5HzRGSziGwVkXuaYxvB5s3cxLWFnYiIiQ91KMYY02yCnvRFJBL4KzAFGA78WESGB3s7wbb2wDv8LblLqMMwxphm1RzX6Z8EbFXV7QAiMhu4CNhQ1wq7sjYy8+/O82h7eaL48eHOLIovYH1MadUyPzvUhfSochYkFFRNm1yUyNCyGJ5Jyq2aNrg8hmmFHXkjMZ9dUeUARCL8PK8rqzuU8GlcYdWyPyroRBdvBC90yiMzsoTT4k4Iyg4wxpjWSlQ1uAWKXAacp6o3uOMzgZNV9bYay90I3OiOjgC+C2ogbVc3IDvUQbQSti+q2b6oZvui2hBV7diUFZqjpl/bnU1H/WdR1eeB5wFEZKWqjm+GWNoc2xfVbF9Us31RzfZFNRFZ2dR1muNEbgbQx2c8BdjbDNsxxhjTRM2R9FcAg0Wkv4jEANOBd5thO8YYY5oo6M07qlohIrcBHwCRwEuqur6B1Z4PdhxtmO2LarYvqtm+qGb7olqT90XQT+QaY4xpveyOXGOMCSOW9I0xJoyEPOm3xS4bgkVEXhKRTBH5zmdaVxFZLCJb3Nd2f5uwiPQRkU9FZKOIrBeRO93p4bgvYkVkuYisdffFg+70/iKyzN0Xc9yLJMKCiESKyBoRWeiOh+W+EJGdIrJORNIqL9X05zcS0qTfVrtsCKJXgPNqTLsH+FhVBwMfu+PtXQXwC1UdBkwAbnW/B+G4L0qBs1R1FDAaOE9EJgCPA0+5++IgcH0IY2xpdwIbfcbDeV+cqaqjfe5TaPJvJNQ1/aouG1S1DKjssiEsqOrnQG6NyRcBr7rvXwUubtGgQkBV96nqavf9YZwfeG/Cc1+oqlb2NRLtDgqcBbzhTg+LfQEgIinA+cAL7rgQpvuiDk3+jYQ66fcG0n3GM9xp4ewYVd0HTjIEeoQ4nhYlIqnAGGAZYbov3OaMNCATWAxsA/JUtcJdJJx+J08Dvwa87ngy4bsvFPhQRFa53diAH7+RUD8YvVFdNpjwICKJwJvAz1U1X8L0WcWq6gFGi0gS8BYwrLbFWjaqliciFwCZqrpKRM6onFzLou1+X7gmqupeEekBLBaRTf4UEuqavnXZcLQDItITwH3NDHE8LUJEonES/r9Vdb47OSz3RSVVzQOW4JznSBKRykpauPxOJgLTRGQnTtPvWTg1/3DcF6jqXvc1E6cycBJ+/EZCnfSty4ajvQtc7b6/GngnhLG0CLed9kVgo6r+yWdWOO6L7m4NHxGJA87BOcfxKXCZu1hY7AtVvVdVU1Q1FSc3fKKqMwjDfSEiCSLSsfI9cC5Oz8RN/o2E/I5cEZmK89+7ssuGR0IaUAsSkf8AZ+B0FXsAuB94G5gL9AV2A5eras2Tve2KiEwCvgDWUd12ex9Ou3647YsTcE7IReJUyuaq6kMiMgCnttsVWANcpaqldZfUvrjNO79U1QvCcV+4n/ktdzQKeF1VHxGRZJr4Gwl50jfGGNNyQt28Y4wxpgVZ0jfGmDBiSd8YY8KIJX1jjAkjlvSNMSaMWNI3xpgwYknfGGPCyP8HfLkInvE7Z7gAAAAASUVORK5CYII=\n",
"text/plain": [
"