def updatepLvlGrid(self):
        '''
        Update the grid of persistent income levels.  Currently only works for
        infinite horizon models (cycles=0) and lifecycle models (cycles=1).  Not
        clear what to do about cycles>1 because the distribution of persistent
        income will be different within a period depending on how many cycles
        have elapsed.  This method uses a simulation approach to generate the
        pLvlGrid at each period of the cycle, drawing on the initial distribution
        of persistent income, the pLvlNextFuncs, and the attribute pLvlPctiles.

        Parameters
        ----------
        None

        Returns
        -------
        None
        '''
        orig_time = self.time_flow
        self.timeFwd()
        LivPrbAll = np.array(self.LivPrb)

        # Simulate the distribution of persistent income levels by t_cycle in a lifecycle model
        if self.cycles == 1:
            pLvlNow = drawLognormal(self.AgentCount,mu=self.pLvlInitMean,sigma=self.pLvlInitStd,seed=31382)
            pLvlGrid = [] # empty list of time-varying persistent income grids
            # Calculate distribution of persistent income in each period of lifecycle
            for t in range(len(self.PermShkStd)):
                if t > 0:
                    PermShkNow = drawDiscrete(N=self.AgentCount,P=self.PermShkDstn[t-1][0],X=self.PermShkDstn[t-1][1],exact_match=False,seed=t)
                    pLvlNow = self.pLvlNextFunc[t-1](pLvlNow)*PermShkNow
                pLvlGrid.append(getPercentiles(pLvlNow,percentiles=self.pLvlPctiles))

        # Calculate "stationary" distribution in infinite horizon (might vary across periods of cycle)
        elif self.cycles == 0:
            T_long = 1000 # Number of periods to simulate to get to "stationary" distribution
            pLvlNow = drawLognormal(self.AgentCount,mu=self.pLvlInitMean,sigma=self.pLvlInitStd,seed=31382)
            t_cycle = np.zeros(self.AgentCount,dtype=int)
            for t in range(T_long):
                LivPrb = LivPrbAll[t_cycle] # Determine who dies and replace them with newborns
                draws = drawUniform(self.AgentCount,seed=t)
                who_dies = draws > LivPrb
                pLvlNow[who_dies] = drawLognormal(np.sum(who_dies),mu=self.pLvlInitMean,sigma=self.pLvlInitStd,seed=t+92615)
                t_cycle[who_dies] = 0
                for j in range(self.T_cycle): # Update persistent income
                    these = t_cycle == j
                    PermShkTemp = drawDiscrete(N=np.sum(these),P=self.PermShkDstn[j][0],X=self.PermShkDstn[j][1],exact_match=False,seed=t+13*j)
                    pLvlNow[these] = self.pLvlNextFunc[j](pLvlNow[these])*PermShkTemp
                t_cycle = t_cycle + 1
                t_cycle[t_cycle == self.T_cycle] = 0

            # We now have a "long run stationary distribution", extract percentiles
            pLvlGrid = [] # empty list of time-varying persistent income grids
            for t in range(self.T_cycle):
                these = t_cycle == t
                pLvlGrid.append(getPercentiles(pLvlNow[these],percentiles=self.pLvlPctiles))

        # Throw an error if cycles>1
        else:
            assert False, "Can only handle cycles=0 or cycles=1!"

        # Store the result and add attribute to time_vary
        self.pLvlGrid = pLvlGrid
        self.addToTimeVary('pLvlGrid')
        if not orig_time:
            self.timeRev()
Update the grid of persistent income levels.  Currently only works for
        infinite horizon models (cycles=0) and lifecycle models (cycles=1).  Not
        clear what to do about cycles>1 because the distribution of persistent
        income will be different within a period depending on how many cycles
        have elapsed.  This method uses a simulation approach to generate the
        pLvlGrid at each period of the cycle, drawing on the initial distribution
        of persistent income, the pLvlNextFuncs, and the attribute pLvlPctiles.

        Parameters
        ----------
        None

        Returns
        -------
        None