PB_final.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Macroeconomics - Final Assignment\n",
    "\n",
    "### Paula Beck\n",
    "04.03.2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all necessary packages\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "\n",
    "# set seed for random numbers\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #1 Define matrices $A$ and $B$\n",
    "$ A z_t = B$ ; $z_t = [g_t, y_t, \\pi_t, i_t]'$\n",
    "\n",
    "$ A = \\begin{pmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    - \\bar{g} & 1 & 0 & (1-\\bar{g}) \\sigma^{-1} \\\\\n",
    "    \\kappa \\sigma \\frac{\\bar{g}}{1-\\bar{g}} & -\\kappa(\\frac{\\sigma}{1-\\bar{g}} + \\eta) & 1 & 0 \\\\\n",
    "    0 & -\\phi_y & - \\phi_{\\pi} & 1 \n",
    "  \\end{pmatrix} $\n",
    "\n",
    "$ B = \\begin{pmatrix}\n",
    "      \\rho_g g_{t-1} + \\epsilon_t^9\\\\\n",
    "      E_t y_{t+1} + (1-\\bar{g}) \\sigma^{-1} E_t \\pi_{t+1} - \\bar{g} E_t g_{t+1} \\\\\n",
    "      \\beta E_t \\pi_{t+1} \\\\\n",
    "      0\\\\\n",
    "      \\end{pmatrix}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup parameters in a dictionary\n",
    "params = {'sigma': 1, 'eta':2, 'kappa': 0.3, 'beta': 0.995, 'phi_pi':1.5, 'phi_y':0.1, 'rho_g':0.7, 'gbar':0.25}\n",
    "\n",
    "# define matrix A\n",
    "a31 = params['kappa']*params['sigma']*params['gbar']/(1-params['gbar'])\n",
    "a32 = -params['kappa']*(params['sigma']/(1-params['gbar']) + params['eta'])\n",
    "A = np.array([[1,0,0,0],\n",
    "                [-params['gbar'],1,0, (1-params['gbar'])*(params['sigma']**(-1))],\n",
    "                [a31,a32,1,0],\n",
    "                [0,-params['phi_y'],-params['phi_pi'],1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #2 Define msv function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define function\n",
    "def NKM(input, A = A, params=params):\n",
    "    '''\n",
    "    Returns the difference between input and the implied model which can be used to solve for the msv solution.\n",
    "\n",
    "    Arguments: \n",
    "        input: vector of length 8\n",
    "        A: matrix A\n",
    "        params (dict): parameters of the model\n",
    "\n",
    "    Returns:\n",
    "        dif: vector of differences between input and results of model\n",
    "    '''\n",
    "    # input: 8 elements\n",
    "    cg = input[0:4] # unpack first 4 elements und define them as cg\n",
    "    ce = input[4:8] # unpack next 4 elements und define them as ce\n",
    "\n",
    "    eg, ey, epi, ei = cg * params['rho_g'] # unpack expectations based on equation 6, shock is zero \n",
    "\n",
    "    # define B for gt_1 = 1 and calculate the implied values\n",
    "    B1 = np.array([[params['rho_g']],\n",
    "                    [ey + (1-params['gbar'])*(params['sigma']**(-1)) * epi - params['gbar']*eg],\n",
    "                    [params['beta']* epi],\n",
    "                    [0]])\n",
    "    # calculate the implied values by inv(A)*B\n",
    "    implied1 = np.linalg.inv(A) @ B1\n",
    "\n",
    "    # define B for epsilon t_1 = 1, gt_1=0 and calculate the implied values\n",
    "    B2 = B1.copy() # same as before\n",
    "    B2[0,0] =  1 # change upper left element\n",
    "    implied2 = np.linalg.inv(A) @ B2 # change implied values corresponding to the shock\n",
    "    \n",
    "    # put implied values for cg and ce in one vector\n",
    "    implied = np.append(implied1,implied2)\n",
    "\n",
    "    # calculate the difference between the input values and the implied values by the model\n",
    "    dif = input - implied\n",
    "    # function returns this difference\n",
    "    return dif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #3 MSV solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cg: [0.7        0.08116102 0.03677436 0.06327765] \n",
      " Ce: [1.         0.13059284 0.05620618 0.09736856]\n"
     ]
    }
   ],
   "source": [
    "# make a random guess about the 8 values of the model\n",
    "guess = np.random.randn(8)\n",
    "\n",
    "# use fsolve to calculate the minimum state variable solution \n",
    "msv = fsolve(NKM, guess) #  first input is the function, second is the guess where fsolve should start\n",
    "print('Cg:',msv[0:4], '\\n', 'Ce:', msv[4:8]) # first 4 elements are the values for Cg, the next 4 element correspond to Ce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #4 Calculate impulse responses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 30 # number of periods\n",
    "epsg = 0.01 # shock to g\n",
    "\n",
    "# unpack and reshape msv solution into vectors, Cg and Ce \n",
    "cg = msv[0:4].reshape(4,1)\n",
    "ce = msv[4:8].reshape(4,1)\n",
    " \n",
    "# start a list in which to save zt\n",
    "# first entry is Ce * error\n",
    "z = [ce * epsg]\n",
    "\n",
    "# loop through N periods\n",
    "for i in range(N):\n",
    "    zt = cg * z[-1][0] # new zt is Cg * gt-1\n",
    "    # z[-1] is last element from list (period t-1), [0] picks the first element from array -> gt-1\n",
    "    z.append(zt) # append new zt to list\n",
    "\n",
    "# unpack elements of zt using list comprehensions\n",
    "gov_spending = [x[0][0] for x in z]\n",
    "output = [x[1][0] for x in z]\n",
    "inflation = [x[2][0] for x in z]\n",
    "interest = [x[3][0] for x in z]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #5 Plot impulse responses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEdCAYAAAC/nDhCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABJt0lEQVR4nO3deXwU9fnA8c+TzQEJNwTkEDlVREAOAUUFFRHQivWo2FbU2nrVo61WsaL11h7219b7qEVt1SqKUrUKoniicojIqYgckXAfgQCBJM/vj+9sWDab7CbZZJjs83699pXdmfnOPDPZ7z4z3/nOjKgqxhhjTBCk+R2AMcYYkyhLWsYYYwLDkpYxxpjAsKRljDEmMCxpGWOMCQxLWsYYYwLDkpYxxpjAsKRljDEmMNL9DsAYY0xsItISeBA4EtgKfARkq+q1fsblJ0taxhhz4JoM/EdVzxeR9sA3wKU+x+SrhJoHReQiEVER6VbbAVWViNwmInV2L6qIbRF+7RGRb0XkHhFpUFdxpAIROVNEfuN3HFUR/j6KSJ3sENZkeXVdd2qLiPxMRL7x6uLWGszngNoeInIi0FhVHwJQ1e+BjcCsWl7uwSIySUS2iUiBiLwiIh2TVVZEOojIAyIyU0R2et/fTonGZ+e0qu9c4BjgNOBt4CbgT75GVP+cCQQqaQXMk7jvcGCJSDvgceAT4CRgeA1md6Btj37AnPAHEWkNNAK+rq0Fikg28C5wOHAhcAHQHXhPRHKSVLYb8CNgC/BhVWO05sHqm6eqy7z300SkO3CJiFyrqqV+BmZMIlQ1D8jzO44a6g6EgKdV9aOazOgA3B4bgDNFJIQ7wHgQ+EJr9y7nvwC6AIeFf99EZD6uWfIy4C9JKPuBqrbxxv8cGFGVAKt1pBXRJHG4iLwtIoUiskpELvbGXyAiS0Rkh4i8JyJdY5Tt5Y3bKSL5InKHiKRFTDdRRFbEWPYMEZkRJ75DRWSyiKwXkd1ebC9FNqGISB8RmSIiW0Rkl4h8LCLHV2d7eOYCDYFWUbHEXU68eBPdZhHzG+kdeu/yDtNfFZHDIsaH59ddRN7w/k8rReTW6Pkla1smMp+o6Sfi9tbay76m2BWJrmNFRKSbiDwrIt95ZZeLyCMi0jyBslVZh84JbNu46+Bt28kissmbbqmI3BQnzpHech+M9f2ImG6/5rCqfC8qmF+m951c6cX6joh0FZEiqUYzbwLf44nADO/jdC/2iZXML6F6FjH9MNn/VEDka2LEdNX+LRGRLuJ+674X17y5XlwdzwBeBLYBi4BXgWJgekIbr/rOAD6N2CFHVb8DPgbGJKNsTXfqa9o8+BLwBq4ZZw7wlIjcA1wBjAcuBg4DnotR9lXgHa/sc8AtwK01jCfsdaC9F8epXixFeOsrIv1wzQktcHsHZwObgHdEpH81l9kJ9wXbFB5QheVUGm+EV4mzzURkJO5/sgM4z5vnkcBH4k7kRpqMO5w/05v37bhEESlZ2zLRdQy7E3gTt7d5jPf6YTXWMVo73N70r7w47gBO9pYVT1XWodJtm8g6iMhAYCbQFfg1rin6L0CHigIUkXHAFOAPqnpVNX8gEvleRC9XcD+y1wIPAT8A5gPTgEzcTl3CEvwf3wlc473/Je47cmcls63qd3Au+7574dfvvHGLvTir/VsiIs1wP+jNvJhO9ubxkqruBTJU9XRVPQzXTN4beKaS+YmIpCfwClUSVk9gQYzhC4EjKlufGpZNnKrGfQEXAQp08z7f5n0eFzFNc9yewCagScTwa7xpD4kqOz5qGU8A24Fm3ueJwIoYscwAZkR8vs2tRtnnVt78z6hkfabjvnSZEcNC3rBXE9wWh+GaV5sDP/PW/aqqLifBeBPaZt6w2bjD8fSIYZ2BvcBfouZ3cdT8vgKmJntbJjKfCuY9EciLMTzuOlZhGenAcV58fSuZLqF1qMK2TeT/9AGwGtfFOd7y0oEbvPI/T3Ddo+tOQrFXMK9LvWUPiPG/KgWaVvH/ktD/GHcOS4FhceaXcD2rZHx3YDMwCZBEv/+VzG+0F9P5QHbkunrjf+PN5wvc794xceY3zJtfvNeMSuaxB7gvxvC7gOI4y69yWeDnXkydEv1u1PRI63/hN6q6BViPOzwsiJhmiff34KiyL0Z9fgF3kvHIGsa0CVgO3CcivxB3rqmMiDQEhuKOEkvDex+A4I5iTkhwOUtwFWgz8A/gMVV9sBrLqTTeKJVuM3EnO/vhusgWhyfSfYfoQ6PKvxH1eQEQ2dMnWduyKutYqWqsY3T5TBH5nbjm6124/2H4ZHBlzYtVXYcKt20i6yDupPYQ4N+qujPOsgD+D3dEdI6qPpnA9NWKvRI3AC+q6uyo4cuAb1V1G4CIpInIBO/ILKaa/o8rUKPvoLjm49dx63OBqmoSfktm4RLzc0AhsC5ypKr+RVV7qGpfVR2mqjPjzG8OcHQCr8vizCfWObMK/19JLJuQmiatLVGf91QwDCC6O/i6Cj7Ha96plLr0fQpuT+1e4Gtx5y2u8CZpgdsTugX3gxX5ugponkj7Pa6p6mjc3tI7wJVe00xYQstJIN5I8bZZc9wXJD9G2bVeTJE2R30uIuL/lKxtWcV1jKeq6xjtXtwe9b9wzW0DgbO8cRVeslCNdahs2yayDs1x9TPRjgHn45ph3klw+spU+r2IJiJdcE2Y0ckO3M7qnIjPvYDzvO1ZkZr+j8upyXfQS0STcNvgDFXd5Y2q6W9JDm7H/wrcDkpNzqmDa0qdl8BrGRXbQuzt25zyv+3JLJswP3sPtsHt+UR+Bvje+7sb1xYerSUR541iUdXlwDhvb64P7gv0sLgT+R/gmiseooL2YU3sPMAC3ddD5l1c+/2fRORlVS3EXb2e0HIqi1dV/xdRJN4224Lb0zkoxuIOIs52qyDGpGzLKqxjPDVdx7HAM6p6V3iAiDRKZMF1vA5bcNs20Z24k4GpwP9EZLSq7qhCPDXVzvu7NnKguLs59AFe8z4fgTt3mCEi84C3VHV8jPkl/XsMNfr/PYTbQR2iqpHruJVq/paISBNck9/VqvrfihYsrhPb00Br3NHYL2IczYYNBd6rZD3C3sc1JcayEHduKtoRuA4hlalJ2YT5eZ3Wj6I+j8XtKYRP5K0E2ohIWW887x8Yt4dYmDrz2Hetz5FeQvkQ96Wdq6qzo19VXRFVLQJ+i/tiXekNq/JyYsUbNUml28xb5hzg3MiTrSJyCHAs7staLcnalgmsY6QiXI/MyPI1Xcds3J5wpIvjlNlPFdchVvm46+A1CX4E/NRrhopnIe6HqDvwlog0rkpMNbTB+xtdN2/BHU3MAVDVRcArwD2qelQFCatWv8fe/BP+/4nIr4FLgLGq+lWMOKv7W3ICcAiwNE64jwITVfVQXBPsvytpWk1G8+AUYLB39AyAuAt/h3jjKlOTsgnz80jrF96h8yxcT56fA7ep6lZv/Eu4nkD/FpG/4E6k3oS7IrxCItIb+BvwH9xhcAjXeaIY1yMK3Jf1A+BtEfkHrhmiFa4dPVRRZaqMqk4RkVnA9SLyoNeEEHc5CcYbFm+bgfuheAN4XUQexp3zuh3Xs/H+qqxTsrZlFdcx0iKghdeEMxvY7f1w1GQd3wIuFJGvvFjOwv0QJmNbVEUi63A97gd6pojcj2sq7AIcpapXR89QVReLyDDc3vZbIjJSVbdXI7aq+gaXNH8vIoW4I//zcNtW2b/nYH/cNownad9jqN7/T0SOBf6MO4raLCKDI0ZvUNVvqf5vSTjRPyMif8IdPXYATlLVn3nLzwUG405BoKrTvHzVH1cf9uP9r6u80x3lCdwR6GsiMgH3/7sT1yHosfBEIjIU1wnlZ6r6TFXKeuXP8d6Ge1iOEpENuO1a+U6JJtaT5yJi9x6M7u2yAvhX1LBh3rTDo8oeiatcu3DNCncCaVFlz8QdRewCvsRdhDaDynsPtsYdTn8N7MS1z78PnBo17x64jgzrcXv0ebi9gdFV2RZR40Z4436d6HISibcq28ybfiSuq/QuXCV/DXfBX/T8ov9/E4nosZmsbZnofGKsRw7wPPuai1Ykuo6VzLOVF+sW7/Vv3N6nAhdVUi7RbZHQtk10HYC+wH9xTVG7cB2AbqxsebijrTxv3k0qWafbiN17MG7sMeZ1uLc9CnGtJHey72LY8DQhoADISfB3J5Htk2jvwYTrWYy6Hus1MQm/JT/FHR0VeDEtAiZEjO8PLI0qMxU4K5HtV90XrtPNy15c23GXPXSKmmYYMepMImW96SrarjPixRfutllnROQ24Pe4axCK40xusG1mgkdEOuP2+m9R1Ye9YR1wzWitfQ0uIMRd5/Wcuuu0wsOmAY+o6iv+ReYvu42TMaZGROQYYBTwKe5Iox+u2exL3H0Bw/KB2SKyFHhdVa+r61gDZhXQTkQy1F1sDO482CofY/KdJS1jTE1l427T8xsgA/gW16PuT7r/dVYleOdnTHyqukFEPsc1Uz4hIqfgLgWYU2nBeq7OmweNMcYkxrsI+mncudidwKWq+rm/UfnLkpYxxpjA8PM6LWOMMaZKLGkZY4wJDEtaxhhjAsOSljHGmMCwpGWMMSYwLGkZY4wJDEtaxhhjAsOSljHGmMCwpGWMMSYwLGkZY4wJDEtaxhhjAsOSljHGmMCwpGWMMSYwLGkZY4wJDEtaxhhjAsOSljHGmMCwpGWMMSYwLGkZY4wJDEtaJvBEZIWIDPfe/05EnvQ7JlM/iEgbEflARLaLyP0icpuI/KsG81soIsOSF2HqsaTlEZGxIvKZiBSKyHrv/ZUiIn7H5hcRmSEiP/c7jqpQ1XtUNVAxm7oXuaMTx6XARqCJql5XxWVMFJG7Ioepak9VnVGV+Zj9WdICROQ64G/An4CDgDbA5cAQILMO40ivq2UZYxJyCLBIVdXvQIxHVVP6BTQFCoGzE5juGWADsBKYgEv6WcBW4MiIaXOBXUBr7/PpwDxvuk+A3hHTrgBuBOYDRUA3QIELgVW4vbybI6a/DXgJ+BewHfgKOBS4CVgPrAZGRMX9DyAf+B64Cwh54y4CPgL+DGwBvgNGeePuBkqA3cAO4MEY26SBF8cmb91mAW28cTOAe4HPgW3Aa0CLiLKDvW2xFfgSGBYxbgZwJ/Cxt45TgVYR4y/w/gebgJu9bTg8Yvv8y3vfKc62bAg87a37YuAGIM/v76S9av8V/s7EqQMTgb3AHq8ODI/8fnnTvASs9b7jHwA9veGXRpX9b+RyvfdZwF+BNd7rr0CWN24YkAdch6vX+cDFfm+3A+FlR1pwDO7L81qc6R7AJYAuwFBgHO5LVAS8ApwfMe2PgPdVdb2I9AOeAi4DWgKPAVNEJCti+vOB04BmQLE37DjgMOBk4FYR6REx/Q+AZ4HmwBfA27gE2h64w1tG2NPePLsBfYERQGTz2SBgKdAK+CPwDxERVb0Z+BC4SlUbqepVMbbJhd42Odhbt8txyTpsHPAzoJ0Xw98BRKQ98AYugbYArgdeFpHciLI/Bi4GWuOOdq/3yh4BPIJLXO285XaIEVukirbl73GJrQtwCvDTOPMx9VNFdeAi4N/AH7068E6Msv8DuuO+p3O96VHVx6PK/iBG2ZtxO29HAX2Agbid4bCDcPWrPXAJ8JCINK/ZqgafJS33Rd2oquFkgYh8IiJbRWSXiJwgIiHgPOAmVd2uqiuA+3E/nADPsX/S+rE3DOAXwGOq+pmqlqjq07gjqsER0/9dVVerauQP/u2quktVv8QdifSJGPehqr7txfwS7sjuPlXdC7wAdBKRZiLSBhgF/EpVC1V1PfB/wNiIea1U1SdUtQSX4NrimkcTsReXNLp56zZHVQsixj+rqgtUtRC4BfiRty1/Crypqm+qaqmqTgNmA6Mjyv5TVb/2tsmLuIoNcA7wuqp+4O0w3AKUxomzom35I+AeVd2iqnl4SdWknGrXAVV9yvtNKMIdhfURkaYJLvcnwB2qul5VNwC3s+83BVz9ukNV96rqm7gjtsMSnHe9ZedQXBNTKxFJDycuVT0WQETycIm9FW5vf2VEuZW4PSCAd4GGIjII11RwFDDZG3cIcKGIXB1RNhN3lBC2OkZcayPe7wQaRXxeF/F+Fy7plkR8xpu+HZAB5Ef0J0mLWl7ZclR1pzdd5LIq8yzuKOsFEWmGayq82UueRC1npRdLK9w2OVdEIvc+M4D3YsXF/uvfLnK+qlooIpvixJnQvIj9fzD1X7XqgLcDdjdwLm7HMbzz1ArXXBhPO8r/pkT+LmyK3Jmm/O9ASrIjLZiJO/IZU8k0G3F7PYdEDOuIO0eEqpbijgbOxx1lva6q273pVgN3q2qziFe2qj4fMa/aOsm7GrdurSKW3URVeyZYvtK4vD3A21X1COBY3Lm7cRGTHBzxviNuG2704no2apvkqOp9CcSUHzlfEcnGHe1VRz77Ny0eXNGExsTwY9zvxnBcM14nb3h4DzFevV5D+d+UNUmMr15K+aSlqltxh+UPi8g5ItJIRNJE5Cggx5umBJeU7haRxiJyCPAb3JFF2HO4JsSfsK9pEOAJ4HIRGSROjoicJiKN62Dd8nGdGO4XkSbeenUVkaEJzmId7nxPTCJyooj08vY4C3BJqSRikp+KyBFeYrkDmORty38BPxCRU0UkJCINRGSYiMQ7NwUwCThdRI4TkUxvvtX9Hr8I3CQizb3zbLHO2xlTkca4ncJNQDZwT9T4SusP8DwwQURyRaQVcCv7/6aYGFI+aQGo6h9xSegGXE+ddbjODDfiergBXI3rZbgc19voOVwHi/A8PvPGt8OdnA0Pn407r/UgrnfSMlyPpboyDtccuchb/iRcm30i/gacIyJbRCTW+Z6DvPkV4Hrfvc/+le5ZXA+stbiehtcAqOpq3B7q73C9MVcDvyWB76OqLgR+idv++d465SW4PtHu8Mp+B7zjrUtRNedlUs8zuCa973H169Oo8f8AjvDOj78ao/xduHO583G9gOd6w0wlRNUuPzDJJyIzcF2DA3N3ChG5AhirqokeiRpj6pgdaZmUJSJtRWSI12x6GO6amMnxyhlj/GO9B00qy8Q1A3fGXeT8AvCwnwEZYypnzYPGGGMCw5oHjTHGBEbSmgdFZCSut1kIeDL6mhvvbul/w931YCdwkarO9cY9hbvGZ72qHhlRpgXwH9z1DyuAH6nqlsriaNWqlXbq1Ck5K2VMAubMmbNRVXPjTxkcVo+MHxKpS0lJWt51Og/h7t+WB8wSkSmquihislG4e3R1x93r6xHvL7hu0Q/iupBGGg9MV9X7RGS89/nGymLp1KkTs2fPrtkKGVMFIrIy/lTBYvXI+CGRupSs5sGBwDJVXa6qe3AntKPvMDEGeEadT4FmItIWQFU/ADbHmO8Y3L3A8P6emaR4jTHGBFCyklZ79r9vWx777stXlWmitfHu6hC+u0Pr6gb4+uuv89hjj8Wf0BhTob0FW1n++J8pWDTP71BMikpW0or1dN/obomJTFO9hYtcKiKzRWT2hg0bYk4zadIkbr31Vqy3pDE1ISy597esf+9NvwMxKSpZSSuP/W822oHyN35MZJpo68JNiN7f9bEmUtXHVXWAqg7IzY19Du+YY45h/fr1LF++PM4ijTEVyWjSlOxO3ShYMNfvUEyKSlbSmgV0F5HO3k1MxwJToqaZAozzbho7GNgWbvqrxBTcgwbx/sZ7UGOFjjnmGABmzpxZ3VkYY4CmPfuxbaElLeOPpCQt75kvV+GeoLsYeFFVF4rI5SJyuTfZm7ibzS7D3fn8ynB5EXke94iQw0QkT0Qu8UbdB5wiIt/geiYm8uiKmHr27Enjxo355JNP4k9sjKlQkyP7sWv1d+zdVunVJ8bUiqRdp+U9WfPNqGGPRrxX3N25Y5U9v4Lhm3CPSK+xUCjEoEGD+O6775IxO2PiquG1izHLVnTtooh0wu0wLvVm/6mqXu6V6Y+7rKQhro5eqzU4udukZz8kPZ0dy5fSvO/g+AWMSaKUuvfga6+9RnZ2tt9hmBRQk2sX45St7NrFb1X1qBjhPAJcint0xpvASCIen1NVLQcPY8SC7YSyGlR3FsZUW0rdxskSlqlDNbl2sbKyVbp20ZtfE1Wd6R1dPROvTDxpGRmWsIxvUipplZSUMHbsWB5+2G7kbWpdTa5drKxsZdcudhaRL0TkfRE5PmIZkQ/JrPD6yEQuHSmbyctP8+V1F1Y6jTG1IaWSVigUYsGCBbzxxht+h2Lqv5pcu1idaxrzgY6q2hf3FO7nRKRJVeaVyKUjYbu+X8n3k5+leMf2OGEZk1wplbTAdX2fOXMmpaWlfodi6reaXLtYWdmY1y6qapHXcQlVnQN8CxzqzatDnDiqrEnPfqBKweIvazorY6okJZPWli1b+Prrr/0OxdRvNbl2sbKyMa9dFJFcrwMHItIF17ljuTe/7SIy2OutOI4aXO8Y1vTIfgAULPyiprMypkpSqvcgwLHHHgvAJ598wuGHH+5zNKa+UtViEQlfuxgCngpfu+iNfxTXk2807trFncDFlZX1Zn0f8KJ3LeMq4Fxv+AnAHSJSDJQAl6tq+CbUV7Cvy/v/qEHPwbCs1m3JbNXGLjI2dS7lktahhx7K0KFDycrK8jsUU8/V8NrFcmW94TGvXVTVl4GXK5jXbODIWOOqS0RoddwphBpYj1xTt1IuaaWlpTFjxgy/wzAm8I76v2f9DsGkoJQ7pxVWXFxMcXGx32EYY4ypgpRMWl988QXNmjXjnXfe8TsUYwJrz9bNfHDKEax6/gm/QzEpJCWTVvfu3dm1a5fd8d2YGsho2pyiDWvZ9tVsv0MxKSQlk1ajRo3o3bu33fHdmBoQEZoc2Y8C60Fo6lBKJi1w12t99tlnlJSU+B2KMYHV5Ii+bF8yn9K9e/0OxaSIlE1axx57LNu3b2fhwoXxJzbGxNT0yH6U7tnDjm8WxZ/YmCRI2aQ1dOhQfv/739OsWTO/QzEmsJodNYj2Z1+IZGT4HYpJESl3nVbYwQcfzG233eZ3GMYEWnbHLvT580S/wzApJGWPtACKiop48803KSgo8DsUYwJLVSnasNbvMEyKSOmk9cUXX3Daaafx+uuv+x2KMYG15N4bmDGsG2qdmkwdSOmkNXDgQNq1a8fLL8e8ZZsxJgHN+gykZGchm2d96HcoJgWkdNJKS0vjhz/8If/73//YuXOn3+EYE0i5J44mlJ3Dmv++4HcoJgWkdNICOOuss9i1axdvv/2236EYE0jp2Tm0Hn4Ga/83ya7XMrUu5ZPWCSecQIsWLZg2bZrfoRgTWO1OH8veLZvY9Ml0v0Mx9VzKdnkPS09PZ9asWXTq1MnvUIwJrFYnnErfh16ixaChfodi6rmUT1oAXbp08TsEYwItlJVF29Hn+B2GSQEp3zwYNn78eCZMmOB3GMYEVvHOQpY9eDebZr7ndyimHrOk5VmxYgWPP/643UDXmGpKy8hkxT//yqrnHvM7FFOPJS1pichIEVkqIstEZHyM8SIif/fGzxeRfvHKishtIvK9iMzzXqOTFW+0s846iw0bNvDxxx/X1iKMqdfSMjI4aNQ5rJ/+X4p3FvodjqmnkpK0RCQEPASMAo4AzheRI6ImGwV0916XAo8kWPb/VPUo7/VmMuKNZdSoUWRlZfHKK6/U1iKMqffanj6Wkl07Wf+u3WXG1I5kHWkNBJap6nJV3QO8AIyJmmYM8Iw6nwLNRKRtgmVrXePGjRkxYgSvvPIKqlrXizemXmhx9HFktWlHvl1obGpJspJWe2B1xOc8b1gi08Qre5XXnPiUiDSPtXARuVREZovI7A0bNlR3HRg3bhwjRoxgx44d1Z6HMalMQiHaj/kJpKXZzp+pFclKWhJjWPQ3tqJpKiv7CNAVOArIB+6PtXBVfVxVB6jqgNzc3IQCjuWcc87hySefpHHjxtWehzFhtXSet4WITBORb7y/zb3hp4jIHBH5yvt7UkSZGd68wueGW9fmeh82/g/0f+RlRGJVbWNqJllJKw84OOJzB2BNgtNUWFZV16lqiaqWAk/gmhJr3aeffsrKlSvrYlGmnqrF87zjgemq2h2Y7n0G2Aj8QFV7ARcCz0Yt6ycR54bXJ29Nywsnq13fr6rNxZgUlaykNQvoLiKdRSQTGAtMiZpmCjDO27scDGxT1fzKynrnvMJ+CCxIUrwV2rRpE0OHDuW+++6r7UWZ+q22zvOOAZ723j8NnAmgql+oanhHcSHQQESyamnd4lrz+n9474TOFCz60q8QTD2VlKSlqsXAVcDbwGLgRVVdKCKXi8jl3mRvAsuBZbijpisrK+uV+aPX3DEfOBH4dTLirUzLli0ZN24cEydOZP36Wt0hNfVbbZ3nbePt7OH9jdXUdzbwhaoWRQz7p9c0eItU0G6XrHPDALnHjyA9pxHLHryrRvMxJlrSrtNS1TdV9VBV7aqqd3vDHlXVR733qqq/9Mb3UtXZlZX1hl/gTdtbVc8IV9badt1117F7924eeuihulicqZ9q6zxv5QsV6Qn8AbgsYvBPvGbD473XBbHKJuvcMEBG0+YcMu5q1r71Mtu/Xhi/gDEJsjtixHD44YczZswYHnzwQQoL7SJJUy21cp4XWBduNvf+ljUHiEgHYDIwTlW/DQ9X1e+9v9uB56ijc8OdL/k1oYbZLHvo7vgTG5MgS1oV+O1vf4uqsmBBrZ9GM/VTrZzn9f5e6L2/EHgNQESaAW8AN6lq2W1dRCRdRFp57zOA06mDc8MAmc1bcsgFv2T99P+yZ+vmulikSQF2l/cKDBkyhLy8PLKzs/0OxQSQqhaLSPhcbQh4Knye1xv/KO4872jced6dwMWVlfVmfR/woohcAqwCzvWGXwV0A24RkVu8YSOAQuBtL2GFgHdw55TrRNcrxtP559eR2axFXS3S1HNS3y4AHDBggM6ePTv+hAkqLS1l3bp1tG3bNv7EJiWJyBxVHeB3HMmU7HqkqpQW7SbUoGHS5mnqn0TqkjUPxnHmmWdy+umn29X9xlSTlpby2Y9PYuHvr/Y7FFMPWNKK4+yzz2bu3Lk88USdtagYU69IWhqNDz2S7195ml15dtG+qRlLWnGMGzeOk08+meuvv57Vq1fHL2CMKafL5TciaWksvd8etGpqxpJWHCLCE088QUlJCZdffrk1ExpTDQ3bdqDLZTey5tV/kf+/l/0OxwSYJa0EdO7cmXvvvZcVK1awadMmv8MxJpC6XX0LTXsPYMVT/2c7f6barMt7gq666iouu+wysrJ8u52bMYGWlpFBv0cnk9mshd0B3lSbHWklKC0tjaysLLZv326dMoyppoZtOxBqmE3xzkI2fvKu3+GYALKkVUX//Oc/ufTSS5k0aZLfoRgTWEvu/S2zf3Ya279Z5HcoJmAsaVXRlVdeyYABA7jkkkuYO3eu3+EYE0jdr7mV9JxGfPmrn1C6Z4/f4ZgAsaRVRenp6UyePJnmzZtz6qmnsmTJEr9DMiZwsnIPote9T1KwaB5f/vYiSouL/Q7JBIQlrWro0KED06ZNIxQKcfHFF1tPKGOqoc2IMRx2433kT3mepX8YH7+AMVjvwWrr3r07U6dOpXHjxtYTyphq6nr5jaQ3bkrrYaP9DsUEhB1p1UDv3r3p3LkzpaWl/OlPf2LLli1+h2RM4Bzyk8tp2L4jWlLCquefsKZCUylLWkmwYMECbr75ZgYPHsz8+fP9DseYQNrwwdss+N2lfPmbCygpKvI7HHOAsqSVBL1792bq1KkUFBQwaNAgHn/8cTvPZUwVtT5xNIeP/yP5/32BT84cSMGSr/wOyRyALGklybBhw5g3bx7HHXccl112Gddee63fIRkTOF0u+y0DnnqDoo3r+GTMAPJeftrvkMwBxpJWErVp04a33nqLO++8kzPOOAOAYmufN6ZKWp84muPf+orWJ59Bo649/A7HHGCs92CShUIhJkzY9/iFW265hVmzZnH33XczaNAgHyMzJjiyWubS7+GXyj4v/P3VhBo0pPPPryMrt42PkRm/2ZFWLevUqRPz589n8ODBjBkzxjpqGFNFqkpx4XaWP3k/753QmUV3/Ird69b4HZbxiSWtWnbZZZexfPly7rrrLt5//3369OnDvffe63dYxgSGiNDnzxMZ+s4S2p5+HiufeZAZJ3Rm7dRX/Q7N+MCSVh1o1KgRN998M8uXL+cPf/gDp5xyCgCzZ8/mnHPOYdKkSWzevNnnKI05sOV07k6fP/2Toe99Q/tzLqZx954AbHj/LRZMuIJNn75v13ilADunVYdatGjBDTfcUPZ51apVfPDBB7z88stub7JPH0466SRuvfVWmjZt6mOkxhy4sg/uTK+7Hy37vGP5Ur5/5RlW/ftRQjmNaN7vWFoMPIGuV4xHQiEfIzW1IWlHWiIyUkSWisgyESl3IzFx/u6Nny8i/eKVFZEWIjJNRL7x/jZPVrwHgrPOOovvv/+ejz76iDvuuIPmzZvzzDPPkJ2dDcCtt97Kaaedxk033cRzzz3HrFmz2Lhxo89Rm0TVdZ0QkZu86ZeKyKkRw/uLyFfeuL9LPbvvWOeLr+Xk2evp++CLtP/hOHavX0PeyxPLEtai26/li2vOZ9kDd7H27cns+HYJJbt2+hy1qa6kHGmJSAh4CDgFyANmicgUVY18WM4ooLv3GgQ8AgyKU3Y8MF1V7/Mq7njgxmTEfKDIyMhgyJAhDBkyhAkTJlBcXEx6uvu3ZGZmkpeXx7Rp09i7dy8AXbt2ZdmyZQBMmDCB/Px82rRpQ+vWrWnTpg1dunQp66W4detWcnJyyMjI8GflUlhd1wkROQIYC/QE2gHviMihqlrizfdS4FPgTWAk8L/a3QJ1Kz07h7annUvb084FoKRod9m4kt272PrFp+T/94WyYS0GD2Pw8+8BrmciqmS2zCWzhXvldDmMJj16A7C3YBvpOY3sqO0AkazmwYHAMlVdDiAiLwBjgMgKOgZ4Rt2tIj4VkWYi0hboVEnZMcAwr/zTwAzqWdKKFk5Y4JLShAkT2LNnD0uXLmX58uX7Xfe1cOFCPv/8c9avX182/KSTTmL69OkAHH300SxbtowGDRqQk5NDTk4Op512Gg8//DAA5513Hrt27aJBgwZkZWWRmZnJcccdx8UXX1y2/HBM4deAAQMYPnw4JSUlPPzww6SlpREKhUhLS0NEOOqoozj66KPZvXs3L774IiKy36tPnz707NmTwsJC3njjjbLhQFn5rl27sm3bNt59992y4WH9+vWjY8eObN68mY8++qjc9hswYADt2rVj/fr1fPbZZ+XGDx48mNzcXPLz85kzZ0658UOGDKF58+asXr2aL7/8smz4oEGDyM3NTeRfGFbXdWIM8IKqFgHficgyYKCIrACaqOpMb17PAGdSz5JWtFBWg7L3ve59HIDiwh3s+GYRhcuXkt6kWdn4rfM+Y+fKZezdtu/eoe3PGkef+92FzdOPbkPpniJCOY1Iz25EqGE2Hc65mG5XT6B0zx7mXDqGtMys/V6tTzqdNiPGULyzkGUP3kVaKB1JT0fSQkh6Oi0Hn0izvoPYW7CV7195FtLSkFAIkTQkLY3mA4bQqFsP9mzdzPrpr7s6EPFq3v9Ysjt0omjTBjbPfM8Nh33j+x1DgzbtKNqwli1zZ+7bMN50zfsPIatlLrvXfs+2r8rXgxYDjyejaXN25a2kYEn5Hs8tBw8jvVFjdq5aHvNBnq2GnEyoQUN2LP+awu++3vd/adCQVkNOTuA/WLFkJa32wOqIz3m4Pcd407SPU7aNquYDqGq+iLSOtXARuRS3J0nHjh2ruQoHrszMTHr16kWvXr32Gz558mQASktL2bp1K+vWrdvv9lHjx48nPz+fgoICCgsLKSws5NBDDy0bv2XLFjZu3Mju3bvZvXs3e/bsoXHjxmXj//rXv7J7925KSkrKhl199dUMHz6cvXv3cs0115SL9aabbuLoo49m+/btXHjhheXG33PPPfTs2ZP169dz3nnnlRv/wAMPcNVVV7Fq1SrOOuuscuMnTpzIhRdeyJIlSxgzZky58ZMmTeLss8/miy++KLvAO9Lbb7/NiBEj+Pjjjzn33HPLjf/kk0845phjmD59elnyBpg6dWpZB5oE1XWdaI87koqe117vffTwcup7PUrPaUSzowbS7KiB+w0f8trnAJTu3cverZso2rSBUIOGgOtuf9gN91K8o4C927dRsrOQkl07adDObR8tKWbPlk2U7ina75XdqTttgJLC7Xz35P1ocTFE1M3Dx/+RZn0HsWfTBhbdXr4e9bzzYRp168Gu71cy//ry9ajPX54lu0MnCpct5oury9ej/o+9SoMRY9i2YC5zLy9fjwb+ezpZx57E5tkfMe/qseXGHzv5M5odNZCNH03jq5t+UW78CdMW0ahbD9ZNe43Fd/2m3PgTP1lNw7YdyH/9P3zzf7eWDW/YoRMnfvhduemrIllJK1YbefTN9yqaJpGylVLVx4HHAQYMGJByN/1LS0ujRYsWtGjRYr/hl1xySaXlpk6dWun4HTt2AK7ilpSUUFxcXHbUk5WVxcaNGykpKaG0tLTs1ahRI8B1Olm2bBmqut+rVatWALRv354FCxaUJdnw37Zt2wLQrVs35s2bV+4ejuEf0969e8c8UurcuTPgjqhmz5693zhVpXv37oA7Ip01a1a58ocffjgAp59++n7jw+WqoK7rRI3nlfL1KCODrNyDyMo9qGyYiND5kl9XWCbUMLss6cWSlXsQo752T2bW0lK0pAQtLUHSXFNjdscuDJ+zwY0rLQFVtLSUDO9IsFG3Ixg6Y5lLeF4dQrUsxia9+nP82wvcwsLjgYZeUm0+4DiOe2OeF82+f2l2x64A5B4/giH/LV+PcrocBkCbEWfSpGffcuMbdugEQLszfkzzo48vv94t3b7UweddQu6wUWXD0zIyK9xWiUpW0soDDo743AGIvvqvomkyKym7TkTaenuUbYH1SYrXVIGIlDUPRg5r2bJlhWVCoRBdu3atcHxmZiY9e/ascHzDhg3p06dPheMbNWpEv379KhzftGlT+vfvX+H4WEk+UqtWrcoSbDXVdZ2oaF553vvK4jB1QNJc0x/sO8csoRCZLSr+noWyssg5pOJ6lJ6dQ+NDK65HGY2bkHFExfUoo2lzmjatuH9bZotWlcaXldum0juUNGjTjgZt2lU4vjqS1XtwFtBdRDqLSCbuhPCUqGmmAOO8HlODgW1eM0dlZacA4WPjC4HXkhSvMbWtruvEFGCsiGSJSGdc547PvfltF5HBXq/BcVg9MgEmyXqEhoiMBv4KhICnVPVuEbkcQFUf9SrMg7ieSzuBi1V1dkVlveEtgReBjsAq4FxVrfQqXBHZAKysYHQrIMh9xoMcf5Bjh8rjP0RVy/XSqOs6ISI3Az8DioFfqer/vOEDgIlAQ1wHjKs1TsWPU4/ibY8DXZBjh/odf8y6FClpSSsIRGS2qg7wO47qCnL8QY4dgh9/sgV5ewQ5drD47TZOxhhjAsOSljHGmMBItaT1uN8B1FCQ4w9y7BD8+JMtyNsjyLFDisefUue0jDHGBFuqHWkZY4wJMEtaxhhjAiMlkla8R0QcaETkKRFZLyILIoYF5jEtInKwiLwnIotFZKGIXOsND8Q6iEgDEflcRL704r/dGx6I+GuT1aW6Y/UotnqftGTfYx5GAUcA54t7jMOBbCLugtNI4UdSdAeme58PVMXAdaraAxgM/NLb5kFZhyLgJFXtAxwFjPTuWBGU+GuF1aU6Z/UohnqftIh4RISq7gHCj3k4YKnqB0D0nT/G4B5Fgff3zLqMqSpUNV9V53rvtwOLcXcWD8Q6qLPD+5jhvZSAxF+LrC7VIatHsaVC0qro8Q9Bs98jKYCYj2k50IhIJ6Av8BkBWgcRCYnIPNwNaaepaqDiryVWl3xi9WifVEhaNX70iakeEWkEvIy7D16B3/FUhaqWqOpRuLuiDxSRI30O6UBgdckHVo/2lwpJK5FHRATBOnGPokAC8JgWEcnAVbR/q+or3uBArQOAqm7FPR14JAGMP8msLtUxq0flpULSSuQREUEQmMe0eHcv/wewWFX/EjEqEOsgIrki0sx73xAYDiwhIPHXIqtLdcjqUQWinyxbH1/AaOBr4FvgZr/jSSDe54F89j0q/RKgJa6nzTfe3xZ+x1lJ/Mfhmo3mA/O81+igrAPQG/jCi38BcKs3PBDx1/K2sbpUd7FbPYrxsts4GWOMCYxUaB40xhhTT1jSMsYYExiWtIwxxgSGJS1jjDGBYUnLGGNMYFjSMsYYExiWtIwxxgSGJS1jjDGBYUnLGGNMYFjSMsYYExiWtIwxxgSGJS1jjDGBYUnLGGNMYFjSCjgRuUhEvhKRnSKyVkQeCT/DJoGyK0RkeBJjSer8jElEVb53IjJDRH5e2zHFWG4nEVERSU/yfFVEuiVzngc6S1oBJiLXAX8Afgs0BQYDhwDTvIf0GWOSSERCdbispCa4+sKSVkCJSBPgduBqVX1LVfeq6grgR7jE9VMRmSgid0WUGSYied77Z4GOwH9FZIeI3BCxN3ipiKwRkXwvMYbLV2l+tb8VjNmf1/LwkYj8WUS2iMh3IjLKG3c3cDzwoPcdfdAbfriITBORzSKyVER+FDG/iV7rxZsiUgicKCLtRORlEdngzf+aiOkHishsESkQkXUiEn7i8Afe363eso+JEfttIjJJRP4lIgXARd78ZorIVq8+PhjeIRWR8Dy/9OZ5njf8dBGZ55X5RER6J3Uj+83vp1vaq9pPBR0JFAPpMcY9jXti60Tgrojhw4C8iM8rgOERnzvhnpT6PJAD9AI2hKep6vzsZa+6eEV+74CLcE8p/gUQAq4A1kDZA29nAD+PKJsDrAYuBtKBfsBGoKc3fiKwDRiC28nPBuYAtwKZQBdgOXCqN/1M4ALvfSNgsPc+XLfK1deIWG7zYj/TW1ZDoD+uBSXdm8di4FcRZRToFvG5H7AeGOSt/4Xe9sny+/+UrJcdaQVXK2CjqhbHGJfvja+u21W1UFW/Av4JnF+DeRlT11aq6hOqWoLbgWsLtKlg2tOBFar6T1UtVtW5wMvAORHTvKaqH6tqKW5HLldV71DVPaq6HHgCGOtNuxfoJiKtVHWHqn5axdhnquqrqlqqqrtUdY6qfurFtgJ4DBhaSflfAI+p6meqWqKqTwNFuMRXL1jSCq6NQKsK2r3beuOra3XE+5VAuxrMy5i6tjb8RlV3em8bVTDtIcAgryltq4hsBX4CHBQxzeqo6dtFTf879iXFS4BDgSUiMktETq9i7JHLQkQOFZHXvU5WBcA9VL5DeghwXVR8B1OP6rAlreCaiduDOityoIjkAKOA6UAhrjkjLLIigmtaiOXgiPcdcc0r1GB+xhwoor+jq4H3VbVZxKuRql5RQZnVwHdR0zdW1dEAqvqNqp4PtMZ1kprk1clE60b0dI8AS4DuqtoElyClkvKrgbuj4stW1ecTXP4Bz5JWQKnqNlxHjAdEZKSIZIhIJ+AlIA94FpgHjBaRFiJyEPCrqNmsw7XJR7tFRLJFpCeurf8/3vDqzs+YA0X0d/R14FARucCrQxkicrSI9Kig/OdAgYjcKCINRSQkIkeKyNEAIvJTEcn1mhK3emVKcOeGS6l6/WgMFAA7RORw3Dm6ytbnCeByERkkTo6InCYijau43AOWJa0AU9U/4va8/oz7Yn+G29M6WVWLcInrS9yJ2KnsSz5h9wITvGaE6yOGvw8swx2t/VlVp3rDqzs/Yw4UfwPO8XoW/l1VtwMjcOek1uCaFv8AZMUq7J0n+wFwFPAdrhn+SdwlJ+A6SC0UkR3essaq6m6vmfJu4GOvfiR6jul64MfAdlxCiq5ztwFPe/P8karOxp3XehDYgqvHFyW4rEAI96gxBu9I7Tsgo4IOHsYY4ys70jLGGBMYlrSMMcYEhiUtU0ZVV6iqWNNg7fM6zywVkWUiMj7GeBGRv3vj54tIv3hlReRcEVkoIqUiMiBqfjd50y8VkVNrd+2MqT2WtIypY+LuX/cQ7tKEI4DzReSIqMlGAd2916W4rs/xyi7AXQLxQeSMvPFjgZ64jgIPSx3eQ8+YZKp3N2Rs1aqVdurUye8wTAqZM2fORlXNrUKRgcAy724KiMgLwBhgUcQ0Y4Bn1PWU+lREmolIW9ytfGKWVdXF3rDo5Y0BXvB6lH4nIsu8GGZWFKDVI+OHROpSvUtanTp1Yvbs2X6HYVKIiKysYpH27H/ngzzcveLiTdM+wbKxlhd5O6HwvPYjIpfijuro2LGj1SNT5xKpS9Y8aEzdi3VHg+hrTyqaJpGy1Vkeqvq4qg5Q1QG5uVU5cDSm7qRM0vr4zzCt3OluY3yRx/63yurAvltlxZsmkbLVWV5CSravZ93E89m59J3qFDemxlImaW1YCF/Vm7tvmYCbBXQXkc7es5HGAlOippkCjPN6EQ4GtqlqfoJlo00BxopIloh0xnXu+Lw6gUuDJhTOe4nd38yoTnFjaixlklbrXrB9Dezc5HckJtV5lxRcBbyNez7Si6q6UEQuF5HLvcnexD2naRnu9j1XVlYWQER+KO6hnMcAb4jI216ZhcCLuI4ebwG/9G5HVGVpGQ3IaNODorwvqlPcmBqrdx0xKtKml/u7fgF0quxpNMbUAVV9E5eYIoc9GvFegV8mWtYbPhmYXEGZu3H3vquxrIP7sWvJ1PgTGlMLUuZIK5y01s33Nw5jgi6rfV9KCtZSXLA2/sTGJFnKJK2cNnDQUWD3BzamZrIOOZrM9kdRumOD36GYFJQyzYMicNksv6MwJvgadBlChxvtnJbxR8ocaRljjAm+lEpay6fD3w+Dzcv8jsSYYNv85m3k/bG/32GYFJRSSatBU9iyHNZ95XckxgSbpGeyJ28upbu2+R2KSTEplbRyjwDEkpYxNZXVwT0ppShvnr+BmJSTUkkrIxtadnfXahljqi+zQ18A9thFxqaOpVTSAmh9pB1pGVNT6U3aEGrSlqK8uX6HYlJMynR5D+s+CrJbueu1yj92yBiTqMbHXkqokd0N3tStKh1p1YdHhPe9CE5/yBKWMTXVYvRtND0h5p2mjKk1CSet+vSI8NIS2LMjGXMyJrWVFu2gtMgqk6k7VTnSKntEuKruAcKP+Y5U9ohwVf0UCD8ivMKyqrpYVZfGWF7ZI8JV9Tvc3a4HVmntYlCFv3SE6RNqOidjUlvxltWsuKEJO2b9y+9QTAqpStKq6PHfiUyTSNnqLA9wjwkXkdkiMnvDhsrvhyYCzTpZD0JjairUrANpDZraY0pMnapK0jogHxEOVX9MeLgHod0815jqExEyO/S1bu+mTlUlaQX2EeHR2vSCXZthR34y5mZM6srq0Jc9+V+hJcV+h2JSRFWSVmAfER6t7Nladr2WMTWS2aEvunc3e9ct8TsUkyISvk5LVYtFJPyY7xDwVPgR4d74R3FPUx2N6zSxE7i4srLgHhEOPADk4h4RPk9VT/XmHX5EeDE1eER4tIP6wNBboXnnZMzNmNTVsNtQWp7zAGl2vZapI6L17MTOgAEDdPbs2X6HYVKIiMxR1QHxpwwOq0fGD4nUpZS7jVPYri2wxuqkMTVWvGU1u76Z4XcYJkWkbNL66I/w1FAo2et3JMYE25ap97D2iTFoaanfoZgUkLJJq00vKNkDm772OxJjgq1B52PQ3QUUrZrldygmBaRu0jrS/bUehMbUTPaRZ0Aokx1zX/A7FJMCUjZptTocMhvD8nf8jsSYYAtlNyO752gK5/4HLU1KB19jKpSySSuUCUecBYtegeLdfkdjTLA16jeWkoJ8uzuGqXUp9zytSMfdAMf8CtIb+B2JMcGWfeQZHPz778ho2cnvUEw9l9JJq+WhfkdgTP2QltmQNEtYpg6kbPNg2MYlMPlCKFzvdyTGBFvxtjWsffwMdi6e6ncoph5L+aRVWgLzn4MFL/odiUkltfQU8BYiMk1EvvH+NveGdxKRXSIyz3s9WhvrFMpuye5vP2THbHu+lqk9KZ+0Wvd09yKc/2+/IzGpohafAj4emK6q3YHp3uewb1X1KO91ea2sV0YWOX3OpnD+q5Tu2VUbizDGkhZAr5+4WzptjPX8ZGOSr1aeAu79fdp7/zRwZi2vRzk5/ceiRdvZuejNul60SRGWtIBe54GkuWZCY+pAbT0FvI33KCC8v60jpussIl+IyPsicnysoKryBPCKNOx+IqHGbSi0C41NLUnp3oNhjdtB7x9Dg6Z+R2JSRF0/BTwf6Kiqm0SkP/CqiPRU1YL9ZqL6OPA4uLu8x5lnTJIWoumJv6lOUWMSYknLc+Y//Y7ApJCaPAU8s5Ky60Skrarme02J6wFUtQgo8t7PEZFvgUOBWnnOQbPhN9TGbI0BrHlwP6Ulrgu8MbWstp4CPgW40Ht/IfAagIjkeh04EJEuuM4dy2tv9aB07252ffthbS7CpChLWhHevAb+cQIUF/kdianPVLUYCD/JezHwYvgp4OEngeOeAr4c9xTwJ4ArKyvrlbkPOEVEvgFO8T4DnADMF5EvgUnA5aq6uTbXcdu7fyb/bydQ9P382lyMSUH25OII306Df42GMx6HvhcnOTBTb9mTi8srKdzMqts6kd1jJG1+ZhdBmsQk/cnF9fGCyEhdhkPbfvD+3e5ZW8aY6gnltKDp0Gso/HISe/IXxi9gTIISTlr19YLISCJw0h2wbSXMfaq2l2ZM/db0xF8jGdlsmXq336GYeqQqR1r19oLISF1HwMHHwpJX/YzCmOAL5bSkyQlXsXfdEnSvnSg2yVGVpHVAXhAJybkoct+84Ef/gZ+8UaPZGGOA5qNuo/31s5GMLL9DMfVEVZKWXxdE9gV+AzwnIk1iTaiqj6vqAFUdkJubG2e28TU6CNJCUFQAe3fWeHbGpKy0jAZIWholO7dSXLDW73BMPVCVpFWTCyIrK7vOa0Ik+oJIVd3kvZ8DhC+IrBOFG+Bvh8KnD9TVEo2pn3RvEXn3HMHmKeX6bhlTZVVJWvX+gshIOblw8GD45M+we2tdLdWY+kcysmjUbyw7Zv+LvRu+9TscE3AJJ61UuCAy2rDfu4Q18291uVRj6p+mJ/8WCWWwecqN1LdrQ03dsouL43hpLCx7Gy6fA827JG22ph6xi4sTs2XqPWx5/WZyf/oMjQdekNR5m/oh6RcXp6JT/gBp6fDVf/yOxJhgazb8Rhp0PZ6iFTP9DsUEmN3lPY5mh8Dlc6HpwfGnNcZUTNJCHHTFW6RlZvsdigkwO9JKQDhhrV8I373nbyzGBFk4Ye1Zs4Dtn9rzgEzVWdJKkCq8fgW8dD4U5PkdjTHBtvXdP7HhhUvZvXKW36GYgLGklSAROOMJKN4Nr4yD0mK/IzImuFqe9VdCTduy4ZmfUFpU6Hc4JkAsaVVBq8PgtAdg5Yfw3m1+R2NMcIWym9P6p8+wd+My1j35Q0r37PI7JBMQlrSqqM8F0Pdn8NEfYNHLfkdjTHA17D6M3B8/xa6v32Hbu3/2OxwTENZ7sBpOfxja9YPDzvA7EmOCrfGgiwg160DDrif4HYoJCDvSqoa0EAy4DEIZsGOtayosLfE7KmOCKfuw4Uh6JiU7NrLxxSvtHJeplCWtGlo8GT64GyZfCCV7/Y7GmODa/d0nFHz8GGsfO42Swk1+h2MOUJa0aujoK2D4fbDgP/D0ybDpa78jMiaYcnqdQesLnmX3d5+Qd28vdi56y++QzAHIklYSDLkOznoWNiyGR/vDt+/4HZExwdRowI9pf93npOW0YO2joyiY+Q+/QzIHGEtaSdJrLFw5H3r9GNof7YbZeS5jqi6rw1G0v342zU6dQM6RrreTdYk3YZa0kqhxWzjjMWjQFEr2wD+Og7d+A9tW+R2ZMcGSltGAFqfdSahxLqpK/gPDWPv4GLuDhrGkVVv2FEJuD5j1CPz9MJh8Eaxf4HdUxgRQaTHZPU9n9/IPWXP/QPIfHsnubz/yOyrjE0tataRhczjzKbhmKRx9JSx+FR7pCys+8DsyY4JFQhk0H3kLHW9bSYsf3EdR3lzW/O14di58E8AeKpli7OLiWta0I4y8H074HXz5LHQc4oa/dxvkfwE9fgiHnQ7ZrXwN05gDXlqDxjQ75UaaDL2a7TOfpEE3d0HytnfvZ+dXr5HT7zyye4wivVUXRMTnaE1tsaRVR7JbwjG/2vc5s7FrLvzG7SzS6nA47Acw/B73WdXdpNcYs7+0zGyaDr2m7HMopyUlu7awadLVbAJCTdrSsMdIWv/kKcAdiVkSqz+qlLREZCTwNyAEPKmq90WNF2/8aGAncJGqzq2srIi0AP4DdAJWAD9S1S3euJuAS4AS4BpVfbtaa3kAGnIdHPsbyJ8D306H1Z/Ars37xj/UCzKzXTJr1QNyD4c2vaFFV/9iNsljdSl5Gg++mMaDL2bP2sXs/mYGu5Z/CMVFZePX3D+Q0qJCMtv2JOOgI8hs04PMtkeS2e5IH6M21ZVw0hKREPAQcAqQB8wSkSmquihislFAd+81CHgEGBSn7HhguqreJyLjvc83isgRwFigJ9AOeEdEDlXVetORXATaDXCvSFoK3Ua4675WfghfPe+GD/wljPqrq4+P9IUm7SA71x3FZbeCrsOh43FufN5nkNUYsppAZiPIyHavtFCdr6aJYnWpdmQe1IPMg3rQ5Pgr9hue3fMHFOXNpej7eRR++TKoktP/fNpc+BwAeff1QbIaEWqUW/Zq0G0Y2T1GoKWl7F42g7QGTUhr0BTJykEysknLaoSErKHKD1XZ6gOBZaq6HEBEXgDGAJEVbQzwjLozo5+KSDMRaYvb86uo7BhgmFf+aWAGcKM3/AVVLQK+E5FlXgwzq76awSJpMPIv+z4XbYdNSyGrqfu8dye07Qvbv3dNjDs3ws5NLjl1PA62r3F354g26q8u8W1YBE+PgFAmpGdBKAvSM+HE26H7KDfPqTdAWoZLchJyf4+9zl2Dtm4+fPo3IM3FGn4NvtodGeZ/AfOecUlZ0gBx7wdd5c7xff85LJrsNX96rTYiMPgayGkNqz6GZZHHAd40x/7GXU6w4v3YT5A+fjykN3AXd6/+pPz4oRNcPF+/CWvm7D8uLR1OuMm9XzwZ1kX19Oz946Qe5VpdqkPNR91a9r50zy6KNy5zX35AS/aS2a43xQVrKN60nKKVn1GyYyNNS4td0iraQf6D5StTs5G/p8Xo2yguWEvefb2RjAZIepZ7hTJpdtL1NBrwY4q3rGbD8z+HtHQkLVT2t8lxV9Dw0JPYu3E5W966w42LqExNjvk5WR37s2f91xR8+BDhSiTe38bH/JzMg3qwJ38h2z9/puxcQtn4Yy8lo2UnivLmUThv0r7AvemaHHcl6U3bsnvl5+xc8Pr+KydC06HXEsppwe7lH7NzybTy63/y9aRlNWLX1++x+9vyvcuajfgdEspg56K3KFr5WdnwtIbNaDrs2oT/d7FUJWm1B1ZHfM7D7QHGm6Z9nLJtVDUfQFXzRaR1xLw+jTGvckTkUuBSgI4dOya4OsGR1Xj/o7GGzeGcf+8/TWkJhPebG7WBcVOhqMC99uyAvbvh4GPd+Iwcd/6spMhdT1ayxx2dpTdw40v2wu5tULp333xLS918wCXI72a4I8LSEvcXhd7nu/FbV7hOJ6g7Nxcef+R5LmmtXwSfPeCGgZsGoM9PXdJaMxs++qO3YhEdwwb8wiWtlR/CB/dEbSSFY37t1mH5O/DJ/eW349AJ7u/Xb8Ccx/cfl95w/6QVProNO/iYpCatA7Iu1fd6BJCW2ZDMdr3KPksog9bjnt1vGlUtu5GoZDak7dUzKN1dQOnubeieQkr37KRBp8Fe+Uxy+pyNFhdBcRHqvSQzx82rtJjSXdvQ0uKyyqQlxZTs3AJA6e4Cdn/zHhquTLgKk91jJFkd+1OybQ07Pn/axaRKuFI1PPxUMg/qwd6N31Lw/t/2xe1VmOwjRpPRshN78hew9Z37witWto45fc4mvWlbilbNZuvUuyJXHoDGR48rS1pb37q93HZsctwVLml98y5b376r3PimJ//WJa3Fb5XFB5DeolONkxaqmtALOBfXfh7+fAHwQNQ0bwDHRXyeDvSvrCywNWoeW7y/DwE/jRj+D+DseHH2799fjalMaWnsV6LjowGzNcF6pAGpS1aPTCJKS0srfFU0vjKJ1KWqHGnlAQdHfO4ArElwmsxKyq4Tkbbq9gzbAuursDxjqixeR7I66GhmdcnUC/F6ZdZGr82qXFw8C+guIp1FJBN3YndK1DRTgHHiDAa2qWuuqKzsFOBC7/2FwGsRw8eKSJaIdMadkP68iutnzIHI6pIx1ZTwkZaqFovIVcDbuK62T6nqQhG53Bv/KPAmrovuMlw33YsrK+vN+j7gRRG5BFiFa/7Am/eLuBPMxcAvNYHeTnPmzNkoIisrGN0K2JjoOh+Aghx/kGOHyuM/pCozCkJdilOPINj/zyDHDvU7/rh1STSFboEiIrNVdUD8KQ9MQY4/yLFD8ONPtiBvjyDHDha/3XvQGGNMYFjSMsYYExiplrQejz/JAS3I8Qc5dgh+/MkW5O0R5NghxeNPqXNaxhhjgi3VjrSMMcYEmCUtY4wxgZESSUtERorIUhFZ5t39+oAmIk+JyHoRWRAxrIWITBORb7y/zf2MsTIicrCIvCcii0VkoYhc6w0PxDqISAMR+VxEvvTiv90bHoj4a5PVpbpj9Si2ep+0ZN+jHEYBRwDni3tUw4FsIjAyalj4sRPdcfehO5B/MIqB61S1BzAY+KW3zYOyDkXASaraBzgKGOndlSIo8dcKq0t1zupRDPU+aRHxGAhV3QOEH+VwwFLVD4DNUYPH4B43gff3zLqMqSpUNV+9Bxaq6nZgMe6u4oFYB+/end497cnwXkpA4q9FVpfqkNWj2FIhaVX0iIeg2e+xE0DrONMfEESkE9AX+IwArYOIhERkHu6ms9NUNVDx1xKrSz6xerRPKiStWLcZtn7+dUBEGgEvA79S1QK/46kKVS1R1aNwd0QfKCL2bHarS76werS/VEha9eWxDOvEPW4C2f+xEwckEcnAVbR/q+or3uBArQOAqm7FPQF4JAGMP8msLtUxq0flpULSSuQxEEFQ0WMnDjgiIrgHDS5W1b9EjArEOohIrog08943BIYDSwhI/LXI6lIdsnpUgXhPiawPL9wjHr4GvgVu9jueBOJ9HsgH9uL2bi8BWuJ62nzj/W3hd5yVxH8crtloPjDPe40OyjoAvYEvvPgXALd6wwMRfy1vG6tLdRe71aMYL7uNkzHGmMBIheZBY4wx9YQlLWOMMYFhScsYY0xgWNIyxhgTGJa0jDHGBIYlLWOMMYFhScsYY0xg/D+KDcdkhp7kUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(2,2) # define figure and axis, including the number of rows and columns for the subfigures\n",
    "cmap = cm.get_cmap('gnuplot') # choose a colormap\n",
    "# for each axis, define the variable, all in dashed lines, and choose a color from the colormap\n",
    "ax[0,0].plot(gov_spending, linestyle = 'dashed', color=cmap(0))\n",
    "ax[0,0].set_title('Government spending') # set title for each subfigure\n",
    "ax[1,0].plot(output, linestyle = 'dashed', color=cmap(0.25) )\n",
    "ax[1,0].set_title('Output')\n",
    "ax[0,1].plot(inflation, linestyle = 'dashed', color=cmap(0.5) )\n",
    "ax[0,1].set_title('Inflation')\n",
    "ax[1,1].plot(interest, linestyle = 'dashed', color=cmap(0.75) )\n",
    "ax[1,1].set_title('Interest rate')\n",
    "fig.suptitle(f'Impulse Responses to a shock in $g_t$ of size $\\epsilon_0^g = {epsg}$', fontsize=16) # title of figure\n",
    "fig.tight_layout() # makes sure the labels from the figures do not overlap\n",
    "plt.show() # show the figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #6 Interpretation of impulse responses\n",
    "\n",
    "Generally, the autocorrelation in government spending shapes the agent's expectations such that government spending today affects to expected future government spending.\n",
    "Looking at the graphs, one first notcies that government spending reacts one-to-one to a shock in government spending in the first period, which it should be by construction, and then steadily decreases to zero. Output reacts positive to a shock in government spending which is not only influenced by the relation of government spending and output but also by how the other variables react to the shock. However, as we see, output is positively influenced which also makes sense intuitively, i.e if consumption (in form of government spending) increases, total output should increase as well. Especially, because of price stickiness output is mainly demand-driven. Next, also inflation reacts positively. A higher demand results in higher prices which raises inflation naturally according to the Phillips curve. Finally, also the interest rate is increased by a shock in government spending. The interest rate, as determined by the Taylor rule, is a function of output and inflation. So when both, output and inflation, are increased, the interest rate should react in the same way. Also, as we set $\\phi_{\\pi}$ larger than $1$ we ensure that the Taylor principle holds, the nominal interest rates increases more than inflation such that monetary policy can steer the economy, whereby we rule out sunspot equilibria. All in all we should also note that all variables are converging back to their steady state level after around 10 periods such that the impact of a shock in government spending is fading out quite fast."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #7 Linear Time Iteration: Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a function using the matrices A, M and D as input as well es a first guess an F\n",
    "def lin_time_iter(A, M, D, Fguess):\n",
    "    '''\n",
    "    Performs linear time iteration given the arguments. \n",
    "    Returns F and Q.\n",
    "    '''\n",
    "    # make a while loop as before until the difference between F and Fnew is small enough\n",
    "    F = Fguess\n",
    "    Fnew = np.linalg.inv(A - M @ F) @ D\n",
    "    while np.max(abs(F - Fnew)) > 1e-10:\n",
    "        F = Fnew\n",
    "        Fnew = np.linalg.inv(A - M @ F) @ D\n",
    "    # save the latest F when the loop stops\n",
    "    F_final = Fnew\n",
    "    # Calculate Q\n",
    "    Q = np.linalg.inv(A - M @ F_final)\n",
    "    # function return the solution for F and Q\n",
    "    return F_final, Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #8 Define Matrices A, M and D\n",
    "\n",
    "$ A z_t = M E z_{t+1} + D z_{t-1} + u_t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = A # same as before\n",
    "\n",
    "M = np.array([np.zeros(4),\n",
    "            [-params['gbar'],1, (1-params['gbar'])/params['sigma'],0],\n",
    "            [0,0,params['beta'],0],\n",
    "            np.zeros(4)])\n",
    "\n",
    "D = np.zeros((4,4))\n",
    "D[0,0] = params['rho_g']\n",
    "\n",
    "# guess a 4x4 matrix\n",
    "F_init = np.random.randn(4,4)\n",
    "\n",
    "# call the function and save results F and Q\n",
    "F, Q = lin_time_iter(A=A, M=M, D=D, Fguess=F_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.7        0.08116102 0.03677436 0.06327765]\n",
      "[1.         0.11594431 0.0525348  0.09039664]\n"
     ]
    }
   ],
   "source": [
    "print(F[:,0]) # should be equal to cg\n",
    "print(Q[:,0]) # should be equal to ce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #9 Interpretation of $F$ and $Q$\n",
    "\n",
    "We should first note that the first columns of $F$ and $Q$ are similar in the way that their elements have the same ratios, i.e. the second element in $F$ is $0.7$ times the second element in $Q$ and so on. The first column of $F$ tells us how $z_t$ reacts to a change in $g_{t-1}$. The first column of $Q$ shows the reaction of $z_t$ to a shock in government spending, the second column would show the reaction to a shock in output and so on. Similarly, the vectors $C_g$ and $C_{\\epsilon g}$ show those reactions and are the same as the first columns of $F$ and $Q$, respectively."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "98f8afcc7b498ff748bcf3669bc6fa2980d8228748497c4de9780cb8988f727b"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 64-bit ('compMacro': conda)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}