problem_sets/problem_set_3.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Assignment 3\n",
        "\n",
        "Jesse Perla, UBC"
      ],
      "id": "49a45037-1672-445f-a8fa-5274a8ea0e82"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Student Name/Number: (doubleclick to edit)"
      ],
      "id": "e4f30de8-4075-42fd-9463-22fe385c2d95"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Instructions\n",
        "\n",
        "-   Edit the above cell to include your name and student number.\n",
        "-   Submit just this `ipynb` to Canvas. Do not rename, it associates\n",
        "    your student number with the submission automatically."
      ],
      "id": "36421aa7-eb5e-4d9d-bb15-adc81f73a353"
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "using Distributions, Plots, LaTeXStrings, LinearAlgebra, Statistics, Random"
      ],
      "id": "6e7c826c"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Question 1\n",
        "\n",
        "Following the notes on [AR(1)\n",
        "processes](https://jlperla.github.io/ECON408/lectures/stochastic_dynamics.html#/title-slide)\n",
        "rather than plotting the distribution as normal instead lets see what\n",
        "the stationary distribution looks like with simulation.\n",
        "\n",
        "### Part (a)\n",
        "\n",
        "From $X_0 = 1.0$ simulate up to $T=1000$ using the process\n",
        "\n",
        "$$\n",
        "X_t = a X_{t-1} + b + c W_t\n",
        "$$\n",
        "\n",
        "Where $a=0.9, b=0.1, c=0.5$."
      ],
      "id": "6b2e1ca0-4c28-4ee1-82c9-e27167123301"
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "0.5"
            ]
          }
        }
      ],
      "source": [
        "T = 1000\n",
        "X_0 = 1.0\n",
        "a = 0.9\n",
        "b = 0.1\n",
        "c = 0.5\n",
        "# Add code here"
      ],
      "id": "366844e3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (b)\n",
        "\n",
        "On the same graph plot the histogram of those simulated values (i.e.,\n",
        "$\\{X_0, \\ldots X_T\\}$) then plot the density of the stationary\n",
        "distribution calculated in closed from in [those\n",
        "notes](https://julia.quantecon.org/introduction_dynamics/ar1_processes.html#equation-mu-sig-star)\n",
        "(i.e. create a normal distribution with $\\mu^* = b/(1-a)$ and\n",
        "$v^*  = c^2/(1 - a^2)$\n",
        "\n",
        "Hint: `histogram(X, normed=true)` normalizes the empirical draws so they\n",
        "are a proper PMF."
      ],
      "id": "45ac1e90-16d4-4bc3-9399-58639921ad9c"
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Add code here"
      ],
      "id": "890dfaa6"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (c)\n",
        "\n",
        "What happens if you discard the first 200 observations from that\n",
        "simulation (i.e., $\\{X_{199}, \\ldots X_T\\}$))?"
      ],
      "id": "37d91567-1d16-40be-83a8-f1dc6d99d685"
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Add code here"
      ],
      "id": "4fe9e39e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Do these line up approximately? Explain why it may be better or worse?"
      ],
      "id": "239581cc-954b-47be-8eca-cd16d58afc2a"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "(double click to edit your answer)"
      ],
      "id": "cb44b1d3-e12b-4660-8ee3-e3d75e0a07bd"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reusable functions, do not modify without clearly annotating your\n",
        "changes."
      ],
      "id": "8601ded9-495f-4f50-9b34-df6fe238b642"
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "gini (generic function with 1 method)"
            ]
          }
        }
      ],
      "source": [
        "function lorenz(v)  # assumed sorted vector\n",
        "    S = cumsum(v)  # cumulative sums: [v[1], v[1] + v[2], ... ]\n",
        "    F = (1:length(v)) / length(v)\n",
        "    L = S ./ S[end]\n",
        "    return (; F, L) # returns named tuple\n",
        "end\n",
        "# Assumes that v is sorted!\n",
        "gini(v) = (2 * sum(i * y for (i,y) in enumerate(v))/sum(v)\n",
        "           - (length(v) + 1))/length(v)"
      ],
      "id": "9fa6df05"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Question 2\n",
        "\n",
        "You can create a Pareto distribution with tail parameter `a` and draw\n",
        "from it with"
      ],
      "id": "2cc544fc-9f0e-45b3-8ee9-053da85db356"
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "3-element Vector{Float64}:\n",
              "  1.3742764221615884\n",
              " 13.637082808351632\n",
              "  1.0737073421937346"
            ]
          }
        }
      ],
      "source": [
        "a = 1.5\n",
        "d = Pareto(a)\n",
        "rand(d, 3)"
      ],
      "id": "31f6493b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The gini coefficient of a Pareto distribution is given by\n",
        "\n",
        "$$\n",
        "G = \\frac{1}{2a - 1}\n",
        "$$\n",
        "\n",
        "### Part (a)\n",
        "\n",
        "Draw `N=1000` observations from the Pareto distribution for `a = 1.5`\n",
        "and calculate the gini coefficient. Compare to the closed-form solution"
      ],
      "id": "66dcf2b1-8325-4f80-bb8f-5d6108343f4c"
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "1.5"
            ]
          }
        }
      ],
      "source": [
        "N = 1000\n",
        "a = 1.5\n",
        "# your code here"
      ],
      "id": "9fa68379"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (b)\n",
        "\n",
        "Given a pdf $f(x)$ the cumulative distribution function is given by\n",
        "$F(x) = \\int_{-\\infty}^x f(y) dy$. The complementary cumulative\n",
        "distribution function (CCDF) is given by $F^c(x) = 1 - F(x)$. It starts\n",
        "at $1$ at the minimum of the support of the distribution and then drops\n",
        "to $0$ at the maximum of the support.\n",
        "\n",
        "Take the following code which calculates the CCDF for a distribution and\n",
        "a range of values `x`.\n",
        "\n",
        "For this, plot the relationship between the `log(x)` and\n",
        "`log(ccdf(d,x))` for the distribution below (this is called a [log-log\n",
        "plot](https://en.wikipedia.org/wiki/Log%E2%80%93log_plot)). Describe\n",
        "what you know about this relationship?"
      ],
      "id": "d75b0689-d88c-4610-b312-4fbd0cd6f86b"
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "1.0:0.09090909090909091:10.0"
            ]
          }
        }
      ],
      "source": [
        "ccdf(d, x) = 1 - cdf(d, x)\n",
        "\n",
        "N = 1000\n",
        "a = 1.5\n",
        "d = Pareto(a)\n",
        "max_x = 10.0\n",
        "num_points = 100\n",
        "x = range(support(d).lb, max_x, num_points)  # from the lower bound of support to max_x\n",
        "\n",
        "# edit your code here"
      ],
      "id": "100e9337"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (b)\n",
        "\n",
        "For the above case, numerically calculate the slope of this line\n",
        "(easiest is just rise over run). Compare it to the tail parameter `a`."
      ],
      "id": "dc919283-51af-43f9-883f-304b6337878e"
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [],
      "source": [
        "# your code here"
      ],
      "id": "a21bf42b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (c)\n",
        "\n",
        "Do we will do the [empirical\n",
        "counterpart](https://en.wikipedia.org/wiki/Empirical_distribution_function).\n",
        "Given unweighted vector $\\{X_n\\}_{n=1}^N$ of observations, we can define\n",
        "the empirical CDF as\n",
        "\n",
        "$$\n",
        "\\hat{F}(x) = \\frac{\\text{number of observations } X_n \\leq x}{N}\n",
        "$$\n",
        "\n",
        "With the equivalent CCDF as $1 - \\hat{F}(x)$.\n",
        "\n",
        "The code below calculates the empirical CDF for some simulated data.\n",
        "Adapt it to (1) draw `N` draws from the Pareto distribution above with\n",
        "`a = 1.5`; (2) plot the same log-log plot as above; and (3) overlay the\n",
        "plot with the theoretical log-log plot from the previous part to compare\n",
        "the two."
      ],
      "id": "83f64b83-0227-4de8-8b0a-39ff40b0e197"
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/html": [
              "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoUAAAHMCAIAAAAs9EW6AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd3wUdf7H8c/M7iYhHUhoAYMBEpEOAURQWhRRRIqIh2c5CweW8zy7d3r+Tk/xDgsgenqCJ6LYUA9FepMmIr2ISAskJEAgbbOb7O7M9/fHehGREiDJbHk9//Cxu8zufjLOd977ne/MdzSllAAAAEvpVhcAAADIYwAAAgB5DACA9chjAACsRx4DAGA98hgAAOuRxwAAWI88BgDAeuQxAADWC5Q8LiwstLoEAAAsExB5XFBQcMkll9T+97pcrtr/0kDDSvB4PD6fz+oqLMZmoJRyu91WV2Ext9vNDMoWtoWAyGOfz2fJDtE0zdr/0kDDSlBKsQ9iMxBWAmtARCxdCQGRxwAAhDnyGAAA65HHAABYjzwGAMB65DEAANarzjwuLi7+9ZlpZWVlHo+nGr8FAIDQU6U8nj59ep8+fRo0aPCPf/zjpAvs3LmzXbt2HTt2bNas2bx58/wvlpWVDR48uFWrVikpKad6IwAAkCrmcePGjf/yl7/06dPnVNfL/+EPf7jpppv27t07ffr0W2+91d8hfuWVV5RSubm5mzdvfvHFF7du3VqdhQMAEEKqlMf9+/fPysqKi4s76b8WFBQsWrRozJgxItK3b9/69esvXrxYRGbMmDFmzBhN0xo3bjxkyJAZM2ZUY90AAIQS+/l/xIEDB+rWrZuYmOh/2qJFi+zsbBHJzs5u0aKF/8W0tLTT94/Ly8sXLlxY+bRr164JCQnnXxsAAGfLVFLsERHxmuL0KRGpMMTlExFx+6TckA71teSoav7Sashjt9sdERFR+TQyMtLlcimlysvLK1+PiooqKys71Sd4PJ6SkpLXXnut8pX777+/S5cu51/b6Z2mpPDBSqioqNB13eFwWF2IldgM/PNXh/nMqS6XyzRNTdOsLuQcOX2a2yc5LslxafvLNKdPK/MpESn2aErEbWgVhhhKSn0iIiUeMUWrTNkSr/jPRi7yRMbaPXZdRCTBIZqIpklihIiILio+QkQkQpex6Ub/Rmcxs6au69HR0adfphryuFGjRseOHTNNU9d1ETly5Ejjxo01TWvYsGFBQUFaWlrli6f6hIiIiAYNGnz66afn8O3l5eUHDx48t8pdLhe7oRpdCRdeeGHgt22Hw0Eei0hsbKzVJVhJKaXrekxMjNWFWEnTtOjo6EBrs06vLMg19zsl26ny3HK0XBlKij3iNcXplXJD3IZy+SRSF5sudSO0xtGSGqulxkpClFbPJjF2EZHECNE0ibZr/sXiHSIi8Q6x6RKpS7RdE5EYh0To4nQ6rWoL557HHo/H6/XGxMSkpqYmJSWtWrWqV69eJSUl69at69Gjh4j06tVr8eLF3bp1E5HFixffd9991Vb1cSZMmDBp0qRGjRqdw3uVUoG25dW+mlsJu3btmjVr1uWXX14THw4g2FUYkutSB8skp0zlueWwW7l84jakqELKfMrtkxKvlHrlWIW6vJGekSCtErR+TbQYh8Q7xKFLrEOibFLHpvlzNARUKY8XLFjwySefrFy5ctOmTQcPHvz973/fuXPnqVOnTp8+fcWKFTab7YknnrjtttvGjBkza9asESNGpKamisjDDz985ZVX+ny+H3/8saSkZPjw4TXxB3i93rvvvvuJJ56oiQ/H+Rg8eLDX67W6CgAW8I+/FnpUiUdKvFLsUYfckueSw26V55bcMpXtVGVeaRarNY2RxtFasxhJjtJiHZIQIXVsEm3X60ZKHbtE2yS5jhZTDUdyg0CV/srU1NSsrKysrCz/0wYNGohIVlZWy5Yt/a+MGTMmIyNj1apVd99994gRI/wvdunSZeXKlZ9//nmPHj0mTpx4/BgzACDYHXbLsQp1rEKOVUihRx0rl1yX+myfyncrmybxDi0hQuIjJCFCEiO05ChpWEfrnKQ1qKM1i5ELYrVEMuGXqpTH6enp6enpJ7zYsmXLyjwWkb59+/bt2/eEZS666KLHHnvsPEsEAASCIo+MXWn8UKQOuaXEq+pGaMl1pH6k1I3U6kVKvUhpEq19OUDPSAj3ccBzEx5HAQAA5yS3TC3PV+sK1I5itemoDG2u/bWvrVE0vdvqRx4DAH7mM2XjMbXxqNpZrLYcU1sKpW9jrVuyNqCpnpEozWLo+9YU8hgAIEcr5G8bbOsKjW2F6qJErUuS1jpRuyJFv7yRFmmzurjwQB4DQLj7dJ/5t/Vmr2R5qbveob5OAFuCPK5Zc+fO3b9/f+XTjh07+i/IPgfr16+fNm3aK6+8csLrn376aXZ29gMPPFDFzykvL7/22msXLFjw63/atm3bV199VVpa2qJFi6FDh8bHx+fk5Hz11VcikpCQcMEFF7Rv375ywoTvvvtu/fr1le+99tprTzPlC4BAo0Q+3mP+N1utK1DxEfJ/XfT+SUZMwE0HEkZC4iLqADZx4sSPP/545/8UFBSc80fFxcX9+ix3EWnQoMEFF1xQ9c8xDGPp0qW/fv2FF1647LLLCgoKGjRosGbNmjZt2ojItm3bnnjiiT179qxevfqZZ55JSUl58803/ct/9tlnkydPrvzTTnXvLwABqMwnIxYZk7abQ5prc66yfXud/bpUnSC2Fv3jGjds2LCxY8ce/8q6desSEhJ27969devWgQMHXnzxxQsXLty6deuAAQNat24tIlu2bNF1/fDhwxs3buzdu3fnzp1FJCEhwZ+RTqdz2bJlXbp0mTlzZufOnRs1alTZZ62oqJg9e3Z2dvaFF144cODAyMjILVu2rFy50uv19urVq1OnTqcqcs2aNU8//fTGjRszMjL8rzz44IP+B8nJyePGjfM/Xr58ef/+/Tt06NC9e3cR6d69+/jx46t3dQGoBX9bb0Tqsvhqu4NOWcDgf0VtUP/jf/rGG28MHTr0yy+/LCkpueyyyx566KFPP/20uLj4sssuy8nJEZEZM2YMGzZs2rRpIjJ06FD/zN6bN29+8sknRSQvL+/mm28eNWpUUVGRy+X66quv3njjDREpLS3t2rXrzJkz7Xb7/Pnzd+3aJSJ/+9vfnE6nx+MZOXLkzJkzT1Xh7NmzBw8eXBnGIlJ5b67jXXbZZddcc81HH31U+Xf5fD6fz2cYRjWtKgA1osQru0rUmsNq9gE1ZoWx+KAa391GGAeU0Owf/+5rI6estu/ToolMvdzW9FcXA/zpT3+qnBTl3//+9w033CAil1xyyaRJk0Rk+/bthw8f9kfvzp0758yZc9ddd4lIkyZN3n77bRFp06bNgw8+OGzYsOM/s7i4+J133mnWrJmIbNu2zf/ipEmTLr744vfee+/4JT/++GMRKSsra9GixeTJk081a+muXbtOGsC/lpGRsWfPHv/jd99915/xKSkpW7ZsqcrbAdSOfaVq/VG16ahafVitPaKi7Vr9KKkfKclRWuNomT/QXjfS6hLxS6GZx4930EtrfeJkTaRx9EnGX1566aUTjleLSLt27fwPGjZs6J/u2//46NGj/seVp3117959x44dJ3RAGzZs6A/j423atOmEKdJM0/zTn/40Z86cunXrappWVFR0quITExNLSkpO/wf6FRYWVt417JZbbqkcTgYQIKb9aD7yrREfoXWur3Wsrz3WQe/RQKsTmjv7kBKa/4vSA362Nv+9KUVE07TjH1ce0y4tLa18EBMTY7P94vqDk94cMCEhobi4+PhX5s6du2bNmu3bt9tstmXLlo0ePfpU9XTs2HHChAmVN808lYqKioULFz7yyCNn+PMA1C5DyQe7zVn71Q9Fan+ZWnaNvV29QN8N4gSMHtS48vLy4v+p+knIn376aV5enoi8+uqrAwcOrMpbhg8f/uabb+bm5opISUlJaWlpeXm5P1+9Xu+ECRNO895Ro0Z5PJ777rvP/zsgPz9/zJgx/n8yDKOwsPDQoUOrVq0aMmRITEzMzTffXMW/AkCNKvbI1/nq5a3m5V/6pvxgjkzTPu5vO3yTgzAORqHZPw4cDRo0ePnll19++WX/0+HDh7/88stJSUkJCQn+V+rXr5+YmOh/XK9evbp16/ofX3nllUOHDj18+HBqaur7778vIlFRUQ0bNhQRh8Nx/AVOcXFx9erVE5EBAwY8+OCDl112md1ut9vts2bNGjRo0H/+85+0tLTY2NibbrrJf0Ra1/VfXx8VGxu7bNmyBx98MCUlJSoqKi4u7o477hCR6OhowzAyMzOjo6ObN2/eu3fve++913+8OjEx0TTNmlpxAE7raIVM3m5O3Gp0rK91qq893dnWr4lmI4WD2c8HSC2Un5/fs2fP3bt3n8N7n332WV3XQ+z+x0888URERMTTTz99xgPIJ2UYxvHHt094erZvP2eDBw++//77+/fvf/4fVaMqKip0XT/pEED4cDqdsbGxVldhJaWUy+WqvHQwYH17RM3KNjcdU98dUSPS9Ifb69U4oXRZWVl0dLhPCGJhW6B/HNDOIYxF5IQ0PdtwrZYwBlAtcsrUwlw1a7/aXaL2laoL47SRLfT72uhdkzROkA4x5HEguv32288tiQGEjPd2mfevNhIitB4NtdvT9VYJ0jyWWzuEMvI4ELVs2dLqEgDUqnUFak+pOuSWgnJ1pFy+O6LchizhNOlwQh4DgAWUyA9Fam2B+uawWpSrYhzSrq6WHCWNorW0OBnVQr+0YXgP5IYf8hgAatsLm8wJ24ykKK1rknZJA+0PbfSMgJ81ATWNPAaAWrXlmHpxi7FxmL3JyWb0Q9gKhTzeunWrf5ZmBJSDBw9aXQIQcPLdctNS4x/dbIQxThD0edyrV6/NmzefWx5X14W2Qa3mVkJGRsZJb9gMhCGfKd8eUW/uMOfnmmNa225L5wIKnCjo87hPnz59+vQ5t/cyB4KwEoCaV+yRGxb5Cj2S1UT7YYQjLqznnsEpBX0eA0CAG7bQ1zVZey7TpnOIGqdGHgNADTpQpn4oloVXM7c0zoAxDACoKaaSx9ea91xMxxhnRh4DQPUr9croFUaj97zHKtQf27CnxZlxvBoAqpnXlBGLfMlR2g8jHNz1AVVEHgNANXt/t+n2yX962xg0RtWRxwBQnbYXqX9uNv/VkzDG2WFUAwCqzXcF6oqvjFtb6Zc2JI1xdugfA0D1mHNA3bnceOsy28BmhDHOGnkMANXA5ZN7Vxnv9LZlpRDGOBfkMQCcr0NueW6j0a2BRhjjnDF+DADnZUGu6vSZV4mM68oeFeeO/jEAnKN9peofm83Ps83ZA+yd6tMzxnnh1xwAnDWfKQ+uMXp+YTSL1bZf7yCMcf7oHwPAWVtXoObnqM3D7fWZfgvVhP4xAJwdQ8m4TebINJ0wRjUijwHgLBx2S5/ZPqdPPdye/SeqE9sTAFRVqVe6z/Jd3Uyfd5U90mZ1NQgtjB8DQFW9scPMaqI93oGeDKofeQwAVZLtVOM3G0uuYbeJGsGvPAA4M5dP+s42nu5sa53IpU2oEeQxAJzZh3vMdvW0Ma3ZZ6KmcOAFAM5gR5F6fK2xlCPVqEn81gOA0znoUgPmGhN72C7iSDVqEnkMAKfkM2XEIuOR9voNaewtUbPYwgDg5Mp8cu8qI9ou91zMrhI1jo0MAE7ioEv9cbVR6JFPsxg2Rm1gOwOAX5j2o/nYWiPKprWpK89m2uIcVheE8EAeA8DPthepv3xnLrnGnpHA2VuoVRyvBoCfTfvR/E0LjTBG7aN/DAA/+bFYvbnD/P56jlDDAvSPAUBEZI9T6/mFb/KltoZ1rC4FYYn+MYBwV+yRKT+Y/9pu/2tn229a0EuBNdjyAISv7UXqzuVG+sfeXSXyYqbBdcawEP1jAGFqfq6642vj0Q76C90c9SKUy2VaXRHCGnkMIBzllKm7Vxr/vsx2VVNNRJSyuiCEPQ7OAAhHb+4wB12g+cMYCAT0jwGEnT2l6o0d5uZhXNeEAEL/GEB42XJMXTrL949uXNeEwEIeAwgvz28yH25vu7UVez8EFrZIAGFk8nZzyzH1hzbs+hBwGD8GEC4Ou+XZDcaGYQ4HcYzAw1YJICy89YN5+Ze+Ma1tjRg2RkCifwwg9O0uUY98aywYaO+SxAVOCFD0jwGEuB+K1WVf+iZfaiOMEcjIYwAh7qUt5gNtuVEEAh3HqwGEsmMV8tEe88BvmPoDgY4fjABClqnkhkW+sa31WOIYAY88BhCavKYMW2hE27W/dbFZXQtwZhyvBhCaVh5SB11q1bV2O/0OBAO2UwAhaEW+uuNr45H2OmGMYEH/GECoyXPJrcuMyT1t3E4RQYQ8BhA6fKb8aY3x6T51Yxr3NkaQIY8BhI7Pss01h9X26+3xnFCNYMPQCoDQsfaIGtpcJ4wRjMhjACFiX6masVvdkMZhagQl8hhAKKgwZOA84/+66Glx5DGCEnkMIBR8tNdMi5Pb09mnIVix7QIIBTP3qltasUNDEGPzBRDc3D554Btja6Ea0JQdGoIYmy+AIFbkkY6f+Y6Wy6Zh9sQIq6sBzgPXHwMIYrOyzcwkbVof7hiBoEf/GEAQW3VIXZHCCdUIBeQxgGA1P1fNPqC6JpPHCAUcrwYQlP71vfnKVvOd3rY2dcljhALyGEDwOeSWv3xnrLjWflEiYYwQUdU8fu211yZPnqxp2h/+8IfRo0ef8K933nlndnZ25dPf/e53o0aN+uabb5588snKFz/44IP69euff8UAcOks3z0X64QxQkmV8njZsmXPP//8woULlVJZWVkXX3xxr169jl/ggQceKC8vFxG3252VlfX888+LyJEjR5xO56uvvupfJj4+vrqLBxCOXv/ebBwt/9eFc6oRUqqUx2+99dbo0aMzMjJE5M4775wyZcoJedymTRv/g/feey89PT0zM9P/ND4+vkuXLtVaMICwpkT+vtFcPogwRqip0vnVO3fubNu2rf9x27Ztd+7ceaolp0yZctddd1U+Xb16dYMGDTp06PDaa6+d/itM0yw8jmmaVSkMQLi5Z6XRqb5cyE0jEHKq1D8uLi6OiYnxP46NjS0qKjrpYnv37l29evVHH33kf9qpU6e1a9empqauX79+xIgR9erVu/HGG0/6Ro/Hk5ub27lz58pXXn311d69e5/F33FOysrKavorAh8roaKiQtd1hyOsb5kbLJvBd0f1uQfsawd6nM5q/mSllNvtVkpV8+cGFZfLZZqmpoX1b50aagu6rkdHR59+mSrlcXJycmUGFxUVJScnn3SxqVOnXnfddUlJSf6nTZs29T+49NJL77777i+++OJUeRwREdGsWbPdu3dXpZjqFRsbW/tfGmjCfCU4HA7yWIJhMzhWIU9u8f25k14/ofonxlRK6bpe2fEIT5qmRUdHh3kei3VtoUrHqzt27LhmzRr/4zVr1nTq1OnXy5imOW3atDvuuOOkn3D06NEw39ABnKfRK4yeDbXbM5jFCKGpSv3jsWPH9urVq1OnTkqpadOmrVq1SkR+/PHHK664Yt++ff5l5s2bp5Tq169f5bsmTpwYHR3dtGnTdevW/fvf/160aFEN1A8g9BV75K/rjR+L1VuX2cO974bQVaVfmhdffPEXX3yxcOHCxYsXf/nll/4TrRMSEo4//uzxeCZOnGiz/XzSY+vWrdetWzd16tSjR4+uWrWqW7du1V49gJD38V4z/WNvqUfmD+QOTghlWiCcv5Cfn9+zZ8/aHz92Op2BP2ZW01gJnM8lAbwZzNxr3rfamHWlPTOpZjvGSimXyxXmw2plZWWMH1vYFhiJARC4nlxnzsyq8TAGAgF5DCBA7SpRh92qG7dvQnggjwEEolnZ5i1LjT+2tdmIY4QH7u8EIOD8UKzuXWVOulS/9gL6DAgX5DGAgPPJXnVLK+26VMIYYYTNHUBgKayQN743b0hj74TwwhYPILBsL1IZidK+HuPGCC/kMYDAsq1QZSQQxgg75DGAADI/V720xezXhDxG2OF8LgCBwuWTUYt9H/a39yePEX7oHwMIFDcuNkak6YQxwhP9YwABwWfKkjyz5Nawnkgc4Yz+MYCAcNcKY0BTna4xwhb9YwDWW5qnNh1Va65jj4TwRf8YgPU+3mvedZHuYIeEMMbmD8BihRXy+T41oCnHqhHWyGMAFvuuQGUma2lx5DHCGnkMwGJrj6gWcVYXAViNPAZgpa8OqNe+N29NZ1+EcMfZjACsUVAub+ww//W9+VmWrQN3j0DYI48BWOPO5UZylCy62pbO3SMA8hiAJQ6UqVWHzD0jHbHMxwWICHkMoJYdcssrW423d5p/7WwjjIFK5DGA2lNuyJVzfD0batuud9SPtLoaIJCQxwBqiRK542ujW7I2uaeNEWPgBOQxgFoy54DaVaKWX2snjIFf45o/ALXBVPK3DcYDbfUI9jrAydAyANQ4JfL69+Zht9zYgn0OcHIcrwZQ4x5ZY6w5omZdabO6ECBwkccAataqQ+rzbPXdEHtChNWlAAGMY0cAatan+8yxrXXCGDg98hhAzVqer3o05JRq4AzIYwA16PXvTUPJJQ3IY+AMGD8GUFMqDPnnZvOLK5n9Azgz+scAasqcHLN1orSpSxwDZ0YeA6gRTq+88b05iguOgaqhqQCoEX1m+xpHa8MvZCcDVAnjxwCq35s7zCibTL2cCUCAquKnK4BqVm7Iu7vMR9qzewHOAg0GQDX7fJ8Zqcs1F7B7Ac4CDQZANXtvtzn6Ip2LnICzQh4DqE6lXll7RPVuzL4FODu0GQDVadI288oUvWEdq+sAgg15DKA6rTpk3prOjgU4azQbANVGiWw8Jm2ZkAs4e+QxgGqzIl/FOYSD1cA5YD4QANUjt0z9dqnxn97MAQKcC/rHAKqBz5Sr5xl/7qj3bczBauBckMcAqsGKQ6pupIy+iF0KcI5oPADOV06ZeniNcVVT9ifAuaP9ADhfU35Q3Rtoj3ZgfwKcO9oPgPOy8pB6/Xvj9xfpjBsD54M8BnDufixWIxcbH/Szt6tHHAPnhTwGcI5MJc9sMO/M0PpwTjVw3shjAOfowz3m1kL1SHsuOAaqAfOBADgXXx1Qj3xrzh5gi2YvAlQH+scAzlqpV8asMD6/wtaeYWOgmpDHAM7O6sOq+399N7XUuiQRxkC14UgTgLOwr1Rdv9CY1sfWvwlhDFQn+scAzsLbO807MzTCGKh25DGAqjpWIVN2qhtbsN8Aqh/tCkBV3bvKuK2V1jqRzjFQ/Rg/BnBmXlNuWWbsK1XTerPTAGoE/WMAZzZ+i1lYoZYPstvZZwA1g5+6AM5gd4l6fbu5bihhDNQgmheAM5i03fxdupYcZXUdQEijfwzgDN7fZW4e7rC6CiDE0T8GcDrfHFZNorVGdayuAwh19I8BnNLuEjVsoe9fPbmDE1Dj6B8DOKXrFxnPZdoGp7KjAGoczQzAya3IV8cq5LZ09hJAbaClATi5idvM57uyiwBqCY0NwEkUlMvyfDOrCbsIoJbQ2ACcxL9/MEe20BtwWjVQW8hjACex+pC6vBH3jQBqD3kM4EQlXll1yLy6GfsHoPbQ3gCcaM4BMytFj+KqY6AWkccATjR+s3lLK3YOQK2iyQH4hbk5qtgrVzdj8BioVeQxgJ+5ffLwGmPCJRyqBmobeQzgZ/NzzWaxMpDOMVDryGMAP3tzhzkyjd0CYAEaHoCf5Jap9QXqphbsFgAL0PAA/OStH9Rt6bqdvQJgBVoeABGRwgr5Yr/ZvQEjx4A1yGMAIiIvbzW6JmtDuNUxYBHaHgARkV0l0r8JnWPAMuQxACmo0Jblqa7J5DFgGfIYCHdzDqgecx33t9VTY8ljwDL2Ki43ZcqUL7/8smHDho899ljz5s1P+NfXXntt//79/sddunQZMWKE//HHH3/8wQcf1K1b94EHHmjTpk011Qyg2uS55J5VxvSevv7NI6yuBQhrVeofT5ky5aWXXrrvvvuaNm3ar18/r9d7wgLTp0/3+XxpaWlpaWnJycn+F7/88ssHH3xw9OjRHTt27NevX1FRUTXXDuC8zckxL2+kdU8yrS4ECHdV6h9PnDjx+eef79evX79+/WbNmjVr1qzhw4efsMygQYP69Olz/CsTJkx48sknBwwYMGDAgHnz5r333nv33HNPddUN4Pwdq5BJ28znMpmtGrDemfvHPp9v27ZtXbt29T/t2rXrpk2bfr3YU089NWjQoKeffrq0tNT/yubNm8/4LgAWun+10b+JxmzVQCA4c//42LFjhmEkJCT4n9atW/fw4cMnLHP33Xc3a9ZMKTV+/PhBgwYtXbpURAoKCirflZiYuH79+lN9hcfj2b9/f926df1PdV1/5513Tuht14SysrKa/orAx0qoqKjQdd3hcFhdSG3b69Tm5Tg2XON1OhWbgVLK7XYrpawuxEoul8s0TU0L699nNdQWdF2Pjo4+/TJnzuPExERN05xOp/+zSkpK6tevf8Iyv/3tb/0PunXrVq9evT179rRo0SIhIcHpdPpfP+m7KkVERDRt2rQysO12e1xc3BkLqxaxsbG180WBLMxXgsPhCMM8ViK3LvD9tbOeUven07jCfDNQSum6HhMTY3UhVtI0LTo6OszzWKxrC2c+Xh0REdG8efNt27b5n27durVVq1anWdhut3s8HhHJyMiofNe2bdtO8y4R0XW97v/UWhgDYeudnWYdu4xtzRWPQKCoUmu88847n3nmmby8vJkzZ27ZssV/Mte8efMef/xxETl69Oj06dP37t27a9euMWPGpKWlpaen+9/1wgsv7N+/f9GiRXPnzr355ptr9C8BUEXlhvx1vflGL5se7h0hIIBU6fzqhx9+uLCwsF+/fg0bNvzvf//r778ahuG/8EnTtM8//3zcuHF2u71bt25z5syx2Wwicvvtt+fm5g4cODAhIWHGjBkpKSk1+pcAqIpspxo83xh0gda2LmkMBBAtEM5fyM/P79mz5+7du2v5e51OZ5iPmQkrIfzO5/rbBrPEo8Z3/8U1TmwGSimXyxXm48dlZWWMH1vYFhg9Ar+WdFYAABxmSURBVMLI7hL11g5zZBoNHwg4NEsgjHyWrUakadw3AghA5DEQRt6kcwwEKlomEC6+2G/GO4TOMRCYyGMgXLy9U/2pHZc4AQGKPAbCgtMrK/LN61Jp8kCAonECYeGrA+alDfWYqt7xHEBtI4+BsDA/Vw1O5Vg1ELjIYyD0+UyZudcc2JT2DgQu2icQ+rYXqfgIrfEZ7vYGwErkMRDiXD4ZusCY1IPGDgQ0migQygwlz2002tfTBnNmNRDYONsSCGUTtppL8tSsK2npQKCjlQIha1eJenGLueQaW/1Iq0sBcCYcwgJCk6Hk+Y3m9Rdq6Qlc5gQEAfrHQGiasdvcVaK+uoo2DgQH+sdACPKYMvUH8742TMgFBA3yGAg1SuSxb43ESO3aC2jgQNDgxzMQUrYXqZuXGnEOmd7HFmmzuhoAVUYeA6HjkFvu/Nq4M0Mf25qeMRBkaLRAiDhWIZfO8mUma7en066B4EP/GAgRf99o9G6sTezBQWogKJHHQCjY71Tv7zJ3jHBYXQiAc8RxLSAUfHVADbpAT4iwug4A54o8BkLBPzebtzFsDAQzGjAQ9D7Za9aPkp4NmRcTCGKMHwPBzWfKfauMT7Joy0Bwo38MBLc1R1RavEbnGAh25DEQ3GZlm4OZFxMIfjRjIIiV+WT2AdWvCZ1jIOiRx0AQe/Rbo0cDLTOZPAaCHnkMBCtTydf56r42OmkMhADyGAhKO4rU0IXGhXFau3rEMRAKyGMg+Hx7RPX/yshqos3sbyONgdDANYtAkCnyyMtbzdEX6fe14fc0EDpoz0CQyfrKVy9SHu1A4wVCCv1jIJjsKFKH3fLqpRymBkINP7GBYHLzUuOpzpxQDYQg8hgIGhuOqj2l6qYWNFsgBNGwgeCw4agaNM94vaetDqNMQCgij4FAp0QW5qp7Vxn3ttFvSKPNAqGJX9pAoPvPTvOFTebfuugjCGMgdJHHQEDLd8tfvjPnXGVrzzxcQEjj5zYQ0F7aYtzSSiOMgZBH/xgIXIaSd380v7mOdgqEPvrHQOBaclC1iNdSY+kcA6GPPAYC1/ObjLsvppECYYGmDgSoVYdUYYXcyDnVQHigqQMBasoPZr8mGnNjAmGC80SAQLThqFp9WK0fSgsFwgX9YyDgeEx5aYv5mxZ6lM3qUgDUFvIYCDh/XmuUeuWBtjRPIIxwNAwILB5T3t1lbh7miHVYXQqAWsQPcCCAKJERi4zBqXqDOlaXAqB2kcdAAJmw1Szzqn/1ZNwYCDscrwYChcsn47eYCwbauMYJCEP0j4FAMX2X2SVJa51IGgPhiP4xEBBKvfL4WmPZIJokEKboHwMB4Y+rjREX6m3r0jkGwhQ/xgGLlRtyx9fGjmK1ejDtEQhf9I8BKxWUy+jlhk/JN4PtETRHIIyxAwCs9H/rDUPJ6z1tDtoiEN44PgZY5kCZmrnPXHudvV6k1aUAsBq/yQHLfLJX3ZCmp8RwDhcA8hiwyPdF6sUt5h3ptEEAIuQxYIlij9ywyJjYQ29Xj84xABHyGLBE1le+61K1Yc1pgAB+wu4AqG0f7DYdujybyU0jAPyM86uBWrW9SD24xpw9gDAG8Av0j4Has/aIummJ8XB7vWN9ho0B/AL9Y6CWVBhy+9fGEx31G1vwOxjAichjoJb0nu3LTNZ+QxgDOBnyGKgNk7aZ2aXcMQLAKfFTHahxiw6qV7aaC6+2M2gM4FT4tQ7UrBKv3LzU935fexvubQzg1OgfAzXIY8qT3xndk/U+jQljAKdD/xioQS9sMr8vUp9k0dAAnAG7CaCmrCtQk7YZm4Y54h1WlwIg4HG8GqgRM/eaV831vd7T1jja6lIABAP6x0CNeOhbc2aW/fJGDBsDqBL6x0A1K/LIHV8bdWxySQPCGEBVkcdANbtxsa9upKwbYo+geQGoMo5XA9VsXYH6pL+9Dm0LwNngBzxQnfaVqnqRWiwnVAM4S+QxUG3cPnlvt2pXj2FjAGeNY2pAtbl+kU8TeaWHzepCAAQf8hioHvluWV+gdo90RNOqAJw9jlcD1SCnTA2c6/ttS50wBnBuyGPgvHhNefRbo+vnvjsz9H9250g1gHNU1R/z+/fvf//99zVNGzVqVLNmzU7415ycnLlz5+bm5qalpY0cOTIiIkJEsrOz582bV7nMb3/72+hoZg5EqJmVba45otYNtTeJ5jQuAOeuSv3j3NzczMxMp9NZUlKSmZmZl5d3wgKDBg367rvvYmNjp02bdvnll3u9XhHZvHnzuHHj9vyPz+er/vIBq03ebt7aSieMAZynKvWPX3/99WuvvfbZZ58VkYMHD77xxhtPP/308QusWbMmMjJSRO67774mTZqsX7++e/fuItKqVatx48ZVf9VAAMh3y81LfcVeuaUV4z4AzleV9iMrV67s16+f/3G/fv1Wrlx5wgL+MBYRj8fjcrnq1q3rf5qdnf3oo4+OHz/+4MGD1VQwECiuX+i7qqm+9jq7jb4xgPNWpf7xoUOHkpKS/I+TkpJ+fby60j333DNs2LD09HQRSU5OHjZsWPPmzdevX9+mTZtvvvkmIyPjpO/yeDz5+fmVkW+32x9//PGuXbue3Z9y9srKymr6KwIfK6GiokLXdYfj7KbU+vaoftRtH31huSsk1h+bgVLK7XYrpawuxEoul8s0TU0L6x+YNdQWdF0/4xlUVcrj6Ojo8vJy/2O32x0bG3vSxR577LHdu3dXnsN1ySWXXHLJJf7HSqlJkya9+uqrJ31jREREYmLiE088UVl3t27dYmJiqlLbeTrV3xJWwnwlOByOc8jjPy/2PdtVj4uNqKGqal+YbwZKKV3Xa2e3E7A0TYuOjg7zPBbr2kKV8jg1NXX37t3+x3v27ElNTf31Mk8//fTSpUvnz59/0g06LS1tw4YNp/mKqKiorKysqhQDWEuJ/Hmt4TFl+IUMGwOoNlXaoYwaNertt98uKytzOp1vv/32qFGjROTIkSNPPfWUf4G///3vH3300YwZMwzDKCws9Hg8IrJhwwa32y0i2dnZU6dO7dOnT039EUAtemiNseqwWnYNE38AqE5VyuOhQ4f27ds3NTW1efPmV1555eDBg0WksLDwrbfe8i8wc+bMioqKrKyszMzMzMzMRYsWicinn37aqFGjBg0aZGZmDh8+fPTo0TX3ZwC1495VxrI89WmWnTs4AaheWtXPX/BfVXy2w2xut7tOnTqnXyY/P79nz56Vh8RrjdPpDPMxM2ElnM35XDll6tJZxq6R9oiQO1DNZqCUcrlcYT5+XFZWxvixhW3hLI65nW0S+50xjIFg8advzKwULfTCGEAgYAwMqJJDbll1WGXfSJMBUCP4qQ9UyeZjqkM9YeoPADWEPAaq5MdilZFAGgOoKeQxcGafZ5vPbTKHNKe9AKgpDIYBZ5DvlrErjNkD7J2T6B8DqCn83gfOYFGumZWiE8YAahR5DJzBGzvM37SgpQCoWexlgNOZm6Py3TKgKZ1jADWLPAZOqdyQ+1cbr/e0cZkTgJpGHgOn9PeNRmaS1r8JaQygxnF+NXBy83PVv3eYW4dz4wgAtYE8Bk60q0Q9s8FcmKvmXmVPirK6GgDhgTwGTvTUOrNxtGwaRhgDqD3kMfALH+w2dxSpN3rZ4zhQDaAWcT4X8AufZau/dNIJYwC1jDwGfpbvlsUHzT6NaRcAahv7HeBnwxb4Hm5vqxdpdR0Awg95DIiIHHSpO742nD55uD2NAoAF2PUAsiRPLv/SuDBOW3mtnbk/AFiC86sR7oo88tuvZerltkEXkMUALEP/GGFNiTy4VrulpRDGAKxF/xjhy1Ry1VxfhU8mXWJ1KQDCHv1jhKljFdLpM19ihDbvShXL1cYArEYeIxwty1NtZ3qvv1D/qD+3UgQQEDhejbBzoEzdu8qYfKltaHN+jwIIFOyPEF6UyG3LjNvSdcIYQEBhl4Tw8sIm01TyhzZs+QACC8erEUb+vtF8bbu5cZjdQRwDCDDkMcLFIbe8us3Ydr2D6akBBCC6CQgXqw+bnZI0whhAYCKPERZcPnn0W/OhdjarCwGAkyOPERauW+Dr0UDr24SLjQEEKMaPEfq2Fao9JTJnAFN/AAhc9I8R4g67ZcQi449tdTsbO4AAxi4Koey9XWbf2b5RLfT7uOAYQGDjeDVC1vRd5pPrzHd723o14kA1gEBHHiM0Tdpmvv69OfcqW0YCYQwgCJDHCEFL8tSLW8w119kb1rG6FACoGvIYoebZDebr35vT+tgIYwBBhDxGSPlsn/nhHnPb9fbECKtLAYCzQR4jROS55On1xqxsc/5AwhhA8OEiEISCZXmq82fe1Fht83BHu3qcwAUg+NA/RtA7WiG/X2H8q5ftulR+XwIIVuy/ENzKDek72zeqpU4YAwhq7MIQ3BbkmnUj5KlObMkAghvHqxHEfKY8+Z35cg/uoggg6NGrQLAylIxYZFwQq13OdJgAgh95jGD1wDeG11Qf9rNxG0UAIYDj1Qg+hpK/bzTn5qj1Q+x12IQBhAT6xwgyeS7p/Jlv8zH1zWB7rMPqagCgmtC5QJCZvN3o20R75RLO4QIQUugfI5jsLFbTflS/v4jtFkCoYb+GoOHyyRVzjHHd9NaJnMEFINSQxwgab+80+zTWRrVgowUQghg/RnD4eK/53Ebzv1cybAwgNJHHCALrC9R9q4x5A+0duHcTgBBFHiOguXzy/CbjPzvVuK42whhACCOPEaDKDZm0zZy60+yWrK0bYm9Qx+qCAKAmkccIUA+tMfY7ZVpvW9dkusUAQh95jEA0N0d9uMf8/npHUpTVpQBAreDSEQSc93ebv1vm++8VdsIYQPigf4zAsvqwemiN8fUge6sEDlMDCCPkMQKC15TnN5kz95oHytSMvoQxgLBDHsNKBeWyrkCtPaKm7TK7Jmkf9rOlJ2g6WQwg/JDHsMy+UtXzCyMzWetQTz7uz+XFAMIaeQwL/FCsblxslHjkj231h9tzUiEAkMeoLSVeOeBUKw6p17abB8rU811t3DYRACqRx6hxT6w13txhOnRpFqs1i9Fe72m7pAGDxADwC+QxatamY+qzbLVjBDN7AMDpkMeoKYaSlYfUcxuNOzN0whgATo88RvUwlRR5pLBCFXmksEIKytXbO81Cj1zdTGOcGADOiDxGNRi52JhzwKwbqSVGSN1ISYzQEiPkxhb6zS11O1kMAFVAHuMcOb3yY4lae0SN32LqIntGMkIMAOeOPMZZc/nk5qXGsjwzPUHLSNQ+6GvrnMTZ0gBwXshjVJXXlHGbzPd3m4fc6ppm+uHfOrhmCQCqC3mMM9taqCZtM+fnqi5J2udX2NITNIIYAKoXeYzTUSJ/XG3M2q/+1FZ/oqOeGksQA0CNII9xciVeWXLQfHunKjfU99fbo2xWFwQAIY08xs9cPlmSpxblmisOqb2lqldD/Zpm2m3pNgfXLAFADSOPw9rOYrUyR8/zmntK1Yajak+J6tlQu7KpfnMrvX09zcbBaQCoLeRx+FpXoK5bYPRpoGfUk35NtD+21VsnksEAYA3yOLyUeiWnTB0okw0FatqP5vhu+qCGvthYJvIAAIuRx6Evt0xtKZTnNhobjqpYuzSO1i6I1TrUk1d62K5I0ZxOq+sDAJDHIW9pnrphka9nQ31Yc33R1TpnZgFAYCKPg1VhhRR6lIiUeMRQUuYTjyHlhrgN5TOl1CtK5FiFTN5uTrrUNjKNHAaAgEYeB5Mj5bLxqPr2iPpyv3mgTKJskhghIpIQIbomsXbNoUsdu0TZJEKXGIfoIm9fbrsihXO0ACDQkceBq9gjh93qoEsOlKkfi9U7Pyqfko71JDNZezbT1q8Jk1YCQOggj63h9snyQyrPpQrKpaBcHXbL0Qop9qhijxR5pMijPIbER0j9SO2CWGkSrV0Qqy0YqLdKIIIBIDRVNY8LCgqWL1/eoEGDSy+9VDtZx2z//v1r165t2bJlhw4dKl8sLi5etmxZXFxc7969dZ0hTDlSLi9uMb49rDYfU12TtbR4LSlSmkRr7epJUpRWL1ISIyQhQkuIkAjWFgCEkyrl8caNGwcMGDBw4MDNmze3a9funXfeOWGBL7744q677ho0aNDSpUtvvfXWJ598UkT27t3bq1evvn377t27Nz4+fvbs2WESyT5TNh9Tu0tViUdKveL0itOnCiuk1CvfFahBzbRnMm2d6mvRHJsAAPyPppQ640LDhg3r0aPHww8/XFZW1qpVq6+++qpjx47HL9CuXbtnnnlmyJAhOTk5rVu33rdvX/369ceOHRsdHf3iiy96vd4OHTq89NJLV1111Uk/Pz8/v2fPnrt3766ev6nKnE5nbGzsaRYo9YrLJ2U+VVghLp+4fFLiVT+96JUijyrzidMrTq8UepTzp+iVY+WqdV3t4kQt1iGx9p/6u7EOiXNIUpTWs2FgHXM+40oIeRUVFbquOxwOqwuxEpuBUsrlcsXExFhdiJXKysqio6NPegQ0fFjYFqrUR5s7d+4LL7wgIjExMVlZWXPmzDk+jw8cOLBz586rr75aRJo2bdq2bdslS5Zcf/31c+bMmT59uog4HI5rrrlm7ty5p8rj81FYIaVe5fRJqUdKvVLkUUrE6RWvKSJSYYjLJyJiKCnx/vTLo9gjphIRKS53mJohIh5Tynw/XSZU7JEyn3L5xGNKtE0SIrQ4h0TbJcYhiRFatF2i7ZIQIbF2LTlKu9AhdSMl1q7FOiTWIQkREu/Q4jnaDAA4S2fO48LCQrfb3ahRI//TJk2aHDx48PgF8vLy6tevHxERcfwCSqn8/Pzj37Vq1apTfYXX6z127NjDDz9c+crNN9+ckZFxmqq+PCC3rdRtmsTYJdYucQ5JiFBxDi3Wrvx3BnToEusQEdFEEhwiInG6xDnErovUkTp2idRVeXl5bJ2oOMdPPwYTI1SMXaLt/k6t6Of6G1F5peIc32qBioqKMO8a+vvHpmlaXYiV2AyUUhUVFXZ7WA8jVVRU2Gy2MO8f11Bb0DStMiVP5cwbn81mE5HKvZVhGCdssifsywzD8P8ftdlsp3nXCV9ht9vT0tIqX4mPjz/9YPPgVDmWesJr2nH/rQrN6TRjY4//lnDcCnVdD5Nx/VPR/8fqQqzEGlBKsRL8ayDM87iGNoOqrNUz53F8fHxsbGxubm5CQoKI5OTkdO7c+fgFmjRpcuzYMbfbXadOHf8CKSkpItK4cePc3Nz09PTjXzwpXdfj4+PHjh17xmKql8PhCPM+gbASREzTZPyYzUApxUrwr4Ewz2MLN4Mq/QoYMmTIjBkzROTo0aMLFy687rrrRGTv3r3Lli0TkSZNmnTq1Onjjz8Wke+//37Xrl39+/f3v+v9998XkbKyslmzZvnfBQAAfq1KgyV//etfs7Kyvvnmm507d95xxx3+Lu+8efOmT5++YsUKERk/fvwNN9wwY8aMjRs3vvDCC3FxcSLy0EMP9e3bt0+fPgcPHuzdu/dll11Wo38JAADBq0rXO4lIRUXFpk2bGjZsmJr607BteXl5eXl5YmKi/2lpaen27dtTU1Mrz+ESEZ/Pt3HjxsTExJYtW57mwwP2eqdwwErgeidhM+B6JxHheicRCfzrnUQkMjKyW7dux78SFRUVFfXzfezj4uK6d+9+4qfb7ZmZmedZIgAAIS+sTyb86KOPCgsLra7CSi6Xy3+NeDhbvXr1pk2brK7CYm+99VaYX/GVk5PzxRdfWF2Fxb744oucnByrq7CSaZpvvfWWVd8e1nk8derU2j9IHlAOHjw4efJkq6uw2Ny5c/1nJoazcePGOZ1Oq6uw0ubNmz/88EOrq7DYhx9+uHnzZqursJLT6Rw3bpxV3x7WeQwAQIAgjwEAsF6gTA6nlKr9oVyfz1dSUhLOQ8jFxcWGYYTzGhCR8vJyt9sd5ivBNM2ioiLDMKwuxDJOp9Pr9Yb5ZuD1ep1OZzivhNLSUtM0a2INREZGRkdHn36Zql7vVKOOHDmSmpp6/NnatcPlckVFRYXzDHmmabrd7jC/xqOioqIqU8uGNq538vl8Xq/XP8lg2HK73Q6HI8wn8a6hthAbG7tlyxb/NJenEhB5DABAmAvfriEAAIGDPAYAwHrkMQAA1iOPAQCwHnkMAID1wvq8dj+Px7No0aIffvhB07SePXuG5w0w3G73woULd+7cWb9+/aFDh57+pPxQ5fV6t23btmvXri5dulx44YVWl1N7lixZ8u2336anpw8ZMiQ87+2Tl5e3ceNGm8125ZVXWl2LNUpKShYsWLBnz56mTZsOGTIkDK/7Mk1z6dKlW7duNQyjc+fOvXv3rv0a6B/LgQMH3nzzzaKiosOHD1977bWTJk2yuiILjB49+s0333Q6nXPnzr344ovz8vKsrsgCw4YNu+mmm8aOHbtkyRKra6k9EydO/P3vf69p2rhx4+6//36ry7HAggUL2rdv/9hjjz344INW12KZ4cOHf/jhh263+913383MzAzDyczdbvf48eMLCgqKiopuv/32Rx991IIiFI7zwQcfdOnSxeoqLFBWVlb5uH///hMmTLCwGKt4vV6l1BVXXDFlyhSra6klHo+ncePG3377rVLq8OHDMTEx+fn5VhdV2/z/3+fMmdO2bVura7FM5R7AMIyLLrpo5syZ1tZjrRUrVjRs2LD2v5f+8c+8Xu/y5cs7duxodSEWOH4iN5fLFRcXZ2ExVgnDaYl27drldrv9YzTJyclt2rT55ptvrC6qtoXh//dfO34P4Ha74+PjLSzGWqZpLlmypFOnTrX/1eGyIZaUlLhcrhNejI2NrZwXrVGjRkePHs3IyFi6dGltF1dbCgoKfD7fCS/Wq1fv+Kki//Of/xQUFNx44421W1otqaio+PXMtA6Ho379+pbUY7n8/PykpKTKMePk5OTwHKpApeeeey4lJaVfv35WF2KNTp06bd++vVGjRosXL679bw+XPP7nP//565ubjh079oEHHvA/zs/PdzqdTz311IgRI0J1+HDkyJEHDhw44cXp06d369bN/3j27NlPPvnkggULQvVsjtWrV48ePfqEFy+66KJZs2ZZUo/l6tSpU1FRUfm0vLz8jFPeI4RNmTJl2rRpy5YtC9tZ/Tds2FBeXj5hwoRBgwZt3brVZrPV6tfX/iHyQJadna3run88KdzMmzcvJSVl06ZNVhdisbAaP87JyYmMjKwcO0xLS1uyZImlFVkmzMePlVLTpk1r0aJFdna21YVYz+Px6Lq+f//+Wv7eMP0RdLx9+/YVFRWJiGma7777btu2bcNwPGn58uW33377559/3r59e6trQe1JSUnp1q3b1KlTRWT+/Pkej6dXr15WFwULfPLJJ0899dT8+fMvuOACq2uxRl5e3qFDh/yPp0+fnpSU1KRJk9ouopbzPwD997//bdSoUbNmzZKSknr16rVlyxarK7LADTfcUPc4f/7zn62uyAKPP/748Svhu+++s7qi2rBp06b09PSWLVumpKQsWLDA6nIskJOTc/z/99tvv93qiizQs2fP41fCa6+9ZnVFtW316tUpKSkpKSnJycmdO3desWJF7dfA/RZ/UlxcHB0d7XA4rC4EqG1KqaKiovj4+NoeLQMCTGlpaURERGRkpCXfTh4DAGA9xo8BALAeeQwAgPXIYwAArEceAwBgPfIYAADrkccAAFiPPAYAwHrkMQAA1iOPAQCwHnkMAID1yGMAAKz3/wKSckdIRqPeAAAAAElFTkSuQmCC\" />"
            ]
          }
        }
      ],
      "source": [
        "# your code to adapt\n",
        "N = 1000\n",
        "x = sort(randn(N))  # this draws normals, adapt to our areto\n",
        "F_hat(x) = (1:length(x)) ./ length(x) # the cdf values at the x values are just the counts normalized by the total number of observations  No need to modify.\n",
        "plot(x, F_hat(x), label=\"Empirical CDF\") # adapt to do the log-log plot "
      ],
      "id": "74792fdb"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This shows methods for understanding the tail behavior of distributions.\n",
        "If it approaches a straight line, then it is called a power-law tail.\n",
        "\n",
        "### Part (d)\n",
        "\n",
        "Now lets do the same empirical check on the tail behavior with a\n",
        "LogNormal distribution - which does not have a power-law tail and the\n",
        "[Frechet](https://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution)\n",
        "distribution - which does.\n",
        "\n",
        "First, adapt your code above to (1) draw `N=1000` elements from the\n",
        "`LogNormal(0.0, 0.5)` (2) plot the pdf to get a feel for the\n",
        "distribution on a reasonable range, then finally (3) plot the\n",
        "theoretical and empirical log-log plots as you did before."
      ],
      "id": "d520255b-0661-4ca2-8342-9736468fa351"
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "LogNormal{Float64}(μ=0.0, σ=0.5)"
            ]
          }
        }
      ],
      "source": [
        "# your code to adapt\n",
        "N = 1000\n",
        "d = LogNormal(0.0, 0.5)\n",
        "# add code here to plot the pdf\n",
        "# add code for the log-log plot"
      ],
      "id": "539f0eb3"
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [],
      "source": [
        "# your code to adapt\n",
        "# add code for the log-log plot"
      ],
      "id": "b14af449"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next do the same thing except with a `Frechet(1.5, 1)` (which has a tail\n",
        "parameter of 1.5)."
      ],
      "id": "9a10db8c-a4a9-44d3-aacf-9b7ad5db649b"
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {},
      "outputs": [
        {
          "output_type": "display_data",
          "metadata": {},
          "data": {
            "text/plain": [
              "Frechet{Float64}(α=1.5, θ=1.0)"
            ]
          }
        }
      ],
      "source": [
        "# your code to adapt\n",
        "N = 1000\n",
        "d = Frechet(1.5, 1.0)\n",
        "# add code here to plot the pdf\n",
        "# add code for the log-log plot"
      ],
      "id": "33df2bfe"
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [],
      "source": [
        "# your code to adapt\n",
        "# add code for the log-log plot"
      ],
      "id": "b3928c8d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part (e)\n",
        "\n",
        "Comparing your answers in Part (d) and (c) what do you notice is\n",
        "different theoretically between the log-log plots of distributions with\n",
        "and without power-law tails? Play around with it and see if you feel you\n",
        "could just look at an empirical plot and tell if it is a power-law tail\n",
        "or not?"
      ],
      "id": "8a5368b8-6cc8-4e3c-bacc-976a506cd863"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "(double click to edit your answer)"
      ],
      "id": "7cc0b430-7e7e-4f27-b5ba-52f80b30f0a2"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "name": "julia-1.11",
      "display_name": "Julia 1.11.2",
      "language": "julia",
      "path": "C:"
    },
    "language_info": {
      "name": "julia",
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "version": "1.11.2"
    }
  }
}