{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov examples\n",
    "\n",
    "This notebook aims at illustrating the theory presented in [Lecture 9 Markov Processes](http://randall-romero.com/wp-content/uploads/Macro2-2017b/handouts/Lecture-9-Markov-Process.pdf) of my class on Macroeconomics II. It presents three examples of Markov chain processes:\n",
    " * Unemployment rate\n",
    " * Marital status\n",
    " * Credit ratings\n",
    "\n",
    "For each of these examples, we simulate the Markov chain and the evolution of the distribution of states. \n",
    "\n",
    "*Randall Romero-Aguilar, PhD*\n",
    "\n",
    "October 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from warnings import warn\n",
    "import seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "We first write up a class to represent the Markov chain, and to define how to plot its simulations and distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class MarkovChain:\n",
    "\n",
    "    def __init__(self, P, pi0=None, z=None):\n",
    "        \"\"\"\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        P   transition probability matrix, n.n\n",
    "        pi0 initial distribution, n vector\n",
    "        z   stochastic outcomes, n vector\n",
    "        \"\"\"\n",
    "        P = np.atleast_2d(P)\n",
    "        n, n2 = P.shape\n",
    "        assert n == n2, \"P must be a square array\"\n",
    "        assert np.all(P >= 0), \"P elements must be nonnegative\"\n",
    "        assert np.allclose(P.sum(1), 1), \"The rows of P must sum to 1\"\n",
    "\n",
    "        self.P = P\n",
    "        self.n = n\n",
    "        self.cdf = P.cumsum(1)\n",
    "\n",
    "        if pi0 is None:\n",
    "            self.pi0 = self.stationary_distribution if self.has_unique_stationary_distribution else np.ones((n))/n\n",
    "        else:\n",
    "            pi0 = np.asarray(pi0)\n",
    "            assert pi0.ndim == 1 and pi0.size == n, \"pi0 must be a vector with %d elements\" % n\n",
    "            assert np.all(pi0 >= 0), \"pi0 elements must be nonnegative\"\n",
    "            assert np.allclose(pi0.sum(), 1), \"The elements of pi0 must sum to 1\"\n",
    "            self.pi0 = pi0\n",
    "\n",
    "        if z is None:\n",
    "            self.outcomes = np.arange(n)\n",
    "        else:\n",
    "            z = np.asarray(z)\n",
    "            assert z.ndim == 1 and z.size == n, \"z must be a vector with %d elements\" % n\n",
    "            self.outcomes = z\n",
    "\n",
    "        self.string_outcomes = isinstance(self.outcomes[0], str)\n",
    "\n",
    "    @property\n",
    "    def stationary_distribution(self):\n",
    "        if self.has_unique_stationary_distribution:\n",
    "            A = np.vstack((np.eye(self.n)-self.P.T, np.ones((1,self.n))))\n",
    "            b = np.hstack((np.zeros(self.n), np.ones((1))))\n",
    "            return np.linalg.lstsq(A,b)[0]\n",
    "        else:\n",
    "            warn(\"This Markov process does not have a unique stationary distribution\")\n",
    "\n",
    "    @property\n",
    "    def has_unique_stationary_distribution(self):  #fixme: ugly implementation\n",
    "        pn = self.P\n",
    "        for it in range(50):\n",
    "            if np.any(pn == 0.0):\n",
    "                pn = pn @ self.P\n",
    "            else:\n",
    "                return True\n",
    "\n",
    "        return False\n",
    "\n",
    "    def rand(self, j):\n",
    "        \"\"\"\n",
    "        Draw next state\n",
    "        Parameters\n",
    "        ----------\n",
    "        j current state, integer = 0, 1,...,n-1\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        Integer: next state\n",
    "        \"\"\"\n",
    "        return (np.random.rand() > self.cdf[j]).sum()\n",
    "\n",
    "    def simulate(self, T: int = 120, j: int = 0):\n",
    "        \"\"\"\n",
    "        Simulate the Markov process\n",
    "        Parameters\n",
    "        ----------\n",
    "        T\n",
    "        j\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        Vector: array of simulated values\n",
    "        \"\"\"\n",
    "        assert j in range(self.n), 'Initial state must be a nonnegative integer less than %d' % self.n\n",
    "\n",
    "        y = np.empty((1 + T))\n",
    "        y[0] = j if self.string_outcomes else self.outcomes[j]\n",
    "        for t in range(T):\n",
    "            j = self.rand(j)\n",
    "            y[t + 1] = j if self.string_outcomes else self.outcomes[j]\n",
    "\n",
    "        return y\n",
    "\n",
    "    def unconditional_distribution(self, T=24):\n",
    "        prob = np.empty((T+1, self.n))\n",
    "        prob[0] = self.pi0\n",
    "        for t in range(T):\n",
    "            prob[t + 1] = np.dot(prob[t], self.P)\n",
    "        return prob\n",
    "\n",
    "\n",
    "    def plot(self, ax, T=120, j=0, *args):\n",
    "        y = self.simulate(T, j)\n",
    "        ax.plot(y, '-o', drawstyle='steps', *args)\n",
    "\n",
    "        if isinstance(self.outcomes[0], str):\n",
    "            ax.set_ylim([-0.5, self.n - 0.5])\n",
    "            ax.set_yticks(np.arange(self.n))\n",
    "            ax.set_yticklabels(self.outcomes)\n",
    "        else:\n",
    "            ax.set_yticks(np.sort(self.outcomes))\n",
    "\n",
    "    def plot_probability(self, ax, T=24):\n",
    "        p = self.unconditional_distribution(T).T\n",
    "        lb = np.zeros((T+1))\n",
    "\n",
    "        current_palette = seaborn.color_palette(n_colors=self.n)\n",
    "\n",
    "        tt = np.arange(1+T)\n",
    "        for k in range(self.n):\n",
    "            ax.bar(tt, p[k], bottom=lb, color=current_palette[k], label=str(self.outcomes[k]))\n",
    "            lb += p[k]\n",
    "        ax.set_ylim([0, 1.25])\n",
    "        ax.legend(loc='upper right', ncol=self.n)\n",
    "\n",
    "    def print_stationary(self):\n",
    "        if self.has_unique_stationary_distribution:\n",
    "            for z, p in zip(self.outcomes, self.stationary_distribution):\n",
    "                print(f'{z:6s}   {100*p:9.2f}')\n",
    "    \n",
    "    def __repr__(self):\n",
    "        txt = 'A Markov Chain process with transition probability\\n'\n",
    "        txt += self.P.__str__()\n",
    "        txt += '\\n\\nand possible outcomes\\n'\n",
    "        txt += self.outcomes.__str__()\n",
    "        return txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the simulation horizon to $T=48$ periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T=48"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: unemployment\n",
    "\n",
    "A worker can either be employed or unemployed:\n",
    "\n",
    "* If unemployed, she will get a job with probability $p = 45\\%$\n",
    "* If employed, she will lose her job with probability $q = 5\\%$\n",
    "\n",
    "The worker is unemployed at $t = 0$. Then the Markov chain is:\n",
    "\n",
    "* outcomes  \\{unemployed, employed\\}  or $z = \\begin{bmatrix}0 \\\\ 1 \\end{bmatrix}$\n",
    "* initial probability $\\pi_0 = \\begin{bmatrix}1 \\\\ 0 \\end{bmatrix}$\n",
    "* transition probability $P = \\begin{bmatrix}1-p & p \\\\ q & 1-q \\end{bmatrix} = \\begin{bmatrix}0.55 & 0.45 \\\\ 0.05 & 0.95 \\end{bmatrix}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unemployed       10.00\n",
      "employed       90.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAHVCAYAAACexG8YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X20XVV9L/zvjyS8tCAoiVQJNFgBoQZBQ7QXCxQV8OUC\nvlSlSpVK0V6xalsq9unA4ONtuep4amu5RVTEKgqKiAzFiyIiRS0kFIQKQhGjBrgSAuIbCCTz+eNs\n4iEk5u28zX0+nzEyctZac635O3vunP3dM3PtU621AAAAfdpisgsAAAA2nUAPAAAdE+gBAKBjAj0A\nAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADo2MzJLqA3s2fPbvPmzZvsMgAAGHJXX331Xa21\nOetrJ9BvpHnz5mXJkiWTXQYAAEOuqr6/Ie0suQEAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAP\nAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjAj0A\nAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA\n0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBA\nxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAx2ZOdgEbqqp+1lrb\ndhyvf1aSz7fWzhuvPibTBdfclvdcfFNu//F9eeIO2+TEw/bMUfvtPKbnTMU+pmJNE2Ui6hqW8dDH\n1OljKtakD33oo6+apqMZixYtmuwaNsgpp5zyN4sWLfq7cbz+UUluXrRo0Q2/rt0ZZ5yx6Pjjjx+v\nMsbFBdfclreff33u/sUDSZKf3v9Qvnbz8sx97DZ5yhMeMybnTMU+pmJNE2Ui6hqW8dDH1OljKtak\nD33oo6+ahs0pp5xyx6JFi85YX7tqrY1551X16iR/nmTLJFcm+R9J7k1yWpLnJrknyd8keXeSXZO8\npbV2YVW9NsmLk2yVZLckn2itnTK45s9aa9tWVQ3Oe36SluRdrbVzq+pjSc5rrX1u0P7sJOcm+UKS\nU5McPLjuaa21Dwyu8/4khyT5XpJKcub6ZugXLFjQlixZsvkP0hj4xJU/yOeuvW297a75wY/zwMpV\nj9q/5Ywtst+uO4zJOVOxj8msaecdtsnXTzpkrX1sqg0d782payz6mIrjoY8++piKNelDH/qYmjWN\nx+vsVFRVV7fWFqyv3Zivoa+qvZK8IskBrbV9k6xM8qokv5nkstbaM5L8NMm7kjwvIwH+naMusXDQ\nft8kf1hVa34TLxkce1pG3hy8p6qekORDSY4d1LB9kv+W5KIkr0tyb2tt/yT7J/nTqtpt0O+eSeYn\n+dNB+3V9T8dX1ZKqWrJ8+fJNelzGw+euvS033PGT9bZb2z+EX7d/U86Zin1MZk23//i+dfaxqTZ0\nvJNNr2ss+piK46GPPvqYijXpQx/6mJo1jcfrbM/GYw39c5I8I8nikUnwbJPkziQPJPk/gzbXJ/ll\na+3Bqro+ybxR53+5tbYiSarq/CTPTjJ6SvzZST7ZWluZ5EdV9bUk+w9m+E+rqsdnJPR/prX2UFUd\nmmSfqnrZ4Pztk+ye5MBR17m9qi5d1zfUWjsjyRnJyAz9Jj0q42TvJzwm577+935tmwNOvTS3reWJ\nv/MO26zz3I09Zyr2MZk1PXGHbdZ6/c21IeOdbF5dm9vHVBwPffTRx1SsSR/60MfUrGm8Xmd7NR6f\nclNJPtpa23fwZ8/W2qIkD7Zfre9ZleSXSdJaW5VHvrFYMzCvuV2/pu+PZWR2/9gkHxnV/k2j6tmt\ntfaldVx7KJ142J7ZZtaMR+zbZtaMnHjYnmN2zlTsYyrWNFEmoq5hGQ99TJ0+pmJN+tCHPvqqaboa\njxn6ryT5XFX9Q2vtzqp6XJLtNuL85w3OuS/JUUn+ZI3jlyd5fVV9NMnjMjLTfuLg2FlJrkryf1tr\n3x7suzjJn1XVpYP/EdgjyW2jrvOvSR6f5A+SfGIjv9cuPHwn+MbcIb6x50zFPiaypr8+77o8sHJV\ndp4id99PRF1TeTz00WcfU7EmfehDH1Ojpqn2OjvVjNdNsa9I8vaM/A/Ag0nemOSShz92sqoWJflZ\na+29g+2Hb3h9bZIXZGS9/ZOzETfFjur7/yS5oLV2+mB7i4ys1//vGZmtX56RNwo/ya9uir15cPrH\ne7op9hUf+GaSbNDyCMbXRIzFpvSxsed4TgEwFU3X16cNvSl2XD6HfhCwz11j97ajji9ao/3oz5e/\ns7V2wlquue3g75aRGfkT12xTVb+RkfXxnxx13qqMfKLO36yl1Ef1AwAAPRma3xRbVc9N8p0k72+t\n3TvZ9QAAwESYUr8ptrV2VkbWwW/KuZdk5DPtAQBg2hiaGXoAAJiOBHoAAOiYQA8AAB0T6AEAoGMC\nPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0\nAADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtAD\nAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8A\nAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAA\ndEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQ\nMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDH\nBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T\n6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEyg\nBwCAjgn0AADQMYEeAAA61k2gr6qDq+rz49zH0qqaPZ59jJULrrkt1/zgx7nye3fngFMvzQXX3DbZ\nJbERLrjmthxw6qXZ7aQvbND4TcR4e04BMCw25XV2Y9pPNTMnuwA23gXX3Ja3n399Hli5Kkly24/v\ny9vPvz5JctR+O09maWyAh8fvvgdXJln/+E3EeHtOATAsNvV1dkPbT0XrDfRVNS/J51trTx1s/1WS\nbZMcnOTKJH+QZIckr2ut/VtVzUhy6uD4VklOa619oKoOTnJKkh8l2TfJ+UmuT/LmJNskOaq19t2q\nOivJ/Ul+N8lOSf6itfaImfmqelySM5M8Kckvkhyf5D+T3JTkv7XWllfVFkluTvKsJJXk9CS7Di7x\nltba16tqxySfTDInyVWDdlPeey6+afWT7mH3Pbgy77n4pm6eeMPohjt+kld84JvrbXfND368Ojg/\n7L4HV+avz7sun7zqBxvcfkPHe0Pq2tw+AGC8TdXX2algc2foZ7bWFlbVC5K8I8lzk7wuyb2ttf2r\naqskX6+qLw3aPy3JXknuTnJrkg8Nzn9zkjclecug3bwkByX5nSRfraonr9HvKUmuaa0dVVWHJPnX\n1tq+VfXxJK9K8r5BLd9qrd1VVZ9I8g+ttSuqatckFw/qeEeSK1pr76yqF2bkjcGjVNXxDx/bdddd\n19ZkQt3+4/s2aj/j78h9N/wf/Jo/NDZ1/4aM94bWtTl9AMB4m6qvs1PF5gb68wd/X52REJ4khybZ\np6peNtjePsnuSR5Isri1dkeSVNV3kzwc9K/PyEz/wz7VWluV5L+q6tYkT1mj32cneWmStNYuraod\nq2r7jMzafy4jgf5Pknxk0P65SfauWj0B/5iq2i7JgUleMrjOF6rqnrV9k621M5KckSQLFixo63lM\nxt0Td9gmt63lSfbEHbaZhGpIkj965q75o2du2Ju9A069dK3jt/MO2+Tc1//eBrffkPHe0Lo2pw8A\nGG9T9XV2qtiQm2IfWqPd1qO+/uXg75X51ZuDSvKm1tq+gz+7tda+tEb7JFk1antVHvnmYs3QvOb2\n2pbGtNbaD5P8aDBr/8wkXxwc2yLJ742qaefW2k/Xce0p78TD9sw2s2Y8Yt82s2bkxMP2nKSK2Bgb\nO34TMd6eUwAMi6n4OjveNiTQ/yjJ4wez4FsledF62l+c5M+qalaSVNUeVfWbG1nXH1bVFlX1OxlZ\nJ3/TGscvz8jSmgzW5t/VWvvJ4NiHknw8I7P8Dy80/1KSEx4+uar2Xct1np/ksRtZ56Q4ar+d8/cv\nmZ+dd9gmlZF3nH//kvndrPOa7jZ2/CZivD2nABgWU/F1drxVa+ufoK6qP0/y50m+l+S2JEszctPr\nX7XWlgw+6nFJa23e4GbUdyX57xmZSV+e5Kgk+w3av2hwzctGnf/wtV40uCn2niQLMuqm2DXaPC4j\ny2l2y+Cm2NbadYPrzkqyIsnC1tp3BvtmJzktI+vmZya5vLX2hlE3xc5O8rWMLL95RmvtrnU9FgsW\nLGhLlixZ72MGAACbo6qubq0tWG+7DQn0E2kQ6D/fWjtvE89fkJEbYH9/TAsbEOgBAJgIGxroh+pz\n6KvqpCR/lsEyGgAAGHZTLtC31l67GeeempHPwAcAgGlhQ26KBQAApiiBHgAAOibQAwBAxwR6AADo\nmEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBj\nAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J\n9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQ\nAwBAxwR6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOlattcmuoStVtTzJ9ye7\njlFmJ7lrsotgwhjv6ceYTz/GfPox5tPLxoz3b7fW5qyvkUDfuapa0lpbMNl1MDGM9/RjzKcfYz79\nGPPpZTzG25IbAADomEAPAAAdE+j7d8ZkF8CEMt7TjzGffoz59GPMp5cxH29r6AEAoGNm6AEAoGMC\nPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0\nAADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtAD\nAEDHBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8A\nAB0T6AEAoGMzJ7uA3syePbvNmzdvsssAAGDIXX311Xe11uasr51Av5HmzZuXJUuWTHYZAAAMuar6\n/oa0s+QGAAA6JtADAEDHBHoAAOjY0K6hr6ozk7woyZ2ttaeu5firkrxtsPmzJH/WWvvWBJYIANPe\ngw8+mGXLluX++++f7FJg0my99daZO3duZs2atUnnD22gT3JWkn9O8q/rOP69JAe11u6pqucnOSPJ\nMyeoNgAgybJly7Lddttl3rx5qarJLgcmXGstK1asyLJly7Lbbrtt0jWGdslNa+3yJHf/muPfaK3d\nM9j89yRzJ6QwAGC1+++/PzvuuKMwz7RVVdlxxx0363+phjbQb6TXJfniug5W1fFVtaSqlixfvnwC\nywKA4SfMM91t7r+BaR/oq+oPMhLo37auNq21M1prC1prC+bMWe9n+wMAwIQZ5jX061VV+yT5UJLn\nt9ZWTHY9ADDdzTvpC2N6vaWnvnBMrzdeLrvssrz3ve/N5z//+XHr4+Ffjjl79uyxv/ii7cf4eveO\n7fU20Lbbbpuf/exn43b91772tXnRi16Ul73sZWN63Wk7Q19VuyY5P8kxrbWbJ7seAADYFEMb6Kvq\nk0m+mWTPqlpWVa+rqjdU1RsGTU5OsmOS/11V11bVkkkrFgCYFEuXLs1Tn/qrT7d+73vfm0WLFuXg\ngw/O2972tixcuDB77LFH/u3f/i1JsnLlypx44onZf//9s88+++QDH/hAkpEZ9oMOOigvf/nLs8ce\ne+Skk07K2WefnYULF2b+/Pn57ne/m2RkhvYNb3hDfv/3fz977LHHWmfk77777hx11FHZZ5998qxn\nPSvXXXddVq1ald133z0P38u3atWqPPnJT85dd92V5cuX56UvfWn233//7L///vn617+eJFmxYkUO\nPfTQ7Lfffnn961+f1tq4PpYT7eMf/3gWLlyYfffdN69//euzcuXKbLvttnnb296WZzzjGXnuc5+b\nq666KgcffHCe9KQn5cILL0ySnHXWWTnyyCNz+OGHZ88998wpp5zyqGu31nLiiSfmqU99aubPn59z\nzz03SXLMMcfkc5/73Op2r3rVq3LhhReu83nRWssJJ5yQvffeOy984Qtz5513jstjMbSBvrV2dGvt\nCa21Wa21ua21D7fWTm+tnT44flxr7bGttX0HfxZMds0AwNTx0EMP5aqrrsr73ve+1aHvwx/+cLbf\nfvssXrw4ixcvzgc/+MF873vfS5J861vfyj/+4z/m+uuvz8c+9rHcfPPNueqqq3Lcccfl/e9//+rr\nLl26NF/72tfyhS98IW94wxse9ekm73jHO7Lffvvluuuuy9/93d/lj//4j7PFFlvk1a9+dc4+++wk\nySWXXJKnPe1pmT17dt785jfnrW99axYvXpzPfOYzOe6445Ikp5xySp797GfnmmuuyRFHHJEf/OAH\nE/GwTYgbb7wx5557br7+9a/n2muvzYwZM3L22Wfn5z//eQ4++OBcffXV2W677fK3f/u3+fKXv5zP\nfvazOfnkk1eff9VVV+Xss8/Otddem09/+tNZsuSR87rnn39+rr322nzrW9/KJZdckhNPPDF33HFH\njjvuuHzkIx9Jktx77735xje+kRe84AXrfF589rOfzU033ZTrr78+H/zgB/ONb3xjXB6Pab2GHgBg\nXV7ykpckSZ7xjGdk6dKlSZIvfelLue6663LeeeclGQl1//Vf/5Utt9wy+++/f57whCckSX7nd34n\nhx56aJJk/vz5+epXv7r6ui9/+cuzxRZbZPfdd8+TnvSkfOc733lEv1dccUU+85nPJEkOOeSQrFix\nIvfee2/+5E/+JEceeWTe8pa35Mwzz8yxxx6bZCTc33DDDavP/8lPfpKf/vSnufzyy3P++ecnSV74\nwhfmsY997Fg/RJPmK1/5Sq6++ursv//+SZL77rsvj3/847Plllvm8MMPTzLyuG+11VaZNWtW5s+f\nv3oMk+R5z3tedtxxxyQj43zFFVdkwYJfze1eccUVOfroozNjxozstNNOOeigg7J48eIcccQReeMb\n35g777wz559/fl760pdm5syZ63xeXH755auv88QnPjGHHHLIuDweAj0AMG3NnDkzq1atWr09erZ8\nq622SpLMmDEjDz30UJKRJRTvf//7c9hhhz3iOpdddtnq9kmyxRZbrN7eYostVp+fPPojCtfcXtvS\nmKrKLrvskp122imXXnpprrzyytWz9atWrco3v/nNbLPNNms9bxi11vKa17wmf//3f/+I/e9973tX\nf89jPQYPO+aYY3L22WfnnHPOyZlnnrm6/dqeFxdddNGEjMHQLrkBAFifnXbaKXfeeWdWrFiRX/7y\nl+v9lJnDDjss//Iv/5IHH3wwSXLzzTfn5z//+Ub1+elPfzqrVq3Kd7/73dx6663Zc889H3H8wAMP\nXB3WL7vsssyePTuPecxjkiTHHXdcXv3qV+flL395ZsyYkSQ59NBD88///M+rz7/22msfdZ0vfvGL\nueeeezIsnvOc5+S8885bvSb97rvvzve///0NPv/LX/5y7r777tx333254IILcsABBzzi+IEHHphz\nzz03K1euzPLly3P55Zdn4cKFSUbug3jf+96XJPnd3/3dJOt+Xhx44IE555xzsnLlytxxxx2P+J+a\nsWSGHgCYMib6YyZnzZqVk08+Oc985jOz22675SlPecqvbX/cccdl6dKlefrTn57WWubMmZMLLrhg\no/rcc889c9BBB+VHP/pRTj/99Gy99daPOL5o0aIce+yx2WefffIbv/Eb+ehHP7r62BFHHJFjjz12\n9XKbJPmnf/qnvPGNb8w+++yThx56KAceeGBOP/30vOMd78jRRx+dpz/96TnooIOy6667blSdG2WC\nP2Zy7733zrve9a4ceuihWbVqVWbNmpXTTjttg89/9rOfnWOOOSa33HJL/uiP/ugRy22S5MUvfnG+\n+c1v5mlPe1qqKu9+97vzW7/1W0lG3gTutddeOeqoo1a3X9fz4sUvfnEuvfTSzJ8/P3vssUcOOuig\nsXkA1lDDdsfzeFuwYEFb88YJAGDT3Hjjjdlrr70mu4wJs7mfQ75kyZK89a1vXf2pO2y8s846K0uW\nLHnE/2psjF/84heZP39+/uM//iPbbz92n7+/tn8LVXX1hnxwiyU3AAAdOPXUU/PSl770UevGmTiX\nXHJJnvKUp+RNb3rTmIb5zWWGfiOZoQeAsTPdZuhhXczQAwDdMrnIdLe5/wYEegBg0my99dZZsWKF\nUM+01VrLihUrHnVz9MbwKTcAwKSZO3duli1bluXLl092KTBptt5668ydO3eTzxfoAYBJM2vWrOy2\n226TXQZ0zZIbAADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADo\nmEAPAAAdE+gBAKBjQxvoq+rMqrqzqv5zHcerqv6pqm6pquuq6ukTXSMAAGyuoQ30Sc5KcvivOf78\nJLsP/hyf5F8moCYAABhTQxvoW2uXJ7n71zQ5Msm/thH/nmSHqnrCxFQHAABjY2gD/QbYOckPR20v\nG+wDAIBuzJzsAiZRrWVfW2vDquMzsiwnu+6663jW9GvNO+kLG9Ru6akvHPli0fYbduFF9476eiPP\nmYp9TMWa9KEPffRVkz70oY++apqoPqao6TxDvyzJLqO25ya5fW0NW2tntNYWtNYWzJkzZ0KKAwCA\nDTGdA/2FSf548Gk3z0pyb2vtjskuCgAANsbQLrmpqk8mOTjJ7KpaluQdSWYlSWvt9CQXJXlBkluS\n/CLJsZNTKQAAbLqhDfSttaPXc7wleeMElQMAAONiOi+5AQCA7gn0AADQMYEeAAA6JtADAEDHBHoA\nAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEA\noGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDHBHoAAOiYQA8AAB0T6AEAoGNDHeir6vCquqmq\nbqmqk9ZyfNeq+mpVXVNV11XVCyajTgAA2FRDG+irakaS05I8P8neSY6uqr3XaPa3ST7VWtsvySuT\n/O+JrRIAADbP0Ab6JAuT3NJau7W19kCSc5IcuUabluQxg6+3T3L7BNYHAACbbZgD/c5Jfjhqe9lg\n32iLkry6qpYluSjJm9Z2oao6vqqWVNWS5cuXj0etAACwSYY50Nda9rU1to9OclZrbW6SFyT5WFU9\n6jFprZ3RWlvQWlswZ86ccSgVAAA2zTAH+mVJdhm1PTePXlLzuiSfSpLW2jeTbJ1k9oRUBwAAY2CY\nA/3iJLtX1W5VtWVGbnq9cI02P0jynCSpqr0yEuitqQEAoBtDG+hbaw8lOSHJxUluzMin2Xy7qt5Z\nVUcMmv1lkj+tqm8l+WSS17bW1lyWAwAAU9bMyS5gPLXWLsrIza6j95086usbkhww0XUBAMBYGdoZ\negAAmA4EegAA6JhADwAAHRPoAQCgYwI9AAB0TKAHAICOCfQAANAxgR4AADom0AMAQMcEegAA6JhA\nDwAAHRPoAQCgYwI9AAB0TKAHAICOCfQAANAxgR4AADom0AMAQMcEegAA6JhADwAAHRPoAQCgYwI9\nAAB0TKAHAICOCfQAANCxoQ70VXV4Vd1UVbdU1UnraPPyqrqhqr5dVZ+Y6BoBAGBzzJzsAsZLVc1I\nclqS5yVZlmRxVV3YWrthVJvdk7w9yQGttXuq6vGTUy0AAGyaYZ6hX5jkltbara21B5Kck+TINdr8\naZLTWmv3JElr7c4JrhEAADbLMAf6nZP8cNT2ssG+0fZIskdVfb2q/r2qDl/bharq+KpaUlVLli9f\nPk7lAgDAxhvmQF9r2dfW2J6ZZPckByc5OsmHqmqHR53U2hmttQWttQVz5swZ80IBAGBTDXOgX5Zk\nl1Hbc5PcvpY2n2utPdha+16SmzIS8AEAoAvDHOgXJ9m9qnarqi2TvDLJhWu0uSDJHyRJVc3OyBKc\nWye0SgAA2AxDG+hbaw8lOSHJxUluTPKp1tq3q+qdVXXEoNnFSVZU1Q1JvprkxNbaismpGAAANt7Q\nfmxlkrTWLkpy0Rr7Th71dUvyF4M/AADQnaGdoQcAgOlAoAcAgI4N9ZKb6W7e/Z/YoHZLx7cMAADG\nkRl6AADomEAPAAAdE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAd\nE+gBAKBjAj0AAHRMoAcAgI4J9AAA0DGBHgAAOibQAwBAxwR6AADomEAPAAAdE+gBAKBjQx3oq+rw\nqrqpqm6pqpN+TbuXVVWrqgUTWR8AAGyuoQ30VTUjyWlJnp9k7yRHV9Xea2m3XZI/T3LlxFYIAACb\nb2gDfZKFSW5prd3aWnsgyTlJjlxLu/83ybuT3D+RxQEAwFgY5kC/c5IfjtpeNti3WlXtl2SX1trn\nJ7IwAAAYK8Mc6Gst+9rqg1VbJPmHJH+53gtVHV9VS6pqyfLly8ewRAAA2DzDHOiXJdll1PbcJLeP\n2t4uyVOuAhy5AAAIS0lEQVSTXFZVS5M8K8mFa7sxtrV2RmttQWttwZw5c8axZAAA2DjDHOgXJ9m9\nqnarqi2TvDLJhQ8fbK3d21qb3Vqb11qbl+TfkxzRWlsyOeUCAMDGG9pA31p7KMkJSS5OcmOST7XW\nvl1V76yqIya3OgAAGBszJ7uA8dRauyjJRWvsO3kdbQ+eiJoAAGAsDe0MPQAATAcCPQAAdEygBwCA\njgn0AADQMYEeAAA6JtADAEDHhvpjK9l48+7/xAa1Wzq+ZQAAsIHM0AMAQMcEegAA6JhADwAAHRPo\nAQCgYwI9AAB0TKAHAICOCfQAANAxgR4AADom0AMAQMcEegAA6JhADwAAHRPoAQCgYwI9AAB0TKAH\nAICOCfQAANAxgR4AADo21IG+qg6vqpuq6paqOmktx/+iqm6oquuq6itV9duTUScAAGyqmZNdwHip\nqhlJTkvyvCTLkiyuqgtbazeManZNkgWttV9U1Z8leXeSV0x8tf2ad/8nNqjd0vEtAwBg2hrmGfqF\nSW5prd3aWnsgyTlJjhzdoLX21dbaLwab/55k7gTXCAAAm2WYA/3OSX44anvZYN+6vC7JF9d2oKqO\nr6olVbVk+fLlY1giAABsnmEO9LWWfW2tDatenWRBkves7Xhr7YzW2oLW2oI5c+aMYYkAALB5hnYN\nfUZm5HcZtT03ye1rNqqq5yb5f5Ic1Fr75QTVBgAAY2KYZ+gXJ9m9qnarqi2TvDLJhaMbVNV+ST6Q\n5IjW2p2TUCMAAGyWoQ30rbWHkpyQ5OIkNyb5VGvt21X1zqo6YtDsPUm2TfLpqrq2qi5cx+UAAGBK\nGuYlN2mtXZTkojX2nTzq6+dOeFEAADCGhnaGHgAApoOhnqFnavLLqAAAxo4ZegAA6JhADwAAHRPo\nAQCgYwI9AAB0zE2xTHluogUAWDcz9AAA0DEz9AwdM/oAwHQi0EM2/k2ANw0AwFQh0MMEmYg3DfrQ\nR499TMWa9KEPffRV00T1MVVZQw8AAB0T6AEAoGMCPQAAdEygBwCAjgn0AADQMYEeAAA6JtADAEDH\nBHoAAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjg11oK+qw6vqpqq6papOWsvxrarq3MHxK6tq3sRX\nCQAAm25oA31VzUhyWpLnJ9k7ydFVtfcazV6X5J7W2pOT/EOS/zWxVQIAwOYZ2kCfZGGSW1prt7bW\nHkhyTpIj12hzZJKPDr4+L8lzqqomsEYAANgs1Vqb7BrGRVW9LMnhrbXjBtvHJHlma+2EUW3+c9Bm\n2WD7u4M2d61xreOTHD/Y3DPJTRPwLWyo2UnuWm8rhoXxnn6M+fRjzKcfYz69bMx4/3Zrbc76Gs3c\nvHqmtLXNtK/57mVD2qS1dkaSM8aiqLFWVUtaawsmuw4mhvGefoz59GPMpx9jPr2Mx3gP85KbZUl2\nGbU9N8nt62pTVTOTbJ/k7gmpDgAAxsAwB/rFSXavqt2qasskr0xy4RptLkzymsHXL0tyaRvWNUgA\nAAyloV1y01p7qKpOSHJxkhlJzmytfbuq3plkSWvtwiQfTvKxqrolIzPzr5y8ijfZlFwKxLgx3tOP\nMZ9+jPn0Y8ynlzEf76G9KRYAAKaDYV5yAwAAQ0+gBwCAjgn0naqqw6vqpqq6papOmux6GHtVdWZV\n3Tn4fQkP73tcVX25qv5r8PdjJ7NGxlZV7VJVX62qG6vq21X15sF+4z6Eqmrrqrqqqr41GO9TBvt3\nq6orB+N97uCDHRgiVTWjqq6pqs8Pto35EKuqpVV1fVVdW1VLBvvG9Oe6QN+hqpqR5LQkz0+yd5Kj\nq2rvya2KcXBWksPX2HdSkq+01nZP8pXBNsPjoSR/2VrbK8mzkrxx8G/buA+nXyY5pLX2tCT7Jjm8\nqp6V5H8l+YfBeN+T5HWTWCPj481Jbhy1bcyH3x+01vYd9fnzY/pzXaDv08Ikt7TWbm2tPZDknCRH\nTnJNjLHW2uV59O9FODLJRwdffzTJURNaFOOqtXZHa+0/Bl//NCMv+DvHuA+lNuJng81Zgz8tySFJ\nzhvsN95DpqrmJnlhkg8NtivGfDoa05/rAn2fdk7yw1Hbywb7GH47tdbuSEbCX5LHT3I9jJOqmpdk\nvyRXxrgPrcHSi2uT3Jnky0m+m+THrbWHBk38fB8+70vy10lWDbZ3jDEfdi3Jl6rq6qo6frBvTH+u\nD+3n0A+5Wss+nz8KQ6Kqtk3ymSRvaa39ZGQCj2HUWluZZN+q2iHJZ5PstbZmE1sV46WqXpTkztba\n1VV18MO719LUmA+XA1prt1fV45N8uaq+M9YdmKHv07Iku4zanpvk9kmqhYn1o6p6QpIM/r5zkuth\njFXVrIyE+bNba+cPdhv3Idda+3GSyzJy78QOVfXwhJuf78PlgCRHVNXSjCyXPSQjM/bGfIi11m4f\n/H1nRt64L8wY/1wX6Pu0OMnug7vit8zIb7i9cJJrYmJcmOQ1g69fk+Rzk1gLY2ywlvbDSW5srf1/\now4Z9yFUVXMGM/Opqm2SPDcj9018NcnLBs2M9xBprb29tTa3tTYvI6/dl7bWXhVjPrSq6jeraruH\nv05yaJL/zBj/XPebYjtVVS/IyLv6GUnObK39z0kuiTFWVZ9McnCS2Ul+lOQdSS5I8qkkuyb5QZI/\nbK2teeMsnaqqZyf5tyTX51fra/8mI+vojfuQqap9MnIz3IyMTLB9qrX2zqp6UkZmbx+X5Jokr26t\n/XLyKmU8DJbc/FVr7UXGfHgNxvazg82ZST7RWvufVbVjxvDnukAPAAAds+QGAAA6JtADAEDHBHoA\nAOiYQA8AAB0T6AEAoGMCPQAAdEygBwCAjv3/nsQ0olwuoTsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2c7fe052390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "outcomes = ['unemployed', 'employed']\n",
    "\n",
    "P = [\n",
    "    [0.55, 0.45],\n",
    "    [0.05, 0.95],\n",
    "]\n",
    "\n",
    "pi0 = [1, 0]\n",
    "\n",
    "\n",
    "worker = MarkovChain(P, pi0, outcomes)\n",
    "\n",
    "fig, (ax0,ax1) = plt.subplots(2, 1, sharex=True, figsize=[12,8])\n",
    "worker.plot(ax0, T, 0)\n",
    "worker.plot_probability(ax1,T)\n",
    "worker.print_stationary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Marital status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAHVCAYAAABfQyW0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2YXWV9L/zvjyRIBIEq8Vy8SMN5iigl8mICCEcKFgEF\nAZ+qiNpqK2J7HXw8bUVFW430UGlpe/RpPdaoPOixvKhVRKVCKVJsq5JEVFCKokYN8RJEQUAUEu/n\nj9ngGCZkMpm5Z8/M53NdudhrrXut9dv7XrPnu2/utadaawEAAPrZZroLAACAuUYIBwCAzoRwAADo\nTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRwAADobP50F9DDLrvs0hYvXjzdZQAAMMutXr36\nB621RZtrNydC+OLFi7Nq1arpLgMAgFmuqr49nnamowAAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBn\nQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4\nAAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAA\ndCaEAwBAZ0I4AAB0Nn9Ld6iq5UnuSbJjkmtba1dNdlFbWM8FST7RWvvwdNbB5Lr0+ltz3hU3Z92d\n92W3nRfmzGP3yckH7j5p7Se6D1OjR/8N6zl6mC2v1Vw9x0QM47U4jK+tc8z897eZbN7y5cu3aIe3\nvOUtRya5v7V2zvLly7+5tQVU1bzly5e3ie7/lre85eQkX1u+fPlXN9VmxYoVy08//fSJnoLOLr3+\n1pz1kRvyw5/cnyS5+6fr869fuz17/MrCPGnXHbe6/UT3YWr06L9hPUcPs+W1mqvnmIhhvBaH8bV1\njpn//jas3vKWt3xv+fLlKzbXrlrbfP6tqjcm+Z0k301ye5LVSfZL8okk9yb53dbaCwZtj0zyx621\n51TVqUnekKSSfLK19rpBm3uS/E2SY5P8cZKfJXl7ku0Hj38zyU+SnJvkyCSPSvKO1tq7qqqS/G2S\nZyT51uDY5z/SSPjSpUvbqlWrNvs8mToXfv47+dgXbx1X2+u/c2fu3/Dzh63fdt42OXDPnbe6/SPt\ns/vOC/Pvr3/GuOrkkY23zyez/3pcI1t6jqm4pmbKz5NzjP8cm7tOJqPPJ/tadB06h9+ZY6uq1a21\npZtrt9k54VX11CQvTHJgkv87ybKNmvxzkkOravvB8ilJLqmq3ZL8RUbC8gFJllXVyYM22ye5sbV2\nSJLrklyS5NWttf2THJ3kviQvT3JXa23Z4JyvqKq9kjw3yT5JliR5RZLDNlH36VW1qqpW3X777Zt7\nmkyxj33x1nz1ez8eV9uxftAnc/0jbVt3532bqY7xGm+fT2b/9bhGtnT9VFxTM+XnyTnGv21z18lk\n9PlkX4uuQ+fwO3PrjGdO+NOTfLS19pMkqarLRm9sra2vqk8leU5VfTjJ8Ulem5HwfU1r7fbBfv+Q\n5IgklybZkOQfB4fYJ8n3WmsrB8f78aD9MUmeUlXPG7TbKcneg2Nc1FrbkGRdVV09VtGttRVJViQj\nI+HjeJ5MsX133TGXvPJpm213+LlX59YxfrB333nhmPtvaftH2me3nRdutj7Gbzx9Ppn91+Ma2dJz\nTNU1NRN+npxj/OcYz3WytX0+Fdei63Bun8PvzK0z3m9H2VyIvSTJCzISvFe21u7OyDSRTfnpIERn\n0G6s41eSV7XWDhj826u1duU462EGO/PYfbJwwbxfWrdwwbyceew+k9J+ovswNXr037Ceo4fZ8lrN\n1XNMxDBei8P42jrHzH9/m+nGMxJ+bZILqurcQfvnJHnXRm2uSfLejEwPuWSw7vNJ3l5VuyT5UZJT\nMzKXe2P/mWS3qlrWWltZVY/JyHSUK5L8QVVd3Vp7oKqemOTWQT2vrKr3J3l8kqOSXDjeJ8zwe/Bu\n6/Hehb2l7Ufv89oPfzn3b/h5dnen97TZmv7rcY1s6TmG7Zoa5tdqLp9jKq+TYbwWZ1v/zdVzDNM1\nNRts6Y2Z306yNslXM7gx88EbIqvq75K8LMnjR01deVGSszIyqn15a+21g/X3tNZ2GHX8ZRkJ6Asz\nEsCPzsiNmf8zI6G/MnJD6MlJfpxf3Jj5tcEhPuDGzOF2yrs+myTj+t+WPQ1rXbPBXH1tezzvufra\nziZb2ocT6fOpvk5ch3OPPh+f8d6YOa7vCW+tnZPknM20OSPJGRutuzBjjFKPDuCD5ZVJDh3jsG8Y\n/NvYGWOsAwCAGcFfzAQAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwI\nBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcA\ngM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDO\nhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRw\nAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA\n6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6EcAAA6EwIBwCAzuZPdwEbq6r/aK0dtgXtlye5\np7X2V1NX1cxx6fW35rwrbs66O+/LbjsvzJnH7pOTD9x90tpP9BzXf+fO3L/h5zn83KvHdY5h1Ou1\nmi3nmA193oPXls3p0eeuQ6bCMP5+GibTFsKran5rbf2o5XmttQ1bEsD5ZZdef2vO+sgNue+BDUmS\nW++8L2d95IYkGfOC3NL2W3OO+zf8fNznGEY9X6vZco6Z3uc9eG3ZnB597jpkKgzj76dhU621Lduh\nanGSTyX5tySHJvlSkv8vyVuSPD7JiwdN35ZkYZL7kvxua+3mqnpZkuOTbJdk+yRnJ3lzku8lOaC1\ntm9V3dNa22FwrjOTvCDJo5J8tLX25sH6Nyb5nSTfTXJ7ktWPNBK+dOnStmrVqi16nsPiws9/Jx/7\n4q3javvgqMTGtp23TQ7cc+etbj+Z59h954X599c/Y8xz9HTKuz6br37vx9l31x0fsd10vlaz5RzD\n0udTbbzXVOK1ncu29r1nPH0+Ve9vrsO5azre3yby+2m6r8WqWt1aW7q5dhMdCf+1JM9PcnqSlUle\nlOS/JTkxyRsyEpCPaK2tr6qjk/x5kt8a7Pu0JE9prf2wqo5McnCS/Vpr39roCRyTZO/B9kpyWVUd\nkeTeJC9McuCg/i8kWb1xgVV1+qC+7LnnnhN8mtPvY1+8ddwX/FgX4mSun8xjrbvzvk2eo6eTDhjf\nJ+XpfK1myzmGpc+n2nivqcRrO5dt7XvPePp8qt7fXIdz13S8v03k99NMuRYnGsK/1Vq7IUmq6itJ\n/qW11qrqhiSLk+yU5H1VtXeSlmTBqH3/ubX2w1HL120cwAeOGfy7frC8Q0ZC+WMyMir+k8H5Lxur\nwNbaiiQrkpGR8Ak9yyGx76475pJXPm2z7Q4/9+rcOsaFt/vOC8fcf0vbT+Y5dtt54ZjH7+1Fh+yZ\nFx2y+Q9p0/lazZZzDEufT7XxXlOJ13Yu29r3nvH0+VS9v7kO567peH+byO+nmXItTvTbUX426vHP\nRy3/PCPB/s+SfLq1tl+S52Rk+smD7t3oWBsvP6iSvLW1dsDg36+11t472DajQ/VUOfPYfbJwwbxf\nWrdwwbyceew+k9K+1zmG0bC+VrPlHHOV15bN6dHnrkOmwjD+fho2U3Vj5k5JHpzI/LIJHuOKJH9W\nVf/QWrunqnZP8kCSa5NcUFXnZqT+5yR511bWOys8eBPCeO8S3tL2vc4xjIb1tZot55irvLZsTo8+\ndx0yFYbx99OwmeiNmZ8YjHKnqi4YLH/4wW1JXpHkfRm5afLqJL/dWls8uDFzaWvtjMG+RyZ5TWvt\nhFHHH31j5quTnDbYdE+Sl7TWvjHqxsxvJ1mb5Kuz9cbMU9712SQZ13QUAACm13hvzNziED4TCeEA\nAPQw3hDuL2YCAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdC\nOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgA\nAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0\nJoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaE\nAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMA\nQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHTWLYRX1Xuqat8J7ru4qm6c7JoAAGA6zO91otba\nab3ONawuvf7WnHfFzVl3533ZbeeFOfPYfXLygbs/Yvvrv3Nn7t/w8xx+7tWbbQ8AwMwwJSPhVbV9\nVX2yqr5UVTdW1SlVdU1VLR1sv6eqzhls/1xV/ZfB+v9rsLyyqs6uqnvGOPa8qjpv0ObLVfXKqXgO\nk+3S62/NWR+5IbfeeV9aklvvvC9nfeSGXHr9rY/Y/v4NP0/G0R4AgJljqkbCj0uyrrV2fJJU1U5J\n/mDU9u2TfK619saq+sskr0jyP5O8PcnbW2sXVdXvb+LYL09yV2ttWVU9Ksm/V9WVrbVvTdFzmRTn\nXXFz7ntgwy+tu++BDXnth7+ci677zsPaPzgCvnH786642Wg4AMAMN1Vzwm9IcnRV/UVVPb21dtdG\n2+9P8onB49VJFg8ePy3JhwaPL9zEsY9J8jtV9cUkn0/yuCR7b9yoqk6vqlVVter222+f+DOZJOvu\nvG/M9RsH7c2t39RxAACYOaZkJLy19rWqemqSZyd5a1VduVGTB1prbfB4wxbWUUle1Vq7YjM1rEiy\nIkmWLl3aHqltD7vtvDC3jhGgd995YS555dMetv7wc68es/1uOy+ckvoAAOhnquaE75bkJ621DyT5\nqyQHjXPXzyX5rcHjF26izRVJ/qCqFgzO9cSq2n5r6u3hzGP3ycIF835p3cIF83LmsftMSnsAAGaO\nqZoTviTJeVX18yQPZGQ++F+NY7//keQDVfXHST6ZZONpLEnynoxMX/lCVVWS25OcPBlFT6UH53GP\n99tRtrQ9AAAzR/1iVsj0q6pHJ7mvtdaq6oVJTm2tnbS1x126dGlbtWrV1hcIAACPoKpWt9aWbq5d\nt+8JH6enJvm7wQj3nUl+b5rrAQCASTdUIby19pkk+093HQAAMJW6/dl6AABghBAOAACdCeEAANCZ\nEA4AAJ0J4QAA0JkQDgAAnQnhAADQmRAOAACdCeEAANCZEA4AAJ0J4QAA0JkQDgAAnQnhAADQmRAO\nAACdCeEAANCZEA4AAJ0J4QAA0JkQDgAAnQnhAADQmRAOAACdCeEAANCZEA4AAJ0J4QAA0JkQDgAA\nnQnhAADQmRAOAACdVWttumuYclV1e5JvT3cdo+yS5AfTXQRd6fO5RX/PPfp87tHnc894+/xXW2uL\nNtdoToTwYVNVq1prS6e7DvrR53OL/p579Pnco8/nnsnuc9NRAACgMyEcAAA6E8Knx4rpLoDu9Pnc\nor/nHn0+9+jzuWdS+9yccAAA6MxIOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAA\ndCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQm\nhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQD\nAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0Nn+6C+hhl112aYsXL57uMgAAmOVWr179g9ba\nos21mxMhfPHixVm1atV0lwEAwCxXVd8eTzvTUQAAoDMhHAAAOhPCAQCgs6GbE15V5yc5IcltrbX9\nxtj+4iSvGyzek+QPWmtf6lgiAEzYAw88kLVr1+anP/3pdJcCbIXtttsue+yxRxYsWDCh/YcuhCe5\nIMnfJXn/JrZ/K8lvtNZ+VFXPSrIiySGdagOArbJ27do85jGPyeLFi1NV010OMAGttdxxxx1Zu3Zt\n9tprrwkdY+imo7TWrk3yw0fY/h+ttR8NFj+XZI8uhQHAJPjpT3+axz3ucQI4zGBVlcc97nFb9X+0\nhi6Eb6GXJ/mnsTZU1elVtaqqVt1+++2dywKATRPAYebb2p/jGRvCq+qojITw1421vbW2orW2tLW2\ndNGizX5fOgAAdDOMc8I3q6qekuQ9SZ7VWrtjuusBgIla/PpPTurx1px7/Bbvc9ppp+WP/uiPsu++\n+275+dasyQknnJAbb7xxi/edsOU7TfLx7prc403AYYcdlv/4j/8Yd/vly5dnhx12yGte85pJr2XJ\n+5ZM6vFueOkNW9T+wef24x//OEcccUSOPvroSa1nS73sZS/LCSeckOc973mTetwZF8Kras8kH0ny\n2621r013PQAw073nPe+Z7hLmjPXr12f+/F/Erw0bNmTevHlbFMDnirPPPntSjvPgazxshm46SlVd\nlOSzSfapqrVV9fKq+v2q+v1BkzcleVyS/11VX6wqf48eAMbp3nvvzfHHH5/9998/++23Xy655JIc\neeSRWbVq5NfpDjvskDe+8Y3Zf//9c+ihh+b73/9+kuQb3/hGDj300CxbtixvetObssMOOzzs2Bs2\nbMiZZ56ZZcuW5SlPeUre9a53dX1uU2nNmjV50pOelNNOOy377bdfXvziF+eqq67K4Ycfnr333jvX\nXXddrrvuuhx22GE58MADc9hhh+Xmm29OklxwwQV5/vOfn+c85zk55phjcs011+Soo47Ki170oixZ\nMjLqPPr1PO+88x56Dd/85jc/tP6cc87JPvvsk6OPPvqhY88WYz23l73sZfnwhz+cf/qnf8oLXvCC\nh9pec801ec5znpMkueiii7JkyZLst99+ed3rfjFDeYcddsib3vSmHHLIIfnsZz+blStX5rDDDsv+\n+++fgw8+OHffffcmr9fWWs4444zsu+++Of7443PbbbdNyXMeupHw1tqpm9l+WpLTOpUDALPKpz71\nqey222755CdHpsHcddddeec73/nQ9nvvvTeHHnpozjnnnLz2ta/Nu9/97vzJn/xJXv3qV+fVr351\nTj311Pz93//9mMd+73vfm5122ikrV67Mz372sxx++OE55phjJvwVbsPmlltuyYc+9KGsWLEiy5Yt\ny4UXXph/+7d/y2WXXZY///M/z/vf//5ce+21mT9/fq666qq84Q1vyD/+4z8mST772c/my1/+ch77\n2MfmmmuuyXXXXZcbb7zxYa/NlVdema9//eu57rrr0lrLiSeemGuvvTbbb799Lr744lx//fVZv359\nDjrooDz1qU+djpdh0q1evfoRn9szn/nMvPKVr8y9996b7bffPpdccklOOeWUrFu3Lq973euyevXq\n/Mqv/EqOOeaYXHrppTn55JNz7733Zr/99svZZ5+d+++/P0960pNyySWXZNmyZfnxj3+chQsXbvJ6\nvf7663PzzTfnhhtuyPe///3su++++b3f+71Jf95DNxIOAEydJUuW5KqrrsrrXve6fOYzn8lOO/3y\n/Optt902J5xwQpLkqU99atasWZNkJEQ+//nPT5K86EUvGvPYV155Zd7//vfngAMOyCGHHJI77rgj\nX//616fuyXS21157ZcmSJdlmm23y67/+6/nN3/zNVFWWLFmSNWvW5K677srzn//87LfffvnDP/zD\nfOUrX3lo32c+85l57GMf+9DywQcfPOaHkyuvvDJXXnllDjzwwBx00EH5z//8z3z961/PZz7zmTz3\nuc/Nox/96Oy444458cQTuzznHjb33ObPn5/jjjsuH//4x7N+/fp88pOfzEknnZSVK1fmyCOPzKJF\nizJ//vy8+MUvzrXXXpskmTdvXn7rt34rSXLzzTdn1113zbJly5IkO+64Y+bPn7/J6/Xaa6/Nqaee\nmnnz5mW33XbLM57xjCl53kM3Eg4ATJ0nPvGJWb16dS6//PKcddZZOeaYY35p+4IFCx766rV58+Zl\n/fr14z52ay1/+7d/m2OPPXZSax4Wj3rUox56vM022zy0vM0222T9+vX50z/90xx11FH56Ec/mjVr\n1uTII498qP3222//S8faePlBrbWcddZZeeUrX/lL69/2trfN6q+23NxzO+WUU/KOd7wjj33sY7Ns\n2bI85jGPSWttk+232267h+aBt9bGPP6mrtfLL7+8y2ttJBwA5pB169bl0Y9+dF7ykpfkNa95Tb7w\nhS+Ma79DDz30oakVF1988Zhtjj322Lzzne/MAw88kCT52te+lnvvvXdyCp8B7rrrruy+++5JRuaB\nT8Sxxx6b888/P/fcc0+S5NZbb81tt92WI444Ih/96Edz33335e67787HP/7xySp72o3nuR155JH5\nwhe+kHe/+9055ZRTkiSHHHJI/vVf/zU/+MEPsmHDhlx00UX5jd/4jYft+6QnPSnr1q3LypUrkyR3\n33131q9fv8nr9YgjjsjFF1+cDRs25Hvf+14+/elPT8nzNhIOANNoIl8puDVuuOGGnHnmmdlmm22y\nYMGCvPOd7xzX19y97W1vy0te8pL89V//dY4//viHTWNJRr7qcM2aNTnooIPSWsuiRYty6aWXTv6T\nGIKvFBzLa1/72rz0pS/N3/zN30x4CsMxxxyTm266KU972tOSjNxg+IEPfCAHHXRQTjnllBxwwAH5\n1V/91Tz96U+fzNJ/yZZ+peDWGs9zmzdvXk444YRccMEFed/73pck2XXXXfPWt741Rx11VFprefaz\nn52TTjrpYftuu+22ueSSS/KqV70q9913XxYuXJirrrpqk9frc5/73Fx99dVZsmRJnvjEJ44Z7CdD\nPdJQ/myxdOnS9uBd3wAwnW666aY8+clPnu4ytthPfvKTLFy4MFWViy++OBdddFE+9rGPTXdZMK3G\n+nmuqtWttaWb29dIOACwWatXr84ZZ5yR1lp23nnnnH/++dNdEsxoQjgAsFlPf/rT86UvfWm6y4BZ\nw42ZANDZXJgKCrPd1v4cC+EA0NF2222XO+64QxCHGay1ljvuuCPbbbfdhI9hOgoAdLTHHntk7dq1\nuf3226e7FGArbLfddtljjz0mvL8QDgAdLViwYNb8GXdg4kxHAQCAzoRwAADoTAgHAIDOhHAAAOhM\nCAcAgM6EcAAA6EwIBwCAzoRwAADoTAgHAIDOhHAAAOhMCAcAgM6GLoRX1flVdVtV3biJ7VVV/29V\n3VJVX66qg3rXCAAAW2PoQniSC5Ic9wjbn5Vk78G/05O8s0NNAAAwaYYuhLfWrk3yw0doclKS97cR\nn0uyc1Xt2qc6AADYekMXwsdh9yTfHbW8drAOAABmhPnTXcAE1Bjr2sMaVZ2ekekq2XPPPae6pk1a\n8r4l42p3w0tvmFB753COmXqOYazJOZzDOWZWTc7hHOPdZxjNxJHwtUmeMGp5jyTrNm7UWlvRWlva\nWlu6aNGibsUBAMDmzMQQflmS3xl8S8qhSe5qrX1vuosCAIDxGrrpKFV1UZIjk+xSVWuTvDnJgiRp\nrf19ksuTPDvJLUl+kuR3p6dSAACYmKEL4a21UzezvSX5753KAQCASTcTp6MAAMCMJoQDAEBnQjgA\nAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0\nJoQDAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0JoQDAEBnQjgAAHQmhAMAQGdCOAAAdDZ0\nIbyqjquqm6vqlqp6/Rjb96yqT1fV9VX15ap69nTUCQAAEzVUIbyq5iV5R5JnJdk3yalVte9Gzf4k\nyQdbawcmeWGS/923SgAA2DpDFcKTHJzkltbaN1tr9ye5OMlJG7VpSXYcPN4pybqO9QEAwFYbthC+\ne5LvjlpeO1g32vIkL6mqtUkuT/KqsQ5UVadX1aqqWnX77bdPRa0AADAhwxbCa4x1baPlU5Nc0Frb\nI8mzk/yfqnrY82itrWitLW2tLV20aNEUlAoAABMzbCF8bZInjFreIw+fbvLyJB9MktbaZ5Nsl2SX\nLtUBAMAkGLYQvjLJ3lW1V1Vtm5EbLy/bqM13kvxmklTVkzMSws03AQBgxhiqEN5aW5/kjCRXJLkp\nI9+C8pWqOruqThw0++Mkr6iqLyW5KMnLWmsbT1kBAIChNX+6C9hYa+3yjNxwOXrdm0Y9/mqSw3vX\nBQAAk2WoRsIBAGAuEMIBAKAzIRwAADoTwgEAoDMhHAAAOhPCAQCgMyEcAAA6E8IBAKAzIRwAADoT\nwgEAoDMhHAAAOhPCAQCgMyEcAAA6E8IBAKAzIRwAADqbP90FzHY3fOs7010CAABDxkg4AAB0JoQD\nAEBnQjgAAHQmhAMAQGdCOAAAdCaEAwBAZ0I4AAB0NnQhvKqOq6qbq+qWqnr9Jtq8oKq+WlVfqaoL\ne9cIAABbY6j+WE9VzUvyjiTPTLI2ycqquqy19tVRbfZOclaSw1trP6qqx09PtQAAMDHDNhJ+cJJb\nWmvfbK3dn+TiJCdt1OYVSd7RWvtRkrTWbutcIwAAbJWhGglPsnuS745aXpvkkI3aPDFJqurfk8xL\nsry19qmND1RVpyc5PUn23HPPKSl2Kvgz9wAAs9+wjYTXGOvaRsvzk+yd5MgkpyZ5T1Xt/LCdWlvR\nWlvaWlu6aNGiSS8UAAAmathC+NokTxi1vEeSdWO0+Vhr7YHW2reS3JyRUA4AADPCsIXwlUn2rqq9\nqmrbJC/Acl5KAAAKmUlEQVRMctlGbS5NclSSVNUuGZme8s2uVQIAwFYYqhDeWluf5IwkVyS5KckH\nW2tfqaqzq+rEQbMrktxRVV9N8ukkZ7bW7pieigEAYMsN242Zaa1dnuTyjda9adTjluSPBv8AAGDG\nGaqRcAAAmAuGbiScLedrDQEAZhYj4QAA0JkQDgAAnQnhAADQmTnhc5A55AAA08tIOAAAdCaEAwBA\nZ0I4AAB0Zk4442IeOQDA5DESDgAAnQnhAADQmekoTAnTVwAANs1IOAAAdCaEAwBAZ6ajMDRMYQEA\n5goj4QAA0JmRcGYsI+cAwExlJBwAADozEs6cYvQcABgGRsIBAKAzI+HwCIycAwBTYehCeFUdl+Tt\nSeYleU9r7dxNtHtekg8lWdZaW9WxRHhEgjsAsDlDNR2lquYleUeSZyXZN8mpVbXvGO0ek+T/SfL5\nvhUCAMDWG7aR8IOT3NJa+2aSVNXFSU5K8tWN2v1Zkr9M8pq+5cHkM3IOAHPPsIXw3ZN8d9Ty2iSH\njG5QVQcmeUJr7RNVJYQzJwnuADCzDVsIrzHWtYc2Vm2T5H8ledlmD1R1epLTk2TPPfecpPJgZhLa\nAWC4DNWc8IyMfD9h1PIeSdaNWn5Mkv2SXFNVa5IcmuSyqlq68YFaaytaa0tba0sXLVo0hSUDAMCW\nGbaR8JVJ9q6qvZLcmuSFSV704MbW2l1JdnlwuaquSfIa344Ck8/oOQBMnaEK4a219VV1RpIrMvIV\nhee31r5SVWcnWdVau2x6KwQ2RWgHgPEbqhCeJK21y5NcvtG6N22i7ZE9agImn9AOwFw2dCEcYFME\ndwBmCyEcmLUmEtoFfQB6EMIBtoLQDsBECOEAnW1pcBf0AWYfIRxgFhL0AYabED7FFv/0wnG1WzO1\nZQBMOkEfYOKEcACGVo+g78MBMB2EcADYAj2Cvg8GMPsJ4QAwC0x10PfhAyaXED5kzCEHgPGbLR8m\nnGPuEcIBAJgxZktw32a6CwAAgLlGCAcAgM5MR5kFzCMHAJhZjIQDAEBnQjgAAHQmhAMAQGfmhM9B\n5pADAEwvI+EAANCZEA4AAJ2ZjsK4mMICADB5jIQDAEBnQjgAAHQ2dNNRquq4JG9PMi/Je1pr5260\n/Y+SnJZkfZLbk/xea+3b3QvlEZm+AgCwaUM1El5V85K8I8mzkuyb5NSq2nejZtcnWdpae0qSDyf5\ny75VAgDA1hmqEJ7k4CS3tNa+2Vq7P8nFSU4a3aC19unW2k8Gi59LskfnGgEAYKsM23SU3ZN8d9Ty\n2iSHPEL7lyf5p7E2VNXpSU5Pkj333HOy6mMKmcICAMwVwzYSXmOsa2M2rHpJkqVJzhtre2ttRWtt\naWtt6aJFiyaxRAAA2DrDNhK+NskTRi3vkWTdxo2q6ugkb0zyG621n3WqjSFj5BwAmKmGbSR8ZZK9\nq2qvqto2yQuTXDa6QVUdmORdSU5srd02DTUCAMBWGaoQ3lpbn+SMJFckuSnJB1trX6mqs6vqxEGz\n85LskORDVfXFqrpsE4cDAIChNGzTUdJauzzJ5Rute9Oox0d3L4pZwxQWAGAYDNVIOAAAzAVDNxIO\nw8TIOQAwFYyEAwBAZ0bCYZIZPQcANsdIOAAAdGYkHKaZkXMAmHuEcJiBBHcAmNlMRwEAgM6MhMMc\nYOQcAIaLEA6MSXAHgKljOgoAAHRmJByYFEbOAWD8hHBgWgjtAMxlQjgwYwjuAMwWQjgwawntAAwr\nIRxgFMEdgB6EcICtILQDMBFCOEBnWxrcBX2A2UcIB5iFBHeA4SaEAzCh0C7oA0ycEA5AF4I+wC8I\n4QDMGj2Cvg8GwGQQwgFgivUI+j4cwMwydCG8qo5L8vYk85K8p7V27kbbH5Xk/UmemuSOJKe01tb0\nrhMAZrJh/GDggwRzyVCF8Kqal+QdSZ6ZZG2SlVV1WWvtq6OavTzJj1prv1ZVL0zyF0lO6V8tADDd\nZsuHiWE8B1NrqEJ4koOT3NJa+2aSVNXFSU5KMjqEn5Rk+eDxh5P8XVVVa631LBQAYDbzYWJq1TBl\n16p6XpLjWmunDZZ/O8khrbUzRrW5cdBm7WD5G4M2P9joWKcnOX2wuE+Smzs8hfHaJckPNtuK2USf\nzy36e+7R53OPPp97xtvnv9paW7S5RsM2El5jrNv4U8J42qS1tiLJiskoarJV1arW2tLproN+9Pnc\nor/nHn0+9+jzuWey+3ybyTrQJFmb5AmjlvdIsm5TbapqfpKdkvywS3UAADAJhi2Er0yyd1XtVVXb\nJnlhkss2anNZkpcOHj8vydXmgwMAMJMM1XSU1tr6qjojyRUZ+YrC81trX6mqs5Osaq1dluS9Sf5P\nVd2SkRHwF05fxRM2lNNkmFL6fG7R33OPPp979PncM6l9PlQ3ZgIAwFwwbNNRAABg1hPCAQCgMyG8\no6o6rqpurqpbqur1010Pk6+qzq+q2wbfZ//gusdW1T9X1dcH//2V6ayRyVVVT6iqT1fVTVX1lap6\n9WC9fp+lqmq7qrquqr406PO3DNbvVVWfH/T5JYMvGGCWqKp5VXV9VX1isKy/Z7GqWlNVN1TVF6tq\n1WDdpL6vC+GdVNW8JO9I8qwk+yY5tar2nd6qmAIXJDluo3WvT/IvrbW9k/zLYJnZY32SP26tPTnJ\noUn+++BnW7/PXj9L8ozW2v5JDkhyXFUdmuQvkvyvQZ//KMnLp7FGJt+rk9w0all/z35HtdYOGPXd\n4JP6vi6E93Nwkltaa99srd2f5OIkJ01zTUyy1tq1efj31p+U5H2Dx+9LcnLXophSrbXvtda+MHh8\nd0Z+Se8e/T5rtRH3DBYXDP61JM9I8uHBen0+i1TVHkmOT/KewXJFf89Fk/q+LoT3s3uS745aXjtY\nx+z3X1pr30tGAluSx09zPUyRqlqc5MAkn49+n9UGUxO+mOS2JP+c5BtJ7mytrR808R4/u7wtyWuT\n/Hyw/Ljo79muJbmyqlZX1emDdZP6vj5U3xM+y9UY63w/JMwSVbVDkn9M8j9aaz8eGShjtmqtbUhy\nQFXtnOSjSZ48VrO+VTEVquqEJLe11lZX1ZEPrh6jqf6eXQ5vra2rqscn+eeq+s/JPoGR8H7WJnnC\nqOU9kqybplro6/tVtWuSDP572zTXwySrqgUZCeD/0Fr7yGC1fp8DWmt3JrkmI/cD7FxVDw5ueY+f\nPQ5PcmJVrcnIVNJnZGRkXH/PYq21dYP/3paRD9oHZ5Lf14XwflYm2XtwN/W2GflLn5dNc030cVmS\nlw4evzTJx6axFibZYG7oe5Pc1Fr7m1Gb9PssVVWLBiPgqaqFSY7OyL0An07yvEEzfT5LtNbOaq3t\n0VpbnJHf3Ve31l4c/T1rVdX2VfWYBx8nOSbJjZnk93V/MbOjqnp2Rj49z0tyfmvtnGkuiUlWVRcl\nOTLJLkm+n+TNSS5N8sEkeyb5TpLnt9Y2vnmTGaqq/luSzyS5Ib+YL/qGjMwL1++zUFU9JSM3Zc3L\nyGDWB1trZ1fVf83ISOljk1yf5CWttZ9NX6VMtsF0lNe01k7Q37PXoG8/Olicn+TC1to5VfW4TOL7\nuhAOAACdmY4CAACdCeEAANCZEA4AAJ0J4QAA0JkQDgAAnQnhAADQmRAOAACd/f91Nu8i+eb6UAAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2c7fdf54eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "outcomes = ['single', 'married', 'divorced']\n",
    "\n",
    "P = [\n",
    "    [0.92, 0.08, 0.00],\n",
    "    [0.00, 0.65, 0.35],\n",
    "    [0.00, 0.15, 0.85]\n",
    "]\n",
    "\n",
    "civil = MarkovChain(P, pi0=None, z=outcomes)\n",
    "\n",
    "fig, (ax0,ax1) = plt.subplots(2, 1, sharex=True, figsize=[12,8])\n",
    "civil.plot(ax0, T, 0)\n",
    "civil.plot_probability(ax1,T)\n",
    "civil.print_stationary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Credit rating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AAA           0.28\n",
      "AA            0.88\n",
      "A             2.78\n",
      "BBB           4.06\n",
      "BB            5.75\n",
      "B            11.25\n",
      "CCC          12.93\n",
      "D            40.29\n",
      "N.R.         21.78\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAHVCAYAAAAHPLatAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXFWd7/3PjyaaBoQoBJWEGEAIMIBBWxmEKOA4REG5\nDEhgHBkPnggPnBEvOMYzo8FnRqLoiI48Z8zgBcZj4sgAchBBnwFGbmqCQSNELoEISbiESBAkSEh+\n54+uhE6nd7q6ui67qj7v1ysvutZee6+1axWVb6+svXdkJpIkSZK2tE2rOyBJkiSVlWFZkiRJKmBY\nliRJkgoYliVJkqQChmVJkiSpgGFZkiRJKmBYliRJkgoYliVJkqQChmVJkiSpwLat7sBAu+yyS06e\nPLnV3ZAkSVKHu+OOO57IzPHD1StVWJ48eTILFy5sdTckSZLU4SLit9XUcxmGJEmSVMCwLEmSJBUw\nLEuSJEkFDMuSJElSAcOyJEmSVMCwLEmSJBUY9tZxEZHAP2XmRyuvPwbskJmzB9U7Avg+8ADQC1yT\nmR+rd4el4Vy1aAUXXn8PK9esZbdxvZx39BSOP3hC27XRDJ1yHmXk57BcfK/KxfFQO6nmPst/BE6M\niAsy84lh6t6cmcdGRC+wKCKuzMxbR99NqTpXLVrBrCsWs3bdegBWrFnLrCsWA9Tti7gZbTRDp5xH\nGfk5LBffq3JxPNRuqgnLLwBzgQ8D/7Oag2bm2oi4E/BTr1H7zs8e4vt3rqiq7qKH1vD8+g2bla1d\nt54Lr79nq1/CzWijGTrlPMrIz2G5+F6Vi+OhTlbtmuWLgb+MiJ2qqRwRLwf2Bn5SRd2ZEbEwIhau\nWrWqyu6om3z/zhXc/cjvq6o7+At4o5Vr1ra8jWbolPMoIz+H5eJ7VS6OhzpZVY+7zszfR8RlwN8A\nW/s0T4uIXwFTgDmZ+WgVx55L/8w1fX19WU1/1H32f/WOfPeDhw5b77A5N7BiiC/c3cb1lqKNZuiU\n8ygjP4fl4ntVLo6HOtVI7oZxEXAGsP1W6tycmQcBBwJnRcTU0XROGqnzjp5C75iezcp6x/Rw3tFT\n2qqNZuiU8ygjP4fl4ntVLo6H2k3VYTkzfwf8O/2Bebi69wIXAH9be9ekkTv+4AlccOKBvKSn/6M9\nYVwvF5x4YF3XwTWjjWbolPMoIz+H5eJ7VS6Oh9pNVcswBvgicM7GFxHxbqAvMz81RN1/AT4WEXsA\nOwNnZuYHau6pVKXjD57AvJ8/BFDVPwmWtY1m6JTzKCM/h+Xie1UujofaybBhOTN3GPDzY8B2A15f\nDVxd+fkm4KYB29by4t0wHgQMypIkSWorPsFPkiRJKmBYliRJkgoYliVJkqQChmVJkiSpgGFZkiRJ\nKmBYliRJkgoYliVJkqQChmVJkiSpgGFZkiRJKmBYliRJkgoYliVJkqQChmVJkiSpgGFZkiRJKrBt\nIw4aEeuBxcAY4AXgUuCizNzQiPak0bhq0QoWPbSG59dv4LA5N3De0VM4/uAJdW/jwuvvYeWatew2\nrreqNka6TyedR9naaNZ76/iVazzK+F51q7J+DjuhDQ2vIWEZWJuZUwEiYlfgO8BOwKcb1J5Uk6sW\nrWDWFYt5fn3/73Er1qxl1hWLAer25bKxjbXr1lfdxkj36bTzKFMbzXxvHb9yjUfZ3qtuVebPYbu3\noepEZtb/oBHPZOYOA17vCSwAdsmtNNjX15cLFy6se3/U3k752u0AfPeDh45on7sf+T37v3rHrdbb\nODs12IRxvdz6iaMa2sZLerbh4Enj6rJPp59HK9to5Xvr+JVrPMpyHmXUiu/pMn4Oy9pGu36uGi0i\n7sjMvuHqNWpmeTOZ+UBEbAPsCjw2cFtEzARmAkyaNKkZ3VEXOG5qdb9BD/WlArByzdqGt1FUXss+\nnX4erWyjle+t41d9eae8V9WcR6do1Hh0+me9lja66XPVCE0JyxUxVGFmzgXmQv/MchP7ow522iGT\nOO2Q4X/5OmzODawY4ktkt3G9DW9jwrjewlmYke7T6efRyjZa+d46fuUaj7KcR6do1Hh0+me9lja6\n6XPVCE25G0ZlGcZ64PFmtCdV67yjp9A7pmezst4xPZx39JSWtjHSfbr5PBrdRlnf22a00QnjV4tO\nOY9O0SmfwzK2oer0zJ49u+4HPf/88z85e/bszwJExHjgm8AVmXnj1vabO3fu7JkzZ9a9P2pvl9+x\nHICT+3av+7H3ffWOTHx5L4tXPMUzz73AhHG9fOpd+9f1Qoha2hjpPt18Ho1uo6zvbTPa6ITxq0Uz\nz+OG3zzO+syGnEczlel7uqyfw2a20Smfq0Y7//zzH5k9e/bc4eo16gK/wbeO+zfgn4a7dZwX+Gko\ntVw4IkntoFO+3zrlPDqF41Gdll7gl5k9w9eSJEmSys0n+EmSJEkFDMuSJElSAcOyJEmSVMCwLEmS\nJBUwLEuSJEkFDMuSJElSAcOyJEmSVMCwLEmSJBUwLEuSJEkFDMuSJElSAcOyJEmSVMCwLEmSJBUw\nLEuSJEkFtq2mUkS8CrgIeCPwR2AZcG5l80XAPsA6YDHwPzLzsYh4E/AF4JVAArcAf5OZz9bzBCRJ\nklSbqxat4MLr72HlmrXsNq6X846ewvEHT2h1t0pl2LAcEQFcCVyamTMqZVPpD8HfAD6Smf+nUn4k\nML5/F74HzMjM2yvH+AvgZYBhWZIkqcWuWrSCWVcsZu269QCsWLOWWVcsBjAwD1DNzPKRwLrM/JeN\nBZl5Z0T8N+D2jUG5Un4jQER8hv5wfXulPIHL69pzSZIkDenuR37PKV+7fat1Fj20hufXb9isbO26\n9Vx4/T2G5QGqCcsHAHeMoHzjtkur6UBEzARmAkyaNKmaXSRJklTguKnVBd3BQXmjlWvW1rM7ba+q\nNcuNlJlzgbkAfX192eLuSJIktbXTDpnEaYcMPwF52JwbWDFEMN5tXG8jutW2qrkbxl3AG0ZQPtw2\nSZIktdh5R0+hd0zPZmW9Y3o47+gpLepROVUTlm8AXhoR/31jQUS8EbgfeHNEHDOgfHpEHAh8FTg9\nIg4ZsO29lbtqSJIkqcWOP3gCF5x4IC/p6Y+DE8b1csGJB7peeZBhl2FkZkbECcBFEfEJ4DlevHXc\nsZXyi+i/ddyvgA9Vbh03A/hCROwKbAB+AlzRmNOQJEnSSB1/8ATm/fwhAL77wUNb3JtyqmrNcmau\nBN5TsHl6wT63A9Nq7JckSZLUcj7BT5IkSSpgWJYkSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSpgWJYk\nSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSpgWJYkSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSqwbaMOHBHr\ngcVAAOuBczLztka1J0lSp7tq0QouvP4eVq5Zy27jejnv6Ckcf/CEutWvtY1FD63h+fUbOGzODVW1\nofbWjM9hmTQsLANrM3MqQEQcDVwAvLWB7UmS1LGuWrSCWVcsZu269QCsWLOWWVcsBhgyeIy0/mja\neH79hqrbUHtrxuewbCIzG3PgiGcyc4fKzycDf5mZx29tn76+vly4cGFD+qP2dcrXbgfgux88tMU9\nkaT6OuVrt3P3I79n/1fvOGzdjbO3g72kZxsOnjRu1PXr2caEcb3c+omjhmxD5dOKz2EZPiMRcUdm\n9g1Xr5Ezy70RcScwFng1MOQ7EhEzgZkAkyZNamB3JEkql+OmVj+zNlTgqGd5PY+1cs3awjZUPq34\nHLbTZ6RZyzAOBS6LiANy0FR2Zs4F5kL/zHID+yNJUqmcdsgkTjukuomiw+bcwIohAsaEcb1D/svb\nSOvXs43dxvUOeXyVUys+h+30GWnK3TAy83ZgF2B8M9qTJKnTnHf0FHrH9GxW1jumh/OOnlKX+s1q\nQ+2tGz8jjZxZ3iQi9gV6gNXNaE+SpE6z8WKoau8qMNL6zWpD7a0bPyONvMBv463joP/2cZ/MzB9s\nbR8v8NNQvMBPkiTVW8sv8MvMnuFrSZIkSeXlE/wkSZKkAoZlSZIkqYBhWZIkSSpgWJYkSZIKGJYl\nSZKkAoZlSZIkqYBhWZIkSSpgWJYkSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSpgWJYkSZIKGJYlSZKk\nAqMKyxGxPiLujIhfRsQvIuLNlfLJEbF2wLbbImJKfbosSZIkNcdoZ5bXZubUzHwdMAu4YMC2pQO2\nXQp8cpRtSZIkSU1Vz2UYOwJP1rBNkiRJKqVtR7l/b0TcCYwFXg0cNWDbXpVtLwO2Aw4Z6gARMROY\nCTBp0qRRdkeSJEmqn3otw9gXmA5cFhFR2bZxGcZewLnA3KEOkJlzM7MvM/vGjx8/yu5IkiRJ9VO3\nZRiZeTuwCzBU4r0aeEu92pIkSZKaoW5hOSL2BXqA1UNsPhxYWq+2JEmSpGao15plgABOz8z1lZUY\nG9csB/A88IFRtiVJkiQ11ajCcmb2FJQvA3pHc2xJkiSp1XyCnyRJklTAsCxJkiQVMCxLkiRJBQzL\nkiRJUgHDsiRJklTAsCxJkiQVMCxLkiRJBQzLkiRJUgHDsiRJklTAsCxJkiQVMCxLkiRJBQzLkiRJ\nUoFhw3JErI+IOyPilxHxi4h4c6V8ckSsHbDttoiYUtl2REQ8Vdn2q4j4/yNi10afjCRJklRP1cws\nr83MqZn5OmAWcMGAbUsHbLsU+OSAbTdXth0ELADOrluvJUmSpCYY6TKMHYEnR7ItIgJ42Vb2kyRJ\nkkpp2yrq9EbEncBY4NXAUQO27VXZ9jJgO+CQAdumVbbtDPyBzWedJUmSpNIbyTKMfYHpwGWV2WJ4\ncRnGXsC5wNwB+21chrE78E3g80MdPCJmRsTCiFi4atWqUZyKJEmSVF8jWoaRmbcDuwDjh9h8NfCW\ngl0Lt2Xm3Mzsy8y+8eOHOqwkSZLUGiMKyxGxL9ADrB5i8+HA0oJdt7ZNkiRJKqWRrFkGCOD0zFxf\nWYmxcc1yAM8DHxiw37QB254atE2SJEkqvWHDcmb2FJQvA3oLtt0E7DSajkmSJEmt5hP8JEmSpAKG\nZUmSJKmAYVmSJEkqYFiWJEmSChiWJUmSpAKGZUmSJKmAYVmSJEkqYFiWJEmSChiWJUmSpAKGZUmS\nJKmAYVmSJEkqYFiWJEmSChiWJUmSpAINDcsRcUJEZETs28h2JEmSpEZo9MzyqcAtwIwGtyNJkiTV\nXcPCckTsABwGnIFhWZIkSW2okTPLxwPXZea9wO8i4vUNbEuSJEmqu0aG5VOB+ZWf51debyEiZkbE\nwohYuGrVqgZ2R5IkSRqZbRtx0IjYGTgKOCAiEugBMiI+npk5sG5mzgXmAvT19eUWB5MkSZJapFEz\nyycBl2XmazJzcmbuDjwIHN6g9iRJkqS6a1RYPhW4clDZfwCnNag9SZIkqe4asgwjM48YouwrjWhL\nkiRJahSf4CdJkiQVMCxLkiRJBQzLkiRJUgHDsiRJklTAsCxJkiQVMCxLkiRJBQzLkiRJUgHDsiRJ\nklTAsCxJkiQVMCxLkiRJBQzLkiRJUgHDsiRJklTAsCxJkiQVGHVYjogTIiIjYt9B5R+OiOciYqfR\ntiFJkiS1Qj1mlk8FbgFmDFG+ADihDm1IkiRJTTeqsBwROwCHAWcwICxHxF7ADsDf0R+aJUmSpLYz\n2pnl44HrMvNe4HcR8fpK+anAPOBmYEpE7DrKdiRJkqSmG21YPhWYX/l5Pi/OIs8A5mfmBuAK4OSi\nA0TEzIhYGBELV61aNcruSJIkSfWzba07RsTOwFHAARGRQA+QEfFtYG/gxxEB8BLgAeDioY6TmXOB\nuQB9fX1Za38kSZKkehvNzPJJwGWZ+ZrMnJyZuwMPAhcBsytlkzNzN2BCRLymHh2WJEmSmmU0YflU\n4MpBZf8BTB6i/Eq2vFuGJEmSVGo1L8PIzCOGKPsK8JUhyj9SazuSJElSq/gEP0mSJKmAYVmSJEkq\nYFiWJEmSChiWJUmSpAKGZUmSJKmAYVmSJEkqYFiWJEmSChiWJUmSpAKGZUmSJKmAYVmSJEkqYFiW\nJEmSChiWJUmSpAKGZUmSJKlA1WE5Ik6IiIyIfQeVfzginouInYbY58sRsSIiDOWSJElqO9uOoO6p\nwC3ADGD2oPIFwAnAtzYWVgLyCcDDwFuAm0bV0ya6atEKLrz+HlauWctu43o57+gpHH/whJbV7/Y2\nFj20hufXb+CwOTdU1YYkSVK9VBWWI2IH4DDgSOBqKmE5IvYCdgDOAz7JgLBcqftr4Lv0B+qb6tPl\nxrpq0QpmXbGYtevWA7BizVpmXbEYYMiQ1uj6trGY59dvqLoNSZKkeorMHL5SxHuBIzPzjIi4DTgn\nM38REX8HBPCPwAPAmzLz8co+lwD/BXwfWAJMzsx1W2unr68vFy5cOKoTGq3D5tzAijVrtyh/Sc82\nHDxp3BblG2c9G1XfNrasP2FcL7d+4qgh25AkSapGRNyRmX3D1at2LfGpwPzKz/Mrr6F/Scb8zNwA\nXAGcXGn8JcA7gasy8/fAz4A/L+jozIhYGBELV61aVWV3GmflEEEZGDK0NaPcNrZUNEaSJEn1Nuwy\njIjYGTgKOCAiEugBMiK+DewN/DgiAF5C/+zyxcB0YCdgcWXbdsCzwA8GHz8z5wJzoX9mefSnNDq7\njesdcmZ5wrhevvvBQ7coL5qJrld929iy/m7jeoc8viRJUr1VM7N8EnBZZr4mMydn5u7Ag8BFwOxK\n2eTM3A2YEBGvoX/m+QMbtwF7AH8eEds16Dzq5ryjp9A7pmezst4xPZx39JSW1LeNkbUhSZJUT9Vc\n4HcqMGdQ2X8AHwauHFR+JfDXwNHABzcWZuYfIuIW4F30X/BXWhsvHKv2jg2Nrm8bI2tDkiSpnqq6\nwK9ZynCBnyRJkjpfvS/wkyRJkrqOYVmSJEkqYFiWJEmSChiWJUmSpAKGZUmSJKmAYVmSJEkqUKpb\nx0XEKuC3re7HALsAT7S6E2oax7v7OObdxzHvLo539xnJmL8mM8cPV6lUYblsImJhNfffU2dwvLuP\nY959HPPu4nh3n0aMucswJEmSpAKGZUmSJKmAYXnr5ra6A2oqx7v7OObdxzHvLo5396n7mLtmWZIk\nSSrgzLIkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIB\nw7IkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7Ik\nSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJU\nwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJUwLAsSZIkFTAsS5IkSQUMy5IkSVIBw7IkSZJUYNtW\nd2CgXXbZJSdPntzqbkiSJKnD3XHHHU9k5vjh6pUqLE+ePJmFCxe2uhuSJEnqcBHx22rquQxDkiRJ\nKmBYliRJkgoYliVJkqQCNa1ZjohvAMcCj2fmAUNs/0vgbysvnwHOysxf1txLtY1169axfPlynnvu\nuVZ3RZIkibFjxzJx4kTGjBlT0/61XuD3LeCrwGUF2x8E3pqZT0bEO4C5wCE1tqU2snz5cl72spcx\nefJkIqLV3ZEkSV0sM1m9ejXLly9njz32qOkYNS3DyMyfAL/byvbbMvPJysufAhNraUft57nnnmPn\nnXc2KEuSpJaLCHbeeedR/Yt3M9YsnwH8sGhjRMyMiIURsXDVqlVN6I4azaAsSZLKYrS5pKFhOSKO\npD8s/21Rncycm5l9mdk3fvyw94WWJEmSmqZhDyWJiIOAS4B3ZObqRrWjcpv8iR/U9XjL5hxTVb0r\nr7ySE088kSVLlrDvvvtuKv/Sl77ErFmzeOyxx9hpp5022+dDH/oQl19+OQ8//DDbbNPA3yNn7zR8\nnREd76mqqtXyntTLgZceWNfjLT59cVX1is65npbsu19dj7ffb5YMW6enp4cDDzyQzKSnp4evfvWr\nvPnNb2bZsmXst99+TJkyhcxk++2355vf/CZTpkzhpptu4rjjjmOPPfZgw4YN7LrrrnznO99h1113\nrWv/AS4+84a6Hu/sfzlq2Dq1vCf19sVTjq3r8T763WuGrVN03vW0/BM31/V4E+dMq6reo48+yrnn\nnsuCBQt46UtfyuTJk7nooosAOPfcc7n33nsZM2YMBx54IP/8z//MK1/5Sn7+85/zsY99jMcee4yI\n4PDDD+crX/kK2223XV3PYfbs2U0/3saxXrduHdtuuy2nn3465557bkP+vvrPG/aq6/HedtTSYetE\nBB/5yEf44he/CMAXvvAFnnnmmS3em43fZXvuuSdr167l2GOP5Qtf+EJd+ztYQxJBREwCrgD+KjPv\nbUQb0tbMmzePww8/nPnz529R/sY3vpErr7xys/INGzZw5ZVXsvvuu/OTn/ykmV1tmpG+J52g6Jzb\nXW9vL3feeSe//OUvueCCC5g1a9ambXvttdembaeffjqf/exnN22bNm0ad955J7/61a944xvfyMUX\nX9yK7jdEre9Ju9vaebezzOSEE07giCOOYOnSpdx999189rOf5bHHHuOYY47hrLPO4v7772fJkiWc\nddZZrFq1iscee4yTTz6Zz33uc9xzzz0sWbKE6dOn8/TTT7f6dOpi41jfdddd/PjHP+baa6/l/PPP\nb3W36ualL30pV1xxBU888cSwdadNm8aiRYtYtGgR11xzDbfeemtD+1ZTWI6IecDtwJSIWB4RZ0TE\nmRFxZqXKp4Cdgf8vIu6MCJ9hraZ55plnuPXWW/n617++WUhaunQpzzzzDP/wD//AvHnzNtvnxhtv\n5IADDuCss87aYlsnqOU9aXdF59xpfv/73/Pyl798RNsyk6effrpwv3ZXy3vSCTrp3G688UbGjBnD\nmWeeuals6tSp3HfffRx66KG8613v2lR+5JFHcsABB3DxxRdz+umnc+ihhwL9M5UnnXQSr3zlK5ve\n/0bbddddmTt3Ll/96lfJzFZ3py623XZbZs6cyZe+9KWq9+nt7WXq1KmsWLGigT2rcRlGZp46zPYP\nAB+oqUfSKF111VVMnz6dffbZh1e84hX84he/4PWvfz3z5s3j1FNPZdq0adxzzz08/vjjm/4JeuO2\n4447jk9+8pOsW7eu5vsxllEt70m7KzrnTrB27VqmTp3Kc889xyOPPMINN7y47GHp0qVMnTqVp59+\nmmeffZaf/exnm7bdfPPNTJ06ldWrV7P99tt31Axrre9Ju9vaebezX//617zhDW+ounzjttNPP73R\nXSuNPffckw0bNvD44493zC8EZ599NgcddBAf//jHq6r/5JNPct999/GWt7ylof3yCX7qOPPmzWPG\njBkAzJgxY9OM6fz585kxYwbbbLMNJ554It/73vcAeP7557n22ms5/vjj2XHHHTnkkEP40Y9+1LL+\nN8JI35NOUHTOnWDjP8f+5je/4brrruN973vfptmljUsOli5dykUXXcTMmTM37bdxGcbDDz/M+9//\n/qr/QmoHtb4n7W5r563O12ljveOOO/K+972Pr3zlK1utd/PNN3PQQQfxqle9imOPPZZXvepVDe1X\nwy7wk1ph9erV3HDDDfz6178mIli/fj0RwXvf+17uu+8+3v72twP9AXnPPffk7LPP5rrrruOpp57i\nwAP7L0J79tln2W677TjmmOouJiy7Wt6Tdld0zp///Oc77taGhx56KE888QRD3Xrz3e9+N+9///uH\n3O/d7343f/EXf9Ho7rVEre9Juxt43u3+L0R/8id/wuWXXz5k+X/9138V7nPHHXdw3HHHNbp7pfDA\nAw/Q09PT9mM92LnnnsvrX//6rf5/Om3aNK655hruvfdeDj/8cE444QSmTp3asD45s6yOcvnll/O+\n972P3/72tyxbtoyHH36YPfbYg3PPPZfZs2ezbNkyli1bxsqVK1mxYgW//e1vmTdvHpdccsmmbQ8+\n+CA/+tGPePbZZ1t9OnVRy3vS7orO+ZZbbml11+ruN7/5DevXr2fnnXfeYtstt9zCXnsNfVX71ra1\nu1rfk3a3tfNuN0cddRR//OMf+dd//ddNZQsWLOC1r30tt912Gz/4wYt3WrruuutYvHgx55xzDpde\neulmy2y+/e1v8+ijjza1782watUqzjzzTM4555yOmwB4xStewXve8x6+/vWvD1t3n332YdasWXzu\nc59rbKcyszR/3vCGN6Ta2913393S9t/61rfmD3/4w83KvvzlL+fkyZNzyZIlm5V/+MMfztmzZ+fL\nX/7yfOqppzbbdsIJJ+T8+fMb3t9mGOl7MmfOnGZ2ryGKzvnMM89sUY/qa5tttsnXve51+brXvS4P\nOuigvOaaazIz88EHH8yxY8duKu/r68uf/vSnmZl544035o477rhp27Rp0/Kee+5p5WnUVS3vSSco\nOu9OsGLFijz55JNzzz33zP333z/f+c535r333ptLlizJo48+Ol/72tfmfvvtl6eccko++uijmZl5\n22235eGHH5777LNP7rvvvjlz5sz8wx/+0OIzqY+NY73//vvnQQcdlBdeeGGuX7++1d2qm+23337T\nz48++mj29vbmpz/96czM/P73v59///d/n5n932XHHHPMprrPPvts7rbbbvnAAw/kggUL8owzzhjy\n+EPlE2BhVpFPI0u03qWvry8XLvTGGe1syZIl7Ldffe87K0mSNBpD5ZOIuCMz+4bb12UYkiRJUgHD\nsiRJklTAsKy6K9PSHkmS1N1Gm0sMy6qrsWPHsnr1agOzJElqucxk9erVjB07tuZjeJ9l1dXEiRNZ\nvnz5kPc3lSRJaraxY8cyceLEmvc3LKuuxowZwx577NHqbkiSJNWFyzAkSZKkAoZlSZIkqYBhWZIk\nSSpgWJYkSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSpgWJYkSZIKGJYlSZKkAoZlSZIkqYBhWZIkSSpQ\nU1iOiG9ExOMR8euC7RERX4mI+yPiVxHx+tF1U5IkSWq+WmeWvwVM38r2dwB7V/7MBP5Xje1IkiRJ\nLVNTWM71dCzqAAAS1UlEQVTMnwC/20qV44DLst9PgXER8epa2pIkSZJapVFrlicADw94vbxSJkmS\nJLWNbRt03BiiLIesGDGT/qUaTJo0qUHdGd7kT/xg2DrL5hyz6ef/vGGvqo77tqOW1lTfNmzDNkbf\nRhn7ZBu2YRvt1SfbaF4bZdWomeXlwO4DXk8EVg5VMTPnZmZfZvaNHz++Qd2RJEmSRq5RYflq4H2V\nu2L8KfBUZj7SoLYkSZKkhqhpGUZEzAOOAHaJiOXAp4ExAJn5L8C1wDuB+4FngffXo7OSJElSM9UU\nljPz1GG2J3B2TT2SJEmSSqJRF/h1vHMf3q6qeosb3A9JkiQ1jmG5RosffGhE9W/+yV9VVe9tR9XS\nG0mSJDWCYbnEDNiSJEmt1ai7YUiSJEltz5nlDuJMtCRJUn0ZlpvkA8+9rdVdGJIBW5IkqZhhuUne\nsd//U1W9xSW/f4bhWpIkdRPDshrOgC1JktqVYblJRnqruW5muJYkSWVhWC6xsq5zLiMDtiRJagTD\nco0mP/edquota2w3VKNawrWBXJKk7mNYLrGnr5pZXcU5SxrbEdXEcC1JUvszLJfYe2ZVNzzlvn+G\nRsKALUlSuRiWO4hrnLuPy0kkSWosw7KkrTJcS5K6mWG5xEZ6u7nvPvi5qup9lGm1dEeq2kgDtoFc\nklRWhuUO8s5fLm11F6SmqSZgG64lSaNlWO4gtVwQ6DpndQtnuyVJtTAsS1KdGMglqfMYliWpjRjI\nJam5DMtd7vtr1lVV7+zKf122IXU+A7kkvciw3OWOuuns4SsB4FMCJdVPMwK5IV5SPRiWO8hIbzUH\nPiVQkjYykEsaSs1hOSKmA18GeoBLMnPOoO2TgEuBcZU6n8jMa0fRV7Upl25IUr9mBHIDvFRfNYXl\niOgBLgbeDiwHFkTE1Zl594Bqfwf8e2b+r4jYH7gWmDzK/ratyc99p6p6yxrbjVHzwSeS1HnKsizG\nWXuVUa0zy28C7s/MBwAiYj5wHDAwLCewY+XnnYCVtXZS5TH25R9pdRckSdpCWWftDf3tr9awPAF4\neMDr5cAhg+rMBn4UEf8D2B74s6EOFBEzgZkAkyZNqrE7apZmXBDosg1JUrfqlNDfSb8k1BqWY4iy\nHPT6VOBbmfnFiDgU+LeIOCAzN2y2U+ZcYC5AX1/f4GOoZGq5INClG5IkqV3VGpaXA7sPeD2RLZdZ\nnAFMB8jM2yNiLLAL8HiNbaoBarmDhiRJUreoNSwvAPaOiD2AFcAM4LRBdR4C3gZ8KyL2A8YCq2rt\nqNpXM9Y5u3RDkiQ1Qk1hOTNfiIhzgOvpvy3cNzLzroj4DLAwM68GPgr8a0R8mP4lGn+dmS6z6EI+\n+ESSJLWrmu+zXLln8rWDyj414Oe7gcNq75q6lWucJUlSWfgEPzVcGZ8S6LINSZJUjW1a3QFJkiSp\nrJxZVumU9cEnzkZLktR9DMsakWbcaq6WCwJd5yxJkhrBsFxik5/7TlX1ljW2G5IkSV3LsKzSKeMF\ngbVw2YYkSe3PsKyOMNJ1zi7bkCRJ1TAsqyN0yoNPnI2WJKlcvHWcJEmSVMCZZXWlWm5PV8alG85E\nS5LUWIZlNVwzbjfXKRcFSpKkcjEsqyt1yhpnSZLUWIZlqUHKuGwDXLohSdJIGJalKpX1MdySJKlx\nDMuStsqZaElSNzMsdxAfj129Wi4IbMY657Iu3ZAkqVsZllU6zbh7hhrL2WhJUqcwLEsN4hpnSZLa\nn2FZamOdsmzDmWhJUlkZlqUqjXSds/dyliSp/RmWpRJx6Ub1nI2WJDWDYVnqMp2ydEOSpGYwLKsj\neAcNDceZaElSLWoOyxExHfgy0ANckplzhqjzHmA2kMAvM/O0WtuTukEZ1zl360y04VqSBDWG5Yjo\nAS4G3g4sBxZExNWZefeAOnsDs4DDMvPJiNi1Hh2W2kUtDz4ZKdc4S5LUWLXOLL8JuD8zHwCIiPnA\nccDdA+r8d+DizHwSIDMfH01H1Rg+9U+qH2ejJanz1BqWJwAPD3i9HDhkUJ19ACLiVvqXaszOzOsG\nHygiZgIzASZNmlRjdySVSbcu3ZAkdZ5aw3IMUZZDHHtv4AhgInBzRByQmWs22ylzLjAXoK+vb/Ax\nJG1FLWucXbpRHs5ES1L51RqWlwO7D3g9EVg5RJ2fZuY64MGIuIf+8Lygxjaluinr3TOasc65jJyJ\nrp4BW5Kaa5sa91sA7B0Re0TES4AZwNWD6lwFHAkQEbvQvyzjgVo7KkmSJDVbTTPLmflCRJwDXE//\neuRvZOZdEfEZYGFmXl3Z9ucRcTewHjgvM1fXq+OSajPSpRsu22hvzkRL0ujUfJ/lzLwWuHZQ2acG\n/JzARyp/JKmuqlm64bKN2hiwJelFPsFPamPdusZZkqRmMSxLqrsyLt3wIsLGcSZaUiczLGtEuvkh\nJmW9g0ajlfER3Gp/BmxJ7cKwLElDcCZakgSGZanrlHGdcxmXbahcnImW1CqGZUl1161LN5yNLhcD\ntqR6MCxLkoThWtLQDMtSg3TrBYHN0glLN5yJbn8GbKnzGZYlbVUz1jh367INdR/DtdR+DMuS1Eac\nje4+BmyptQzLarhuvjezyqMTlm1I1TBcS/VlWJbUlqpbuuGyDWeiVQ0DtlTMsCyVSKdcFFjGezmP\nlDPRUjHDtbqJYVmStBlno9UIBmy1K8OypK7gHTek9mK4VlkYliW1XCcs24DuXbrhTLTKoJZwbSBX\nNQzLkqSmM2CrHRmuu5NhWWpjnXJBYBk1Y9lGt85ES91kpAHbQF4+hmWVjvdlljSYM9FSMQN2YxmW\nJbWlatY5l32NsyS1guu7R8awLEl14tKNcnE2WmqdTgrXhmVJkjBcSxpazWE5IqYDXwZ6gEsyc05B\nvZOA7wFvzMyFtbYnqT669aLATrk93Ug5E91YBmyp89UUliOiB7gYeDuwHFgQEVdn5t2D6r0M+Bvg\nZ6PtqCR1Gh+UIknlV+vM8puA+zPzAYCImA8cB9w9qN7/C3we+FjNPZSq4B00pPpxNrpxnImW2k+t\nYXkC8PCA18uBQwZWiIiDgd0z85qIMCxLaitlXbbhbHT3MWBLrVVrWI4hynLTxohtgC8Bfz3sgSJm\nAjMBJk2aVGN3JDVKt65xVvWciS4Xw7VUX9vUuN9yYPcBrycCKwe8fhlwAHBTRCwD/hS4OiL6Bh8o\nM+dmZl9m9o0fP77G7kiSJEn1V+vM8gJg74jYA1gBzABO27gxM58Cdtn4OiJuAj7m3TAkdbIyLt0o\n67INZ6PLxdloqVhNYTkzX4iIc4Dr6b913Dcy866I+AywMDOvrmcnJUlSeRiu1U1qvs9yZl4LXDuo\n7FMFdY+otR2pEbx7RmO5zrm9lXE22pno9ldNwDZcq4x8gp8ktUgZl21IZeHstcrCsCxJ6kjORncX\nw7UaxbAsSRqVMi7bkKphwFY1DMuSWs41ztVz6UbjOBOt4dQSrg3k7c+wLElquk6ZjTZgq94M1+Vj\nWJaq5B00JEllZMBuLMOypLbk0o3qdMqyDWeipfpxOcnIGJYlSSoxA7baUSeFa8OyJGkzzkZL0osM\ny5IkdRBnoqX6MixLDeIFgeXiGmcNp5tnog3YUjHDsiRpVDpl2YYkDcWwLElSjbp1NtqZaHUTw7Ik\nDcFlG43VrbPR3RquwYCt9mVYlkrEdc6S1M9wrbIwLEuSSq9bZ6Khu2ejR8JwrUYxLEtSnbh0Q2ov\nBmxVw7AsSepI3Tob7Uy0VF+GZamNucZZUj0YsKvjTHR3MixLUou4bKNcunUmWo010oBtIC8fw7Ik\nSRoRZ6LLxYDdWIZlSWojzkaXi7PR1asuYBuum8FwPTKGZanLuM5ZUjtw9rpcujlg1xyWI2I68GWg\nB7gkM+cM2v4R4APAC8Aq4L9l5m9H0VdJktqaM9GNY7gul04K1zWF5YjoAS4G3g4sBxZExNWZefeA\naouAvsx8NiLOAj4PnDLaDkuSqueyjfZnwJZaq9aZ5TcB92fmAwARMR84DtgUljPzxgH1fwq8t9ZO\nSmodl21I6lTORqsatYblCcDDA14vBw7ZSv0zgB8OtSEiZgIzASZNmlRjdyRJ9eJsdHtzJrpxDNfd\nqdawHEOU5ZAVI94L9AFvHWp7Zs4F5gL09fUNeQxJktQ4BuzGMWC3v1rD8nJg9wGvJwIrB1eKiD8D\n/ifw1sz8Y41tSWozLt3oLs5ES/VjuC6fWsPyAmDviNgDWAHMAE4bWCEiDga+BkzPzMdH1UtJklQa\nzkSXiwG7sWoKy5n5QkScA1xP/63jvpGZd0XEZ4CFmXk1cCGwA/C9iAB4KDPfXad+S5LamLPR3aea\ngG24bg7D9cjUfJ/lzLwWuHZQ2acG/Pxno+iXpC7isg1Jgzl7XS4jDdidFMh9gp8kqfScidZwDNdq\nFMOypLZUzWz0ssZ3QyVmwJZUD4ZlSZLUlZyNVjUMy5K6guuiNRxnojUcw3V3MixLklQjA7aGY8Bu\nf4ZlSRqCM9GSWqGWcG0gbyzDsiRJTeJMtMrAcD0yhmVJqhNno9UIBmyVQTcHbMOyJEkdxHCtMuik\ncG1YlqQWcSZaZWHAlooZliWpjRiwVQaGa3UTw7IkdTDDtcrCgK12ZViWJEmlY7hWWRiWJUmbcTZa\n7chwrUYxLEuSRsVwrXY10oBtIO9OhmVJUtMZsNUtDNjtz7AsSSo9w7W6RS3h2kDeWIZlSVJHGmnA\nNpCrWxjIR8awLElSjQzY0tA6KVwbliVJapJawrWBXGotw7IkSR3EQC7Vl2FZkiSNSKMC+UjrD95H\nagTDsiRJanvNmFE3wHenmsNyREwHvgz0AJdk5pxB218KXAa8AVgNnJKZy2rvqiRJUntpRiAvYxud\npKawHBE9wMXA24HlwIKIuDoz7x5Q7Qzgycx8bUTMAD4HnDLaDkuSJKncOilcb1Pjfm8C7s/MBzLz\neWA+cNygOscBl1Z+vhx4W0REje1JkiRJTReZOfKdIk4CpmfmByqv/wo4JDPPGVDn15U6yyuvl1bq\nPDHoWDOBmZWXU4B7ajmRBtkFeGLYWuoUjnf3ccy7j2PeXRzv7jOSMX9NZo4frlKta5aHmiEenLqr\nqUNmzgXm1tiPhoqIhZnZ1+p+qDkc7+7jmHcfx7y7ON7dpxFjXusyjOXA7gNeTwRWFtWJiG2BnYDf\n1dieJEmS1HS1huUFwN4RsUdEvASYAVw9qM7VwOmVn08Cbsha1nxIkiRJLVLTMozMfCEizgGup//W\ncd/IzLsi4jPAwsy8Gvg68G8RcT/9M8oz6tXpJirl8hA1jOPdfRzz7uOYdxfHu/vUfcxrusBPkiRJ\n6ga1LsOQJEmSOp5hWZIkSSpgWB5CREyPiHsi4v6I+ESr+6P6i4hvRMTjlfuBbyx7RUT8OCLuq/z3\n5a3so+onInaPiBsjYklE3BURH6qUO+YdKiLGRsTPI+KXlTE/v1K+R0T8rDLm361cpK4OEhE9EbEo\nIq6pvHbMO1hELIuIxRFxZ0QsrJTV9bvdsDzIgEd5vwPYHzg1IvZvba/UAN8Cpg8q+wTwn5m5N/Cf\nldfqDC8AH83M/YA/Bc6u/H/tmHeuPwJHZebrgKnA9Ij4U+BzwJcqY/4kcEYL+6jG+BCwZMBrx7zz\nHZmZUwfcX7mu3+2G5S1V8yhvtbnM/Alb3vd74CPaLwWOb2qn1DCZ+Uhm/qLy89P0/0U6Ace8Y2W/\nZyovx1T+JHAUcHml3DHvMBExETgGuKTyOnDMu1Fdv9sNy1uaADw84PXySpk63ysz8xHoD1fAri3u\njxogIiYDBwM/wzHvaJV/jr8TeBz4MbAUWJOZL1Sq+P3eeS4CPg5sqLzeGce80yXwo4i4IyJmVsrq\n+t1e6+OuO1lVj+mW1H4iYgfgP4BzM/P3/ZNO6lSZuR6YGhHjgCuB/Yaq1txeqVEi4ljg8cy8IyKO\n2Fg8RFXHvLMclpkrI2JX4McR8Zt6N+DM8paqeZS3OtNjEfFqgMp/H29xf1RHETGG/qD8vzPzikqx\nY94FMnMNcBP969XHRcTGiSK/3zvLYcC7I2IZ/Usoj6J/ptkx72CZubLy38fp/6X4TdT5u92wvKVq\nHuWtzjTwEe2nA99vYV9UR5V1i18HlmTmPw3Y5Jh3qIgYX5lRJiJ6gT+jf636jcBJlWqOeQfJzFmZ\nOTEzJ9P/d/cNmfmXOOYdKyK2j4iXbfwZ+HPg19T5u90n+A0hIt5J/2+jGx/l/Y8t7pLqLCLmAUcA\nuwCPAZ8GrgL+HZgEPAScnJmDLwJUG4qIw4GbgcW8uJbxk/SvW3bMO1BEHET/hT099E8M/XtmfiYi\n9qR/1vEVwCLgvZn5x9b1VI1QWYbxscw81jHvXJWxvbLyclvgO5n5jxGxM3X8bjcsS5IkSQVchiFJ\nkiQVMCxLkiRJBQzLkiRJUgHDsiRJklTAsCxJkiQVMCxLkiRJBQzLkiRJUoH/C1H8RHw0JG6bAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2c7fe8dbf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P = [\n",
    "    [90.34, 5.62, 0.39, 0.08, 0.03, 0, 0, 0, 3.54],\n",
    "    [0.64, 88.78, 6.72, 0.47, 0.06, 0.09, 0.02, 0.01, 3.21],\n",
    "    [0.07, 2.16, 87.94, 4.97, 0.47, 0.19, 0.01, 0.04, 4.15],\n",
    "    [0.03, 0.24, 4.56, 84.26, 4.19, 0.76, 0.15, 0.22, 5.59],\n",
    "    [0.03, 0.06, 0.4, 6.09, 76.09, 6.82, 0.96, 0.98, 8.57],\n",
    "    [0, 0.09, 0.29, 0.41, 5.11, 74.62, 3.43, 5.3, 10.75],\n",
    "    [0.13, 0, 0.26, 0.77, 1.66, 8.93, 53.19, 21.94, 13.12],\n",
    "    [0, 0, 0, 0, 1, 3.1, 9.29, 51.29, 35.32],\n",
    "    [0, 0, 0, 0, 0, 0.1, 8.55, 74.06, 17.29]\n",
    "]\n",
    "\n",
    "P = np.array(P)/100\n",
    "outcomes = ['AAA', 'AA', 'A', 'BBB', 'BB', 'B', 'CCC', 'D', 'N.R.']\n",
    "initial = [1, 0, 0, 0, 0, 0, 0, 0, 0]\n",
    "\n",
    "ratings = MarkovChain(P,initial, outcomes)\n",
    "\n",
    "fig, (ax0,ax1) = plt.subplots(2, 1, sharex=True, figsize=[12,8])\n",
    "ratings.plot(ax0, T, 0)\n",
    "ratings.plot_probability(ax1,T)\n",
    "ratings.print_stationary()"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}