1"""
2Income mobility measures.
3"""
4
5__author__ = "Wei Kang <weikang9009@gmail.com>, Sergio J. Rey <sjsrey@gmail.com>"
6
7__all__ = ["markov_mobility"]
8
9import numpy as np
10import numpy.linalg as la
11
12
13def markov_mobility(p, measure="P", ini=None):
14    """
15    Markov-based mobility index.
16
17    Parameters
18    ----------
19    p       : array
20              (k, k), Markov transition probability matrix.
21    measure : string
22              If measure= "P",
23              :math:`M_{P} = \\frac{m-\sum_{i=1}^m P_{ii}}{m-1}`;
24              if measure = "D",
25              :math:`M_{D} = 1 - |\det(P)|`,
26              where :math:`\det(P)` is the determinant of :math:`P`;
27              if measure = "L2",
28              :math:`M_{L2} = 1  - |\lambda_2|`,
29              where :math:`\lambda_2` is the second largest eigenvalue of
30              :math:`P`;
31              if measure = "B1",
32              :math:`M_{B1} = \\frac{m-m \sum_{i=1}^m \pi_i P_{ii}}{m-1}`,
33              where :math:`\pi` is the initial income distribution;
34              if measure == "B2",
35              :math:`M_{B2} = \\frac{1}{m-1} \sum_{i=1}^m \sum_{
36              j=1}^m \pi_i P_{ij} |i-j|`,
37              where :math:`\pi` is the initial income distribution.
38    ini     : array
39              (k,), initial distribution. Need to be specified if
40              measure = "B1" or "B2". If not,
41              the initial distribution would be treated as a uniform
42              distribution.
43
44    Returns
45    -------
46    mobi    : float
47              Mobility value.
48
49    Notes
50    -----
51    The mobility indices are based on :cite:`Formby:2004fk`.
52
53    Examples
54    --------
55    >>> import numpy as np
56    >>> import libpysal
57    >>> import mapclassify as mc
58    >>> from giddy.markov import Markov
59    >>> from giddy.mobility import markov_mobility
60    >>> f = libpysal.io.open(libpysal.examples.get_path("usjoin.csv"))
61    >>> pci = np.array([f.by_col[str(y)] for y in range(1929,2010)])
62    >>> q5 = np.array([mc.Quantiles(y).yb for y in pci]).transpose()
63    >>> m = Markov(q5)
64    The Markov Chain is irreducible and is composed by:
65    1 Recurrent class (indices):
66    [0 1 2 3 4]
67    0 Transient classes.
68    The Markov Chain has 0 absorbing states.
69    >>> m.p
70    array([[0.91011236, 0.0886392 , 0.00124844, 0.        , 0.        ],
71           [0.09972299, 0.78531856, 0.11080332, 0.00415512, 0.        ],
72           [0.        , 0.10125   , 0.78875   , 0.1075    , 0.0025    ],
73           [0.        , 0.00417827, 0.11977716, 0.79805014, 0.07799443],
74           [0.        , 0.        , 0.00125156, 0.07133917, 0.92740926]])
75
76    (1) Estimate Shorrock1 mobility index:
77
78    >>> mobi_1 = markov_mobility(m.p, measure="P")
79    >>> print("{:.5f}".format(mobi_1))
80    0.19759
81
82    (2) Estimate Shorrock2 mobility index:
83
84    >>> mobi_2 = markov_mobility(m.p, measure="D")
85    >>> print("{:.5f}".format(mobi_2))
86    0.60685
87
88    (3) Estimate Sommers and Conlisk mobility index:
89
90    >>> mobi_3 = markov_mobility(m.p, measure="L2")
91    >>> print("{:.5f}".format(mobi_3))
92    0.03978
93
94    (4) Estimate Bartholomew1 mobility index (note that the initial
95    distribution should be given):
96
97    >>> ini = np.array([0.1,0.2,0.2,0.4,0.1])
98    >>> mobi_4 = markov_mobility(m.p, measure = "B1", ini=ini)
99    >>> print("{:.5f}".format(mobi_4))
100    0.22777
101
102    (5) Estimate Bartholomew2 mobility index (note that the initial
103    distribution should be given):
104
105    >>> ini = np.array([0.1,0.2,0.2,0.4,0.1])
106    >>> mobi_5 = markov_mobility(m.p, measure = "B2", ini=ini)
107    >>> print("{:.5f}".format(mobi_5))
108    0.04637
109
110    """
111
112    p = np.array(p)
113    k = p.shape[1]
114    if measure == "P":
115        t = np.trace(p)
116        mobi = (k - t) / (k - 1)
117    elif measure == "D":
118        mobi = 1 - abs(la.det(p))
119    elif measure == "L2":
120        w, v = la.eig(p)
121        eigen_value_abs = abs(w)
122        mobi = 1 - np.sort(eigen_value_abs)[-2]
123    elif measure == "B1":
124        if ini is None:
125            ini = 1.0 / k * np.ones(k)
126        mobi = (k - k * np.sum(ini * np.diag(p))) / (k - 1)
127    elif measure == "B2":
128        mobi = 0
129        if ini is None:
130            ini = 1.0 / k * np.ones(k)
131        for i in range(k):
132            for j in range(k):
133                mobi = mobi + ini[i] * p[i, j] * abs(i - j)
134        mobi = mobi / (k - 1)
135
136    return mobi
137