{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DiscreteDP Example: Bioeconomic Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Daisuke Oyama**\n",
    "\n",
    "*Faculty of Economics, University of Tokyo*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Miranda and Fackler, <i>Applied Computational Economics and Finance</i>, 2002,\n",
    "Section 7.6.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import quantecon as qe\n",
    "from quantecon.markov import DiscreteDP, backward_induction, sa_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emax = 8             # Energy carrying capacity\n",
    "n = emax + 1         # Number of states, 0, ..., emax\n",
    "m = 3                # Number of areas (actions), 0, ..., m-1\n",
    "e = [2, 4, 5]        # Energy offerings\n",
    "p = [1.0, 0.7, 0.8]  # Survival probabilities\n",
    "q = [0.5, 0.8, 0.7]  # Success probabilities\n",
    "\n",
    "T = 10               # Time horizon\n",
    "\n",
    "# Terminal values\n",
    "v_term = np.ones(n)\n",
    "v_term[0]= 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We follow the state-action pairs formulation approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = n * m  # Number of feasible state-action pairs\n",
    "s_indices, a_indices = sa_indices(n, m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reward vector\n",
    "R = np.zeros(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Transition probability array\n",
    "Q = sp.lil_matrix((L, n))\n",
    "it = np.nditer((s_indices, a_indices), flags=['c_index'])\n",
    "for s, k in it:\n",
    "    i = it.index\n",
    "    if s == 0:\n",
    "        Q[i, 0] = 1\n",
    "    elif s == 1:\n",
    "        Q[i, np.minimum(emax, s-1+e[k])] = p[k] * q[k]\n",
    "        Q[i, 0] = 1 - p[k] * q[k]\n",
    "    else:\n",
    "        Q[i, np.minimum(emax, s-1+e[k])] = p[k] * q[k]\n",
    "        Q[i, s-1] = p[k] * (1 - q[k])\n",
    "        Q[i, 0] = 1 - p[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Discount factor\n",
    "beta = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.6/site-packages/quantecon/markov/ddp.py:416: UserWarning: infinite horizon solution methods are disabled with beta=1\n",
      "  warnings.warn(msg)\n"
     ]
    }
   ],
   "source": [
    "# Create a DiscreteDP\n",
    "ddp = DiscreteDP(R, Q, beta, s_indices, a_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us use the `backward_induction` routine to solve our finite-horizon problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vs, sigmas = backward_induction(ddp, T, v_term)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAEWCAYAAABG5QDSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcnQV97/HP10TKLipRA2FTEYpaUSPYqoh1A0Sw1faC\nW7Uq5V612sWC7bXqdUFq25fVaikFBAXBDS1FrLjhjhIWQUBsZEsCKUFcQK0Q/N0/nmfk5DiTOTPP\nnHMmmc/79ZpXzrP/zpnJd37zrKkqJEmSJM3OvcZdgCRJkrQps6GWJEmSOrChliRJkjqwoZYkSZI6\nsKGWJEmSOrChliRJkjqwodZAkpyQ5A1zsJ5Tk7x1LmoaYFuV5KGzXPb6JE+bYtqTklwz2bxJ/jrJ\nSbOreHSS7JrkjiSLZrHsgUlWD6MuSXPDzN5gmpltZg+dDfUmLMkTk3w9yY+T3Jbka0keN4xtVdXR\nVfWWYax7QpKXJLm7DY2fJLksyaHD3OZsVNVXqmqvKaa9vapeDpBk9/YXxOK52G77i+3O9vO5Lcln\nk+w9m3VV1Y1VtW1V3T0XtfXVuXuSLyb5WZLvTvVLTlpozOzxMLOnrfP6JD9v67wjyflzvY2FwIZ6\nE5Vke+Bc4D3A/YCdgTcDv5jFupJkvvwsfKOqtgV2AE4GPpLkvv0zzVXgbYL+rv18lgG3AKfOdAUj\n+OzOBC4F7g/8DfCxJEuGvE1pXjOzzWzmb2YDPLtt2LetqmeMYHubnfnyH1Iz9zCAqjqzqu6uqp9X\n1flVdTlAkjclOX1i5v6/vJNckORtSb4G/Ax4XZIVvRtI8mdJzmlf/+qwX5Kre/dCJFmcZF2Sx7TD\nH02ytt0L8+UkD5/pm6uqXwKnAFsBD5k4ZJXkmCRrgfe323pFkpXtX//nJNmpb1WHJLk2ya1J3jnx\nSyjJQ5J8IckP2mlnJNmhb9nHJbkqyQ+TvD/Jlu2yUx4+6/vcv9z++6P2r/4nt3U+smf+B7R7cmfU\ncFbVz4APAY9o13OvJMcm+X77nj6S5H7ttInv/cuS3Ah8YZKfh53az++29vN8RU+NW7Xf/x8muQqY\nco9akocBjwHe2P5Mfhy4HHjuTN6ftBkyszGzmWeZrbljQ73p+h5wd5LTkhycSfYIDOBFwFHAdsAJ\nwF5J9uyZ/nyaAOh3JnBkz/AzgVur6pJ2+NPAnsADgEuAM2ZaWBsaLwfuAP6rHf0gmj07uwFHJfld\n4DjgD4GlwA3AWX2r+j1gOU2TdzjwxxObaJfdCfhNYBfgTX3LvqB9bw+h+WX4f2f4Ng5o/92h/av/\nS219L+yZ50jg81W1biYrTrJtW9+l7ahXA88Bnkzznn4IvLdvsSfTvNdnTrLKs4DV7bLPA97efr4A\nb6T5DB7SLvtHGynt4cC1VXV7z7hvt+OlhczMNrPnY2ZPOKP9I+v8JI8a9H2pR1X5tYl+0fxHO5Xm\nP9V64Bzgge20NwGn98y7O1DA4nb4AuD/9a3vdOBv29d7ArcDW7fDpwJvbV8/tG/aGRPLTVLjDu12\n79O/nknmfUn7Pn4E3ApcCDytnXYgcCewZc/8J9McTpsY3ha4C9i9HS7goJ7p/4cmCCfb9nOAS3uG\nrweO7hk+BPh+Ty2r++adqPNXn3v/Z96O2x+4EUg7vAL4wwG/36cC/9N+Pmvb7/dD2mlXA0/tmXdp\n+1ks7qnjwZP9PND8Yrob2K5n+nHAqe3ra/s+x6N6339fjS8CLuwb97aJdfnl10L+wsw2s+dZZrfT\nn0BzZGFr4PVtrTuM+//LpvblHupNWFVdXVUvqaplNIeRdgLeNYNVrOob/hD37MV4PvDJag5T9W93\nJU0YPDvJ1sBh7bIkWZTkHe1hrJ/QBBfAjgPWdGFV7VBVO1bV46vqcz3T1lXV//QM70Szh2OirjuA\nH9CcmzjZe7yhXYYkD0xyVpI1bZ2nT1LjpMt2UVXfpDlce2Cai1MeShOyg/r79vN5UFUdVlXfb8fv\nBnwiyY+S/Ijm+3M38MCeZfu/3xN2Am6rDfcq38A9n+NO/PpnMZU7gO37xt2H5pe5tKCZ2Wb2PMxs\nqupr1ZyC9LOqOo7mD4AnDfTu9Cs21JuJqvouzV/Dj2hH/ZTmr80JD5pssb7hzwJLkuxLE9KTHTqc\nMHEI8XDgqjawoQn1w4Gn0TRSu7fjM8j7mEZ/vTfRhFKzgWQbmgvh1vTMs0vP613bZQDe3q7vkVW1\nPc0hvf4ap1p2tvVOOK3d3ouAj/X9wpmtVcDBbXBPfG1ZVb2fxVT13ATcL8l2PeN25Z7P8WZ+/bOY\nypXAg/vW9ah2vKSWmW1mMz8yezLF3Hz/FxQb6k1Ukr2T/EWSZe3wLjRheWE7y2XAAWnuXXkfmsM4\nG1VVdwEfBd5Jc97bZzcy+1nAM4D/zYYhvh3NVes/oPnl8PaZvK8ZOhN4aZJ9k/xGu61vVtX1PfO8\nLsl928/nNcCHe+q8A/hxkp2B102y/lcmWdZeKPI3PcsOah3wS+DBfeNPpzlP8IXAB3ontBedHDjD\n7UBzPuXbkuzWrmdJksMHWbCqVgFfB45LsmWS3wJe1tYJ8BHg9e3nuIzm3L+p1vU9mp+9N7br+n3g\nkcDHZ/GepM2GmQ2Y2b3mRWa3P29PSLJFu67X0ez5/9os3tOCZkO96bqd5tyubyb5KU0ofwf4C4Cq\n+ixNmFwOXExzu6ZBfIhmT8VHq2r9VDNV1c3AN4DfYcPQ+gDN4aU1wFXc88tizrWHFt9A06zdTHMB\nxhF9s/07zfu/DPgUzTl80Nyu6jHAj9vxZ0+yiQ8B59Ocj/Z9YEYPN2gPvb4N+Fp7WO/x7fhVNBf+\nFPCVifnbXyC3A1fMZDutf6I5DHl+kttpPvf9Z7D8kTR7pm4CPkFzl46JQ7dvpvmeXkfzeXxwmnUd\nQXNR0Q9pzut7Xs3wAh5pM2Rmm9m95ktmbwf8C01erwEOotlz/oMZ1CLuOcle0gglOQW4qar+b8+4\nFwIPr6pp90xJkkbHzNZ0bKilEUuyO83el0dX1XXjrUaStDFmtgYxtFM+kpyS5JYk35liepK8O80N\nyS9Pe4N5aXOW5C00h3nfaTBrvjG3pQ2Z2RrU0PZQJzmA5gKCD1TVIyaZfgjNifKH0Jw39E9VNZPz\nhyRJc8jclqTZGdoe6qr6MnDbRmY5nCa0q6ouBHZIsnRY9UiSNs7clqTZWTzGbe/MhjceX92Ou7l/\nxiRH0Tzph2222eaxe++990gKlKS5dvHFF99aVUvGXccsDZTbZrY0c9++/HLW33XXuMuABObD9XXz\npQ4YKLPH2VAPrKpOBE4EWL58ea1YsWLMFUnS7CTZ6FPLNgdmtjRzSdjtmEHvljg8Nxx/qHVsWMdA\nmT3O+1CvYcMn+Sxjw6clSZLmF3NbkiYxzob6HODF7VXjjwd+3N54XpLm3NJlu5Jk7F+bOHNbkiYx\ntFM+kpwJHAjsmGQ18Ebg3gBVdQJwHs2V4iuBnwEvHVYtksZn6bJdWbtm1fQzjsA8OXw47hKmZG5L\n0uwMraGuqiOnmV7AK4e1fUnzw9o1q2xkNxHmtiTNziZxUaIkSdKwzKcjado02VBLkqQFbT4cSfMo\n2qZtnBclSpIkSZs8G2pJkiSpAxtqSZIkqQPPoZY2U15kI0nSaNhQS5up+XCRDXihjaSp+Ye/Nhc2\n1JIkaSz8w1+bC8+hliRJkjqwoZYkSZI6sKGWJEmSOrChliRJkjrwokRpjnnVuiRJC4sNtTTHvGpd\nkqSFxVM+JEmSpA5sqCVJkqQObKglSZKkDmyoJUmSpA5sqCVJkqQObKglSZKkDmyoJUmSpA5sqCVJ\nkqQObKglSZKkDmyoJUmSpA5sqCVJkqQObKglSZKkDmyoJUmSpA4Wj7sAaa4sXbYra9esGncZkiRp\ngbGh1mZj7ZpV7HbMueMugxuOP3TcJUiSpBHylA9JkiSpA/dQS5K0wHiKnDS3bKglSVpgPEVOmlue\n8iFJkiR1YEMtSZIkdWBDLUmSJHUw1IY6yUFJrkmyMsmxk0y/T5L/SPLtJFcmeekw65EkTc3MlqTZ\nGVpDnWQR8F7gYGAf4Mgk+/TN9krgqqp6FHAg8A9JthhWTZKkyZnZkjR7w9xDvR+wsqqurao7gbOA\nw/vmKWC7JAG2BW4D1g+xJknS5MxsSZqlYTbUOwO9N7lc3Y7r9c/AbwI3AVcAr6mqX/avKMlRSVYk\nWbFu3bph1StJC5mZLUmzNO6LEp8JXAbsBOwL/HOS7ftnqqoTq2p5VS1fsmTJqGuUJDXMbEmaxDAb\n6jXALj3Dy9pxvV4KnF2NlcB1wN5DrEmSNDkzW5JmaZgN9UXAnkn2aC9aOQI4p2+eG4GnAiR5ILAX\ncO0Qa5IkTc7MlqRZGtqjx6tqfZJXAZ8BFgGnVNWVSY5up58AvAU4NckVQIBjqurWYdUkSZqcmS1J\nsze0hhqgqs4Dzusbd0LP65uAZwyzBknSYMxsSZqdoTbUWjiWLtuVtWtWTT+jJEnSZsaGWnNi7ZpV\n7HbMuWOt4YbjDx3r9iVJ0sI07tvmSZIkSZs0G2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtq\nSZIkqQMbakmSJKkDG2pJkiSpAx/sIknSiPhUWWnzZEMtSdKIzIenyoJPlpXmmqd8SJIkSR3YUEuS\nJEkd2FBLkiRJHdhQS5IkSR3YUEuSJEkd2FBLkiRJHdhQS5IkSR3YUEuSJEkd2FBLkiRJHdhQS5Ik\nSR346PFN3NJlu7J2zapxlyFJkrRg2VBv4tauWcVux5w77jK44fhDx12CJEnSWHjKhyRJktSBDbUk\nSZLUgQ21JEmS1IENtSRJktSBDbUkSZLUgQ21JEmS1IENtSRJktSBDbUkSZLUgQ21JEmS1IENtSRJ\nktTBUBvqJAcluSbJyiTHTjHPgUkuS3Jlki8Nsx5J0tTMbEmanYEa6iRnJ3lWkoEb8CSLgPcCBwP7\nAEcm2advnh2A9wGHVdXDgT8YuHJJ0qTMbEkarUHD9n3A84H/SvKOJHsNsMx+wMqquraq7gTOAg7v\nm+f5wNlVdSNAVd0yYD2SpKmZ2ZI0QgM11FX1uap6AfAY4Hrgc0m+nuSlSe49xWI7A6t6hle343o9\nDLhvkguSXJzkxZOtKMlRSVYkWbFu3bpBSpakBcvMlqTRmsnhwPsDLwFeDlwK/BNNWH+2w/YXA48F\nngU8E3hDkof1z1RVJ1bV8qpavmTJkg6bk6SFwcyWpNFZPMhMST4B7AV8EHh2Vd3cTvpwkhVTLLYG\n2KVneFk7rtdq4AdV9VPgp0m+DDwK+N6A9UuS+pjZkjRag+6h/req2qeqjpsI5iS/AVBVy6dY5iJg\nzyR7JNkCOAI4p2+efweemGRxkq2B/YGrZ/wuJEm9zGxJGqFBG+q3TjLuGxtboKrWA68CPkMTuB+p\nqiuTHJ3k6Haeq4H/BC4HvgWcVFXfGbR4SdKkzGxJGqGNnvKR5EE0F6VsleTRQNpJ2wNbT7fyqjoP\nOK9v3Al9w+8E3jmDmiVJkzCzJWk8pjuH+pk0F7UsA/6xZ/ztwF8PqSZJ0uyY2ZI0BhttqKvqNOC0\nJM+tqo+PqCZJ0iyY2ZI0HtOd8vHCqjod2D3Jn/dPr6p/nGQxSdIYmNmSNB7TnfKxTfvvtsMuRJLU\nmZk9haXLdmXtmlXTzyhJszDdKR//2v775tGUI0maLTN7amvXrGK3Y84ddxnccPyh4y5B0hBMd8rH\nuzc2var+dG7LkSTNlpktSeMx3SkfF4+kCknSXDCzJWkMBrnLhyRpE2BmS9J4THfKx7uq6rVJ/gOo\n/ulVddjQKpMkzYiZLUnjMd0pHx9s//37YRciSerMzJakMZjulI+L23+/lGQLYG+avR7XVNWdI6hP\nkjQgM1uSxmO6PdQAJHkWcALwfSDAHkn+pKo+PcziJEkzZ2ZL0mgN1FAD/wA8papWAiR5CPApwHCW\npPnHzJakERq0ob59Iphb1wK3D6GeTYZP3ZI0j5nZkjRC093l4/fblyuSnAd8hOZ8vD8ALhpybfOa\nT92SNN+Y2ZI0HtPtoX52z+v/Bp7cvl4HbDWUiiRJs2VmS9IYTHeXj5eOqhBJUjdmtiSNx6B3+dgS\neBnwcGDLifFV9cdDqkuSNEtmtiSN1r0GnO+DwIOAZwJfApbhBS6SNF+Z2ZI0QoM21A+tqjcAP62q\n04BnAfsPryxJUgdmtiSN0KAN9V3tvz9K8gjgPsADhlOSJKkjM1uSRmjQ+1CfmOS+wBuAc4Bt29eS\npPnHzJakERqooa6qk9qXXwIePLxyJEldmdmSNFoDnfKR5P5J3pPkkiQXJ3lXkvsPuzhJ0syZ2ZI0\nWoOeQ30WcAvwXOB5wK3Ah4dVlCSpEzNbkkZo0HOol1bVW3qG35rkfw2jIElSZ2a2JI3QoHuoz09y\nRJJ7tV9/CHxmmIVJkmbNzJakEdroHuoktwMFBHgtcHo76V7AHcBfDrU6SdLAzGxJGo+NNtRVtd2o\nCpEkdWNmS9J4DHoONUkOAw5oBy+oqnOHU5IkqSszW5JGZ9Db5r0DeA1wVfv1miTHDbMwSdLsmNmS\nNFqD7qE+BNi3qn4JkOQ04FLg9cMqTJI0a2a2JI3QoHf5ANih5/V95roQSdKcMrMlaUQG3UN9HHBp\nki/SXD1+AHDs0KqSJHVhZkvSCE3bUCcJ8FXg8cDj2tHHVNXaYRYmSZo5M1uSRm/aUz6qqoDzqurm\nqjqn/RoomJMclOSaJCuTTLl3JMnjkqxP8rwZ1C5J6mNmS9LoDXoO9SVJHjf9bPdIsgh4L3AwsA9w\nZJJ9ppjveOD8maxfkjQlM1uSRmjQhnp/4MIk309yeZIrklw+zTL7ASur6tqquhM4Czh8kvleDXwc\nuGXgqiVJG2NmS9IIDXpR4jNnse6dgVU9w6tpQv5XkuwM/B7wFO451+/XJDkKOApg1113nUUpkrSg\nmNmSNEIbbaiTbAkcDTwUuAI4uarWz+H230Vzscwvm+toJldVJwInAixfvrzmcPuStNkwsyVpPKbb\nQ30acBfwFe45r+41A657DbBLz/Cydlyv5cBZbTDvCBySZH1VfXLAbUiS7mFmS9IYTNdQ71NVjwRI\ncjLwrRms+yJgzyR70ITyEcDze2eoqj0mXic5FTjXYJakWTOzJWkMpmuo75p4UVXrN3aIr187/6uA\nzwCLgFOq6sokR7fTT5hFvZKkqZnZkjQG0zXUj0ryk/Z1gK3a4dDc7nT7jS1cVecB5/WNmzSUq+ol\nA1UsSZqKmS1JY7DRhrqqFo2qEElSN2a2JI3HoPehliRJkjQJG2pJkiSpAxtqSZIkqQMbakmSJKkD\nG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqQMb\nakmSJKkDG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtq\nSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqQMbakmSJKkDG2pJ\nkiSpAxtqSZIkqQMbakmSJKkDG2pJkiSpAxtqSZIkqYOhNtRJDkpyTZKVSY6dZPoLklye5IokX0/y\nqGHWI0mampktSbMztIY6ySLgvcDBwD7AkUn26ZvtOuDJVfVI4C3AicOqR5I0NTNbkmZvmHuo9wNW\nVtW1VXUncBZweO8MVfX1qvphO3ghsGyI9UiSpmZmS9IsDbOh3hlY1TO8uh03lZcBn55sQpKjkqxI\nsmLdunVzWKIkqWVmS9IszYuLEpM8hSacj5lselWdWFXLq2r5kiVLRlucJGkDZrYkbWjxENe9Btil\nZ3hZO24DSX4LOAk4uKp+MMR6JElTM7MlaZaGuYf6ImDPJHsk2QI4Ajind4YkuwJnAy+qqu8NsRZJ\n0saZ2ZI0S0PbQ11V65O8CvgMsAg4paquTHJ0O/0E4G+B+wPvSwKwvqqWD6smSdLkzGxJmr1hnvJB\nVZ0HnNc37oSe1y8HXj7MGiRJgzGzJWl25sVFiZIkSdKmyoZakiRJ6sCGWpIkSerAhlqSJEnqwIZa\nkiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqS\nJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIk\nSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ\n6sCGWpIkSerAhlqSJEnqwIZakiRJ6sCGWpIkSerAhlqSJEnqwIZakiRJ6mCoDXWSg5Jck2RlkmMn\nmZ4k726nX57kMcOsR5I0NTNbkmZnaA11kkXAe4GDgX2AI5Ps0zfbwcCe7ddRwL8Mqx5J0tTMbEma\nvWHuod4PWFlV11bVncBZwOF98xwOfKAaFwI7JFk6xJokSZMzsyVpllJVw1lx8jzgoKp6eTv8ImD/\nqnpVzzznAu+oqq+2w58HjqmqFX3rOopmbwjAXsA1HcvbEbi14zrmgnVsaD7UMR9qAOvotznVsVtV\nLZmLYuaSmT0Q69iQdcyvGsA6+o0ssxd33MhIVNWJwIlztb4kK6pq+Vytzzo2nzrmQw3WYR2bOjPb\nOhZaHfOhBusYbx3DPOVjDbBLz/CydtxM55EkDZ+ZLUmzNMyG+iJgzyR7JNkCOAI4p2+ec4AXt1eO\nPx74cVXdPMSaJEmTM7MlaZaGdspHVa1P8irgM8Ai4JSqujLJ0e30E4DzgEOAlcDPgJcOq54+c3Yo\nsiPr2NB8qGM+1ADW0c86hszMHoh1bMg67jEfagDr6DeyOoZ2UaIkSZK0EPikREmSJKkDG2pJkiSp\ngwXXUE/3aN0R1XBKkluSfGcc229r2CXJF5NcleTKJK8ZUx1bJvlWkm+3dbx5HHX01LMoyaXt/XbH\nVcP1Sa5IclmSFdMvMbQ6dkjysSTfTXJ1kt8e8fb3aj+Dia+fJHntKGvoqeXP2p/P7yQ5M8mW46hj\nITKzN6hj7LltZk9ag5l9Tw3zIrfHkdkL6hzqNI/W/R7wdGA1zVXtR1bVVSOu4wDgDponjj1ilNvu\nqWEpsLSqLkmyHXAx8JwxfBYBtqmqO5LcG/gq8Jr2KWwjl+TPgeXA9lV16JhquB5YXlVjvSl+ktOA\nr1TVSe1dH7auqh+NqZZFNLdn27+qbhjxtnem+bncp6p+nuQjwHlVdeoo61iIzOxfq2PsuW1mT1rD\n9ZjZk9UzltweV2YvtD3Ugzxad+iq6svAbaPebl8NN1fVJe3r24GrgZ3HUEdV1R3t4L3br7H8lZdk\nGfAs4KRxbH8+SXIf4ADgZICqunOcwQw8Ffj+qJvpHouBrZIsBrYGbhpTHQuNmb1hHWPPbTN7fpqH\nmQ3jze2RZ/ZCa6h3Blb1DK9mDE3kfJNkd+DRwDfHtP1FSS4DbgE+W1VjqQN4F/BXwC/HtP0JBXwu\nycVpHuE8DnsA64D3t4dTT0qyzZhqgeaeyGeOY8NVtQb4e+BG4Gaaey+fP45aFiAzewrjzG0z+9eY\n2ZMbS26PK7MXWkOtPkm2BT4OvLaqfjKOGqrq7qral+apa/slGfkh1SSHArdU1cWj3vYknth+HgcD\nr2wPN4/aYuAxwL9U1aOBnwLjOn91C+Aw4KNj2v59afaK7gHsBGyT5IXjqEWC8ee2mf1rzOw+48zt\ncWX2QmuofWxuj/b8t48DZ1TV2eOupz089UXgoDFs/gnAYe25cGcBv5vk9DHUMfHXNVV1C/AJmsPe\no7YaWN2z5+ljNGE9DgcDl1TVf49p+08DrquqdVV1F3A28DtjqmWhMbP7zKfcNrMbZvakxpnbY8ns\nhdZQD/Jo3QWhvbDkZODqqvrHMdaxJMkO7eutaC4++u6o66iq11fVsqranebn4gtVNfK9kEm2aS82\noj1c9wxg5HcWqKq1wKoke7WjngqM9EKwHkcyptM9WjcCj0+ydfv/5qk0565q+MzsHvMht83sDZnZ\nUxpnbo8ls4f26PH5aKpH6466jiRnAgcCOyZZDbyxqk4ecRlPAF4EXNGeCwfw11V13ojrWAqc1l4N\nfC/gI1U1ttsfzQMPBD7RZACLgQ9V1X+OqZZXA2e0jcy1jO4x07/S/oJ6OvAno972hKr6ZpKPAZcA\n64FLmT+P1d2smdm/Zj7ktpm9ITO7z7hze1yZvaBumydJkiTNtYV2yockSZI0p2yoJUmSpA5sqCVJ\nkqQObKglSZKkDmyoJUmSpA5sqDUWSf4myZVJLk9yWZL92/GvTbL1LNf5piR/OctllyT5ZvvI1if1\nTbsgyTVtnZe1t+ORpAXDzJY2bkHdh1rzQ5LfBg4FHlNVv0iyI7BFO/m1wOnAz0Zc1lOBK6rq5VNM\nf0FVrZirjbU3m09V/XKu1ilJw2Bmm9mannuoNQ5LgVur6hcAVXVrVd2U5E+BnYAvJvkiQJIjk1yR\n5DtJjp9YQZKDklyS5NtJPt+/gSSvSPLp9kleveN3T/KFdi/L55PsmmRf4O+Aw9u9GVv1r28ySU5N\n8u4kX09ybZLn9Ux7XZKL2u28uWfb1yT5AM2TtHZJ8rIk30vyrST/luSfk2yX5Lo0jxgmyfa9w5I0\nYma2ma1p2FBrHM6nCabvJXlfkicDVNW7gZuAp1TVU5LsBBwP/C6wL/C4JM9JsgT4N+C5VfUo4A96\nV57myWqHAs+pqp/3bfs9wGlV9VvAGcC7q+oy4G+BD1fVvpMsA83TpyYOH76zZ/xS4Int9t7Rbv8Z\nwJ7Afm3dj01yQDv/nsD7qurhwF3AG4DH0zwBbe/2c7gduAB4VrvMEcDZVXXXNJ+rJA2DmW1maxqe\n8qGRq6o7kjwWeBLwFODDSY6tqlP7Zn0ccEFVrQNIcgZwAHA38OWquq5d3209y7wYWEUTzJOF2W8D\nv9++/iDNXo5BTHX48JPtIcCrkjywHfeM9uvSdnhbmlC+Ebihqi5sx+8HfGmi/iQfBR7WTjsJ+Cvg\nkzSPj30XBp6FAAAByUlEQVTFgHVK0pwys81sTc+GWmNRVXfT/EV/QZIrgD8CTp2DVV9Bs4dhGXDd\nHKxvOr/oeZ2ef4+rqn/tnTHJ7sBPB1lpVX2tPdx4ILCoqr7TvVRJmh0ze+PMbHnKh0YuyV5J9uwZ\ntS9wQ/v6dmC79vW3gCcn2THJIuBI4EvAhcABSfZo13e/nnVdCvwJcE57+LHf12kOxwG8APjKHLyl\nfp8B/jjJtm19Oyd5wCTzXUTz/u6bZDHw3L7pHwA+BLx/CDVK0kDM7F8xszUl91BrHLYF3pNkB2A9\nsBI4qp12IvCfSW5qz8k7FvgizR6ET1XVvwMkOQo4O8m9gFuAp0+svKq+muZWTJ9K8vSqurVn268G\n3p/kdcA6mkNzgzgjycR5erdW1dOmmrGqzk/ym8A3kgDcAbyQ5rBn73xrkryd5pfQbcB3gR/3bhN4\nK3DmgDVK0jCY2ZjZ2rhU1bhrkBasJNu25ycuBj4BnFJVn2inPQ84vKpeNNYiJUmAma2puYdaGq83\nJXkasCXNlfSfBEjyHuBg4JAx1iZJ2pCZrUm5h1qSJEnqwIsSJUmSpA5sqCVJkqQObKglSZKkDmyo\nJUmSpA5sqCVJkqQO/j8RhXEomqMKwwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11258ab00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "ts = [0, 5]\n",
    "for i, t in enumerate(ts):\n",
    "    axes[i].bar(np.arange(n), vs[t], align='center', width=1, edgecolor='k')\n",
    "    axes[i].set_xlim(0-0.5, emax+0.5)\n",
    "    axes[i].set_ylim(0, 1)\n",
    "    axes[i].set_xlabel('Stock of Energy')\n",
    "    axes[i].set_ylabel('Probability')\n",
    "    axes[i].set_title('Survival Probability, Period {t}'.format(t=t))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}