economics_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\\rightarrow$Run All).\n",
    "\n",
    "Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\", as well as your name and collaborators below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"\"\n",
    "COLLABORATORS = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Introduction\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the notebook to accompany the course Applied Economic Analysis at Tilburg University. The idea is to bring economic concepts \"alive\" by programming them in python. The choice of topics is loosely based on Tirole (2017).\n",
    "\n",
    "The point is not that we go into models in detail. Instead, we sketch the trade offs and then model these in python. We make graphs to explain the intuition. In this way, you learn both python and to use economics in a more \"broad brush\" fashion than \"taking the third derivative\".\n",
    "\n",
    "When we go through the notes below, there will be **Questions** and **Assignments**. Questions we will typically do in class together. Assignments will take more time and you do these at home. If you have questions about these, come back to it in class.\n",
    "\n",
    "**Note** that you only learn to program by actually typing code and doing things \"wrong\" (e.g. you see error messages). Just typing in correct answers may seem useful but it actually is not. Also play around with the code below; change things to see what the effects are; experiment with different solutions than what we did in class etc. \n",
    "\n",
    "In contrast to, say math, programming you can learn by yourself because you can see whether what you tried worked or not. If you believe that the derivative of $x^2$ equals $5x$, no one will correct you. If you want to plot a figure and you get errors, you did something wrong and should keep on trying...\n",
    "\n",
    "## Why this course?\n",
    "\n",
    "In short, what are we trying to achieve with this course:\n",
    "\n",
    "* you learn to talk about economics \"outside your comfort zone\":\n",
    "    * that is, not within the confines of a model\n",
    "    * but also not \"chatting away\" as if you never saw a model before in your life\n",
    "* you learn to program in python\n",
    "    * you only learn python by actually doing it\n",
    "    * not by waiting to type in the perfect answers\n",
    "    * programming is frustrating, but we are around to help you\n",
    "* goal is not to rush through this notebook to make sure we finish everything in time\n",
    "* if you find the pace in class too slow, feel free to work ahead; it is set upin such a way that you can go through things on your own"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Python packages that we will be using\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need the following packages to run the code below. If you get an error saying that the package is not available, you can install it using either `pip install` or `conda install`.\n",
    "\n",
    "If you are using this on the russet server, you can install a package that you need by doing the following.\n",
    "\n",
    "Say, you want to install [wordcloud](https://amueller.github.io/word_cloud/) on the russet server. Then create a new cell in a python notebook and type the following:\n",
    "\n",
    "```\n",
    "%%bash\n",
    "\n",
    "pip install --user wordcloud\n",
    "\n",
    "```\n",
    "\n",
    "It will install the package in your environment on the russet server.\n",
    "\n",
    "E.g. if the following cell gives an error that the package \"wbdata\" is not installed, run the cell below the import statements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy import stats, optimize\n",
    "import random\n",
    "import wbdata as wb\n",
    "\n",
    "plt.style.use('seaborn')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: wbdata in /home/janboone/.local/lib/python3.8/site-packages (0.3.0)\n",
      "Requirement already satisfied: appdirs<2.0,>=1.4 in /home/janboone/anaconda3/lib/python3.8/site-packages (from wbdata) (1.4.4)\n",
      "Requirement already satisfied: requests>=2.0 in /home/janboone/anaconda3/lib/python3.8/site-packages (from wbdata) (2.25.1)\n",
      "Requirement already satisfied: tabulate>=0.8.5 in /home/janboone/anaconda3/lib/python3.8/site-packages (from wbdata) (0.8.9)\n",
      "Requirement already satisfied: decorator>=4.0 in /home/janboone/anaconda3/lib/python3.8/site-packages (from wbdata) (5.0.9)\n",
      "Requirement already satisfied: urllib3<1.27,>=1.21.1 in /home/janboone/anaconda3/lib/python3.8/site-packages (from requests>=2.0->wbdata) (1.26.4)\n",
      "Requirement already satisfied: idna<3,>=2.5 in /home/janboone/anaconda3/lib/python3.8/site-packages (from requests>=2.0->wbdata) (2.10)\n",
      "Requirement already satisfied: chardet<5,>=3.0.2 in /home/janboone/anaconda3/lib/python3.8/site-packages (from requests>=2.0->wbdata) (4.0.0)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in /home/janboone/anaconda3/lib/python3.8/site-packages (from requests>=2.0->wbdata) (2021.5.30)\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "pip install --user wbdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# The market\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is based on Tirole chapters: 1, 2, 6.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## why do we love the market?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many (but not all) economists love the market as an organizing institution. This affection for markets is traced back by some to Adam Smith and his [\"invisible hand\"](https://en.wikipedia.org/wiki/Invisible_hand). Others may not think that the market is so great, but basically distrust other organizing institutions like a government. They will point, for instance, to the fall of communism in the [former Soviet Union](https://en.wikipedia.org/wiki/Revolutions_of_1989).\n",
    "\n",
    "Although you, as an economist, have seen many market models already,\n",
    "it may be illustrative to go back to basics. To understand the\n",
    "advantages of the market, let us consider a very simple economy. We\n",
    "focus on one type of product and there is an exogenous endowment of\n",
    "this product (supply) equal to `number_of_goods`. Further, in this\n",
    "economy there are `number_of_agents` agents who have a valuation for\n",
    "this good which is randomly drawn from a normal distribution. An agent\n",
    "can at max. consume one unit of the product and her utility is then\n",
    "given by her valuation. The value of consuming zero units and the\n",
    "additional value of consuming more than one unit both equal 0.\n",
    "\n",
    "The vector `valuations` contains for each agent her valuation and we sort this vector such that the first agent has the highest valuation.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_agents = 1000\n",
    "number_of_goods = 100\n",
    "\n",
    "valuations = np.array(sorted(pm.Normal.dist(100,20).random(size=number_of_agents),reverse = True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we are using [pymc3](http://docs.pymc.io/notebooks/getting_started) here to generate random numbers from a distribution. There are also other python libraries which can do this, e.g [scipy](https://scipy.org/). We use pymc3 as we will use it later for estimation as well.\n",
    "\n",
    "\n",
    "**Question** What is the economic name for the following expression? To answer this question, you need to understand both how indexing works in python and which economic concept is captured by this expression.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"{0:.2f}\".format(valuations[number_of_goods]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1022894352f41ade043f5914b4f8fbed",
     "grade": true,
     "grade_id": "cell-992dd3c425d326e4",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that we have an omniscient social planner who knows the valuation of every agent in the economy. This planner aims to maximize the (unweighted) sum of agents' utilities.\n",
    "\n",
    "**Question** Calculate the total welfare that this planner can achieve. Denote this value `max_welfare`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "12489fe2b128707cb69c36ef8207e665",
     "grade": true,
     "grade_id": "cell-dc7591c2493f44d5",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, this is the best that we can do. It gives us an upperbound on the welfare that can be achieved. Having an omniscient social planner seems unrealistic, but perhaps there is an institution that can achieve this outcome without omniscient intervention. You guessed it&#x2026;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### market outcome\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compare the maximum welfare that a planner can achieve with the market outcome.\n",
    "\n",
    "****Question**** Define a function `demand(p,valuations)` which has as arguments a price $p$ and a vector of agents' valuations. This function returns the number of agents who are willing to buy the product at price $p$. Since each agent who buys, buys exactly one unit in our set up, this function returns demand at each price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "059102389055e61614dd2bcfe23eb3c6",
     "grade": true,
     "grade_id": "cell-6e537cfe5c558002",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Using `matplotlib` plot (fixed) supply and the function `demand` against price, where we maintain the economic convention of having quantity on the horizontal axis and price on the vertical axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "ef8fcc8d495bb64988ff9b0a8462fae5",
     "grade": true,
     "grade_id": "cell-edc08a2fd2d198d0",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to calculate the equilibrium price, we define a function `excess_demand`. We will then look for the price where `excess_demand` equals 0; this is the equilibrium price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def excess_demand(p,valuations,number_of_goods):\n",
    "    return demand(p,valuations)-number_of_goods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to find the equilibrium price, we use from `scipy.optimize` the function [fsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html). If you want to know more about solve, there are people at [stackoverflow](https://stackoverflow.com/questions/8739227/how-to-solve-a-pair-of-nonlinear-equations-using-python) discussing this function. Stackoverflow is generally a great resource if you are wondering how to solve a problem in python.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "price = optimize.fsolve(lambda x: excess_demand(x,valuations,number_of_goods),120)\n",
    "print(\"{0: .2f}\".format(price[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, now we know the equilibrium price\n",
    "\n",
    "****Question**** Calculate total welfare at this equilibrium price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8d746f8923b01a8dd3c65615bf9eff87",
     "grade": true,
     "grade_id": "cell-ce3db0f6b82b1bd5",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** How does this welfare compare to the maximum welfare that the omniscient social planner can achieve? Recall that this level is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"{0:.2f}\".format(max_welfare)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "4157b950ffa265433f2f439978e2fe39",
     "grade": true,
     "grade_id": "cell-12281671fab9433b",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### elastic demand and supply\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up till now we assumed that supply was inelastic: there was a given endowment of goods and this was auctioned off to consumers. Now we assume that some agents initially own the goods. However, these agents are not necessarily the ones that value the goods the most.\n",
    "\n",
    "In particular, we give `number_of_goods` agents one unit of the good. They become suppliers.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.shuffle(valuations)\n",
    "valuations_supply = valuations[:number_of_goods]\n",
    "valuations_demand = valuations[number_of_goods:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Define a function `supply` which depends on the price and the valuations of the suppliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8e49fe86042374bc904781fabe82be57",
     "grade": true,
     "grade_id": "cell-16c6b97a04b8f81b",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Use `matplotlib` to plot demand and supply in a single figure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "7152e41f2125a49371acbd179f023de5",
     "grade": true,
     "grade_id": "cell-2ef0599165ae9dbc",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Define the function `demand_minus_supply` which looks\n",
    "like `excess_demand` above but now with elastic supply. The function\n",
    "depends on the price, the valuations of people demanding the good and\n",
    "the valuations of people supplying it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1ee31021850576d2da61e42b1bc517f8",
     "grade": true,
     "grade_id": "cell-65682795af4c438c",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Use `fsolve` to find the equilibrium price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "db4bb99e584d0deb9c0ad7247db153d0",
     "grade": true,
     "grade_id": "cell-c2036a71256a3f4f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** How does the equilibrium price here compare to the equilibrium price above with exogenous supply? Is the price here higher? Why (not)? Is welfare higher here than above?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "012c6347e584f4cfa008d0ecd4379bbc",
     "grade": true,
     "grade_id": "cell-cc4eaf62b184e229",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "e7d953fe1e3f83c7051419611398c637",
     "grade": true,
     "grade_id": "cell-6aa67ec501987431",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## why do others not love the market?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the results above look great, the assumptions we made, may not be realistic in every market. To illustrate, without saying so, we assumed above that the market is perfectly competitive and without external effects. Here we program three reasons why the market outcome may not necessarily lead to maximum welfare. First, we look at income inequality and the problem that this causes for the market. Then we consider market power and finally we model external effects.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### income distribution\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In micro economics we usually do not do much with income distributions. Often because models where income distributions play a role are tricky to solve analytically. But here we program/simulate and hence we do not worry about analytical solutions.\n",
    "\n",
    "Now in addition to the valuations introduced above (the utility an agent gets from consuming the good), we need an income distribution. The former determines the willingness to pay (wtp) for an agent, the latter the price an agent can pay. A consumer is willing to buy the product at a price $p$ if both her wtp and her income exceed $p$.\n",
    "\n",
    "First, we randomly draw an income for each agent in the economy.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "incomes = pm.Normal.dist(100,20).random(size=number_of_agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to redefine demand, now denoted `demand_2` which takes into account both whether an agent values the good more than $p$ and whether she can afford $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def afford(p,incomes):\n",
    "    return incomes>p\n",
    "\n",
    "def wtp(p,valuations):\n",
    "    return valuations>p\n",
    "\n",
    "def demand_2(p,valuations,incomes):\n",
    "    return np.sum(afford(p,incomes)*wtp(p,valuations))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Define the function `excess_demand_2` which depends on $p$, agents' valuations, incomes and number of goods (which we assume to be inelastically supplied again)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "6691352c1cdbd89bf7a953b3ea6786ec",
     "grade": true,
     "grade_id": "cell-5776d2bf65e4b258",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Use `fsolve` to determine the equilibrium price in this case. Is this price higher or lower than above? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "169a39bc6ce1158a26c7d091006a9567",
     "grade": true,
     "grade_id": "cell-bf49952904f411cb",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5e4425a5e12b951aede635ef596a969d",
     "grade": true,
     "grade_id": "cell-0d5fd89ae902ceb1",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Calculate welfare in the market equilibrium. How does it compare to `max_welfare`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "61923cf0d787f5bc0047db043a439746",
     "grade": true,
     "grade_id": "cell-a22e02589e9baa2f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5fdd9ad919ee313131e2279541a37063",
     "grade": true,
     "grade_id": "cell-bf1e4efc594edac1",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1aa18650c678dbe953c43cbb439e18bd",
     "grade": true,
     "grade_id": "cell-e5028eb1b241aaee",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Assignment**** Model an economy where an increase in income inequality reduces welfare.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "7408bcc3935bdee87e04bbbe2344973c",
     "grade": true,
     "grade_id": "cell-fa2f71f71129011b",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "3541d96573b52569b6b7cdfae16f1303",
     "grade": true,
     "grade_id": "cell-a17482be644a4ce0",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### market power\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Warning**** We are going to do a couple of things wrong in this section. No need to panic; this actually happens a lot when you are programming. Use your economic intuition to see where the mistakes are and correct them.\n",
    "\n",
    "Suppose that we now give all the products to 1 agent who then owns `number_of_goods` units of this good. To simplify, we assume that this agent values the good at 0.\n",
    "\n",
    "****Question**** Suppose we use the function `demand_minus_supply` defined above to calculate the equilibrium price. Would the equilibrium price increase due to market power? Why (not)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "2dd3769b69b6237dcd894ba86d6f0a66",
     "grade": true,
     "grade_id": "cell-ae595d48aa3c2109",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps a monopolist would not use an auction to sell all the goods. Let's calculate the profits of the monopolist as a function of the price and the valuations of the agents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def profit(p,valuations):\n",
    "    return p*demand(p,valuations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "range_p = np.arange(0,140)\n",
    "\n",
    "plt.plot(range_p, [profit(p,valuations) for p in range_p], label = \"profit\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"$P$\")\n",
    "plt.ylabel(\"$\\pi$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the profit maximizing price is around 80. Recall the equilibrium price under perfect competition above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "price"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Since when does a monopolist charge a lower price than a perfectly competitive market?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "397246e4a54932e39369acda57069376",
     "grade": true,
     "grade_id": "cell-e7432304771a32b7",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "aed50c825612231b67f1b8886c67a7d2",
     "grade": true,
     "grade_id": "cell-ce5155f9fa17d736",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Assignment**** Calculate the profit maximizing price in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "6b03b83814f5ff741888a44eb23c45f0",
     "grade": true,
     "grade_id": "cell-b024a8a8ccd90426",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "9b57ee58e9f9fe67f3a3c62a2ce641c0",
     "grade": true,
     "grade_id": "cell-834900a6e8f5e531",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### merger simulation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is based on Tirole chapter 13.\n",
    "\n",
    "In this section, we model a more standard oligopoly market with\n",
    "Cournot competition. We start with three firms and then calculate what\n",
    "happens if two firms merge such that only two firms are left in the\n",
    "industry. Hence, we first calculate the equilibrium with three firms,\n",
    "denoted by 1, 2 and 3. Then firms 2 and 3 merge so that we are left with 2 firms; denoted by\n",
    "1 and 2.\n",
    "\n",
    "We are interested in the effects of the merger on the equilibrium price.\n",
    "\n",
    "We assume that before the merger each firm has constant marginal costs\n",
    "equal to 0.3. We assume a simple linear (inverse) demand curve of the\n",
    "form $p=1-Q$ where $p$ denotes price and $Q$ total output on the market.\n",
    "Total output equals the sum of each firm's output: $Q= q_1 + q_2+q_3$.\n",
    "\n",
    "The function `reaction` gives the optimal reaction of a firm to the total output `Q_other` from its competitors. In this function, we use the routine [fminbound](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html). Python does not have maximization routines, hence we minimize \"minus profits\" (which is the same from a mathematical point of view). The parameters `0,1` in this routine give the bounds over which we optimize. Since demand is of the form $p(Q)=1-Q$, we know that no firm will choose $q>1$; further we also know that $q \\geq 0$.\n",
    "\n",
    "The fixed point makes sure that for each of the three firms, their output level is equal to its optimal reaction to the output levels of its competitors. If each firm plays its optimal response, given the actions of the other players, we have a Nash equilibrium.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c0 = 0.3\n",
    "vector_c = [c0]*3\n",
    "\n",
    "def p(Q):\n",
    "    return 1 - Q\n",
    "\n",
    "def costs(q,c):\n",
    "    return c*q\n",
    "\n",
    "def profits(q,Q_other,c):\n",
    "    return p(q+Q_other)*q-costs(q,c)\n",
    "\n",
    "def reaction(Q_other,c):\n",
    "    q1 =  optimize.fminbound(lambda x: -profits(x,Q_other,c),0,1,full_output=1)\n",
    "    return q1[0]\n",
    "\n",
    "def fixed_point_three_firms(vector_q,vector_c):\n",
    "    return [vector_q[0]-reaction(vector_q[1]+vector_q[2],vector_c[0]),\n",
    "            vector_q[1]-reaction(vector_q[0]+vector_q[2],vector_c[1]),\n",
    "            vector_q[2]-reaction(vector_q[0]+vector_q[1],vector_c[2])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the reaction functions of a duopoly where both firms have zero marginal costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8c1d2e2a5ead0c1f987ffc405736e315",
     "grade": true,
     "grade_id": "cell-035c561d4e99654c",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the equilibrium output level, price and the Herfindahl index. The Herhindahl index is defined as the sum of squared market shares:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:1}\n",
    "H = \\sum_j \\left( \\frac{q_j}{\\sum_i q_i} \\right)^{2}\n",
    "\\end{equation}\n",
    "\n",
    "If we have $n$ symmtric firms, we have $H = 1/n$. Hence, more competition in the form of more firms in the market leads to a lower Herfindahl index.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_guess_3 = [0,0,0]\n",
    "\n",
    "Q0 = np.sum(optimize.fsolve(lambda q: fixed_point_three_firms(q,vector_c), initial_guess_3))\n",
    "P0 = p(Q0)\n",
    "H0 = 3*(1.0/3.0)**2\n",
    "\n",
    "print(\"Before the merger\")\n",
    "print(\"=================\")\n",
    "print(\"total output: {:.3f}\".format(Q0))\n",
    "print(\"equil. price: {:.3f}\".format(P0))\n",
    "print(\"Herfn. index: {:.3f}\".format(H0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Define a function `fixed_point_two_firms` with the same\n",
    "structure as the function `fixed_point_three_firms` above, except that\n",
    "it derives the equilibrium output levels for a duopoly (two firms).\n",
    "Test this function by showing that each of the two firms produces\n",
    "0.3333 in case both firms have zero costs; use `fsolve` as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "018453f4503aaf72db01e0498e8d53cc",
     "grade": true,
     "grade_id": "cell-2618d929927d2965",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A competition authority (CA) is asked to evaluate the effects\n",
    "of a merger between firms 2 and 3. Firms 2 and 3 claim that by merging\n",
    "they can reduce their constant marginal costs. But it is not clear by\n",
    "how much they will reduce their costs.\n",
    "\n",
    "The CA assumes that the marginal cost level of the merged firm is\n",
    "uniformly distributed between 0 and the current marginal cost level\n",
    "`c0`. The merger will not affect the marginal cost level of firm 1 which\n",
    "does not merge. Firm 1's cost level remains `c0`.\n",
    "\n",
    "The next cell generates a vector of cost levels for the merged firm,\n",
    "denoted `c_after_merger`. Then it calculates the equilibrium output\n",
    "levels for (the non-merging) firm 1 and (the merged) firm 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_after_merger = pm.Uniform.dist(0,c0).random(size = 100)\n",
    "\n",
    "initial_guess = [0.2,0.2]\n",
    "\n",
    "q1_after_merger = [optimize.fsolve(lambda q: fixed_point_two_firms(q,[c0,c]), initial_guess)[0] for c in c_after_merger]\n",
    "q2_after_merger = [optimize.fsolve(lambda q: fixed_point_two_firms(q,[c0,c]), initial_guess)[1] for c in c_after_merger]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Create a dataframe called `df_after_merger` with\n",
    "three columns: `c_merged_firm`, `output_non_merging_firm`,\n",
    "`output_merged_firm` containing resp. the cost level of the merged firm,\n",
    "the output level of firm 1 and the output level of firm 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "7a3f8e887389e661046097c180e47f3a",
     "grade": true,
     "grade_id": "cell-cae9c13d60b629a2",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Add three columns to the dataframe with resp. total\n",
    "equilibrium output on the market, `Q`, equilibrium price, `P` and the\n",
    "Herfindahl index, `H`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0e13ca8ea6579acc4b2e695e918ccac3",
     "grade": true,
     "grade_id": "cell-12866a56fe44ae8f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Make a histogram of the equilibrium price `P` after\n",
    "the merger. Also indicate in the histogram the equilibrium price before\n",
    "the merger `P0`. Label the horizontal axis with $P$.\n",
    "\n",
    "[hint: you may want to use matplotlib's `hist`, `vlines` and `legend` to\n",
    "make this graph (e.g use google to find these functions); but feel free\n",
    "to use something else]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "3ed8d2cecc50f369f6926c8e4fcd9177",
     "grade": true,
     "grade_id": "cell-9aed05a5ad5cc7e4",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Excersise** Explain why sometimes the equilibrium price after\n",
    "the merger exceeds the equilibrium price before the merger and sometimes\n",
    "it is lower than the pre-merger price."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "586ab3657ceada4ab1375a999ca2a0f0",
     "grade": true,
     "grade_id": "cell-d720e381cb3cf3f2",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** What is calculated in the following cell?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(df_after_merger.P < P0)/len(df_after_merger.P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "a79811f5fe2ae30aa854f571a7404e8d",
     "grade": true,
     "grade_id": "cell-299f9ac6b26f5603",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Make a graph with the Herfindahl index on the\n",
    "horizontal axis and the equilibrium price on the vertical axis. This is\n",
    "straightforward for $(H,P)$ after the merger as both values are in the\n",
    "dataframe. Add in another color, the pre-merger combination `(H0,P0)`\n",
    "that we calculated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "3c75667ac0506c8c2748dfd2dabb3dbe",
     "grade": true,
     "grade_id": "cell-b31921136d5ee546",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** What does the figure above illustrate about the relation\n",
    "between the Herfindahl index and the equilibrium price? To illustrate,\n",
    "some people think that lower values of the Herfindahl index are\n",
    "associated with more competitive outcome. Would you agree with this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5aef5249c499644e02ce43dd49d4fc1e",
     "grade": true,
     "grade_id": "cell-db504e62423559ea",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### external effects\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A final reason why people are not always enthusiastic about markets is the presence of external effects. One can think of pollution associated with the production of a good. We model this as follows. Assume a monopolist can produce the product at cost $c q$. But production leads to an external effect equal to $\\gamma q$. Hence, the social cost of production equals $(c+\\gamma)q$\n",
    "\n",
    "We can model this as follows.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_agents = 1000\n",
    "valuations = np.array(sorted(pm.Normal.dist(100,20).random(size=number_of_agents),reverse = True))\n",
    "\n",
    "def demand(p,valuations):\n",
    "    return sum(valuations>p)\n",
    "\n",
    "c = 30\n",
    "γ = 80\n",
    "def costs(q):\n",
    "    return c*q\n",
    "\n",
    "def externality(q):\n",
    "    return γ*q\n",
    "\n",
    "def profit_c(p,valuations):\n",
    "    return p*demand(p,valuations)-costs(demand(p,valuations))\n",
    "\n",
    "def welfare_e(p,valuations):\n",
    "    return np.sum(valuations[:demand(p,valuations)])-costs(demand(p,valuations))-externality(demand(p,valuations))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Show graphically that the welfare maximizing price exceeds the profit maximizing price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b424370f3056e491285ae73e90403eff",
     "grade": true,
     "grade_id": "cell-ba3c1fb66fff8090",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** What is the interpretation of this result? Which policy instrument can the government use here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "f9262a5e319ec8483306a7fdb2bbcac6",
     "grade": true,
     "grade_id": "cell-5228df55ca31a43e",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Asymmetric information\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tirole chapter 15.\n",
    "\n",
    "One of the reasons why markets (or other institutions for that matter) work less well than a naive observer may think is asymmetric information. We consider here both adverse selection and moral hazard. Adverse selection we analyze in the context of insurance and moral hazard in the context of taxation.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## adverse selection\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** What is adverse selection?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "494aedd135688ce5b908c19cee1d618b",
     "grade": true,
     "grade_id": "cell-f5f53a947c76feab",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider an economy with `number_of_agents` agents. Each agent has an endowment/income equal to `income` and faces a potential loss of the size `cost`. Agents differ in the probability $\\pi$ of this loss. We randomly draw 50 values for $\\pi$ assuming it is uniformly distributed on $[0,1]$.\n",
    "\n",
    "Further, agents have a utility function of the form $u(x)=x^{\\rho}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "income = 1.1\n",
    "cost = 1\n",
    "ρ = 0.1\n",
    "def u(x):\n",
    "    return x**ρ\n",
    "\n",
    "number_of_agents = 50\n",
    "\n",
    "π = pm.Uniform.dist(0.0,1.0).random(size = number_of_agents)\n",
    "π.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we assume that $\\rho \\in \\langle 0, 1 \\rangle$, agents are risk averse and would like to buy insurance which covers the loss. We assume that insurance covers the loss completely at a premium $\\sigma$. As we assume that the probability of loss, $\\pi$, is exogenous, there is no reason to have co-payments of any sort.\n",
    "\n",
    "An agent buys insurance if and only if\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:2}\n",
    "u(\\text{income}-\\sigma) > \\pi u(\\text{income}-\\text{cost}) + (1-\\pi) u(\\text{income})\n",
    "\\end{equation}\n",
    "\n",
    "**Question** Define a function `insurance_demand` that returns the number of agents buying insurance as a function of the premium $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0531283149f9c2d3e2643c4e225981d2",
     "grade": true,
     "grade_id": "cell-c5f31d2eac1e1052",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume that this insurance market is perfectly competitive. That is, for each quantity supplied, insurance companies compete down the price such that the premium equals the average (expected) cost of the agents buying insurance.\n",
    "\n",
    "**Question** Explain the code of the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def insurance_supply(Q):\n",
    "    return np.mean(π[-Q:])*cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "682d856d1c6748cbf2ce658d0a7e2e3e",
     "grade": true,
     "grade_id": "cell-1c5a96bc85938b9d",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot demand and supply in one figure. In addition, we plot the marginal costs curve.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "range_Q = np.arange(1,number_of_agents+1,1)\n",
    "range_sigma = np.arange(0,1.01,0.01)\n",
    "plt.plot(range_Q,[insurance_supply(Q) for Q in range_Q],label=\"insurance supply\")\n",
    "plt.plot([insurance_demand(sigma) for sigma in range_sigma],range_sigma,label=\"insurance demand\")\n",
    "plt.plot(range_Q,[π[-Q]*cost for Q in range_Q],label=\"marginal cost\")\n",
    "plt.legend()\n",
    "plt.xlabel('$Q$')\n",
    "plt.ylabel('$\\sigma$')\n",
    "plt.title('Perfectly competitive insurance market')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Interpret this figure. In particular, \n",
    "\n",
    "-   explain why all curves are downward sloping (is supply not usually upward sloping?)\n",
    "-   what is approx. the equilibrium premium $\\sigma$?\n",
    "-   is the market outcome efficient?\n",
    "-   what can we learn from the marginal cost curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "9764a14341eb3c0f607426a9c7250cbf",
     "grade": true,
     "grade_id": "cell-f03a05b583f8d1c4",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assignment** Show graphically the effect of an increase in income on the market outcome. Does the inefficiency increase or decrease with income? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "c0fbc6ea9220168d2667f463f1cf5653",
     "grade": true,
     "grade_id": "cell-36a2a18ba8492e55",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b256b883e450fcd494d776a7015f7950",
     "grade": true,
     "grade_id": "cell-9dfa0a5b5186936a",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## moral hazard: optimal taxation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With moral hazard, agents take hidden actions. The actions that they take are affected by the incentives that they face. We consider this in the context of taxation. \n",
    "\n",
    "People differ in their productivity. For some people it is easy to generate a gross income $x$, for others generating such an income would be very costly in terms of effort. In the real world, such differences in productivity can be caused by IQ, education, health status etc. Here, we simply model this as an effort cost. People with a high effort cost have lower productivity than people with low effort costs. We assume that the effort cost is log-normally distributed. \n",
    "\n",
    "The government uses a linear tax schedule: $\\tau x - \\tau_0$. Hence, when you have a gross income $x$, your net income equals $(1-\\tau)x+\\tau_0$. Where we assume that for the economy as a whole the tax revenue is redistributed among the population. Hence, `number_of_agents` times $\\tau_0$ has to equal the total revenue from the marginal tax rate $\\tau$.\n",
    "\n",
    "Agents maximize their utility by choosing production $x$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:3}\n",
    "\\max_{x \\geq 0} (1-\\tau)x+\\tau_0 - cx^2\n",
    "\\end{equation}\n",
    "\n",
    "where agents differ in $c$ and $c$ is not observable.\n",
    "\n",
    "These two aspects are important: if $c$ were observable or if everyone was symmetric (had the same $c$) taxation would be easy. To see why, first note that income $x$ is apparently observable since taxation depends on it. Hence, the government could say to an agent $c$: I want you to produce income $x$ and you give me a share $\\tau$ of this income. \n",
    "\n",
    "In our set-up with heterogeneity in $c$ and $c$ unobservable, the government cannot force people to generate income $x$ because some of these agents may have such a high $c$ that this is inefficient (or even impossible).\n",
    "\n",
    "Hence, the government sets the tax schedule (in our case linear) and allows each agent to choose her own production level. The higher $\\tau$, the lower an agent's production will be.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_agents = 200\n",
    "effort_costs = pm.Lognormal.dist(mu=0.0,sd=0.5).random(size=number_of_agents)\n",
    "def effort(c,τ):\n",
    "    sol = optimize.minimize(lambda x: -(x*(1-τ)-c*x**2),1)\n",
    "    return sol.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the following welfare function:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:4}\n",
    "W = \\left(\\sum_i ( (1-\\tau)x_i + \\tau_0 - c_i x_i^2)^{\\rho} \\right)^{1/\\rho}\n",
    "\\end{equation}\n",
    "\n",
    "With $\\rho=1$, the social planner just maximizes the sum of utility. With $\\rho<1$, the planner has a taste for redistribution: agents with low utility get a relatively high weight in this welfare function.\n",
    "\n",
    "The function `Welfare` first calculates for a given $\\tau$, what the value of $\\tau_0$ is (using budget balance for the government). Then for these values of $\\tau$ and $\\tau_{0}$, $W$ is calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Welfare(τ,ρ):\n",
    "    τ_0 = np.mean([τ*effort(c,τ) for c in effort_costs])\n",
    "    return (np.sum([((1-τ)*effort(c,τ)+τ_0 - c*effort(c,τ)**2)**ρ for c in effort_costs]))**(1/ρ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Plot `Welfare` as a function of $\\tau$ for $\\rho=1$. What is the welfare maximizing tax rate? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "947bb341f001daa1260382c36639f944",
     "grade": true,
     "grade_id": "cell-fc051edb95253c67",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** What happens to the optimal tax rate as $\\rho<1$ falls?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1d19cf7443a72f8cbfceabac37b82566",
     "grade": true,
     "grade_id": "cell-2e3de756f51128b6",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "ec8c964bafef52863fa30b434696a230",
     "grade": true,
     "grade_id": "cell-8277fe6019dfa66a",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Assignment**** Redefine the function `Welfare` above such that it uses [Rawls' criterion](https://en.wikipedia.org/wiki/A_Theory_of_Justice) of maximizing the utility of the person who is worse off in society. Further, suppose that the government needs $g$ per head to finance a public good. What is the effect of $g$ on the optimal marginal tax rate?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "05bf44b1aa0b2561a37e154dfaaf2ec7",
     "grade": true,
     "grade_id": "cell-4f4524226874762e",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "111b0c29e7263b356c6a6c60df8ab7c4",
     "grade": true,
     "grade_id": "cell-f4b3d5d1a0994d2d",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Financial crisis\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tirole chapters: 11, 12\n",
    "\n",
    "We will look at two aspects of the financial crisis. First, why are financial markets problematic in the first place. Second, many people claim that the crisis was (partly) caused by the bonus contracts used by banks. Why do banks offer their employees such contracts?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Why is there a problem in financial markets?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first problem in financial (and other) markets is limited liability. When banks go bankrupt they \"only\" loose their equity even if the debts that they accumulated exceed their equity. We run some simulations to show that this leads to banks taking excessive risks from a social point of view.\n",
    "\n",
    "Let $x$ denote an investment opportunity: $x$ is a vector with 1000 elements. That is, we assume that there are 1000 states of the world and $x$ gives us the return in each of these states of the world. To find the expected `profit`, we take the average over the return in all these states of the world. However, if $x$ is \"very negative\" (a big loss), the bank goes bankrupt and the owners only loose their `equity`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def profit(x,equity=0):\n",
    "    return np.mean(np.maximum(x,-equity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the function `profit` for the case where `equity`=10 and $x$ is a scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1cb60f744345f3f0793cfc8ec506394d",
     "grade": true,
     "grade_id": "cell-63a92efa63b680e9",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that although we think of $x$ as a vector (and we will use this in a second), python does not know this nor does it care. Hence, we can plot `profit` as a function of the scalar $x$.\n",
    "\n",
    "\n",
    "**Question** Create an investment opportunity `vector_returns` where the returns are normally distributed with mean $-10$ and standard deviation 100. As mentioned above, we assume that there are 1000 states of the world."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5ddd69459130ac7f65c47df21b1f7c8e",
     "grade": true,
     "grade_id": "cell-5c9bdb02e8aa0f44",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Calculate the expected (social) return from this investment `vector_returns`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5b9fa0f91b3c38ab19b7fd297802c7e9",
     "grade": true,
     "grade_id": "cell-68f5cf65a7a7cbe9",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Calculate the expected profits of the `vector_returns`, assuming the bank has no equity. Compare the outcome to the one above. What is the interpretation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "ffea4bbf8f047086a406383bcac18637",
     "grade": true,
     "grade_id": "cell-f049260cce9706f5",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "801bc861e61b47e8a52a193f3818d5c3",
     "grade": true,
     "grade_id": "cell-7978647c2b6acd48",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Explain what is coded in the following code cell. E.g. what is the point of the 10000?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "cd0a140568a6dbeef3c9d9d93eed7f3f",
     "grade": true,
     "grade_id": "cell-051371459833eb4f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v_std = np.arange(0,200,1)\n",
    "v_returns = [pm.Normal.dist(-10,std).random(size=1000) for std in v_std]\n",
    "plt.scatter([np.std(vx) for vx in v_returns],[profit(vx) for vx in v_returns], label=\"no equity\")\n",
    "plt.scatter([np.std(vx) for vx in v_returns],[profit(vx,60) for vx in v_returns], label=\"equity equals 60\")\n",
    "plt.scatter([np.std(vx) for vx in v_returns],[profit(vx,10000) for vx in v_returns], label=\"social value\")\n",
    "plt.xlabel('$\\sigma$')\n",
    "plt.ylabel('return')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****Question**** Explain the economic intuition/interpretation of the graph above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0051308f83af582401a4ad7b44e71827",
     "grade": true,
     "grade_id": "cell-99aaf84ac349eb2f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Why these bonus contracts?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bonus contracts are generally a reaction to asymmetric information. Let's first consider moral hazard.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### moral hazard\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do your employees choose their investment opportunities? There are a number of things that you may worry about. Suppose you would pay everyone a fixed salary that does not depend on performance. Then your employees may randomly pick an investment opportunity without analyzing whether this is the best opportunity. The rest of the time, they simply enjoy the sun. Or they may invest your money in the new webshop of their brother in law. This may increase their status in their family but does not necessarily boost your profits. \n",
    "\n",
    "Since it is hard for banks to monitor exactly what investment opportunities their employees choose and how risky these are, it seems a good idea to give them some incentive to choose the right investment. One way to do this is to reward good outcomes. That is, the higher the return is, the higher their income will be. That is, employees get a bonus for good *outcomes*. Not for investments with a high expected outcome because that is hard to monitor.\n",
    "\n",
    "Tirole (2017) claims that before the financial crisis investments bank offered high bonuses to attract talented employees (page 345). These bonuses led to inefficient risky behaviour by these employees. But why should competition for talent lead to inefficiencies?\n",
    "\n",
    "We follow Bijlsma et al. (2018) to model this question. The structure of this model is comparable to our analysis of moral hazard in optimal taxation above. For a given bonus contract, employees choose the investment project that maximizes their own income. Knowing this, the bank sets the bonus contract to maximize its profits.\n",
    "\n",
    "Assume that there are 3 states of the world: the good state where the bank receives $y_g$ as return on the investment; a bad state where the bank makes a loss $y_b <0$ on the project and an \"average state\" where the bank earns $y_a \\in \\langle 0, y_g \\rangle$. The employee can choose from a set of investment projects that differ in their probabilities over these 3 states of the world. We model this as follows. \n",
    "\n",
    "Projects are indexed with their probability $q \\in [0,1]$ of the \"average state\". For a given $q$, the probabilities of the other two states are given by $q_g = \\alpha (1-q)(1+q)$ where $\\alpha \\in [0,0.5]$ denotes the \"talent\" of the employee and $q_b = 1- q - q_g$ resp. More talented agents (higher $\\alpha$) have a higher probability of the good state and a lower probability of the bad state for a given probability $q_a = q$.\n",
    "\n",
    "The bank cannot observe the project choice $q$ of the employee but it can observe and contract upon the outcome $y_{g,a,b}$. Hence, it can specify a wage for each state $w_{g,a,b}$. The limited liability of the agent is modelled as $w_g,w_a,w_b \\geq 0$. The bank cannot fine ($w<0$) its employee. \n",
    "\n",
    "**Question** What does a contract with a constant wage look like?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "43266eed6a1455e2af05c1523c49d8a0",
     "grade": true,
     "grade_id": "cell-f33b013685a0de5a",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify values for $y_{g,a,b}$, define the functions for $q_g$ and $q_b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_a = 1\n",
    "y_g = 10\n",
    "y_b = -20\n",
    "\n",
    "def q_g(q,ability):\n",
    "    return ability*(1-q)*(1+q)\n",
    "\n",
    "def q_b(q,ability):\n",
    "    return 1 - q - q_g(q,ability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Define a function `q_choice` of the wage vector $w$ and `ability` of the employee. [hint: the vector $w$ only needs to have 2 dimensions]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8104053f5c3a7e6f345fbd191bf21115",
     "grade": true,
     "grade_id": "cell-ece78fcea875fb09",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the function `q_choice`, the function `bank_choices` derives the optimal wage vector $w=[w_a,w_g]$ for a given `ability` of the employee. Then contract derives optimal risk choice $q$ as a function of `ability` and the outside option of the employee. If the outside option of the employee is so high that it is no longer profitable for the bank to hire this employee, the function returns $-1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_guess = [0.5,1.5]\n",
    "\n",
    "def bank_choices(ability):\n",
    "    opt_w = optimize.fmin(lambda x: -(q_g(q_choice(x,ability),ability)*y_g + q_choice(x,ability)*y_a + (1-q_choice(x,ability)-q_g(q_choice(x,ability),ability))*y_b),initial_guess, disp=0)\n",
    "    return [opt_w,q_choice(opt_w,ability)]\n",
    "\n",
    "def contract(ability,outside_option):\n",
    "    q = bank_choices(ability)[1]\n",
    "    profit = q*y_a + q_g(q,ability)*y_g + q_b(q,ability)*y_b\n",
    "    if profit - outside_option >= 0:\n",
    "        q_out = q\n",
    "    else:\n",
    "        q_out = -1\n",
    "    return q_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the optimal $q$ for the bank as a function of the outside option of an employee with `ability` equal to $0.5$. Do this for outside options between 0 and 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "46211f41b7682cb54bc7e295af219721",
     "grade": true,
     "grade_id": "cell-a24c6859bbc9c4f9",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Do we see that as competition for talented employees intensifies, thereby increasing their outside options, banks tend to offer them more risky contracts? Why (not)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0a6a6a2e94deff09b3fb7be6182dc906",
     "grade": true,
     "grade_id": "cell-2a431da00d421f09",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we add the second form of asymmetric information: adverse selection. Before we do this and in order to speed up the code below, we will also solve the problem above analytically.\n",
    "\n",
    "The bank solves the following optimization problem:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:5}\n",
    "\\max_q qy_a + \\alpha (1-q)(1+q)y_g + (1-q-\\alpha (1-q)(1+q)) y_b\n",
    "\\end{equation}\n",
    "\n",
    "**Question** By taking the first order condition for $q$, show that the bank would like to implement:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:6}\n",
    "q = \\frac{y_a-y_b}{2\\alpha(y_g-y_b)}\n",
    "\\end{equation}\n",
    "\n",
    "What is the optimal $q$ for a trader with $\\alpha=0.5$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "02647ef7d601deeda20684ad17207afd",
     "grade": true,
     "grade_id": "cell-1e42d1032ded906f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** For a social planner, the damage in the bad state ($y_b$) may be bigger (i.e. more negative) than for the bank. The bank only loses its equity, the government may need to bail out the bank or there can be a bank run if one bank collapses. What is the effect on the optimal $q$ of lower (more negative) $y_b$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "e2b4bfc36de182d4b4a9c5175b39ca30",
     "grade": true,
     "grade_id": "cell-04293c19ab6da9db",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** A less talented trader has lower $\\alpha$; what is the effect on $q$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8df827ad102eaf64985ea4f123f086ac",
     "grade": true,
     "grade_id": "cell-2c1af1a0d48f543e",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Use the equation for $q$ above to find the $q$ chosen by a worker facing a bonus contract with $w_a,w_g$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "524aa6adf3833d80af12126dc88b2701",
     "grade": true,
     "grade_id": "cell-d01571d5a976bbe8",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### moral hazard and adverse selection\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose now that there are two types of traders: one talented and the other less so. The bank cannot distinguish these types by just observing them. They may have the same education degrees and work experience. But some have a \"knack\" for observing opportunities and taking risks, which the others lack.\n",
    "\n",
    "The bank is willing to pay a lot to remunerate the top traders, but not the average ones. However, the average ones will try to look like the top traders to also earn these stellar incomes. To avoid the average ones to take these risks, the bank needs to pay them enough to stop them from mimicking the top traders.\n",
    "\n",
    "To simplify the analysis here, we make a number of assumptions:\n",
    "\n",
    "-   first, the fraction of top traders equals 0.5\n",
    "-   the outside option of the average traders is so low that top traders are never interested in mimicking average traders; but average traders do want to mimic top traders\n",
    "-   top traders are only paid in terms of bonus payments ($w_a,w_g$); not a fixed income component\n",
    "-   when we optimize over the top trader's wages below, we take the $q$ for the average traders and the profit they generate as given (and equal to their optimal profit for the bank)\n",
    "\n",
    "The last two points can actually be proved (in this sense, they are not assumptions); but we are not going to worry about this here.\n",
    "\n",
    "Since top traders only receive bonus payments, it must be the case that:\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:7}\n",
    "\\text{outside_option}=qw_a + \\alpha_h (1-q)(1+q) w_g\n",
    "\\end{equation}\n",
    "\n",
    "Defining the bonus ratio are $R=w_a/w_g$, we see that\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{eq:8}\n",
    "w_g = \\frac{\\text{outside_option}}{qR+\\alpha_h(1-q)(1+q)}\n",
    "\\end{equation}\n",
    "\n",
    "and hence we find that $w_a = R w_g$.\n",
    "\n",
    "The code below derives the optimal value of this bonus ratio $R$.\n",
    "\n",
    "**Question** Explain why $R$ can be interpreted as the \"riskiness\" of the bonus contract.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "ba5a6ba40e00c4b149aa2ef466d1688f",
     "grade": true,
     "grade_id": "cell-164cd4216160ee33",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "α_l = 0.1\n",
    "α_h = 0.5\n",
    "\n",
    "def profit(R,outside_option):\n",
    "    q = R/(2*α_h)\n",
    "    w_g = outside_option/(R*q + q_g(q,α_h))\n",
    "    mimic_q = R/(2*α_l)\n",
    "    w_a = R*w_g\n",
    "    wage_l = mimic_q*w_a+q_g(mimic_q,α_l)*w_g\n",
    "    profit = 0.5*(q*y_a + q_g(q,α_h)*y_g+(1-q-q_g(q,α_h))*y_b - outside_option) - 0.5*wage_l\n",
    "    return [profit, q]\n",
    "\n",
    "initial_guess = 0.5\n",
    "\n",
    "def outcome_h(outside_option):\n",
    "    wages_h = optimize.fmin(lambda x: -profit(x,outside_option)[0],initial_guess, disp=0)\n",
    "    return profit(wages_h,outside_option)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the optimal value of $q$ as a function of the outside option.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(range_outside_options,[outcome_h(o) for o in range_outside_options])\n",
    "plt.xlabel('outside option')\n",
    "plt.ylabel('probability of safe choice $q_a$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** What does the graph above show? What is the intuition for this?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "c018a99f25d265c9e97f41284fecf425",
     "grade": true,
     "grade_id": "cell-4e0009a614c239e4",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Using python for empirical research\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider two ways in which python can be useful for empirical research. First, the use of API's to download datasets. Second, the use of hacker statistics.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## API's to get data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A good reason to use python for data analysis is the option to get on-line data directly into your notebook without going to the website first to download this data. A number of institutes have such python API's. To illustrate this, we use the Worldbank API as described on [this website](http://wbdata.readthedocs.io/en/latest/).\n",
    "\n",
    "The advantage of doing your analysis in this way is that your research becomes better reproducible. Everyone can run the same code and then go through your code of data cleaning steps to end up with the same data set. If instead you first download the data to your computer, then use excel to clean the data and then start analyzing it (say, with stata), no one will be able to reproduce the exact steps that you have taken.\n",
    "\n",
    "To illustrate the Worldbank API, we will look at the development over time of inequality in gdp per head. So we want measures of gdp per head. The API allows us to search for such indicators in the Worldbank data set. The column on the left gives the name of the variables (that we will use below to download the data); the column on the right explains what the variable provides.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wb.search_indicators(\"gdp per capita\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say that we are interested in \"GDP per capita, PPP (constant 2011 international $)\", we specify this indicator in a dictionary where the key is the \"official name\" of the variable and the value is the way that we want to refer to the variable (in this case: \"GDP\\_per\\_head\").\n",
    "\n",
    "With `get_dataframe` we actually download the data into the dataframe `df_wb`. We reset the index in this case (just see what happens to the dataframe if you don't do this). And we look at the first 5 rows to get an idea of what the data are.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indicators = {\"NY.GDP.PCAP.PP.KD\": \"GDP_per_head\"}\n",
    "df_wb = wb.get_dataframe(indicators, convert_date=True)\n",
    "df_wb.reset_index(inplace = True)\n",
    "df_wb.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** What do the last 10 rows look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1f67f232a1b8e7c6059fee6873788e59",
     "grade": true,
     "grade_id": "cell-a7af556672d16f72",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you like the dataframe that you have downloaded from the web, you can save it with pandas `to_csv`. We save the data in the subdirectory \"data\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wb.to_csv('data/worldbank_data_gdp_per_capita.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the distribution of gdp per head in 1990 with the distribution in 2017. In order to illustrate how we can combine dataframes in pandas, we first define separate dataframes for the years 1990 and 2017.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_1990=df_wb[df_wb['date']=='1990-01-01']\n",
    "df_2017=df_wb[df_wb['date']=='2017-01-01']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** What does the dataframe `df_1990` look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "c84cf1343f4f04354f96430abfe67f30",
     "grade": true,
     "grade_id": "cell-0c59c12f078837e9",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both dataframes have a column `country`. Hence, we can merge the dataframes on this column. There are a number of `how` methods, here we use 'inner' which means that only countries that are present in both datasets will be in `df_merged`. To distinguish the columns, like `GDP_per_head` from the two dataframes, we can provide suffixes. All columns from `df_1990` (except for `country`) will be suffixed with '\\_1990'; and similarly for 2017."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_merged = pd.merge(df_1990, df_2017, on=['country'], suffixes=['_1990', '_2017'], how='inner')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** To see how the suffixes work, check what `df_merged` looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "76e86e98d28f2fe6cb5af429b7233384",
     "grade": true,
     "grade_id": "cell-60669260b6f8340a",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot GPD per head in 1990 against GDP per head in 2017. What do you conclude about the development in inequality in these 27 years?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b31d5b88e7ec639bf91668464e00f98d",
     "grade": true,
     "grade_id": "cell-642e222f3ff57506",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "a832758935e5df5ce7638a8a3ce7bcd2",
     "grade": true,
     "grade_id": "cell-d191195d398f91b2",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may wonder which observations (\"dots\") correspond to which countries. For this we need a plotting library that is more sophisticated on interactions than `matplotlib`. A number of these libraries are available; here we consider [bokeh](https://bokeh.pydata.org/en/latest/docs/user_guide/quickstart.html). If you want to know more about bokeh, there is a [datacamp course](https://www.datacamp.com/courses/interactive-data-visualization-with-bokeh)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_file, show, output_notebook\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.models import HoverTool\n",
    "output_notebook()\n",
    "\n",
    "hover = HoverTool(tooltips=[\n",
    "     ('country','@country'),\n",
    "     ])\n",
    "\n",
    "plot = figure(tools=[hover])\n",
    "plot.circle('GDP_per_head_1990','GDP_per_head_2017',\n",
    "    size=10, source=df_merged)\n",
    "output_file('data/inequality.html')\n",
    "show(plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Hacker statistics\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last chapter of the Datacamp course [Intermediate Python for Data Science](https://www.datacamp.com/courses/intermediate-python-for-data-science) you have seen hacker statistics. The idea is to simulate a random process, say, 10,000 times. Then if you want to know the probability that a certain condition is satisfied (say, number of heads bigger than 10 when you throw a coin 15 times), you calculate the number of outcomes where this condition is satisfied and divide it by 10,000.\n",
    "\n",
    "If you can program, you can recap all the statistics that you were taught (and probably forgot).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### high school puzzles\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember when you did statistics at high school? \"In a dark cupboard there are 10 red and 20 white socks&#x2026;\", \"you throw dice 20 times and &#x2026;\" etc.\n",
    "\n",
    "Then you had to do the most counter-intuitive mental jumping around to \"understand\" why the outcome was \"one third\".\n",
    "\n",
    "When you can program, there is no need for this anymore. Program the probability process, replicate it 10,000 times and see what the frequency of an event is.\n",
    "\n",
    "We will practice this, using the example given in this [Ted talk](https://www.ted.com/talks/peter_donnelly_shows_how_stats_fool_juries#t-479966).\n",
    "\n",
    "-   Experiment 1: you throw a coin a number of times and wait for the pattern \"HTH\" to emerge;\n",
    "-   Experiment 2: you throw a coin a number of times and wait for the pattern \"HTT\" to emerge.\n",
    "\n",
    "In expectation, which experiment is over sooner (i.e. takes fewer throws to see the required pattern)?\n",
    "\n",
    "**Question** Define a function `check` which checks whether a list `coin_throws` ends with the required `pattern` (also a list).\n",
    "\n",
    "[hint: in terms of python, you need to use your knowledge of \"slicing\" lists]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0502165ad810d135a04e5646f1820514",
     "grade": true,
     "grade_id": "cell-25a434bc9a339463",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "**Question** Explain what the following code block does:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes = ['H','T']\n",
    "\n",
    "def waiting_time(pattern):\n",
    "    throws = list(np.random.choice(outcomes,size=3))\n",
    "    while not check(throws,pattern):\n",
    "        throws.append(random.choice(outcomes))\n",
    "    return len(throws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "d623873f2d918af06a846ddbd176fc8c",
     "grade": true,
     "grade_id": "cell-f9179c446988b6ab",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Create an empty list `waiting_times_1`. Then run the first experiment 50,000 times. Record the waiting time for the pattern `['H','T','H']` to appear &#x2013;for each of these 50,000 experiments&#x2013; in the list `waiting_times_1`. Calculate the average waiting time and give a histogram for the waiting times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b8987bb0dd9edc319bb790cc8fdf7c87",
     "grade": true,
     "grade_id": "cell-03286f90c658d468",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Do the same as the previous question for experiment 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "162d038b88728a7c63b845b281831f7a",
     "grade": true,
     "grade_id": "cell-d3c731183a0ab43f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Which experiment leads to the lower waiting time in expectation? Do you have an intuition?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "7f3d8b372bf49bfffc27d31c21c7559e",
     "grade": true,
     "grade_id": "cell-1a494ead7b124ab8",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### statistics\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start simple. Suppose we have two normally distributed variables $x,y$ with $\\mu_x = 10, \\mu_y = 20, \\sigma_x = 3, \\sigma_y =4$. We are interested in the variable $z = x+y$. Suppose you forgot what you know about the sum of normal distributions and wondered what the expectation and standard deviation is of $z$.\n",
    "\n",
    "**Question** Simulate $z$ and calculate its mean and standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "f294d7d17e51d76275c0f596bbbb5e78",
     "grade": true,
     "grade_id": "cell-5f8e0c41f1935559",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, indeed, we know that $\\mu_z = \\mu_x + \\mu_y$ and $\\sigma_z = \\sqrt{\\sigma_x^2 + \\sigma_y^2}$\n",
    "\n",
    "**Question** Let $x$ be normally distributed with $\\mu_x=0,\\sigma_x=1$ and $y$ has a poisson distribution with $\\lambda_y=5$. What is the expectation and standard deviation of $z=xy$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "de62b39c219dac5493257bfb0ff975eb",
     "grade": true,
     "grade_id": "cell-12cc43c446a51e4f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many students find it to hard to think of the distribution of an average or the distribution of a standard deviation.\n",
    "\n",
    "**Question** Consider the following code block and try to understand what it does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 1000\n",
    "sd = 100\n",
    "number_of_samples=250\n",
    "\n",
    "def moments(n):\n",
    "    samples = pm.Normal.dist(mu,sd).random(size=(number_of_samples,n))\n",
    "    mus = samples.mean(axis=1)\n",
    "    std = mus.std()\n",
    "    return [mus,std]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Redefine the function `my_function(n)` such that it goes through the points in the figure below. \n",
    "\n",
    "[hint: you do not need to fit a function, just use your knowledge of statistics]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_function(n):\n",
    "    return 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "d171176ed6a7e89f92e07fe3e2133830",
     "grade": true,
     "grade_id": "cell-53085b93bd704d2d",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We plot the second element of the function `moments` against $n$ and the function `my_function`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "range_n = np.arange(2,1000)\n",
    "\n",
    "plt.plot(range_n,[moments(n)[1] for n in range_n], label='moments')\n",
    "plt.plot(range_n,[my_function(n) for n in range_n], label='my_function')\n",
    "plt.legend()\n",
    "plt.xlabel('$n$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Explain what the following distribution below is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(moments(10)[0],bins=30,label='$n=10$',density=True)\n",
    "plt.hist(moments(100)[0],bins=30,alpha=0.6,label='$n=100$',density=True)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "11a43cee0c5499907af09ae9b398377f",
     "grade": true,
     "grade_id": "cell-2c871de3f6e3514f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Compare the probability that the green, blue distribution resp. exceeds 1015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "50ce92ddddd8175673e3345a000ac76b",
     "grade": true,
     "grade_id": "cell-396020e48bc6f3ab",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Suppose you have a sample of 100 observations. The average of these observations equals 1020. Your hypothesis is that these observations were drawn from a normal distribution with mean 1000 and standard deviation 100. Would you reject this hypothesis?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "5d927a3fe01bd40b392538a7652aa786",
     "grade": true,
     "grade_id": "cell-e2b7e3e40e911207",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8115c99662703477a2fbbd91ffe27ed7",
     "grade": true,
     "grade_id": "cell-f9738869567aff41",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you like this approach, see [this video](https://www.youtube.com/watch?time_continue=1&v=ssVsVhZEQ9M) for more examples. There is also a free book (in the form of jupyter notebooks) to [recap your statistics](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Regulation in health care markets\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tirole: chapters 15, 17\n",
    "\n",
    "In this section, the main question is: does government regulation have an effect on markets and can we measure/quantify this?\n",
    "\n",
    "For this we consider the effect of an increase in the deductible $d$ in Dutch basic health insurance.\n",
    "\n",
    "Some institutional background:\n",
    "\n",
    "-   we focus on the basic health insurance market (i.e. we ignore the supplementary health insurance market)\n",
    "-   basic health insurance is mandatory in the Netherlands\n",
    "-   for people below the age of 18, health care is free of charge\n",
    "-   for people older than 18: pay the first $d$ euros of treatments per year yourself, treatments above $d$ are free of charge\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## simple theory\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following simple theoretical framework. People get\n",
    "offered at max. one treatment per year. They decide whether either to\n",
    "accept this treatment or to go without treatment.\n",
    "\n",
    "The figure below plots costs $c$ of treatment vs. values $v$ of treatments. Each\n",
    "treatment is a point in this figure; a combination of $c$ and $v$.\n",
    "\n",
    "The red/blue line is the out-of-pocket payment by an agent facing two deductible levels: 365 and 170 euro resp. Treatments above the red line are always accepted. The value exceeds the out-of-pocket payment for both deductibles. Treatments below the blue line are always rejected: even with the low deductible, the value is below the out-of-pocket payment. Treatments in the yellow area are accepted with the low deductible but are rejected with the high deductible. Hence, to quantify the effect of an increase in the deductible, we want to know the probability that treatments fall in the yellow area. The more treatments in the yellow area, the bigger the fall in health care expenditure in response to an increase in $d$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deductible(c,d):\n",
    "    return min(c,d)\n",
    "\n",
    "range_c = np.arange(0,500,0.1)\n",
    "range_v170 = [deductible(c,170) for c in range_c]\n",
    "range_v365 = [deductible(c,365) for c in range_c]\n",
    "\n",
    "plt.plot(range_c,range_v365,'-', color = 'r', linewidth = 2, label = '$d=365$')\n",
    "plt.plot(range_c,range_v170,'-', color = 'b', linewidth = 2, label = '$d=170$')\n",
    "plt.legend()\n",
    "plt.fill_between(range_c, range_v170, range_v365, facecolor='yellow')\n",
    "plt.ylim(0,500)\n",
    "plt.xlabel('Cost')\n",
    "plt.ylabel('Value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## some data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the effect of an increase in deductible, we compare health care expenditures in the years 2011 (deductible was 170 euro) and 2014 (deductible was 365 euro). We use data from [Vektis](https://www.vektis.nl/intelligence/open-data). Download from this website the 'csv' files for 2011 and 2014. To use the code below, download these csv-files into the sub-directory \"data\" (i.e. \"data\" is sub-directory of the directory in which this notebook resides). \n",
    "\n",
    "When you open the csv files, you can see that it uses \";\" as separator between columns. Hence, we use pandas' `read_csv` statement where we specify the separator as ';'. The data contain total cost per postal code area for a number of cost categories. The expenditures are grouped by sex and age.\n",
    "\n",
    "The function `get_data_into_shape` does a number of things:\n",
    "\n",
    "-   not all health care cost categories in the data \"count\" as far as the deductible is concerned. Hence, we select the ones that fall under the deductible and sum these as the relevant total expenditure under the deductible.\n",
    "-   the cost categories are in Dutch, hence we translate the labels into English\n",
    "-   we drop variables that we do not need for the analysis here\n",
    "-   we calculate cost per head per postal code area\n",
    "-   we also introduce the log of health care costs per head\n",
    "-   we turn the variable `sex` into a category with two values ('M' for males, 'V' for females)\n",
    "-   we drop the age category '91+' and turn the remaining ages into integers\n",
    "-   finally, the function returns this new dataframe.\n",
    "\n",
    "The function illustrates the data manipulation you can do with pandas. We use the same function for the 2011 and 2014 data, which means that we ignore some columns in the 2014 data that are not available in the 2011 data.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2014 = pd.read_csv('data/Vektis Open Databestand Zorgverzekeringswet 2014 - postcode3.csv', sep = ';')\n",
    "\n",
    "cost_categories_under_deductible = ['KOSTEN_MEDISCH_SPECIALISTISCHE_ZORG', 'KOSTEN_MONDZORG', 'KOSTEN_FARMACIE', 'KOSTEN_HULPMIDDELEN', 'KOSTEN_PARAMEDISCHE_ZORG_FYSIOTHERAPIE', 'KOSTEN_PARAMEDISCHE_ZORG_OVERIG', 'KOSTEN_ZIEKENVERVOER_ZITTEND', 'KOSTEN_ZIEKENVERVOER_LIGGEND', 'KOSTEN_GRENSOVERSCHRIJDENDE_ZORG', 'KOSTEN_OVERIG']\n",
    "\n",
    "def get_data_into_shape(df):\n",
    "    df['health_expenditure_under_deductible'] = df[cost_categories_under_deductible].sum(axis=1)\n",
    "    df = df.rename({\n",
    "        'GESLACHT':'sex',\n",
    "        'LEEFTIJDSKLASSE':'age',\n",
    "        'GEMEENTENAAM':'MUNICIPALITY',\n",
    "        'AANTAL_BSN':'number_citizens',\n",
    "        'KOSTEN_MEDISCH_SPECIALISTISCHE_ZORG':'hospital_care',\n",
    "        'KOSTEN_FARMACIE':'pharmaceuticals',\n",
    "        'KOSTEN_TWEEDELIJNS_GGZ':'mental_care',\n",
    "        'KOSTEN_HUISARTS_INSCHRIJFTARIEF':'GP_capitation',\n",
    "        'KOSTEN_HUISARTS_CONSULT':'GP_fee_for_service',\n",
    "        'KOSTEN_HUISARTS_OVERIG':'GP_other',\n",
    "        'KOSTEN_MONDZORG':'dental care',\n",
    "        'KOSTEN_PARAMEDISCHE_ZORG_FYSIOTHERAPIE':'physiotherapy',\n",
    "        'KOSTEN_KRAAMZORG':'maternity_care',\n",
    "        'KOSTEN_VERLOSKUNDIGE_ZORG':'obstetrics'\n",
    "    }, axis='columns')\n",
    "    df.drop(['AANTAL_VERZEKERDEJAREN',\n",
    "             'KOSTEN_HULPMIDDELEN',\n",
    "             'KOSTEN_PARAMEDISCHE_ZORG_OVERIG',\n",
    "             'KOSTEN_ZIEKENVERVOER_ZITTEND',\n",
    "             'KOSTEN_ZIEKENVERVOER_LIGGEND',\n",
    "             'KOSTEN_GRENSOVERSCHRIJDENDE_ZORG',\n",
    "             'KOSTEN_OVERIG',\n",
    "             'KOSTEN_EERSTELIJNS_ONDERSTEUNING'],inplace=True,axis=1)\n",
    "    df.drop(df.index[[0]], inplace=True)\n",
    "    df['sex'] = df['sex'].astype('category')\n",
    "    df['age'] = df['age'].astype('category')\n",
    "    df['costs_per_head']=df['health_expenditure_under_deductible']/df['number_citizens']\n",
    "    df['log_costs_per_head']=np.log(1+df['health_expenditure_under_deductible']/df['number_citizens'])\n",
    "    df = df[(df['age'] != '90+')]\n",
    "    #df['age'] = df['age'].astype(int) # if this line gives an error, use the following:\n",
    "    df['age'] = df['age'].astype(\"Float64\").astype(int)\n",
    "    return df\n",
    "\n",
    "df_2014 = get_data_into_shape(df_2014)\n",
    "df_2014.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Create `costs_per_sex_age` which contains the average health care expenditure (averaged over postal code areas) for each combination of sex and age in the data for 2014. \n",
    "\n",
    "[hint: use `pandas` method `groupby`]\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "48cb63b1f29a31ad9adf272990148583",
     "grade": true,
     "grade_id": "cell-cac2f8b2969ceb62",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## matplotlib\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the average health care expenditure per head against age (on horizontal axis) for males and females.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "7a7d684648fce28d3bbe5b45c1b6b0f3",
     "grade": true,
     "grade_id": "cell-930df9bee293b57d",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Can you interpret how these costs evolve with age and sex?\n",
    "\n",
    "**Question** How can a graph like this help us to determine the effect of $d$ on health care expenditure?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## reversing the probability distributions\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we used `pymc3` to generate vectors of productivities, valuations, incomes etc. using probability distributions. Here we go the \"other way around\". We have here distributions of health care expenditures per head and we want to identify the distributions where these come from. To illustrate this, consider the distribution of (average) costs for 30 year old males. Since, health care costs have a skewed distribution, we actually plot the distribution of log costs.\n",
    "\n",
    "**Question** Plot the log health care cost distribution for 30 year old males.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0bba992bcf08c26ccda3b9e3268f0abb",
     "grade": true,
     "grade_id": "cell-9a0bcb2119fff275",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution looks (sort of) normal. Hence, we assume that for each age and sex category `log_costs_per_head` are normally distributed. This implies that `costs_per_head` have a log-normal distribution.\n",
    "\n",
    "**Question** Why are there no zeroes in this figure? Lots of people have no health care expenditure in a year? How does our `log_costs_per_head` deal with zero costs? Should there be a spike at 0?\n",
    "\n",
    "We focus here on health care costs for women. Clearly, a similar analysis can be done for men. In fact, it is also possible to combine men and women into one analysis with gender fixed effects.\n",
    "\n",
    "Here we focus on women and introduce age-fixed effects. We assume that observed costs $z$ are `log_costs_per_head` which are normally distributed with a mean $\\mu$ and standard deviation $\\sigma$ which both vary with age. We do not know these means and standard deviations `μ[age], σ[age]` but assume they are drawn from prior distributions. A [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) for $\\mu$ and a half-normal distribution for $\\sigma$.\n",
    "\n",
    "This is called a Bayesian analysis which you probably never saw before. Do not worry about this; you do not need to understand the details of this analysis. If you find it interesting, there are some references at the end.\n",
    "\n",
    "\n",
    "The next cell does the following:\n",
    "* we create the vector `log_costs_per_age_female` so that we can plot the prediction of the model against the observed average values\n",
    "* `log_costs_per_head` is the vector with values that we wish to explain with our model (the \"y-values\")\n",
    "* we create an index variable `age`; just evaluate `age` if you want to know what it looks like.\n",
    "* then we create our model:\n",
    "    * $\\mu$ is a vector of fixed effects with a different value for each age category. We \"assume\" (prior) that $\\mu$ is drawn from a normal distribution with expectation 8 and standard deviation 3\n",
    "    * $\\sigma$ is also a vector with fixed effects with a different value for each age category. As a standard deviation is always positive, we assume that the distribution of $\\sigma$ is half-normal.\n",
    "    * finally, we specify our variable `z` as being normally distributed with paramters $\\mu[age],\\sigma[age]$ for an individual of age `age`. We specify that the observed values (that we want to \"match\") are `log_costs_per_head`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_costs_per_age_female_2014 = df_2014[df_2014['sex']=='V'].groupby(['age'])['log_costs_per_head'].mean()\n",
    "\n",
    "log_costs_per_head_2014 = df_2014[df_2014['sex']=='V'].log_costs_per_head.values\n",
    "age_2014 = df_2014[df_2014['sex']=='V'].age.values\n",
    "\n",
    "\n",
    "with pm.Model() as model_2014:\n",
    "    \n",
    "    μ = pm.Normal('μ', 8, 3, shape=len(set(age_2014)))\n",
    "    σ = pm.HalfNormal('σ', 4, shape=len(set(age_2014)))\n",
    "    z = pm.Normal('z', μ[age_2014], σ[age_2014], observed=log_costs_per_head_2014)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model has our prior information (the distributions for $\\mu$ and $\\sigma$ that we specify) and the data. It combines these to generate the posterior distribution (using Bayes' formula).\n",
    "\n",
    "The posterior distribution is generated in the form of samples. Usually, you can easily generate 5000 samples, but since we are all doing this together on the server, it would slow things quite a bit. So we choose a lower number of samples here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 500\n",
    "\n",
    "with model_2014:\n",
    "    trace_2014 = pm.sample(n_samples,step = pm.Metropolis(),start = pm.find_MAP())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the posterior distributions for each of the parameters. Here, we plot it only for $\\mu$ for each age.\n",
    "\n",
    "The posterior shows the uncertainty that is left for each parameter after taking into account the priors and the data. The uncertainty can be reduced by taking more samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.plot_posterior(trace_2014, varnames=['μ'],ref_val = log_costs_per_age_female_2014.values)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figures above compare the estimated distribution for $\\mu$ for each age category with the observed average expenditure for this age category in the data. The observed average expenditure falls in the middle of this distribution for each age category. This suggests that the estimation makes some sense (although more tests should be done). \n",
    "\n",
    "The \"fun\" of Bayesian analysis is that we have a (posterior) distribution of the population parameters $\\mu$. With \"classic\" econometrics this is not possible, since $\\mu$ is not random in that analysis.\n",
    "\n",
    "To elaborate a little, in classical econometric analysis we consider the probability $p(\\text{data}|\\theta)$ where $\\theta$ denotes the parameters of the underlying distribution. Hence, we do $p(m>10|\\mu=9)$. With Bayes, we do $p(\\theta|\\text{data})$. The idea is that the data is fixed (we have observed it) but we face uncertainty over the values of the parameters $\\theta$. The graphs above plot for each age the distribution of $\\mu$.\n",
    "\n",
    "Now we can plot the average $\\mu$ for each age and the observed average expenditure per age category in a graph. \n",
    "\n",
    "To get an idea what `trace_2014`['$\\mu$'] looks like:\n",
    "\n",
    "**Question** Use the `shape` attribute on the trace-vector.\n",
    "\n",
    "[hint: what is the difference between an attribute and a method?]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "d337122449d1ce1a384a165c653b2596",
     "grade": true,
     "grade_id": "cell-750bb7ef3857e5ab",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method `np.mean()` takes the mean of a vector. As you see above, the trace has two dimensions: samples and age. We want to take the mean over the first dimension to get average expenditures across age categories.\n",
    "\n",
    "If `x` is a two dimensional vector, we take the mean over the first dimension by `np.mean(x,axis=0)`.\n",
    "\n",
    "**Question** Plot in one figure the predicted means per category (using the trace) and the observed means per category (using `log_costs_per_age_female_2014`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "fcc3dafd56a3aff1f4c4cc84e625fb72",
     "grade": true,
     "grade_id": "cell-f610a10af3a8812f",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the effect of the deductible, we compare the average $\\mu$ for 17 year olds with the average $\\mu$ for 19 year olds.\n",
    "\n",
    "**Question** Use indexing on `np.mean(trace_2014['μ'],axis=0)` to calculate the average difference in predicted expenditures for 17 and 19 year olds.\n",
    "\n",
    "[hint: python starts indexing at 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "19bac89a277af90cd104f2d67cc9591c",
     "grade": true,
     "grade_id": "cell-9d4b8675b922cda9",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is positive: the $\\mu$ parameter in the normal distribution of log expenditures is higher for 17 year olds than for 19 year olds. Now you may think that although the average $\\mu$ is higher for 17 year olds than for 19 year olds, perhaps it is still likely that $\\mu_{19}>\\mu_{17}$. To verify this, we plot these $\\mu$ distributions for 17 and 19 year olds to see whether they overlap.\n",
    "\n",
    "**Question** Use \"slicing\" to plot the histograms for $\\mu$ for 17 and 19 year olds.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "eec6331151a85338a7f6464e8699baec",
     "grade": true,
     "grade_id": "cell-0eeea179cc34d794",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Does this mean that expenditure of every 17 year old exceeds the expenditure of every 19 year old?\n",
    "\n",
    "**Question** Focus on 19 year olds. What is the probability that $\\mu_{19}>6.45$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8cc471e8f5f53eaa650f92a7c82469c2",
     "grade": true,
     "grade_id": "cell-45fffbb335301eb3",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(trace_2014['μ'][:,19]>6.45)/len(trace_2014['μ'][:,19])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One explanation for the figure is that 17 year olds do not face a deductible, while 19 year olds do. But an other explanation is that health care expenditure simply differs by age (irrespective of a deductible). In order to control for the age effect, we re-do the analysis above for 2011. Also in 2011 we can take the difference in means for 17 and 19 year olds. If there is (only) a biological reason for the different expenditures between 17 and 19 year olds, the difference in 2011 should be the same as the difference in 2014. \n",
    "\n",
    "If, however, the difference in expenditures is caused by the deductible, we expect a bigger difference in 2014 than in 2011 as the deductible was higher in 2014 than in 2011. In terms of our model above: the yellow area is the additional effect due to the higher deductible in 2014.\n",
    "\n",
    "Hence, we do the same analysis as above for 2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2011 = pd.read_csv('data/Vektis Open Databestand Zorgverzekeringswet 2011 - postcode3.csv', sep = ';')\n",
    "\n",
    "df_2011 = get_data_into_shape(df_2011)\n",
    "df_2011.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Specify the model for 2011, call this `model_2011`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "2ca62090c1e6d045d437aefbf6e44e9a",
     "grade": true,
     "grade_id": "cell-6af493086ffedb93",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Sample from `model_2011` and call the trace `trace_2011`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b956932c577aeae7b43e506aa41aa843",
     "grade": true,
     "grade_id": "cell-05c6508591f4529a",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that in 2014 the difference in mean log expenditure equals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(trace_2014['μ'],axis=0)[17]-np.mean(trace_2014['μ'],axis=0)[19]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Show that this difference is smaller in 2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "e194a05819c9b6883a47d07bdde26e93",
     "grade": true,
     "grade_id": "cell-1467054a0864f722",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the distributions (for 2011 and 2014) of the difference between $\\mu$ for 17 and 19 year olds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "a59e7fbbb4d895686e3a37133c18394c",
     "grade": true,
     "grade_id": "cell-40e3f8493be198a3",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Calculate the probability that in fact the parameter difference $\\mu_{17}-\\mu_{19}$ is bigger in 2014 than in 2011."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "9c25a9016c2a0e35d1583c8e8297ee8a",
     "grade": true,
     "grade_id": "cell-f3ee376c6123c607",
     "locked": false,
     "points": 1,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence we find that government regulation in terms of a deductible has an effect on health care expenditures. People spend less on health care if they face a higher deductible.\n",
    "\n",
    "One can extend the analysis above to make it more convincing. Think of things like\n",
    "\n",
    "-   add more years (using the Vektis website)\n",
    "-   add year dummies to distinguish year effects (like changes in treatments covered by basic insurance) from changes in deductible\n",
    "-   add the level of the deductible in the estimation.\n",
    "\n",
    "If you like this way of modelling with pymc3, you can look at [this video](https://www.youtube.com/watch?v=TMmSESkhRtI&t=9076s)\n",
    "\n",
    "A simple introduction to Bayesian estimation with python can be found [here](https://www.youtube.com/watch?v=TpgiFIGXcT4&t=6s). Allen Downey also has [free books](http://www.allendowney.com/wp/books/) on statistics with python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10"
  },
  "org": null,
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}