{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multicriteria approach to Miller-Orr cash management model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Loading dataset of cash flows and visualizing some properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we import the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading the data set and dividing training and test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cf=np.loadtxt(\"datasets/datasets/cash.csv\",delimiter=\";\")\n",
    "f = np.array(cf,dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(800, 200)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = int(0.8*len(f))\n",
    "f_train = f[:n] \n",
    "f_test = f[n:] \n",
    "len(f_train),len(f_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10694.88625, 129146.59506757549)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu = np.mean(f_train)\n",
    "sigma = np.std(f_train)\n",
    "mu,sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of the evolution of daily cash flow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fa51513f090>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEZCAYAAADCJLEQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXm8HUWV+L8neQlbIiTIToAAQYjiAgjixkNIDI4CgyDI\ngKgZxhEXHOfnAs5AMiKCo4OIwqgsBkZZRpRFEROBoI5CWBUISCJGSUJYQkLAQMxyfn9UF7duveq+\nfff73jvfz+d+bnd1d3V1dVedOqdOVYmqYhiGYRi9xohuJ8AwDMMwUpiAMgzDMHoSE1CGYRhGT2IC\nyjAMw+hJTEAZhmEYPYkJKMMwDKMnMQFlGD2MiHxPRL5Yx/l/LyKPi8gqEXm9iCwSkUPamUbDaBcm\noAyjRYjI8SJyt4g8LyJLReQmEXlLk9Fq9ivLV4FTVPUVqnp/A9cbRs9gAsowWoCIfBo4DzgL2BqY\nAHwLOLyDaRBgJ2B+p+5pGO3EBJRhNImIbA7MxGku16nqi6q6XlV/qqqfy87ZX0R+KyIrMu3qAhEZ\nFcRxnog8KSLPicjvRWRycIvxIvKTzGx3h4jsmkjDRsDzwEjgdyKyIHWOiHxdRJZkv/NEZHR27HYR\nOSrbfouIbBCRd2X7h4jIfa3LMcMohwkow2ieA4GNgR8XnLMOOBXYMjv/EOAUABF5J/A2YJKqbg4c\nAzybXSfAccAMYBywEPhSHLmqrlHVMdnua1V1UiINXwD2B16X/fYH/i07Nhfoz7YPAh4D3h7szy14\nNsNoCyagDKN5tgSeUdUNeSeo6r2qOk9VN6jqn4Hv4Cp+gLXAWGAvERmhqn9Q1WX+UuBHqnq3qq4H\nvg+8vsF0Hg/8h6o+o6rP4LS+E7NjvwzS8zbgy8H+QcDtDd7TMBrGBJRhNM9y4JUiklueRGSPzEz3\nhIg8h9OCtgRQ1VuBb+L6rJ4UkW+LyNjg8ieD7ReBMTTG9sCfg/2/ZGEAvwX2EJGtcQLwcmCCiGwJ\nvBEnwAyjo5iAMozm+S2wBvj7gnMuwjkv7J6Z8b5AUP5U9QJV3Q+YDOwBfKYN6VwK7BLs75SFoaqr\ngXuATwEPqOpa4DfAvwILVfVZDKPDmIAyjCZR1eeAM4BvicgRIrKpiIwSkcNE5NzstDE4J4bVIrIn\n8FEy928R2U9EDsicJlYDLwHrs+ukhUm9Evg3EXmliLwyS/MVwfHbgY9RMefNBT6OmfeMLmECyjBa\ngKr+F/BpnNPBUzjz2SlUHCf+H64PaBWu/+mq4PJXZGHPAouAZ4D/9FEzcBxT0bimomNnAXcDv89+\nd2dhnttxgtSb834JbIaZ94wuId1csFBETgNOADYADwAfwhWIq4GdcYX1faq6Mjj/w7jW5SdVdXYW\nvi/wPZwn1U2qemoWvhHOlr4Prp/g2KyDGhE5CWdmAThLVS9v8+MahmEYddA1DUpEdgFOBvZR1b1x\n4zeOAz4PzFHVPYBbsn2ycSHH4mz004ALs4GJ4Oz70zPX2kkiMi0Lnw4sz8LPA87N4hqPM2/sn/3O\nFJEt2vrAhmEYRl1008S3Cudeu6mI9AGb4jpsDwdmZefMAo7Mto8ArlTVtaq6CDce5AAR2Q4Yq6rz\nsvMuD64J47oWN/YE4J3AbFVdmWlnc3BCzzAMw+gRuiagMq+gr+Fs9UuBlao6B9hGVb1b7ZPANtn2\n9sDiIIrFwA6J8CVZONn/49n91gHPZW6zeXEZhmEYPUI3TXy74Vxad8EJjDEickJ4jroOMpvo0jAM\nYxjS18V77wf8RlWXA4jIj3BTwCwTkW1VdVlmvnsqO38JbgJOz444zWdJth2H+2t2ApZmZsTNVXW5\niCyhMq0LWby3xgkUEROOhmEYDaCqTQ+R6GYf1CPAm0Rkk8zZ4VDcQMYbgZOyc04Crsu2bwCOE5HR\nIjIRmATMy6aEWZWNIxHc1C3XB9f4uI7GOV0AzAamisgWIjIOmAL8PJVIVbWfKmeeeWbX09ArP8sL\nywvLi+Jfq+iaBqWqvxORy3FjMTYA9+LGgowFrhGR6WRu5tn580XkGpwQW4ebOdrnxCk4N/NNcG7m\nN2fhlwBXZDM7L8d5CaKqz2aLwN2VnTdTM1d2wzAMozfopokPVf0K8JUo+FmcNpU6/2zg7ET4PcDe\nifA1ZAIucewy4LI6k2wYhmF0CJtJwihFf39/t5PQM1heVLC8qGB50Xq6OpNEryMiavljGIZRHyKC\nDnInCcMwDMPIxQSUYRiG0ZOYgDIMwzB6EhNQhmEYRk9iAsowDMPoSUxAGYZhGD2JCSjDMAyjJzEB\nZRiGYfQkJqAMwzCMnsQElGEYhtGTmIAyDMMwehITUEZPMmUKPPRQt1NhGEY3MQFl9CS/+IX7GYYx\nfDEBZRiGYfQkJqCMmqxZA+vWdf6+0vRk/YZhDGZMQBk12WknOP74bqfCMIzhhgkooyZPPQW/+123\nU2EYxnCjqwJKRLYQkR+KyMMiMl9EDhCR8SIyR0QeFZHZIrJFcP5pIrJARB4RkalB+L4i8kB27Pwg\nfCMRuToLv0NEdg6OnZTd41ER+UDnntowDMMoQ7c1qPOBm1R1L+C1wCPA54E5qroHcEu2j4hMBo4F\nJgPTgAtFXu6luAiYrqqTgEkiMi0Lnw4sz8LPA87N4hoPnAHsn/3ODAWhYRiG0X26JqBEZHPgbap6\nKYCqrlPV54DDgVnZabOAI7PtI4ArVXWtqi4CFgIHiMh2wFhVnZedd3lwTRjXtcAh2fY7gdmqulJV\nVwJzcEJv0LNiBfz+991OhWEYRvN0U4OaCDwtIpeJyL0i8l0R2QzYRlWfzM55Etgm294eWBxcvxjY\nIRG+JAsn+38cnAAEnhORLQviGvR88pPwutd1OxWGYRjN09fle+8DfFxV7xKRr5OZ8zyqqiKiXUld\nxowZM17e7u/vp7+/v2tpKcNLL7UnXu3qWzAMo5eZO3cuc+fObXm83RRQi4HFqnpXtv9D4DRgmYhs\nq6rLMvPdU9nxJcCE4PodsziWZNtxuL9mJ2CpiPQBm6vqchFZAvQH10wAbk0lMhRQgwETJIZhdJq4\n8T5z5syWxNs1E5+qLgMeF5E9sqBDgYeAG4GTsrCTgOuy7RuA40RktIhMBCYB87J4VmUegAKcCFwf\nXOPjOhrndAEwG5iaeRGOA6YAP2/HcxqNYwN1DWN4000NCuATwPdFZDTwR+BDwEjgGhGZDiwC3geg\nqvNF5BpgPrAOOEX1ZX3hFOB7wCY4r8Cbs/BLgCtEZAGwHDgui+tZEfki4LW3mZmzhGEYhtEjdFVA\nqervgDcmDh2ac/7ZwNmJ8HuAvRPha8gEXOLYZcBl9aTXMAzD6BzdHgdlGIZhGElMQBk9i/VBGcbw\nxgRUD/Od78Cf/9ztVBiGYXQHE1AdYMUKN+FqvXzkI3D++bXPCzGtwxjuLFoEjz3W7VQYrcAEVAeY\nMgV2GBLzVBhG7/P618OrXtXtVBitwARUB1iypDsL/hnGcOSvf7XyNlQwAWX0LGauNBrBvpuhgwko\nwzAMoycxAdUBrEVnGJ3DytvQwQSUYRiG0ZOYgOpxbHZywzCGKyagjJ7FTDWGMbwxAdUBrKI1jM5h\n5W3oYAKqAzRjphvOJr7h/OxG45iAGjqYgDIMwzB6EhNQHWAotOi6oc0MhXwzOo99N0MHE1A9jpm5\nDMMYrpiAMgyjJ1m2DJ59ttupMLqJCSjDMHqSHXaAQw7pdiqMbtJ1ASUiI0XkPhG5MdsfLyJzRORR\nEZktIlsE554mIgtE5BERmRqE7ysiD2THzg/CNxKRq7PwO0Rk5+DYSdk9HhWRD3TqeQ3DKMeGDfDk\nk/VfZ31QQ4euCyjgVGA+4HtbPg/MUdU9gFuyfURkMnAsMBmYBlwo8vKneBEwXVUnAZNEZFoWPh1Y\nnoWfB5ybxTUeOAPYP/udGQrCVtNMgRnOfVBW0XQPVfjmN7udisaw72bo0FUBJSI7Au8CLgb8Z3U4\nMCvbngUcmW0fAVypqmtVdRGwEDhARLYDxqrqvOy8y4NrwriuBbzB4J3AbFVdqaorgTk4oWcMYe6+\nu9spGDy88AJ84hPdToUx3Om2BnUe8BlgQxC2jap6xf5JYJtse3tgcXDeYmCHRPiSLJzs/3EAVV0H\nPCciWxbENegZzhpXES+8AG98Y7dTMXjole+oEW3INKihQ1+3biwi7waeUtX7RKQ/dY6qqoh0tajM\nmDHj5e3+/n76+/s7ev9eqSgGOxs21D7H6G0WLoTddjMB1IvMnTuXuXPntjzergko4M3A4SLyLmBj\n4BUicgXwpIhsq6rLMvPdU9n5S4AJwfU74jSfJdl2HO6v2QlYKiJ9wOaqulxElgD9wTUTgFtTiQwF\nVKNYgeo+JugHJ2HZmTQJfvMbOPDA7qXHSBM33mfOnNmSeLtm4lPV01V1gqpOBI4DblXVE4EbgJOy\n004Crsu2bwCOE5HRIjIRmATMU9VlwCoROSBzmjgRuD64xsd1NM7pAmA2MFVEthCRccAU4Odte9gu\nMXNm6yrmRuP5y18av7aVgt0EVGP0Wr6tXt3tFBidpNt9UCG+KJwDTBGRR4F3ZPuo6nzgGpzH38+A\nU1RfLj6n4BwtFgALVfXmLPwSYEsRWQB8iswjUFWfBb4I3AXMA2ZmzhJDihkzYN267qZh553httu6\nmwbovYq21/H51e18sz6o4U03TXwvo6q3A7dn288Ch+acdzZwdiL8HmDvRPga4H05cV0GXNZ4qjtD\noxVEtyuWkOef73YKOpcf994L48bBxImduV+78H12vfQdlcUE1NChlzQoo4W0ugXcTKHvhQqjUxXt\nvvvCYYd15l7txOeXOZcY3cQEVAfoRgXtK5Z6K5jVq3tnvNBg7YMaCpV6r2hQvdC4MbqHCagO0I0F\nCxvVoM4918YLNctQEFCmQRm9gAmoIUqjLeAXX2x9WnqhFdzJTv+hUKkPZg2qF743ozWYgOoA3Sgw\njbaAm62QnnoK/u3fmoujHXRSaHS7Um8FpkEZvYAJqCFKoy3gZiukG2+EL32pOqwXWrSDSYN64gn4\nh39oTVoapVc0qEbohe/NaA0moHqcRiuIRp0keqlCaqSiUXVCMqaTFW6zAuqXv4Qf/KA1aWmUXtGg\nzMQ3vDEBNURpVGPIO79sPKnKodEKoxFhsnQpHH74wPDBJKB6gcGsQRlDBxNQHaCbbuadNvH1auvV\n3MzrYzDPJGEMHUxA9TiDzcTX7Qol7/6DSYPqtlCAxr8fw2glJqCGKI22gNuhQTUqtFrZ/9BJjSB1\nj6VLYfny6rBly+Dzn29/ehqhV0x81gfVWe64A9au7XYqKpiAGqI0WsEMdg3KEz9HJzWB1L122AEO\njWaYvPFGNzC6F2nGSeKnP4U3vam16TE6w4EHwlVX1XfNrbfCZW2a1dQEVI/T7EwS9VYwg92kk/fc\nveBm/uST1fu9IsxTNKNB/eQncOedrUlHvXn0q18N/m+429SrQX3sY/DhD7cnLT0xm/lQZ7g7SXTy\n+devd/8bNsDIkZXwXuiDisN7WUANVjfzt7+9PekwuoNpUEOUTjtJvO99sGJF6yvdRx+t73z/vF5Q\neXpBg4rvXSavVOGRR5pPU700I9C73W9lNEe976+dDS0TUEOUTjtJ/O//wv33N3ZtHr/7HbzqVfVd\nkyegut0H1Wga5syBvfZqLj2NMFg1KKN5eqmBYQKqx2nWzbyTThIirTHx+TQ0sry3F0x5Airv+e65\np/575dFKE98LLzSfnkboFS8+T6+MyzIGYhqUUTeNmvjKmqdSxAKq1bNZlKERE58q7LcfrFtXCZs/\nv/E05KW/rIAKr++WBtMrGpTHBFTn6KU8NgHVAbo5m3kzwkG1NcKi0TQ0cu9aJj5VOP98eOCBgcfC\nyvjVr25Mg4vjKQrvZfNVr2hQPo96JT3DgV7K464JKBGZICK3ichDIvKgiHwyCx8vInNE5FERmS0i\nWwTXnCYiC0TkERGZGoTvKyIPZMfOD8I3EpGrs/A7RGTn4NhJ2T0eFZEPdOq5O0UrTHyTJsE//VP5\na2MNqhtmxjwTXxjnpz5VPf4o9PxrBWUEVFnhF+fFypXw7LONpasemtGg2lHBtUKj27DBTcTbCp5/\nvrcbGM1gThKOtcC/qOqrgTcBHxORvYDPA3NUdQ/glmwfEZkMHAtMBqYBF4q8nDUXAdNVdRIwSUSm\nZeHTgeVZ+HnAuVlc44EzgP2z35mhIMxj9Wq48sr0sXPOgX//9/oyoJ3kFegrroAtt8y/Ljz/j3+E\nX/+6/D3jD7VZT8JmNKj4nkXC0pv2vKBqtoItY+LbbDN48MH643rrWyuOIwsXtq9yaKfGsmZN+XNj\nDaoZAfXrX8NBBzV+fchf/9qaeIxiuiagVHWZqt6fbb8APAzsABwOzMpOmwUcmW0fAVypqmtVdRGw\nEDhARLYDxqrqvOy8y4NrwriuBQ7Jtt8JzFbVlaq6EpiDE3qF/PjHcPzx6WMzZ8JZZ9WKof3Egimu\nYH796+IWeLMzMKQ0qG4IqFp9UGHcsdbVrKAq248XT31U65qFC+Ghh+CZZ9z+n//cWPpqMW9ee/t8\nNt7YrXlVD60QmK2cwqdX+uaGOj3RByUiuwBvAO4EtlFVP+b+SWCbbHt7YHFw2WKcQIvDl2ThZP+P\nA6jqOuA5EdmyIK4mnqGZq1tPnnColc48zSN0IMij0yY+VfjFL6rDannxpYhNfM1Whs06SeTF9Yc/\nlLtPMzz3HBxwQPsniy1r4vR51CoTX1nWrSvW9DoloG680Qn0TtJLfVBdn0lCRMbgtJtTVfV5CUqt\nqqqIdDW7ZsyY8fL22rX9QP+Ac1ThxRfbc/9mK8l6B4fG5/v9v/zFtdh33nngNam4VRuv5Hbfvdx5\nCxfClClpr7d63MxbqUGJlBdQZShKQzsqkjif2lVZjaizadwKgVnPtSedBDffnK/ldqoSf+ih+kyi\nraCxPqi5zJgxt+Vp6aqAEpFROOF0hapelwU/KSLbquqyzHz3VBa+BJgQXL4jTvNZkm3H4f6anYCl\nItIHbK6qy0VkCdWSZgJwayqNoYDKW+V09uyip+wOrZrqKLx+xYpiARWev2FD45XK4sUD750iFkLh\nvVph4mukMqxHQDWrdbezFd9uDareZ2+Fia+ea++7r2IKX73a/V75yoHpaTdb1OwZ7yxPPAHbb5/K\ny35mzOh/eW/mzJktuV83vfgEuASYr6pfDw7dAJyUbZ8EXBeEHycio0VkIjAJmKeqy4BVInJAFueJ\nwPWJuI7GOV0AzAamisgWIjIOmAL8vHaa0+G1zBXdnIuv2XFQ4f4f/gAf/GD+tSLV9y1TqVx0keu/\nS1FvZfTEE65igcY0qFaZ+PJoZBxUpzUoH2erHEbyKKtBlTHxiZQbaJ269hvfgE9/emB4+Nz/8A+w\n1VbpuNopqH7xi+4IqKJ3/sc/di4d0N0+qLcAJwAHi8h92W8acA4wRUQeBd6R7aOq84FrgPnAz4BT\nVF/OylOAi4EFwEJVvTkLvwTYUkQWAJ8i8whU1WeBLwJ3AfOAmZmzREM0U4h/+ctij6AycavCD3+Y\nvq5ZE19YAK+7DmbNIpdQQK1fX64Q/9d/QaCkFqalFsceC//8z5X7F8UV7sdefK12O89LQ5mGS14a\nnnoKHn+8+TTF+Gf3edJODaoep4Va39Jf/lI+jpAvfxnOO29gePiu/vSngcfD/KnleNQoU6Y4qwWk\nrQXd4G9/q97/0Y/c0Id20TUTn6r+mnwBeWgqUFXPBs5OhN8D7J0IXwO8Lyeuy4DLyqYXyrV46+Wg\ng+Dss+G00xqPY9EiOOaYdF9MvZVikYmvDGELvIyAmjjR9SUVxVXrXp7QkaOWBtUuL768a4pMf43w\n1rfCggWti8/jn90Lj1bmQchdd8Gll7qlOYqoZ6BurfumvsOXXqodV+q6sBHztrfBZz/bnrW9/PM/\n9xyMH9/6+FMU5WMsoN773vampSe8+HqZk0+utGJ6ld/9bmBYGS++W27J1yz8+i5lWtD+mkZMfNts\nk3+sXkKz0YYNznzoTRJFnf55Jr5WOjWkTFplteMUZVzUG8E/sxdQrdagfHwrVsCqVeWvK+PFV2R+\nzrs2zwGhloCKNcywf6qV+PLaCUeJMnnc6dV2TUDV4OKL3bgQKK9BrVlTn1dfuG5RI/gJRVN9F0WV\n4KGHOlNRiP84/QqZ9Qio0HOvrAY1enTteMsSCoH16+GUUyrmm05pUCleeiltosnrhwnTWKSRtYLw\nnUFrNKgi/L3y8qTWdUXf0uWXl4sjfKZGNajYLLzddsX3bpa8dLaSMnkca1DtxgRUE6QEA8A73wmv\ne11lv1ZlUiSgylQQqf6CZgfJxvEUERb88L5l+nI22qhcvEXH/X8soAD6+orjCM+N09usB6Tnu99N\nh+d9F2XysF437Tw++9lqLTYWUK3WoHz8ZQVUK+fiS5UT1co3ElLWxOe1wFaPVYq/7U5oUHljCENM\nQPUwcYUydizcfvvAQnPfffX1DzSrQeUVvPDfEz9DvF/kTJBXWYWVaEqDKqpUvAaVOqeWkI0LVEpA\n+bwt48W3dCnMndu4k0Tee6y3UKcEVFxphM/a6KS24MzDfmaK8D7+m2pWg1q40DnCxPHXq0HVMj+V\nSacXuvF9awmXIg3queeK09Uo8XvopIAqGpRvAmqQ8eSTzRfiVAuuHlImqbLaRy2NKWX+ifHeZKGA\nuuOOypiposLr758qgP5+H/wgvPa1+ceLBJQPKzLx+QL5wQ/CwQc3ZuIrOjfWdmqZX0MBFZuSPGHD\nYs89y6czZrPN0vdulQb1zW/Cv/5rZb+sgPrMZ9x/K+fiyxNQKS2+Xg2q3nT94z86D7g84ve+Zo3z\n1G3nkJW8MYQh1gfVg9SqfIq0jjKkWt7e5txqE5//wPMqviKBlffhTppUudafH66nVFR4fTpSWoC/\n7pZbqpfH+Otf3b1iTaceE9+aNZUB1v5cP7ddvZXh6acXd/jHlUqRADzqqGoHlbz3FD7r44/DDTeU\nS2vMmDHptOVpUD/5Sf58lCnib7tsH9RXv5q+rpnGoK9cYw0h7Af961+dy3jZPqjnn3f/9bqBX3IJ\nfOc7+cdTAqrVK1bHlDHx+YZku8bHxdQUUCJyS5mw4YZ/USkBVS9hIX7gAVfZfeAD5a9vxMTnr4kL\na5EGVavCzvPcK7rOF4aUgMqryMeMceOxmjHx3XCD6ysMz6113zy+/GU3JY2nViu3yIT44x9XtuM+\nqPD8+B5HHFEurTF5AirPSWLWrPwZ/VP4d3LllS7NoQaV913cdVdlu5Vz8eVpUKrw9a+7gd5HHw1b\nb12/gCqTrv7+6nsXfV9x+Xz729vnuekpY+Lz9V6rTZp55AooEdkkm1h1q2yNJv/bhSYnVh2shJWC\nt1uHrtWNEgqo174W/u7vyvdhvfnNcNVVbjslTPI8wspqUClvtzzCSjQlLG+7baCJoEhAFeXrggUD\nnyHMx1iryhPYoSYWp8nHMWdO/gh6H1/R0ux5y5CUGefln/Fvf6u2/7fKScILqHD8GuSb+OqZ4Pac\ncyrp9EKnjIlv//0HhrVSg0qNkfuXf3HC949/dMfzBNRb3uI0+jwT3/nnp5dRWb/e9VeHmnaZhlso\nLMK+whSqA71y66GMBpVXb7SLos/8I8DdwKuAe4LfDcA325+03uN92ZDf8OMVqT3Td63KJDaD/OUv\n5Tupf/tb198D+Sa+k092afi3f6scv/de919r1u9wv9Zzxi7LcRzveEe1hhDGWY8GBa6yKaNBeRPf\nY4+l40oJqLgynDrVTWDrp1GK0wHVo+lrve+y7uyhiW/aNNe699QSFE8/XZ8wiTWBPA0qHDgaml1T\n3HLLwD7AZr342qFBhe8h1ZAJ7/mb3ziN2WuRcb596lPVTiHxvcPvpB4NCmq/z6uuanxc4YsvVhpZ\nReU8Lw/BleFWz6iRW5RU9euqOhH4jKpODH6vVdVhJaDmzKneD1vLI0bUrrjjsT7LlrmP7aST3H4s\noNasKbe8RUzeOKiLL3bbv/xl5SM/+GD3X4+Jrx4NKq+Ax/fzcaa8g4pazbUElC9IPm9POSU/zWVN\nfIcm5jfx491C80stT8mynnJh2u6+u/pYLSFYdhokb7LxFYu/36mnVvbDhTj9s5x6atpxJWTkyIEm\n1rJ9UPH9WuHFl2fCCr+zWEAtWjTwnrfcUhknmHKSSL2bUEAdd1ztNKc0lTwPuu22cw2Sp5/Oj8+z\naFF6WqhDD628T5+uj34U3v3u6vOKBNSxxxYvhtoINQ0FqvoNEXmziBwvIh/wv9Ymo7eJW0RvelNl\nu8ycYqNGVe8vXer+8wYW1uuC66nlJJGa8y++T/wsKa2k6P7+fvHSG6k0QrUJKy9tqYK8bl19AiqV\nFr+fJzQ3bKgWDKnWodf8fEs6RS0BlVfphhpUTFkTX61FNH2+ewEbv+MXX3RxxO/VD14voq+vks74\nXdX7jff3u/+8Sj3ULvOIzZieVL+Q/584sVhr82UqJaCeeaaSv6GAuvrq6nukSAmo2NP1vvucSXLZ\nMlcnlVnRe889YZ99Bob/4Q8DtcEf/hB++tPa6fI8+qj7/9//rZ2OspRxkvgf4KvAW4E3Br9hRegK\nHnqoNaJB5VVYnkY1qFpOEn/9a+17NyOgzjrLuRXH9y1ytGiFBpXy4vPx1RJQKQ0qvO8ba3zpXoMK\ntepaGtRRR7n/MgIqb0YSb7Ksha+0nngi/f7yNKj4uPcq9c/y8MO17x0KqJSJb9EimDChOA5/P78E\ny4YNro8nXryxDHlu1GWHZKTw7ycloLbaqmJW999jaBYN7xc62YRpPDuYeTSeTWKffeCww9z2Y49V\ntLlnnoHvfS+d3jVrXGMk/iZDZxmfriJNMPUt+Wd8X3L208Yo0w7bF3iLqp6iqp/wv9YlYXCQVzmX\n0aBiV+f4xcdxv/RS2vyz6ab5Hx4UO0lA+wXU//2fq3TKpMvjnzPV59FMH1QzAqrIyy7GV1D1aFC3\n3+7+ywioopnui4gbONtvXzH1epYurWiA//d/TqOP05InoDzXXpsvLGoJKKgInrIsXAh77w1veEN9\n10Hl3efNWUr1AAAgAElEQVRpzGX6oGL8+w+/ofA7fDJbG9x/3+GacmG8r3mNayT85CdprR6qNagL\nLnD/qTGUZ58NH/pQfppThOPhUg0+T5EG1Y5BvGUE1INAm2eaGrwUOUk89JB7kb5Q33abK5hFY418\nhZr6AF58EX71q/y0pApVGHfKESFOeyygwgqpHpNMWRNfrEGF96gloP7lXyrnrV9fveRInoCKqacP\nKkUjGpSnjIBqdJaI1NxtXkAdeSS8/vWwww4Vc8zMma5PNNVYgspzxs9y9NHVA3Ghkm99fcV9UI3g\n+4PzNEtvCkyRZ+JLNeZSHq8pQg3Ke9OGFbs376caYLEwVIX3vMcJqloC6pOfdP8pAVVGu4XqdxkK\nKFXXhxkLqBdfLNag2jGIt4yA2gqYLyKzReTG7NfgsMChx4gR+S/mNa9xrVL/AR5yiBucl+dqDRX3\n9bxCUdRKqWXiS31UtTSosN+lHgFV1sS3bp0rZP65wueOK44nnqhOp7d5r18PS5ZUxxuOU8tLl9/P\nq7Aa1aB8wReptKBT+Gd98cX0qsz1alBh3qUEwN13u8HI11+fngEf6hdQPp0pRo6srUHVIr5f/P3H\n77NoVvFafVB5GlTREIJQQL3//W67SECFruL+Hv6Z/LRJfX3pmdlTeZYSUHGj5oknajtQhPH88pew\n004Dy86mm1Zmv+glDWoGcCRuHaavZb+EI+XwZMSI4nmyVqwY2BqLhcCPf1yp4Hx/lR/PkJoBwROb\nbGqZ+BoRUCGh2fGFF4ormZSgSd1v/XrYZJNyGlRoHgnTuX79QE3JVxxr1uRrcz5tzWhQPg9iE5+v\nJBYtqr7/NddUtn/9a/d/5ZWVQcMhf/d3FSFchte/vrIdfifh86UGe+66a2U7Fjb++fJMfODy6f/+\nr9JI8Pk2cmT+zCXxt3P44fmd/L6fBQZWgvG7K5o2LDTxrVhRaTyUbdjdeuvAOP17zvPiKyOgfF54\nIaI60GszPG/ZsoHxp87zzhjbb1/b5Bfmo4/fl6lnn630efpGYs8IKFWdm/q1PimDE9Xi/oe//W1g\nRRe/yNmzK+MX8jQn3/cUFuzvf7/6nJSACs9fty6/D8hTJKBCIbLrrvD3f59/bhhPLRPfppumTQf+\nXD92ZPToyvG4Ao5be17z8K3SPFJmtHr6oHz+xSa+iRPd9oIFFZdtqDaHxR5ge0dLbt5/f2W8Whl8\nR/teezlzsiesFGNNE9yCe5644vn4x91/0fIxqm4BxY9+1O3751m1qmKS9t9hSkCtXg033lj5nl/z\nmsoxEbj55sp++F2lBsmXEVDr17vF/448sjq9eRqU55BDBoaVFVC+7IQDdf09/LfsG6XhuwPXeDn0\n0Mp5ocBO9RP5d+Xd2cP75xG+d/8t+7hvuqmiHfr87ZSAqjlNqYi8APjXNRoYBbygqq9ofXIGH0cf\nXSyg1qypfnGrV7u522L8R5UnIP77v93/jTe6MRh//evAubnCQuJV8TPPrIT5vpqQejSoUEA9/XSx\nJ1nohpwSnL/4hfOQW7euokEdcgh88YuVc++5pzrO0aMrhTTM85Qg8abJ+N2kNCjv9h/HV0aDyhNQ\nvrK56qrqPE21eL3QSM1A0AiPPFL9zX0zGLWYcmIJ05RnxvXf59/+Bv/zP3DCCZVjPp823bQ6juuv\nr5zj8yfVcPINHX9s661h883TjYu4EqwloFQHjqXy7yzW4lKaVMikSQNneHnxRffc4bWhlhlrUPH9\nTj21Mpel/2Zi124fj09vONg3pdGmGhMvvAD/+Z8Dwz3he4/Nyttu6xo98+alLR2pOFpFGQ1qjKqO\nVdWxwCbAUcCFrU/K4KRIOIFrfTz8cGVBs7vvdrM/pNhjj4EF8OqrnVkmDL//fvjEJ6o/VKguVJdc\n4v7Dec3CsUOecH/duoGVdYgXDr6/JJ4JOyQc3JwSUFOmwMc+5u6/dq17nltvhZ/9LD/OjTaq5ENo\n5li5cqA3mD8eT+Ka6oOKr/V5kiroqtWVjd/O+w7i95kSUF/5SvpaTz0zWHszX2jKu+iiynZqLFcZ\nAXX55c4stnr1wPceC6iUY8cjj1THHwoHP2zDfxu+0QIDG2FF/bdQLaD+9jenBfhv0afTr8+1dm2+\nl2dqyiD/fCEvvujSWtbEF9/rG9+oTIzrTXyxWVfVNc58nsVaZEz4Lfr3vXKlW/sr5sEHnat/aEnx\n23/6k/sfPdo1GMK8uv/+YoetVlHXjF6qukFVrwOmtSk9HUVEponIIyKyQEQ+14o4V62q/kDmzXMf\n1vjxbr9oPq0FCwZWoC+95ApY2EL/xjfS12/YULvVXySgwjEXKXxle/jh7j9VYFPkmfhWrXKFYcmS\nSuErcgqYMwfGjXPboVA5/viB87d5W3ksOOI+mH/8x4FmL58nBxwwMA3f/a4rsLfd5iqH977Xhed5\n8cX9k/X0KXm23bb8ud75IdQ+wjSkOvxDAeVNejGXXOK08eefh1dEthP/TjfbzDUMUjO7e237uuvc\nfyig/Dt//HE3bVcooOLvOa7oizQon9dTp1bHFQqoUJiEJr4UeQJq002ry1FYkY8a5cr8r3/tKvmQ\n2MTnBVSq4h81qnJeLW+50DHHz+ywYkX63L33dmWpSPtZv959N9tvXwk75hg3gW27KTNQ973B7xgR\nOQeoY0Hz3kRERuLmFJwGTAbeLyJ7NRuvL4Axfs0Zv6SDJ3zpkG5hf+Qj5SaP3XXX2vOjxX1Op57q\nCsqaNdXmwEsuGdh6jTWKsgLqsceql7IIvZfiglEkoLyTxEYbVZ+Xau3maVD77Ve9f911roCGA0aL\n+p585/VPflIdHnoYhnmc6vOpl0bnV/OaTp6A8oOQw28u/j5DlixxAmrs2Opw/z5Hj3aWgkWLYIst\n0nH4mS3CPAorz0svze+HPeGEgRpp3KgKBVTsgZrqC44nGE6Ny/JCLGUxSGlQYQNl9Gg3zdZpp7mB\nuyE+Pf5bzfO0U3XvqGjWlVrUmiOvyJXeC6gdujBFeBkN6j3Au7PfVOB5oMHJ/XuK/YGFqrpIVdcC\nV9HEc910k/uPC7gvML4VHFemoQcVuI8vniuraK2hmIULi49femn1/qOPOm3Az83nERnYARv3CfT1\nDezQTXHJJbDLLm77sccq8a5YMbBglHGr9nF5YqHe11cR6GXybvHi6vdQNEuBb32H01/FFVdYgeTN\ngl4PW2/t/utdVtyPlQlb3Cl3+Himkzx+9CPXLxgLKO/d5r+PBQsaX/4DXHq9eSlk9OiBlfM551Tv\nh99sKPhSC4uuW1d9fl5jwpfhVINs1Sr3/jdsqGh9YZruuKMyLilPQHnyZiL3AspTqw8qfj9HHVV7\nPGCRgPrmN10ZiRvTnaBMH9QHVfVD2e9kVf2SqjYxqXvPsAMQTqm5mCaWEXnhBfexXH+9U6v9h+/N\nIW99q/uPC1jKC6fW9C/18rWvFR//5392/WJ+dLpPV958b95U8bOfuRnKi9hpp+r9sMN+6dKBGtR9\n98GrXlVZtM8TmpX8Sr2euGN86lRXqHfeuTLT++67VzzrYpYurRZQ8cDTWhx5ZGXi31YRmvW8gKql\nHcekPNq8BhUOaM0bJpH3HcYmPk8427s3xdbLtdemZ42HcoI0rLDDAatLlgwUCIsXV/eXrV5dyetU\nnKn7r1jh5rfbsMF9t+96lwv3C0hef33F+SUeoxWb3X7+8/QzQdq5JY/YJHz66ZUuhjxSDik77+z6\niq+/3s3iXqRBhfkeLzbZDGVMfBNE5Mci8nT2u1ZEdmxdErpGCR8tgBmMGDGDz31uBjD35dCwkH7h\nC84VVNW1ML/wBef1ApUXV7bAfuxj7j9sMYXzdPnVTL3AK8NJJ7m0XXCBM5PEiyF6bSP8sEeMyO+c\nX7iwehLSM85In7frrvAf/5Gfrqefdi36sGU6f74r5DsWfGETJri0ffCDzkwVtii/+EU3qzJUWpJv\nfSvMnVts5ohbnSHxeDNwQun889326tVuGEDRhKVFAmyvhGE5NBmFa4+FY3HCmaN9P0tI+A3tvrv7\n9xr+t75VOea11riCzOuPycure++FAw9022GFWKv17vnmN4tnwy4joMJG1WmnVbZTAipm9er08u9e\nM8oz/e6+e2U8nddEUqvfxvlWVrtWzV/FOGV2/MpXXJr/6Z9c2dh333xzfH+/i8OXjTD/Dj20uq8w\nbmxWp3EubsjsDH7/+xn5J9ZJGRPfZbg1oLbPfjdmYYOdJUDYRpyA06IiZjB69AzOOWcG0P9yaPhh\nnHVWtcAKPwZfyed5vMUahP9Arr22EjZ5cmXb3yevFRszaVKl0H/843DFFW5hthRjxlT6WGITn9fC\nvvIV1xIMzU1bbpnuJ/HjTfLYZRdXSP3cdJ6NNy4ed7Pppi4/x42D3/+++tiGDRUNz+fRHnu41l/R\nmKg889nuu1ebP17/ejfocfr0ignNV1zhIFyoOFDAwFZt2JE9YYLTYL3pcv786saB14TidxIKmVTH\neahBnXyyS8/99ztT3eTJlXv4vJ46FXbbrXJNOA4mNOVuthnMmFF9rylT3P+pp7p3EvZB1RJQXrC/\n+92uIQFu8Pp551WfV6+AClm8eKCA+tKXqvdXrXJ55gdQe7bc0uVZSkDtv79rCPiZ572FJDWDSEr4\n+sZULT7ykep9v6JC3Ijcay/XeFq9Gr797UpZyMuXV72qOl/DcrDVVtVlJh6nV00/XkB95CMzik6s\ni1JTHanqZaq6Nvt9D0gowoOOu4FJIrKLiIwGjsUJ4gGkXm6Ri3XYafqhD7mPMHzxp59e6aCObb++\nMOd9UN5dvayASrUI8xg7ttJKik18vvXn0xdWfiNGpL2Y1q8v1hx9q378+GrvsVGjBgoo34r19/YC\nKjZPrV9fKZRxmosI4w9bq2PHVptixo1z/XjhOjn+2UWq+x/C6WrCdG68cXWaRo92FY43u40alR5L\nE7+TWvPFhe+or6/SGo7XtQrzOtT0/L0OO6yShk02ceHveU/lvD33rHw3fX2uIgufb999B6Yt5Kij\nXB7vvHPFDLb33gP7Gst8y3nl5oknBgooX4a9yXjVKpfu8FvwccbfuGfjjSvH/JAJcP3IcaMkJaDy\nzKjhYGWf7vB6/y7j7z+vLynPGjJmTHW3Q/jsr3xltVD2af3P/xw4pirUrsr2aZahjIBaLiInishI\nEekTkROAGosP9z6qug74OPBzYD5wtaomp1lMvdwiDzbfogKncVx1VbW55cwz4dxz3XbcKvMFLO+D\n8qavUEDF66/EFVNZxo6t3Dc28fkP11c8YQVUJKBi19oQL6BGj672rgsrU3CryXrTp4/fC6h43FR4\nT7+EQJm1k8KCGb6Tvr5q1/SU6TF89tDMF+bRn/5UKcRjx1a/F1/JeGHe11cdZ54G5SuuMWMqS3iE\nhN+caqW/xQtuf30ooP793+HEE912ajyPNzGGx+68s/KtxDP3v/e9bpl0cPNQphg5sqJx+cpt9eqB\nZSDl4RqT965ffHGggPLn/uM/uv88AeXzPeU95wVUrEGtWDGwcZYSUHnejuFsEZ7QbOgF2Ny5lXy5\n8sr8fqy8fAmnGRs5stKQPvlk11gJy8IuuziPxP/3/wb21YWzy9fTKK5FGQH1YeB9wDLgCeAYoM7J\n3HsTVf2Zqr5KVXdX1S/nnecLyg47VEwg4fopManxSGHhGj268iHkmfjiD8qPyveeNGEh8hWAJ/xA\nahXqUHMZM6Z6otMwDXFHcVj5jhxZ/Rz+3NQceSG+hTlqVHV+hpPHgtM2w+MjR7r9ceOc8AoJBVSq\nEXHggdWmMU+Yn+G76+tz07x86lNuPzUZaerZ/bWevfaCz2Uj7V7xiup88a7N/hlDl2K/7+MO34mv\nPGbMSHdgh/cvs2SEJ9UI8Wnw7z8WXrGA8prEhRdWzs1bWyv+1m+91VXAsYAq09DIO2fNmoF54OM/\n7ji3Sq438aU0qJEj0wJqo43cMd8H5Rsbzz03sI5ICaiwAffqV1e2w3eXWqbn7/6u0jfq83qvvfId\ngfIavGH90NdXqZdOPdV96/FyP77sxPkc7ndUQGVu2O9R1a2y3xGq+pda1w0l/FiLBx6o9NHEH3FI\nuLKsJ1Z7d9vNfdh5Jr68wrnNNs4RIxYQIaE5MU9ALV3qPLrCxcX6+vI1KL/tP74iE18451mRgPLp\nHDWqok35uMPKoK+v+qMfMaKiQcWsX1+sOYmkXaDDPIs1qH32cf0hzzyTHsycNwbG59FrXwtf/nIl\nL2IN6qCD3L9/V3191X1KeRpUSJjP3jkg1qDynAQOOQRe97qBcaU0qNDc6Bk9Ol9Abb11Jb6NNkp7\nfcbPdPDBLr4yAiruuM/LnwsvHOhdFjbGNtusokHF5nOf76khFbEG9fWvu/AHHxzoFFFLgwrriPjd\nxc82evTAslXkKp6XL3191ZP7+nLgz8/7ZnpGQInIriJyXubJNyyX2/hEtjzjuHGVD6qo4k0JKP/B\neRfmrbdOT96aZ+ILC9M++1RXcHGFUkaD2m47Vyj9ta94RcU7zsfl4zv66EoFlmpBxxpUKKCKWr0+\nnaNHu/i//e1KmosEFDjPwTe/eWCcm2zi8vbhh/NHx8cOERdfXKxBebbcMu1MkSegfN7Ghfy88yrH\nPvaxylREoRDI06DyZmYP89mnOXz3RbOMnHZatfmoSEDF7/+3v62etdzfO0x/+F3eckulle/7wvLK\nUpl1teLvouy0UC+9NFBA+YG748dX553XoFKEAmr9eud04PENJd+oaVRAeeIBwK0SUOG2LwfxMim1\n4uuagAKuA/4EXEBluY0aI2uGDnmFuqhv5y1vGfhi99jDfWBxKy6uREMBERKb/lIalK88y2hQ8bXH\nHFNdAYat9fe+txKnL0Rx5Zg3eaSPP9VxGmpQUK2dhRrEyJEDr582LW2/99rDnnum0yQyUMgcdli+\nBhV6UOZRxsQX0t+f1pRraVCxk4RPp0j19+DPj81E8becV5n768qY+OIGVaxBhfH4a33jw3uTljHd\nheeFy3LE77LMDPTgvrXwW/fmYJ/+O++snOsdQ/LiCQVUmOdeQJ18svuvJaDC4QZlNKhw/6c/zTeh\nQv67juuRejUo77QVxtNpAfWSqn5DVW8Nltu4vfZlQ5uiymf33Qe+2O23TxeeSy91LuXe5TUsNCFx\nhRBrMFD5uOrpg4q98lImvhEjKpVLSkAVaVA+/lR/kE9nWDH7/1CDOuSQ6mcqaiWHFVaegIoL0MYb\npyuE3XYb6OqcoqyASmk8YVioQYVu+7WcJML4wvNrmfjyKp96THyxoE0JqFiz9+82r781L33+vNDL\nMn6XZQVUGF8ooFINh732yk/jxz9ebeILK2rvJej7mYr6oF7zmsqKBZCuX2IB5fNnwQI3drBI0Der\nQcWz28Tfbzc1qAtEZIaIHCgi+/hf65IwOMmr+OtZSwice/FRRw1cgqPIxAfV6nyegJo0Kb0QXkjc\nIZ4y8dUSUL6AQnUr7sgjK3GkBFRcEYbCMhRQ++3X2EefMnmIVNvdYaDbtz82dmw5z7E8ARWb+FIC\nIbw2FNT33lvRFPKcJOJGQrwdO0mUWT4kvL6Mia8ZDSplDQjJc2oI30ncD1mPi3NYpmLnH5+2c85x\nv1QaP/axiun7ggsG9rn6yn7UKGfmGzPG9WN7hxuoaFDjxuU3LP17C7+tjTaq5E/Yf1svcVeBr0Pi\nbzd2h/fHfRraJaDKOCG/GjgROBgIP5mD06cPD1IV7h/+UN9aQkXUMvGlWqixgCozc3aqMoPqyjA0\nsfn/1PICIm72dh/X979fWb4jNW4sfkaflrgPKryvv08Zyq5P4z2xPKHprAy1nCTKXhsKqG22qXhs\nNqpBxefWK6DKmPjyNKhUH1S8eF+jGlSYrqOPdpW7HyRddgLjMM2hgIotCccf795D0dpnfjDrunXV\n7zysqL3p+TWvqRaq3pEizoNUwyjWoOpZfymvwRwLqHr7oFICquzMIWUoo0EdA0xU1YNU9WD/a10S\nBh/33pv25tpjD9f3AfWZGlL4AuJXE41bqqkWqi8Q9bRg8iqaek184X+q4gwrjk99yrUk40IZmvji\nmREaaZWVKcA77ODSnXLrbVZA5a01VOvauB8y/M9zkkgJk/Dcer7HdmlQ/lhZAZU3BCPubwtnOCjy\nro2JnQ7CuON7FVW6fsxePDt63vfj39tvf1vRoIrGfOX1QdUjoPLO7etza33NmlXtZl5WQPm0NVvf\n5VFGQD0ANDj149DkDW9wAzbzBtlB8y/MfwDeBblIg/KFwhfOema9LuqDqsfEF8aVij8UUK9+tWtJ\nxuenTHx+iY16Bhx78vqgQryJxt97v/0qsymU6bzfYgvnRp4XPxRrLkWVTFz552lQcXg8WPbkk50T\nTCtMfP4/zzTs793fX5mA18cTLhcea+gp8kx8Ybpix5FGBJTIwL61OM1F34KfgDfUoN7znvw6wL8H\nP13RtGnVy7OH9w3Pjxs/rRJQ73+/m5+zESeJdguoMsV+HPCIiNwF+Ik1VFUPb0+ShgbNmvji1nRc\nOItMfI0IqLiiCT36RoyoFJh6NSgfFpr44oo3PjfUoN7//uprUtflUUZA+X1/7yuvrNj0y9znySfL\ntZqhvAYV74f5VaYPKhbmfgaHvAo/psiLL66w8zSov//7ylLuKcEycmRtDSpvCEYcT5jOekx8Ydp9\nGmNhkBLWMeEaW+H5tQSUJ7WK9KhRriEXThQdv/tWm/ga0aD8cf9fZrmceigjoPwydgoI8HbguPzT\nDWidic9TZOLzjB7tRn+X6dj3lNWgYkGZ0qBSrS4fFs50kFfoQ++zuA+qyMSyfLlz7ffLinvCOfI8\neUIgfv5U+lLEnfIpE1y7NCi/vHuc1rzB3jG1vPje8pbKcux9fW5G8FjY5GlQIamyED5LvX1QRRpU\nqM3WIuUxG/epFJn4/HWhgArfWV4dUKZuGDXKDfYVSZv46hVQRRqUJ9UHVUuD+va33ewZfgHPehoI\nZahZBNXNo74Kt2DhLOAdwEWtTcbgpKjiaZWTRF5/hBdQs2dXrhk50i1hkTfdSYqyXnwxqdZ7qnLy\nx6ZNq6xp4+OL57XL8+KL0xBXvOPHpysQP41T0bsoaimX1dRS8aXuW0uDquUYM2JEZT69DRsq8xfm\njYPKW/m2Fj6umTOrJ8oNF6yrpUGF1BJQZUx8Bx2Ur4mFeb733tWLMhaRElC+Im/ExBcyYkR+fpcR\nUKl89OmcNs156DaqQYXpCic9TmlQtRox73iHm7uxnrTUQ262i8irMvfyh4GvA38BRFX7VfWCvOsM\nR6s0qLx/L6D8MgdQ+WjOPrt6gtMi8jSo2MQXE1esF1+cngw0rIS8QPJhb3hDtSAK03LCCdXTMNXy\nDEql8cgjB4bVMmu1UkCVoYyJL/x/3evcIox5GjZUnmX33auXEW+mDyqmHg0qdd8yGlSYN3Pn5mtQ\nefkWklppukiDqmXi22MPmD7dbV+UaK7XY+JLUeQkccYZzmmoWQ1q5Mhq0/u3vlXpky1r4qt1XrMU\nmfgeBn4CvNPPvScin25PMgYnRR/aVluVb8mlqNWaTpn4/LFRo2qvoOmpV4Pygw1jM54vrHnx58UX\nFsRQQIULIkL186QqoLKCoVYfVFz51YuP78gjK95lRaa+oulp8vqgDjgg/1yobmyEk9u2Q0A1qkHV\n0tCh2MTX318RWmUEVMoLNNUXGmtQeSa+r361YmLdaSc3xinui8oTIM0KqEaEQiotcTre9ja39leZ\ne8T50w0BdRTwfuCXInIz8L+4Pigjo+hD+9Wv0kKkLLVcjlNxNzL+oGwfFFQ/b9jyKrpvWMBTLdaQ\n1AwInl13deatVtu4YwFV1F9VT3xTp5a7PnyPtfod8+KLw/PMUs2Mg4qpR0DV0qAa8eK77bb8fCnb\ngCmjQeXlZZzmcOl4f32zfVB56fVpOfLIymrYtSgrQFL38Kswh3Rdg1LV64DrRGQMcATwL8BWInIR\n8GNVnZ13reE0qGbIE1B+/73vrc+RII88U00tE9+HPuRMdG98Y30t7SKK+rKg0oGbqmyK4i8jdFIV\nUTMaVOo+qYo6tYRDXlxF6QnjzktDvQKqSMB2wkmijBdfrXTG16auSwmo+Lq8IRHhdfEM/82Y+Iqm\nOvL/kyblr7EVU0aDSt0jL/74vI73QXlU9QVV/b6qvhu3LPp9wOfbk5zBRbOOEEXktZ59+Ic/DHPm\npK+ph3o0qJC+vsr6NUUeTikzSq2O11rT1XTKxNeMBhV7XIX/IaEGFbeai/KwiLzzy36vZe7XCieJ\nWoK3zDioMJ74vFR6U+elTHx55kXvEp7ScEKrwj77tF+DqodGNag8uq5BpVDVZ4HvZD+jjdQy96Ro\nxsRXayaJomvLai+18JVb2eXsy5LSLvLu3ayA8qTyZMqUgcuthxrU0UdXu+PnOcjEiKSfsZGKrCx5\nGlTqG6wloPKeq8xUR6nrW2Hiy9OkfOMpFiDTp/PyWlf+Wr9Sb0yjfVDNUG/jpNa3E5f9nhBQRueo\n1QeVohV9UJ5Qg8qLt8wYkRR5x3w8RcvE513faGXcrj6oVHomTx6o9cYa1NvfXtnP01JSpMaeNapB\nlals8tKW+hby+qDqTUdZE1/Z76PISSK+dzydWFxe/Oq2YXyt1qDi+OvB33PmzOLz6tWguukkYdSg\nmya+omvqocjNvJZQrNf05Kll4qsloFIcdVR6QlooFjpFJr5m+qBSJr4UneyDKluJLFlS+5yivsuY\nPA2qFmVMfGXjKqtB+QmW87S3PA0qRTNefKmpjprBx3HGGcXx1mqUxud1vQ+qHYjIf4rIwyLyOxH5\nkYhsHhw7TUQWiMgjIjI1CN9XRB7Ijp0fhG8kIldn4XeIyM7BsZNE5NHs94EgfKKI3Jldc5WI1DH3\nQoV2Cqh6TXz//M+VhdHqIZ4rrx4TXyNaUhFem2hkYtjTT4fbb699XrtNfPWa14o8Pb2ps5YprNk0\nxNeC+ewAABlaSURBVDz+eO1z8jSoss4hZTTVMhpUK/ug7r+/ouHm3Tuek7CIVo+DaoaycZT9duK8\n6+Zkse1gNvBqVX0d8ChwGoCITAaOBSYD04ALRV7+3C4CpqvqJGCSiEzLwqcDy7Pw84Bzs7jGA2cA\n+2e/MwNBeC7wteyaFVkcPUFeYa2lzVx0UcUGXg/xXHn1aFBFNGJeKjuPV9kxXmXohJNEEUUa1NSp\nboBprcZJ3AeVd/5BBznPr1p89atumEQReRpU6rkb1aDiefHyhGCrBNTrXufGAgHssotr9MWkVivO\no5mZJIo0mGb6Rj2nnw5f+EJ+3LXuMaQFlKrOUVX/SHcCftKbI4ArVXWtqi4CFgIHiMh2wFhVnZed\ndzng5wk4HDcFE8C1wCHZ9juB2aq6UlVXAnOAwzKBdzDww+y8WUFcdT5HI1eVi7MRJ4lG8PHG5rEy\nndiNkhff3ntXpifKY9Ei+OhH67ufn1U7dW+/3+hcfDGtNvHttlu595ByrojTf/755dYI23FHeOtb\ni89ppQaVR966XLH21UonCc/GG1fPEOGfwQuOZjSoevr4wnt7GimPcRxf+tLAwfCp88oylPugPgxc\nmW1vD9wRHFsM7ACszbY9S7Jwsv/HAVR1nYg8JyJbZnEtTsQ1HlgZCMgwrrropImvXZ5ZPt56limI\nKVMBhefk5dvYsW5l0iJ23rn4eMwLL1Q/W9k+qIsuqsx1Vw+t1KA8ZTTZXXaB3/wG3vzm9jVmQuLn\nLCoLjWpQOzRUKptzksijlQKqTL3R6nLeLgHiaVcfVNsElIjMAbZNHDpdVW/MzvkC8DdV/UG70hFR\nt0iZMWPGy9v9/f309/e3MDkDadTE1yzxEh2t8mbrNnmOEzGxgEqZd8rQDgFVS4PKE7rtfgczZtS/\narCnzPfll+uoRZnnLDsOKo94lvNOCqh2aFDNnhezcuVcYC5BddkS2iagVHVK0XER+SDwLiomOXDa\nzIRgf0ec5rOEihkwDPfX7AQsFZE+YHNVXS4iS4D+4JoJwK3As8AWIjIi06J2zOJIMqMgx3vJSaJZ\nYg2qVwRMq6nHxNdM/GE8ed/JpEn1a4RFeGHc7saM58wzK9v1zuxfj+ed/8+Lpx0mvpgJE+Dqq2vP\ndhIyfXp6UdMy2kzoKNRJJ4lGGTOmH+h/WUDNrOXPXpJuefFNAz4DHKGqLwWHbgCOE5HRIjIRmATM\nU9VlwCoROSDrQzoRuD64JlsDlaOBW7Lt2cBUEdlCRMYBU4Cfq6oCt+GWsie79rq2PGgTdMrE5yky\n8dUqwPWa+LpJWRNfs/GXieeBB+CnPy0fZy1e+1rXR9cpDaosH/0onHtudVg9aSsSUKm4GnWSKGLE\nCDe7fj0a1LRp8N3vDgyvJaBe+UrYZpva8ddDuzWoodYHdQEwGpiTOen9VlVPUdX5InINMB9YB5yS\nCRSAU4DvAZsAN6nqzVn4JcAVIrIAWE62mKKqPisiXwTuys6bmTlLAHwOuEpEzgLuzeLoKTqpQX3m\nMwPt/b1SubWavOfqhoBqxJ2+FjvvDI89Vp2WTlBUsU2aBJ/9bHWYT9t22w1cF6weUhpU3nl5YfXk\nUz19UHnUEgK77VZ8fi+a+IaUgMrcu/OOnQ2cnQi/B9g7Eb4GeF8cnh27DLgsEf4nILFoQX10Uito\np9nmK18ZGFZPIWjUhNNt3v72Sj9Hq/K1Xi++VhD3s3XKxNcMPp+WLi1/TVkNKkWzThLxNe0UUCE/\n+EG1FyqUMy82c89GGFICaqjQjT6oTlU6Q8VJoohwYK+I84JrpuLx8UDn3tODD8Kee6bT0CsaVIpG\nTHx5xzrhJOGppw8qj3oq8/e/v3r/t7+FV72q/nu2W4MadF58RnN0yzMrdf+ie+61V7kC06t9UCEH\nHti6+DsloPyM8iGdcDNvllYJqLJxtcrEFzvVNEIz2sab3tTYde0ufxdc4Po/W40JqCZo9Uu/5Zb8\nsTfd1KCKmD+/ues7TbvT1Q4TX71p7rSQbIR60lbkJDFmDCxfXj6OVFg9+Zs3eLgePvKR1s6GUoZW\na1Bx/+nBB9eXnrKYgGqCVguooqmKuumZ1ayJ76abBi4z0S26IaCapVHz2VAx8eWxYIFzKJg1q7Hr\nmxFQzfDmN7tfJ2m1gJo8Ob9x2kp6uI01vOm2BtLK+x92WKVv59Of7h1h1Q56QXspq2030pfRKlrh\nibb77uX7oFI0YgodP779+daOst+OPqi99mosLfVgAqoJOukk0WmacZKYPDnfVPm1r8FWWzWermYZ\nTia+WtddfDE8/3xjaYrplpNEvXGlrqvn+o03hkceaex+3aRX+oDrxUx8PUq3BZRn331hp53qu+be\ne3sn/TGdElDdfP6ymsHo0ZX1jTpNK8fydFJADVba7cXXLkxADRI6XYj8/e6+u/5r2zEAdTDT6UI/\n1Pqg6p1Jot54e01AtXOVhMGGmfiaYLiY+IYS3XDT7xa9kIY8BquJb7AyWDUoE1BNMJQF1FBlMAqo\nwfAtdMOLr9m4ejVfu5musrP/dwoz8TXBLru0f50Vz+abw7apxUvaRK8W3sFCmH+91irtBVrZB9XJ\ncjHUedvb3CrOvYIJqCa49972xR0X4I03hieeaN/9at1/qDAYNaihSCtNfEXjB436EBk4WW03MQHV\nBOPGtS/uLbdsX9y1uOYaN8ZkKNINAZVaE2io0U03c3Dj7NaurS8Nw43B2GgyAdWjjB3bPdPQMcfU\nPmew0mkBtWgRbLppZ+45mGilia/WscFGu75RE1CGYQCVyqAVK+YOhoqlm04Sjdy/HWnodXp5bsY8\nBmGSDaNxBmMf1FDSDjytNvGl8uhLX+o9r7QytOt9D0aBbALKGFYMxkI6GOhEH1TRHHip+2+9tb3v\nwY4JKMNoA90c5zNUK+W3vS1fENYbPhwZjN+FCShjWDEYTXz10o1Kudt9UN26RzswJ4kKXRVQIvKv\nIrJBRMYHYaeJyAIReUREpgbh+4rIA9mx84PwjUTk6iz8DhHZOTh2kog8mv0+EIRPFJE7s2uuEpEm\nF/o2BguDUUANxoqlFq3ugxo3rv5FAIdivhYxGJ+3awJKRCYAU4A/B2GTgWOBycA04EKRl7P1ImC6\nqk4CJonItCx8OrA8Cz8PODeLazxwBrB/9jtTRDbPrjkX+Fp2zYosDmMYMBgF1GC6d1laLaAefBAe\neKC5+wxltt0WJk7sdirqp5sa1H8Bn43CjgCuVNW1qroIWAgcICLbAWNVdV523uXAkdn24YBfU/Na\n4JBs+53AbFVdqaorgTnAYZnAOxj4YXberCAuw2gJw61irNfE11fHAJcyebn99u5npHnoIfjNb7qd\nivrpyjgoETkCWKyqv5fqr2974I5gfzGwA7A22/YsycLJ/h8HUNV1IvKciGyZxbU4Edd4YKWqbkjE\nZRgtYbgJqHq491545Su7nYrhRb3mz16hbQJKROYAqWkcvwCcBkwNT29XOiLq7j6eMWPGy9v9/f30\n9/e3MDlGp+mU4OjmoMhed5J4wxval44QayR0jrlz5zJ37tyWx9s2AaWqU1LhIvIaYCLwu0x72hG4\nR0QOwGkzE4LTd8RpPkuy7Tic7NhOwFIR6QM2V9XlIrIE6A+umQDcCjwLbCEiIzItascsjiShgDIG\nP8OhD2qoYV5tvU/ceJ85c2ZL4u14O09VH1TVbVR1oqpOxAmafVT1SeAG4DgRGS0iE4FJwDxVXQas\nEpEDsj6kE4HrsyhvAE7Kto8Gbsm2ZwNTRWQLERmHc8j4uaoqcBvgZ5w7CbiurQ9t9AzDQUB14962\nNprRDnphLr6XP21VnS8i1wDzgXXAKZlAATgF+B6wCXCTqt6chV8CXCEiC4DlwHFZXM+KyBeBu7Lz\nZmbOEgCfA64SkbOAe7M4DKMnqbeC3nNPOPnk9qTFMDpJ1wWUqu4a7Z8NnJ047x5g70T4GuB9OXFf\nBlyWCP8TcECDSTYGMYNRg9p3X9hnn/Lnb7YZfOc7rbt/GUyDah3D7XmLsJkkjGHDNtvA29/emXu1\nspLZcUe4557WxdcOBuOUQiYIep+ua1CG0SmWLevMfU44Yegu+NgNhpsgGYzCvl2YgDKMFnPFFd1O\nQeexStVoByagDMPoWaZPhx2G2TD64aYxFmECyjCMpmmXBnXxxY1fOxgr+i23hP3373YqegcTUIZh\nDEt6UYA99VRvpqtbmBefYRhNMxj6oL7/fdh4426nopgRI0xAhZiAMgxjWHD88TB2bGXfBEHvYwLK\nMIymGQwalDH4MAFlGMaQJKUhhYJ0333h29/uXHqM+jEBZRjGsGT0aPinf+p2KowiTEAZhtE0ZuIz\n2oG5mRuGMWw47zx45plup8IoiwkowzCaZrBoUCec0O0UGPVgJj7DMIYk5kY++DEBZRhG0wwWDcoY\nXJiAMgxjyLL55m52BmNwYn1QhmE0Ta9qUPPmwbp13U6F0Shda1uIyCdE5GEReVBEzg3CTxORBSLy\niIhMDcL3FZEHsmPnB+EbicjVWfgdIrJzcOwkEXk0+30gCJ8oIndm11wlIqM68cyGYXQOEdh+e9hp\np26nxGiUrggoETkYOBx4raq+BvhqFj4ZOBaYDEwDLhR5uavzImC6qk4CJonItCx8OrA8Cz8PODeL\nazxwBrB/9jtTRDbPrjkX+Fp2zYosDsMwGuT44+GYY7qdCmOo0S0N6qPAl1V1LYCqPp2FHwFcqapr\nVXURsBA4QES2A8aq6rzsvMuBI7Ptw4FZ2fa1wCHZ9juB2aq6UlVXAnOAwzKBdzDww+y8WUFchmE0\nwCGHwDXXdDsVxlCjWwJqEvD2zCQ3V0T2y8K3BxYH5y0GdkiEL8nCyf4fB1DVdcBzIrJlQVzjgZWq\nuiERl2EYhtEjtM1JQkTmANsmDn0hu+84VX2TiLwRuAbYtV1pCejRrlzDMAwjpm0CSlWn5B0TkY8C\nP8rOu0tENojIK3HazITg1B1xms+SbDsOJzu2E7BURPqAzVV1uYgsAfqDayYAtwLPAluIyIhMi9ox\niyPJjBkzXt7u7++nv78/71TDMIxhydy5c5k7d27L4xXtgn+oiHwE2F5VzxSRPYBfqOpOmZPED3BO\nDTsAvwB2V1UVkTuBTwLzgJ8C31DVm0XkFGBvVf2oiBwHHKmqx2VOEncD+wAC3APso6orReQa4FpV\nvVpE/hu4X1X/O5FO7Ub+GIbRHCIwaxZ84AO1zzVaj4igqk3P5dGtcVCXApeKyAPA34APAKjq/Ex4\nzAfWAacEEuIU4HvAJsBNqnpzFn4JcIWILACWA8dlcT0rIl8E7srOm5k5SwB8DrhKRM4C7s3iMAzD\nMHqIrmhQgwXToAxjcGIaVHdplQZlk4AYhmEYPYkJKMMwhhzvehccdFC3U2E0i5n4CjATn2EYRv2Y\nic8wDMMY0piAMgzDMHoSE1CGYRhGT2ICyjAMw+hJTEAZhmEYPYkJKMMwDKMnMQFlGIZh9CQmoAzD\nMIyexASUYRiG0ZOYgDIMwzB6EhNQhmEYRk9iAsowDMPoSUxAGYZhGD2JCSjDMAyjJzEBZRiGYfQk\nJqAMwzCMnqQrAkpE9heReSJyn4jcJSJvDI6dJiILROQREZkahO8rIg9kx84PwjcSkauz8DtEZOfg\n2Eki8mj2+0AQPlFE7syuuUpERnXiuQ3DMIzydEuD+grw76r6BuCMbB8RmQwcC0wGpgEXiohflfEi\nYLqqTgImici0LHw6sDwLPw84N4trfBb3/tnvTBHZPLvmXOBr2TUrsjiMAubOndvtJPQMlhcVLC8q\nWF60nm4JqCcALyy2AJZk20cAV6rqWlVdBCwEDhCR7YCxqjovO+9y4Mhs+3BgVrZ9LXBItv1OYLaq\nrlTVlcAc4LBM4B0M/DA7b1YQl5GDFb4KlhcVLC8qWF60nr4u3ffzwK9F5Ks4IXlgFr49cEdw3mJg\nB2Bttu1ZkoWT/T8OoKrrROQ5Edkyi2txIq7xwEpV3ZCIyzAMw+gR2iagRGQOsG3i0BeATwKfVNUf\ni8gxwKXAlHalJUA7cA/DMAyjFahqx3/AqmBbgOey7c8Dnw+O3QwcgBN0Dwfh7wcuCs55U7bdBzyd\nbR8H/Hdwzbdx/VsCPA2MyMIPBG7OSafaz372s5/96v+1QlZ0y8S3UEQOUtXbgXcAj2bhNwA/EJH/\nwpndJgHzVFVFZJWIHADMA04EvhFccxLONHg0cEsWPhs4W0S2wAmlKcDnsrhuA44Brs6uvS6VSFWV\nVLhhGIbRfiTTFDp7U5H9gG8BGwEvAqeo6n3ZsdOBDwPrgFNV9edZ+L7A94BNgJtU9ZNZ+EbAFcAb\ngOXAcZmDBSLyIeD07LZnqeqsLHwicBWuP+pe4ARVXdvepzYMwzDqoSsCyjAMwzBqYTNJJBCRadlA\n4QUi8rlup6fdiMgEEblNRB4SkQdFxGun40VkTjbQeXZmLvXXJAdUDxVEZGQ2kPzGbH9Y5oWIbCEi\nPxSRh0VkvogcMIzz4rSsjDwgIj/IJgkYFnkhIpeKyJMi8kAQVvez5024kEs3nCR6+QeMxI2/2gUY\nBdwP7NXtdLX5mbcFXp9tjwH+AOyFG0D92Sz8c8A52fbkLF9GZfm0kMzpZKj8gE8D3wduyPaHZV7g\nxgl+ONvuw41fHHZ5kT3PY8BG2b7vvx4WeQG8DdeN8kAQVs+ze2vdPGD/bPsmYFrRfU2DGsj+wEJV\nXaSuX+oq3ADiIYuqLlPV+7PtF4CHcU4q4SDocEBzakD1/h1NdBsRkR2BdwEX4xxsYBjmRTbzyttU\n9VIAVV2nqs8xDPMCWIUbj7mpiPQBmwJLGSZ5oaq/ws26E1LPs9eacCGJCaiBvDzwN8MP8B0WiMgu\nuJbSncA2qvpkduhJYJtsO28Q9FDhPOAzwIYgbDjmxUTgaRG5TETuFZHvishmDMO8UNVnga8Bf8EJ\nppWqOodhmBcB9T57HF5zkgQTUAMZtl4jIjIGN13Uqar6fHhMnU5elDdDIt9E5N3AU+q8SpPDDIZL\nXuBMevsAF6rqPsBfcWMVX2a45IWI7AZ8Cmey2h4YIyInhOcMl7xIUeLZG8IE1ECWABOC/QlUS/0h\nSTaj+7XAFarqx4U9KSLbZse3A57KwuM82pHKfIqDnTcDh4vIn4ArgXeIyBUMz7xYDCxW1buy/R/i\nBNayYZgX+wG/UdXlqroO+BFukP9wzAtPPWVicRa+YxRemCcmoAZyN2629F1EZDRu9okbupymtpJN\noHsJMF9Vvx4c8oOgoXpA8w3AcSIyOhtTNgnX+TnoUdXTVXWCqk7EzUZyq6qeyPDMi2XA4yKyRxZ0\nKPAQcCPDLC+AR4A3icgmWXk5FJjP8MwLT11lIvueVmWeoIKbcCE5ScLLdNs7pBd/wGE4T7aFwGnd\nTk8HnvetuP6W+4H7st803EDmX+Bm+pgNbBFcc3qWP48A7+z2M7QpXw6i4sU3LPMCeB1wF/A7nNaw\n+TDOi8/iBPQDOKeAUcMlL3DWhKXA33B99B9q5NmBfbP8Wwh8o9Z9baCuYRiG0ZOYic8wDMPoSUxA\nGYZhGD2JCSjDMAyjJzEBZRiGYfQkJqAMwzCMnsQElGEYhtGTmIAyjCYRkS2zpTnuE5EnRGRxtv28\niHyzTff8uIh8MBG+S7gkQgvus5GI/FJErK4wOk63lnw3jCGDqi7HTbCLiJwJPK+q/9Wu+2Wj8KcD\nb2zXPTyqukZEfoWbdfpH7b6fYYRYq8gwWo8AiEh/sODhDBGZlWkji0TkKBH5qoj8XkR+li3h4Bd0\nmysid4vIzX6us4i3AI+omxPOX/M7EbkfOOXlRDht6v+3d+cgUmVRGMf/n8K0GIzJCEYDgijYjAsK\nKmNiYKRooDgM4hYpGBkbKxi6JiLtiokoCo6IwoA0Loi2SwsKgkzkEgkKMyrtN8G9NbwpW7tpuosS\nvh8U9bjvcl9VQXHq3Hqcc0PSvfpYVsdPSFrbmHdG0hpJvZLu1OzvoaRZdcol4PcJ+JwivikBKqJz\nZgIrKH10TgPXbM8D/gZW1YK9B4F1thcDfcCeYdZZTqkZ2dIH7LS9oG3ea2Cl7UWUuoIH6vgxYCv8\n1/NpGXAZ2AHst72QUpKmVST5AaWIbkRHZYsvojMMXLE9JGmQ0l31aj33mNLGYTbQC1wvu3hMptQ/\na/cz0A+lJTswzXZ/PXeKUksS4AfgkKT5wFBdH9s3JB2R9BOwHjhXX9dNYHdt2Hje9vM6/4OkSZKm\n2P5nvD6QiJEkQEV0zkcA258lfWqMf6Z8FwU8sT2abGXYXlVt47uAl7Y3SZoMNIPLSUo16d+o2ZTt\ns5JuA6uBPyRtt/1nY90U7oyOyhZfRGd8LaA0PQOmS1oKpUeXpLnDzPsLmAFg+y3wVtKv9dzGxrwf\ngVf1eDMlI2s5TmnAZ9tP6/Vm2n5h+yBwEfiljvcAQ7Y/jOI9RIybBKiI8efG83DH8GU2YtufKFtu\n++oNDwOU/4fa9VMa6LVsAw5LGmhb+wiwpa41B3jfuNgbSj+jvsY6GyQN1nV6KVkWlDsUb3397UZM\njLTbiPjO1NvM7wNLbH8c4xpTgUfAQtvvRpi7F7hr+8JYrhUxVsmgIr4zLr8qj/L/7bxRk9TqBntg\nFMGph3LX4Lc7n0ZMgGRQERHRlZJBRUREV0qAioiIrpQAFRERXSkBKiIiulICVEREdKUEqIiI6Er/\nAroXw+alw6rIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa54026a850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(1,1000,num=1000)\n",
    "plt.plot(x, f, label=\"Cash flow\")\n",
    "plt.xlabel(\"Time (days)\")\n",
    "plt.ylabel(\"Amount\")\n",
    "plt.title(\"Cash flow\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram and best fit normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fa515035690>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEZCAYAAAAZnxsyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2clXWd//HXW+40RUfUuJPEciipMLSE7nS8XaQEy1Zw\nTcHYVZcwWWsTtnaD3a1Faw3NJH+loaUoqctOhTdoO6mljHiDN4BCNqZj4C1imgry+f1xfQcOhzNn\nDjBn5pzh/Xw8zmOu872+3+/1uQ7MfM51Xd/reykiMDMzq0a7dHYAZmZm28tJzMzMqpaTmJmZVS0n\nMTMzq1pOYmZmVrWcxMzMrGo5iVmXI+kxSUd0dhydSdLnJD0j6TVJh3TA9uZK+o9tqN8S3zpJH5HU\nJOmYcsZoXZOTmFWVQn/sJE2UdHfL+4j4UETc1UY/gyVtlNRVfwe+B0yOiN4RsbQDthfpVaqW+PaM\niIe3o70Z4CRm1ae9/9ipHfva3KnUrRz9lrhtAe8BlnVWDMVUenxWXZzErCvYIqmlo7Wj0/LhkpZI\nelXSaknfS9VajtTWplNuI5T5Zmq/RtLVkvbM6fcMSU9LejGnXst2Zki6UdLPJL0KTJD0MUn3SnpF\n0nOSfiCpR05/GyX9o6SV6bTav0t6X2qzVtL1ufXz9rFgrJJ6Aa8B3YClkla20v6DkhZJeil9LtNz\nPq9iMX8/be9VSY9IGprTbR9Jv0r7cp+k9xbYbpvxSeolabak5vT6vqSead1vJX0+LX8yfYaj0/tj\nJD1UaH+t63ISs2qUf/SU/z43qV0CfD8i9gLeC/wilX86/dwrnXJbDJwJTADqUt09gMsA0h/rHwKn\nAv2BvYABedsdA/wibes64B3gPGAf4OPAMcDkvDbHA8OBkcAFwI/TNt4DfDgtF1Iw1oh4KyL2SHWG\nRURtfkNJvYE7gIVpXw4C7kyrN7QWs6S/SZ9bbdrHvwVebukWGA/MAPYGVgHfzt92KfEB3wAOBw5J\nr8OBb6Z1DWmfAY4EngKOyHnfUKA/68KcxKzaCFiQjhRekfQKWXJp7RTj20CtpH0j4o2UrFr6yXca\n8N8R0RQRrwPTgfHp1OAXgPqI+H1ErAf+rcA2fx8R9QAR8WZEPBgRjRGxMSKeBv4f2R/aXBdFxF8i\nYhnwKHBL2v464BayBFdIa7GW8jv9WeC5iPh+RLydtt+Y4i4W83qgN3CwpF0i4omIWJ3WBXBzRCyJ\niHeAa4GPlBBLIX8H/HtEvBgRLwIzgdPTurty4vk08F85748Efrud27Qq5SRm1SaAsRGxd8uL7Eih\ntWtbk4AhwHJJjZI+U6Tv/sDTOe//BHQH+qZ1z24KIuKvwEt57Z/NfSNpSDq99ud0ivHbZEc4udbk\nLP+1wPs9KKxYrG0ZRHYEs5ViMUfEb8iOTH8IrJF0RTqqa21fWou9LQPYet9ajnrvBYZIejdZkrwG\nGCRpH+BjbD5NbDsJJzHrClodnBERqyLi7yJiP+BC4EZJu1H4yO05YHDO+/eQnV5bDfwZ2H/TBrM+\n8hNSfp9zyAYvHJROv32D9vuday3WNQVrb+lPZKcgCykac0T8ICI+Cgwl+3Lwz9scedsK7dtzaftv\nAA8AU4FH01Hx74GvAqsi4mVsp+IkZl2apC9K2i+9fZUs0WwEXkg/35dTfR7wT8qG3+8BfAe4PiI2\nAjcBJ0r6eBpkMIO2RzbuQTaI4Q1JHwD+sZSQW1nOVyzWtvwK6C/pvDSIorekw4vEHACSPpoGwPQA\n3gDeJLvu11as22oe8E1J+0ral+zU7c9y1v8W+DKbTx02AFPwqcSdkpOYdQXFht3/DfCYpNeA7wPj\n0+CCN8hOlf0uXVs7HLiK7I/lXWSn294AzgWIiMfT8vVkRwWvAc8DbxWJ4Wtk13fWkV1buj6vTqGY\n89e3tl+txlqkb9K+/AU4DjiR7AjzSTYPligUc4s9U9nLQBPwIvDdIrEWuxWi2Lr/BJYAj6TXklTW\n4rdkybbl1OFdwO74VOJOSeV8KKakUcBssuG0P4mICwvUuRQ4geyXcGJEPFSsraQ+wA3AAWS/SKdE\nxNq0bjrwJbJvh1+JiNtT+WHAXGBXYGFEnJfKDyD7Y7Av2S/mFyOiud0/COty0tHPK2Sn3Z5uq76Z\nlUfZjsTSiK7LgFFk589PlXRwXp3RZH8EaoGzyM7Ht9V2GrAoIoaQDQueltoMBcal+qOAyyW1nOKY\nA0xK26lNCRKyWQPmRsQhwL+TjXQyK0jSiZLeJWl3sv87jziBmXWucp5OPJzsQmtTuvh6PTA2r84Y\n4GqANPS5RlK/NtpuapN+npSWxwLzImJ9RDSR3acyQlJ/oHfLEGKy0UwtbQ4GfpOWGwrEZ5ZrDNCc\nXu8juy/KzDpROZPYQOCZnPfPprJS6gwo0rZvRLSMwFrD5iHFA9hyiHNuX7nlzTl9LQVOTsufA3pL\n2rutHbOdU0T8QxrWXxMRx0VEwdkwzKzjlDOJlXqxrZRRTSrUX2QX9Hbkot7XgCMlPUh2138zm0db\nmZlZhetexr6byW6qbDGIvJtBC9TZP9XpUaC8ZcDFGkn9ImJ1OlX4fBt9NZNzf09uXxHxZ9KRWLpQ\nf3KaKWELkjy7tpnZdoiIskyy3aKcR2JLyAZRDE731YwD6vPq1ANnAEgaCaxNpwqLta0nmzOO9HNB\nTvl4ST0lHQjUAo1pWpx16f4WkU1fsyBtc5+caXqmA1e2tjMRUfGvb33rW50eQ1eJsxpidJyOs9Jf\nHaFsR2IRsUHSFOA2smHyV0bEcklnp/VXRMRCSaMlrQJeJ5vUtNW2qetZwHxJk0hD7FObZZLmk802\nsIHsWUUtn+JksiH2u5ENsb81ldcB/5WOtFpuoDQzsypRztOJRMQtZJOY5pZdkfd+SqltU/nLwLGt\ntPkO2cwF+eUPkM0Inl9+E9lMDGZmVoU8Y0cXUldX19khlKQa4qyGGMFxtjfHWX3KOmNHVyEp/DmZ\nmW0bSUQVD+wwMzMrKycxMzOrWk5iZmZWtZzEzMysajmJmZlZ1XISMzOzquUkZmZmVctJzMzMqpaT\nmJmZVS0nMTMzq1pOYmZmVrWcxMzMrGqV9VEsZlYeU6dNZe2ba0uqW7NrDbNnzS5zRGadw0nMrAqt\nfXMtg08aXFLdpgVNZY3FrDOV9XSipFGSVkhaKemCVupcmtYvlTS8rbaS+khaJOlJSbdLqslZNz3V\nXyHp+JzywyQ9mtZdklN+kKS7JT2Utn9C+38KZmZWLmVLYpK6AZcBo4ChwKmSDs6rMxo4KCJqgbOA\nOSW0nQYsioghwJ3pPZKGAuNS/VHA5ZJanmMzB5iUtlMraVQq/ybw84gYDowHLm/fT8HMzMqpnEdi\nhwOrIqIpItYD1wNj8+qMAa4GiIjFQI2kfm203dQm/TwpLY8F5kXE+ohoAlYBIyT1B3pHRGOqd01O\nmz8De6XlGqB5x3fbzMw6SjmviQ0Ensl5/ywwooQ6A4EBRdr2jYg1aXkN0DctDwDuK9DX+rTcojmV\nA/wXcK+kc4HdgWNK2TEzM6sM5UxiUWK9Uh5drUL9RURIKnU7hVwM/CQivi9pJPBz4IOFKs6YMWPT\ncl1dHXV1dTuwWTOzrqehoYGGhoYO3WY5k1gzMCjn/SC2PCIqVGf/VKdHgfKWU31rJPWLiNXpVOHz\nbfTVnJbzywE+AXwLICLuk7SrpH0j4sX8nclNYmZmtrX8L/gzZ84s+zbLeU1sCdkgisGSepINuqjP\nq1MPnAGQjoTWplOFxdrWAxPS8gRgQU75eEk9JR0I1AKNEbEaWCdpRBrocTrwv6nNCuDYtP2DgV0L\nJTAzM6tMZTsSi4gNkqYAtwHdgCsjYrmks9P6KyJioaTRklYBrwNnFmubup4FzJc0CWgCTkltlkma\nDywDNgCTI6LlVONkYC6wG7AwIm5N5f8MXCnpn8hOV7YkRzMzqwLa/HfeWiMp/DlZJZk4deI23ew8\nd/bcssZjVogkIqKUcQ/bzXMnmplZ1XISMzOzquUkZmZmVctJzMzMqpaTmJmZVS0nMTMzq1pOYmZm\nVrWcxMzMrGo5iZmZWdVyEjMzs6rlJGZmZlXLSczMzKqWk5iZmVUtJzEzM6taTmJmZla1nMTMzKxq\nlTWJSRolaYWklZIuaKXOpWn9UknD22orqY+kRZKelHS7pJqcddNT/RWSjs8pP0zSo2ndJTnlF0t6\nKL2ekPRK+38KZmZWLmVLYpK6AZcBo4ChwKmSDs6rMxo4KCJqgbOAOSW0nQYsioghwJ3pPZKGAuNS\n/VHA5ZJanig6B5iUtlMraRRARJwfEcMjYjjwA+Cm9v8kzMysXMp5JHY4sCoimiJiPXA9MDavzhjg\naoCIWAzUSOrXRttNbdLPk9LyWGBeRKyPiCZgFTBCUn+gd0Q0pnrX5LTJ9XfAvB3ZYTMz61jlTGID\ngWdy3j+bykqpM6BI274RsSYtrwH6puUBqV6hvnLLm/PjkHQAMBj4TRv7ZGZmFaR7GfuOEuup7Sqo\nUH8REZJK3U4x44FfRESrfc2YMWPTcl1dHXV1de2wWTOzrqOhoYGGhoYO3WY5k1gzMCjn/SC2PCIq\nVGf/VKdHgfLmtLxGUr+IWJ1OFT7fRl/NablQXy3GAZOL7UxuEjMzs63lf8GfOXNm2bdZztOJS8gG\nUQyW1JMsUdTn1akHzgCQNBJYm04VFmtbD0xIyxOABTnl4yX1lHQgUAs0RsRqYJ2kEWmgx+k5bZD0\nAWDviLivPXfezMzKr2xHYhGxQdIU4DagG3BlRCyXdHZaf0VELJQ0WtIq4HXgzGJtU9ezgPmSJgFN\nwCmpzTJJ84FlwAZgcs7pwcnAXGA3YGFE3JoT6jg8oMPMrCqpyGUgSyQVu1xm1uEmTp3I4JMGl1S3\naUETc2fPLWs8ZoVIIiJKGfew3Txjh5mZVS0nMTMzq1pOYmZmVrWcxMzMrGo5iZmZWdVyEjMzs6rl\nJGZmZlXLSczMzKqWk5iZmVUtJzEzM6taTmJmZla1nMTMzKxqOYmZmVnVchIzM7Oq5SRmZmZVy0nM\nzMyqVlmTmKRRklZIWinpglbqXJrWL5U0vK22kvpIWiTpSUm3S6rJWTc91V8h6fic8sMkPZrWXZK3\n/VMkPS7pMUnXtu8nYGZm5VS2JCapG3AZMAoYCpwq6eC8OqOBgyKiFjgLmFNC22nAoogYAtyZ3iNp\nKDAu1R8FXC6p5Ymic4BJaTu1kkalNrWp/Sci4kPAee3+QZiZWdmU80jscGBVRDRFxHrgemBsXp0x\nwNUAEbEYqJHUr422m9qknyel5bHAvIhYHxFNwCpghKT+QO+IaEz1rslp8w/AZRHxaorhxfbZdTMz\n6wjlTGIDgWdy3j+bykqpM6BI274RsSYtrwH6puUBqV6hvnLLm3P6qgXeL+keSfdK+pvSds3MzCpB\n9zL2HSXWU9tVUKH+IiIklbqdQnoABwFHAoOAuyR9uOXILNeMGTM2LdfV1VFXV7cDmzUz63oaGhpo\naGjo0G2WM4k1kyWGFoPY8oioUJ39U50eBcqb0/IaSf0iYnU6Vfh8G301p+X8csiO9hZHxDtAk6Qn\nyZLaA/k7k5vEzMxsa/lf8GfOnFn2bZbzdOISskEUgyX1JBt0UZ9Xpx44A0DSSGBtOlVYrG09MCEt\nTwAW5JSPl9RT0oFkpwobI2I1sE7SiDTQ43Tgf1ObBUBd2v6+wBDgqfb6AMzMrLzKdiQWERskTQFu\nA7oBV0bEcklnp/VXRMRCSaMlrQJeB84s1jZ1PQuYL2kS0AScktoskzQfWAZsACZHRMupxsnAXGA3\nYGFE3Jra3CbpeEmPA+8AX4uIV8r1mZiZWfvS5r/z1hpJ4c/JKsnEqRMZfNLgkuo2LWhi7uy5ZY3H\nrBBJREQp4x62m2fsMDOzquUkZmZmVctJzMzMqpaTmJmZVS0nMTMzq1pOYmZmVrWcxMzMrGo5iZmZ\nWdVqM4lJGiPJyc7MzCpOKclpHLBK0kWSPlDugMzMzErVZhKLiNOA4WQT485Nz906S1LvskdnZmZW\nREmnCdPztW4EbiB7yOTngIckfaWMsZmZmRVVyjWxsZL+B2gge87XxyLiBGAYcH55wzMzM2tdKY9i\n+Tzw/Yi4K7cwIt6Q9PflCcvMzKxtpZxOXJOfwCRdCBARd5QlKjMzsxKUksSOK1A2upTOJY2StELS\nSkkXtFLn0rR+qaThbbWV1EfSIklPSrpdUk3Ouump/gpJx+eUHybp0bTukpzyiZJekPRQen2plP0y\nM7PK0GoSk/SPkh4F3p8SQMurCXikrY4ldQMuA0YBQ4FTJR2cV2c0cFBE1AJnAXNKaDsNWBQRQ4A7\n03skDSW7HWBoane5pJaHsc0BJqXt1EoalcoDmBcRw9Prqrb2y8zMKkexI7HrgBOBeuCzaflE4LA0\n7L4thwOrIqIpItYD1wNj8+qMAa4GiIjFQI2kfm203dQm/TwpLY8lS0jrI6IJWAWMkNQf6B0Rjane\nNTltlF5mZlaFiiWxSMngy8BrwLr0Ckl9Suh7IPBMzvtnU1kpdQYUads3Itak5TVA37Q8INUr1Fdu\neXNOXwGcLOkRSb+QtH8J+2VmZhWi2OjEecBngAfI/tjnO7CNvgu1KaSUIyEV6i8iQlKp2ynkl8B1\nEbFe0llkR3bH7EB/ZmbWgVpNYhHxmfRz8Hb23QwMynk/iC2PiArV2T/V6VGgvDktr5HULyJWp1OF\nz7fRV3Na3qqviHg5p/xK4KLWdmbGjBmbluvq6qirq2utqpnZTqmhoYGGhoYO3WarSUzSocUaRsSD\nbfS9hGwQxWDgObJBF6fm1akHpgDXSxoJrI2INZJeKtK2HpgAXJh+Lsgpv07SxWSnC2uBxnS0tk7S\nCKAROB24NO1jv4hYndqPAZa1tjO5SczMzLaW/wV/5syZZd9msdOJF1P8lOBRxTqOiA2SpgC3Ad2A\nKyNiuaSz0/orImKhpNGSVgGvA2cWa5u6ngXMlzQJaAJOSW2WSZpPlog2AJMjoiX+ycBcYDdgYUTc\nmsq/ImlMqv8SMLHYPpmZWWXR5r/z1hpJ4c/JKsnEqRMZfNLgkuo2LWhi7uy5ZY3HrBBJRERZR4AX\nO514dET8RtLJFB5UcXM5AzMzM2tLsdOJRwK/Ibs3rNBhiJOYmZl1qmKjE7+Vfk7ssGjMzMy2QSmP\nYtlX0g/S3IIPSrpE0j4dEZyZmVkxpUwAfD3ZvVifB74AvED2cEwzM7NOVcrzxPpFxH/kvP9PSePK\nFZCZmVmpSjkSu13SqZJ2Sa9xwO3lDszMzKwtxYbY/4XNoxKnAj9Ly7uQ3Zj81fKGZmZmVlyx0Yl7\ndGQgZmZm26qUa2JI2ptsLsJdW8oi4q5yBWVmZlaKNpOYpH8AvkI2Q/xDwEjgXuDo8oZmZmZWXCkD\nO84je9JyU0QcBQwHXi1rVGZmZiUoJYm9GRF/BZC0a0SsAN5f3rDMzMzaVso1sWfSNbEFwCJJr5A9\nAsXMzKxTtZnEIuJzaXGGpAZgT+DW1luYmZl1jFJHJx4GfIrsvrF7IuLtskZlZmZWglImAP43sqci\n9wH2BX4q6V9L6VzSKEkrJK2UdEErdS5N65dKGt5WW0l9JC2S9KSk2yXV5KybnuqvkHR8Tvlhkh5N\n6y4pEMPJkjZKOrSU/TIzs8pQysCOLwIfi4hvRcS/kQ2xP72tRpK6AZcBo4ChwKmSDs6rMxo4KCJq\ngbOAOSW0nQYsioghwJ3pPZKGAuNS/VHA5ZJanig6B5iUtlMraVRODL3JRmDeV8JnYWZmFaSUJNYM\n7Jbzflfg2RLaHQ6sioimiFhPNhv+2Lw6Y4CrASJiMVAjqV8bbTe1ST9PSstjgXkRsT4imoBVwAhJ\n/YHeEdGY6l2T0wbgP4BZwFtAWR+jbWZm7avVJJaeIfYDsnvCHpc0V9Jc4DFKu09sIPBMzvtnU1kp\ndQYUads3Itak5TVA37Q8gC2Ta25fueXNLX2l04cDI2JhWlfoCdZmZlahig3seIDsj/oSsuH1LX/g\nGyjtj32pCaGUox8V6i8iQtJ2JZ50qvFiYMI2xmJmZhWi2ATAc1uWJfUChqS3K9IpvrY0k01V1WIQ\nW5+GzK+zf6rTo0B5c1peI6lfRKxOpwqfb6Ov5rScX94b+CDQkC6d9QPqJZ0YEQ/m78yMGTM2LdfV\n1VFXV1don83MdloNDQ00NDR06DYVUfxARlId2bWnp1PRe4AJEfHbNtp1B54AjgGeAxqBUyNieU6d\n0cCUiBgtaSQwOyJGFmsr6SLgpYi4UNI0oCYipqWBHdeRXU8bCNxBNmgkJC0mm/+xEfg1cGlEbHGv\nm6T/A75aKIFJirY+J7OONHHqRAafNLikuk0Lmpg7e25Z4zErRBIRUdYzXKXcJ3YxcHxEPJGCGkI2\n0KLocPSI2CBpCnAb0A24MiWhs9P6KyJioaTRklaRPaPszGJtU9ezgPmSJpHNHHJKarNM0nxgGbAB\nmJyTeSaT3SawG7AwP4GZmVl1KuVI7JGIGNZWWVfmIzGrND4Ss2pQKUdiD0j6CfBzsoEPp5EN9jAz\nM+tUpSSxc4ApZNeUAO4GLi9bRGZmZiUqmsTSAIulEfEB4L87JiQzM7PSFJ2xIyI2AE9IOqCD4jEz\nMytZKacT+5DN2NFINoIQsvuMx5QvLDMzs7aVksS+mX7mjjDxUD0zM+t0rSYxSbuRDeo4CHgEuKrE\nmTrMzMw6RLFrYlcDh5ElsNHA9zokIjMzsxIVO514cER8GEDSlcD9HROSmZlZaYodiW1oWUijFM3M\nzCpKsSOxYZJey3m/W877iIg9yxiXmZlZm4o9iqVbRwZiZma2rYre7GxmZlbJnMTMzKxqOYmZmVnV\nchIzM7OqVdYkJmmUpBWSVkq6oJU6l6b1SyUNb6utpD6SFkl6UtLtkmpy1k1P9VdIOj6n/DBJj6Z1\nl+SUnyPpEUkPSbpX0iHt/ymYmVm5lC2JSeoGXAaMAoYCp0o6OK/OaOCgiKgFzgLmlNB2GrAoIoYA\nd6b3SBoKjEv1RwGXS2qZ73EOMCltp1bSqFR+bUQMi4jhwHfw42bMzKpKOY/EDgdWRURTmnPxemBs\nXp0xZNNbERGLgRpJ/dpou6lN+nlSWh4LzIuI9RHRBKwCRkjqD/SOiMZU75qWNhGRex/cHsCLO77b\nZh1jlw3v0ONNT2dqO7dSZrHfXgOBZ3LePwuMKKHOQGBAkbZ9I2JNWl4D9E3LA4D7CvS1Pi23aE7l\nAEiaDJwP7A58ooT9MqsIfZ96ni+dexUbenXn1f32Yt1+e7Lu3XvS/P4BPPSZQzs7PLMOUc4kVurj\nWtR2FVSov4gISTv0WJiIuJzs1OOpwFXAUYXqzZgxY9NyXV0ddXV1O7JZsx325yH9+fat/8K7Xn2D\nPV9Yx17Pr2PPF9YR3TxeyzpHQ0MDDQ0NHbrNciaxZmBQzvtBbHlEVKjO/qlOjwLlzWl5jaR+EbE6\nnSp8vo2+mtNyob5y3QD8qLWdyU1iZh1u3TrYs8BMbxJv1OzOGzW7s7q2f8fHZZYj/wv+zJkzy77N\ncn5lW0I2iGKwpJ5kgy7q8+rUA2cASBoJrE2nCou1rQcmpOUJwIKc8vGSeko6EKgFGiNiNbBO0og0\n0OP0ljaSDsqJ5TNkj50xqyyNjfCBD0Bzoe9eJbr7bvjGN9ovJrMKUbYjsYjYIGkKcBvQDbgyIpZL\nOjutvyIiFkoaLWkV8DpwZrG2qetZwHxJk4Am4JTUZpmk+cAyshn4J0dEy6nGycBcYDdgYUTcmsqn\nSDqW7LrZCy3bN6sYv/0tfOELcNVVMHBg2/Vb86EPwTnnwD77wPnnt198Zp1Mm//OW2skhT8n63C3\n3gpnnAHz5sExx2yxauLUiQw+aXBJ3TQtaGLu7Lnwpz/BJz8J3/0ujB/f/vGa5ZFERJQy7mG7lfOa\nmJltrwUL4Oyz4X//Fz7+8fbp8z3vgV//Go49Fvr2haMKjmEyqyoexmRWiQYMyI7E2iuBtRg2DG64\nASZOhDfeaN++zTqBj8TMKsTUaVNZ++baLQuvLly38YHGkk8nbuWoo+Dhh+Fd79q+9mYVxEnMrEKs\nfXNtyYnpnsZ7dmxje++9Y+3NKoRPJ5qZWdXykZhZZ3vuOVi4sLOjMKtKPhIz62znnQdNTZ0bw513\nwle/Cr6VxKqMk5hZZ/rVr7JBFp09m8ZHPwp33AGXX965cZhtIycxs87y+uswZQrMmQO77da5sey1\nF8yfD9/6VnZ606xK+JqYWWeZMQM+/ens5uMyalzcyMSpE0uqe1rtgRz3ta/BddeVNSaz9uIkZtYZ\nNmyAJ56AK68s+6be1tslD93/xfr1HPfr38FvfgNHH13ewMzagZOYWWfo3h3q8x/q0Pne7tEDrrkG\n9tuvs0MxK4mTmJlt6YgjOjsCs5J5YIeZmVUtJzEzM6taTmJmHeX+++H55zs7CrMupexJTNIoSSsk\nrZR0QSt1Lk3rl0oa3lZbSX0kLZL0pKTbJdXkrJue6q+QdHxO+WGSHk3rLskpP1/S42nbd0h6T/t/\nCrbTe+01OPlkeOyxzo5k27z1VjajyFtvdXYkZgWVNYlJ6gZcBowChgKnSjo4r85o4KCIqAXOAuaU\n0HYasCgihgB3pvdIGgqMS/VHAZdLanmq6BxgUtpOraRRqfxB4LCIOAS4EbiofT8FM7KbiI8+uvqG\nrffqBU8/nT0N2qwClftI7HBgVUQ0RcR64HpgbF6dMaSnJkXEYqBGUr822m5qk36elJbHAvMiYn1E\nNAGrgBGS+gO9I6Ix1bumpU1ENETEm6l8MbB/++y6WfLgg3DttfC973V2JNvnkktg9mz44x87OxKz\nrZQ7iQ0Ensl5/2wqK6XOgCJt+0bEmrS8BuiblgekeoX6yi1vLhAHwCTA04lb+3nnHTjrLLjwQth3\n386OZvsccACcf352WtGswpT7PrFSp8RW21VQof4iIiTt8NTbkr4IHAr8U6H1M2bM2LRcV1dHXV3d\njm7Sdga/+x3U1MCECZ0dyY756ldh2DD45S/hxBM7OxqrUA0NDTQ0NHToNsudxJqBQTnvB7HlEVGh\nOvunOj0JbnTlAAAUf0lEQVQKlDen5TWS+kXE6nSqsGXIV2t9NbPlacLcvpB0LPAvwBHp1OVWcpOY\nWcmOOAJuvRVUyve0CtarF/zwh7B0aWdHYhUs/wv+zJkzy77Ncp9OXEI2iGKwpJ5kgy7y59qpB84A\nkDQSWJtOFRZrWw+0fLWdACzIKR8vqaekA4FaoDEiVgPrJI1IAz1Ob2mTRkP+CDgxIl5s5/03y6aY\n6gqOPTY7IjOrIGX97YqIDZKmALcB3YArI2K5pLPT+isiYqGk0ZJWAa8DZxZrm7qeBcyXNAloAk5J\nbZZJmg8sAzYAkyM2PeVvMjAX2A1YGBG3pvKLgN2BG9NAxqcjomWgiJmZVbCyf0WMiFuAW/LKrsh7\nP6XUtqn8ZaDg8ysi4jvAdwqUPwB8uED5cUXCNzOzCuYZO8zamx8qadZhnMTM2tMf/wiHHAIvv9zZ\nkZTfnXfCunWdHYXt5JzEzNrTV76SDX7o06ezIym/efPgm9/s7ChsJ+ckZtZe6uth5crsxuCdwUUX\nwfz52cTGZp3EScysPbzxRjajxQ9/CD17dnY0HaNPn2wqrbPOgg0bOjsa20k5iZm1h4svhhEj4Jhj\nOjuSjnXaabDPPnDppZ0die2knMTM2sOUKfCDH3R2FB1Pgssvh5/9zEdj1im6yFQCZp2spqbtOl3V\nkCGwZAl069bZkdhOyEdiZrbjnMCskziJmZlZ1XISM9tescNPADKzHeQkZrY9Hnssm9XdiWxrb70F\nzzzTdj2zduCBHWbbKgK+/GU45ZTqf05YnsbFjUycOrGkujW71jB71uytV9TXw6xZsHhx13kMjVUs\n/w8z21Y//zn85S9wzjmdHUm7e1tvM/ikwSXVbVrQVHjFF74Ac+ZkN36fd167xWZWiJOY2bZYuxa+\n/nVYsMAj8lojZUnsk5+Ez38eBg1qu43ZdvI1MbNSRcCXvpQdaYwY0dnRVLb3vx/OPTebENmsjMp+\nJCZpFDCb7OnMP4mICwvUuRQ4AXgDmBgRDxVrK6kPcANwAOnJzhGxNq2bDnwJeAf4SkTcnsoPI3uy\n865kT3Y+L5UfkbbxYWB8RNzU/p+CdQkSnHEGnHBCyU2mTpvK2jfXllS38YHGkk/lVYVp02DYMLjn\nHvjUpzo7GuuiyprEJHUDLiN7CnMzcL+k+ohYnlNnNHBQRNRKGgHMAUa20XYasCgiLpJ0QXo/TdJQ\nYBwwFBgI3CGpNiIi9TspIholLZQ0KiJuBZ4GJgBfK+dnYV3ESSdtU/W1b64tOTHd03jPdgRUwXr1\nyhLYvvt2diTWhZX7dOLhwKqIaIqI9cD1wNi8OmOAqwEiYjFQI6lfG203tUk/W/6yjAXmRcT6iGgC\nVgEjJPUHekdEY6p3TUubiHg6Ih4FNrbjfpsZwH77dbkRnFZZyp3EBgK5N4w8m8pKqTOgSNu+EbEm\nLa8B+qblAaleob5yy5sLxGFmZlWm3NfESr0TtJSvairUX0SEpLLfcTpjxoxNy3V1ddTV1ZV7k1YJ\nmpthoL/vmJWioaGBhoaGDt1muZNYM5A7vnYQWx4RFaqzf6rTo0B5c1peI6lfRKxOpwqfb6Ov5rRc\nqK9crSbD3CRmO4mrr4bZs+HBB31KrL0sXw7r1nl0ZxeV/wV/5syZZd9muU8nLgFqJQ2W1JNs0EV9\nXp164AwASSOBtelUYbG29WSDMUg/F+SUj5fUU9KBQC3QGBGrgXWSRkgScHpOmxaitCNC2xk8/DB8\n7WvZjc1OYO2nqSkbHNPU1NmRWBdR1iQWERuAKcBtwDLghohYLulsSWenOguBpyStAq4AJhdrm7qe\nBRwn6Ung6PSeiFgGzE/1bwEmp5GJpH5/AqwkGzByK4Ckj0l6BvgCcIWkR8v2gVh1eOUVOPnk7CGX\nH/xgZ0fTtZxwAkyfDp/9bHZEZraDyn6fWETcQpZQcsuuyHs/pdS2qfxlsqH3hdp8B/hOgfIHyO4F\nyy+/ny1PQdrObONG+OIXYcwYGD++s6Ppms49F1asyD7f+nrPr2g7xDN2mOV67DF45x246KLOjqTr\nkuDSS7PP+fzzOzsaq3JOYma5hg2DW26BHj06O5KurXt3uOGGbH5Fsx3gJGaWzwM5OkZNDYwb19lR\nWJVzEjMzs6rlJGY7t9de6+wIzGwHeFiQ7bxuuQUmTYL77/esHNuhXZ4Cne/xx+F974Ndd92x4Gyn\n4SRmO5+NG+Hb34Yf/Qjmz3cC207t8hTofD/4ASxenN1k7nv0rAROYrZzefVVOP10ePllWLIE+vff\n5i526meElducOXDllVBXB9/8ZnZP2S6+6mGtcxKznUdENlPERz4CN94IPXtuVzc79TPCyk2Cv//7\nLImdfjr86lcwd66Plq1VTmK285Dgppvg3e/u7EisLQcdBHffDbNmwcqVTmLWKicx27m0ksB8irAC\nde+enVI0K8JJzAyfIiy3soxkNMNJzLqi116DefO4ef511H9ocElNfHRVXu0+kvHXv4Zjj4VevXYo\nLqt+TmLWNURAYyP8+MfZda+jjuLJml19dNUVbdyYDfaYPh3+7d9g9Gh417s6OyrrJB67atXvrbfg\n0EPhtNOyAQHLl8PNN7Ni/36dHZmVwy67ZPf3feMbcMUVMGAAnHoq3H57Z0dmnaCsSUzSKEkrJK2U\ndEErdS5N65dKGt5WW0l9JC2S9KSk2yXV5KybnuqvkHR8Tvlhkh5N6y7JKe8l6YZUfp+kA9r/U7Cy\n69ULrroKnnwSpk2Dfk5eXZ6UTR68aFH2737kkdkoRtvplO10oqRuwGVkD69sBu6XVJ/zdGYkjQYO\niohaSSOAOcDINtpOAxZFxEUpuU0DpkkaCowDhgIDgTsk1aYnO88BJkVEo6SFkkalJztPAl5K2x8H\nXAhU7ZMQGxoaqKur6+ww2rRdca5bB8uXM3PuFfyx18aSmuzIda6mh5sY/JHta9uRqiXOv77615Lr\nbtcgkHPOab3S009nN7WXcF9gl/4d6qLKeU3scGBVRDQBSLoeGAssz6kzBrgaICIWS6qR1A84sEjb\nMcCRqf3VQANZIhsLzIuI9UCTpFXACElPA70jojG1uQY4Cbg19fWtVH4TWeKsWh3xH7vUoejFRpiV\nFOfixXD99bBsWfZ6+WU4+GD677MbMf2YkmLdketc1ZIcqiXOv64rPYm1+yCQCy/Mnl02bBi8973Z\n3IzvfS8cfzz06bNF1WpJDtUSZ0coZxIbCDyT8/5ZYEQJdQYCA4q07RsRa9LyGqBvWh4A3Fegr/Vp\nuUVzKt9i+xGxQdKrkvpExMul7GBHeumll3jnnXeK1nn99dd5/vnnAdhnn33o1q1bu8dR6lD0+dPn\nb/Ftevc336LfK6+y29vraXxsJVfdczvvens9q2t6c/XLr3DIRw/Zon1t8xret/oFnutTw3PHjeCl\n3rsTEo0PNHJKO++TVa+Sjtp6wl6fP4Y3732Qkbu8xbuXNrLfute4+ZYbWL33XltUffi+h+nzwF08\nu/Ip+r3/IN7q0Z03e/Tg7R7deX6vPXinwO/UttwSsC33I/pWg9KUM4lFifVKeQKhCvUXESGp1O1U\ntZumnsvhv7t7q/InBvbj1sM+BMBD9z3Ei+tfpHev3hz41NMcveSBreqvGNiPWz76oS3Kli5Zyvi+\n7+YzSx6FyD5sRfaxLh/Un18ePmxT3ZZTdLX3reTYH9/BLu9sZJcNG7Of72zkiU8MYeHUz2z1bbr2\n3ic58q5HeHP3XXlq/et8pOdfeWvvXdnlQ/vy2m1/2ioxrmcwK9LyHukFHkVoW9qWo7afP/4g/f/1\neJ5O73cF8ls2rW3i0LfeovaFFzlorx70fPNtevx1PT3ffJtrZ53GKwP23qrfU06/Am78VXZzdrdu\nm18LFsABW15mX/vmWqYuf5g9X1xHtDx8VRCI+q+PYd1+e26q2/JF8Mw7fk/N629s0c/qV9Yx9dnH\nWbvH7ixdsnSLL4ET7ryXvfPqA8w9eiTsO7DrJcaIKMsLGAncmvN+OnBBXp0fAeNz3q8gO7JqtW2q\n0y8t9wdWpOVpwLScNreSHb31A5bnlJ8KzMmpMzItdwdeaGVfwi+//PLLr21/lSvHtLzKeSS2BKiV\nNBh4jmzQxal5deqBKcD1kkYCayNijaSXirStByaQDcKYACzIKb9O0sVkpwlrgcZ0tLYuDRxpBE4H\nLs3r6z7gC8CdhXYkIvy8ejOzClS2JJauMU0BbgO6AVdGxHJJZ6f1V0TEQkmj0yCM14Ezi7VNXc8C\n5kuaBDRBdokkIpZJmg8sAzYAk9PIRIDJwFxgN2BhGpkIcCXwM0krgZeo4pGJZmY7I23+O29mZlZd\ndooZOyQdLqlR0kOS7pf0sZx17XaDtKQJ6SbsJyWdkVN+oKTFqc31knoUifVcScslPSbpwkqNM9X/\nqqSNkvrklFVMnJK+mz7LpZJulrRXzrqKiXN7qYTJBNphG4Mk/Z+kx9P/ya+k8g6ZdKC1z7aVWLsp\n+x3/ZQXHWCPpxvT/cpmkERUa5/T0b/6opOtSvxUXJ0BZL7hVyovsXrK/ScsnAP+XlocCDwM9yAYq\nrWLz0WkjcHhaXgiMSsuTgcvT8jjg+rTcB/gDUJNefwD2SuvmA6ek5TnAOa3EeRSwCOiR3u9XiXGm\n9YPIBsb8EehTiXECxwG7pOVZwKxKjHM7/093S3EPTvvxMHBwGX53+gEfSct7AE8ABwMXAV9P5Rd0\n4GdbUyTW84Frgfr0vhJjvBr4UlruDuxVaXGmbT0F9ErvbyAbO1BRcW6Kt73/01fiC5jH5j8mpwI/\nT8tbjJgkjVYkG/WYO6JxPPCjnDojcv4TvpDT75ycNj9K7QS8wOY/pluMvMyLcz5wdIHyioozrf8F\nMIwtk1jFxZnT/nOV+u++nf+nP86WI3i3GJ1bxt+lBWQz6awgu2cTskTXMkq47J9tK3HtD9xB9kXw\nl6ms0mLcC3iqQHmlxdmH7MvK3qmPX5J9IayoOFteO8XpRLJf8P+W9Cfgu2QfOmQ3SOfeCJ17s3VJ\nN0gDr0rap0hffchGXW4s0Fe+WuCIdHjdIOmjlRinpLHAsxHxSN6qioozz5fIvglWepylam2igLJR\nNlp4OLCY4pMOlPuzLeT7wD8DuXOSVVqMBwIvSPqppAcl/VjS7pUWZ2STPfw38Cey0eFrI2JRpcXZ\noss8ikXSIrJvB/m+AXwF+EpE/I+kvwWuIvtmUW5RoOxastsHHs0r/wbZv8feETFS2XW7+cB7yxwj\nbHuc04Hjc8o66haEbYnzXyKi5drIN4C3I+K6cgeYFIqzGrexiaQ9yKZmOy8iXpM2/5NHdO6kA5I+\nCzwfEQ9JqitUp7NjTLoDhwJTIuJ+SbPJvmBvUglxSnofMJXs1OCrwC8kfTG3TiXE2aLLHIlFxHER\n8eECr3qyc7L/k6reSDavI2TfDAbldLM/WeZvTsv55S1t3gMgqTvZ9Y+XCvQ1KJW9DNRIavmsLwAa\nWonzWeDmtD/3Axsl7VtJcZKdKz8QWCrpj2mbD0jqW0lx5iSwicBo4LScvjojzv1TeXsptN1nW6m7\nQ5QNSLkJ+FlEtNyXuUbZPKdI6g8830pc7f3ZFtrHTwBj0v/HecDRkn5WYTGSyp9Nv9uQ/S06FFhd\nYXF+FPh9RLyUjpJuJjt9XWlxZoqda+wqL+BB4Mi0fAxwf1puuSDZk+wP8x/YfEFyMdmMH2LrC5It\nM36MZ8sLkk+RXYzcu2U5rZsPjIvN53hbG4hwNjAzLQ8B/lSJcebFnHtNrKLiBEYBjwP75pVXVJzb\n+X+6e4p7cNqPcg3sENmk2d/PK7+IzbPoTGPri/xl/WyLxHskm6+JVVyMwF3AkLQ8I8VYUXEChwCP\nkd1XK7LBKF+utDg3xdve/+kr8UX2zWJx+qDvBYbnrPsXstE0K0gjGFP5YcCjad2lOeW9yP44rSSb\n6WNwzrozU/lKYEJO+YFp+yvJRvr0aCXOHsDP0nYfAOoqMc68mJ8iJbFKizOtfxp4KL0ur8Q4d+D/\n9QlkF+BXAdPL9LvzKbLrTA/nfI6jyP7Y3AE8CdxOzh+ajvhsi8R7JJtHJ1ZcjGQJ4n5gKdkRzl4V\nGufXyb4APkqWxHpUYpwR4ZudzcysenWZa2JmZrbzcRIzM7Oq5SRmZmZVy0nMzMyqlpOYmZlVLScx\nMzOrWk5iZmUmqZ+yR7GskrRE0q8l1W5HP03KeexNkXrzlD1+Zmqap+/k7YvcrPJ1mbkTzSqRskkG\n/wf4aUSMT2XDyCZPXbmN3QVtzFOZpgX6aETUpvc/pYPnWTTrSD4SMyuvo8gmH/5/LQUR8UhE3CNp\nd0l3SHpA0iOSxgCk8l9Lejg9UPBvc/o7N6f++wts73ZgoLKHQ34qlSn1e0yaPf0RSVdK6inpY5Ju\nSuvHSnpDUndJu0r6Q1k+EbN25CRmVl4fIptCrJA3gc9FxGHA0WSPv4BsWqfmiPhIZJMu35bT5oVU\nfw7wtQJ9ngj8ISKGR8Q9qSwk7Qr8lOy5esPIzsL8I9m8oh9J9T5NNkXQ4WTz3d23zXtr1sGcxMzK\nq9ipvF2A/5K0lOyJ3gMkvRt4BDhO0ixJn4qIdTltbk4/HySb/DdfodONAt4P/DEiVqWyq4EjIuId\n4A+SPgB8DLgYOIJszsS7S9lBs87kJGZWXo+TTYJayGnAvsChETGc7NEWu0bESrKHTz4K/Kekf81p\n81b6+Q7bdk07P5nmJru7yB5Xsx64k+yIzEnMqoKTmFkZRcRvgF6S/qGlTNKwdL1qT7KHOb4j6Sjg\ngLS+P/BmRFwLfI8soe1QGGSz3Q9ODzwEOB1oSMt3kz0E8fcR8SKwD9njQh7fwe2alZ1HJ5qV3+eA\n2ZIuILsO9keypHEt8EtJjwBLgOWp/oeB70raSHZ0dE6BPoPWT1VuVR4Rb0k6k+wpvd2BRrJnnJGW\n3012RAbZY0L65vdhVon8KBYzM6taPp1oZmZVy0nMzMyqlpOYmZlVLScxMzOrWk5iZmZWtZzEzMys\najmJmZlZ1XISMzOzqvX/ATsrvTRMUB4MAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa51511f350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_bins = 25\n",
    "n, bins, patches = plt.hist(f, num_bins, normed=1, facecolor='green', alpha=0.5)\n",
    "y = mlab.normpdf(bins, mu, sigma)\n",
    "plt.plot(bins, y, 'r--')\n",
    "plt.xlabel('Cash flow')\n",
    "plt.ylabel('Probability')\n",
    "plt.title(r'Histogram of cash flow')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Definition of the class policy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we construct the class policy defined by its particular cost structure summarized as: \n",
    "\n",
    "$v$: the holding cost per money unit of a positive cash balance at the end of the day; \n",
    "\n",
    "$u$: the shortage cost per money unit of a negative cash balance at the end oh the day, \n",
    "\n",
    "$\\gamma_{0}^{+}$: the fixed cost of transfer into account; \n",
    "\n",
    "$\\gamma_{0}^{-}$: the fixed cost of transfer from account; \n",
    "\n",
    "$\\gamma_{1}^{+}$: the variable cost of transfer into account; \n",
    "\n",
    "$\\gamma_{1}^{-}$: the variable cost of transfer from account. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A policy in the Miller-Orr cash management model is characterized by a set of bounds or control limits $\\{h,z,l\\}$ that determines the transfer to be done. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class policy(object):\n",
    "    \n",
    "    def __init__(self,gzeropos,gzeroneg,gonepos,goneneg,v,u,h,z,l):\n",
    "        self.gzeropos = gzeropos\n",
    "        self.gzeroneg = gzeroneg\n",
    "        self.gonepos = gonepos\n",
    "        self.goneneg = goneneg\n",
    "        self.v = v\n",
    "        self.u = u\n",
    "        self.h = h\n",
    "        self.z = z\n",
    "        self.l = l\n",
    "        self.daily_cost = []\n",
    "        self.cost_list =[]\n",
    "        \n",
    "    def transfer(self,ob):\n",
    "        \"\"\"Determines transfer x from parameters, opening balance (ob) and high, z, low\"\"\"\n",
    "        if ob > self.h:\n",
    "            x = self.z-ob\n",
    "        elif ob < self.l:\n",
    "            x = self.z-ob\n",
    "        else:\n",
    "            x = 0\n",
    "        return x\n",
    "    \n",
    "    def trans_cost(self,x):\n",
    "        \"\"\"Computes the cost of transfer x\"\"\"\n",
    "        cost = 0\n",
    "        if x<0:\n",
    "            cost = self.gzeroneg-self.goneneg*x\n",
    "        elif x>0:\n",
    "            cost = self.gzeropos+self.gonepos*x\n",
    "        return cost\n",
    "    \n",
    "    def holding_cost(self,final):\n",
    "        \"\"\"Computes the holding cost\"\"\"\n",
    "        cost = 0\n",
    "        if final>=0:\n",
    "            cost = self.v*final\n",
    "        else: \n",
    "            cost = -self.u*final\n",
    "        return cost\n",
    "                \n",
    "    def cost_calc(self,cf,ob):\n",
    "        \"\"\"Computes a vector of daily cost of cash flows and an opening balance\"\"\"\n",
    "        inibal = ob\n",
    "        bal = ob\n",
    "        if len(self.daily_cost)>0:\n",
    "            del self.daily_cost[:]\n",
    "        for element in cf:\n",
    "            trans = self.transfer(inibal)\n",
    "            bal = inibal+trans+element\n",
    "            self.daily_cost.append(self.trans_cost(trans)+self.holding_cost(bal))\n",
    "            inibal = bal\n",
    "        return self.daily_cost\n",
    "    \n",
    "    def cost_std(self,daily_cost):\n",
    "        \"\"\"Computes std dev of daily cost\"\"\"\n",
    "        #Check that before using this function daily_cost\n",
    "        return np.std(np.array(daily_cost,dtype=float))\n",
    "    \n",
    "    def cost_var(self,daily_cost):\n",
    "        \"\"\"Computes variance of daily cost\"\"\"\n",
    "        #Check that before using this function daily_cost\n",
    "        return st.pvariance(np.array(daily_cost,dtype=float))\n",
    "    \n",
    "    def up_dev(self,daily_cost,target):\n",
    "        \"\"\"Computes upside std dev above a given target of daily cost\"\"\"    \n",
    "        dc =np.array(daily_cost)\n",
    "        subset = dc[dc[:] > target]\n",
    "        length = len(daily_cost)\n",
    "        result = np.sqrt(np.sum((subset-target)**2/length))\n",
    "        return result\n",
    "    \n",
    "    def space(self,vh,vz,vl,cf,method):\n",
    "        \"\"\"Returns a list of daily_cost vectors from the combination of feasible h,z,l\"\"\"\n",
    "        method_name = str(method)\n",
    "        space={}\n",
    "        n = 0\n",
    "        for i in vh:\n",
    "            for j in vz:\n",
    "                for k in vl:\n",
    "                    if (i>=j and j>=k):\n",
    "                        n = n + 1\n",
    "                        self.h = i\n",
    "                        self.z = j\n",
    "                        self.l = k\n",
    "                        space[n]=[i,j,k,np.mean(self.cost_calc(cf,self.z)),self.risk_method(method_name,self.daily_cost)]\n",
    "        return space\n",
    "    \n",
    "    def space_up_dev(self,vh,vz,vl,cf,target):\n",
    "        \"\"\"Returns a list of daily_cost vectors from the combination of feasible h,z,l\"\"\"\n",
    "        space={}\n",
    "        n = 0\n",
    "        for i in vh:\n",
    "            for j in vz:\n",
    "                for k in vl:\n",
    "                    if (i>=j and j>=k):\n",
    "                        n = n + 1\n",
    "                        self.h = i\n",
    "                        self.z = j\n",
    "                        self.l = k\n",
    "                        space[n]=[i,j,k,np.mean(self.cost_calc(cf,self.z)),self.up_dev(self.daily_cost,target)]\n",
    "        return space\n",
    "    \n",
    "    def risk_method(self, argument, daily_cost):\n",
    "        \"\"\"Dispatch method for selcting cost_std or cost_var\"\"\"\n",
    "        method_name = str(argument)\n",
    "        # Get the method from 'self'. Default to a lambda.\n",
    "        method = getattr(self, method_name, lambda: \"nothing\")\n",
    "        # Call the method as we return it\n",
    "        return method(daily_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can initialize an object pol of the class policy with a particular cost structure and $h = 100000$, $z = 50000$ and $l = 25000$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pol = policy(100,5,0.05,0.001,0.00027,0.3,100000,50000,25000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then if the initial balance is 125000 the transfer is the amount necessary to reach the targe level that is -75000 with cost 80"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-75000"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pol.transfer(125000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pol.trans_cost(-75000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average daily cost over the entire data set with initial balance the target z is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7480.6350171600006"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pol.cost_calc(f,pol.z)\n",
    "np.mean(pol.daily_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Deriving the efficient frontier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we initialize an instance of the class policy with arbitrary $h,z,l$ values, for example, optimal values from Miller-Orr model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(442537.58082518424, 319707.98703182873, 258293.19013515097)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gzeropos = 5\n",
    "v = 0.00027\n",
    "l_miller = np.std(f_train)*2\n",
    "z_miller = l_miller+((3*gzeropos*np.std(f_train)**2)/(4*v))**(0.333333)\n",
    "h_miller = 3*z_miller-2*l_miller\n",
    "h_miller,z_miller,l_miller"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pol_miller = policy(50,5,0.1,0.001,0.00027,0.3,h_miller,z_miller,l_miller)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2641.4121707881218, 6885.6872754047972)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(pol_miller.cost_calc(f_test,z_miller)),np.std(pol_miller.cost_calc(f_test,z_miller))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3338.7059763075904, 9500.3402519796746)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_miller =np.mean(pol_miller.cost_calc(f_train,z_miller))\n",
    "y_miller =np.std(pol_miller.cost_calc(f_train,z_miller))\n",
    "x_miller,y_miller"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define two functions to transform a dict in an array and to derive the efficient frontier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def transform(d):\n",
    "    \"\"\"Transforms a dictionary in an array\"\"\"\n",
    "    S = np.zeros((len(d),5))\n",
    "    for i in range(len(d)):\n",
    "        S[i,:]=d[i+1][:]\n",
    "    return S\n",
    "\n",
    "def frontier(S,col):\n",
    "    \"\"\"Accepts a list of policies in increasing order of cost and returns the efficient frontier\"\"\"\n",
    "    P=[]\n",
    "    i = 0\n",
    "    P.append(S[i])\n",
    "    j = 1\n",
    "    while j<len(S):\n",
    "        if S[j][col]< S[i][col]:\n",
    "            i=j\n",
    "            P.append(S[j])\n",
    "        j = j+1 \n",
    "    return np.array(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define the limits to search for better policies and create the ranges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(310000, 830000, 190000, 710000, 130000, 650000)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma1 = np.std(f_train)\n",
    "sigma2 = 3*np.std(f_train)\n",
    "h_ini = int(np.round(h_miller-sigma1,-4))\n",
    "h_end = int(np.round(h_miller+sigma2,-4))\n",
    "z_ini = int(np.round(z_miller-sigma1,-4))\n",
    "z_end = int(np.round(z_miller+sigma2,-4))\n",
    "l_ini = int(np.round(l_miller-sigma1,-4))\n",
    "l_end = int(np.round(l_miller+sigma2,-4))\n",
    "h_ini,h_end,z_ini,z_end,l_ini,l_end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "140608"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vh = range(h_ini,h_end,10000)\n",
    "vz = range(z_ini,z_end,10000)\n",
    "vl = range(l_ini,l_end,10000)\n",
    "len(vh)*len(vz)*len(vl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function space produce a complete search over the input space formed by vh, vz y vl. Here we use standar deviation as a measure of risk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pol = policy(50,5,0.05,0.001,0.00027,0.3,h_miller,z_miller,l_miller)\n",
    "sp = pol.space(vh,vz,vl,f_train,'cost_std')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A the efficient frontier is stored in P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spa = transform(sp)\n",
    "S = spa[np.argsort(spa[:,3])]\n",
    "P = frontier(S,col=4)\n",
    "len(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cost index $\\theta_1$ and risk index $\\theta_2$ are computed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(754.01176111249993, 1226.5445175, 3312.9392483749189, 4050.2542982082164)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaxCost = np.max(P[:,3])\n",
    "MinCost = np.min(P[:,3])\n",
    "MaxRisk = np.max(P[:,4])\n",
    "MinRisk = np.min(P[:,4])\n",
    "MinCost,MaxCost,MinRisk,MaxRisk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot the efficient frontier in the space C-R and the values obtained from Miller-Orr optimal values for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEZCAYAAACuIuMVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHppJREFUeJzt3X+cV1W97/HXGwjFxLykIaAm2RDyKI8/KuqW+fWeBDID\nrZtgRWORXS/3qKeoI+o5wa3s0L1ZqY+rj2s/ZMzkXE4lB48cYDKnsptNPzApIsEbJihjWYoefwTy\nuX/sNbodhxkGvuv7ne933s/H4/uY9V177b3XcuN8Zq299tqKCMzMzHIZVu8KmJlZc3OgMTOzrBxo\nzMwsKwcaMzPLyoHGzMyycqAxM7OsHGisqUn6rKQ/SHowfT9b0gOSdkg6QdKvJL1tL47zuKRjctc3\nnestkjalc86swfneL2lN7vPY0CU/R2ONTNIW4BXAs6XsGyLiIklHAxuBoyLikVT+PuBvI+LWmle2\nOP9S4IGI+Ic+ytwOrIiIazKc/xjg/wEjImJ3tY9v1psR9a6A2X4K4MyI+F4v244GHikFGaW8DTWs\n377YYx1TG4j9/wtR+7l/cRBpmAOW9cdDZ9aUJL0dWAuMT0NQNwM7gOHALyVtSuW2SPrrlB4u6TJJ\nm9PQ2s8kTUjbdkt6VUofIOkLku6XtF3SdZIOTNsqkrZK+rikLkkPSjovbfso8D7g71Kd/qWXet8H\nvAq4NdVhpKSONAT4I+DfgYmS/qOkn0p6VFKnpDeXjtEh6dOS7kzHWCPp5WnzD9LPR9O2N0k6T9IP\nS/tPltQu6RFJGyW9t7RtaWrvKklPAJX9u1I2FDjQWDN40V/nEfFd4B3AgxExOiLeFxGj0+bjI6Kl\nu2j6AHwcmAO8IyIOAT4MPNXL+ZYArwb+Kv2cAHyqtH0scAgwHpgH/C9JL4uI64FvAp9PdZrVS72P\nBX5P0Us7JCL+kjZ9APgIcDBFsLkN+DIwBvgicJuk/1A61LnAeRTDiiOBT6T8U9LPl6Xj31U+v6SX\nAu3ATcDh6b/HtZKO63Hsz0TEwcCPevnvY/YCDjTW6ASskPTn0mdeadtAfAS4PCI2AUTEPRHxpxec\nrBi6Oh/4eEQ8GhFPAP9I8Qu5207g0xHxbET8G/AE8JoedR6IAJZGxG/SMNU04LcR8c2I2B0R/0Rx\nL2pmqfwNEbE5Ip4GlgMn7OW5zwR+FxFt6dh3A98B3lsqsyIifgwQEc8MsC02BPkejTW6AGbt4R7N\nQB0J3NdPmcOBg4Cfp9slUPzyLv/R9kiP+xZPUvRE9scDpfR4il5P2f0pv9v2UvqpAZz/lcBUSX8u\n5Y0AbkzpALbu5bHMAAcas7IHKIbC+pos8EeKX9xTIuKhfTjHvt7EL++3DXh3j+2vBP6tCuf/PfD9\niJg2gLqZ9clDZ9YMqjKDCvgq8BlJr1bheEljygVST+UrwJclHQ4gaYKkvf3F3EVxs3+gym1cBUyS\ndK6kEZJmA5OBf91D+bI/ALuBY/ew/bZ07A9Iekn6vEHS5H6Oa7ZHDjTWDG5Ns7i6P98ubev5F3xf\nf9F/keJ+xlrgMYqAcmAv+10CbAbukvQYxc3zSXt5jq8BU9K9pO/0Ua6n546Z7hudCSyg6GF9gmLy\nwJ96K09pwkNEPAlcAfxI0p8kTe2x/XGKe0BzKHpOD1HcgxrZ81hmeyvbA5uSvg68E3g4Il6X8sYA\n/4eim78FOCciHk3bLqWY5fMscFFErE35JwNLKf6HXxURF6f8AyjGjU8CHgFmR8T9WRpjZmb7LGeP\n5gZgRo+8hUB7REwCbk/fkTQFmA1MSftc2/1gGnAdMC9NR22R1H3MeRQ3XVuALwGfz9gWMzPbR9kC\nTUT8EPhzj+yZQFtKtwFnpfQsYFlE7IyILRTDElMljQNGR0RnKndjaZ/ysb4N/HXVG2FmZvut1vdo\nxkZEV0p3UTzYBsW0zPKUya0UD8H1zN+W8kk/HwCIiF3AYz1v3JqZWf3VbTJAWqvJNxXNzJpcrZ+j\n6ZJ0RERsT8NiD6f8bcBRpXJHUvRktqV0z/zufY4GHpQ0gmJJjRc8xQ0gycHMzGwfRERVprPXukez\nEmhN6VZgRSl/TlpAcCLQAnRGxHZgh6SpaXLAXOBfejnWf6aYXNCriGjaz6JFi+peB7fP7RtqbRsK\n7aumbD0aScuAU4HDJD1AsejgEmB5WotqC3AOQERskLSc4onsXcD8eL6l8ymmN4+imN68OuV/DfhG\nWoX3EV641pSZmQ0S2QJNRJy7h01v30P5zwGf6yX/58Dresl/hhSozMxs8PLKAA2uUqnUuwpZuX2N\nq5nbBs3fvmpq+lc5S4pmb6OZWbVJIhp0MoCZmQ0xDjRmZpaVA42ZmWXlQGNmZlk50JiZWVYONGZm\nlpUDjZmZZVXrRTXNzCyDH9x2G2uvvpoRzzzDrgMOYNpFF/G2d76z3tUCHGjMzBreD267jTUXX8wV\n9933XN7lKT0Ygo2HzszMGtzaq69+QZABuOK++2i/5po61eiFHGjMzBrciGee6TV/+NNP17gmvXOg\nMTNrcLsOOKDX/GcPPLDGNemdA42ZWYObdtFFXH7ssS/Iu+zYYzn9wgvrVKMX8urNZmZN4Ae33Ub7\nNdcw/OmnefbAAzn9wgv3ayJANVdvdqAxM7MX8WsCzMysYTjQmJlZVg40ZmaWlQONmZll5UBjZmZZ\nOdCYmVlWdQk0ki6WtF7SryRdnPLGSGqXdK+ktZIOLZW/VNImSRslTSvln5yOs0nSVfVoi5mZ9a3m\ngUbSa4GPAG8A/go4U9KxwEKgPSImAben70iaAswGpgAzgGsldc/tvg6YFxEtQIukGTVtjJmZ9ase\nPZrJwE8i4umIeBb4PvAeYCbQlsq0AWel9CxgWUTsjIgtwGZgqqRxwOiI6EzlbiztY2Zmg0Q9As2v\ngFPSUNlBwBnAkcDYiOhKZbqAsSk9Htha2n8rMKGX/G0p38zMBpGav/gsIjZK+jywFvh34G7g2R5l\nQlLV1o1ZvHjxc+lKpUKlUqnWoc3MmkJHRwcdHR1Zjl33tc4kXUHRM7kYqETE9jQsdkdETJa0ECAi\nlqTyq4FFwP2pzHEp/1zg1Ii4oMfxvdaZmdkANfxaZ5JekX4eDbwbuBlYCbSmIq3AipReCcyRNFLS\nRKAF6IyI7cAOSVPT5IC5pX3MzGyQqPnQWfItSS8HdgLzI+IxSUuA5ZLmAVuAcwAiYoOk5cAGYFcq\n391FmQ8sBUYBqyJidW2bYWZm/an70FluHjozMxu4hh86MzOzocOBxszMsnKgMTOzrBxozMwsKwca\nMzPLyoHGzMyycqAxM7OsHGjMzCwrBxozM8vKgcbMzLJyoDEzs6wcaMzMLCsHGjMzy8qBxszMsnKg\nMTOzrBxozMwsKwcaMzPLyoHGzMyycqAxM7OsHGjMzCwrBxozM8vKgcbMzLJyoDEzs6zqEmgkXSrp\n15LWS7pZ0gGSxkhql3SvpLWSDu1RfpOkjZKmlfJPTsfYJOmqerTFzMz6VvNAI+kY4HzgpIh4HTAc\nmAMsBNojYhJwe/qOpCnAbGAKMAO4VpLS4a4D5kVEC9AiaUYNm2JmZnuhHj2aHcBO4CBJI4CDgAeB\nmUBbKtMGnJXSs4BlEbEzIrYAm4GpksYBoyOiM5W7sbSPmZkNEjUPNBHxJ+BK4PcUAebRiGgHxkZE\nVyrWBYxN6fHA1tIhtgITesnflvLNzGwQGVHrE0o6Fvhb4BjgMeCfJX2gXCYiQlJU65yLFy9+Ll2p\nVKhUKtU6tJlZU+jo6KCjoyPLsRVRtd/ne3dCaTZwekR8JH2fC7wJ+E/AaRGxPQ2L3RERkyUtBIiI\nJan8amARcH8qc1zKPxc4NSIu6HG+qHUbzcwanSQiQv2X7F897tFsBN4kaVS6qf92YANwK9CayrQC\nK1J6JTBH0khJE4EWoDMitgM7JE1Nx5lb2sfMzAaJmg+dRcQvJd0I/AzYDfwCuB4YDSyXNA/YApyT\nym+QtJwiGO0C5pe6KPOBpcAoYFVErK5hU8zMbC/UfOis1jx0ZmY2cI0+dGZmZkOIA42ZmWXlQGNm\nZlk50JiZWVYONGZmlpUDjZmZZeVAY2ZmWTnQmJlZVg40ZmaWlQONmZll5UBjZmZZOdCYmVlWDjRm\nZpaVA42ZmWXlQGNmZlk50JiZWVYONGZmlpUDjZmZZeVAY2ZmWTnQmJlZVg40ZmaWlQONmZll5UBj\nZmZZ1TzQSHqNpHWlz2OSLpI0RlK7pHslrZV0aGmfSyVtkrRR0rRS/smS1qdtV9W6LWZm1r+aB5qI\n+G1EnBgRJwInA08CtwALgfaImATcnr4jaQowG5gCzACulaR0uOuAeRHRArRImlHb1piZWX/qPXT2\ndmBzRDwAzATaUn4bcFZKzwKWRcTOiNgCbAamShoHjI6IzlTuxtI+ZmY2SNQ70MwBlqX02IjoSuku\nYGxKjwe2lvbZCkzoJX9byjczs0FkRL1OLGkk8C7gkp7bIiIkRbXOtXjx4ufSlUqFSqVSrUObmTWF\njo4OOjo6shxbEVX7fT6wE0uzgP8aETPS941AJSK2p2GxOyJisqSFABGxJJVbDSwC7k9ljkv55wKn\nRsQFPc4T9WqjmVmjkkREqP+S/avn0Nm5PD9sBrASaE3pVmBFKX+OpJGSJgItQGdEbAd2SJqaJgfM\nLe1jZmaDRF16NJJeStEjmRgRj6e8McBy4GhgC3BORDyatl0GfBjYBVwcEWtS/snAUmAUsCoiLurl\nXO7RmJkNUDV7NHUbOqsVBxozs4FrlqEzMzMbAhxozMwsKwcaMzPLqt9Ak26498w7M091zMys2exN\nj+Yrkl7X/SU9r/KpfFUyM7Nm0u+sM0mvAr4FvA84BfggcGZEPJa/evvPs87MzAau5tObJb2G4mHI\n+4F3R8ST1Th5LTjQmJkNXE0CjaT1PbJeATwK/IViObLjq1GB3BxozMwGrpqBpq9FNd/VS173b+yq\nnNzMzJrfHicDRMSW9P6X4cD2lH4VxfthHq1J7czMrOHtzayz7wC7JL0a+N/AUcDNWWtlZmZNY28C\nze6I2AW8G7gmIj4JjMtbLTMzaxZ7E2j+Iul9FNOa/zXlvSRflczMrJnsTaD5MPBm4IqI+F16ruam\nvNUyM7Nm4dcEmJnZi9RkerOkf46I9/byPA000HM0ZmZWX309sDk+Ih6UdEwvmyMi7s9ZsWpxj8bM\nbOBq0qOJiAfTzy09Ti7gHIrlaMzMzPq0x8kAkg6WtEDStZLmSxom6Wzg18D7a1dFMzNrZH0NnX0H\n2AH8GJhG8aDm08BFEXF3zWq4nzx0ZmY2cLVaVPOe7hv+koYDDwGvjIinqnHiWnGgMTMbuGoGmr6e\no3m2OxERzwLbGi3ImJlZ/fXVo3kWKL93ZhTQHWgiIg7JXLeqcI/GzGzgatKjiYjhETG69BlRSu9X\nkJF0qKRvSfqNpA2SpkoaI6ld0r2S1ko6tFT+UkmbJG2UNK2Uf7Kk9WnbVftTJzMzy2NvlqDJ4Spg\nVUQcBxwPbAQWAu0RMQm4PX1H0hRgNjAFmAFcm6ZYA1wHzIuIFqBF0ozaNsPMzPpT80Aj6WXAKRHx\ndYCI2BURjwEzgbZUrA04K6VnAcsiYmd6pmczMFXSOGB0RHSmcjeW9jEzs0GiHj2aicAfJN0g6ReS\nviLppcDYiOhKZbqAsSk9Htha2n8rMKGX/G0p38zMBpG+XuWc85wnAX8TET+V9GXSMFm3iAhJVbuD\nv3jx4ufSlUqFSqVSrUObmTWFjo4OOjo6shy75qs3SzoC+HFETEzf3wpcSvGa6NMiYnsaFrsjIiZL\nWggQEUtS+dXAIoolcO5I93mQdC5wakRc0ON8nnVmZjZAtXqOJouI2A48IGlSyno7xbI2twKtKa8V\nWJHSK4E5kkZKmgi0AJ3pODvSjDUBc0v7mJnZIFGPoTOAC4FvShoJ3Ad8CBgOLJc0D9hCsXAnEbFB\n0nJgA7ALmF/qoswHllI847MqIlbXshFmZtY/v/jMzMxepKGHzszMbGhxoDEzs6wcaMzMLCsHGjMz\ny8qBxszMsnKgMTOzrBxozMwsKwcaMzPLyoHGzMyycqAxM7OsHGjMzCwrBxozM8vKgcbMzLJyoDEz\ns6wcaMzMLCsHGjMzy8qBxszMsnKgMTOzrBxozMwsKwcaMzPLyoHGzMyycqAxM7OsHGjMzCyrugQa\nSVsk3SNpnaTOlDdGUrukeyWtlXRoqfylkjZJ2ihpWin/ZEnr07ar6tEWMzPrW716NAFUIuLEiHhj\nylsItEfEJOD29B1JU4DZwBRgBnCtJKV9rgPmRUQL0CJpRi0bYWZm/avn0Jl6fJ8JtKV0G3BWSs8C\nlkXEzojYAmwGpkoaB4yOiM5U7sbSPmZmNkjUs0fzXUk/k3R+yhsbEV0p3QWMTenxwNbSvluBCb3k\nb0v5ZmY2iIyo03nfEhEPSTocaJe0sbwxIkJSVOtkixcvfi5dqVSoVCrVOrSZWVPo6Oigo6Mjy7EV\nUbXf5/tWAWkR8ARwPsV9m+1pWOyOiJgsaSFARCxJ5VcDi4D7U5njUv65wKkRcUGP40e922hm1mgk\nERE9b3Hsk5oPnUk6SNLolH4pMA1YD6wEWlOxVmBFSq8E5kgaKWki0AJ0RsR2YIekqWlywNzSPmZm\nNkjUY+hsLHBLmjg2AvhmRKyV9DNguaR5wBbgHICI2CBpObAB2AXML3VR5gNLgVHAqohYXcuGmJlZ\n/+o+dJabh87MzAauoYfOzMxsaHGgMTOzrBxozMwsKwcaMzPLyoHGzMyycqAxM7OsHGjMzCwrBxoz\nM8vKgcbMzLJyoDEzs6wcaMzMLCsHGjMzy8qBxszMsnKgMTOzrBxozMwsKwcaMzPLyoHGzMyycqAx\nM7OsHGjMzCwrBxozM8vKgcbMzLJyoDEzs6wcaKpozZo1TJv2HqZNew9r1qypd3XMzAaFugUaScMl\nrZN0a/o+RlK7pHslrZV0aKnspZI2SdooaVop/2RJ69O2q+rRjm5r1qzh7LNbaW+fSXv7TM4+u9XB\nxsyM+vZoLgY2AJG+LwTaI2IScHv6jqQpwGxgCjADuFaS0j7XAfMiogVokTSjhvV/gSuvvJ6nnvo8\n0Aq08tRTn+fKK6+vV3XMzAaNugQaSUcCZwBfBbqDxkygLaXbgLNSehawLCJ2RsQWYDMwVdI4YHRE\ndKZyN5b2MTOzQaJePZovAZ8EdpfyxkZEV0p3AWNTejywtVRuKzChl/xtKb8uFiz4KKNGXQJ8Angz\nw4Yt4NRTT6pXdczMBo0RtT6hpDOBhyNinaRKb2UiIiRFb9v2xeLFi59LVyoVKpVeT7tfpk+fzuWX\nX8inPnUlu3d/id274YorLuH1r38906dPr/r5zMyqqaOjg46OjizHVkTVfp/v3QmlzwFzgV3AgcAh\nwHeANwCViNiehsXuiIjJkhYCRMSStP9qYBFwfypzXMo/Fzg1Ii7ocb6oVRunTXsP7e0zKe7TALRx\n+ukrWbv22zU5v5lZtUgiItR/yf7VfOgsIi6LiKMiYiIwB/heRMwFVvL8b+hWYEVKrwTmSBopaSLQ\nAnRGxHZgh6SpaXLA3NI+ZmY2SNR86KwX3d2NJcBySfOALcA5ABGxQdJyihlqu4D5pS7KfGApMApY\nFRGra1jvF1mw4KPceWcrTz1VfB816hIWLGjreyczsyZX86GzWqvl0BkUz9N0T2tesOCjvj9jZg2p\nmkNnDjRmZvYiDX2PxszMhhYHmjrwmmhmNpR46KzGutdEK5arKSYM3HJLm+/lmNmg4qGzBlavNdHc\nizKzenGgGQKeX1l6Iu3tD3LGGe/niiuuqHe1zGyI8NBZjdVj6KxYsWAicBNQnHfYsI+xatUyD9mZ\nWa+qOXQ2GB7YHFKmT5/OLbe0lZ61qdX9mR9RBJli8YXdu4thPAcaM8vNPZohYM2aNZxxxvvZvftK\n4AjgemADo0f/O2960xv8YKmZvYgnA9iATJ8+nU9/+mNIfwN8AJgI/JHHH/+M3wZqZtm5RzOEnHRS\nhXXrPkSxTulMnu/dPMiJJw7nF7+4s671M7PBwz0a2yeHHfby0rf1FPdrJgKwbt0Gz0QzsywcaIaQ\n598COhG4gWIY7SbgLcBr+Pu/X8KrX32in7Uxs6ry0NkQ07269F13/ZTHH59AEWRuogg6bcCpwE+A\nJ2ltfRdLly6tX2XNrG68evMAOND07vmZaC3ABRT3bQS0A+dTTIdez7hxE3jta1/rmWlmQ4zv0dh+\n656JBr8t5f6EIsh0D6eN4qGHLvPMNDPbLw40Q9jll1/OZz+7gGHDPkZx3+ZJnn+w83fAF6j1mmxm\n1nwcaIa4yy+/nFWrlnH66b9j3LiDgY31rpKZNRnfo7EXOO+882hrWwF8hGJywBcAv87AbKjxZIAB\ncKAZuO6ZaX/8YxcwgsMOe7knA5gNMQ40A+BAY2Y2cJ51ZmZmDaPmgUbSgZJ+IuluSRsk/WPKHyOp\nXdK9ktZKOrS0z6WSNknaKGlaKf9kSevTtqtq3RYzM+tfzQNNRDwNnBYRJwDHA6dJeiuwEGiPiEnA\n7ek7kqYAs4EpwAzgWknd3bnrgHkR0QK0SJpR29bUX0dHR72rkJXb17iauW3Q/O2rproMnUXEkyk5\nEhgO/JliOeG2lN8GnJXSs4BlEbEzIrYAm4GpksYBoyOiM5W7sbTPkNHs/9jdvsbVzG2D5m9fNdUl\n0EgaJuluoAu4IyJ+DYyNiK5UpAsYm9Ljga2l3bcCE3rJ35byzcxsEKnLq5wjYjdwgqSXAWskndZj\ne0jyVDEzsyZQ9+nNkv4BeIriCcFKRGxPw2J3RMRkSQsBImJJKr8aWATcn8ocl/LPBU6NiAt6HN8B\ny8xsH1RrenPNezSSDgN2RcSjkkYBpwP/nWL54FaKhbZagRVpl5XAzZK+SDE01gJ0pl7PDklTgU5g\nLnB1z/NV6z+UmZntm3oMnY0D2iQNo7hH9I2IuF3SOmC5pHnAFuAcgIjYIGk5sAHYBcwvPYE5H1gK\njAJWRcTqmrbEzMz6VfehMzMza24NvzKApC2S7pG0TlJnyhvww5+DhaSvS+qStL6U1zQPs+6hfYsl\nbU3XcJ2kd5S2NVr7jpJ0h6RfS/qVpItSfsNfwz7a1hTXr9kfJu+jffmvX0Q09IfixSljeuT9D+Dv\nUvoSYElKTwHuBl4CHEPxTM6werehR91PAU4E1u9je7p7qZ3AG1N6FTCj3m3ro32LgI/3UrYR23cE\ncEJKH0zxZrnjmuEa9tG2Zrp+B6WfI4C7gLc2w7Xrp33Zr1/D92iSnjf8B/Lw5xtrUsO9FBE/pHiA\ntaxpHmbdQ/vgxdcQGrN92yPi7pR+AvgNxSSWhr+GfbQNmuf6NfXD5HtoH2S+fs0QaAL4rqSfSTo/\n5Q304c/Bbig8zHqhpF9K+lppaKKh2yfpGIre209osmtYattdKasprp+a/GHyPbQPMl+/Zgg0b4mI\nE4F3AP9N0inljVH07fqa8dBQsyH2oj2N6DqKd0mfADwEXFnf6uw/SQcD3wYujojHy9sa/Rqmtn2L\nom1P0ETXLyJ2R7EO45HA29TLw+Q08LXrpX0VanD9Gj7QRMRD6ecfgFsohsK6JB0BkLp5D6fi24Cj\nSrsfmfIGu4G0Z2vKP7JH/qBtZ0Q8HAnwVZ4fzmzI9kl6CUWQ+UZEdD8P1hTXsNS2m7rb1mzXDyAi\nHgNuA06mSa5dWal9r6/F9WvoQCPpIEmjU/qlwDRgPc8//AkvfvhzjqSRkiaSHv6sba33yYDaExHb\ngR2SpkoSxcOsK3oedLBI//N2O5viGkIDti/V52vAhoj4cmlTw1/DPbWtWa6fpMO6h430/MPk62iC\nawd7bl93EE3yXL9azXbI8aHo7t2dPr8CLk35Y4DvAvcCa4FDS/tcRnFTayMwvd5t6KVNy4AHgb8A\nDwAf2pf2UPwltj5tu7re7eqjfR+muJl4D/DL9A92bAO3763A7vRvcl36zGiGa7iHtr2jWa4f8Drg\nF6l99wCfTPkNf+36aV/26+cHNs3MLKuGHjozM7PBz4HGzMyycqAxM7OsHGjMzCwrBxozM8vKgcbM\nzLJyoDHLQNIRkv5J0ua0Dt9tkloGeIzLctXPrJb8HI1ZlaWnpf8vcENEXJ/yjgcOiYg7B3CcxyNi\ndKZqmtWMezRm1Xca8JfuIAMQEfdExJ2S/md6YdQ9ks6BYgkXST9IL51aL+mtkpYAo1LeN+rVELNq\ncI/GrMpUvHnymIj4eI/89wD/BZgOHA78FJgKvB84ICI+J2kYxcupnnCPxprFiHpXwKwJ7emvt7cA\nN0fx193Dkr4PvIFiYdevp5WRV0TEL2tUT7Oa8NCZWfX9mmLRwd70fJNhRPHW0VMollpfKmluzsqZ\n1ZoDjVmVRcT3gANKb3ztngzwKDA7veXwcOBtQKeko4E/RMRXKZbhPzHttlOSRx2s4fkfsVkeZwNf\nlnQJ8DTwO+BjwMEUy7EHxTLtD0v6IPBJSTuBx4EPpmNcD9wj6ecR4V6ONSxPBjAzs6w8dGZmZlk5\n0JiZWVYONGZmlpUDjZmZZeVAY2ZmWTnQmJlZVg40ZmaWlQONmZll9f8BGwDkxQNxXrcAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa510346290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(P[:,3],P[:,4])\n",
    "plt.plot(x_miller,y_miller,'ro')\n",
    "plt.xlabel('Cost')\n",
    "plt.ylabel('Risk')\n",
    "plt.title(r'Efficient frontier')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now in the normalized space $\\theta_1$-$\\theta_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAEeCAYAAABVItqsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHKVJREFUeJzt3Xu4XHV97/H3hyTARm6m0RwNAY6GChyLUsrVetwKJJGm\nUIhCA2poPULbhx6tu4dUWyXloDRIODxVixAQYp9qTlsIDUj3kAPuQjjcJQlo4BAumgSk3C9lA4F8\nzx9r7TCZ7Mus2bNm1pr5vJ5nnqzbrPXdG+azf+u31vqNIgIzszzt0O4CzKzzOWjMLHcOGjPLnYPG\nzHLnoDGz3DlozCx3DpoOJOk8SU9LeiKdP1HSBkkvSfqwpAck/dc69vOypH3zrjc91kckPZwe8/gW\nHO80SZW8j2MJ+T6a8pH0OPBu4K2qxVdGxH+XtDfwIDA9Ip5Nt38E+FJEXNfyYpPjXwVsiIivjbLN\nTcC1EfHtHI6/L/AoMDEitjR7/za2ie0uwBoSwJyIuHmYdXsDz1aFjNJlP29hfY0Yscb0ZyDG/1dR\n43x/shNpBwdWNj516iCSjgFuBN6bnoL8EHgJmACskfRwut3jko5OpydI+qqk9emp1T2SpqXrtkh6\nXzq9k6QLJf1C0q8kXSJp53Rdr6SNkr4s6SlJT0g6PV13BnAqcHZa078MU/cjwPuA69IadpQ0kJ4C\n3gb8B/CfJR0l6W5JL0i6S9KRVfsYkHSupFXpPiqSfi1dfUv67wvpuiMknS7p1qr37y9ppaRnJT0o\n6dNV665Kf94bJL0C9I7vv1QXigi/SvYCHgOOHmHdx0hOU6qXbQHeV/P+T6TT/wNYC+yXzh8ETK59\nH/C/gGuBPYFdgRXAN9N1vcBmYCFJqH2SJBz2SNdfCZxbx8/0iar5AeBx4ACSP4hTgeeB09L53wee\nA95Ztf3DwAxgZ+AnwPnpun3Sn2WHqv2fDtyaTr8D2ADMT/f9YeBp4IB0/VXAC8CR6fxO7f5/oGwv\nt2jKScC1kp6ven2+al0W/w34y4h4GCAi1kbEc9scLDl1+QLw5Yh4ISJeAc4n+bAP2UwSJm9FxL8C\nrwAfqKk5iwCuioh1kZymzAQeioh/iIgtEbGMpC/q+Krtr4yI9RHxGvCPJIFRz7HnAI9FxNJ036uB\na4BPV21zbUTcDhARr2f8Wbqe+2jKKYATYvg+mqz2Ah4ZY5t3AbsA96bdJZB8eKv/UD0b2/ZbvErS\n8hmPDVXT7wV+WbP+F+nyIb+qmh7McPx9gMMlPV+1bCLwg3Q6gI117suG4aCxDSSnG6N1Fj9D8sE9\nMCKebOAYjXbiVr9vE3BSzfp9gH9twvF/CfxbRMzMUJtl4FOn8mrKFRTgcuB/SpqhxEGSJldvkLZU\nlgAXS3oXgKRpkur9YD5F0tmbVfXPeAPw65LmSZoo6RRgf+D6Ebav9jRJH837R1j/43Tfn5E0KX0d\nKmn/MfZrdXLQlNd16VWcodfVVetq/4KP9hf9IpL+jBuBF0kCZedh3rcAWA/cIelFYCXw63Ue4wrg\nwLQv6ZpRtqu1dZ9pv9EcoI+khfXnJJf4nxtu+3Q60ve+CnwDuE3Sc5IOr1n/Mkkf0O+TtJyeJOmD\n2rF2X9aYlt6wJ+n7wO8A/x4RvzHCNn9LctXiVeD0iLivZQWaWS5a3aK5Epg90kpJxwEzImI/4Azg\nklYVZmb5aWnQRMStJPdCjOR4YGm67Z3AnpKmtqI2M8tP0fpoprHtJc2NJJdfzazEihY0sH0Pvzvh\nzEquaPfRbAKmV83vlS7bjiQHkFmbRESmS/5Fa9GsAD4HIOkI4IWIeGqkjXt6ptLf39/25zjqeZ1z\nzjltr6HTay5bvWWqefXqt6cb0dIWjaQfkTz0N0XSBuAcYBJARFwaETdIOk7SepKH8v5gtP0tX76U\nWbNm5V22WVe74AK44gpYswZ23nns7YfT0qCJiHl1bHNWvftzyJjl64ILYMkSGBhoPGSgeKdOHau3\nt7fdJWRWtprLVi8Uu+bqkJk2bXz7Ku1QnpKirLWbFd1oISOJKHlnsJm1WTNbMkMcNGa2VR4hAw4a\nM0vlFTLgoDEz8g0ZcNCYdb28QwYcNGZdrRUhAw4as67VqpABB41ZV2plyICDxqzrtDpkwEFj1lXa\nETLgoDHrGu0KGXDQmHWFdoYMOGjMOl67QwYcNGYdrQghAw4as45VlJABB41ZRypSyICDxqzjFC1k\nwEFj1lGKGDLgoDHrGEUNGXDQmHWEIocMOGjMSq/oIQMOGrNSK0PIgIPGrLTKEjLgoDErpTKFDDho\nzEqnbCEDHRA0lUqFmTPnMnPmXCqVSrvLMctVGUMGSv6VuP39/Zx44nwGBxcB0NOzgOXLlzJr1qw2\nV2fWfEUJmUa+ErfUQXPssSexcuXxwPx06VKOPXYFN954dTtLM2u6ooQM+Lu3zTpSkUKmURPbXcB4\n9PWdwapV8xkcTOZ7ehbQ17e0vUWZNVEnhAyU/NQpIqhUKixefBmQBI/7Z6xTFDVkuq6Ppqy1m42l\nqCED7qMx6whFDplGOWjMCqQTQwbaEDSSZkt6UNLDkhYMs36KpH5JqyU9IOn0Vtdo1g6LFnVmyECL\n+2gkTQAeAo4BNgF3A/MiYl3VNguBnSLiK5KmpNtPjYg3a/blPhrrGIsWweWXlyNkytBHcxiwPiIe\nj4jNwDLghJptngR2T6d3B56tDRmzTlKmkGlUq++jmQZsqJrfCBxes80S4GZJTwC7ASe3qDazluuG\nkIHWt2jqOdf5KrA6It4LfBj4rqTd8i3LrPUuuKA7QgZa36LZBEyvmp9O0qqpdhTwDYCIeETSY8AH\ngHtqd7Zw4cKt0729vfT29ja3WrOclOnq0sDAAAMDA+PaR6s7gyeSdO4eDTwB3MX2ncEXAS9GxF9L\nmgrcCxwUEc/V7MudwVZKZQqZ4TTSGdzSFk1EvCnpLKACTACuiIh1ks5M118KfBO4UtIaklO7s2tD\nxqysyh4yjfIjCGYt0ikhU4bL22ZdqVNCplEOGrOcdXvIgIPGLFcOmYSDxiwnDpm3OWjMcuCQ2ZaD\nxqzJHDLbc9CYNZFDZngOGrMmcciMzEFj1gQOmdE5aMzGySEzNgeN2Tg4ZOrjoDFrkEOmfg4aswY4\nZLJx0Jhl5JDJzkFjloFDpjEOGrM6OWQa56Axq4NDZnwcNGZjcMiMn4PGbBQOmeZw0JiNwCHTPA4a\ns2E4ZJrLQWNWwyHTfA4asyoOmXw4aMxSDpn8OGjMcMjkzUFjXc8hkz8HjXU1h0xrOGisazlkWsdB\nY13JIdNaDhrrOg6Z1nPQWFdxyLSHg8a6hkOmfRoOGkm7pv9OkjSheSWZNZ9Dpr0aChpJZwNfl3QR\nsAfwvaZWZdZEDpn2m9jg++5MX5uBU/ApmBWUQ6YYGg2IV4HTI+KtiPghcEu9b5Q0W9KDkh6WtGCE\nbXol3SfpAUkDDdZoXc4hUxyKiNE3kN4PfIrkFGkDcHtErG7oYElfzkPAMcAm4G5gXkSsq9pmT+A2\nYFZEbJQ0JSKeGWZfMVbt1r0cMvmRREQoy3vqadF8ArgauAk4GjhP0r2STmugxsOA9RHxeERsBpYB\nJ9RscypwdURsBBguZMxG45ApnnqCZgdg14i4CbguIuYARwFbJP1xxuNNI2kVDdmYLqu2HzBZ0k8k\n3SPpsxmPYV3MIVNM9XQGXwZ8UdKFwMuSngUeBe4CTsp4vHrOdSYBv0nSetoFuF3SHRHxcMZjWZdx\nyBTXmEGTdoRcLOm7JKdRRwK/CzwL/Cjj8TYB06vmp5O0aqptAJ6JiEFgUNItwIeA7YJm4cKFW6d7\ne3vp7e3NWI51CodMfgYGBhgYGBjXPsbsDG4mSRNJOoOPBp4gaRXVdgbvD3wHmAXsRHIZ/ZSI+HnN\nvtwZbIBDptUa6Qxu9D6ahkTEm5LOAirABOCKiFgn6cx0/aUR8aCkfmAtsAVYUhsyZkMcMuXQ0hZN\nM7lFYw6Z9sjr8vbQzg8cZllvloOZNYtDplyy3Bn8j5IWKLGLpG8Df5NXYWYjcciUT5agOZzkKtHt\nJJ24T5LcT2PWMg6ZcsoSNG8Cg0APsDPwaERsyaUqs2E4ZMorS9DcBbwG/BbwUeBUSf+US1VmNRwy\n5Vb3VSdJh0bE3TXLPhsRf59LZWPX46tOXcIhUyy5XnUC7pX0WUlfTw+2N/D/shzMLCuHTGfIEjR/\nR/L4wanp/CvAd5tekVnKIdM5stwZfHhEHCzpPoCIeE7SpJzqsi7nkOksWVo0b1QPQi7pXSSPCJg1\nlUOm82QJmm8Dy4F3S/omySh45+dSlXUth0xnyvSsk6QDSIaKEHBT9VPXrearTp3HIVMOjVx1qmfM\n4L6q2SAJmaFpIuKiLAdsFgdNZ3HIlEdew0TsRhIqHwAOBVaQhM0ckpv4zMbFIdP5stywdytwXES8\nnM7vBtwQER/Nsb7R6nGLpgM4ZMon7xv23k3yhXFDNqfLzBrikOkeWYLmB8BdkhZK+muSITaX5lNW\n+VUqFWbOnMvMmXOpVCrtLqdwHDLdJetVp0NIHqgM4JaIuC+vwuqopbCnTpVKhRNPnM/g4CIAenoW\nsHz5UmbNmtXmyorBIVNuuVx1KqoiB83MmXNZufJ4YH66ZCnHHruCG2+8up1lFYJDpvxyHZxc0s7A\nXGDfqvdFRJyb5YDWvRwy3SvLs07/ArwA3EsyLo2NoK/vDFatms/gYDLf07OAvr7u7s5yyHS3LJe3\nH4iID+ZcT92KfOoEST/N4sWXAUnwdHP/jEOms+TaRyPpMuA7EbG2keKarehBYwmHTOfJO2jWATOA\nx4DX08UREQdlqrJJHDTF55DpTHl/U+UnM9ZjXcwhY9V8eduaziHT2XJ5BEHSbem/r0h6ueb1UqPF\nWmdyyNhw3KKxpnHIdIe8H6o0G5FDxkbjoLFxc8jYWBw0Ni4OGatH3UGTPrldu2xOc8uxMnHIbMtD\ng4wsyw17PwXmR8T96fw84M8i4rAc6xutHncGt5FDZlvdNDRI3ncGvw/4Z5Jvqvwo8DlgTkS8mLXQ\nZnDQtI9DZnvdNDRIrncGR8SjaSvmWuAXwKyIeDVjjVZyDhlrxJhBI+n+mkWTSfp27kxbFW151sla\nzyEzMg8NMrp6vtdp32EWD71JEfF4pgNKs4GLgQnA5RGxaITtDgVuB06OiGuGWe9TpxZyyIytW4YG\nybuP5tNAJSJekvQ14GDgvIj4aYYCJwAPAccAm4C7gXm133iZbrcSeBW4MiK2O9F10LSOQ8aq5X1n\n8NfTkPlt4Gjg+8AlWQ4GHAasj4jHI2IzsAw4YZjt/pSk4/npjPu3JnPIWDNkCZq30n/nAEsi4npg\nx4zHmwZsqJrfmC7bStI0kvAZCjE3W9rEIWPNkiVoNqWj7J0C/DgdrDzrncX1hMbFwF+k50Xi7e/6\nthZyyFgzZRn46mRgNvCtiHhB0nuAszMebxMwvWp+OkmrptohwDJJAFOAT0raHBErane2cOHCrdO9\nvb309vZmLMeG45CxagMDAwwMDIxrHy0dJkLSRJLO4KOBJ4C7GKYzuGr7K4HrfNWpdRwyzdWJV6Jy\nuWFP0m0R8RFJr7D9qU9ExO71Hiwi3pR0FlAhubx9RUSsk3Rmuv7SDLVbkzlkmqv2sYRVq+Z37GMJ\nY2m4RaPk3OaUiFjW3JLqPr5bNE3kkGm+Tn0sIa+hPHeV1Cfpu5L+RNIOkk4EfgbMa7RYKw6HjOWt\nns7gHwAvAXcAxwKnk3xT5akRsTq/0qwVHDL58WMJb6vnEYS1Q88zpXfsPgnsExGDLahvtLp86jRO\nDpn8uTM4fU8dQXNfRBw80ny7OGjGxyFjjcoraN4ieeZoSA8w1JrJdNWpmRw0jXPI2Hjkcnk7IiY0\nXpIVjUPG2sGDk3cRh4y1i4OmSzhkrJ0cNF3AIWPt5qDpcA4ZKwIHTQdzyFhROGg6lEPGisRB04Ec\nMlY0DpoO45CxInLQdBCHjBWVg6ZDOGSsyLKMGWwFtWgRXH65Q8aKyy2aknPIWBk4aErMIWNl4aAp\nKYeMlYmDpoQuuMAhY+XioCkZX12yMnLQlIhDxsrKQVMSDhkrMwdNCThkrOwcNAXnkLFO4KApMIdM\nZ6lUKsycOZeZM+dSqVTaXU5LNfzd2+3W6V+34pDpLJVKhRNPnM/g4CIg+dbK5cuXlvIL5XL5Xqei\n6uSgcch0npkz57Jy5fHA/HTJUo49dgU33nh1O8tqSCNB41OngnHIdK+OPrWKiFK+ktI7y6JFETNm\nRGzc2O5KrNn6+/ujp2dqwFUBV0VPz9To7++ve32RpJ+9TJ9XnzoVhFsyna9SqbB48WUA9PWdsU3/\nTJlOrXL5SlzLn0OmO8yaNauUnb/N4KBpM4eMQdLCWbVqPoODyXxPzwL6+pa2t6gm8qlTGzlkrNpo\np1ZF4svbJeKQsbIqxeVtSbMlPSjpYUkLhll/mqQ1ktZKuk3SQa2uMW8OGes2LW3RSJoAPAQcA2wC\n7gbmRcS6qm2OBH4eES9Kmg0sjIgjhtlXKVs0DhkruzK0aA4D1kfE4xGxGVgGnFC9QUTcHhEvprN3\nAnu1uMbcOGSsW7U6aKYBG6rmN6bLRvJ54IZcK2oRh4x1s1Zf3q77XEfSx4E/BD6SXzmt4ZCxbtfq\noNkETK+an07SqtlG2gG8BJgdEc+PtLOFCxdune7t7aW3t7dZdTaNQ8bKbmBggIGBgXHto9WdwRNJ\nOoOPBp4A7mL7zuC9gZuBz0TEHaPsq/CdwQ4Z60SFfwQhIt6UdBZQASYAV0TEOklnpusvBb4OvBO4\nRBLA5og4rJV1NoNDxuxtvmEvBw4Z62RluLzd8Rwy1mydME6NWzRN5JCxZiviEKB+1qmNHDKWhyKO\nU+NTpzZxyJiNzuPRjJNDxvLUKePU+NRpHBwy1gpFG6fGfTQt5JCxbuU+mhZxyJhl46DJyCFjlp2D\nJgOHjFljHDR1csiYNc5BUweHjNn4OGjG4JAxGz8HzSgcMlZEZXzI0vfRjMAhY0VUhIcsfcNekzhk\nrKiK8JClb9hrAoeMWfP5ocoqDhkrurI+ZOlTp5RDxsqi3Q9Zuo+mQQ4Zs/q5j6YBDhmz/HV10Dhk\nrOzKck9N1546OWSs7Np1T437aOrkkLFO0K57atxHUweHjHWjtp9iRUQpX0np2SxaFDFjRsTGjZnf\nalY4/f390dMzNeCqgKuip2dq9Pf3N7xdvdLPXqbPa9ecOrklY52onntqtj/F+nMmT76WQw75UEP3\n4biPZgQOGetm2wZNBfgMcCHQWAdyI0HT8Y8gOGSs22372ML3SEImad0MDsLixZflfqWqozuDHTJm\nMGvWLJYvT65ITZ78dFtq6NhTJ4eM2faace+N+2hSDhmzkY33oUwHDQ4Zs3o1GjhdHzQOGbP6jOcU\nqquDxiFjVr/xPL5QikcQJM2W9KCkhyUtGGGbv03Xr5F08Fj7dMiYFVtLg0bSBOA7wGzgQGCepANq\ntjkOmBER+wFnAJeMts+yhMzAwEC7S8isbDWXrV5oX819fWew445fAo4EjmTHHb9EX98ZuR2v1S2a\nw4D1EfF4RGwGlgEn1GxzPLAUICLuBPaUNHW4nZUlZMAfglYoW73Q7ponAX+UviaNufXQg5mNaHXQ\nTAM2VM1vTJeNtc1ew+2sLCFjVjSLF1/GG298i6SPZj5vvPGtrVeghjPUeZz062TX6qCpt+e5tqNp\n2Pc5ZMxaY/Hiy9IrVPPH3HY4Lb3qJOkIYGFEzE7nvwJsiYhFVdt8DxiIiGXp/IPAxyLiqZp9lfNy\nmVkHKPpDlfcA+0naF3gCOAWYV7PNCuAsYFkaTC/Uhgxk/0HNrH1aGjQR8aaks0ieVZ8AXBER6ySd\nma6/NCJukHScpPXAfwB/0Moazaz5SnvDnpmVR6GHicjj5r68jVWzpNPSWtdKuk3SQe2os6qeMX/H\n6XaHSnpT0kmtrG+EWur5/6JX0n2SHpA00OISh6tnrP8vpkjql7Q6rfn0NpRZXc/3JT0l6f5Rtqn/\ns5d17M9WvUhOrdYD+5Jc5F8NHFCzzXHADen04cAdJaj5SGCPdHp2O2uup96q7W4GrgfmluB3vCfw\nM2CvdH5KCWpeCJw/VC/wLDCxjTV/FDgYuH+E9Zk+e0Vu0TT15r4WGbPmiLg9Il5MZ+9khHuEWqSe\n3zHAnwL/DLRn1KRt1VPzqcDVEbERICKeaXGNteqp+Ulg93R6d+DZiHizhTVuIyJuBZ4fZZNMn70i\nB01Tb+5rkXpqrvZ54IZcKxrdmPVKmkbyoRh6FKTdnXr1/I73AyZL+omkeyR9tmXVDa+empcA/0XS\nE8Aa4Istqq1RmT57RR4zuKk397VI3ceW9HHgD4GP5FfOmOqp92LgLyIiJIntf9+tVk/Nk4DfBI4G\ndgFul3RHRDyca2Ujq6fmrwKrI6JX0vuBlZI+FBEv51zbeNT92Sty0GwCplfNTydJzdG22Std1i71\n1EzaAbwEmB0RozVP81ZPvYeQ3NMESd/BJyVtjogVrSlxO/XUvAF4JiIGgUFJtwAfAtoVNPXUfBTw\nDYCIeETSY8AHSO49K6Jsn712dpKN0Rk1EXiEpANtR8buDD6C9ncG11Pz3iQdg0eU4Xdcs/2VwElF\nrxnYH/g/JJ2wuwD3AwcWvOaLgHPS6akkQTS5zb/rfamvM3jMz15hWzRRwpv76qkZ+DrwTuCStJWw\nOSIOK3C9hVLn/xcPSuoH1gJbgCUR8fMi1wx8E7hS0hqSvtOzI+K5dtUs6UfAx4ApkjYA55A+4t3I\nZ8837JlZ7op81cnMOoSDxsxy56Axs9w5aMwsdw4aM8udg8bMcuegMbPcOWjMLHcOmi4l6T9JWiZp\nffqE848l7ZdxH3tI+uNR1t+WcX8LJfVlec8o+5KkMyV9IX1I0drIQdOF0qewlwM3R8SMiPgt4Csk\nz9hk8U7gT0ZaGRFZn0xv5m3qXyQZ7+cnwKeauF9rgIOmO30ceCMitn5jWESsjYhVkr4s6f70tXVM\nFEnvSFs9q9N1JwPnA+9Ph8xcVHsQSa+k/+4raZ2ky9JhKiuSdk7X/aWkhyTdSvK08tB7PyPpznTf\n35O0Q7r80HToyJ3Smh6QdGDNcScBcyJiNbAPsEcTf3fWgMI+VGm5+iBwb+1CSYcAp5OMCLcDcKek\nf0s/sLOBTRHxO+m2u5O0GD4YESONF1vdQpkBnBIRZ0j638BcJd/ZdQrJEA6TgJ8C9yj5PvaTgaMi\n4i1JfwecBvx9RNwtaQVwHtCTLqt9YPITwMuS5gO/C9yU5Zdjzeeg6U4jnaL8NnBNJOO4IOkakrFj\nV5M8CX2hpL8Brk9bP5MzHPOxiFibTt9LMgTBlPR4rwGvpQEikqA4hCR0IAmUX1Xt61yScVoGSYYZ\nrXUkyRPS10v6NHB7hjotBz516k4/I/kg1wq2HTVN6TIiGZ3uYJKxXc6T9DWy9am8XjX9Fm//kas9\n3tC/SyPi4PS1f0ScW7XdFOAdwK4kIVTrPcCjknYC3hMRqyWdIOm9Geq1JnLQdKGIuBnYSdIXhpal\no/6tBn5PUo+kdwC/B9yarn8P8FpE/ANwIUnovAzsNo5SbkmPt7Ok3YA5JOF1E/ApSe9Kjz1Z0t5V\n77sU+Cvgh8B2fUMk3yDwOnAScFE6aPZ82j8MadfyqVP3OhG4WMl3DL0GPAb8GXAVcFe6zZKIWJNO\n/wbwLUlbgM3AH0XEc0q+m+p+ktHWar+vKEaYBoiIuC/tr1kD/PvQcdNBof4KuDHtBN5McnXrl5I+\nB7weEcvSdf9XUm9EDFTt+0ckIfNKRFwCkA4oZW3iga+sK0g6B7g8Ito5pnTX8qmTdTxJ7ya5dP7x\ndtfSrdyiMbPcuUVjZrlz0JhZ7hw0ZpY7B42Z5c5BY2a5c9CYWe4cNGaWOweNmeXOQWNmufv/hsCP\nui4SwMIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa51783c490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta1 = (P[:,3]-MinCost)/(MaxCost-MinCost)\n",
    "theta2 = (P[:,4]-MinRisk)/(MaxRisk-MinRisk)\n",
    "x_line,y_line = 0.01*np.array(range(0,100,1)),0.01*np.array(range(0,100,1))\n",
    "plt.scatter(theta1,theta2)\n",
    "plt.axis([0,1,0,1])\n",
    "plt.plot(x_line, y_line)\n",
    "plt.xlabel('Cost index $\\\\theta_1$')\n",
    "plt.ylabel('Risk index $\\\\theta_2$')\n",
    "plt.title(r'Efficient frontier')\n",
    "plt.axes().set_aspect('equal')\n",
    "#plt.savefig('EF.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A new data array is formed with $\\theta_1$ and $\\theta_2$ to store results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "new_col1 = theta1[...,None]\n",
    "new_col2 = theta2[...,None]\n",
    "PA = np.concatenate((P, new_col1, new_col2), 1)\n",
    "#np.savetxt('MillerTable.txt',PA,fmt='%.2f',delimiter =';')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best policy minimizes distance to the ideal point $(0,0)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def find_policy(P,col1=5,col2=6):\n",
    "    \"\"\"Find the policy tha minimizes distance to the ideal in the space col1-col2\"\"\"\n",
    "    min_dist = 1\n",
    "    best = P[0,:]\n",
    "    for element in range(len(P)):\n",
    "        dist = np.sqrt(P[element,col1]**2+P[element,col2]**2)\n",
    "        if dist < min_dist:\n",
    "            best = P[element,:]\n",
    "            min_dist = dist\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  8.20000000e+05   5.30000000e+05   2.30000000e+05   8.44900947e+02\n",
      "   3.67882650e+03   1.92344732e-01   4.96242761e-01]\n"
     ]
    }
   ],
   "source": [
    "best = find_policy(PA)\n",
    "print(best)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Evaluate the best policy over the test set and the training set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test = policy(50,5,0.05,0.001,0.00027,0.3,best[0],best[1],best[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(652.30655200000012, 3067.421072463615)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(test.cost_calc(f_test,best[1])),np.std(test.cost_calc(f_test,best[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(844.9009473750001, 3678.8265045417697)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(test.cost_calc(f_train,best[1])),np.std(test.cost_calc(f_train,best[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Considering the test set error as a new measure of risk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minimizing cost-risk square error and cost-risk understimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Cost-risk square error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we search for feasible policies in the training set and the test set. Again, initilization of $h,z,l$ is arbitrary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pol = policy(50,5,0.05,0.01,0.00027,0.3,h_miller,z_miller,l_miller)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sp1 = pol.space(vh,vz,vl,f_train,'cost_std')\n",
    "sp2 = pol.space(vh,vz,vl,f_test,'cost_std')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "S1 = transform(sp1)\n",
    "S2 = transform(sp2)\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\theta_1$ and $\\theta_2$ are obtained for the training set and the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(939.08371201250009,\n",
       " 3593.3192204624997,\n",
       " 516.48381760000007,\n",
       " 2508.8660864000003)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaxCost1 = np.max(S1[:,3])\n",
    "MinCost1 = np.min(S1[:,3])\n",
    "MaxRisk1 = np.max(S1[:,4])\n",
    "MinRisk1 = np.min(S1[:,4])\n",
    "MaxCost2 = np.max(S2[:,3])\n",
    "MinCost2 = np.min(S2[:,3])\n",
    "MaxRisk2 = np.max(S2[:,4])\n",
    "MinRisk2 = np.min(S2[:,4])\n",
    "MinCost1,MaxCost1,MinCost2,MaxCost2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_theta1 = (S1[:,3]-MinCost1)/(MaxCost1-MinCost1)\n",
    "train_theta2 = (S1[:,4]-MinRisk1)/(MaxRisk1-MinRisk1)\n",
    "test_theta1 = (S2[:,3]-MinCost2)/(MaxCost2-MinCost2)\n",
    "test_theta2 = (S2[:,4]-MinRisk2)/(MaxRisk2-MinRisk2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we obtain $\\theta_3$ as the normalized distance between predicted points ($\\theta_1^*,\\theta_2^*$) from the training set and real points ($\\theta_1,\\theta_2$) from the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.80903436572271636, 0.0011035970604048393, 1.0)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dist = np.sqrt((train_theta1-test_theta1)**2+(train_theta2-test_theta2)**2)\n",
    "MaxDist = np.max(Dist)\n",
    "MinDist = np.min(Dist)\n",
    "theta3 = (Dist - MinDist)/(MaxDist-MinDist)\n",
    "MaxDist,MinDist,np.max(theta3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\theta_3$ is added to the initial C-R space and a new efficient frontier is obtained in the space Cost-Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "new_col3 = theta3[...,None]\n",
    "S = np.concatenate((S1, new_col3), 1)\n",
    "P2=frontier(S[np.argsort(S[:,3])],col=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "new_theta1 = (P2[:,3]-np.min(P2[:,3]))/(np.max(P2[:,3])-np.min(P2[:,3]))\n",
    "new_theta3 = (P2[:,5]-np.min(P2[:,5]))/(np.max(P2[:,5])-np.min(P2[:,5]))\n",
    "col1=new_theta1[...,None]\n",
    "col3=new_theta3[...,None]\n",
    "PA = np.concatenate((P2, col1, col3), 1)\n",
    "np.savetxt('MillerTable.txt',PA,fmt='%.2f',delimiter =';')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find the most robust policy to prediction errors as the policy closer to the ideal point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  8.20000000e+05,   5.40000000e+05,   1.60000000e+05,\n",
       "         1.01215038e+03,   4.70321518e+03,   2.84824488e-02,\n",
       "         1.28074870e-01,   1.21949126e-01])"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best2=find_policy(PA,col1=6,col2=7)\n",
    "best2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAEeCAYAAABVItqsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHk5JREFUeJzt3XucHHWd7vHPI0EyIiAxwMEQYCUgsIqicvPGCCYTXQyG\nqCyim7AeYffI4jnkHKKikuMmanIMy2EXWW4m42VlLxA3spghCwwCy02WBFgTToKgCagrQZTLgAn5\nnj+qhnSauXT1dFV3dT/v16tfqeqqrvrOQD/zq19V/UoRgZlZnl7R7ALMrP05aMwsdw4aM8udg8bM\ncuegMbPcOWjMLHcOmjYkaYGkX0t6PJ2fKWmjpN9JeoukByW9p4btPC3pwLzrTff1Tknr033OKGB/\np0vqy3s/lpCvoykfSY8CewMvVry9NCLOkbQ/sA6YHBGb0/UfBv57RPyg8GKT/S8DNkbEF0dY50bg\n+xHx1zns/0Dgp8C4iNjW6O3b6MY1uwCrSwAnRcRNQyzbH9hcETJK3/tJgfXVY9ga05+BGPtfRY3x\n88lGpFc4sLLxoVMbkfQ+4AbgdekhyN8BvwN2AtZIWp+u96ikE9PpnSR9XtKG9NDqx5Impcu2SXp9\nOr2LpK9L+pmkX0q6VNL4dFm3pE2SzpX0K0mPS5qTLjsT+BhwXlrTPw9R98PA64EfpDW8UlJ/egh4\nO/As8AeS3iHpHklPSbpb0nEV2+iX9GVJt6Xb6JP02nTxj9J/n0qXHStpjqRbKz5/qKRVkjZLWifp\nIxXLlqU/7/WSngG6x/ZfqgNFhF8lewGPACcOs+x4ksOUyve2Aa+v+vwJ6fT/Au4HDk7njwAmVH8O\n+Cvg+8BrgFcDK4CvpMu6gS3AfJJQez9JOOyRLl8KfLmGn+mEivl+4FHgMJI/iPsAvwFOT+f/GHgS\n2LNi/fXAFGA8cDPw1XTZAenP8oqK7c8Bbk2ndwU2ArPTbb8F+DVwWLp8GfAUcFw6v0uz/x8o28st\nmnIS8H1Jv6l4fbJiWRb/FTg/ItYDRMT9EfHkDjtLDl0+BZwbEU9FxDPAV0m+7IO2kITJixHxQ+AZ\n4A1VNWcRwLKIWBvJYco04KGI+G5EbIuIq0n6omZUrL80IjZExPPAP5AERi37Pgl4JCJ6022vBq4F\nPlKxzvcj4g6AiHgh48/S8dxHU04BnBxD99FktR/w8Cjr7AW8Crg37S6B5Mtb+Ydqc+zYb/EcSctn\nLDZWTL8O+HnV8p+l7w/6ZcX0QIb9HwAcI+k3Fe+NA76VTgewqcZt2RAcNLaR5HBjpM7iJ0i+uIdH\nxC/q2Ee9nbiVn3sMOKVq+QHADxuw/58Dt0TEtAy1WQY+dCqvhpxBAa4E/lLSFCWOkDShcoW0pXIF\ncJGkvQAkTZJU6xfzVySdvVlV/ozXA4dIOk3SOEmnAocC1w2zfqVfk/TRHDTM8n9Jt/1xSTunr6Mk\nHTrKdq1GDpry+kF6FmfwdU3Fsuq/4CP9Rb+QpD/jBuC3JIEyfojPzQM2AHdK+i2wCjikxn1cBRye\n9iVdO8J61V7aZtpvdBIwl6SF9T9JTvE/OdT66XSkn30OWAjcLulJScdULX+apA/oj0laTr8g6YN6\nZfW2rD6FXrAn6ZvAHwH/GRFvGmadi0nOWjwHzImI+wor0MxyUXSLZikwfbiFkj4ATImIg4EzgUuL\nKszM8lNo0ETErSTXQgxnBtCbrnsX8BpJ+xRRm5nlp9X6aCax4ynNTSSnX82sxFotaODlPfzuhDMr\nuVa7juYxYHLF/H7pey8jyQFk1iQRkemUf6u1aFYAfwIg6VjgqYj41fCr786CBQuafh9HLa8LLrig\n6TW0e81lq7dMNa9evX26HoW2aCR9j+Smv4mSNgIXADsDRMRlEXG9pA9I2kByU94ZI21vwYLzOP/8\n8/Mu26yjLV4MV10Fa9bA+PGjrz+UQoMmIk6rYZ2za92eQ8YsX4sXwxVXQH9//SEDrXfo1La6u7ub\nXUJmZau5bPVCa9dcGTKTJo1tW6UdylNSlLV2s1Y3UshIIkreGWxmTdbIlswgB42ZvSSPkIGSB820\nabPo6/MTM8waIa+QgZL30cBE4EVmz57BsmXLml2SWWllCZkO7KP5OvBX9PYuZ+HChc0uxqyU8mzJ\nDCp5i2aw9l4mTPhLNm/e0NSazMqmnpDpwBaNmdWriJbMoFa7qTKj3vTfczj33POaWolZmRQZMlDy\nQ6dx4/amq2s88+ad6dsRzGo01pCp59Cp1EFT1trNmqURLRn30ZjZsIo+XKrkoDHrAM0MGXDQmLW9\nZocMOGjM2lorhAyUPGgWLlzIa187hde+doqvDDar0iohAyU/6wS7AHul7/yaBQu+6NPcZuR8g2Sn\nnd6G3YGL03fOYfz4cQwMbG5mWWZNl3dLpp6gKfmVwRcDs1+ae/75uc0rxawFtNLhUqVS99FU6+ra\npdklmDVNq4YMlL5Fc84O0+ef7/udrDO1cshAyftoFixYwIUXLgXg3HPPcEewdaTCb5DstM7gstZu\n1ijNaMn4XiezDtLqh0uVHDRmJVSmkAEHjVnplC1kwEFjViplDBlw0JiVRllDBtogaPr6+pg2bZYf\nJmdtrcwhAyU/vb1y5UpmzpzNwMAiALq65rF8eS89PT1Nrs6scVotZDruOpqpU09h1aoZbL/fqZep\nU1dwww3XNLM0s4ZptZABX0dj1lZaMWTqVep7nebOPZPbbpvNwEAy39U1j7lze0f+kFkJtFPIQMkP\nnSKCvr4+liy5HEiCx/0zVnatHjId10dT1trNhtPqIQMl6aORNF3SOknrJc0bYvlESSslrZb0oKQ5\nRddo1gyLFrV+yNSr0BaNpJ2Ah4D3AY8B9wCnRcTainXmA7tExOckTUzX3ycitlZtyy0aaxuLFsGV\nV5YjZMrQojka2BARj0bEFuBq4OSqdX5BMhgw6b+bq0PGrJ2UKWTqVfRZp0nAxor5TcAxVetcAdwk\n6XFgN+CjBdVmVrhOCBkovkVTy7HO54HVEfE64C3AJZJ2y7css+ItXtwZIQPFt2geAyZXzE8madVU\negewECAiHpb0CPAG4MfVG5s/f/5L093d3XR3dze2WrOclOHs0qD+/n76+/vHtI2iO4PHkXTungg8\nDtzNyzuDLwR+GxH/W9I+wL3AERHxZNW23BlspVSmkBlKyz/XKSK2Sjob6AN2Aq6KiLWSzkqXXwZ8\nBVgqaQ3Jod151SFjVlZlD5l6+YI9s4K0S8iU4fS2WUdql5Cpl4PGLGedHjLgoDHLlUMm4aAxy4lD\nZjsHjVkOHDI7ctCYNZhD5uUcNGYN5JAZmoPGrEEcMsNz0Jg1gENmZA4aszFyyIzOQWM2Bg6Z2jho\nzOrkkKmdg8asDg6ZbBw0Zhk5ZLJz0Jhl4JCpj4PGrEYOmfo5aMxq4JAZGweN2SgcMmPnoDEbgUOm\nMRw0ZsNwyDSOg8ZsCA6ZxnLQmFVxyDSeg8asgkMmHw4as5RDJj8OGjMcMnlz0FjHc8jkz0FjHc0h\nUwwHjXUsh0xxHDTWkRwyxXLQWMdxyBRvzEEj6RBJuzSiGLO8OWSaY1w9H5L0FWBv4B5gCvAC8IUG\n1mXWcA6Z5qkraIAbgPXAbsB3gLc2rCKzHDhkmqveQ6cngKMiYh3waWBb40oyayyHTPMpIkZeQToI\n+DCwB7ARuCMiVhdQ24gkxWi1mzlkGk8SEaEsn6mlRXMCcA1wI3AisEDSvZJOr6NGJE2XtE7Seknz\nhlmnW9J9kh6U1F/PfswcMq2jlqB5BfDqiLgR+EFEnAS8A9gm6c+z7EzSTsDfANOBw4HTJB1Wtc5r\ngEuAD0bEG0laU2aZOGRaSy1BcznQLelfgQ9JOgk4CLgbeHXG/R0NbIiIRyNiC3A1cHLVOh8DromI\nTQAR8UTGfViHc8i0nlHPOqUdIRdJuoTkMOo44IPAZuB7Gfc3iaSfZ9Am4JiqdQ4GdpZ0M8lZrf8b\nEd/OuB/rUA6Z1lTz6e20BdKXvupVS+/tziSny08EXgXcIenOiFg/hv1aB3DItK56r6Op12PA5Ir5\nySStmkobgSciYgAYkPQj4M0k1+3sYP78+S9Nd3d3093d3eByrSwcMvnp7++nv79/TNsY9fR2I0ka\nBzxE0lp5nKSf57SIWFuxzqEkHcY9wC7AXcCpEfGTqm359LYBDpmi5XV6e3Djhw/xXneWnUXEVuBs\nksOvnwB/HxFrJZ0l6ax0nXXASuB+kpC5ojpkzAY5ZMqh5haNpAeBbwOLgS5gEcnVwcfmV96I9bhF\n0+EcMs2Ra4uG5OzQZOAOkkOeX5BcT2NWOIdMuWQJmq3AAElrZjzw04jwPU5WOIdM+WQJmruB54G3\nA+8GPibpH3OpymwYDplyytJHc1RE3FP13ieadTGd+2g6j0OmNeTdR3OvpE9I+lK6s/2B/5dlZ2b1\ncsiUW5ag+QbJ7QcfS+efIbn50SxXDpnyy3Jl8DERcaSk+wAi4klJO+dUlxngkGkXWVo0v0+HeQBA\n0l54ZD3LkUOmfWQJmr8GlgN7p4OT3w58NZeqrOM5ZNpLpnud0kGqTgAE3Fh5j1LRfNapfTlkWls9\nZ51qGTN4bsVskITM4DQRcWGWHTaKg6Y9OWRaXz1BU0tn8G4kofIG4ChgBUnYnERyEZ9ZQzhk2leW\nC/ZuBT4QEU+n87sB10fEu3Osb6R63KJpIw6Z8sj7gr29gS0V81vS98zGxCHT/rJcR/Mt4G5J15Ic\nOn0I6M2lKusYDpnOkPWs09tIbqgM4EcRcV9ehdVQiw+dSs4hU065nHVqVQ6acnPIlFdeZ50GNz4e\nmAUcWPG5iIgvZ9mhmUOm82Tpo/ln4CngXpJxacwyc8h0pixBMykienKrxNqeQ6ZzZTm9/W+Sjsit\nEmtrDpnOluWCvbXAFOAR4IX07YiIpoSPO4PLwyHTXnLtDAben7EeM4eMAT69bTlyyLSnXG5BkHR7\n+u8zkp6uev2u3mKtvTlkrNKoQRMR70z/fXVE7Fb12j3/EkfX19fHtGmzmDZtFn19fc0up+M5ZKxa\n6Q+d+vr6mDlzNgMDiwDo6prH8uW99PT4THwzOGTaX0fegjBt2ixWrZoBzE6X9DJ16gpuuOGaZpbX\nkRwynSG3YSKUmFxfWdYJHDI2kiynt38IvDGvQuo1d+6Z3HbbbAYGkvmurnnMnevRK4rkkLHRZLlg\nrxe4JCJaYvjOytPbfX19LFlyOZAEj/tniuOQ6Ty59tFIeojkyuCfAc+mb/vK4A7mkOlMeV8ZPNhM\nGPx2Z9qRtReHjGWRdYS9t7B9hL1bI2JNXoXVUItbNE3ikOlsuQ5OLukzwHeAvYB9gO9IOidbiVZ2\nDhmrR5Y+mgeAYyPi2XR+V+DOiHhTjvWNVM+QLRp3DOfHIWOQ/+NWALYNM10zSdMlrZO0XtK8EdY7\nStJWSafUuu3Bq4RXrZrBqlUzmDlztm9JaBCHjI1Fls7gpcBdVY9b+WaWnUnaCfgb4H3AY8A9klZU\nP8M7XW8RsJIMnc5Lllye3oqQXCU8MJC851bN2DhkbKxqChpJAv4JuAV4F0ln8Jw6HrdyNLAhIh5N\nt3s1cDKwtmq9v0j3d1TG7VuDOWSsEbK0aK6PiDeSDE5er0nAxor5TcAxlStImkQSPieQBE3Np5Z8\nlXBjOWSsUWrqo0l7Xe+VdPQY91dLaFwEfDbdp8hw6NTT08Py5clNlVOnrvBd3GPgkLFGytKiORb4\nuKSxXBn8GFB5c+ZkklZNpbcBVydHa0wE3i9pS0SsqN7Y/PnzX5ru7u6mu7ubnp4eh8sYOWSsUn9/\nP/39/WPaRk2nt9M+mncDP69eNtjfUtPOpHHAQ8CJwOPA3cBp1Z3BFesvBX4QEdcOscwX7OXAIWOj\nyfsWhG+kfTR1i4itks4G+oCdgKsiYq2ks9Lll41l+zY2DhnLS1vcvW1j55CxWvnubauLQ8ayKOru\n7Uql+Kb7toThOWSsCLU8buU8eKnT96iIeHTwBZyVb3lj59sShueQsaKMeugk6b6IOLJ6eqj5ItV6\n6OTBy4fmkLF6FXFTpbUBh4wVLUsfTSn5toQdOWSsGWo5dHoReC6d7QIGKhZ3RURTwirLWSd3Bicc\nMtYIHfkAOauNQ8YaxX00NiSHjDWbg6bNOWSsFTho2phDxlqFg6ZNOWSslTho2pBDxlqNg6bNOGSs\nFTlo2ohDxlqVg6ZNOGSslbX9LQidYNEiuPJKh4y1LrdoSs4hY2XgoCkxh4yVhYOmpBwyViYOmhJa\nvNghY+XioCkZn12yMuqYoOnr62PatFlMmzartGMGO2SsrDpiPJrBAcoHBhYBySh7ZXsut0PGWoUH\nvhpG2Qcod8hYK/HAV23IIWPtoCOuDC7rAOUOGWsXHXHoBOUboNwhY63KfTRtwiFjrcx9NG3AIWPt\nyEHTQhwy1q46Omha6SI+h4y1s47to9l+Ed+7gLuA55g9+4MsW7asUSXWzCFjZeI+mgyWLLk8DZlV\nwALgQnp7l7Nw4cJC63DIWCfo2KBJ3AVcTHLF8GzgYi68cGlhe3fIWKfo2KCZO/dM4Lmm7d8hY52k\n8KCRNF3SOknrJc0bYvnpktZIul/S7ZKOyKOOnp4eZs/+IHAO0Ju+zuHcc8/IY3c7cMhYx4mIwl7A\nTsAG4EBgZ2A1cFjVOscBe6TT04E7h9lWNMKCBQtiwoSDYsKEg2LBggUN2eZIFi2KmDIlYtOm3Hdl\nlov0u5fpu1/oWSdJxwEXRMT0dP6zaWJ8bZj19wQeiIj9hlgWRdbeCG7JWDsow1mnScDGivlN6XvD\n+SRwfa4VFcQhY52s6Lu3a26CSHov8KfAO/MrpxgOGet0RQfNY8DkivnJJK2aHaQdwFcA0yPiN8Nt\nbP78+S9Nd3d3093d3ag6G8YhY2XX399Pf3//mLZRdB/NOOAh4ETgceBu4LSIWFuxzv7ATcDHI+LO\nEbbV8n00DhlrR/X00RTaoomIrZLOBvpIzkBdFRFrJZ2VLr8M+BKwJ3CpJIAtEXF0kXU2gkPGbLuO\nvdcpTw4Za2dlOOvU0hpxN7dDxuzl3KJJ7Xg39y3ANvbddyJLl15c87CfDhnrBB7KcwySR7II+CHw\nKuDrAIwbN5frrvvuqGHjkLFO0fKdwa3vLuAI4M8YfAbU1q3JkBIjBY1Dxmxk7qNJ1Xs3t0PGbHQ+\ndKowZ84cenv/kVoPnRwy1oncR9MACxcu5Gtfu4Tnn3+RAw54HZdc8jWHjFkFB01BHDLWyXwdTQEc\nMmbZOWgycMiY1cdBM4rBq4UPOeRbLF78JPvu+ynOOKP5z4EyKxP30Yxg+9XC1wEHAO8GPgdAV9c8\nli/vrfmqYbN24c7gBkuuFv4c8HbgU8C7GLyQD3qZOnUFN9xwTa41mLUadwY32KOPngwcks492cxS\nzErNQTOMxYvh2WdnMX78e0gexfIHVD6apatrXno18Xat9Cxvs1biQ6chVJ5devDBPpYsuRyA449/\nK7fc8u9AcstCZf/M9v6cRYD7cKx9uY+mAeo9hZ3058zAfTjW7txHM0a+TsYsHx4mIjXWkJk790xu\nu202AwPJfNKH09vQGs3KyodONK4l09e3vT+nug/HrF24j6YOPlwyy8Z9NBk5ZMyK0bFB45AxK05H\nBo1DxqxYHRc0Dhmz4nVU0DhkzJqjY4LGIWPWPB0RNA4Zs+Zq+6BxyJg1X1sHjUPGrDW0bdA4ZMxa\nR1sGjUPGrLW0XdA4ZMxaT1sFTauGjIf4tE7XNkHTyiEzc+ZsVq2awapVM5g5c/aoYeNgsrYTEYW+\ngOnAOmA9MG+YdS5Ol68BjhxmnRi0aFHElCkRmzZFy5k69ZSAZQGRvpbF1KmnDLv+ypUro6trn/Qz\ny6Kra59YuXLliPtYuXJlTJ16Skydesqo62aR13at3NLvXrbvfdYPjOUF7ARsAA4EdgZWA4dVrfMB\n4Pp0+hjgzmG2FRGtHTIRlUFzc01BU0Qw1WLlypXxylfu2fDt5unmm29udgmZlbHmeoKm6EOno4EN\nEfFoRGwBrgZOrlpnBskzTYiIu4DXSNpnqI216uFSpblzz6Srax5wEcM9pmUsliy5PH3ywmwgeQrD\n4Ch/Y93u73//noZvN0/9/f3NLiGzMtU8eEhfj6KDZhKwsWJ+U/reaOvsN9TGdt21tUMGoKenh+XL\ne3n96x9i6tQVoz6CZXswDf/8KLOiVfY11qPowclrHXuzepjAIT/36U+PrZii9PT08IlPnMr8+fNr\nWnf58t6KsYdHD6Y8BkWfO/dMbr55Flu39jZ0u1ZOO7ac52T+fKFjBks6FpgfEdPT+c8B2yJiUcU6\nfwv0R8TV6fw64PiI+FXVtso52LFZG4iMYwYX3aL5MXCwpAOBx4FTgdOq1lkBnA1cnQbTU9UhA9l/\nUDNrnkKDJiK2Sjob6CM5A3VVRKyVdFa6/LKIuF7SByRtAJ4FziiyRjNrvNI+bsXMyqOlrwyWNF3S\nOknrJc0bZp2L0+VrJB1ZdI1D1DNizZJOT2u9X9Ltko5oRp0V9Yz6O07XO0rSVkmnFFnfMLXU8v9F\nt6T7JD0oqb/gEoeqZ7T/LyZKWilpdVrznCaUWVnPNyX9StIDI6xT+3cv64U3Rb1o4MV9LVbzccAe\n6fT0ZtZcS70V690EXAfMKsHv+DXAfwD7pfMTS1DzfOCrg/UCm4FxTaz53cCRwAPDLM/03WvlFk1D\nL+4ryKg1R8QdEfHbdPYuhrlGqCC1/I4B/gL4J+DXRRY3jFpq/hhwTURsAoiIJwqusVotNf8C2D2d\n3h3YHBFbC6xxBxFxK/CbEVbJ9N1r5aBp6MV9Baml5kqfBK7PtaKRjVqvpEkkX4pL07ea3alXy+/4\nYGCCpJsl/VjSJwqrbmi11HwF8IeSHie5x+8zBdVWr0zfvaJPb2fR0Iv7ClLzviW9F/hT4J35lTOq\nWuq9CPhsRIQk8fLfd9FqqXln4K3AicCrgDsk3RkR63OtbHi11Px5YHVEdEs6CFgl6c0R8XTOtY1F\nzd+9Vg6ax4DJFfOTSVJzpHX2S99rllpqJu0AvgKYHhEjNU/zVku9byO5pgmSvoP3S9oSESuKKfFl\naql5I/BERAwAA5J+BLyZZESAZqil5ncACwEi4mFJjwBvILn2rBVl++41s5NslM6occDDJB1or2T0\nzuBjaX5ncC0170/SMXhsGX7HVesvBU5p9ZqBQ4F/JemEfRXwAHB4i9d8IXBBOr0PSRBNaPLv+kBq\n6wwe9bvXsi2aKOHFfbXUDHwJ2BO4NG0lbImIo1u43pZS4/8X6yStBO4HtgFXRMRPWrlm4CvAUklr\nSPpOz4uIJ5tVs6TvAccDEyVtBC4gOSSt67vnC/bMLHetfNbJzNqEg8bMcuegMbPcOWjMLHcOGjPL\nnYPGzHLnoDGz3DlozCx3DpoOJem/SLpa0ob0Dud/kXRwxm3sIenPR1h+e8btzZc0N8tnRtiWJJ0l\n6VPpTYrWRA6aDpTehb0cuCkipkTE24HPkdxjk8WewH8bbmFEZL0zvZGXqX+GZLyfm4EPN3C7VgcH\nTWd6L/D7iHjp0ZMRcX9E3CbpXEkPpK+XxkSRtGva6lmdLvso8FXgoHTIzEXVO5H0TPrvgZLWSro8\nHaayT9L4dNn5kh6SdCvJ3cqDn/24pLvSbf+tpFek7x+VDh25S1rTg5IOr9rvzsBJEbEaOADYo4G/\nO6tDy95Uabl6I3Bv9ZuS3kbydLCjSf4I3SXplvQLOx14LCL+KF13d5IWwxsjYrjxYitbKFOAUyPi\nTEl/D8xS8syuU0mGcNgZ+Hfgx5IOAz4KvCMiXpT0DeB04NsRcY+kFcACoCt9r/qGyROApyXNBj4I\n3Jjll2ON56DpTMMdorwLuDaScVyQdC3J2LGrSe6E/rqkrwHXpa2fCRn2+UhE3J9O30syBMHEdH/P\nA8+nASKSoHgbSehAEii/rNjWl0nGaRkgGWa02nEkd0hfJ+kjwB0Z6rQc+NCpM/0HyRe5WrDjqGlK\n3yOS0emOJBnbZYGkL5KtT+WFiukX2f5Hrnp/g//2RsSR6evQiPhyxXoTgV2BV5OEULV9gZ9K2gXY\nNyJWSzpZ0usy1GsN5KDpQBFxE7CLpE8NvpeO+rca+JCkLkm7Ah8Cbk2X7ws8HxHfBb5OEjpPA7uN\noZQfpfsbL2k34CSS8LoR+LCkvdJ9T5C0f8XnLgO+APwd8LK+IZInCLwAnAJcmA6aPZvmD0PasXzo\n1LlmAhcpecbQ88AjwP8AlgF3p+tcERFr0uk3Af9H0jZgC/BnEfGkkmdTPUAy2lr184pimGmAiIj7\n0v6aNcB/Du43HRTqC8ANaSfwFpKzWz+X9CfACxFxdbrs3yR1R0R/xba/RxIyz0TEpQDpgFLWJB74\nyjqCpAuAKyOimWNKdywfOlnbk7Q3yanz9za7lk7lFo2Z5c4tGjPLnYPGzHLnoDGz3DlozCx3Dhoz\ny52Dxsxy56Axs9w5aMwsdw4aM8vd/wcgh0qUVL/arAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e598c4350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_line,y_line = 0.01*np.array(range(1,100,1)),0.01*np.array(range(1,100,1))\n",
    "plt.scatter(new_theta1,new_theta3)\n",
    "plt.plot(x_line, y_line)\n",
    "plt.axis([0,1,0,1])\n",
    "plt.xlabel('Cost index $\\\\theta_1$')\n",
    "plt.ylabel('Error index $\\\\theta_3$')\n",
    "plt.title(r'Efficient frontier')\n",
    "plt.axes().set_aspect('equal')\n",
    "plt.savefig('EF.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test2 = policy(50,5,0.05,0.001,0.00027,0.3,best2[0],best2[1],best2[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((652.30655200000012, 3067.421072463615),\n",
       " (475.81300040000002, 2734.1533723922753))"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.mean(test.cost_calc(f_test,best[1])),np.std(test.cost_calc(f_test,best[1]))\n",
    "b = np.mean(test2.cost_calc(f_test,best2[1])),np.std(test2.cost_calc(f_test,best2[1]))\n",
    "a,b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 3.1 Cost-risk under estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cost-risk underestimation, real(test)-pred(train) Manhattan distance, we can use weights for risk preferences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1216362176267416, -0.60402388820089636)"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Man_Dist = (train_theta1-test_theta1)+(train_theta2-test_theta2)\n",
    "MaxDist = np.max(Man_Dist)\n",
    "MinDist = np.min(Man_Dist)\n",
    "theta3U = (Man_Dist - MinDist)/(MaxDist-MinDist)\n",
    "MaxDist,MinDist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  7.30000000e+05,   5.10000000e+05,   1.30000000e+05,\n",
       "         1.01862009e+03,   4.92593692e+03,   2.99658315e-02,\n",
       "         1.81584203e-01,   1.11063806e-02])"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_col3 = theta3U[...,None]\n",
    "SU = np.concatenate((S1, new_col1, new_col2, new_col3), 1)\n",
    "P2U=frontier(SU[np.argsort(SU[:,5])],col=7)\n",
    "best2U=find_policy(P2U,col1=5,col2=7)\n",
    "best2U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEeCAYAAAAuBhhgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGyxJREFUeJzt3XuYVPWd5/H3RzDKqtEQ74gBb1ETL0TxkqwzHRUaXQMi\nRkWdoJNVks1tN8xqZpJMeJzsKLMh4cmYcdQ4wpgZcScRBxO1JZqeoFFRVNQNGIgSERNXAY0aNCjf\n/eOchqLsS1V3Vf1OVX9ez1OP5/Lrc77Vdn34nd+5lCICM7OUtktdgJmZg8jMknMQmVlyDiIzS85B\nZGbJOYjMLDkH0SAk6ZuSXpL0Qj4/WdIaSb+XdLSkpyT9SQXbeU3SqHrXm+/rY5JW5vuc2ID9XSCp\no977sYx8HVHrkbQa2BN4p2TxjRHxRUn7AyuAkRGxLm//a+C/R8TtDS822/9cYE1EfL2XNvcAt0XE\n39dh/6OAZ4ChEbG51tu3vg1NXYDVRQBnRMS93azbH1hXEkLKl/2ygfX1R4815u+BGPi/qhrgz2cb\nkbZzoFXHh2aDiKRTgbuBffNDnH8Ffg8MAZZJWpm3Wy3plHx6iKS/krQqP3R7RNKIfN1mSQfk0ztI\n+pak30j6naRrJO2Yr2uT9LykL0t6UdILki7K110KnA9cltf0793U/WvgAOD2vIb3SOrMDzHvB94A\nRkv6qKSHJb0iaYmkE0u20SnpCkn35dvokPT+fPXP8/++kq87QdJFkhaX/PyhkhZJWidphaRPlqyb\nm7/fOyS9DrQN7P/UIBQRfrXYC3gWOKWHdX9KdhhUumwzcEDZz5+cT/9P4Ang4Hz+SGB4+c8B3wFu\nA3YDdgYWAn+br2sDNgEzyULvNLLw2DVffyNwRQXv6eSS+U5gNXAY2T+oewEbgAvy+fOA9cD7Stqv\nBA4CdgR+BlyZr/tA/l62K9n+RcDifHonYA0wLd/20cBLwGH5+rnAK8CJ+fwOqf8Gmu3lHlFrEnCb\npA0lr0+XrKvGfwW+GhErASLiiYhYv83OskOjS4AvR8QrEfE6cCVZGHTZRBY270TEncDrwAfLaq5G\nAHMjYnlkh0Hjgacj4l8iYnNEzCcbC5tY0v7GiFgVEW8C/4csUCrZ9xnAsxExL9/248CtwCdL2twW\nEQ8ARMRbVb6XQc9jRK0pgEnR/RhRtfYDft1Hmz2A/wQszYdrIPtwl/5Dty62HTf5A1nPaSDWlEzv\nCzxXtv43+fIuvyuZ3ljF/j8AHC9pQ8myocA/59MBPF/htqwbDiLryxqyw5neBrNfJvtgHx4Rv+3H\nPvo7yFz6c2uBs8rWfwC4swb7fw74j4gYX0VtVgUfmrWumpwBAr4P/I2kg5Q5UtLw0gZ5T+d6YI6k\nPQAkjZBU6Qf3RbLB6GqVvsc7gEMkTZU0VNK5wKHAj3toX+olsjGiA3tY/5N82xdK2j5/jZV0aB/b\ntQo5iFrX7flZqK7Xj0rWlfcAeusRfJtsPOVu4FWywNmxm5+7HFgFPCjpVWARcEiF+7gBODwfy7q1\nl3bltmwzH7c6A5hB1kP7C7JLGNZ31z6fjvxn/wD8L+B+SeslHV+2/jWyMajzyHpevyUbA3tP+bas\nfwp1QaOkCcAcsjMr34+IWT20Gws8AJwTEdX84ZpZARWmRyRpCHA1MAE4HJgq6bAe2s0C7sJdYrOW\nUJggAo4DVkXE6ojYBMwHJnXT7gvAD8mO682sBRQpiEaw7enY5/NlW+RX9E4CrskXFee40sz6rUhB\nVEmozAG+EtnAlvChmVlLKNJ1RGuBkSXzI3n3RWLHAPPzi+Z2B06TtCkiFpY2kuSeklkiEVF1B6FI\nPaJHgIMljZL0HuBcsvuVtoiIAyJidESMJhsn+mx5CJW0barXN77xjeQ1tHK9rrkxr/4qTI8oIt6W\n9Hmgg+z0/Q0RsVzS9Hz9tUkLNLO6KUwQAUR2M+SdZcu6DaCIuLghRZlZ3RXp0GxQa2trS11CVZqt\nXnDNRVaoK6trRVK04vsyKzpJRJMPVpvZIOUgMrPkHERmlpyDyMyScxCZWXIOIjNLzkFkZsk5iMws\nOQeRmSXnIDKz5BxEZpacg8jMknMQmVlyDiIzS85BZGbJOYjMLDkHkZkl5yAys+QcRGaWnIPIzJJz\nEJlZcg4iM0vOQWRmyTmIzCw5B5GZJecgMrPkHERmlpyDyMyScxCZWXIOIjNLzkFkZsk5iMwsOQeR\nmSXnIDKz5AoXRJImSFohaaWky7tZP0nSMkmPSVoq6eQUdZpZ7SgiUtewhaQhwNPAqcBa4GFgakQs\nL2mzU0S8kU8fASyIiIPKthNFel9mg4UkIkLV/lzRekTHAasiYnVEbALmA5NKG3SFUG5n4OUG1mdm\ndVC0IBoBrCmZfz5ftg1JZ0paDtwJfLFBtZlZnRQtiCo6noqI2yLiMOATwE31LcnM6m1o6gLKrAVG\nlsyPJOsVdSsiFksaKun9EbGudN3MmTO3TLe1tdHW1lbbSs2Mzs5OOjs7B7ydog1WDyUbrD4FeAFY\nwrsHqw8EnomIkPQR4N8i4sCy7Xiw2iyB/g5WF6pHFBFvS/o80AEMAW6IiOWSpufrrwWmAJ+StAl4\nHTgvWcFmVhOF6hHVintEZmm0yul7MxuEHERmlpyDyMyScxCZWXIOIjNLzkFkZsk5iMwsOQeRmSXn\nIDKz5BxEZpacg8jMknMQmVlyDiIzS85BZGbJOYjMLDkHkZkl5yAys+QcRGaWnIPIzJJzEJlZcg4i\nM0vOQWRmyTmIzCw5B5GZJecgMrPkHERmlpyDyMyScxCZWXIOIjNLzkFkZskNOIgkHSJph1oUY2aD\n09D+/JCkvwX2BB4GDgLeAr5Ww7rMbBDpVxABdwMrgV2AHwAfqVlFZjbo9PfQ7GVgbESsAD4HbK5d\nSWY22Cgiem8gHQicDewKrAEeiIjHG1Bbv0mKvt6XmdWeJCJC1f5cJT2ik4EfAfcApwDflLRU0gXV\n7qwSkiZIWiFppaTLu1l/gaRlkp6QdL+kI+tRR6N0dHQwfvwUxo+fQkdHR+pyzJKoJIi2A3aOiHuA\n2yPiDOCjwGZJn61lMZKGAFcDE4DDgamSDitr9gzwJxFxJPA3wHW1rKGROjo6mDx5GosWTWTRoolM\nnjzNYWSDUiVBdB3QJumnwJmSzgAOBJYAO9e4nuOAVRGxOiI2AfOBSaUNIuKBiHg1n30I2K/GNTTM\n7NnXsXHjLGAaMI2NG2cxe3bT5qpZv/V51iwfbJkj6Xtkh2knAp8A1gE317ieEWTjUF2eB47vpf2n\ngTtqXIOZNVjFp+/zHkpH/qqXikeYJX0c+HPgY/Urp75mzLiU++6bxsaN2fywYZczY8a8tEWZJdDf\n64jqZS0wsmR+JFmvaBv5APX1wISI2NDdhmbOnLlluq2tjba2tlrWWRPt7e0sWDBvy+HYjBnzaG9v\nT1yVWeU6Ozvp7Owc8Hb6PH3fSJKGAk+TnZ17gWwcampELC9psz9wL3BhRDzYw3Z8+t4sgXqevu/a\nweHdLGurdoe9iYi3gc+THf79ErglIpZLmi5pet7sr4H3AddIekzSklrWYGaNV3GPSNJTwE3A3wHD\ngFlkV1efUL/y+sc9IrM06t4jIjt7NRJ4gOyQ6bdk1xOZFcq6dXDOOfDqq323tWKoJojeBjaS9YZ2\nBJ6JCN9jZoWybh2ceiqMHg3vfW/qaqxS1QTREuBN4FjgJOB8Sf9Wl6rM+qErhMaPh6uuAlV9gGCp\nVDNGNDYiHi5b9mcRcVNdKhsAjxENPg6hYmjEGNFSSX8m6a/zHe4P/KraHZrVmkOo+VUTRP9AdnvH\n+fn868D3al6RWRUcQq2hmiurj4+IMZIeA4iI9ZK2r1NdZn1yCLWOanpEf8wf0wGApD3wkxktEYdQ\na6kmiP4eWADsmT88/37gyrpUZdYLh1Drqepes/whZScDAu4pvQesSHzWrHU5hIqtv2fNKnlm9YyS\n2SALoa5pIuLb1e603hxErckhVHz9DaJKBqt3IQudDwJjgYVkYXQG2UWOZnXnEGpt1VzQuBg4PSJe\ny+d3Ae6IiJPqWF+/uEfUWhxCzaMRFzTuCWwqmd+ULzOrG4fQ4FDNdUT/DCyRdCvZodmZgJ9ranXj\nEBo8qj1rdgzZDa8B/DwiHqtXYQPhQ7Pm5xBqTnU7a9aMHETNzSHUvOp51qxrBzsCU4BRJT8XEXFF\ntTs164lDaHCqZozo34FXgKVkzyUyqymH0OBVTRCNiAh/143VhUNocKvm9P0v8u8TM6sph5BVc0Hj\ncuAg4FngrXxxREThwsmD1c3DIdRa6j5YDZxW7cbNeuMQsi4+fW9JOIRaU91u8ZB0f/7f1yW9Vvb6\nfX+KtcHNIWTl3COyhnIItbZG3PRqNiAOIeuJg8gawiFkvakoiJQZWe9irDU5hKwv1fSI7qxbFday\nHEJWiYqCKB/5XSrpuDrXYy3EIWSVqubK6qfJrqz+DfBGvthXVlu3HEKDUyOurO664bXrE+4/LeuW\nQ8iqVe0TGo9m6xMaF0fEsnoVNhDuEaXjEBrc6n4dkaQvAT8A9gD2An4g6YvV7tBal0PI+quaMaIn\ngRMi4o18fifgwYg4oo719Yt7RI3nEDJo3JXVm3uYrhlJEyStkLRS0uXdrD9U0gOS3iz7FlpLxCFk\nA1XNYPWNwENlXyf0T7UsRtIQ4GrgVGAt8LCkhRGxvKTZOuAL+f4tMYeQ1ULFV1YDPwQuBjaQhcFF\nEfGdGtdzHLAqIlZHxCZgPjCptEFEvBQRj7Dtlz1aAg4hq5VqekR3RMSHyR6eXy8jgDUl888Dx9dx\nf9ZPDiGrpaJdWe0R5ibgELJaq6ZHdAJwoaR6Xlm9Fii9uXYkWa+oajNnztwy3dbWRltb20DqspxD\nyEp1dnbS2dk54O1UdPo+HyM6CXiufF1ErB5wFVv3MxR4GjgFeAFYAkwtG6zuajsTeC0iZnezzqfv\n68AhZH2p61dO50H0ZD5GVFeSTgPmAEOAGyLiSknTASLiWkl7Aw8D7yW7hOA14PCIeL1kGw6iGnMI\nWSXqGkT5DuYB34uIJdXupNEcRLXlELJKNSKIfPf9IOQQsmo08u77Uv60tzCHkDVKJV8ndBlsGZQe\nm19suDqfn17f8iwVh5A1UiXXEU0tmf6rsnX+9tcW5BCyRvO3eNg2HEKWgoPItnAIWSp9njWT9A7w\nh3x2GLCxZPWwiKhmwLshfNaseg4hq4W6nTWLiCH9K8mahUPIUvOh2SDnELIicBANYg4hKwoHUWId\nHR2MHz+F8eOn0NHR0bD9OoSsSKr6OqFm0SyD1R0dHUyePI2NG2cBMGzY5SxYMI/29u4uYq8dh5DV\nS6Menm81NHv2dXkITQP2ZuPG0Zx//ufq2jNyCFkROYgKoYMsjD7D+vVfZ/LkaXUJI4eQFZUPzRLa\nemg2GvgMsDdwHfACY8YM4dFH76vZvtatg1NOgfZ2h5DVjw/NmlB7ezsLFsxj+PCXgCfJekUTgc+w\nbNkva9YrcghZ0blHVAAdHR2cfvoFbN48m1r3ihxC1kjuETWx9vZ2jjrqw9S6V+QQsmbhHlFBbNsr\nmpYvnce4cQu5++4fVb29rhAaPx5mzXIIWWO4R9TktvaKunQA/8jSpcuq7hWV9oQcQtYM3CMqkK1n\n0S4E5gHfAqq70NGHY5ZS3R+e30yaNYggC6Pzz/8c69d/nWoHrh1ClpoPzVpEe3s7xxxzFNUOXHdd\nrOgQsmbkICqgGTMuZbvt5gJdt39MY/Pm7zB79nXdtvcV09bsHEQFVM3AtUPIWoHHiAqqkoFrh5AV\njQerS7RCEEH5wPW21xbdfPOPHEJWOA6iEq0SRADjx09h0aKJlJ5BO+KI9zFkyB0OISscB1GJVgqi\njo4OJk48jz/+cSjZ4dkOSB/i7LN35pZbRjuErFB8+r5Ftbe386EPHUUWQtOA84iADRv+wiFkLcNB\n1AR23/39wA753LPApTz6aPW3fpgVlYOoCVxyyeeQPgT8DDiOej/J0azRPEZUcF2n6A8++Fl++tNx\nbNjw7jNo/bk736wePEbUgkqvE7rlltEce+xR+ZoOYApwFQ8++HDDv4rIrNbcIyqo7i5W3PYM2jT6\ne4e+Wb24R9RCerpietszaM+y9Uxa9t1oPd2LVolUX/RoBgULIkkTJK2QtFLS5T20+W6+fpmkMY2u\nsd76um0jO4NWW123kyxaNJFFiyZ6ENwaLyIK8QKGAKuAUcD2wOPAYWVtTgfuyKePBx7sYVvRjF5+\nOeLooyMuuyxi8+bu29x1110xbNheATMCdg+YGzA3hg3bK+66665+7XfcuLPy7UT+mhvjxp01gHdi\ng1X+2av681+kHtFxwKqIWB0Rm4D5wKSyNhPJBkaIiIeA3STt1dgy66PSG1i7voJo3LhnGTPmg4wZ\ncyPjxi2s8fjQkyxdusyHadYwQ1MXUGIEsKZk/nmyXk9fbfYDXqxvafVV7V307e3tNR2UnjHjUu67\nbxobN0L2QLbrWb/+uyxaBPfdN82D4FZ3ReoRVXqaq/xj2tynx4CpU9PeRb+1l7WQ4cNvA75LrQbB\nzSpRpB7RWmBkyfxIsh5Pb232y5e9y8yZM7dMt7W10dbWVosa62LuXNhnn7R30Xf1srK7/dPVYc2l\ns7OTzs7OAW+nMNcRSRoKPA2cArwALAGmRsTykjanA5+PiNMlnQDMiYgTutlWFOV9NZutD2SbBfj6\nJKtOSzwGRNJpwByyM2g3RMSVkqYDRMS1eZurgQnAG8DFEfFoN9txEA1AR0fHlsOxGTMudQhZxVoi\niGrFQWSWhq+sNrOm5SAys+QcRGaWnIPIzJJzEJlZcg4iM0vOQWRmyTmIzCw5B5GZJecgMrPkHERm\nlpyDyMyScxCZWXIOIjNLzkFkZsk5iMwsOQeRmSXnIDKz5BxEZpacg8jMknMQmVlyDiIzS85BZGbJ\nOYjMLDkHkZkl5yAys+QcRGaWnIPIzJJzEJlZcg4iM0vOQWRmyTmIzCw5B5GZJecgMrPkHERmllxh\ngkjScEmLJP1K0t2Sduuh3T9JelHSk42u0czqozBBBHwFWBQRhwD35PPduRGY0LCqGqSzszN1CVVp\ntnrBNRdZkYJoIjAvn54HnNldo4hYDGxoVFGN0mx/cM1WL7jmIitSEO0VES/m0y8Ce6UsxswaZ2gj\ndyZpEbB3N6u+WjoTESEpGlOVmaWmiGJ83iWtANoi4neS9gF+FhGH9tB2FHB7RBzRw/pivCmzQSgi\nVO3PNLRH1IeFwDRgVv7f2/q7of78IswsnSKNEV0FjJP0K+DkfB5J+0r6SVcjSTcDvwAOkbRG0sVJ\nqjWzminMoZmZDV5F6hH1W7NcDClpgqQVklZKuryHNt/N1y+TNKbRNXZTT681SzpU0gOS3pQ0I0WN\n5Sqo+YL89/uEpPslHZmizpJ6+qp3Ul7vY5KWSjo5RZ1lNfX5t5y3GyvpbUln9brBiGj6F/B3wGX5\n9OXAVT20OwkYAzyZoMYhwCpgFLA98DhwWFmb04E78unjgQcT/14rqXkP4Fjgm8CMAvwtVFLzicCu\n+fSElL/nCuvdqWT6CGBV0X/HJe3uBX4MTOltmy3RI6I5LoY8juwPaHVEbALmA5PK2mx5HxHxELCb\npJTXU/VZc0S8FBGPAJtSFNiNSmp+ICJezWcfAvZrcI2lKqn3jZLZnYGXG1hfdyr5Wwb4AvBD4KW+\nNtgqQdQMF0OOANaUzD+fL+urTcoPSSU1F021NX8auKOuFfWuonolnSlpOXAn8MUG1daTPmuWNIIs\nnK7JF/U6GF2k0/e9aoGLISutqfzSg5TvpYi/x75UXLOkjwN/DnysfuX0qaJ6I+I24DZJJwE3AR+s\na1V9lFNBmznAV/LPo3j33/U2miaIImJcT+vyAei9Y+vFkP+vgaVVai0wsmR+JNm/JL212S9flkol\nNRdNRTXnA9TXAxMiIuW9i1X9jiNisaShkt4fEevqXl33Kqn5GGB+lkHsDpwmaVNELOxug61yaNZ1\nMSQM8GLIOnoEOFjSKEnvAc4lq7vUQuBTAJJOAF4pOeRMoZKauxTlItI+a5a0P3ArcGFErEpQY6lK\n6j0w71Ug6SMACUMIKqg5Ig6IiNERMZpsnOizPYVQ1w80/QsYDvwU+BVwN7Bbvnxf4Ccl7W4GXgDe\nIjvGvbjBdZ4GPE12xuEv82XTgeklba7O1y8DPlKA322vNZMdLq8BXiU7EfAcsHPBa/4+sA54LH8t\nKXi9lwFP5bUuBsYW/e+irO2NwFm9bc8XNJpZcq1yaGZmTcxBZGbJOYjMLDkHkZkl5yAys+QcRGaW\nnIPIzJJzEJlZcg4i65akvSXNl7RK0iOSfiLp4Cq3saukz/ay/v4qtzezVg9fU2a6pEskHViLbVr/\nOYjsXfL7mhYA90bEQRFxLPCXVP94lfcB/62nlRFR7V3vtbwN4EtkzyL6GXB2Dbdr/eAgsu58HPhj\nRFzXtSAinoiI+yR9WdKT+etLXesl7ZT3mh7P150DXAkcmD/idFb5TiS9nv93lKTlkq6T9JSkDkk7\n5uu+KulpSYspefSFpAslPZRv+x8lbZcvH5s/VnWHvKanJB1ett/tgTMi4nHgA8CuNfzdWT80zWNA\nrKE+DCwtXyjpGOAisif0bQc8JOk/8g/0BGBtRPyXvO17yXocH46Inp69XdrDOQg4NyIulXQLMEXZ\nd92dCxxF9kjSR4FHJB0GnAN8NCLekfQPwAXATRHxsKSFZI+uHZYv+2XZfk8GXpM0DfgEcE81vxyr\nPQeRdaenQ6D/DNwaERsBJN1K9hzwx4EngG9Jugr4cd57Gl7FPp+NiCfy6aVkz0PePd/fm8CbecCI\nLEiOIQslyALndyXbuoLsURUbyR5XWu5E4IaI+LGkTwIPVFGn1YEPzaw7/5fsg14u2Pa5Q8qXEREr\nyb+YAPimpK9T3ZjOWyXT77D1H8ny/XX9d15EjMlfh0bEFSXtdgd2Inu+87Bu9rUP8IykHYB9IuLx\n/Jsy9q2iXqshB5G9S0TcC+wg6ZKuZfkTDR8HzpQ0TNJOZF9SsDhfvw/wZkT8C/AtslB6DdhlAKX8\nPN/fjpJ2Ac4gC7d7gLMl7ZHve3j+sLMu1wJfA/6V7JuDy60jC76zgG/nX1AwjeI83G3Q8aGZ9WQy\nMCf/zqo3gWeB/wHMBZbkba6PiGX59BHA/5a0mewbPT4TEeuVfW/Yk2Rfk1T+/VfRwzRkjx9/LB8v\nWkb2+N8l+Yrlkr4G3J0PUm8iOzv3nKRPAW9FxPx83S8ktUVEZ8m2byYLodcj4hoAScuwZPxgNDNA\n0jeA70dEymeED1o+NLNBT9KeZJcGfDx1LYOVe0Rmlpx7RGaWnIPIzJJzEJlZcg4iM0vOQWRmyTmI\nzCw5B5GZJecgMrPkHERmltz/ByIYHWdS+ZIYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f94535b8710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_line,y_line = 0.01*np.array(range(0,20,1)),0.01*np.array(range(0,20,1))\n",
    "plt.scatter(P2U[:,5],P2U[:,7])\n",
    "plt.plot(x_line, y_line)\n",
    "plt.axis([-0.1,0.4,-0.1,0.4])\n",
    "plt.xlabel('Cost index $\\\\theta_1$')\n",
    "plt.ylabel('Error index $\\\\theta_3$')\n",
    "plt.title(r'Efficient frontier')\n",
    "plt.axes().set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((708.44105390000016, 3332.536183964849),\n",
       " (476.56300040000002, 2740.2201476429686),\n",
       " (943.54962525000008, 4940.2050685173981))"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test2U = policy(100,5,0.05,0.001,0.00027,0.3,best2U[0],best2U[1],best2U[2])\n",
    "a = np.mean(test.cost_calc(f_test,best[1])),np.std(test.cost_calc(f_test,best[1]))\n",
    "b = np.mean(test2.cost_calc(f_test,best2[1])),np.std(test2.cost_calc(f_test,best2[1]))\n",
    "c = np.mean(test2U.cost_calc(f_test,best2U[1])),np.std(test2U.cost_calc(f_test,best2U[1]))\n",
    "a,b,c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Cost risk-under estimation does not produce better results\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}