Add reshuffling method [WIP]

Documentation/CHANGELOG.md

@@ -15,6 +15,7 @@ Release Date: TBD
 ### Major Changes
 
 - Adds `HARK.core.AgentPopulation` class to represent a population of agents with ex-ante heterogeneous parametrizations as distributions. [#1237](https://github.com/econ-ark/HARK/pull/1237)
+- Adds methods to improve efficiency of Monte Carlo simulations. [#1244] (https://github.com/econ-ark/HARK/pull/1244)
 
 ### Minor Changes
 

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -1561,6 +1561,8 @@ def prepare_to_calc_EndOfPrdvP(self):
     "T_cycle": 1,  # Number of periods in the cycle for this agent type
     "PerfMITShk": False,
     # Do Perfect Foresight MIT Shock: Forces Newborns to follow solution path of the agent he/she replaced when True
+    "reshuffle": False, # Whether to use reshuffling method for Monte Carlo Simulation
+    "perf_reshuffle": False, #reshuffle must be set to true to use this. Whether to have reshuffling method be perfectly across both newborns and nonnewborns. 
 }
 
 
@@ -1792,9 +1794,20 @@ def sim_death(self):
         # they die.
         # See: https://github.com/econ-ark/HARK/pull/981
 
-        DeathShks = Uniform(seed=self.RNG.integers(0, 2**31 - 1)).draw(
-            N=self.AgentCount
-        )
+        if self.reshuffle ==True:
+
+            DeathShks_dstn = Uniform(seed=self.RNG.integers(0, 2**31 - 1))._approx_equiprobable(
+                N=self.AgentCount
+            )
+
+            DeathShks = DeathShks_dstn.draw(
+                N=self.AgentCount, exact_match = self.reshuffle
+            )
+        else:
+
+            DeathShks = Uniform(seed=self.RNG.integers(0, 2**31 - 1)).draw(
+                N=self.AgentCount
+            )
         which_agents = DeathShks < DiePrb
         if self.T_age is not None:  # Kill agents that have lived for too many periods
             too_old = self.t_age >= self.T_age
@@ -2227,6 +2240,7 @@ def reset_rng(self):
                 dstn.reset()
 
     def get_shocks(self):
+
         """
         Gets permanent and transitory income shocks for this period.  Samples from IncShkDstn for
         each period in the cycle.
@@ -2240,34 +2254,82 @@ def get_shocks(self):
         -------
         None
         """
+
         NewbornTransShk = (
             self.NewbornTransShk
         )  # Whether Newborns have transitory shock. The default is False.
 
         PermShkNow = np.zeros(self.AgentCount)  # Initialize shock arrays
         TranShkNow = np.zeros(self.AgentCount)
+
         newborn = self.t_age == 0
+
+        if self.reshuffle ==True: # when true permanent, transitory, and  unemployment shocks are perfectly distributed across population
+
+            def get_mult_fac(prb):
+                i = 0
+                while isinstance(prb,float):
+                    prb = prb*10
+                    i += 1
+                    if ( abs(prb - round(prb)) )<(1e-10):
+                        break
+                return i
+
+            if self.UnempPrb>0:
+
+                potential_fac = 1/self.UnempPrb
+
+                if (potential_fac).is_integer() ==True:
+                    fac = potential_fac
+                else:
+                    i_1 = get_mult_fac(self.UnempPrb)
+                    fac = 10**i_1
+            else:
+                fac = 1.0
+
+
+            lcm = (self.PermShkCount * self.TranShkCount) * fac # minimum multiple for both the newborns and and oldborns individuals
+            if (self.AgentCount/lcm).is_integer() == False: # check if Agentcount is appropriate to implement reshuffling
+                raise Exception("AgentCount must be a multiple of " + str( lcm))
+
+
+            def check_and_convert_to_int(val):
+
+                if abs(round(val) - val) < 1e-6:
+                    return abs(round(val))
+
+            if self.perf_reshuffle == True:# when true, permanent and transitory shocks are evenly distributed across both newborns and non newborns.
+                Min_AgentCount = check_and_convert_to_int(lcm/(1-self.LivPrb[0])) # total number of agents
+                if (self.AgentCount/Min_AgentCount).is_integer() == False: # check if Agentcount is appropriate to implement perfect reshuffling
+                    raise Exception("AgentCount must be a multiple of " +str(Min_AgentCount))
+
         for t in range(self.T_cycle):
             these = t == self.t_cycle
 
             # temporary, see #1022
             if self.cycles == 1:
                 t = t - 1
-
+
+
+            if self.reshuffle == True:
+                not_newborn = self.t_age != 0
+                these = np.logical_and(not_newborn, these) 
+
             N = np.sum(these)
+
             if N > 0:
                 # set current income distribution
                 IncShkDstnNow = self.IncShkDstn[t]
                 # and permanent growth factor
                 PermGroFacNow = self.PermGroFac[t]
                 # Get random draws of income shocks from the discrete distribution
-                IncShks = IncShkDstnNow.draw(N)
-
+                IncShks = IncShkDstnNow.draw(N , exact_match = self.reshuffle)
+           
                 PermShkNow[these] = (
                     IncShks[0, :] * PermGroFacNow
                 )  # permanent "shock" includes expected growth
                 TranShkNow[these] = IncShks[1, :]
-
+                
         # That procedure used the *last* period in the sequence for newborns, but that's not right
         # Redraw shocks for newborns, using the *first* period in the sequence.  Approximation.
         N = np.sum(newborn)
@@ -2276,23 +2338,29 @@ def get_shocks(self):
             # set current income distribution
             IncShkDstnNow = self.IncShkDstn[0]
             PermGroFacNow = self.PermGroFac[0]  # and permanent growth factor
-
+            
             # Get random draws of income shocks from the discrete distribution
-            EventDraws = IncShkDstnNow.draw_events(N)
+            EventDraws = IncShkDstnNow.draw_events(N,exact_match = self.reshuffle)
+            #EventDraws = IncShkDstnNow.draw(N,exact_match = self.reshuffle)
+
             PermShkNow[these] = (
                 IncShkDstnNow.atoms[0][EventDraws] * PermGroFacNow
             )  # permanent "shock" includes expected growth
             TranShkNow[these] = IncShkDstnNow.atoms[1][EventDraws]
-        #        PermShkNow[newborn] = 1.0
+
+
+        # PermShkNow[newborn] = 1.0
         #  Whether Newborns have transitory shock. The default is False.
         if not NewbornTransShk:
             TranShkNow[newborn] = 1.0
-
+
+
         # Store the shocks in self
         self.EmpNow = np.ones(self.AgentCount, dtype=bool)
         self.EmpNow[TranShkNow == self.IncUnemp] = False
         self.shocks["PermShk"] = PermShkNow
         self.shocks["TranShk"] = TranShkNow
+
 
     def calc_bounding_values(self):
         """

HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py

@@ -922,3 +922,64 @@ def test_calc_jacobian(self):
         self.assertAlmostEqual(CJAC_Perm.T[30][29], -0.06120, places=HARK_PRECISION)
         self.assertAlmostEqual(CJAC_Perm.T[30][30], 0.05307, places=HARK_PRECISION)
         self.assertAlmostEqual(CJAC_Perm.T[30][31], 0.04674, places=HARK_PRECISION)
+
+#%%  Test Reshuffling Monte Carlo Simulation Methods
+
+
+class test_reshuffling_methods(unittest.TestCase):
+    def test_reshuffling(self):
+
+        dict_harmenberg['UnempPrb'] = .1
+        dict_harmenberg['T_sim'] = 500
+        dict_harmenberg['reshuffle'] = True
+        Agent = IndShockConsumerType(**dict_harmenberg)
+        Agent.track_vars = ['aNrm']
+        Agent.solve()
+
+        Agent.neutral_measure = True
+        Agent.update_income_process()
+
+        Agent.perf_reshuffle = False
+        Agent.initialize_sim()
+        Agent.simulate()
+
+        Agg_A = []
+        for t in range(Agent.T_sim):
+            A = np.mean(Agent.history['aNrm'][t])
+
+            Agg_A.append(A)
+
+        Agg_A = np.array(Agg_A) 
+
+
+        self.assertAlmostEqual(np.var(Agg_A[100:]),2.338245625592009e-06, places=HARK_PRECISION)
+
+    def test_perf_reshuffling(self):
+
+        dict_harmenberg['UnempPrb'] = .1
+        dict_harmenberg['T_sim'] = 500
+        dict_harmenberg['reshuffle'] = True
+
+        Agent = IndShockConsumerType(**dict_harmenberg)
+        Agent.track_vars = ['aNrm']
+        Agent.solve()
+
+        Agent.perf_reshuffle = True
+        Agent.AgentCount  = 40000
+
+
+        Agent.initialize_sim()
+        Agent.simulate()
+
+        Agg_A = []
+        for t in range(Agent.T_sim):
+            A = np.mean(Agent.history['aNrm'][t])
+
+            Agg_A.append(A)
+
+        Agg_A = np.array(Agg_A) 
+
+
+        self.assertAlmostEqual(np.var(Agg_A[100:]),3.144866247778965e-08, places=HARK_PRECISION)
+
+

HARK/distribution.py

@@ -814,18 +814,39 @@ def dim(self) -> int:
         """
         return self.atoms.shape[:-1]
 
-    def draw_events(self, N: int) -> np.ndarray:
+    def draw_events(self, N: int,
+                    exact_match: bool = False) -> np.ndarray:
         """
         Draws N 'events' from the distribution PMF.
         These events are indices into atoms.
         """
-        # Generate a cumulative distribution
-        base_draws = self._rng.uniform(size=N)
-        cum_dist = np.cumsum(self.pmv)
+
+        if exact_match == True :            
+
+            events = np.arange(self.pmv.size)  # just a list of integers
+            cutoffs = np.round(np.cumsum(self.pmv) * N).astype(
+                int
+            )  # cutoff points between discrete outcomes
+            top = 0
 
-        # Convert the basic uniform draws into discrete draws
-        indices = cum_dist.searchsorted(base_draws)
+            # Make a list of event indices that closely matches the discrete distribution
+            event_list = []
+            for j in range(events.size):
+                bot = top
+                top = cutoffs[j]
+                event_list += (top - bot) * [events[j]]
+
+            # Randomly permute the event indices
+            indices = self._rng.permutation(event_list)
+
+        else:
+            # Generate a cumulative distribution
+            base_draws = self._rng.uniform(size=N)
+            cum_dist = np.cumsum(self.pmv)
 
+            # Convert the basic uniform draws into discrete draws
+            indices = cum_dist.searchsorted(base_draws)
+
         return indices
 
     def draw(