In [1]:
from __future__ import division
import numpy as np
from quantecon.mc_tools import mc_compute_stationary



def set_tran_matrix(N=6,p=1/3,epsilon=0.01):  #setting up the transition matrix A
	"""
	N :the number of players
	p : the turning point as to which equilibrium is prefered
	epsilon : tthe possibility that a man becomes "irrational"

	if the ratio of players playing 0 exceeds p, we end up the equilibrium(0,0) in an unperturbed game
	I think it's better if I can compute the value of p in this program.For now,we assume it's exogenous.
	"""
	A = np.zeros((N+1,N+1))
	for n in range(1, N):
		A[n, n-1] = \
	    	(n/N) * (epsilon * (1/2) +
                     (1 - epsilon) * (((n)/(N) < p) + ((n)/(N) == p) * (1/2))
                     )
		A[n, n+1] = \
            ((N-n)/N) * (epsilon * (1/2) +
                         (1 - epsilon) * ((n/(N) > p) + (n/(N) == p) * (1/2))
                         )
		A[n,n] = 1- A[n,n-1] - A[n,n+1]

	A[0,0] = 1 - epsilon/2
	A[0,1] = epsilon/2

	A[N,N-1] = epsilon/2
	A[N,N] = 1 - epsilon/2
	
	return A

In [2]:
P=set_tran_matrix()

In [3]:
mc_compute_stationary(P)


Out[3]:
array([  2.45031189e-05,   7.38787504e-07,   1.84696875e-08,
         2.46262500e-06,   3.67546782e-04,   2.92567238e-02,
         9.70348006e-01])

In [4]:
P=set_tran_matrix(N=25)

In [5]:
mc_compute_stationary(P)


Out[5]:
array([[  2.28974045e-01,  -6.60219850e-26],
       [  2.87655835e-02,  -8.29420748e-27],
       [  1.73460805e-03,  -5.00146613e-28],
       [  6.68274458e-05,  -1.92641995e-29],
       [  1.84698971e-06,  -5.29103886e-31],
       [  3.89816906e-08,  -8.50800510e-33],
       [  6.52958230e-10,   2.07333066e-33],
       [  8.90489368e-12,   1.92757373e-33],
       [  9.96225258e-14,   1.68367666e-33],
       [  5.32184164e-16,   2.97148825e-31],
       [  5.68020111e-16,   9.48767575e-29],
       [  6.09076129e-16,   2.57463425e-26],
       [  6.55492459e-16,   5.97744273e-24],
       [  7.05784387e-16,   1.18951110e-21],
       [  7.70086588e-16,   2.02896608e-19],
       [  8.72133584e-16,   2.96093783e-17],
       [  3.66805592e-15,   3.68266643e-15],
       [  2.88337478e-13,   3.87979740e-13],
       [  2.54097246e-11,   3.43146526e-11],
       [  1.86284758e-09,   2.51580584e-09],
       [  1.11211918e-07,   1.50193609e-07],
       [  5.26932652e-06,   7.11631623e-06],
       [  1.90653814e-04,   2.57481260e-04],
       [  4.94870987e-03,   6.68331792e-03],
       [  8.20661053e-02,   1.10831689e-01],
       [  6.53246198e-01,   8.82220243e-01]])

In [6]:
P=set_tran_matrix(N=26)

In [7]:
mc_compute_stationary(P)


C:\Users\Atsushi\Anaconda\lib\site-packages\quantecon\mc_tools.py:142: UserWarning: Elements of your invariant distribution were negative; Re-trying with additional precision
  "Re-trying with additional precision")
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-7-86586f448331> in <module>()
----> 1 mc_compute_stationary(P)

C:\Users\Atsushi\Anaconda\lib\site-packages\quantecon\mc_tools.pyc in mc_compute_stationary(P, precision, tol)
    143 
    144         if precision is None:
--> 145             stationary_dists = mc_compute_stationary(P, precision=18, tol=tol)
    146 
    147         elif precision is not None:

C:\Users\Atsushi\Anaconda\lib\site-packages\quantecon\mc_tools.pyc in mc_compute_stationary(P, precision, tol)
    127 
    128         with mp.workdps(precision):
--> 129             eigvals, eigvecs = mp.eig(mp.matrix(P), left=True, right=False)
    130 
    131             for ind, el in enumerate(eigvals):

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\matrices\eigen.pyc in eig(ctx, A, left, right, overwrite_a)
    753             A[y,x] = 0
    754 
--> 755     hessenberg_qr(ctx, A, Q)
    756 
    757     E = [0 for i in xrange(n)]

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\matrices\eigen.pyc in hessenberg_qr(ctx, A, Q)
    476             totalits += 1
    477 
--> 478             qr_step(ctx, n0, n1, A, Q, shift)
    479 
    480             if its > maxits:

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\matrices\eigen.pyc in qr_step(ctx, n0, n1, A, Q, shift)
    370                 y = Q[k,j+2]
    371                 Q[k,j+1] =           c  * x +          s  * y
--> 372                 Q[k,j+2] = -ctx.conj(s) * x + ctx.conj(c) * y
    373 
    374 

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\ctx_mp_python.pyc in __mul__(self, other)

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\libmp\libmpf.pyc in python_mpf_mul(s, t, prec, rnd)
    866         bc += int(man>>bc)
    867         if prec:
--> 868             return normalize1(sign, man, sexp+texp, bc, prec, rnd)
    869         else:
    870             return (sign, man, sexp+texp, bc)

C:\Users\Atsushi\Anaconda\lib\site-packages\sympy\mpmath\libmp\libmpf.pyc in _normalize1(sign, man, exp, bc, prec, rnd)
    228     bc = prec
    229     # Strip trailing bits
--> 230     if not man & 1:
    231         t = trailtable[int(man & 255)]
    232         if not t:

KeyboardInterrupt: 

In [ ]: