Recitation 5: Simulating from the life-cycle model

Part 1

Review the details of the simple life-cycle savings model. Today we’re only going to be thinking about the income process, but for the sake of completeness I include below all of the code to parameterize and solve that model.

using Distributions,Random
using LinearAlgebra
Φ(x) = cdf(Normal(),x)

function tauchen(ρ,ση,Kϵ)
    sd = ση/sqrt(1-ρ^2)
    grid = range(-3sd,stop=3sd,length=Kϵ)
    Π = zeros(Kϵ,Kϵ)
    Δ = grid[2]-grid[1]
    for j=1:Kϵ
        Π[1,j] = Φ((grid[1] + Δ/2 - ρ*grid[j])/ση)
        Π[end,j] = 1 - Φ((grid[end] - Δ/2 - ρ*grid[j])/ση)
        for k=2:(Kϵ-1)
            Π[k,j] = Φ((grid[k] + Δ/2 - ρ*grid[j])/ση) - Φ((grid[k] - Δ/2 - ρ*grid[j])/ση)
        end
    end
    return Π,grid
end

u(c,σ) = c^(1-σ) / (1-σ)
function solve_max(V,t,iϵ,ia,pars)
    (;agrid,ϵgrid,Π,σ,Ka,r,β) = pars
    cash = exp(pars.μ[t] + ϵgrid[iϵ]) + agrid[ia]
    amax = 0
    vmax = -Inf
    loop = true
    a = 1
    while loop && a<Ka
        c = cash - agrid[a] / (1+r)
        if c>0
            #@views v = u(c,σ) + β * dot(Π[:,iϵ],V[:,a,t+1])
            v = u(c,σ)
            for iϵ′ in axes(V,1)
                v += β * Π[iϵ′,iϵ] * V[iϵ′,a,t+1]
            end
            if v>vmax
                vmax = v
                amax = a
            end
        else
            loop = false
        end
        a += 1 #<- move one up the grid space
    end
    return amax,vmax
end

function iterate!(V,A,t,pars)
    for ia in axes(V,2), iϵ in axes(V,1)
        A[iϵ,ia,t],V[iϵ,ia,t] = solve_max(V,t,iϵ,ia,pars)
    end
end
function terminal_values!(V,pars)
    (;σ,ψ,agrid) = pars
    for ia in axes(V,2), iϵ in axes(V,1)
        V[iϵ,ia] = ψ * u(agrid[ia],σ)
    end
end

function backward_induction!(V,A,pars)
    (;ψ,σ,T,agrid) = pars
    # set the values at T+1 (bequest motives)
    @views terminal_values!(V[:,:,T+1],pars)
    for t in reverse(1:T)
        iterate!(V,A,t,pars)
    end
end

pars = (;
    T = 45, β = 0.95, σ = 2,ρ = 0.9,ση = 0.1, μ = fill(2.,45), ψ = 5., r = 0.05
)
Ka = 100
Kϵ = 5
agrid = LinRange(0,pars.μ[1] * pars.T,Ka) #
Π,ϵgrid = tauchen(pars.ρ,pars.ση,Kϵ)
pars = (;pars...,Ka,agrid,Π,ϵgrid,Kϵ)
V = zeros(pars.Kϵ,pars.Ka,pars.T+1)
A = zeros(Int64,pars.Kϵ,pars.Ka,pars.T)
backward_induction!(V,A,pars)
@time backward_induction!(V,A,pars)

  0.011214 seconds

Part 2: The Stationary Distribution of \(\varepsilon\)

Suppose that \(\pi^*\) is the stationary distribution of \(\varepsilon\). It must obey:

\[ \pi^* = \Pi \pi^* \]

where the \(k\)th column of \(\Pi\) gives the transition probabilities given that the current value is \(\varepsilon_{k}\). This is equivalent to \(\pi^*\) being the eigenvector of \(\mathbf{I} - \Pi\) that corresponds to an eigenvalue of 0.

function stat_dist(Π)
    Kϵ = size(Π,1)
    vals, vec = eigen(I(Kϵ) .- Π)
    return vec[:,1] ./ sum(vec[:,1])
end

π_stat = stat_dist(Π)
pars = (;pars...,π_stat)

(T = 45, β = 0.95, σ = 2, ρ = 0.9, ση = 0.1, μ = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0  …  2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0], ψ = 5.0, r = 0.05, Ka = 100, agrid = LinRange{Float64}(0.0, 90.0, 100), Π = [0.8490507777857362 0.019473727871012682 … 7.346962855655732e-17 3.459030953951908e-30; 0.15094537665867613 0.8961919626850798 … 7.260018586910001e-7 1.237828285826989e-15; … ; 1.2212453270876722e-15 7.260018586308092e-7 … 0.8961919626850798 0.15094537665867613; 0.0 1.1102230246251565e-16 … 0.019473727871012647 0.8490507777857362], ϵgrid = -0.6882472016116855:0.34412360080584276:0.6882472016116855, Kϵ = 5, π_stat = [0.03046350803405268, 0.23613279404893608, 0.46680739583402275, 0.23613279404893592, 0.03046350803405253])

