FIX: markov: Respect dtype of P in cdfs (#592)

Author: oyamad

quantecon/markov/core.py

@@ -273,7 +273,6 @@ def get_index(self, value):
         else:
             raise ValueError('invalid value')
 
-
     def _get_index(self, value):
         """
         Return the index of the given value in `state_values`.
@@ -420,7 +419,7 @@ def stationary_distributions(self):
     def cdfs(self):
         if (self._cdfs is None) and not self.is_sparse:
             # See issue #137#issuecomment-96128186
-            cdfs = np.empty((self.n, self.n), order='C')
+            cdfs = np.empty((self.n, self.n), order='C', dtype=self.P.dtype)
             np.cumsum(self.P, axis=-1, out=cdfs)
             self._cdfs = cdfs
         return self._cdfs
@@ -431,7 +430,7 @@ def cdfs1d(self):
             data = self.P.data
             indptr = self.P.indptr
 
-            cdfs1d = np.empty(self.P.nnz, order='C')
+            cdfs1d = np.empty(self.P.nnz, order='C', dtype=data.dtype)
             for i in range(self.n):
                 cdfs1d[indptr[i]:indptr[i+1]] = \
                     data[indptr[i]:indptr[i+1]].cumsum()

quantecon/markov/tests/test_core.py

@@ -10,7 +10,9 @@
 import numpy as np
 from scipy import sparse
 import itertools
-from numpy.testing import assert_allclose, assert_array_equal, assert_raises
+from numpy.testing import (
+    assert_allclose, assert_array_equal, assert_array_less, assert_raises
+)
 from nose.tools import eq_, ok_, raises
 
 from quantecon.markov import (
@@ -339,6 +341,33 @@ def test_simulate_for_matrices_with_C_F_orders():
     assert_array_equal(computed_F, sample_path)
 
 
+def test_simulate_issue591():
+    """
+    Test MarkovChasin.simulate for P with dtype=np.float32
+    https://github.com/QuantEcon/QuantEcon.py/issues/591
+    """
+    num_states = 5
+    transition_states = 4
+
+    transition_seed = 2
+    random_state = np.random.RandomState(transition_seed)
+    transitions = random_state.uniform(0., 1., transition_states)
+    transitions /= np.sum(transitions)
+    P = np.zeros((num_states, num_states), dtype=np.float32)
+    P[0, :transition_states] = transitions
+    P[1:, 0] = 1.
+    mc = MarkovChain(P=P)
+
+    simulate_seed = 22220
+    ts_length = 10000
+    seq = mc.simulate(
+        ts_length=ts_length, init=0, num_reps=1, random_state=simulate_seed
+    )
+    max_state_in_seq = np.max(seq)
+
+    assert_array_less(max_state_in_seq, num_states)
+
+
 def test_mc_sample_path():
     P = [[0.4, 0.6], [0.2, 0.8]]
     Ps = [P, sparse.csr_matrix(P)]