{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Bioeconomic Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miranda and Fackler, Applied Computational Economics and Finance, 2002, 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": [
    "s_bar = 8  # Energy carrying capacity\n",
    "n = s_bar + 1  # Number of states, 0, ..., s_bar\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.empty(n)\n",
    "v_term[0], v_term[1:] = 0, 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We follow the state-action pairs formulation approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "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": false
   },
   "outputs": [],
   "source": [
    "# Reward vector\n",
    "R = np.zeros(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "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",
    "    else:\n",
    "        Q[i, np.minimum(s_bar, 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": false
   },
   "outputs": [],
   "source": [
    "# Discount factor\n",
    "beta = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "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": false
   },
   "outputs": [],
   "source": [
    "vs, sigmas = backward_induction(ddp, T, v_term)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAEZCAYAAACHP6MmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYJXV97/H3B4hL2BS3ERBwJxoFSUSMBjtujBLFuNwA\nUaPGyM2j0SQ3RjQxjLkmYHJv4pbEi0EMiKJxRWMUNUzcFQXFBQQkjmwZF0QRlwB+7x9Vw5xTc7qn\nprtPn5rp9+t5+pmz1KnzOdXT3/72r35VlapCkiRJ0sJ2mnUASZIkaXtg4yxJkiT1YOMsSZIk9WDj\nLEmSJPVg4yxJkiT1YOMsSZIk9WDjrAUl+cckf7oM6zk1yV8sR6Ye7/WzJHdb5Gv/M8nD53nuoUku\nnLRskhcnOXlxiVdOkrsk+UGSLOK1D0ty+TRySVoe1uyx56zZ1uxlZ+O8HWqLwSeSXJvkO0k+luSX\npvFeVfV7VfWX01j3Jkl+O8mNbXG4Nsl5SY5cwiqncnLyqvp4Vf3CPM+dWFXPAUiyf/uLYFl+vtpf\nYD9tt893kpyd5N6LWVdVXV5Ve9TiT+A+7+vaz/3vSa5P8tUkj1jke0g7FGv2Vlmz5zHlmv2NJD9q\nc/4gyQcW+R6rio3zdibJ7sB7gVcBtwX2AV4G/HSR69vmv2Kn5JNtcbgN8AbgbUn27C6UZOce65r1\nZwpNsVrOHK+oqj2AfYFvAaduc6h+224p3gJ8HtgL+DPg7UluN+X3lAbNmm3NZrg1u4Aj2+/jHlW1\ndsrvt0Owcd7+3AuoqnpbNX5aVR+uqi8DJDkhyembFu7+JZ3knCQvT/LxJNcDL0xy7ugbJPnDJO9u\nb9+8u64dRXzsyHI7J/lWkoPb+29LcnWS7yVZn+Q+i/yMbwBuDdx9066mJH+S5Or2OZL8bpJL2r/m\n353kzp11HJnk622+vx7JfLckH2lf960kb0qyR+e1hyb5SpLvJjklyS3a186726vd7qe1d/+j/ffa\n9q/4w9t13Xdk+Tu0I7Pb1FhW1U+ANwO/2K4nSY5PcmmSbyc5M8lt2uc2fe+flWQD8JEJ/x/unOQ9\nbb6Lkzx7JOOtkrwxyTVJvgw8cL5cSe4JPABY1/6ffCdwAfCkbfl80g7Imo01m4HV7NFNsS2fRzbO\n26OLgZvaH461m37gOrq7Zrr3nwo8G9gdeB1wryR3H3n+GOCMCet9C3DsyP21wLer6gvt/fcDdwfu\nCJw3zzoWlGQX4HeB64BL2ofXALcB9gOek2aO2l8BTwbuDHwTOLOzqicAh7RfRyV51qa3aF+7BvgF\nmtGAdZ3XHgs8qv0s96YZPd2kz+6yw9t/N/0V/1GabffUkWWOAT5cVd/tsb6bJdkN+C2a7QvwfODx\nwK8CewPfA/5hQp4DgSMmfIa30my/NcBTgL9KMtc+tw64a/t1BPDbC0S7L3BZVV0/8tgX28el1cya\nbc0eYs3e5IwkG5N8IMn9e3+wVczGeTtTVdcBDwV+BpwMfKv96/MO27CaN1bVRVX1s6r6AfAemqKw\naeTw3jS7FrveDDw+ya3a+8fQFJdN2d5YVT+qqhuAvwAOSrObso8HJ7kGuAr4TeAJ7WcFuAk4oapu\nqKqf0hTJU6rqi+17vbh9/X4j6zupqr5fVVcAr9z0+arq61X1kaq6sS2Afwc8rJPlNVV1VVVdC/zl\nptcuwuhf8qcx/gvsacDp9PfCdvtcDOwKPKN9/DjgT6vq6pHt/uRsnqtXNNvux+222xwuuQvwYOBF\n7bb9IvBPwNPbRZ4CvLzdjlcCr14g327A9zuP/YDmF720almzrdkMs2bTfr4DgP2B9cAHJ4zmq8PG\neTtUVV+rqmdV1X40u3/2pik0fXV3Xb2FzYXmWODd7e6l7vt+Hfgq8Lgkt6b5q/nNAEl2SnJSu/vp\nWuA/aQrA7Xtm+lRV7VVVd6yqX6mqc0ae+3ZbYDbZG9gwkut64Ls0cwc3uWLk9ob2NSS5Y5K3JLmi\nzfmmCRknvnYpquqzwPXtrsN704yMnLUNq/ibdvvsXVVPqKpvtI/vD7yr3TV3Dc335wbgTiOvvYLJ\n7gxcU1U/GnlsA5u3495suS3m80OgW3D3pBmFklY1a7Y1e4A1m6r6VDt16CdVdRJwLc1IuBZg47yd\nq6qLgTfSzp8Crgd+fmSR7jwy2HLX1YeAOyQ5CDiatrDO40yaQn0U8JWquqx9/FjgccDDqzlY5ACa\nv96XY/5UN+9VNMUHgCS7ArdjvGDcZeT2/u1rAE6kGfm5b5vzqRMyzvfaxebd5J9pRi2eBry9qv57\nG9c7yTeBx7QFeq+qum1V7VpVV/fIcxWwV7v9NtkPuLK9fTVbbov5fAW4W2ddB7WPS2pZs63ZDKNm\nT7LcB0jukGyctzNJ7p3kj5Ls096/C83Iw6faRb4AHJ7m3I97AsdvbZ1VdSPwL8Df0Bz1/aEFFj8T\neDTwe4wX691pjhL/XvtDfSJTOsUQzWjLM5PcP8ktaea/fbqqRkdlXpjkNu32eT6b59PtRjM6el27\nDV84Yf3PTbJPkr2Al7DlXLyt+TZNob975/EzgN+gme922ugT7cEfh7Pt/h/NHLf92vXcIcnjR1c9\n4TUBaHeJfhI4Mckt2/ltv8Pm3ZFvA17cbsd9gefNF6KqLqH5v3dCu64n0jQG71jEZ5J2GNZswJo9\nahA1u/3/9itJfq5d1wtp/pj5xCI+06pi47z9uQ54EPCZJNfR/BBdAPwxQFV9mObggQuAc9ly3tt8\nhfEtwCOAt1XVz+Zbvqr+i6bgH9a+zyan0fwlfSXw5TbXVFTVR4CXAu9s3++uNKMuo5nfQ3NqtPNo\ntsEb2udeBvwSzS6p97JlY1c0v1zOBi6lOdhlvnOiTtyWVfXj9jWfaHfHHdo+fkWbp6rq45uWb39R\n/AD40ra8T+tVNJ/17CTfp9nuh27ltaOPHUOz/a6i2RYvHdnl+jKa7+l/Ah+g84tjgqNpjuL+Hs3n\nf1Jt44E00g7Imm3NHjWUmr078I/ANTQj/48G1lbV9xZ4jYDUos+p3WPlySnArwMbq2ri0ZpJXg08\nhmZ31TNq89G+0g6n/Zm4sqr+fOSx3wLuU1VLvtqXtFTWbWkza7a6pt04P5RmF8tpkwpwkscAz6uq\nI5M8CHhVVR02tUDSDCU5gGb04gFVteBBG9KsWLelhjVbk0x1qka7a2OhYf+jaHclVNVngD2T3GmB\n5aXtUpoLElwA/LUFWENm3Zas2ZrfrOc478P4aXauZPz0NNIOoar+vJoT65806yzSElm3tcOzZms+\ns26cJUmSpO3CLjN+/ysZP+fgvmw+H+GYJNObjC1JK6CqdoRzpPaq29ZsSdu7STV7JUacFzqh+lm0\nl4pMchhwbVVtnG9FVbWkrxNOOGHJ69iRvtwebgu3x8ptj+3MstTtIWz3HenL7bFjbo/2p2WJXycs\nwzqWI8dyfA1pe2xpqiPOSd4MzAG3S/JNmk9yC5pzIp5cVe9P8tgkl9Kc1uiZ08wjafVZs+YANm5c\n+rE9L3vZy5YhzfBZtyVpflNtnKvq2B7LzHtlG0laqqZpXuqI77r2aym2j1ka1m1Jmt+s5zivqLm5\nuVlHGBS3x2Zui3HLsT2Wa6R3GOZmHWBV8udynNtjnNtj1NysAwzM3NTWPNULoCynJLW9ZJUESVj6\nSO9yGE6O2jEODuzFmi31Y63sGk6OSTV7VY04S5IkbbJj7RnTSnDEWdJUOIrS5YizNDTDqFNDyADm\n6Jpcs70AiiRJktSDjbMkSZLUg3OcpR2Mc/YkSZoO5zhLO5hhzNmDIc1TG0oO5zhLwzKMejmEDGCO\nLs+qIUmSBsA9Y9peOeIs7WCGMYICQxo1GEoOR5ylhnVqaBnAHF2eVUOSJElaNBtnSZIkqQcbZ0mS\nJKkHDw6UlokHu0iStGPz4EBpmXiwS5c5xnlwoLSJ9XJoGcAcXR4cKEmSJC2ajbMkSZLUg42zJEmS\n1IONsyRJktSDjbMkSZLUg42zJEmS1IONsyRJktSDjbMkSZLUg42zJEmS1IONsyRJktSDjbMkSZLU\ng42zJEmS1IONsyRJktSDjbMkSZLUwy6zDiAt1Zo1B7Bx44ZZx5AkSTu4VNWsM/SSpLaXrFpZSYAh\n/N8wxzhzjAtVlVmnWCnWbC3Euj20DGCOrsk126kakiRJUg9O1ZAkaZVwapu0NE7V0HbPXX5d5hg3\nnBxO1dCsWS+7hpBjCBnAHF1O1ZAkSZIWzcZZkiRJ6mHqjXOStUkuSnJxkhdNeH6PJGcl+UKSLyV5\nxrQzSZIms2ZL0vymOsc5yU7AxcAjgKuAc4Gjq+qikWVeDOxRVS9Ocnvga8CdqurGzrqcL6eJnLPX\nZY5xw8kx9DnO1uwdn/Wyawg5hpABzNE1mznOhwKXVNWGqroBOBM4qrNMAbu3t3cHvtstwJKkFWHN\nlqQFTLtx3ge4fOT+Fe1jo14L3CfJVcAXgRdMOZMkaTJrtiQtYAjncT4COL+qHp7k7sCHkty/qn7Y\nXXDdunU3356bm2Nubm7FQkrStlnffu1wrNmSdkDr6VOzpz3H+TBgXVWtbe8fD1RVvWJkmfcBJ1bV\nJ9r7HwFeVFWf66zL+XKayDl7XeYYN5wc28EcZ2v2Ds562TWEHEPIAOboms0c53OBeyTZP8ktgKOB\nszrLbAAeCZDkTsC9gMumnEuStCVrtiQtYKpTNarqpiTPA86madJPqaoLkxzXPF0nAy8H3pjkgvZl\nf1JV10wzlyRpS9ZsSVqYl9zWds9dj13mGDecHEOfqrGcrNnDZL3sGkKOIWQAc3R5yW1JkiRp0Wyc\nJUmSpB6GcDo6bcfWrDmAjRs3zDqGJEnS1DnHWUsyjPlyQ8gA5ugyxzjnOGv2hlGzYUg/l7PPMYQM\nYI4u5zhLkiRJi2bjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmS\nJPXglQMlSZoyr7Iq7Ri8cqCWZBhXoRpCBjBHlznGeeXA1WwYtRKG9PNgjiFlAHN0eeVASZIkadFs\nnCVJkqQebJwlSZKkHmycJUmSpB5snCVJkqQebJwlSZKkHmycJUmSpB5snCVJkqQebJwlSZKkHmyc\nJUmSpB52mXUALc6aNQewceOGWceQJElaNVI1hOuBb12S2l6yroRkONdyn32OIWQAc3SZY1yoqsw6\nxUqxZo+zZneZY1gZwBxdk2u2UzUkSZKkHmycJUmSpB5snCVJkqQebJwlSZKkHmycJUmSpB5snCVJ\nkqQebJwlSZKkHmycJUmSpB5snCVJkqQept44J1mb5KIkFyd50TzLzCU5P8mXk5wz7UySpMms2ZI0\nv6lecjvJTsDFwCOAq4BzgaOr6qKRZfYEPgk8uqquTHL7qvrOhHV5+dYRXr51aBnAHF3mGDf8S25b\ns6fHmt1ljmFlAHN0zeaS24cCl1TVhqq6ATgTOKqzzLHAO6rqSoBJBViStCKs2ZK0gF6Nc5J3Jjmy\nHY3YFvsAl4/cv6J9bNS9gL2SnJPk3CRP28b3kCSNsGZL0nT0Lar/QDPKcEmSk5Lcexkz7AIcAjwG\nWAu8NMk9lnH9krTaWLMlaQp26bNQVX0Y+HA7t+2Y9vblwOuBN7W79Ca5Ethv5P6+7WOjrgC+U1U/\nAX6S5KPAQcCl3ZWtW7fu5ttzc3PMzc31iS9JM7C+/Vp51mxJ2lbr6VOzex8cmOR2wFOBp9EcNHIG\n8FDgflU1N89rdga+RnOgydXAZ4FjqurCkWUOBF5DM3JxS+AzwG9W1Vc76/JAkxEeaDK0DGCOLnOM\nW9mDA63Zw2LN7jLHsDKAObom1+xeI85J3gXcGzgdeFxVXd0+9dYkn5vvdVV1U5LnAWfTTAs5paou\nTHJc83SdXFUXJfkgcAFwE3BytwBLkvqzZkvSdPQacU7y2Kp6f+exW1bVT6eWbMsMjl6McPRiaBnA\nHF3mGLdyI87W7OGxZneZY1gZwBxdSzsd3csnPPappQWSJE2JNVuSpmDBqRpJ1tCciujWSR5A82cA\nwB7Az085myRpG1izJWm6tjbH+QjgGTRHVv/tyOPXAS+ZUiZJ0uJYsyVpivrOcX5SVb1jBfIslMH5\nciOcLze0DGCOLnOMW9E5ztbsgbFmd5ljWBnAHF2Ta/aCjXOSp1bVm5L8LyZ8iqr62wkvmwqL8DiL\n8NAygDm6zDFu+o2zNXu4rNld5hhWBjBH1+JOR7dr++9uyx9IkrTMrNkda9YcwMaNG2YdQ9IOovcF\nUGbN0Ytxjl4MLQOYo8sc41b2AiizNpSaba3sMse4IeQYQgYwR9ciRpyTvHqh56vq+UuNJUlaHtZs\nSZqurU3V+PyKpJAkLQdrtiRNkVM1tlPufhxaBjBHlznGOVVjRjkYyvffHKPMMawMYI6uxU3VeGVV\n/UGS9zL5CO3HL2NCSdISWLMlabq2NlXj9Pbf/zPtIJKkJbNmS9IU9Z6qkeQWwIE0oxhfq6r/nmaw\nCe8/iN1+Q+Hux6FlAHN0mWPcyk7VsGbfnIOhfP/NMcocw8oA5uha3Hmcm5cmRwKvA75O84numuS4\nqvq35Q0pSVoqa7YkTUffS25fBPx6VV3a3r878K9VdeCU841mGMToxVA4ijK0DGCOLnOMW9FLbluz\nN+dgKN9/c4wyx7AygDm6ljDiDFy3qQC3LgOuW5Zc2xmvQiVpO2DNlqQp2NpZNZ7Y3vxckvcDb6P5\nM+ApwLlTzjZITdM8jL+EJGmUNVuSpmtrI86PG7m9EXhYe/vbwK2nkkiStFjWbEmaIi+Asu05GM6I\nszmGkwHM0WWOcV4AZUY5GMr33xyjzDGsDGCOrqWdVeNWwO8A9wVutenxqnrWsuWTJC0La7YkTcdO\nPZc7HVgDHAH8B7AvHmgiSUNlzZakKeh7Orrzq+oBSS6oqvsn+TngY1V12PQj3pzB3X5jzDGsDGCO\nLnOMW9HT0VmzN+dgKN9/c4wyx7AygDm6JtfsviPON7T/XpvkF4E9gTsuVzRJ0rKyZkvSFPQ9j/PJ\nSW4LvBQ4C9itvS1JGh5rtiRNgWfV2PYcDGUXgjmGlAHM0WWOcZ5VY0Y5GMr33xyjzDGsDGCOriVM\n1UhyuySvSXJeks8neWWS2y1/SEnSUlmzJWk6+s5xPhP4FvAk4MnAd4C3TiuUJGlJrNmSNAV9z6rx\n5ar6xc5jX6qq+00t2ZYZ3O03xhzDygDm6DLHuBU9q4Y1e3MOhvL9N8cocwwrA5ija2ln1Tg7ydFJ\ndmq//gfwweUNKElaJtZsSZqCBUeck1xH0/YH2BX4WfvUTsAPq2qPqSfcnMXRizHmGFYGMEeXOcZN\nf8TZmj0xB0P5/ptjlDmGlQHM0bWIS25X1e7TCyRJWk7WbEmarr7ncSbJ44HD27vrq+p904kkSVoq\na7YkLb++p6M7CXgB8NX26wVJTpxmMEnS4lizJWk6+p5V4wLg4Kr6WXt/Z+D8qrr/lPONZnC+3Bhz\nDCsDmKPLHONW9Kwa1uzNORjK998co8wxrAxgjq6lnVUD4DYjt/dceiBJ0hRZsyVpmfWd43wicH6S\nc2j+FDgcOH5qqSRJS2HNlqQp2OqIc5r9XB8HDgPeCbwDeHBV9boKVZK1SS5KcnGSFy2w3AOT3JDk\niT2zS5I6rNmSND195zgv6opTSXYCLgYeAVwFnAscXVUXTVjuQ8CPgTdU1TsnrMv5cmPMMawMYI4u\nc4xb0TnO1uzNORjK998co8wxrAxgjq6lzXE+L8kDF/GuhwKXVNWGqroBOBM4asJyvw+8HfjWIt5D\nkjTOmi1JU9B3jvODgKcm+QZwPe2fAz2O0N4HuHzk/hU0hflmSfYGnlBVv5Zk7DlJ0qJYsyVpCvo2\nzkdMMcMrgdF5dCuyK1OSdmDWbEmaggUb5yS3Av4ncA/gS8ApVXXjNqz/SmC/kfv7to+N+mXgzPaA\nltsDj0lyQ1Wd1V3ZunXrbr49NzfH3NzcNkSRpJW0vv1aOdZsSVqs9fSp2QseHJjkrcANwMeAxwAb\nquoFfSO0J93/Gs2BJlcDnwWOqaoL51n+VOC9HmjShzmGlQHM0WWOcdM/ONCavSVrdpc5xg0hxxAy\ngDm6JtfsrU3VuM+mI7OTnEJTRHurqpuSPA84m+ZAxFOq6sIkxzVP18ndl2zL+iVJY6zZkjRFWxtx\nPq+qDpnv/kpy9KLLHMPKAOboMse4FRlxtmZvmYOhfP/NMcocw8oA5uiaXLO31jjfRHNEdrMGuDXw\nIzYfob3HFJLOl8UiPMYcw8oA5ugyx7gVaZyt2VvmYCjff3OMMsewMoA5uhYxVaOqdp5eIEnScrJm\nS9J09b0AiiRJkrSq2ThLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhL\nkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuS\nJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5Ik\nST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJPdg4S5IkST3YOEuSJEk92DhLkiRJ\nPdg4S5IkST3YOEuSJEk9TL1xTrI2yUVJLk7yognPH5vki+3Xx5Pcb9qZJEmTWbMlaX5TbZyT7AS8\nFjgCuC9wTJIDO4tdBhxeVQcBLwdeP81MkqTJrNmStLBpjzgfClxSVRuq6gbgTOCo0QWq6tNV9f32\n7qeBfaacSZI0mTVbkhYw7cZ5H+DykftXsHCRfTbwb1NNJEmajzVbkhawy6wDbJLk14BnAg+db5l1\n69bdfHtubo65ubmp55KkxVnffu2YrNmSdizr6VOzU1VTi5DkMGBdVa1t7x8PVFW9orPc/YF3AGur\n6uvzrKummbWvJMDsc4A5hpUBzNFljnGhqjLrFAuxZk+TOcaZY1gZwBxdk2v2tKdqnAvcI8n+SW4B\nHA2cNRYr2Y+mAD9tvgIsSVoR1mxJWsBUp2pU1U1JngecTdOkn1JVFyY5rnm6TgZeCuwF/EOaoYEb\nqurQaeaSJG3Jmi1JC5vqVI3l5G6/LnMMKwOYo8sc44Y/VWM5WbO7zDHOHMPKAOboms1UDUmSJGmH\nYOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg\n4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDj\nLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMs\nSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJkiT1YOMsSZIk9WDjLEmSJPVg4yxJ\nkiT1YOMsSZIk9TD1xjnJ2iQXJbk4yYvmWebVSS5J8oUkB087kyRpMmu2JM1vqo1zkp2A1wJHAPcF\njklyYGeZxwB3r6p7AscBr5tWnvXr109r1dup9bMOMCDrZx1gYNbPOsDArJ91gBVhzR669bMOMDDr\nZx1gQNbPOsDArJ/amqc94nwocElVbaiqG4AzgaM6yxwFnAZQVZ8B9kxyp2mEsQh3rZ91gAFZP+sA\nA7N+1gEGZv2sA6wUa/agrZ91gIFZP+sAA7J+1gEGZv3U1jztxnkf4PKR+1e0jy20zJUTlpEkTZ81\nW5IW4MGBkiRJUg+pqumtPDkMWFdVa9v7xwNVVa8YWeZ1wDlV9db2/kXAw6pqY2dd0wsqSSugqjLr\nDAuxZkvSZpNq9i5Tfs9zgXsk2R+4GjgaOKazzFnAc4G3tkX72m4BhuH/wpGkHYA1W5IWMNXGuapu\nSvI84GyaaSGnVNWFSY5rnq6Tq+r9SR6b5FLgeuCZ08wkSZrMmi1JC5vqVA1JkiRpR7FqDg7sc1L/\n1SDJvkn+PclXknwpyfNnnWkIkuyU5LwkZ806y6wl2TPJvyS5sP1/8qBZZ5qVJH+Y5MtJLkhyRpJb\nzDrTamHN3sy6vSVr9mbW7M1Womavisa5z0n9V5EbgT+qqvsCDwaeu4q3xagXAF+ddYiBeBXw/qr6\nBeAg4MIZ55mJJHsDvw8cUlX3p5nadvRsU60O1uwtWLe3ZM3ezJrNytXsVdE40++k/qtCVf1XVX2h\nvf1Dmh+wVX0O1iT7Ao8F/mnWWWYtyR7Ar1bVqQBVdWNV/WDGsWZpZ2DXJLsAPw9cNeM8q4U1e4R1\ne5w1ezNr9hamXrNXS+Pc56T+q06SA4CDgc/MNsnM/R3wQsAJ/3BX4DtJTm13g56c5NazDjULVXUV\n8H+Bb9Jc5OPaqvrwbFOtGtbseVi3AWv2KGt2a6Vq9mppnNWRZDfg7cAL2hGMVSnJkcDGdjQn7ddq\ntgtwCPD3VXUI8CPg+NlGmo0kt6EZ5dwf2BvYLcmxs02l1cy6bc2ewJrdWqmavVoa5yuB/Ubu79s+\ntiq1uzDeDpxeVe+ZdZ4Zewjw+CSXAW8Bfi3JaTPONEtXAJdX1efa+2+nKcqr0SOBy6rqmqq6CXgn\n8CszzrRaWLM7rNs3s2aPs2ZvtiI1e7U0zjef1L89wvJompP4r1ZvAL5aVa+adZBZq6qXVNV+VXU3\nmv8X/15VT591rllpL2RxeZJ7tQ89gtV7AM43gcOS3CpJaLbFqjzoZgas2VuybmPN7rJmj1mRmj3t\nKwcOwnwn9Z9xrJlI8hDgt4AvJTmfZo7YS6rqA7NNpgF5PnBGkp8DLmOVXuCiqj6b5O3A+cAN7b8n\nzzbV6mDNHmfd1lZYs1m5mu0FUCRJkqQeVstUDUmSJGlJbJwlSZKkHmycJUmSpB5snCVJkqQebJwl\nSZKkHmycJUmSpB5snLWikvxpki8n+WKS85I8sH38BUlutch1npDkjxb52tsn+XSSz7fnSh197pwk\nF7U5z0/ytsW8hyRtr6zZ0rhVcQEUDUOSw4DHAgdX1Y1J9gJu0T79B8DpwE9WONYjgQuq6jnzPH9M\nVZ2/nG9K7BaHAAADmUlEQVSYJOUJ1CUNnDW7Yc3WKEectZLuDHynqm4EaK8n/19Jfh/YGzgnyUcA\nkhyT5IL266RNK0iyth1p+EKSD3XfIMnvJvnXJLfsPL5/ko+0oyYfSrJvkoOAVwBHtSMUt+yujwk/\nI0lOTfKqJJ9IcmmSJ44898dJPtvmO2HkvS9K8s9JvgTsm+R3knytHTk5Ocmrk+yW5LIkO7ev2330\nviStMGu2NVsdNs5aSWcD+7UF6e+THA5QVa8BrgTmquoRSe4MnATMAQcDD0zy+CS3p7l85m9U1cHA\nU0bWnSTPpRkdeUJV/bTz3q8BTq2qg4A3A6+pqi8Cfw68taoOmfAagDe1Bfq8JK8YeXxNVT0EeBxN\nISfJo4B7VtWhwAOAX07y0Hb5ewKvrar7ATcCfwYcCjwEOLDdDj8EzgGObF9zNPCOqrpp65tWkpad\nNduarQ6namjFVNX1SQ4BfhV4OHBmkuOr6jQg7RfAA4FzquoagCRnAIcDPwP+o6q+2a7v2pHVPx34\nJk0BnlS0Hgz8Rnv7dNrC2cOx8+z2e3eb4cIkd2wfezTwqCTntZ9lV5rieznwjao6t13uUGB9VX2/\n/Xz/0i4HcArwQuAs4JnAs3vmlKRlZc22ZmtLNs5aUe08sY8CH213gT0dOG3Copnw2EKPX0Az0nEX\n4BuT3nrbkm71/X46YZkAJ1bV68dWkOwPXN9nvVX1ySQHJHkYsFNVfXURmSVpWVizF16vNXv1caqG\nVkySeyW5x8hDBwMb2ts/APZob38WODzJXu1csWOA9cCngV9tixpJbjuyrvOB44Cz2t2GXZ9s1wPw\nVOBjfWNvwzIfBJ6VZNc2395J7jBhPefSfL49k+wCPKmzvtNpdk2+oWdGSVp21uybWbN1M0ectZJ2\nA16TZE+aOWOXApuOjH498IEkV7Zz5l5MU3gB3ldV7wNI8hzgXUkCfAs4YtPK27/8/xh4X5JHbdpt\n2Ho+cGr7/Ldpdqn18aYkP6Ypot+uqkez5UhIte//oSQHAp9q4nEdTcH/2ehrquqqJH9F88vmGuAi\n4Psj6zsD+N/AmT0zStI0WLOxZmtcPMOKtPKS7NrOH9wZeBdwSlW9p33uycDjquq3ZxpSkgRYs7WZ\nI87SbKxL8kjglsDZIwX41cBamiPNJUnDYM0W4IizJEmS1IsHB0qSJEk92DhLkiRJPdg4S5IkST3Y\nOEuSJEk92DhLkiRJPdg4S5IkST38f62CocR5TAMwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119f64438>"
      ]
     },
     "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)\n",
    "    axes[i].set_xlim(0-0.5, s_bar+0.5)\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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}