180202_MSM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# INPUT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## input from xls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5P/DPMzPJZGdJAoQ1Yd8RiOyyKwoqVqvV1l2L\nW7/Vn7YWbW1dUKlVWy21St23Wq1aUUQFBAVkERCQPSxhzwIBsk9mOb8/7jL3zp7Mncz2vF8vXszc\nuXPn3mTyzJlznvMcEkKAMcZY4jNF+wQYY4y1Dg74jDGWJDjgM8ZYkuCAzxhjSYIDPmOMJQkO+Iwx\nliQ44DPGWJLggM8YY0mCAz5jjCUJS7RPQCsvL08UFhZG+zQYYyyubNq06aQQIj/YfjEV8AsLC7Fx\n48ZonwZjjMUVIjoUyn7cpcMYY0mCAz5jjCUJDviMMZYkOOAzxliS4IDPGGNJggM+Y4wlCQ74jDGW\nJDjgs5jRaHfiv5uOgpfdZCwyYmriFUtu85fsxuvflSI/24pJfYNOGmSMNZMhLXwiaktE/yWi3US0\ni4jGElF7IlpKRCXy/+2MeC2WuI6daQAAlJ9tjPKZMJaYjOrSeQ7AF0KI/gCGAdgFYC6A5UKIPgCW\ny/cZ88vhdAEA7v9wG7YfOxvls2Es8YQd8ImoDYCJAF4BACFEkxDiDIDZAN6Qd3sDwGXhvhZLbA6X\nu+9+5Z6KKJ4JY4nJiBZ+EYBKAK8R0Q9E9DIRZQLoKIQ4Ie9TBqCjAa/FEliTw6Xe3nTodBTPhLHE\nZETAtwAYAeCfQojhAOrg0X0jpLQLn6kXRDSHiDYS0cbKykoDTofFK7vTHfBX7OH3AmNGMyLgHwVw\nVAixXr7/X0gfAOVEVAAA8v8+v6MLIRYKIYqFEMX5+ZyZkcy0XTqMMeOFHfCFEGUAjhBRP3nTNAA7\nASwCcIO87QYAn4T7Wiyx2Z0c8BmLJKPy8P8PwDtElArgAICbIH2YvE9EtwA4BOAqg16LJSiHpkuH\nMWY8QwK+EGILgGIfD00z4vgsOXCXDmORxaUVWMzQZukwxozHAZ/FDIdLH/Cd3OJnzFAc8FnMcHgM\n2tq5T58xQ3HAZzGj0e7U3W/igM+YoTjgs5hh8+jDP3yqPkpnwlhi4oDPYoZnls7Ff18dpTNhLDFx\nwGcxIdiiJxf+7Vu88V1p65wMYwmKAz6LCW+vPxzw8d1lNfjToh2tdDaMJSYO+CwmbDhY5fexUGfg\nVtU1YV9FTdjnUl7diPc3Hgn7OIzFGg74LOZ0yknT3W/UDOYGWhhlxGNLMf3Zb8N+/Rtf+x73/3cb\nztQ3hX0sxmIJB3wWc4Z0bQMAmD6gAwB9umadzeHzOSdrbYa9/q4T1QCA+iZnkD0Ziy8c8FnMMRHQ\npW06ctJTAOgDPhH5fI5nDn9L3fbWRvW2EvgZSxQc8FnMGdq1LUwmQEnceUczoOtZfgGQavAY0Rrf\nU1aDL3eUq/e/2F4W9jEZiyUc8FlMcGoC+R2TesFMBJcc8XdrWtq+6uuMfGwpLvhr+H33dU367qIP\nNh0N+5iMxRIO+Cwm1DRKwXZYt7YwmQgmIiix/ZJhndX93ttwBOXVjfrn+unXby7PWj7ZVqOWi2As\nNnDAZzGhW/sMAMD0/tJALRHgkiO+NkVy8Y8nMPqJ5RE5B8/yzDU2B4QQOFlrw6Of7sTZentEXpex\n1sIBn8WEnnmZAIDrxxUCAMwmd5fOugPeOfrKIK2RFTWbnNIxf3lekbrt0Kl6LNtZjlfXHMSCFSWG\nvRZj0cABn8UEpW/eYpKycEyaPnxfbHYp0H+y5bhh5/D17goAwOUjurpfx+FCrdxlxCtysXjHAZ/F\nBKcc3M1ywCciBGq82xxSazwj1WzYOby9TsoGSrWYMLhLDgDgo81HMW/xLgCA1WLcazEWDRzwWUxQ\n+utNpLTwAxdUU0opp6Xo38LKN4RwpKeY8dsZ/QEAL317QN1uRNkGxqKJAz6LCQ6PLh2lD99f0Fda\n+E0O/eMOl//nBJNill67Y06azw+OZbsqWnRcxmIFB3wWE5ocLqSYCSZNl45LuD8I7ju/r27/RrkP\nX5mI9emvJqiPFT3weYvOYUiXNjivTx7MJsK4Xrk+93FxPz6LYxzwWUywOVxINbvfjiYCXEKoqZKp\nFpPH/lILX8mdz7SG379udwqkyOfgr4RDzwdb9mHCWCzggM9igs3hhDXFHbSVmbZKX73VM+DLLXxl\n3dsUc/hvZbvTpXbrePpZcbewj89YtHHAZzGhyauFT3C5oGnh61vwjR4tfM+A/+SSXT7LMPhjd7qw\nu6wGK3ZXqtu0HzINmuJsRhVqY6y1ccBnMcHmcMGqybghklI1/XbpePThW8yE2yb1VB9/6ZsDWLv/\nVMiv/728AMukfvnqtq9/M1m9raRpAsDe8ho0OVxYvstdaI2xeGBYwCciMxH9QESfyffbE9FSIiqR\n/29n1GuxxNPkcOla1CYinK5rwsS/rADgo0tH/iDYUyalSqaYTbhnmn5gV5k5GwrlG8NdU3qr27Td\nO788ryfunNwLAPDR5mPo+4cluOWNjSgp51RNFj+MbOHfDWCX5v5cAMuFEH0ALJfvM+aTzeHSteLN\nJkJJRa16P9Viwvu3jcUvRneX95cCtFI6Wcrw0R/z5tc3Blw6UUuZ5GXWDNZaTNpvHIQJffIAAK9r\nFlM/dKqeM3dY3DAk4BNRVwCzALys2TwbwBvy7TcAXGbEa7HEZHM4dTNZPZNkUi0mjCpqj3vl9Ewl\nLVNhMZmQajapM3UVV720FjuPB1/I5Nu9Ut+99kMjJ82Cfh2z8cIvRgDwPdP21jc34p/f7A96fMZi\ngVEt/L8BuB+A9q+woxDihHy7DEBHX08kojlEtJGINlZWVvrahSWBNftOYW+Zu3vE5BHxlS4d5VuA\nZ2XLFDOBiPDQrAFex575/KqgA61vrTsEQN+qt5hN+PL/TcTMIQW6c/C0qoTftyw+hB3wiehiABVC\niE3+9hHS1Eef33uFEAuFEMVCiOL8/Hxfu7Akoa1r7znRVQm2SgveJQScLgGLiXDn5F5q3ry2G0ir\nujH80saeZRwUfTpkh31sxlqDESs8jAdwKRHNBJAGIIeI3gZQTkQFQogTRFQAgOelM5989YGv2KNv\nNaeape4UpeXvFALLd5XD4RLo0i5d3c/fYuZ1NifgJy5rSzEESudXzsFTVV2T/ycxFkPCbuELIR4Q\nQnQVQhQCuBrA10KIawEsAnCDvNsNAD4J97VYYlIqZd40vtDvPkrKphLwXS6BOW9JXyrzsqzuY/kZ\nQK0LsCqWkvHTrX06euVn+d1PeHxJvXBQJwDSoiyMxYNI5uHPB3A+EZUAmC7fZ8yLEqTzs61+91Em\nZSldOtrSydqAf/5An0NF+PMXu/0eWwn4N44r8ltSAQAK2qRjWNc26v0Xrxvpd1/GYpGhi3YKIVYC\nWCnfPgVgmpHHZ4lJCfjmAME2xaK08KX7LiEwqqg9Nhyswjnd2qr7XVXcDRcOKsCwR7/SPX/zodN+\nj62kePoblFWkWkz45FcTUFljU9fVvXNyL10JZcZiGc+0ZVG3eJvUJbLzhP/0SeWjgIik9W6FgM3u\nRFFepi4Vk4jQJiNFvb933kXItlrUlv+Jsw3494bDun57ZdZusICvyM+2YnCXNvJzzHC6RLPKODAW\nLYa28Blrifs/3AYACBQzM63ut6qZCPsqarH16Nmgx061mJCfbYVTAC+s3IenvtgDAFi5pwJzLxqA\nBz7ahrE9pQlV2uJtodKmiaYbuPoWY5HAAZ/FjHmXDVZv3ziuEK9/V4rz+uRh/hVD0Sbd3Wp3uASW\nbC8LeKyBBTnqNwazieB0udRgDwBf7ijHjuPVOHq6QV0kPTut+X8OyrcCm8PJAZ/FPA74LKq0Bc60\nQV0JvsU92qNL23Sv5wXz0Z3j1MlWZhOpVTW1PFM4R/ZofrknpYX/9rpDGN69Hb4vrcI90/sGeRZj\n0cEBn0XVp9uO+9x+2fAuWPjtAVw2vHPA5/92Rj+f29NSzEiTu2hSzCbYfayI7lmeISctxWufYJQW\n/tNf7VW33T2tT8BsH8aihQdtWVQdPlXvc3uv/CzsmXcReuRmBnz+8O5tAz4OAFlWC2oD5OGHY0BB\njte2e9/fGpHXYixcHPBZVNWEGYgzU4N/Sc1Jt+D7Uv9pmQDw/1rYDdOvk/f03Y9/ONaiYzEWaRzw\nWVT1lme2rv7dlBY9PyOEgdJQumrunt6nRa8faO4AY7GGAz6LqvomB/p2zELXdhkten6GNXgLP1sT\n8HvmeXcRvXJDcYteGwBMnlXeZPsrfRdxYyyaOOCzqGqwO9XB1ZbICOG5mVb3Pp5LJQLAlH4dWvz6\nnpQxhenPfmPYMRkzCgd81ipuef17/OVL73o2jeEGfGvw52rz43eX6ZckfHT2IL+t9Jb44fAZAIDg\nibcsBnHAZ61i+e4K/GOF98pQDXZXswO+sgIV4C6qFoi/bwF5Wam4fmxhs147kA9uH4uCNmmGHY8x\no3EePouKs/V2fPbjcdQ02NG1XfMmVikrUAEIKd89w08mz8Y/nN+s1w1kWv8OOLewPdLlD5cJvfMM\nOzZjRuGAz1rd44t34l+rDqr3p/Y3rg/dF21xtTdvHoWTtTZ8vdu49Xi2PzJDnYClfP54LsHIWCzg\nLh3Wqs7UN+mCPQB0yPFfB9+fKf3y1dZ0MErAv3RYZ0zsm4/LR3TFgp+PCPKs0GVZLUiRu5Z+dm43\nAFLuP2Oxht+VrFU98fkur23aGjqheu2mUSHvq9TlaUlxtOb65Xk98fzyfWifmRrx12KsuTjgs4jT\n1p4/UtXg9XiWtfkBvzmm9OuAB2f2x5Uju0X0dQC5Hn96CnyU7mEs6jjgs4jTLg6y9sApr8ezItzy\nNpkIcyb2iuhraCnlmBmLNdyHzyLO7qM0sVZWCLNl44nZRAhyyYxFBQd8FnE1Nrvu/rzLBuOnI7uq\n93NaoW+9NZlNBBcvechiEAd8FnGPfLpTdz8vK1XNZgGAthmJNcBpJoKDu3RYDOKAzyLuZI1+ZalM\nq0XXjdOSLJ1YZjIRD9qymMQBn0XcsTP6zJy0FDP6dnTXkfdV0CyeWUwEFxfTYTEosf7SWMxxOF04\nelof8O0Ol272a6IxmQgO7sNnMSixRstYzKlrcuru52dbcW5RewDAd3On4ky93dfT4pqZwIO2LCaF\n3cInom5EtIKIdhLRDiK6W97enoiWElGJ/H+78E+XxZtXV+vLKMy/fIhahqBz23QM7Oy9Jmy8s5hM\nurkHjMUKI7p0HADuE0IMBDAGwF1ENBDAXADLhRB9ACyX77Mk89zyEgDAz4q7Ictqwcgeif+5bzJJ\nk83WHzilm2XMWLSFHfCFECeEEJvl2zUAdgHoAmA2gDfk3d4AcFm4r8XiV8/8TGx/ZEbCpWD6YjYR\nNpRW4WcL16Hogc+jfTqMqQwdtCWiQgDDAawH0FEIcUJ+qAxARyNfi8WXZBrENJs4F4LFJsPemUSU\nBeBDAPcIIaq1jwnpe63Pv3gimkNEG4loY2VlpVGnw2LMaHmgNhmYEzcBicU5QwI+EaVACvbvCCE+\nkjeXE1GB/HgBAJ8rTgghFgohioUQxfn5+UacDosRdTYH8rOtuHRYZxQXJlHA90g55QFcFivCTssk\naY25VwDsEkI8q3loEYAbAMyX//8k3Ndi8WHFngrY7C7c/vYmANLM2mTiGfDPNti5Pj6LCUb8JY4H\ncB2AH4loi7ztQUiB/n0iugXAIQBXGfBaLA7c9Nr3uvsrDFxOMB54BvyquiYO+CwmhB3whRCrAfjr\ntZwW7vFZbJry9Ep0bZeOJy8fArtToCgv0+++xYWJn4qp5TloW1XXFKUzYUwvub5rM8McPFmHgyfr\nMOHPKwAApfNn+d33/hn9W+u0YoIyaJtltaDW5uCAz2IG54+xkDmcLoyf/zX+98Mxr8eaHP7LQ+Zn\nN3+R8nhmkrt0am0OANzCZ7GDAz4L2felp3HsTAPu+c8Wr8ea5HrAn249rtu+4cFpSE81t8r5xQqL\nRx9+fZMjSmfCmB4HfBayU3U2v4855TX91h/Ur1mbkWQZOoD3oO28xbtQOHcxqhu9C8Wt3FOBwrmL\nseXImdY6PZbEOOCzkHkuZKKlrPDUMy9Ltz09Jbla94B7Dd+bxxfptt8hp6lqvbxKKi532T/WRP7E\nWNLjgM9C1hign/7tdYcBSP3VZhPhu7lT8fGd4xK67r0/DXapJHRhXoZ+u0epaCDxFn9hsY3fbSwo\nIQTeXX844MDsX5ftBQDsOlGN3vlZ6Nw2HcO7J1c6psImB/zMVH13lq9KoYer6lvlnBgDOOCzECzf\nVYEHP/4Rzy7dG3Tfk7U2dGqT1gpnFbuUFr7nDON/rTroVWaBF0phrYkDPguqqj60tMIjVfU402BH\n24zEWpS8uZSumywfA9b7Kmp193Oz3DNwD52qi+yJsaTHAZ8FlZfluyzAAxf1x4d3jFXv3/z69zh0\nqh7ZacmXmaPVaJe6vjKt3gPWJ2v1A982TTfZqRbm67+17hDWHTgVfEeW9JL7L5OFxFfffVqKCbdN\n6gUAyLZaUGNzoERuvdY0JnfeeaPdu4WfnWZBTaMDNod74LaiphHbjp5V79sDjJEE8tD/tgMIPNuZ\nMYBb+CwEt7+92Wvbi9eOVG+/fvMo3WOzz+kc8XOKZUofvnYOwuM/GQJA/+H50Wb9jOW5H/3YCmfH\nkhkHfNYi2tZrr3x94bSp/ZN7cTOlha+dg9A7X5qfoF35Ky9LKjkxoEBayP3gyeb34TucLftWwJIT\nB3zWItrkkjRNYCvMzfCxd3J58vIh6N0hC23S3YPXSr69NkvHLgfr+87vq3u+0yWw64Ru0Ti/6nzk\n9jNJdaMdZxu8ZzcnM+7DZy1i0azjpw34pac4r/zCwQW4cHABAOC7uVNhtZhQLY9ruIR3wNembzpd\nAs8vL8Fzy0vw0Z3jMCLIXIarF64z+vQTxtCHvwIAPHvVMByuqsc90/sGeUbi4xY+a5ZOOWl4+sph\nGN6trc/Hrz63WyufUWzr3DYduVlWmEn6gHS6gONnGnD5C2tQdrYRgD6bx+50Ye1+KeNGeTyQUL8J\nJLN739+Kvy0rifZpxAQO+Cwgz4lCD108ED8d2RVE+pIJSiomt6J8U9ZEcbkE/rXqADYfPoMXVu4H\nAGRoZuSuLjmJDaVVAKQ1gQPZfPh0s86hvLpR/VaR6Hz97G5/axPqmxxJPdmNAz4LaNmucgCA1WLC\nwSdnYtbQAt87yn9DyVYKOVRKTSGnEF718bUD4Le+uVG9XR+kf/7V1Qd197Upn57O1tsx+onleOqL\n3SGfczzbU17jte2LHWUY+McvceFz37b6+QghYmKAnQM+C2hvmfSH85sL+nm16nXkhwLtkszcXToC\n5dX6rpq0FP2foTLAWxekjr7waKi+tfYQnl2612dgGfao1J+9Yk9ls847Xh073eD3sb3ltdh1ohrC\n8wcYQQu+3ofev1+iZnBFCwd8FpCSU/6LMd0D7vfenDG4bWJPZCdh/ftQKKtguYTwqiDq+a1IydX/\neLP3ymIKh9OFxT+e0G175qu9eH55CdZ6zLrVBpnKACWuE0lDkMB60XOr8PXuilY6G+AZuQ5VtH/+\nHPBZQBtLTyMz1azrZ/ZlUOc2eGDmgMDfApKYtoU/rleeuv3aMd1htfjuBiupqPXbDbD9uHuwdvGv\nJwBwB7kUs/7Puv9DX6i3kyVNcc2+k0H32V3m3e0TaQdaMNfCSBzwWUDHzzagbYbvWjosdEoL3+kS\navlkADAF+YC0+Sm3oA3cuZn6NYM/3nwMj3220+/g5KznV6Fw7mKvQm6J5JMt+qU23711tNc+FdXB\ns6CMdsOrG1r9NbU44DO/hBCoqLbh4mF+BmpZyJReHCH0g7El5VLQ3fzQ+br9x/RsDwB++3yVLJRH\nZw9CB49F4v+z8QheWX0QV720FgCQmWrGJcPc5S52yN8OVjSzS8PpEq3a720kX2d9PIS0VyNEu99e\niwM+8+tUXROanC50yknu+vZGUPrtl+4s1/Ud18vBoH1mKjI0fflZVmmWrr8Wfq0c8Kf066B+e/C0\n8ZCUtplqMaFtunfJas/KncH0evBz3PaW9zKN8aBfp2xcM0o/R2R9K1UY9bXSWbRwwGd+vbn2EACg\nV35WkD1ZMErXzYbSKl0/rjZUF+a6axIpwf+bvb6zamrlmbtKSqe/lSRfWX0QDqfQzYxWHGrGrGil\nZf/VzvKQnxNN6SlmzBpSgANPzMSBJ2YiL8uKyf06AAAuGNgRI7q3VctYv77mYESD/+fb9YPrH246\nGrHXCoYDPvPr+eXS7MTxvfOC7MmC8be2r1JADQDe0FQdzZe7aR7wU0FTGaBVMnz8zSWqszlgd7mQ\nYjbhP3PG6B77YkcZPtniPxNISztmUBtkQlhrq6prwrajZ9T7b60tRYPdifaZqTCZSP0GNGNQJ2z9\n0wVYeH0xxvbKhVNIXVQPf7oTN7/+fcTOz7Nc+H0fbMXfl5dEpXss4gGfiC4koj1EtI+I5kb69Zhx\nesiF0JJxIXKjWXz8DB+7bDCevnKoej9f0xffLsiqYTa7E0TShDgAGNxFqrh5zSgpffYnw7sAkFr+\njXYXLCbC6J65uHBQJ91xnvpiT0jnf/yMu79bKf0QK37zwVZcumANfpBnHj/0yQ4AQKmPFcSUgnYm\nIriEUIvPGVGEbtnOcp8zmZVf/bWa1OZnlu7FmfrWz5iKaMAnIjOAfwC4CMBAANcQ0cBIviZrnsK5\ni1E4d7HXVPRTtTZevcpARISJffPV+/ee3xfXjenhNwNq1lD3IKvn76a60Y4tR88izWJW02DfuGkU\nHrl0EJ74yWCUzp+Fv/7sHGRZLXjxmwMAgP2V0uDwn68YikdnD1KPFWqphRNn3ROZPGf4RpsS2H/y\nwne67Z4T3LRMRBCi5Wmqn249jsK5i7FXntG77sAp3PrmRjz9lf4DdMHXJXji893qa2odOV2Pxz7b\nicK5i3HXu95rTkRCpFv4owDsE0IcEEI0AXgPwOwIvyYLkfYrpecSeSPnLQPAq1cZSbtU5MV+SlQ8\ndPFAnNcnDz3au8tMV9bYdIFp6MNf4du9leipWYcgN8uKG8YV6uZB5Galqt0vyuBvm4wUXD+2EB/e\nMQ4AUFFjC6m2jHaAtyTG0jm1pSm+3u0eY/CsA6WlfGs9E+J6zZ7+798/AAAu+KtUpkH5/eyv0H+r\nePqrvertu6b01j126YI1eEX+8Fy8Td/PHymRDvhdABzR3D8qb1MR0Rwi2khEGysrk2Pad7R9uaMM\nFdWNeHvdIXWb0l8PhN7qY82TYnL/uXlOjlLcMqEIb90yGiYT4fczBwAAJj+9EsMekUojaFfMygwy\nGa59pvsD5ukrh+keG9nDXXZ5x/HgFTebnO7g2dzsHqN9s7dSN/bQTfPhOO+zXertQJ9jSsB/eVXL\nvq0oXWkKdykH3y86pEsbdAyS7dYa6ZtRH7QVQiwUQhQLIYrz8/ODP4GFxeZw4ra3NuHqhevUvk4A\n2KpZW1X7R/DbGf1a9fwSWa2mNo5nwPClXaZ3d8/Ff1+l3laqavqjTMjqlJOmGxxWjO+dCwBqt0Qg\nLV1v12hCCNzw6gbc/d4W3TZFYZ77W0+gFcSU7pWPf5A+OFJD+H1o9e2Yrbt/uKrea/uBSvc3od9d\n2B8AMKqwvd9j+svIMlKkA/4xANrk167yNhYlSheN5xTvy4e7v3i9/p0U8B+c2d/rayhruSNV7jTI\nUAJMikcqpd3pwt5ydxC59/zApai7tksP+PijswcDAJZsLwMgdYH4+3bncMVGwP+fpmX/xfYyVDfa\ndRPZQv12mu5RsK7J4WpW1oy2Nb7j+Fm8/l0pAKglrwFg6jPfqLdHyxPpBhToPygA4O/XDAfQOt2n\nkQ743wPoQ0RFRJQK4GoAiyL8miwAh9P3mzpLHpx1uQTKq6Wv7FePClwwjTWPdrWrUMpIe9YlKinX\n9537LVUtUwZ7/Q1MFsl5/0oQuv3tTejz+yU+97XL75tJffPRv5N30Got2m6sV9ccxNCHv8LKPZVq\nRpm2Sqa/cRJA3/EyqLOU4dQU4MNix/GzWCIXq9tdVq0bx5j1/Grdvp5jIuN756pdeA9fOgh/vmII\n3pszBhmpZnx5z0Q17bm2MfJZOxEN+EIIB4BfAfgSwC4A7wshdgR+FgtVVV0T5i/ZHXBwypO/lppy\njOOabAyufGksZaIPAKT66cPX8kyhnPn8Kt197UQtX26eUATA3YL0ZDIRstMsastyaYBJVUrL+Ux9\nE3aX1UQlh3xVSaUu82bDQXeX1qFT9ZgxqKP6zfXl64vx7FXn+D2Wdm7JtAEdAfif1QxIQf2Odzaj\npLwGF/5N+j3k+Mlg6/ng5/hOU7xtct8O6m0iws/O7Y4xPXOx89EL0a9TtrriWWvMb4h4H74Q4nMh\nRF8hRC8hxOORfr1k8sdPtuPFb/Y3qybKCT/1Q95Zfxi1NoeaY/3WLaO48qXBtN0AofxsUy0mvx+6\ny+6dFHR+xICCHJTOn4XpAzv63aem0YHXvytFvWZ8wdfg4bdy/7Iy1nOkyn+9+UhwuQSue2WDbtzJ\nU+8O7hnhBW3TAnab9e2YjQ/vGIcXrx2pzn+w2YN3B137ynr19nMeH6TaD9YFK/YBAGYO6YTrx/UI\neEyrxYxUswmrQ6jwGa6oD9qyllNqdDSnrXXnO/p83ws0wWDwn77E5sPSjMU+HaL3tZ25aVMvFfMu\nG6wLbkbYWOpeLnHRVqnS5L83HMaox5dBCKG+Ly6Vi7Ad9DGpKZIcIXyL1S4GP7AgJ+j+I3u0w4WD\nO6kD6P5WDOv7B3c3l9LdObJHO4wpylW3j+jeFhM03xq+kxtOz1893G/5a61+nbLRITvyNas44Mcx\npU84hN4BVaPHjML/m9pHd//fGw4D0C+szYyhNOqbk/nkmQ0CAD8d2dWoU0LfjtIHx/Wasr1KdsvD\ni3agosamBjkAKC6U0jn9lSJYsbsCP3tprdcyjuEKNmj80MUD1Xz83MzUZn07VQK+r282p2ptulRY\nAOiYY8XQDZnQAAAbDUlEQVSHd4xDWopJnUU7ons7tMtMxY3jCnX7WkL84/z4znF47LLBIZ9zS3HA\nj2Nl8h9ic97cl2myca4f2wNDurbxuV96Cgd8o5FcKm3WkNDLTT90yUD8eqo7U6p3hyykGfi7uX1S\nL939VLPJa+lEpauhb8cstTSBdtzo0Kk6NQNp/pLdWH+wSp3ZaxTPFn4fj284t0woUgd0m5vPrvw8\n7/tgm9dj2pISCmXshIjUXH/l28Wdk90/z+kD/HelebKYTerPNpI44MexXSekCTPmZgR8pzZTxE/g\n+MOsASG3TFjoWjIkkpOWgrs0Ad/oRUs8M31SLSbYnS5sO3pGHcRcVSL13+8tr8UlmpIPh+Vqm5P+\nshLnPbUCf126V1083OgBSG12Wec2aWinKUnx7W+nAHAH3ebWxVH68LceOeP1WGWtO+Ar6w4M69bW\na7+xvaTuHe3ciWF+GlPRxH/VCSDYqkla764/rN5WBrW2PXyBbp/cLF7hKhKUdMZQUjK1tBk9Rtc2\nslrMmKfpSki1mNDkcGH5LncigNK6zpWrTyqu+dc6XWmC5zSztd9a657FbQTtUo8dctJw7wXSHIQ2\n6SnoLqdkes5bCJW/8t/XvbIeN7++EQCw6v4pui4jxZyJPQG4B4y1M6i3+PgAiTbOu4tT2q/UrhBT\n5DxT6ZQ3cE5aCpbdOxHTn5XqgnRvHzjdj7XMs1edg61HzwSdYu9J22V3boCZmi01Tm6dts1IQarZ\nhLMNdrylKbuh1Hl57aZzAUgDlpsOncaxMw0459GlPo/ZIzcD972/Ffsqa/HJXePDPkdtl860/h0w\npmcuSufP0u0zqLPUoh5d1LyfkbYrxeZwqoOsq0rcWTP52Va0kSuYarN/fjujH34+qrtuJnO/jtnY\nU16jfhjEEg74capB008Zah7+sTPuVLrbJvbE9WML1fu9NVk5I7p7f2Vl4cu0WnQLmLfE835y6sNh\nkWv85GVZYXM41SwdT0rm1tNXDsOUp1cGPOZra0qNPEVdl46/2d+d2qR5fQiESgnSDU1SwPf8m0pL\nMeMHOVNJ262WYjbpyjkAwN+uPgd7ymowumcuYg136cQpbcncUKaTO11CXdh5Yt98PDBzgN+uBc6/\njz3L7p2IrX+8QFcZ0ihd26Xjl+cV4V/XFwecEKa8X4ryMvHgzP4hHz9YNc6tR84ELcimzIK9fmwP\nv0s6huPG8YUA3DVxfJVWvmKElB3lWa7a04CCHF1yRCzhFn6cqtCkylWHUIOj14Ofq7fvnsb1ceJN\n7wjOizCZCL+fJS1TkRpCzjiAkHLLFSfrbAFzzGf/Yw065aRh3YPT/O7z5Q6p3o9niqRRlCUlL12w\nBqXzZ+nqHimevnIozuneFpdqBq7jDbfw41SZpgWy+0Q1hBB49NOd2H7sbIBnSfwN8r77y9F4+fpi\nw86RxR/P+vDKQOjci/Qtem0/9pyJPfHET4ao91+5Qf8emvmcvtaMljKuVBZgsRIA6rjHTeOLAu7X\nUp6zlrdqlkz8jTxATES4bkwPtS8/HnELP05p/0Dq7U5U1tjw6pqDeHNtKfY9MTPgc3v5maUZbv8y\ni3+epTdmDOqEz7adwJUek7205Z0n9c3H2J65sJgJo4vao4dHjZ+TtTY4nC6fqb52Td98ZY1Nt8yj\nVmWN9I02WAXQlvKcmausUrVn3oXN+jYT67iFH6fKzzbCRNJA27vrD2OnnJOvzWb4bNtxnzW2c9Li\nt4XCWtcTlw/BqvunINejnr42eKenmmEyEa4q7qYG+wdn9sfDl7hXM/1w81Gfx9cmH3whd9v4Umuz\nw2IiXfkEI/XUpGZqa1MlUrAHOODHrbLqRuRnW9XBrhtfc091V6bG/+rdH3DDqxt4BSsWsmn93ZUd\nS+fPQk5aim5FKcWUfu7FijJ8DP7PmdgLN44vwvzLpa6e3334IwCptPC9/9mi9sVru5Ae+t92v+fV\naHeFtGiMEbQrwSUaDvhxqry6EZ385HM7PAL8nxZJFQYL2qThwzvGRvzcWPwKVoFTka35lhioFazU\n3lH8eclufPTDMazZdxLfl1bhNx9sVR9TVuDypdHuNLSkhC87H50BAFgut/CvPrdboN3jEgf8OHXw\nZB06+An4dqfAnz5xt5aU2bXv/nIMRvYwfuIOSxzKJD7PWvy+KP3emQFmDvfukK0GciGE2iVztsGO\nK19ci+/lKp09cjMCZuA02l0RD/iepUa081YSBQf8OHS6rglHTzf4XbBi9b5KvOFjansel0xgQfx8\ntLTK2cOXDgq67ye/Go9l90702/BQTOwjdf/UNznVOR7LPdZwyM1MDVhhU5oBG9lwRUS6wnbFCdg4\n4oAfh7QtD1+ldpUMA09cAZMFM7V/R5TOn4VObYKXf0gxm0KaH6B0/9Q0OtRywp96zOYd0qUNDp2q\n99vKb7S7YG2F96/2W8TNEwoj/nqtjQN+HNKWf70gwGpGnrgCJosGZb3kWpsd/kYIBndpA4dL4MRZ\n390otTY70lIi//7dUHpKva1dPzdRcASIQ0oO/uM/GazLY2YsFikVPktP1uN/W7zr9PTvlK3uU2fz\nLm1cdrYR6w5UeW2PBG39+0iUcIg2DvhxqKS8FiaSanu0zwzcL7/g51KxrUDrezIWScpi3yv2+F57\nudHuRLrcmm6we5cJeWX1AQCtU254ptyHH4kidbGAo0AcKqmoQff2GUhLMev6Wp+6Yqhuv4/vHIfh\n3aW0uLwgHwyMRUqWVerDVypeeqYTT+qbr+by13ssXlLdaFdXdhvUOfg6teF6+sqhWPvAVHXt3kTD\nAT8Onaxt0tVUV7IXrvLIG+7bMRvt5Lofd0zWL2XHWGtRumvONEhZOP+8doT62JOXD8HciwaoCQWe\nAX/ow1+pA7x/veqciJ+r1WJGQZvIlG+IBYk3KpEEHE4XMjQDShsenA67vMjzO7eOxi9eXg9AmgFJ\nRC2uEc6YEZQFRr7cIaUR56Sn4KXrRuL55SW4fEQXWC1mtfRyQ4DlCXv6WZmKhY4DfhxyuAQsmuXc\ntNX7Bnd2r6PJde1ZLPCsf9OlbTp65WdhhmZyl9Klo83F91yhLdRZwMw/7tKJI59sOYbCuYux7ehZ\ndZUiT9ZWSF1jLBy+ZsxmpEgfCo9+thM/Hj2LjaVV6iLqAPxW0WTNE1Z0IKK/ENFuItpGRB8TUVvN\nYw8Q0T4i2kNEM8I/VXb3e1vU2/4Ce1qKGQMLcnBJgg46sfgUbF1b7eLslyxYjZ++uFbXn98tQmWR\nk024XTpLATwghHAQ0Z8BPADgd0Q0EMDVAAYB6AxgGRH1FUL476BjzdIhQIvn87vPa8UzYSy4Yd3a\nouTxi/w+bjIRnrpiKO7/cJu6bcRj7gXSrxvbI6LnlyzCauELIb4SQiiJs+sAKKskzAbwnhDCJoQ4\nCGAfgFHhvBbT23GsOtqnwFizpJhNSAkw27tnfqbP7V3apqv58Sw8Rnb43gxgiXy7C4AjmseOyttY\nC3kOYE3un+9nT8bi00A/efbv3Do64RYiiZagAZ+IlhHRdh//Zmv2+T0AB4B3mnsCRDSHiDYS0cbK\nSu/VmZjE5lFUiv8AWKLJSLXg11N7e233tQALa5mgffhCiOmBHieiGwFcDGCacDdDjwHQzgLqKm/z\ndfyFABYCQHFxMReG8UMpmJaXlYqTtU1ei5wwlgiuHdsDz3+9T73/7i9HczqmgcLN0rkQwP0ALhVC\n1GseWgTgaiKyElERgD4ANoTzWsnsQGUtth49CwC4a0pv3DapJw9isYTUITsND108ENlpFmRZLRjX\nKy/ap5RQws3SWQDACmCpPMlnnRDidiHEDiJ6H8BOSF09d8VKhk6tzYFMeQZqvJj6zDfq7faZqbhp\nfFEUz4axyLplQhFumcDv8UgIK+ALIbw73NyPPQ7g8XCOb7QDlbWY+sw3eOqnQ3FVcXyuV5mfxRNQ\nGGMtk1TTMpUlAbcfOxvlMwnd7jJ9+mWw5eQYY8yfpAn4Qgg8uURa+q81q+E98fkuvLW2tMXPP3Sq\nXnefp5gzxloqaYqnrdzrTvlMMbde//3Cb6XFG64bW9ii59/21ibdfaXyIGOMNVfStPD/tqxEve10\ntU7259HT9UH3aWhy4kBlbdD9zCbClH482Yox1nJJ08IvP+teq9IRJOAfO9OA8fO/xkd3jsMIecWo\nFr1mdWPQfQb88QsAwMrfTEZhnvfU8m7t01Hcoz2euXJYi8+DMcaAJGrhn6y1qbeDtfAXyQstf7Dx\nSMD9gvFXwtgXf+t1NtpdSEsxwWSihFxUmTHWepIm4HfRlFcN1sI/XFUHAPj3hvACvr0Zs2G1+zba\nnXhh5T7U2hxotDu5jAJjzBBJ06XTKz8LdTYnTtc3wekKHIjtTmP6+M822NXbr64+iKq6JvxmRj91\n22nN6j4mzUSwjzYfw1Nf7EGq2QSb3eVzwQjGGGuupGnhH6mqR1FeBpwugX+s2I/qRrvffbOsxnwO\nagP+o5/txIIV+3SP/0Nz/7NtUjeSyyXw4Mc/AgDmLd6FJqdL/cbBGGPhSIoW/v7KWpRU6DNhrnjh\nOwDA0nsn6bZX1TXh9e9KAQAdc8LLef9g49GAj2urO6w9cAoAMPP5VV77NTm4phxjLHxJ0cI/UuWd\nHllS4f0hAACvrD6g3naE2bWjBHEtbV+90nV0bmE7NNpdeOmb/Th6usHrOU/9dGhY58EYY0CCB3wh\nBArnLsZHm6XKzC/8YoTP/T7/8QQ+3CS1xrXjrKfqmmBztLzmW7sM70lS+zU59412JzpkW/F96WkA\nwJNLdqPW5vB6TvvM1BafA2OMKRI64K/ZJ7WwF22V+sd7d8jyud+d72zGfR9sBQB8X1qle+y9Fmbq\nbDpUhdP13uMEjyzaCUD6MHrv+yOoqLF57aN1z/Q+LXp9xhjzlNAB/8RZffdIWx9lCbRdLEIIZKTq\nM2KOn/XuYgnFFf9cC8C79o3SzVOqqZFz47hCr+f36ZCFt28ZjV9P5YDPGDNGQgf8Ro9lAdtmeHeN\nbDvqnvB0pt6u9tsrNWtSfSy6XGtzwBVieYbMVN8plSXlNertP8waoHvs4UsGYum9kzChTx5PtmKM\nGSZhA35DkxPPLdur25ZqMXkVH9tb7u5TL69phNMlMLZnLjb9QVrZceUe/Tq7/910FIP/9CWeXao/\ntj+lp7wHjF9ZfRBz5KJov585ABazCT88dD66tJUmhxW0bb1qnoyx5JGwAf+RT3fgZG2T1/b/3TVe\nd/90vXufu97ZjA2lVXAJAYvcstc+/sX2MvxG7uv/fPuJZp/T2J65SDWb8NhnO9Vt14+Tlipsl5mK\nIrmWzpii3GYfmzHGgknYPPwfPRY5+a08w7XIo0DZYU0LfH+lNMFp/UFp4PbyEV2w/oB7EPcP//tR\nvZ0SoE6Ov+qX+dlWmEwANIk/2rIJz119DvZX1qGNj+wexhgLV8IG/EqP7Jebxhf63M9f0TJAWlC5\nosZd8bJJMyYQqG/9iCaXfvqADgAIy3aVIy/Lika7/7IOuVlW5PIShoyxCEnYgO9ZETPNTwGy3WU1\nPrcDQHqKGXangNMlYDYRqhvdOfI+xnJVds0Hw8s3nIs6mwP7Kmq9Uj4Hdc4JdAmMMWaohO3Dt8lB\nd0BBDnY8MqNZ2S6d20jrxiopm00O71a5mXwf74fDp9EkP2/eZYMBAJlWC4Z1a4t0j4ydiX15QRPG\nWOtJ2BZ+TpoFtTYHLhlWgEyPYmj/ur4YnXLScMmC1T6f++39UwBALXb2bUklZgzqhCyrRZ0J6+sD\n5MVv9mP+kt1QHhrTUz/4euKMu3uoS9t0nlTFGGtVCdnCr29y4PjZRvTvlI3bJvbyevz8gR0xpGsb\n3bY7Jrv3UzJ0bhgrZdDkpEmDqEK4u4k8W/hNDhfmy4ukK71JnmvnZqW5P3junNKL69wzxlpVQgZ8\npUrl7rIamAN05VwwsKN6e3yvPADA7ZPcgf+iIQUA3IG+rsmdXmPyCPjPLN3jdfwUj47+PprSDp3b\ncK49Y6x1JWSXTju52NiAgsCDokrQnjW0AON756J0/izd40qZhfom7wJqnlmZ+yu8a9ZbPFr40wZ0\nxMd3jsOZejsm84LkjLFWlpABX8mS+fs1wwPup3Sx3DiuEORjEFYN+HbvgJ/q0R3jazlDXytVDQ9j\nUXTGGAuHIV06RHQfEQkiytNse4CI9hHRHiKaYcTrhKqsWhoczUkP/Hn2x0sG4qkrhqK4h+8gnJEq\nPb+hybtk8bd79SUXfAV8pe+fMcZiQdgtfCLqBuACAIc12wYCuBrAIACdASwjor5CiJYXl2+Gv3wp\n9aenB1kLNictBVed283v44G6dDx5TvRSUjsZYyxWGNHC/yuA+wFoZzrNBvCeEMImhDgIYB+AUQa8\nVrMEC/hBny8H/Ec+3Ympz6wMuG+hR8kGrnLJGIs1YQV8IpoN4JgQYqvHQ10AaFcOOSpv83WMOUS0\nkYg2VlZW+tolJDWNdmyVyyR0a5+Oqf07qOmVLaUtjXygsg7ZaRZM7d/B574ul8DAghx8dOc4AL67\neBhjLJqCRkQiWkZE2338mw3gQQB/DOcEhBALhRDFQoji/PyWZ67MeXMTZv9jDSqqG1Hd4EDXduGn\nPXoO5NY0OtAxJw33nt8XgBTkG+1OPLlkF06cbUR2mgU58kCwr9m5jDEWTUH78IUQ031tJ6IhAIoA\nbJUDY1cAm4loFIBjALSd413lbRGjrCQ16onlAHwvdmIEE0HN7Xe4BBZtOY6XvpEWPp8+oCOy5YFa\ne5gLoDPGmNFa3OchhPhRCNFBCFEohCiE1G0zQghRBmARgKuJyEpERQD6ANhgyBmHKNsamYxTIsAi\nB3ynS8Ca4v4R5qRbkCW/LnfpMMZiTUSiohBiBxG9D2AnAAeAu1orQ0cRqYlNJiJNC9+FGk0FzYxU\nMzJSzbhxXCEuPadzRF6fMcZayrCAL7fytfcfB/C4UcdvLl8TqQw5LtxdOtWNDlQ32tXHrBYziAgP\nXzooIq/NGGPhSIhaOp617wHvla2MQkTYdOg0AOCJz3fpBmdTLQnx42SMJaiEiFAbDuoXFnnmymEB\ni6Y1h1UO4rPlLppv91aqK1rVNDp0HzapYaaBMsZYJCVEhBqpKY3wx4sH4oqRXQ079pf3TMTbt4xW\nC60dOFmHc+XX69cxSx/wuYXPGIthCVE8LdViwraHL0BGijnsyVaeCvMyUZiXiTfXlgIABnfJwZxJ\nPfHy6oPo3j4DRzXr11o54DPGYlhCBHwg8oXKlHo6913QDylybeRjZxpxpt49aGvhcgqMsRiWMAE/\n0mwOKeCnp5hhluvcv/jNft0+KdzCZ4zFMI5QIVLSPIXw35LnQVvGWCzjCBUiJca7hPBa3hAABnXO\nwfmaJRMZYyzWcJdOiJQ0T5cQPlv4i399XmufEmOMNQu38EM0uEsbAEBuptWwHH/GGGtN3MIP0W8v\n6IcLB3XCwM6BF0ZnjLFYxS38EFnMJl6AnDEW1zjgG+DWCUXRPgXGGAuKA74BIlSYkzHGDMUBv4Ve\nu/FcXDOqOwCpLDJjjMU6HrRtoSn9O2Bsr1xkp1lwx+Re0T4dxhgLigN+GNJSzHhw5oBonwZjjIWE\nu3QYYyxJcMBnjLEkwQGfMcaSBAd8xhhLEhzwGWMsSXDAZ4yxJMEBnzHGkgQHfMYYSxIkhIj2OaiI\nqBLAIYMOlwfgpEHHihfJds3Jdr1A8l1zsl0v0LJr7iGEyA+2U0wFfCMR0UYhRHG0z6M1Jds1J9v1\nAsl3zcl2vUBkr5m7dBhjLElwwGeMsSSRyAF/YbRPIAqS7ZqT7XqB5LvmZLteIILXnLB9+IwxxvQS\nuYXPGGNMI24CPhF1I6IVRLSTiHYQ0d3y9vZEtJSISuT/22me8wAR7SOiPUQ0Q7P9GiL6kYi2EdEX\nRJQXjWsKprnXTES58v61RLTA41gj5WveR0TPE8XewoxGXS8RZRDRYiLaLR9nfrSuKRgjf8eaYy4i\nou2teR2hMvg9nUpEC4lor/y7viIa1xSMwdccXuwSQsTFPwAFAEbIt7MB7AUwEMBTAObK2+cC+LN8\neyCArQCsAIoA7AdghrToSwWAPHm/pwA8HO3rM+iaMwFMAHA7gAUex9oAYAwAArAEwEXRvr5IXS+A\nDABT5NupAFbF4vUa/TuWH78cwLsAtkf72iJ9vQAeATBPvm1S/qZj7Z+B7+uwY1fctPCFECeEEJvl\n2zUAdgHoAmA2gDfk3d4AcJl8ezaA94QQNiHEQQD7AIyCFPAIQKbcys0BcLzVLqQZmnvNQog6IcRq\nAI3a4xBRAYAcIcQ6Ib1T3oT75xQzjLpeIUS9EGKFfLsJwGYAXVvlIprJqGsGACLKAnAvgHmtcOot\nYuT1ArgZwJPyfi4hRExO0DLwmsOOXXET8LWIqBDAcADrAXQUQpyQHyoD0FG+3QXAEc3TjgLoIoSw\nA7gDwI+QflgDAbwS+bMOT4jX7E8XSNevOCpvi1lhXq/2OG0BXAJgucGnaDgDrvkxAM8AqI/E+Rkt\nnOuVf68A8BgRbSaiD4go5PdFtIRzzUbErrgL+HIr5kMA9wghqrWPya3XgGlHRJQC6Yc2HEBnANsA\nPBCZszVGuNccb4y6XiKyAPg3gOeFEAcMP1EDGfC+PgdALyHEx5E7S+MY8Du2QPrW9p0QYgSAtQCe\njsS5GiUWYldcBXz5gj8E8I4Q4iN5c7ncZaF0XVTI248B6KZ5eld52zkAIITYL/+Q3wcwrhVOv0Wa\nec3+HIO+S0P5WcQcg65XsRBAiRDib8afqXEMuuaxAIqJqBTAagB9iWhlZM44PAZd7ylI32SU538A\nYEQETtcQBl1z2LErbgK+3Gf1CoBdQohnNQ8tAnCDfPsGAJ9otl9NRFYiKgLQB9LA5TEAA4lIKTR0\nPqQ+tZjTgmv2Sf7aWE1EY+RjXh/sOdFg1PXKx5oHoA2Ae4w+TyMZ+Dv+pxCisxCiENKA314hxGTj\nzzg8Bl6vAPApgMnypmkAdhp6sgYx8H0dfuwyaiQ60v8gvYkFpK8xW+R/MwHkQuqfLQGwDEB7zXN+\nDyk7Zw80WRqQRr93ycf6FEButK/PwGsuBVAFoBZSX/1AeXsxgO3yz2MB5El3sfTPqOuF9A1GyL9j\n5Ti3Rvv6Iv071jxeiNjN0jHyPd0DwLfysZYD6B7t62uFaw4rdvFMW8YYSxJx06XDGGMsPBzwGWMs\nSXDAZ4yxJMEBnzHGkgQHfMYYSxIc8BljLElwwGeMsSTBAZ8xxpLE/wfNftmTvIVE2QAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1f12574e5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from pandas import DataFrame as df\n",
    "\n",
    "GLD = pd.read_excel('data_GVZ_GLD.xlsx')\n",
    "date_GLD = GLD.iloc[:,3]\n",
    "GLD = GLD.loc[:,'GLD']\n",
    "GLD_d = GLD - np.mean(GLD)\n",
    "kbar = 5\n",
    "plt.plot(date_GLD,GLD_d)\n",
    "data = GLD_d\n",
    "\n",
    "input_param = [50, 1.99,0.945585657754, 5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TEST My own function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## - T_mat_temp \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def T_mat_temp(kbar):\n",
    "    A = np.fromfunction(lambda i, j: i ^ j, (2**kbar, 2**kbar), dtype=int)\n",
    "    A_template = A.astype(float)\n",
    "    return A_template"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MSM_starting_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "expected an indented block (<ipython-input-4-c8ea5f34c5dc>, line 29)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-4-c8ea5f34c5dc>\"\u001b[1;36m, line \u001b[1;32m29\u001b[0m\n\u001b[1;33m    function [startingvals,LLs,output_parameters] =...\u001b[0m\n\u001b[1;37m                                                      ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m expected an indented block\n"
     ]
    }
   ],
   "source": [
    "def MSM_starting_values(data,startingvals,kbar,A_template):\n",
    "    # options=optimset('fminbnd');\n",
    "    xtol = 1e-3\n",
    "    LLs = []\n",
    "    output_parameters=[]\n",
    "    \n",
    "    # implicit  booleanness of the empty list is quite pythonic.\n",
    "    if not startingvals:\n",
    "        print('No starting values entered: Using grid-search')\n",
    "        # A grid search used to find the best set of starting values in the\n",
    "        # event none are supplied by the user\n",
    "        b = [1.5, 3, 6, 20]\n",
    "        lb = len(b)\n",
    "        g = [.1, .5, .9]\n",
    "        lg = len(g)\n",
    "        sigma = np.std(data)*np.sqrt(252)\n",
    "        output_parameters=np.zeros((lb*lg,3))\n",
    "        LLs=zeros((lb*lg,1))\n",
    "        m0_lower = 1.2\n",
    "        m0_upper = 1.8\n",
    "        \n",
    "        index = 1\n",
    "        \n",
    "        for i in range(0,lb):\n",
    "            for j in range(0,lg):\n",
    "                \n",
    "    \n",
    "## ติดค้างก่อย\n",
    "function [startingvals,LLs,output_parameters] =...\n",
    "    \n",
    "                temp_in = [b(i);1.5;g(j);sigma];\n",
    "                [m0,LL]=fminbnd('MSM_likelihood',m0_lower,m0_upper,options,kbar,data,A_template,[b(i);g(j);sigma]);\n",
    "                parameters=[b(i);g(j);m0];\n",
    "                output_parameters(index,:)=parameters';\n",
    "                LLs(index)=LL;\n",
    "                index=index+1;\n",
    "            end\n",
    "    end\n",
    "    [LLs,index]=sort(LLs);\n",
    "    startingvals=[output_parameters(index(1),:),sigma]';\n",
    "    output_parameters=output_parameters(index,:);\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MSM_likelihood \n",
    "\n",
    "output form matlab using demean GLD\n",
    "\n",
    "[parameters,LL,exitflag,output]=fmincon('MSM_likelihood',startingvals,...\n",
    "    [],[],[],[],LB,UB,[],options,kbar,data,A_template);\n",
    "   \n",
    "   param = [50,1.99,0.945585657754,5]\n",
    "   LL = 56886.6637963524\n",
    "   exitflag = 4\n",
    "   output = \n",
    "   \n",
    "  [LL,LLs] = MSM_likelihood(parameters,kbar,data,A_template);\n",
    "  LL = 56886.6637963524\n",
    "  LLs = \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MSM_likelihood(input_param, kbar, data, A_template, estim_flag):\n",
    "    if len(input_param) == 1:\n",
    "        input_param = [estim_flag[0], input_param, estim_flag[1], estim_flag[2]]\n",
    "\n",
    "    sigma = input_param[3] / np.sqrt(252)\n",
    "    k2 = 2 ** kbar\n",
    "\n",
    "    A = transition_mat(A_template, input_param, kbar)\n",
    "    g_m = gofm(input_param, kbar)\n",
    "    T = len(data)\n",
    "    pi_mat = np.zeros((T + 1, k2))\n",
    "    LLs = zeros(1, T)\n",
    "    pi_mat[0, :] = (1 / k2) * np.ones((1, k2))\n",
    "\n",
    "    # Likelihood Algorithm\n",
    "\n",
    "    pa = (2 * np.pi) ** (-0.5)\n",
    "    s = np.matlib.repmat(sigma * g_m, T, 1)\n",
    "    w_t = np.matlib.repmat(data, 1, k2)\n",
    "    w_t = np.divide(pa * np.e ** (-0.5 * np.power(np.divide(w_t, s), 2)), s)\n",
    "    w_t = w_t + 10 ** -16\n",
    "\n",
    "    for t in range(0, T):\n",
    "        piA = np.matmul(pi_mat[t, :], A)\n",
    "        C = np.multiply(w_t[t, :], piA)\n",
    "        ft = sum(C)\n",
    "\n",
    "        # stop div by zero if prob are too low\n",
    "        if ft == 0:\n",
    "            pi_mat[t + 1, 0] = 1\n",
    "            print(\"prob=0\")\n",
    "        else:\n",
    "            pi_mat[t + 1, :] = np.divide(C, ft)\n",
    "\n",
    "        LLs[t] = np.log(np.dot(w_t[t, :], piA))\n",
    "\n",
    "    LL = sum[LLs]\n",
    "\n",
    "    if np.isinf(LL):\n",
    "        print('Log-likelihood is inf. Probably due to all zeros in pi_mat.')\n",
    "\n",
    "    return LL, LLs\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - transition_mat\n",
    "\n",
    "\n",
    "Calculate transition matric using template\n",
    "\n",
    "\n",
    "$$ \\gamma_{k}=1-(1-\\gamma_{1})^{(b^{k-1})} $$\n",
    "\n",
    "<font color = red>equation (3.2) [Page 21 in MSM book]</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " for i in range(0,kbar):\n",
    "        print(i)\n",
    "kbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def transition_mat(A,input_param,kbar):\n",
    "    b = input_param[0]\n",
    "    gamma_kbar = input_param[2]\n",
    "    print(gamma_kbar)\n",
    "\n",
    "    gamma = np.zeros((kbar,1))\n",
    "    gamma[0]=1-(1-gamma_kbar)**(1/(b**(kbar-1))) #TOWARN : Equation does not conform bbut math is ok \n",
    "\n",
    "    def bitget(number,position):\n",
    "    \n",
    "        bi_number = bin(number)\n",
    "        bi_number = bi_number[2:]\n",
    "        if len(bi_number) >= position :\n",
    "        \n",
    "            fn_output = bi_number[-position]\n",
    "        else :\n",
    "            fn_output = 0\n",
    "        \n",
    "        return fn_output\n",
    "\n",
    "    for i in range(1,kbar):\n",
    "        gamma[i,0]= 1-(1-gamma[0])**(b**i)\n",
    "        \n",
    "\n",
    "    gamma = gamma/2\n",
    "    gamma = np.append(gamma,gamma,axis=1)\n",
    "    gamma[:,0]=1-gamma[:,0]\n",
    "    print(gamma)\n",
    "    kbar1 = kbar+1\n",
    "    kbar2 = 2**kbar\n",
    "    A2 = A\n",
    "\n",
    "    prob = np.ones((kbar2,1))\n",
    "\n",
    "    for i in range(0,kbar2):\n",
    "        for m in range(1,kbar+1):\n",
    "            prob[i,0]=prob[i,0]*gamma[kbar1-m-1, int(bitget(i,m))]\n",
    "    \n",
    "    for i in range(0,2**(kbar-1)):\n",
    "        for j in range(i,2**(kbar-1)):\n",
    "            \n",
    "            A2[kbar2-i-1,j] = prob[kbar2-int(A2[i,j])-1,0]\n",
    "            A2[kbar2-j-1,i] = A2[kbar2-i-1,j]\n",
    "            A2[j,kbar2-i-1] = A2[kbar2-i-1,j]\n",
    "            A2[i,kbar2-j-1] = A2[kbar2-i-1,j]\n",
    "            \n",
    "            A2[i,j] = prob[int(A[i,j]),0]\n",
    "            A2[j,i] = A2[i,j]\n",
    "            A2[kbar2-i-1,kbar2-j-1] = A2[i,j]\n",
    "            A2[kbar2-j-1,kbar2-i-1] = A2[i,j]\n",
    "    return A2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - gofm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gofm(input_param, kbar):\n",
    "    m0 = input_param[1]\n",
    "    m1 = 2-m0\n",
    "    kbar2 = 2**kbar\n",
    "    g_m1 = list(range(0,kbar2))\n",
    "    for i in range(0,kbar2): \n",
    "        g = 1\n",
    "        for j in range(0,kbar): #not req -1\n",
    "            if g_m1[i]&2**j != 0:\n",
    "                g = g*m1\n",
    "            else :\n",
    "                g=g*m0\n",
    "        g_m1[i]=g\n",
    "    \n",
    "    g_m = np.sqrt(g_m1)\n",
    "    return g_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0,\n",
       " 2,\n",
       " 4,\n",
       " 6,\n",
       " 8,\n",
       " 10,\n",
       " 12,\n",
       " 14,\n",
       " 16,\n",
       " 18,\n",
       " 20,\n",
       " 22,\n",
       " 24,\n",
       " 26,\n",
       " 28,\n",
       " 30,\n",
       " 32,\n",
       " 34,\n",
       " 36,\n",
       " 38,\n",
       " 40,\n",
       " 42,\n",
       " 44,\n",
       " 46,\n",
       " 48,\n",
       " 50,\n",
       " 52,\n",
       " 54,\n",
       " 56,\n",
       " 58,\n",
       " 60,\n",
       " 62,\n",
       " 64,\n",
       " 66,\n",
       " 68,\n",
       " 70,\n",
       " 72,\n",
       " 74,\n",
       " 76,\n",
       " 78,\n",
       " 80,\n",
       " 82,\n",
       " 84,\n",
       " 86,\n",
       " 88,\n",
       " 90,\n",
       " 92,\n",
       " 94,\n",
       " 96,\n",
       " 98,\n",
       " 100,\n",
       " 102,\n",
       " 104,\n",
       " 106,\n",
       " 108,\n",
       " 110,\n",
       " 112,\n",
       " 114,\n",
       " 116,\n",
       " 118,\n",
       " 120,\n",
       " 122,\n",
       " 124,\n",
       " 126,\n",
       " 128,\n",
       " 130,\n",
       " 132,\n",
       " 134,\n",
       " 136,\n",
       " 138,\n",
       " 140,\n",
       " 142,\n",
       " 144,\n",
       " 146,\n",
       " 148,\n",
       " 150,\n",
       " 152,\n",
       " 154,\n",
       " 156,\n",
       " 158,\n",
       " 160,\n",
       " 162,\n",
       " 164,\n",
       " 166,\n",
       " 168,\n",
       " 170,\n",
       " 172,\n",
       " 174,\n",
       " 176,\n",
       " 178,\n",
       " 180,\n",
       " 182,\n",
       " 184,\n",
       " 186,\n",
       " 188,\n",
       " 190,\n",
       " 192,\n",
       " 194,\n",
       " 196,\n",
       " 198,\n",
       " 200,\n",
       " 202,\n",
       " 204,\n",
       " 206,\n",
       " 208,\n",
       " 210,\n",
       " 212,\n",
       " 214,\n",
       " 216,\n",
       " 218,\n",
       " 220,\n",
       " 222,\n",
       " 224,\n",
       " 226,\n",
       " 228,\n",
       " 230,\n",
       " 232,\n",
       " 234,\n",
       " 236,\n",
       " 238,\n",
       " 240,\n",
       " 242,\n",
       " 244,\n",
       " 246,\n",
       " 248,\n",
       " 250,\n",
       " 252,\n",
       " 254,\n",
       " 256,\n",
       " 258,\n",
       " 260,\n",
       " 262,\n",
       " 264,\n",
       " 266,\n",
       " 268,\n",
       " 270,\n",
       " 272,\n",
       " 274,\n",
       " 276,\n",
       " 278,\n",
       " 280,\n",
       " 282,\n",
       " 284,\n",
       " 286,\n",
       " 288,\n",
       " 290,\n",
       " 292,\n",
       " 294,\n",
       " 296,\n",
       " 298,\n",
       " 300,\n",
       " 302,\n",
       " 304,\n",
       " 306,\n",
       " 308,\n",
       " 310,\n",
       " 312,\n",
       " 314,\n",
       " 316,\n",
       " 318,\n",
       " 320,\n",
       " 322,\n",
       " 324,\n",
       " 326,\n",
       " 328,\n",
       " 330,\n",
       " 332,\n",
       " 334,\n",
       " 336,\n",
       " 338,\n",
       " 340,\n",
       " 342,\n",
       " 344,\n",
       " 346,\n",
       " 348,\n",
       " 350,\n",
       " 352,\n",
       " 354,\n",
       " 356,\n",
       " 358,\n",
       " 360,\n",
       " 362,\n",
       " 364,\n",
       " 366,\n",
       " 368,\n",
       " 370,\n",
       " 372,\n",
       " 374,\n",
       " 376,\n",
       " 378,\n",
       " 380,\n",
       " 382,\n",
       " 384,\n",
       " 386,\n",
       " 388,\n",
       " 390,\n",
       " 392,\n",
       " 394,\n",
       " 396,\n",
       " 398,\n",
       " 400,\n",
       " 402,\n",
       " 404,\n",
       " 406,\n",
       " 408,\n",
       " 410,\n",
       " 412,\n",
       " 414,\n",
       " 416,\n",
       " 418,\n",
       " 420,\n",
       " 422,\n",
       " 424,\n",
       " 426,\n",
       " 428,\n",
       " 430,\n",
       " 432,\n",
       " 434,\n",
       " 436,\n",
       " 438,\n",
       " 440,\n",
       " 442,\n",
       " 444,\n",
       " 446,\n",
       " 448,\n",
       " 450,\n",
       " 452,\n",
       " 454,\n",
       " 456,\n",
       " 458,\n",
       " 460,\n",
       " 462,\n",
       " 464,\n",
       " 466,\n",
       " 468,\n",
       " 470,\n",
       " 472,\n",
       " 474,\n",
       " 476,\n",
       " 478,\n",
       " 480,\n",
       " 482,\n",
       " 484,\n",
       " 486,\n",
       " 488,\n",
       " 490,\n",
       " 492,\n",
       " 494,\n",
       " 496,\n",
       " 498,\n",
       " 500,\n",
       " 502,\n",
       " 504,\n",
       " 506,\n",
       " 508]"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# f=lambda x: x*2\n",
    "# L=list(map(f,g_m1))\n",
    "\n",
    "# # TODO: use Lambda function Reduce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TEST full"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'A' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-abbfcdb3177c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mestim_flag\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mA_template\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mT_mat_temp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkbar\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mLL\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mLLs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mMSM_likelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0minput_param\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkbar\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mA_template\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mestim_flag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-5-9ed3d395b195>\u001b[0m in \u001b[0;36mMSM_likelihood\u001b[1;34m(input_param, kbar, data, A_template, estim_flag)\u001b[0m\n\u001b[0;32m      6\u001b[0m     \u001b[0mk2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m2\u001b[0m \u001b[1;33m**\u001b[0m \u001b[0mkbar\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m     \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtransition_mat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA_template\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minput_param\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkbar\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      9\u001b[0m     \u001b[0mg_m\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgofm\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0minput_param\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkbar\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m     \u001b[0mT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-6-1f4177f32498>\u001b[0m in \u001b[0;36mtransition_mat\u001b[1;34m(A_template, input_param, kbar)\u001b[0m\n\u001b[0;32m     37\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkbar\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     38\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 39\u001b[1;33m             \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mprob\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     40\u001b[0m             \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     41\u001b[0m             \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkbar2\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'A' is not defined"
     ]
    }
   ],
   "source": [
    "estim_flag = []\n",
    "A_template = T_mat_temp(kbar)\n",
    "LL, LLs = MSM_likelihood(input_param, kbar, data, A_template, estim_flag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TEST part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OPTIMIZOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### fminbnd \n",
    "\n",
    "* [scipy ref](https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.optimize.fminbound.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### fmincon \n",
    "This has various issues due to the typr of problems as [This link](https://stackoverflow.com/questions/49926/open-source-alternative-to-matlabs-fmincon-function).\n",
    "\n",
    "- [stack overflow scipy minimize](https://stackoverflow.com/questions/27851087/matrix-constraints-with-scipy-optimize-minimize)\n",
    "- [scipy ref](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df_A = df(A)\n",
    "df_A.to_excel('transition_mat2.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### - Test gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gamma_kbar = input_param[2]\n",
    "\n",
    "gamma = np.zeros((kbar,1))\n",
    "gamma[0]=1-(1-gamma_kbar)**(1/(b**(kbar-1))) #TOWARN : Equation does not conform bbut math is ok \n",
    "\n",
    "for i in range(1,kbar):\n",
    "    gamma[i,0]= 1-(1-gamma[0])**(b**i)\n",
    "\n",
    "gamma = gamma/2\n",
    "gamma = np.append(gamma,gamma,axis=1)\n",
    "gamma[:,0]=1-gamma[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.955052</td>\n",
       "      <td>0.044948</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.934117</td>\n",
       "      <td>0.065883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.904506</td>\n",
       "      <td>0.095494</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.863833</td>\n",
       "      <td>0.136167</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.810361</td>\n",
       "      <td>0.189639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.744521</td>\n",
       "      <td>0.255479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.670998</td>\n",
       "      <td>0.329002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.600000</td>\n",
       "      <td>0.400000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0         1\n",
       "0  0.955052  0.044948\n",
       "1  0.934117  0.065883\n",
       "2  0.904506  0.095494\n",
       "3  0.863833  0.136167\n",
       "4  0.810361  0.189639\n",
       "5  0.744521  0.255479\n",
       "6  0.670998  0.329002\n",
       "7  0.600000  0.400000"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df(gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### - Test prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.69315661e-01]\n",
      " [  1.12877107e-01]\n",
      " [  8.30185683e-02]\n",
      " [  5.53457122e-02]\n",
      " [  5.80998131e-02]\n",
      " [  3.87332087e-02]\n",
      " [  2.84874020e-02]\n",
      " [  1.89916014e-02]\n",
      " [  3.96227820e-02]\n",
      " [  2.64151880e-02]\n",
      " [  1.94277755e-02]\n",
      " [  1.29518503e-02]\n",
      " [  1.35963574e-02]\n",
      " [  9.06423824e-03]\n",
      " [  6.66654294e-03]\n",
      " [  4.44436196e-03]\n",
      " [  2.66894686e-02]\n",
      " [  1.77929790e-02]\n",
      " [  1.30863351e-02]\n",
      " [  8.72422338e-03]\n",
      " [  9.15835621e-03]\n",
      " [  6.10557081e-03]\n",
      " [  4.49050972e-03]\n",
      " [  2.99367315e-03]\n",
      " [  6.24579552e-03]\n",
      " [  4.16386368e-03]\n",
      " [  3.06242789e-03]\n",
      " [  2.04161859e-03]\n",
      " [  2.14321316e-03]\n",
      " [  1.42880877e-03]\n",
      " [  1.05085665e-03]\n",
      " [  7.00571100e-04]\n",
      " [  1.78757231e-02]\n",
      " [  1.19171487e-02]\n",
      " [  8.76479430e-03]\n",
      " [  5.84319620e-03]\n",
      " [  6.13396401e-03]\n",
      " [  4.08930934e-03]\n",
      " [  3.00759486e-03]\n",
      " [  2.00506324e-03]\n",
      " [  4.18322721e-03]\n",
      " [  2.78881814e-03]\n",
      " [  2.05111289e-03]\n",
      " [  1.36740859e-03]\n",
      " [  1.43545327e-03]\n",
      " [  9.56968847e-04]\n",
      " [  7.03829019e-04]\n",
      " [  4.69219346e-04]\n",
      " [  2.81777567e-03]\n",
      " [  1.87851711e-03]\n",
      " [  1.38160700e-03]\n",
      " [  9.21071332e-04]\n",
      " [  9.66905476e-04]\n",
      " [  6.44603651e-04]\n",
      " [  4.74091457e-04]\n",
      " [  3.16060971e-04]\n",
      " [  6.59408059e-04]\n",
      " [  4.39605372e-04]\n",
      " [  3.23319844e-04]\n",
      " [  2.15546562e-04]\n",
      " [  2.26272542e-04]\n",
      " [  1.50848361e-04]\n",
      " [  1.10945570e-04]\n",
      " [  7.39637133e-05]\n",
      " [  1.19418639e-02]\n",
      " [  7.96124259e-03]\n",
      " [  5.85531449e-03]\n",
      " [  3.90354300e-03]\n",
      " [  4.09779022e-03]\n",
      " [  2.73186015e-03]\n",
      " [  2.00922157e-03]\n",
      " [  1.33948105e-03]\n",
      " [  2.79460192e-03]\n",
      " [  1.86306794e-03]\n",
      " [  1.37024448e-03]\n",
      " [  9.13496322e-04]\n",
      " [  9.58953520e-04]\n",
      " [  6.39302347e-04]\n",
      " [  4.70192467e-04]\n",
      " [  3.13461645e-04]\n",
      " [  1.88241300e-03]\n",
      " [  1.25494200e-03]\n",
      " [  9.22981556e-04]\n",
      " [  6.15321038e-04]\n",
      " [  6.45940504e-04]\n",
      " [  4.30627003e-04]\n",
      " [  3.16716455e-04]\n",
      " [  2.11144303e-04]\n",
      " [  4.40517077e-04]\n",
      " [  2.93678051e-04]\n",
      " [  2.15993588e-04]\n",
      " [  1.43995725e-04]\n",
      " [  1.51161208e-04]\n",
      " [  1.00774139e-04]\n",
      " [  7.41171078e-05]\n",
      " [  4.94114052e-05]\n",
      " [  1.26077795e-03]\n",
      " [  8.40518634e-04]\n",
      " [  6.18182512e-04]\n",
      " [  4.12121675e-04]\n",
      " [  4.32629581e-04]\n",
      " [  2.88419720e-04]\n",
      " [  2.12126204e-04]\n",
      " [  1.41417469e-04]\n",
      " [  2.95043765e-04]\n",
      " [  1.96695843e-04]\n",
      " [  1.44665360e-04]\n",
      " [  9.64435730e-05]\n",
      " [  1.01242776e-04]\n",
      " [  6.74951842e-05]\n",
      " [  4.96411868e-05]\n",
      " [  3.30941245e-05]\n",
      " [  1.98738223e-04]\n",
      " [  1.32492148e-04]\n",
      " [  9.74449891e-05]\n",
      " [  6.49633261e-05]\n",
      " [  6.81960164e-05]\n",
      " [  4.54640110e-05]\n",
      " [  3.34377554e-05]\n",
      " [  2.22918369e-05]\n",
      " [  4.65081685e-05]\n",
      " [  3.10054457e-05]\n",
      " [  2.28038065e-05]\n",
      " [  1.52025377e-05]\n",
      " [  1.59590429e-05]\n",
      " [  1.06393619e-05]\n",
      " [  7.82501092e-06]\n",
      " [  5.21667395e-06]\n",
      " [  7.96859197e-03]\n",
      " [  5.31239465e-03]\n",
      " [  3.90714653e-03]\n",
      " [  2.60476436e-03]\n",
      " [  2.73438205e-03]\n",
      " [  1.82292137e-03]\n",
      " [  1.34071758e-03]\n",
      " [  8.93811721e-04]\n",
      " [  1.86478783e-03]\n",
      " [  1.24319188e-03]\n",
      " [  9.14339612e-04]\n",
      " [  6.09559741e-04]\n",
      " [  6.39892516e-04]\n",
      " [  4.26595010e-04]\n",
      " [  3.13751016e-04]\n",
      " [  2.09167344e-04]\n",
      " [  1.25610049e-03]\n",
      " [  8.37400328e-04]\n",
      " [  6.15889068e-04]\n",
      " [  4.10592712e-04]\n",
      " [  4.31024534e-04]\n",
      " [  2.87349689e-04]\n",
      " [  2.11339220e-04]\n",
      " [  1.40892814e-04]\n",
      " [  2.93949159e-04]\n",
      " [  1.95966106e-04]\n",
      " [  1.44128654e-04]\n",
      " [  9.60857694e-05]\n",
      " [  1.00867168e-04]\n",
      " [  6.72447786e-05]\n",
      " [  4.94570191e-05]\n",
      " [  3.29713461e-05]\n",
      " [  8.41294555e-04]\n",
      " [  5.60863037e-04]\n",
      " [  4.12502123e-04]\n",
      " [  2.75001415e-04]\n",
      " [  2.88685974e-04]\n",
      " [  1.92457316e-04]\n",
      " [  1.41548018e-04]\n",
      " [  9.43653454e-05]\n",
      " [  1.96877422e-04]\n",
      " [  1.31251615e-04]\n",
      " [  9.65326044e-05]\n",
      " [  6.43550696e-05]\n",
      " [  6.75574921e-05]\n",
      " [  4.50383281e-05]\n",
      " [  3.31246752e-05]\n",
      " [  2.20831168e-05]\n",
      " [  1.32614458e-04]\n",
      " [  8.84096386e-05]\n",
      " [  6.50232966e-05]\n",
      " [  4.33488644e-05]\n",
      " [  4.55059808e-05]\n",
      " [  3.03373206e-05]\n",
      " [  2.23124155e-05]\n",
      " [  1.48749437e-05]\n",
      " [  3.10340682e-05]\n",
      " [  2.06893788e-05]\n",
      " [  1.52165718e-05]\n",
      " [  1.01443812e-05]\n",
      " [  1.06491836e-05]\n",
      " [  7.09945574e-06]\n",
      " [  5.22148969e-06]\n",
      " [  3.48099313e-06]\n",
      " [  5.62026219e-04]\n",
      " [  3.74684146e-04]\n",
      " [  2.75571745e-04]\n",
      " [  1.83714497e-04]\n",
      " [  1.92856456e-04]\n",
      " [  1.28570971e-04]\n",
      " [  9.45610512e-05]\n",
      " [  6.30407008e-05]\n",
      " [  1.31523819e-04]\n",
      " [  8.76825463e-05]\n",
      " [  6.44885366e-05]\n",
      " [  4.29923577e-05]\n",
      " [  4.51317337e-05]\n",
      " [  3.00878225e-05]\n",
      " [  2.21289153e-05]\n",
      " [  1.47526102e-05]\n",
      " [  8.85929927e-05]\n",
      " [  5.90619952e-05]\n",
      " [  4.34387663e-05]\n",
      " [  2.89591775e-05]\n",
      " [  3.04002376e-05]\n",
      " [  2.02668251e-05]\n",
      " [  1.49057931e-05]\n",
      " [  9.93719538e-06]\n",
      " [  2.07322868e-05]\n",
      " [  1.38215245e-05]\n",
      " [  1.01654198e-05]\n",
      " [  6.77694653e-06]\n",
      " [  7.11417942e-06]\n",
      " [  4.74278628e-06]\n",
      " [  3.48821242e-06]\n",
      " [  2.32547494e-06]\n",
      " [  5.93366557e-05]\n",
      " [  3.95577704e-05]\n",
      " [  2.90938486e-05]\n",
      " [  1.93958991e-05]\n",
      " [  2.03610734e-05]\n",
      " [  1.35740490e-05]\n",
      " [  9.98340707e-06]\n",
      " [  6.65560472e-06]\n",
      " [  1.38857998e-05]\n",
      " [  9.25719989e-06]\n",
      " [  6.80846188e-06]\n",
      " [  4.53897459e-06]\n",
      " [  4.76484202e-06]\n",
      " [  3.17656135e-06]\n",
      " [  2.33628928e-06]\n",
      " [  1.55752618e-06]\n",
      " [  9.35332147e-06]\n",
      " [  6.23554765e-06]\n",
      " [  4.58610476e-06]\n",
      " [  3.05740317e-06]\n",
      " [  3.20954498e-06]\n",
      " [  2.13969665e-06]\n",
      " [  1.57369866e-06]\n",
      " [  1.04913244e-06]\n",
      " [  2.18883839e-06]\n",
      " [  1.45922559e-06]\n",
      " [  1.07322754e-06]\n",
      " [  7.15485024e-07]\n",
      " [  7.51088828e-07]\n",
      " [  5.00725885e-07]\n",
      " [  3.68272603e-07]\n",
      " [  2.45515069e-07]]\n"
     ]
    }
   ],
   "source": [
    "prob = np.ones((kbar2,1))\n",
    "\n",
    "for i in range(0,kbar2):\n",
    "    for m in range(1,kbar+1):\n",
    "        prob[i,0]=prob[i,0]*gamma[kbar1-m-1, int(bitget(i,m))]\n",
    "print(prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df_prob = df(prob)\n",
    "df_prob.to_excel('prob.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00838513],\n",
       "       [ 0.00838513],\n",
       "       [ 0.00838513],\n",
       "       [ 0.00838513],\n",
       "       [ 0.00838504],\n",
       "       [ 0.00838504],\n",
       "       [ 0.00838504],\n",
       "       [ 0.00838504],\n",
       "       [ 0.00838028],\n",
       "       [ 0.00838028],\n",
       "       [ 0.00838028],\n",
       "       [ 0.00838028],\n",
       "       [ 0.00838019],\n",
       "       [ 0.00838019],\n",
       "       [ 0.00838019],\n",
       "       [ 0.00838019],\n",
       "       [ 0.00831208],\n",
       "       [ 0.00831208],\n",
       "       [ 0.00831208],\n",
       "       [ 0.00831208],\n",
       "       [ 0.00831199],\n",
       "       [ 0.00831199],\n",
       "       [ 0.00831199],\n",
       "       [ 0.00831199],\n",
       "       [ 0.00830727],\n",
       "       [ 0.00830727],\n",
       "       [ 0.00830727],\n",
       "       [ 0.00830727],\n",
       "       [ 0.00830719],\n",
       "       [ 0.00830719],\n",
       "       [ 0.00830719],\n",
       "       [ 0.00830719],\n",
       "       [ 0.00794822],\n",
       "       [ 0.00794822],\n",
       "       [ 0.00794822],\n",
       "       [ 0.00794822],\n",
       "       [ 0.00794814],\n",
       "       [ 0.00794814],\n",
       "       [ 0.00794814],\n",
       "       [ 0.00794814],\n",
       "       [ 0.00794362],\n",
       "       [ 0.00794362],\n",
       "       [ 0.00794362],\n",
       "       [ 0.00794362],\n",
       "       [ 0.00794354],\n",
       "       [ 0.00794354],\n",
       "       [ 0.00794354],\n",
       "       [ 0.00794354],\n",
       "       [ 0.00787897],\n",
       "       [ 0.00787897],\n",
       "       [ 0.00787897],\n",
       "       [ 0.00787897],\n",
       "       [ 0.00787889],\n",
       "       [ 0.00787889],\n",
       "       [ 0.00787889],\n",
       "       [ 0.00787889],\n",
       "       [ 0.00787441],\n",
       "       [ 0.00787441],\n",
       "       [ 0.00787441],\n",
       "       [ 0.00787441],\n",
       "       [ 0.00787434],\n",
       "       [ 0.00787434],\n",
       "       [ 0.00787434],\n",
       "       [ 0.00787434],\n",
       "       [ 0.0070083 ],\n",
       "       [ 0.0070083 ],\n",
       "       [ 0.0070083 ],\n",
       "       [ 0.0070083 ],\n",
       "       [ 0.00700823],\n",
       "       [ 0.00700823],\n",
       "       [ 0.00700823],\n",
       "       [ 0.00700823],\n",
       "       [ 0.00700424],\n",
       "       [ 0.00700424],\n",
       "       [ 0.00700424],\n",
       "       [ 0.00700424],\n",
       "       [ 0.00700417],\n",
       "       [ 0.00700417],\n",
       "       [ 0.00700417],\n",
       "       [ 0.00700417],\n",
       "       [ 0.00694724],\n",
       "       [ 0.00694724],\n",
       "       [ 0.00694724],\n",
       "       [ 0.00694724],\n",
       "       [ 0.00694717],\n",
       "       [ 0.00694717],\n",
       "       [ 0.00694717],\n",
       "       [ 0.00694717],\n",
       "       [ 0.00694322],\n",
       "       [ 0.00694322],\n",
       "       [ 0.00694322],\n",
       "       [ 0.00694322],\n",
       "       [ 0.00694315],\n",
       "       [ 0.00694315],\n",
       "       [ 0.00694315],\n",
       "       [ 0.00694315],\n",
       "       [ 0.00664313],\n",
       "       [ 0.00664313],\n",
       "       [ 0.00664313],\n",
       "       [ 0.00664313],\n",
       "       [ 0.00664306],\n",
       "       [ 0.00664306],\n",
       "       [ 0.00664306],\n",
       "       [ 0.00664306],\n",
       "       [ 0.00663928],\n",
       "       [ 0.00663928],\n",
       "       [ 0.00663928],\n",
       "       [ 0.00663928],\n",
       "       [ 0.00663922],\n",
       "       [ 0.00663922],\n",
       "       [ 0.00663922],\n",
       "       [ 0.00663922],\n",
       "       [ 0.00658525],\n",
       "       [ 0.00658525],\n",
       "       [ 0.00658525],\n",
       "       [ 0.00658525],\n",
       "       [ 0.00658519],\n",
       "       [ 0.00658519],\n",
       "       [ 0.00658519],\n",
       "       [ 0.00658519],\n",
       "       [ 0.00658144],\n",
       "       [ 0.00658144],\n",
       "       [ 0.00658144],\n",
       "       [ 0.00658144],\n",
       "       [ 0.00658138],\n",
       "       [ 0.00658138],\n",
       "       [ 0.00658138],\n",
       "       [ 0.00658138],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039463],\n",
       "       [ 0.00039441],\n",
       "       [ 0.00039441],\n",
       "       [ 0.00039441],\n",
       "       [ 0.00039441],\n",
       "       [ 0.0003944 ],\n",
       "       [ 0.0003944 ],\n",
       "       [ 0.0003944 ],\n",
       "       [ 0.0003944 ],\n",
       "       [ 0.0003912 ],\n",
       "       [ 0.0003912 ],\n",
       "       [ 0.0003912 ],\n",
       "       [ 0.0003912 ],\n",
       "       [ 0.00039119],\n",
       "       [ 0.00039119],\n",
       "       [ 0.00039119],\n",
       "       [ 0.00039119],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00039097],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037407],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037385],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.00037081],\n",
       "       [ 0.0003706 ],\n",
       "       [ 0.0003706 ],\n",
       "       [ 0.0003706 ],\n",
       "       [ 0.0003706 ],\n",
       "       [ 0.00037059],\n",
       "       [ 0.00037059],\n",
       "       [ 0.00037059],\n",
       "       [ 0.00037059],\n",
       "       [ 0.00032984],\n",
       "       [ 0.00032984],\n",
       "       [ 0.00032984],\n",
       "       [ 0.00032984],\n",
       "       [ 0.00032983],\n",
       "       [ 0.00032983],\n",
       "       [ 0.00032983],\n",
       "       [ 0.00032983],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032964],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032696],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00032677],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031265],\n",
       "       [ 0.00031247],\n",
       "       [ 0.00031247],\n",
       "       [ 0.00031247],\n",
       "       [ 0.00031247],\n",
       "       [ 0.00031246],\n",
       "       [ 0.00031246],\n",
       "       [ 0.00031246],\n",
       "       [ 0.00031246],\n",
       "       [ 0.00030993],\n",
       "       [ 0.00030993],\n",
       "       [ 0.00030993],\n",
       "       [ 0.00030993],\n",
       "       [ 0.00030992],\n",
       "       [ 0.00030992],\n",
       "       [ 0.00030992],\n",
       "       [ 0.00030992],\n",
       "       [ 0.00030975],\n",
       "       [ 0.00030975],\n",
       "       [ 0.00030975],\n",
       "       [ 0.00030975],\n",
       "       [ 0.00030974],\n",
       "       [ 0.00030974],\n",
       "       [ 0.00030974],\n",
       "       [ 0.00030974]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# TEST with hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n",
      "[[ 0.75     0.25   ]\n",
      " [ 0.625    0.375  ]\n",
      " [ 0.53125  0.46875]]\n"
     ]
    }
   ],
   "source": [
    "kbar = 3\n",
    "input_param = [2, 1, 0.5, 1]\n",
    "A = T_mat_temp(kbar)\n",
    "TT = transition_mat(A,input_param,kbar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.24902344,  0.21972656,  0.14941406,  0.13183594,  0.08300781,\n",
       "         0.07324219,  0.04980469,  0.04394531],\n",
       "       [ 0.21972656,  0.24902344,  0.13183594,  0.14941406,  0.07324219,\n",
       "         0.08300781,  0.04394531,  0.04980469],\n",
       "       [ 0.14941406,  0.13183594,  0.24902344,  0.21972656,  0.04980469,\n",
       "         0.04394531,  0.08300781,  0.07324219],\n",
       "       [ 0.13183594,  0.14941406,  0.21972656,  0.24902344,  0.04394531,\n",
       "         0.04980469,  0.07324219,  0.08300781],\n",
       "       [ 0.08300781,  0.07324219,  0.04980469,  0.04394531,  0.24902344,\n",
       "         0.21972656,  0.14941406,  0.13183594],\n",
       "       [ 0.07324219,  0.08300781,  0.04394531,  0.04980469,  0.21972656,\n",
       "         0.24902344,  0.13183594,  0.14941406],\n",
       "       [ 0.04980469,  0.04394531,  0.08300781,  0.07324219,  0.14941406,\n",
       "         0.13183594,  0.24902344,  0.21972656],\n",
       "       [ 0.04394531,  0.04980469,  0.07324219,  0.08300781,  0.13183594,\n",
       "         0.14941406,  0.21972656,  0.24902344]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.],\n",
       "       [ 1.,  0.,  3.,  2.,  5.,  4.,  7.,  6.],\n",
       "       [ 2.,  3.,  0.,  1.,  6.,  7.,  4.,  5.],\n",
       "       [ 3.,  2.,  1.,  0.,  7.,  6.,  5.,  4.],\n",
       "       [ 4.,  5.,  6.,  7.,  0.,  1.,  2.,  3.],\n",
       "       [ 5.,  4.,  7.,  6.,  1.,  0.,  3.,  2.],\n",
       "       [ 6.,  7.,  4.,  5.,  2.,  3.,  0.,  1.],\n",
       "       [ 7.,  6.,  5.,  4.,  3.,  2.,  1.,  0.]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}