{
 "metadata": {
  "name": "",
  "signature": "sha256:48a49797a062d3c22eaaebbfe4ed0ccbe380f194e3d651a6a837cab4107ce770"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import numpy as np\n",
      "import quantecon as qe\n",
      "from mc_tools import mc_compute_stationary, mc_sample_path\n",
      "from quantecon.mc_tools import MarkovChain\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Scaled Gaussian Elimination"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ScaledGE(A):\n",
      "    A = np.array(A, dtype=float)\n",
      "    n = A.shape[0]\n",
      "    x = np.zeros(n) \n",
      "#Reduction stage\n",
      "# 1st: scaling column vector for i-th column\n",
      "# 2nd: reduction.Obtain a submatrix for each i\n",
      "    for i in range(n-1):\n",
      "      \n",
      "         A[i+1:n, i] /= -A[i,i]\n",
      "         A[i+1:n, i+1:n] += np.dot(A[i+1:n, i:i+1], A[i:i+1, i+1:n])\n",
      "            \n",
      "#Backward Substitution        \n",
      "    x[n-1]=1\n",
      "    for i in range(n-2, -1, -1):\n",
      "      x[i] = np.sum((x[j] * A[j, i] for j in range(i+1, n)))\n",
      "        \n",
      "    x /= np.sum(x)\n",
      "\n",
      "    return x    \n",
      "      "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "GTH algorithm"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def GTH(A):\n",
      "    A = np.array(A, dtype=float)\n",
      "    n = A.shape[0]\n",
      "    x = np.zeros(n)\n",
      "#Reduction stage\n",
      "# 1st: scaling column vector for i-th column.Scaling is conducted by sum of other components in each row.\n",
      "# 2nd: reduction.Obtain a submatrix for each i\n",
      "# Key point is sum of row vector equals to 0. Keep an attention to this point.\n",
      "    for i in range(n-1):\n",
      "        scale = np.sum(A[i, i+1:n])\n",
      "        \n",
      "        A[i+1:n, i] /= scale\n",
      "        A[i+1:n, i+1:n] += np.dot(A[i+1:n, i:i+1], A[i:i+1, i+1:n])\n",
      "\n",
      "    x[n-1] = 1\n",
      "    for i in range(n-2, -1, -1):\n",
      "        x[i] = np.sum((x[j] * A[j, i] for j in range(i+1, n)))\n",
      "\n",
      "    x /= np.sum(x)\n",
      "\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Apply the algorithms"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q=0.5\n",
      "\n",
      "for i in range(8,18):\n",
      "     e=1e-1**i\n",
      "\n",
      "     P =np.array([[-(q+e), q     , e   ],\n",
      "             [q     , -(q+e), e   ],\n",
      "             [e     , e     , -2*e]])\n",
      "    \n",
      "     print('e = 1e-{0}'.format(i))\n",
      "     print(ScaledGE(P))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "e = 1e-8\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-9\n",
        "[ 0.33333334  0.33333334  0.33333333]\n",
        "e = 1e-10\n",
        "[ 0.33333332  0.33333332  0.33333335]\n",
        "e = 1e-11\n",
        "[ 0.33333332  0.33333332  0.33333335]\n",
        "e = 1e-12\n",
        "[ 0.33333579  0.33333579  0.33332842]\n",
        "e = 1e-13\n",
        "[ 0.33329879  0.33329879  0.33340243]\n",
        "e = 1e-14\n",
        "[ 0.33342217  0.33342217  0.33315567]\n",
        "e = 1e-15\n",
        "[ 0.33342217  0.33342217  0.33315567]\n",
        "e = 1e-16\n",
        "[ 0.32152035  0.32152035  0.3569593 ]\n",
        "e = 1e-17\n",
        "[ nan  nan  -0.]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q=0.5\n",
      "\n",
      "for i in range(8,18):\n",
      "     e=1e-1**i\n",
      "\n",
      "     P =np.array([[-(q+e), q     , e   ],\n",
      "             [q     , -(q+e), e   ],\n",
      "             [e     , e     , -2*e]])\n",
      "    \n",
      "     print('e = 1e-{0}'.format(i))\n",
      "     print(GTH(P))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "e = 1e-8\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-9\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-10\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-11\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-12\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-13\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-14\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-15\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-16\n",
        "[ 0.33333333  0.33333333  0.33333333]\n",
        "e = 1e-17\n",
        "[ 0.33333333  0.33333333  0.33333333]\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "KMR model analysis"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Parameter\n",
      "p = 1/3\n",
      "N = 3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def KMR_2x2_P_sequential(N, p, epsilon):\n",
      "          P = np.zeros((N+1, N+1), dtype=float)\n",
      "          \n",
      "          P[0, 0], P[0, 1] = 1 - epsilon * (1/2), epsilon * (1/2)    \n",
      "       \n",
      "          P[N, N-1], P[N, N] = epsilon * (1/2), 1 - epsilon * (1/2) \n",
      "     \n",
      "          \n",
      "          for x in range(1, N):\n",
      "             \n",
      "                if x/N >p:\n",
      "                    P[x, x-1] = (x/N)*epsilon * (1/2)\n",
      "                    P[x, x+1] = (1-(x/N))*(1 - epsilon*(1/2))\n",
      "                    P[x, x] = 1 - P[x, x-1] - P[x, x+1] \n",
      "                elif x/N<p:\n",
      "                    P[x, x-1] =(x/N)*(1 - epsilon*(1/2))\n",
      "                    P[x, x+1] =(1-(x/N))*epsilon * (1/2)\n",
      "                    P[x, x] = 1 - P[x, x-1] - P[x, x+1]     \n",
      "                else:\n",
      "                    P[x, x-1] = (x/N)*(1/2)\n",
      "                    P[x, x+1] = (1-(x/N))*(1/2)\n",
      "                    P[x, x] = 1 - P[x, x-1] - P[x, x+1]           \n",
      "          return P\n",
      "               \n",
      "               "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(2,17):\n",
      "     epsilon=1e-1**i\n",
      "\n",
      "     kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)\n",
      "\n",
      "     print('epsilon = 1-{0}'.format(i))\n",
      "     print(kmr_seq)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "epsilon = 1-2\n",
        "[[ 0.995       0.005       0.          0.        ]\n",
        " [ 0.16666667  0.5         0.33333333  0.        ]\n",
        " [ 0.          0.00333333  0.665       0.33166667]\n",
        " [ 0.          0.          0.005       0.995     ]]\n",
        "epsilon = 1-3\n",
        "[[  9.99500000e-01   5.00000000e-04   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-04   6.66500000e-01   3.33166667e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-04   9.99500000e-01]]\n",
        "epsilon = 1-4\n",
        "[[  9.99950000e-01   5.00000000e-05   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-05   6.66650000e-01   3.33316667e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-05   9.99950000e-01]]\n",
        "epsilon = 1-5\n",
        "[[  9.99995000e-01   5.00000000e-06   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-06   6.66665000e-01   3.33331667e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-06   9.99995000e-01]]\n",
        "epsilon = 1-6\n",
        "[[  9.99999500e-01   5.00000000e-07   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-07   6.66666500e-01   3.33333167e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-07   9.99999500e-01]]\n",
        "epsilon = 1-7\n",
        "[[  9.99999950e-01   5.00000000e-08   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-08   6.66666650e-01   3.33333317e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-08   9.99999950e-01]]\n",
        "epsilon = 1-8\n",
        "[[  9.99999995e-01   5.00000000e-09   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-09   6.66666665e-01   3.33333332e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-09   9.99999995e-01]]\n",
        "epsilon = 1-9\n",
        "[[  9.99999999e-01   5.00000000e-10   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-10   6.66666666e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-10   9.99999999e-01]]\n",
        "epsilon = 1-10\n",
        "[[  1.00000000e+00   5.00000000e-11   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-11   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-11   1.00000000e+00]]\n",
        "epsilon = 1-11\n",
        "[[  1.00000000e+00   5.00000000e-12   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-12   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-12   1.00000000e+00]]\n",
        "epsilon = 1-12\n",
        "[[  1.00000000e+00   5.00000000e-13   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-13   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-13   1.00000000e+00]]\n",
        "epsilon = 1-13\n",
        "[[  1.00000000e+00   5.00000000e-14   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-14   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-14   1.00000000e+00]]\n",
        "epsilon = 1-14\n",
        "[[  1.00000000e+00   5.00000000e-15   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-15   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-15   1.00000000e+00]]\n",
        "epsilon = 1-15\n",
        "[[  1.00000000e+00   5.00000000e-16   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-16   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-16   1.00000000e+00]]\n",
        "epsilon = 1-16\n",
        "[[  1.00000000e+00   5.00000000e-17   0.00000000e+00   0.00000000e+00]\n",
        " [  1.66666667e-01   5.00000000e-01   3.33333333e-01   0.00000000e+00]\n",
        " [  0.00000000e+00   3.33333333e-17   6.66666667e-01   3.33333333e-01]\n",
        " [  0.00000000e+00   0.00000000e+00   5.00000000e-17   1.00000000e+00]]\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Compute the stationary distribution for each KMR by using the qe code."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(2,17):\n",
      "     epsilon=1e-1**i\n",
      "\n",
      "     kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)\n",
      "     kmr_seq_m=qe.MarkovChain(kmr_seq)\n",
      "     kmr_seq_s=kmr_seq_m.stationary_distributions\n",
      "     print('epsilon = 1-{0}'.format(i))\n",
      "     print(kmr_seq_s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "epsilon = 1-2\n",
        "[[  4.92538049e-03   1.47761415e-04   1.47761415e-02   9.80150717e-01]]\n",
        "epsilon = 1-3\n",
        "[[  4.99250376e-04   1.49775113e-06   1.49775113e-03   9.98001501e-01]]\n",
        "epsilon = 1-4\n",
        "[[  4.99925004e-05   1.49977501e-08   1.49977501e-04   9.99800015e-01]]\n",
        "epsilon = 1-5\n",
        "[[  4.99992500e-06   1.49997750e-10   1.49997750e-05   9.99980000e-01]]\n",
        "epsilon = 1-6\n",
        "[[  4.99999250e-07   1.49999775e-12   1.49999775e-06   9.99998000e-01]]\n",
        "epsilon = 1-7\n",
        "[[  4.99999925e-08   1.49999978e-14   1.49999978e-07   9.99999800e-01]]\n",
        "epsilon = 1-8\n",
        "[[  4.99999993e-09   1.49999998e-16   1.49999998e-08   9.99999980e-01]]\n",
        "epsilon = 1-9\n",
        "[[  4.99999999e-10   1.50000000e-18   1.50000000e-09   9.99999998e-01]]\n",
        "epsilon = 1-10\n",
        "[[  5.00000000e-11   1.50000000e-20   1.50000000e-10   1.00000000e+00]]\n",
        "epsilon = 1-11\n",
        "[[  5.00000000e-12   1.50000000e-22   1.50000000e-11   1.00000000e+00]]\n",
        "epsilon = 1-12\n",
        "[[  5.00000000e-13   1.50000000e-24   1.50000000e-12   1.00000000e+00]]\n",
        "epsilon = 1-13\n",
        "[[  5.00000000e-14   1.50000000e-26   1.50000000e-13   1.00000000e+00]]\n",
        "epsilon = 1-14\n",
        "[[  5.00000000e-15   1.50000000e-28   1.50000000e-14   1.00000000e+00]]\n",
        "epsilon = 1-15\n",
        "[[  5.00000000e-16   1.50000000e-30   1.50000000e-15   1.00000000e+00]]\n",
        "epsilon = 1-16\n",
        "[[  5.00000000e-17   1.50000000e-32   1.50000000e-16   1.00000000e+00]]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Next results are obtained by our algorithms."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(2,17):\n",
      "     epsilon=1e-1**i\n",
      "     kmr_seq=KMR_2x2_P_sequential(N, p, epsilon)   \n",
      "     I=np.identity(kmr_seq.shape[0]) # used for SGE algorithm.\n",
      "     kmr_seq_SGE=ScaledGE(I-kmr_seq)\n",
      "     kmr_seq_GTH=GTH(I-kmr_seq)\n",
      "\n",
      "     print('epsilon = 1-{0}'.format(i))  \n",
      "     print(kmr_seq_SGE)\n",
      "     print(kmr_seq_GTH)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "epsilon = 1-2\n",
        "[  4.92538049e-03   1.47761415e-04   1.47761415e-02   9.80150717e-01]\n",
        "[  4.92538049e-03   1.47761415e-04   1.47761415e-02   9.80150717e-01]\n",
        "epsilon = 1-3\n",
        "[  4.99250376e-04   1.49775113e-06   1.49775113e-03   9.98001501e-01]\n",
        "[  4.99250376e-04   1.49775113e-06   1.49775113e-03   9.98001501e-01]\n",
        "epsilon = 1-4\n",
        "[  4.99925004e-05   1.49977501e-08   1.49977501e-04   9.99800015e-01]\n",
        "[  4.99925004e-05   1.49977501e-08   1.49977501e-04   9.99800015e-01]\n",
        "epsilon = 1-5\n",
        "[  4.99992500e-06   1.49997750e-10   1.49997750e-05   9.99980000e-01]\n",
        "[  4.99992500e-06   1.49997750e-10   1.49997750e-05   9.99980000e-01]\n",
        "epsilon = 1-6\n",
        "[  4.99999250e-07   1.49999775e-12   1.49999775e-06   9.99998000e-01]\n",
        "[  4.99999250e-07   1.49999775e-12   1.49999775e-06   9.99998000e-01]\n",
        "epsilon = 1-7\n",
        "[  4.99999925e-08   1.49999977e-14   1.49999978e-07   9.99999800e-01]\n",
        "[  4.99999925e-08   1.49999978e-14   1.49999978e-07   9.99999800e-01]\n",
        "epsilon = 1-8\n",
        "[  4.99999997e-09   1.49999998e-16   1.49999998e-08   9.99999980e-01]\n",
        "[  4.99999993e-09   1.49999998e-16   1.49999998e-08   9.99999980e-01]\n",
        "epsilon = 1-9\n",
        "[  4.99999937e-10   1.49999994e-18   1.50000000e-09   9.99999998e-01]\n",
        "[  4.99999999e-10   1.50000000e-18   1.50000000e-09   9.99999998e-01]\n",
        "epsilon = 1-10\n",
        "[  4.99999938e-11   1.49999994e-20   1.50000000e-10   1.00000000e+00]\n",
        "[  5.00000000e-11   1.50000000e-20   1.50000000e-10   1.00000000e+00]\n",
        "epsilon = 1-11\n",
        "[  4.99999938e-12   1.49999994e-22   1.50000000e-11   1.00000000e+00]\n",
        "[  5.00000000e-12   1.50000000e-22   1.50000000e-11   1.00000000e+00]\n",
        "epsilon = 1-12\n",
        "[  4.99933333e-13   1.49993333e-24   1.50000000e-12   1.00000000e+00]\n",
        "[  5.00000000e-13   1.50000000e-24   1.50000000e-12   1.00000000e+00]\n",
        "epsilon = 1-13\n",
        "[  5.00600178e-14   1.50060018e-26   1.50000000e-13   1.00000000e+00]\n",
        "[  5.00000000e-14   1.50000000e-26   1.50000000e-13   1.00000000e+00]\n",
        "epsilon = 1-14\n",
        "[  5.00600178e-15   1.50060018e-28   1.50000000e-14   1.00000000e+00]\n",
        "[  5.00000000e-15   1.50000000e-28   1.50000000e-14   1.00000000e+00]\n",
        "epsilon = 1-15\n",
        "[  4.29061342e-16   1.42906134e-30   1.50000000e-15   1.00000000e+00]\n",
        "[  5.00000000e-16   1.50000000e-30   1.50000000e-15   1.00000000e+00]\n",
        "epsilon = 1-16\n",
        "[ nan  nan  nan  nan]\n",
        "[  5.00000000e-17   1.50000000e-32   1.50000000e-16   1.00000000e+00]\n"
       ]
      }
     ],
     "prompt_number": 18
    }
   ],
   "metadata": {}
  }
 ]
}