We can make sure we’re doing this correctly by checking:

π_stat .- Π * π_stat

5-element Vector{Float64}:
  0.0
  2.7755575615628914e-17
  0.0
  0.0
 -3.469446951953614e-18

Part 3: Simulating the Income Process

In principle, one can simulate draws from any distribution by simulating draws from the uniform distribution using rand() and using the inverse cdf. To see this, note that for any random variable:

\[ F_{X}(X) \sim \text{unif}[0,1] \]

and hence if \(U\sim\text{unif}[0,1]\):

\[ F_{X}^{-1}(U) \sim X \]

Here we will make use of the Distributions package and the Categorical distribution for simulation.

# this function simulates on path
function simulate_shocks!(ε,π_stat,Π)
    ε_current = rand(π_stat)
    for t in eachindex(ε)
        ε[t] = ε_current
        ε_current = rand(Π[ε_current])
    end
end

function simulate_shocks!(ε,pars)
    (;Π,π_stat) = pars
    π_stat_dist = Categorical(π_stat)
    Π_dist = [Categorical(Π[:,k]) for k in axes(Π,2)]
    for n in axes(ε,2)
        @views simulate_shocks!(ε[:,n],π_stat_dist,Π_dist)
    end
end

nsim = 500
ε_sim = zeros(Int64,pars.T,nsim)
simulate_shocks!(ε_sim,pars)

Part 4: Holding Simulation Error Fixed

Now suppose we plan to use the simulated draws of \(\varepsilon\) to calculate some moments to match to the data. For example, we could try to match the lagged covariance to lagged covariance in the data. In practice we won’t have to simulation to do this, but for the sake of argument, let’s pretend we do. We could write the function to calculate moments as:

function moments_ε(ε,pars)
    simulate_shocks!(ε,pars)
    ε_value = pars.ϵgrid[ε]
    return [var(ε_value[1,:]) ; cov(ε_value[2,:],ε_value[1,:]) ; cov(ε_value[3,:],ε_value[1,:])]
end

moments_ε (generic function with 1 method)

But notice what happens when we call this function multiple times:

@show moments_ε(ε_sim,pars);
@show moments_ε(ε_sim,pars);
@show moments_ε(ε_sim,pars);

moments_ε(ε_sim, pars) = [0.080886562598882, 0.07397874696761952, 0.06613875118658376]
moments_ε(ε_sim, pars) = [0.08498739584432013, 0.07629638223816074, 0.07038339837569878]
moments_ε(ε_sim, pars) = [0.08589157261892205, 0.08071616918046617, 0.0733280244699928]

The moments keep jumping around! This is going to be a problem if we want to match the simulated moments to ones we have calculated in the data. The simulation error is moving the goalposts around in ways that we can’t predict. The solution is to fix simulation error by fixing the seed of our pseudorandom number generator. To see how, let’s rewrite the function:

function moments_ε(ε,pars ; seed = 202404)
    Random.seed!(seed)
    simulate_shocks!(ε,pars)
    ε_value = pars.ϵgrid[ε]
    return [var(ε_value[1,:]) ; cov(ε_value[2,:],ε_value[1,:]) ; cov(ε_value[3,:],ε_value[1,:])]
end

@show moments_ε(ε_sim,pars);
@show moments_ε(ε_sim,pars);
@show moments_ε(ε_sim,pars);

@show moments_ε(ε_sim,pars ; seed = 12345);
@show moments_ε(ε_sim,pars ; seed = 12345);
@show moments_ε(ε_sim,pars ; seed = 12345);

moments_ε(ε_sim, pars) = [0.07921632739162526, 0.0694721020989348, 0.06356291530429277]
moments_ε(ε_sim, pars) = [0.07921632739162526, 0.0694721020989348, 0.06356291530429277]
moments_ε(ε_sim, pars) = [0.07921632739162526, 0.0694721020989348, 0.06356291530429277]
moments_ε(ε_sim, pars; seed = 12345) = [0.08676489821748767, 0.08202520831135969, 0.07516960236262014]
moments_ε(ε_sim, pars; seed = 12345) = [0.08676489821748767, 0.08202520831135969, 0.07516960236262014]
moments_ε(ε_sim, pars; seed = 12345) = [0.08676489821748767, 0.08202520831135969, 0.07516960236262014]

This is something that is very important to keep in mind when you are using the simulated method of moments to estimate models.