def calcValueAtBirth(cLvlHist,BirthBool,PlvlHist,MrkvHist,DiscFac,CRRA):
    '''
    Calculate expected value of being born in each Markov state using the realizations
    of consumption for a history of many consumers.  The histories should already be
    trimmed of the "burn in" periods.

    Parameters
    ----------
    cLvlHist : np.array
        TxN array of consumption level history for many agents across many periods.
        Agents who die are replaced by newborms.
    BirthBool : np.array
        TxN boolean array indicating when agents are born, replacing one who died.
    PlvlHist : np.array
        T length vector of aggregate permanent productivity levels.
    MrkvHist : np.array
        T length vector of integers for the Markov index in each period.
    DiscFac : float
        Intertemporal discount factor.
    CRRA : float
        Coefficient of relative risk aversion.

    Returns
    -------
    vAtBirth : np.array
        J length vector of average lifetime value at birth by Markov state.
    '''
    J = np.max(MrkvHist) + 1 # Number of Markov states
    T = MrkvHist.size        # Length of simulation
    I = cLvlHist.shape[1]    # Number of agent indices in histories
    u = lambda c : CRRAutility(c,gam=CRRA)

    # Initialize an array to hold each agent's lifetime utility
    BirthsByPeriod = np.sum(BirthBool,axis=1)
    BirthsByState = np.zeros(J,dtype=int)
    for j in range(J):
        these = MrkvHist == j
        BirthsByState[j] = np.sum(BirthsByPeriod[these])
    N = np.max(BirthsByState) # Array must hold this many agents per row at least
    vArray = np.zeros((J,N)) + np.nan
    n = np.zeros(J,dtype=int)

    # Loop through each agent index
    DiscVec = DiscFac**np.arange(T)
    for i in range(I):
        birth_t = np.where(BirthBool[:,i])[0]
        # Loop through each agent who lived and died in this index
        for k in range(birth_t.size-1): # Last birth event has no death, so ignore
            # Get lifespan of this agent and circumstances at birth
            t0 = birth_t[k]
            t1 = birth_t[k+1]
            span = t1-t0
            j = MrkvHist[t0]
            # Calculate discounted flow of utility for this agent and store it
            cVec = cLvlHist[t0:t1,i]/PlvlHist[t0]
            uVec = u(cVec)
            v = np.dot(DiscVec[:span],uVec)
            vArray[j,n[j]] = v
            n[j] += 1

    # Calculate expected value at birth by state and return it
    vAtBirth = np.nanmean(vArray,axis=1)
    return vAtBirth
Calculate expected value of being born in each Markov state using the realizations
    of consumption for a history of many consumers.  The histories should already be
    trimmed of the "burn in" periods.

    Parameters
    ----------
    cLvlHist : np.array
        TxN array of consumption level history for many agents across many periods.
        Agents who die are replaced by newborms.
    BirthBool : np.array
        TxN boolean array indicating when agents are born, replacing one who died.
    PlvlHist : np.array
        T length vector of aggregate permanent productivity levels.
    MrkvHist : np.array
        T length vector of integers for the Markov index in each period.
    DiscFac : float
        Intertemporal discount factor.
    CRRA : float
        Coefficient of relative risk aversion.

    Returns
    -------
    vAtBirth : np.array
        J length vector of average lifetime value at birth by Markov state.