diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index 6c596b4f9..ff03fcc53 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -2496,7 +2496,8 @@ def constructLognormalIncomeProcessUnemployment(self):
                 ShkPrbsRet      = np.array([1.0])
             IncomeDstnRet = DiscreteDistribution(ShkPrbsRet,
                                                  [PermShkValsRet,
-                                                  TranShkValsRet])
+                                                  TranShkValsRet],
+                                                 seed=self.RNG.randint(0, 2**31-1))
 
         # Loop to fill in the list of IncomeDstn random variables.
         for t in range(T_cycle): # Iterate over all periods, counting forward
@@ -2518,7 +2519,9 @@ def constructLognormalIncomeProcessUnemployment(self):
                 ).approx(PermShkCount, tail_N=0)
                 ### REPLACE
                 ###REPLACE
-                IncomeDstn.append(combineIndepDstns(PermShkDstn_t,TranShkDstn_t)) # mix the independent distributions
+                IncomeDstn.append(combineIndepDstns(PermShkDstn_t,
+                                                    TranShkDstn_t,
+                                                    seed = self.RNG.randint(0, 2**31-1))) # mix the independent distributions
                 PermShkDstn.append(PermShkDstn_t)
                 TranShkDstn.append(TranShkDstn_t)
         return IncomeDstn, PermShkDstn, TranShkDstn
diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py
index 934a0e8ef..27c635932 100644
--- a/HARK/ConsumptionSaving/ConsPortfolioModel.py
+++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py
@@ -353,7 +353,7 @@ def getRisky(self):
 
         mu = np.log(RiskyAvg / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd)))
         sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd))
-        self.RiskyNow = Lognormal(mu, sigma).draw(1, seed=self.RNG.randint(0, 2**31-1))
+        self.RiskyNow = Lognormal(mu, sigma, seed=self.RNG.randint(0, 2**31-1)).draw(1)
         
         
     def getAdjust(self):
@@ -370,7 +370,7 @@ def getAdjust(self):
         -------
         None
         '''
-        self.AdjustNow = Bernoulli(self.AdjustPrb).draw(self.AgentCount, seed=self.RNG.randint(0, 2**31-1))
+        self.AdjustNow = Bernoulli(self.AdjustPrb, seed=self.RNG.randint(0, 2**31-1)).draw(self.AgentCount)
         
         
     def getRfree(self):
diff --git a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
index 374803b10..1c219047e 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
@@ -93,4 +93,4 @@ def test_simulation(self):
 
         self.assertAlmostEqual(
             np.mean(self.agent.history['mLvlNow']),
-            1.2063544810887799)
+            1.2043946738813716)
diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
index 9353cba49..8c046f193 100644
--- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
+++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
@@ -23,11 +23,11 @@ def test_getShocks(self):
         self.agent.getShocks()
 
         self.assertEqual(self.agent.PermShkNow[0],
-                         0.9686755529883603)
+                         1.0427376294215103)
         self.assertEqual(self.agent.PermShkNow[1],
-                         1.0050166461586711)
+                         0.9278094171517413)
         self.assertEqual(self.agent.TranShkNow[0],
-                         1.2093790234554653)
+                         0.881761797501595)
 
     def test_ConsIndShockSolverBasic(self):
         LifecycleExample = IndShockConsumerType(
@@ -102,13 +102,11 @@ def test_simulated_values(self):
         self.agent.initializeSim()
         self.agent.simulate()
 
-        print(self.agent.aLvlNow)
-
         self.assertAlmostEqual(self.agent.MPCnow[1],
                                0.5711503906043797)
 
         self.assertAlmostEqual(self.agent.aLvlNow[1],
-                               0.21251164208206255)
+                               0.18438326264597635)
 
 
 class testBufferStock(unittest.TestCase):
diff --git a/HARK/distribution.py b/HARK/distribution.py
index 542ac559a..c8e2b8b4e 100644
--- a/HARK/distribution.py
+++ b/HARK/distribution.py
@@ -824,7 +824,7 @@ def addDiscreteOutcome(distribution, x, p, sort = False):
 
     return DiscreteDistribution(pmf,X)
 
-def combineIndepDstns(*distributions):
+def combineIndepDstns(*distributions, seed=0):
     '''
     Given n lists (or tuples) whose elements represent n independent, discrete
     probability spaces (probabilities and values), construct a joint pmf over
@@ -897,4 +897,4 @@ def combineIndepDstns(*distributions):
     P_out = np.prod(np.array(P_temp),axis=0)
 
     assert np.isclose(np.sum(P_out),1),'Probabilities do not sum to 1!'
-    return DiscreteDistribution(P_out, X_out)
+    return DiscreteDistribution(P_out, X_out, seed = seed)