Replicating Maliar, Maliar, and Valli (2010, JEDC) to solve Krusell and Smith (1998, JPE) model using Julia

By Shunsuke Hori

Overview of the notebook

This notebook solves the model of Krusell and Smith (1998, JPE) and succesfully replicating the result of Maliar, Maliar, and Valli (2010, JEDC)

The solution strategy is as follows

  1. Solve individual problem by Euler equation method with 2D interpolation
    • Agents are boundedly rational. In the code, they take into account the information about the mean of capital
    • Aggregate law of motion is approximated by log-linear relation, i.e. $\log(K_{t+1})=B1+B2\log(K_{t})$ for good aggregate state and $\log(K_{t+1})=B3+B4\log(K_{t})$ for bad aggregate state
    • If specified, Howard's policy iteration is used
  2. Compute the path of aggregate capital using the policy function obtained by VFI
    • Simulation is based on Monte Carlo. That is, aggregate technology shocks and idiosyncratic employment shocks are drawn for many agents and many periods. Then, using the LLN, the aggregate capital is computed by aggregating all agents for all period.
  3. Update the coefficient of aggregate capital law of motion, $B1$, $B2$, $B3$ and $B4$, by regression
  4. Check convergence of $B1$, $B2$, $B3$ and $B4$

Additionally, this notebook also includes the result with value function iteration as a solution method for individual utility maximization problem.

NOTE: Regarding interpolation, Krusell and Smith uses various interpolation scheme depending on the purpose, including polynomial interpolation. Maliar, Maliar, and Valli uses spline interpolation in their paper. This notebook only uses linear interpolation because I could not find interpolation package for polynomial and spline 2D-interpolation.

Code to solve models

First thing to do is import some packages

In [1]:
using Interpolations # to use interpolation 
using QuantEcon  # to use a function, `gridmake`
using Plots      # to plot the result
pyplot()
#plotlyjs()
using Optim      # to use minimization routine to maximize RHS of bellman equation
using GLM        # to regress
using JLD        # to save the result

Model Setup

Functions in this cell are prepared for model parameters and initial guess of the solutions

  • Types
    • TransitionMatrix: collection of transition matrix
    • KSParameter: collection of model parameters, functional forms, and grids
    • KSSolution: collection of solution, which is guess at first
  • Methods
    • create_transition_matrix: construct an instance of type TransitionMatrix
    • KSParameter: construct an instance of type KSParameter
    • place_polynominal_grid: create polynominal grid
    • r: compute interest rate ( = marginal productivity of capital)
    • w: compute wage rate ( = marginal productivity of labor)
    • KSSolution: construct an instance of type KSSolution
      • Guess can be the load of previous result
In [2]:
"""
Collection of transition matrix
"""
immutable TransitionMatrix
    P::Array{Float64,2}       # 4x4 
    Pz::Array{Float64,2}      # 2x2 aggregate shock
    Peps_gg::Array{Float64,2} # 2x2 idiosyncratic shock conditional on good to good
    Peps_bb::Array{Float64,2} # 2x2 idiosyncratic shock conditional on bad to bad
    Peps_gb::Array{Float64,2} # 2x2 idiosyncratic shock conditional on good to bad
    Peps_bg::Array{Float64,2} # 2x2 idiosyncratic shock conditional on bad to good
end

"""
Collection of model parameters
"""
immutable KSParameter
    u::Function
    beta::Float64
    alpha::Float64
    delta::Float64
    theta::Float64
    l_bar::Float64
    k_grid::Vector{Float64}
    K_grid::Vector{Float64}
    z_grid::Vector{Float64}
    eps_grid::Vector{Float64}
    s_grid::Array{Float64,2}
    k_size::Int64
    K_size::Int64
    z_size::Int64
    eps_size::Int64
    s_size::Int64
    ug::Float64
    ub::Float64
    TransMat::TransitionMatrix # bunch of transition matrix
    mu::Float64
end

""" 
Create transition matrices for aggregate shock,
idiosyncratic shock, and shock state

##### Arguments
- `ug` : unemployment rate in good state
- `ub` : unemployment rate in bad state
- `zg_ave_dur` : average duration of good state
- `zb_ave_dur` : average duration of bad state
- `ug_ave_dur` : average duration of unemployment in good state
- `ub_ave_dur` : average duration of unemployment in bad state
- `puu_rel_gb2bb` : prob. of u to u cond. on g to b relative to that of b to b 
- `puu_rel_bg2gg` : prob. of u to u cond. on b to g relative to that of g to g
"""
function create_transition_matrix(ug::Float64,ub::Float64,
        zg_ave_dur::Float64,zb_ave_dur::Float64,
        ug_ave_dur::Float64,ub_ave_dur::Float64,
        puu_rel_gb2bb::Float64,puu_rel_bg2gg::Float64)
    
    # probability of remaining in good state
    pgg = 1-1/zg_ave_dur
    # probability of remaining in bad state
    pbb = 1-1/zb_ave_dur
    # probability of changing from g to b
    pgb = 1-pgg
    # probability of changing from b to g
    pbg = 1-pbb  
    
    # prob. of 0 to 0 cond. on g to g
    p00_gg = 1-1/ug_ave_dur
    # prob. of 0 to 0 cond. on b to b
    p00_bb = 1-1/ub_ave_dur
    # prob. of 0 to 1 cond. on g to g
    p01_gg = 1-p00_gg
    # prob. of 0 to 1 cond. on b to b
    p01_bb = 1-p00_bb
    
    # prob. of 0 to 0 cond. on g to b
    p00_gb=puu_rel_gb2bb*p00_bb
    # prob. of 0 to 0 cond. on b to g
    p00_bg=puu_rel_bg2gg*p00_gg
    # prob. of 0 to 1 cond. on g to b
    p01_gb=1-p00_gb
    # prob. of 0 to 1 cond. on b to g
    p01_bg=1-p00_bg
    
    # prob. of 1 to 0 cond. on  g to g
    p10_gg=(ug - ug*p00_gg)/(1-ug)
    # prob. of 1 to 0 cond. on b to b
    p10_bb=(ub - ub*p00_bb)/(1-ub)
    # prob. of 1 to 0 cond. on g to b
    p10_gb=(ub - ug*p00_gb)/(1-ug)
    # prob. of 1 to 0 cond on b to g
    p10_bg=(ug - ub*p00_bg)/(1-ub)
    # prob. of 1 to 1 cond. on  g to g
    p11_gg= 1-p10_gg
    # prob. of 1 to 1 cond. on b to b
    p11_bb= 1-p10_bb
    # prob. of 1 to 1 cond. on g to b
    p11_gb= 1-p10_gb
    # prob. of 1 to 1 cond on b to g
    p11_bg= 1-p10_bg
    
    #   (g1)         (b1)        (g0)       (b0)
    P=[pgg*p11_gg pgb*p11_gb pgg*p10_gg pgb*p10_gb;
        pbg*p11_bg pbb*p11_bb pbg*p10_bg pbb*p10_bb;
        pgg*p01_gg pgb*p01_gb pgg*p00_gg pgb*p00_gb;
        pbg*p01_bg pbb*p01_bb pbg*p00_bg pbb*p00_bb
        ]
    Pz=[pgg pgb;
        pbg pbb]
    Peps_gg=[p11_gg p10_gg
              p01_gg p00_gg]
    Peps_bb=[p11_bb p10_bb
              p01_bb p00_bb]
    Peps_gb=[p11_gb p10_gb
              p01_gb p00_gb]
    Peps_bg=[p11_bg p10_bg
              p01_bg p00_bg]
    TransMat=TransitionMatrix(P,Pz,Peps_gg,Peps_bb,Peps_gb,Peps_bg)
    return TransMat
end

"""
Creates KSParameter instance
"""
function KSParameter(;
            beta::Float64=0.99,
            alpha::Float64=0.36,
            delta::Float64=0.025,
            theta::Float64=1.0,
            k_min::Float64=1e-16,
            k_max::Float64=1000.0,
            k_size::Int64=100,
            K_min::Float64=30.0,
            K_max::Float64=50.0,
            K_size::Int64=4,
            z_min::Float64=0.99,
            z_max::Float64=1.01,
            z_size::Int64=2,
            eps_min::Float64=0.0,
            eps_max::Float64=1.0,
            eps_size::Int64=2,
            ug::Float64=0.04,
            ub::Float64=0.1,
            zg_ave_dur::Float64=8.0,
            zb_ave_dur::Float64=8.0,
            ug_ave_dur::Float64=1.5,
            ub_ave_dur::Float64=2.5,
            puu_rel_gb2bb::Float64=1.25,
            puu_rel_bg2gg::Float64=0.75,
            mu::Float64=0.0
            )
    if theta == 1.0
        u = (c) -> log(c)
    else
        u = (c) -> (c^(1.0-theta)-1.0)/(1.0-theta)
    end
    l_bar=1/(1-ub)
    k_grid=place_polynominal_grid(k_min,k_max,k_size,degree=7.0)   # individual capital grid
    K_grid=collect(linspace(K_min,K_max,K_size))   # aggregate capital grid
    z_grid=collect(linspace(z_max,z_min,z_size))   # aggregate technology shock
    eps_grid=collect(linspace(eps_max,eps_min,eps_size))  # idiosyncratic employment shock
    s_grid=gridmake(z_grid,eps_grid)               # shock grid
    # collection of transition matrices
    TransMat=create_transition_matrix(ug,ub,
        zg_ave_dur,zb_ave_dur,
        ug_ave_dur,ub_ave_dur,
        puu_rel_gb2bb,puu_rel_bg2gg)

    ksp=KSParameter(u,beta,alpha,delta,theta,l_bar,k_grid,K_grid,z_grid,eps_grid,s_grid,
    k_size,K_size,z_size,eps_size,z_size*eps_size,ug,ub,TransMat,mu)

    return ksp
end

function place_polynominal_grid(k_min::Float64,k_max::Float64,k_size::Int64;degree::Float64=0.7)
    grid=Array{Float64}(k_size)
    for (i, x) in enumerate(linspace(0,0.5,k_size))
        grid[i]=(x/0.5)^degree*(k_max-k_min)+k_min
    end
    return grid
end
"""
Compute interest rate given aggregate capital, labor, and productivity

##### Arguments
- `alpha` : capital share
- `z` : aggregate shock
- `K` : aggregate capital
- `L` : aggregate labor
"""
r(alpha::Float64,z::Float64,K::Float64,L::Float64)=alpha*z*K^(alpha-1)*L^(1-alpha)

"""
Compute wage given aggregate capital, labor, and productivity

##### Arguments
- `alpha` : capital share
- `z` : aggregate shock
- `K` : aggregate capital
- `L` : aggregate labor
"""
w(alpha::Float64,z::Float64,K::Float64,L::Float64)=(1-alpha)*z*K^(alpha)*L^(-alpha)    

"""
Collection of KS solution
"""
type KSSolution
    k_opt::Array{Float64,3}
    value::Array{Float64,3}
    B::Vector{Float64}
    R2::Vector{Float64}
end

"""
Create KSSolution instance
"""
function KSSolution(
        ksp::KSParameter;
        load_value::Bool=false,
        load_B::Bool=false,
        filename::String="result.jld")
    if load_value || load_B
        result=load(filename)
        kss_temp=result["kss"]
    end
    if load_value
        k_opt=kss_temp.k_opt
        value=kss_temp.value
    else
        k_opt=ksp.beta*repeat(ksp.k_grid,outer=[1,ksp.K_size,ksp.s_size])
        k_opt=0.9*repeat(ksp.k_grid,outer=[1,ksp.K_size,ksp.s_size])
        value=ksp.u.(0.1/0.9*k_opt)/(1-ksp.beta)
    end
    if load_B
        B=kss_temp.B
    else
        B=[0.0, 1.0, 0.0, 1.0]
    end
    kss=KSSolution(k_opt,value,B,[0.0,0.0])
    return kss
end
Out[2]:
KSSolution

Shock Generation

The functions in this cell are used to draw aggregate and idiosyncratic shocks

  • Methods
    • get_shock: translate an index shock combination grid to shock value
    • generate_shocks: draw aggregate and idiosyncratic shocks, which is the main function in this cell
      • draw_eps_shock (local function inside generate_shocks): draw idiosyncratic shocks given aggregate shocks
In [3]:
""" 
Translate shock index into shock value

##### Arguments
- `s_grid` : shock  grid
- `s_i` : shock index
"""
get_shock(s_grid::Array{Float64,2},s_i::Integer) = 
    s_grid[s_i,1], s_grid[s_i,2] 

"""
Generate aggregate and idiosyncratic shock

##### Arguments
- `ksp` : instance of KSParameter type
- `z_shock_size` : size of aggregate shock
- `population` : size idiosyncratic shock in one period
"""
function generate_shocks(ksp::KSParameter;
        z_shock_size::Int64=11000,population::Int64=10000)
    
    
    """ 
    Draw idiosyncratic shock given previous idiosyncratic shock and 
    transition matrix.
    The transition matrix must be consistent with aggregate shock

    ##### Arguments
    - `eps_shock` : preallocated vector. current period shock is stored in it
    - `eps_shock_before` : previous period idiosyncratic shock
    - `Peps` : transition matrix of the current period
    """
    function draw_eps_shock(eps_shock_before::Vector{Float64},
                              Peps::Array{Float64,2})
    
        eps_shocks = similar(eps_shock_before)
        # loop over entire population
        for i=1:length(eps_shocks)
            rand_draw=rand()
            eps_shocks[i]=ifelse(eps_shock_before[i]==1.0,
                Float64(Peps[1,1]>rand_draw),  # if employed before
                Float64(Peps[2,1]>rand_draw))  # if unemployed before
        end
        return eps_shocks
    end
    
    # unpack parameters
    Peps_gg=ksp.TransMat.Peps_gg
    Peps_bg=ksp.TransMat.Peps_bg
    Peps_gb=ksp.TransMat.Peps_gb
    Peps_bb=ksp.TransMat.Peps_bb
    
    # draw aggregate shock
    z_shock=simulate(MarkovChain(ksp.TransMat.Pz,ksp.z_grid),z_shock_size)
    
    ### Let's draw individual shock ###
    eps_shock=Array{Float64}(z_shock_size,population) # preallocation
    
    # first period
    rand_draw=rand(population)
    if z_shock[1]==ksp.z_grid[1] # if good
        eps_shock[1,:].=Int64.(rand_draw.>ksp.ug) # if draw is higher, become employed
    elseif z_shock[1]==ksp.z_grid[2] # if bad
        eps_shock[1,:].=Int64.(rand_draw.>ksp.ub) # if draw is higher, become employed
    else
        error("the value of z_shock[1] (=$(z_shock[1])) is strange")
    end
    
    # from second period ...
    for t=2:z_shock_size
        if z_shock[t]==ksp.z_grid[1] && z_shock[t-1]==ksp.z_grid[1]  # if g to g
            eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_gg)
        elseif z_shock[t]==ksp.z_grid[1] && z_shock[t-1]==ksp.z_grid[2]  # if b to g
            eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_bg)
        elseif z_shock[t]==ksp.z_grid[2] && z_shock[t-1]==ksp.z_grid[1]  # if g to b
            eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_gb)
        elseif z_shock[t]==ksp.z_grid[2] && z_shock[t-1]==ksp.z_grid[2]  # if b to b
            eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_bb)
        else
            error("the value of z_shock[t] (=$(z_shock[t])) is strange")
        end
    end
    
    # adjustment
    for t=1:z_shock_size
        n_e=Int64(sum(eps_shock[t,:]))
        er_ideal = ifelse(z_shock[t] == ksp.z_grid[1],
                    1.0-ksp.ug, 1.0-ksp.ub)
        gap = Int64(round(er_ideal*population)) - n_e
        if gap > 0
            for j=1:gap
                become_employed_i = 
                    sample(find(x-> isapprox(x,ksp.eps_grid[2]), eps_shock[t,:]))
                eps_shock[t,become_employed_i] = ksp.eps_grid[1]
            end
        elseif gap < 0
            for j=1:(-gap)
                become_unemployed_i = 
                    sample(find(x-> isapprox(x,ksp.eps_grid[1]), eps_shock[t,:]))
                eps_shock[t,become_unemployed_i] = ksp.eps_grid[2]
            end
        end 
    end
            
    return z_shock, eps_shock    
end
Out[3]:
generate_shocks

Functions for Euler equation method

Functions in this cell is used to solve individual problem by Euler equation method

  • Methods

    • compute_Kp_L: compute approximated $K'$ and $L$, which depend on current aggregate shock
    • compute_expectation_FOC: compute expectation term in Euler equation
    • euler_method!: find optimal policy by Euler equation methods, namely, iteration of Euler equation

      $$ \left(c\right)^{-\theta}=\beta E\left[\left(c'\right)^{-\theta}(1-\delta+r')\right] $$

In [4]:
"""
Compute next period aggregate capital and labor

##### Arguments
- `K` : Current aggregate capital
- `s_i` : current shock index
- `B` : coefficient of ALM for capital
"""
function compute_Kp_L(K::AbstractFloat,s_i::Integer,B::Vector{Float64},ksp::KSParameter)
    Kp,L=ifelse(ksp.s_grid[s_i,1]==ksp.z_grid[1],
            (exp(B[1]+B[2]*log(K)),ksp.l_bar*(1-ksp.ug)),
            (exp(B[3]+B[4]*log(K)),ksp.l_bar*(1-ksp.ub)))
    Kp = ifelse(Kp<ksp.K_grid[1],ksp.K_grid[1],Kp)
    Kp = ifelse(Kp>ksp.K_grid[end],ksp.K_grid[end],Kp)
    return Kp, L
end

"""
Compute expectation term in Euler equation

##### Arguments
- `kp` : next period individual capital holding
- `Kp` : next period aggregate capital
- `s_i` : current state of shock
- `ksp` : KSParameter instance
"""
function compute_expectation_FOC(kp::Float64,
                                  Kp::Float64,
                                  s_i::Int64,
                                  ksp::KSParameter)
    alpha, theta, delta, l_bar, mu, P =
        ksp.alpha, ksp.theta, ksp.delta, ksp.l_bar, ksp.mu, ksp.TransMat.P
    expec = 0.0
    for s_n_i = 1:ksp.s_size
        zp, epsp = ksp.s_grid[s_n_i,1], ksp.s_grid[s_n_i,2]
        Kpp, Lp = compute_Kp_L(Kp,s_n_i,kss.B,ksp)
        rn=r(alpha,zp,Kp,Lp)
        kpp=interpolate((ksp.k_grid,ksp.K_grid),
                kss.k_opt[:,:,s_n_i],Gridded(Linear()))
        cp = (rn+1-delta)*kp+
                   w(alpha,zp,Kp,Lp)*(epsp*l_bar+mu*(1.0-epsp))-kpp[kp,Kp]
        expec = expec + P[s_i,s_n_i]*(cp)^(-theta)*(1-delta+rn)
    end 
    return expec
end

"""
Solve individual problem by Euler equation method

##### Arguments
- `ksp` : KSParameter instance
- `kss` : KSSolution instance
- `n_iter` : maximum number of iteration of Euler equation method
- `tol` : tolerance of policy function convergence
- `update_k` : degree of update of policy function
"""
function euler_method!(ksp::KSParameter,
                        kss::KSSolution;
                        n_iter::Integer=30000,
                        tol::AbstractFloat=1e-8,
                        update_k::AbstractFloat=1e-8)
    println("solving individual problem by Euler equation method")
    alpha, beta, delta, theta, l_bar, mu = 
        ksp.alpha, ksp.beta, ksp.delta, ksp.theta, ksp.l_bar, ksp.mu
    k_grid, k_size = ksp.k_grid, ksp.k_size
    K_grid, K_size = ksp.K_grid, ksp.K_size
    s_grid, s_size = ksp.s_grid, ksp.s_size
    k_min, k_max = minimum(k_grid), maximum(k_grid)
    counter=0
    k_opt_n=similar(kss.k_opt)
    while true
        counter += 1
        dif_k=0.0
        for s_i = 1:s_size
            z, eps = s_grid[s_i,1], s_grid[s_i,2]
            for K_i = 1:K_size
                K = K_grid[K_i]
                Kp, L = compute_Kp_L(K,s_i,kss.B,ksp)
                for k_i = 1:k_size
                    k=k_grid[k_i]
                    wealth = (r(alpha,z,K,L)+1-delta)*k+
                                w(alpha,z,K,L)*(eps*l_bar+mu*(1.0-eps))
                    expec=compute_expectation_FOC(kss.k_opt[k_i,K_i,s_i],Kp,s_i,ksp)
                    cn=(beta*expec)^(-1.0/theta)
                    k_opt_n[k_i,K_i,s_i] = wealth-cn
                    k_opt_n[k_i,K_i,s_i] = ifelse(k_opt_n[k_i,K_i,s_i]>k_max,k_max,k_opt_n[k_i,K_i,s_i])
                    k_opt_n[k_i,K_i,s_i] = ifelse(k_opt_n[k_i,K_i,s_i]<k_min,k_min,k_opt_n[k_i,K_i,s_i])
                    
                end
            end
        end
        dif_k=maxabs(k_opt_n-kss.k_opt)
        kss.k_opt.=update_k.*k_opt_n.+(1-update_k).*kss.k_opt
        if dif_k<tol
            println("Euler method converged with $counter iterations")
            break
        end
        if counter >=n_iter
            println("Euler method failed to converge with $counter iterations (dif = $dif_k)")
            break
        end
    end
end
Out[4]:
euler_method!

Simulate path of aggregate capital

  • Mathods
    • simulate_aggregate_path!: simulating aggregate path of capital using shocks drawn
In [5]:
"""
Simulate aggregate capital's path using policy function 
and generated aggregate and idiosyncratic shock

##### Arguments
- `ksp` : KSParameter instance
- `z_shocks` : aggregate shocks
- `eps_shocks` : idiosyncratic shocks
- `k_population` : initial capital holding of all agents
"""
function simulate_aggregate_path!(ksp::KSParameter,kss::KSSolution,
        z_shocks::Vector{Float64},eps_shocks::Array{Float64,2},
        k_population::Vector{Float64},K_ts::Vector{Float64})
    
    println("simulating aggregate path ... please wait ... ")
    
    T=length(z_shocks)   # simulated duration
    N=size(eps_shocks,2) # number of agents
    
    
    # loop over T periods
    for (t,z) = enumerate(z_shocks)
        K_ts[t]=mean(k_population) # current aggrgate capital
        
        # s_i_base takes 1 when good and 2 when bad 
        s_i_base=ifelse(z==ksp.z_grid[1],1,2)        
        
        # loop over individuals
        for (i,k) in enumerate(k_population)
            eps = eps_shocks[t,i]   # idiosyncratic shock
            s_i=s_i_base+2*(1-Int64(eps))  # transform (z,eps) to s_i
            # obtain next capital holding by interpolation
            
            itp_pol=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,s_i],Gridded(Linear()))
            k_population[i]=itp_pol[k,K_ts[t]]
        end
    end
    
    return nothing
end
Out[5]:
simulate_aggregate_path!

Obtaining aggregate law of motion coefficient

The functions in this cell are used to obtain the coefficient of approximate aggregate capital law of motion (ALM)

  • Methods
    • regress_ALM: regress aggregate capital law of motion and obtain approximate ALM
    • find_ALM_coef!: main iteration. update ALM coefficient until convergence
In [6]:
"""
Obtain new aggregate law of motion coefficients 
using the aggregate capital flaw

##### Arguments
- `ksp` : KSParameter instance
- `z_shock` : aggregate shocks
- `K_ts` : aggregate capital flaw
- `n_discard` : number of discarded samples
"""
function regress_ALM!(ksp::KSParameter,kss::KSSolution,
                        z_shock::Vector{Float64},K_ts::Vector{Float64};
                        n_discard::Int64=100)
    z_grid=ksp.z_grid
    n_g=count(i->(i==z_grid[1]),z_shocks[n_discard+1:end-1])
    n_b=count(i->(i==z_grid[2]),z_shocks[n_discard+1:end-1])
    B_n=Vector{Float64}(4)
    x_g=Vector{Float64}(n_g)
    y_g=Vector{Float64}(n_g)
    x_b=Vector{Float64}(n_b)
    y_b=Vector{Float64}(n_b)
    i_g=0
    i_b=0
    for t = n_discard+1:length(z_shocks)-1
        if z_shocks[t]==z_grid[1]
            i_g=i_g+1
            x_g[i_g]=log(K_ts[t])
            y_g[i_g]=log(K_ts[t+1])
        else
            i_b=i_b+1
            x_b[i_b]=log(K_ts[t])
            y_b[i_b]=log(K_ts[t+1])
        end
    end
    resg=lm(hcat(ones(n_g,1),x_g),y_g)
    resb=lm(hcat(ones(n_b,1),x_b),y_b)
    kss.R2=[r2(resg),r2(resb)]
    B_n[1],B_n[2]=coef(resg)
    B_n[3],B_n[4]=coef(resb)
    dif_B=maximum(abs(B_n-kss.B))
    println("difference of ALM coefficient is $dif_B and B = $B_n")
    return B_n, dif_B
end

function find_ALM_coef!(
                       ksp::KSParameter,
                       kss::KSSolution,
                       z_shocks::Vector{Float64},
                       eps_shocks::Array{Float64,2};
                       tol_ump::AbstractFloat=1e-8,
                       max_iter_ump::Integer=100,
                       Howard_on::Bool=true,
                       Howard_n_iter::Integer=20,
                       tol_B::AbstractFloat=1e-8,
                       max_iter_B::Integer=20,
                       update_B::AbstractFloat=0.3,
                       T_discard::Integer=100,
                       print_skip_VFI::Integer=10,
                       method::Symbol=:Euler,
                       update_k::AbstractFloat=0.3)

    K_ts=similar(z_shocks)
    # populate initial capital holdings
    k_population=37.9893*ones(size(eps_shocks,2))
    counter_B=0
    while true
        counter_B=counter_B+1
        println(" --- Iteration over ALM coefficient: $counter_B ---")

        # solve individual problem
        if method == :VFI
            solve_bellman!(ksp,kss,
                           tol=tol_ump,
                           max_iter=max_iter_ump,
                           Howard=Howard_on,
                           Howard_n_iter=Howard_n_iter,
                           print_skip=print_skip_VFI)
        elseif method == :Euler
            euler_method!(ksp,kss,
                        n_iter=max_iter_ump,
                        tol=tol_ump,
                        update_k=update_k)
        end
        

        # compute aggregate path of capital
        simulate_aggregate_path!(ksp,kss,z_shocks,eps_shocks,k_population,K_ts)
        
        # obtain new ALM coefficient by regression
        B_n,dif_B = regress_ALM!(ksp,kss,z_shocks,K_ts,n_discard=T_discard)

        # check convergence
        if dif_B < tol_B
            println("-----------------------------------------------------")
            println("ALM coefficient successfully converged : dif = $dif_B")
            println("-----------------------------------------------------")
            kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B
            break
        elseif counter_B == max_iter_B
            println("----------------------------------------------------------------")
            println("Iteration over ALM coefficient reached its maximum ($max_iter_B)")
            println("----------------------------------------------------------------")
            kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B
            break
        end
        
        # Update B
        kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B
    end
    return K_ts
end
Out[6]:
find_ALM_coef! (generic function with 1 method)

Plotting

The function is used to plot the path of true path of aggregate capital and approximated one

  • Methods
    • plot_ALM: plot the path of true path of aggregate capital and approximated one
In [7]:
"""
Plot true and approximated ALM of capital 

##### Arguments
- `ksp.z_grid` : aggregate shock grid
- `z_shocks` : aggregate shock
- `kss.B` : ALM coefficient
- `K_ts` : actual path of capital
"""
function plot_ALM(z_grid::Vector{Float64},
                  z_shocks::Vector{Float64},
                  B::Vector{Float64},
                  K_ts::Vector{Float64};
                  T_discard=1000)

    K_ts_approx = similar(K_ts) # preallocation

    # compute approximate ALM for capital
    K_ts_approx[T_discard]=K_ts[T_discard]

    for t=T_discard:length(z_shocks)-1
        if z_shocks[t] == z_grid[1]
            K_ts_approx[t+1]=exp(B[1]+B[2]*log(K_ts_approx[t]))
        elseif z_shocks[t] == z_grid[2]
            K_ts_approx[t+1]=exp(B[3]+B[4]*log(K_ts_approx[t]))
        end
    end

    plot(T_discard+1:length(K_ts),K_ts[T_discard+1:end],lab="true",color=:red,line=:solid)
    plot!(T_discard+1:length(K_ts),K_ts_approx[T_discard+1:end],lab="approximation",color=:blue,line=:dash)
    title!("aggregate law of motion for capital")
    
end
Out[7]:
plot_ALM

Implementation by Euler method

First, construct model instance ksp and initial guess of the solution inside kss

(Grid size inconsistency is also checked, which may return error when exiting result is loaded by load_value=true)

In [8]:
# instance of KSParameter
ksp=KSParameter(k_size=100,
                k_min=1e-16)
# instance of KSSolution
kss=KSSolution(ksp,load_value=false,load_B=false)
if size(kss.k_opt,1) != length(ksp.k_grid)
    error("loaded data is inconsistent with k_size")
end
if size(kss.k_opt,2) != length(ksp.K_grid)
    error("loaded data is inconsistent with K_size")
end

Let's draw the shock for stochastic simulation of aggregate law of motion

In [9]:
# generate shocks
srand(123)
z_shocks,eps_shocks =generate_shocks(ksp;
        z_shock_size=11000,population=5000);

Now, the following cell solves the model with Euler equation method

In [10]:
# find ALM coefficient
@time K_ts = find_ALM_coef!(ksp,kss,z_shocks,eps_shocks,
               tol_ump=1e-8,max_iter_ump=10000,
               Howard_on=true,Howard_n_iter=20,
               tol_B=1e-8, max_iter_B=50,update_B=0.3,
               T_discard=1000,print_skip_VFI=10,
               method=:Euler, update_k=0.7);
 --- Iteration over ALM coefficient: 1 ---
solving individual problem by Euler equation method
Euler method converged with 1868 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.4556643772779668 and B = [0.455664,0.877296,0.438162,0.880418]
 --- Iteration over ALM coefficient: 2 ---
solving individual problem by Euler equation method
Euler method converged with 2005 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.01194706494540887 and B = [0.137682,0.963462,0.119502,0.96685]
 --- Iteration over ALM coefficient: 3 ---
solving individual problem by Euler equation method
Euler method converged with 1479 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.012429438071340365 and B = [0.149423,0.960332,0.135204,0.962694]
 --- Iteration over ALM coefficient: 4 ---
solving individual problem by Euler equation method
Euler method converged with 1516 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.007084614314802812 and B = [0.147807,0.960707,0.134031,0.962974]
 --- Iteration over ALM coefficient: 5 ---
solving individual problem by Euler equation method
Euler method converged with 1159 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.004062551435912554 and B = [0.146911,0.960937,0.133124,0.963222]
 --- Iteration over ALM coefficient: 6 ---
solving individual problem by Euler equation method
Euler method converged with 1044 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.002383876077763747 and B = [0.146451,0.961055,0.132628,0.963361]
 --- Iteration over ALM coefficient: 7 ---
solving individual problem by Euler equation method
Euler method converged with 1065 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0014303729909590501 and B = [0.146213,0.961116,0.132363,0.963436]
 --- Iteration over ALM coefficient: 8 ---
solving individual problem by Euler equation method
Euler method converged with 1050 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0008768443163971185 and B = [0.146088,0.961148,0.132222,0.963476]
 --- Iteration over ALM coefficient: 9 ---
solving individual problem by Euler equation method
Euler method converged with 1017 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0005485109058530613 and B = [0.146023,0.961164,0.132146,0.963498]
 --- Iteration over ALM coefficient: 10 ---
solving individual problem by Euler equation method
Euler method converged with 974 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.00034954334404080356 and B = [0.145988,0.961172,0.132105,0.96351]
 --- Iteration over ALM coefficient: 11 ---
solving individual problem by Euler equation method
Euler method converged with 926 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0002264320715512913 and B = [0.14597,0.961177,0.132083,0.963517]
 --- Iteration over ALM coefficient: 12 ---
solving individual problem by Euler equation method
Euler method converged with 874 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.00014875045461057446 and B = [0.14596,0.961179,0.132071,0.96352]
 --- Iteration over ALM coefficient: 13 ---
solving individual problem by Euler equation method
Euler method converged with 820 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 9.885788332636425e-5 and B = [0.145955,0.96118,0.132065,0.963522]
 --- Iteration over ALM coefficient: 14 ---
solving individual problem by Euler equation method
Euler method converged with 765 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 6.631382514887818e-5 and B = [0.145952,0.96118,0.132062,0.963524]
 --- Iteration over ALM coefficient: 15 ---
solving individual problem by Euler equation method
Euler method converged with 708 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.48072540583655e-5 and B = [0.145951,0.961181,0.13206,0.963524]
 --- Iteration over ALM coefficient: 16 ---
solving individual problem by Euler equation method
Euler method converged with 651 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.0443087124187862e-5 and B = [0.14595,0.961181,0.132059,0.963525]
 --- Iteration over ALM coefficient: 17 ---
solving individual problem by Euler equation method
Euler method converged with 594 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.076821253785277e-5 and B = [0.145949,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 18 ---
solving individual problem by Euler equation method
Euler method converged with 537 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.4209519261981773e-5 and B = [0.145949,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 19 ---
solving individual problem by Euler equation method
Euler method converged with 480 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 9.741695670306694e-6 and B = [0.145949,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 20 ---
solving individual problem by Euler equation method
Euler method converged with 424 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 6.687466508698003e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 21 ---
solving individual problem by Euler equation method
Euler method converged with 370 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.594411339209348e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 22 ---
solving individual problem by Euler equation method
Euler method converged with 317 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.1576579837266916e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 23 ---
solving individual problem by Euler equation method
Euler method converged with 266 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.1703967988118134e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 24 ---
solving individual problem by Euler equation method
Euler method converged with 219 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.4916336017467557e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 25 ---
solving individual problem by Euler equation method
Euler method converged with 176 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.024884198480569e-6 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 26 ---
solving individual problem by Euler equation method
Euler method converged with 137 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.039676926667848e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 27 ---
solving individual problem by Euler equation method
Euler method converged with 103 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.834097304673435e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 28 ---
solving individual problem by Euler equation method
Euler method converged with 74 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.319196261453161e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 29 ---
solving individual problem by Euler equation method
Euler method converged with 50 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.2794120549396446e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 30 ---
solving individual problem by Euler equation method
Euler method converged with 37 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.5634870192959838e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 31 ---
solving individual problem by Euler equation method
Euler method converged with 27 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.0716623247142287e-7 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 32 ---
solving individual problem by Euler equation method
Euler method converged with 20 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.337721727451729e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 33 ---
solving individual problem by Euler equation method
Euler method converged with 18 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 5.00449526541491e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 34 ---
solving individual problem by Euler equation method
Euler method converged with 15 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.408943094473926e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 35 ---
solving individual problem by Euler equation method
Euler method converged with 13 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.3183727299036505e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 36 ---
solving individual problem by Euler equation method
Euler method converged with 10 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.5783928553059212e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 37 ---
solving individual problem by Euler equation method
Euler method converged with 8 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.0743600203921844e-8 and B = [0.145948,0.961181,0.132058,0.963525]
 --- Iteration over ALM coefficient: 38 ---
solving individual problem by Euler equation method
Euler method converged with 7 iterations
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.295628601244886e-9 and B = [0.145948,0.961181,0.132058,0.963525]
-----------------------------------------------------
ALM coefficient successfully converged : dif = 7.295628601244886e-9
-----------------------------------------------------
16333.330134 seconds (20.04 G allocations: 9.153 TB, 27.86% gc time)

Let's compare the true aggreate law of motion for capital and approximated one with figure and regression

In [11]:
plot_ALM(ksp.z_grid,z_shocks,
         kss.B,K_ts)
Out[11]:
In [12]:
#kss.B  # Regression coefficient
println("Approximated aggregate capital law of motion")
println("log(K_{t+1})=$(kss.B[1])+$(kss.B[2])log(K_{t}) in good time (R2 = $(kss.R2[1]))")
println("log(K_{t+1})=$(kss.B[3])+$(kss.B[4])log(K_{t}) in bad time (R2 = $(kss.R2[2]))")
Approximated aggregate capital law of motion
log(K_{t+1})=0.14594821741362846+0.9611811624862514log(K_{t}) in good time (R2 = 0.9999988558436316)
log(K_{t+1})=0.13205800455894173+0.9635249205659238log(K_{t}) in bad time (R2 = 0.9999972131809477)

The approximated law of motion of capital is very close to the true one, which implies that assuming agents are partially rational is not bad idea since the difference of their actions are negligible.

Note: The mean of capital, about 40, is sufficiently close to Maliar, Maliar, Valli but higher than Krusell-Smith.

In [13]:
save("result_Euler.jld","kss",kss)
In [14]:
# Compute mean of capital implied by regression
mc=MarkovChain(ksp.TransMat.Pz)
sd=stationary_distributions(mc)[1]
logKg=kss.B[1]/(1-kss.B[2])
logKb=kss.B[3]/(1-kss.B[4])
meanK_reg=exp(sd[1]*logKg+sd[2]*logKb)
meanK_sim=mean(K_ts[1001:end])
println("mean of capital implied by regression is $meanK_reg")
println("mean of capital implied by simulation is $meanK_sim")
mean of capital implied by regression is 40.049373762168585
mean of capital implied by simulation is 40.160799927814296

Figures in Krusell-Smith

Let's prepare some functions for figures in Krusell-Smith

In [15]:
function plot_Fig1(ksp,kss,K_ts)
    B=kss.B
    K_min, K_max = minimum(K_ts), maximum(K_ts)
    K_lim=collect(linspace(K_min,K_max,100))
    Kp_g=exp(B[1]+B[2]*log(K_lim))
    Kp_b=exp(B[3]+B[4]*log(K_lim))
    
    p=plot(K_lim,Kp_g,linestyle=:solid,lab="Good")
    plot!(p,K_lim,Kp_b,linestyle=:solid,lab="Bad")
    plot!(p,K_lim,K_lim,color=:black,linestyle=:dash,lab="45 degree",width=0.5)
    title!(p,"FIG1: Tomorrow's vs. today's aggregate capital")
    p
end

function plot_Fig2(ksp,kss,K_eval_point)
    k_lim=collect(linspace(0,80,1000))
    itp_e=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,1],Gridded(Linear()))
    itp_u=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,3],Gridded(Linear()))
    
    kp_e(k)=itp_e[k,K_eval_point]
    kp_u(k)=itp_u[k,K_eval_point]
    
    p=plot(k_lim,kp_e.(k_lim),linestyle=:solid,lab="employed")
    plot!(p,k_lim,kp_u.(k_lim),linestyle=:solid,lab="unemployed")
    plot!(p,k_lim,k_lim,color=:black,linestyle=:dash,lab="45 degree",width=0.5)
    title!(p,"FIG2: Individual policy function \n at K=$K_eval_point when good state")
    p
end
Out[15]:
plot_Fig2 (generic function with 1 method)

Now, plot the replication figure:

Figure 1

In [16]:
plot_Fig1(ksp,kss,K_ts)
Out[16]:

Figure 2

In [17]:
plot_Fig2(ksp,kss,40)
Out[17]:

Both figures are replicated well.

Note that the mean of capital is approximately 40 in this replication, which is different from Krusell-Smith but same as Maliar, Maliar, Valli. Therefore, Figure 1 is plotted around 40 and policy function for Figure 2 is evaluated at K=40

Solution with value function iteration

In this section, each agent's utility maximization problem is solved by value function iteration.

Let's prepare functions for value function iteration

Bellman equation and VFI

The functions in this cell are used to solve individual household problem by VFI

  • Methods
    • rhs_bellman: evaluate RHS of bellman equation given the guess of value function
    • compute_expectation_VFI: compute expectation term of Bellman equation given current shock state
    • solve_bellman_once!: maximize RHS of Bellman equation for an agent of state (k_i,K_i,s_i)
    • solve_bellman!: maximize RHS of Bellman equation for all state until convergence
    • iterate_policy!: iterating policy and compute values under the policy. used for Howard
In [18]:
"""
Compute right hand side of bellman equation

##### Arguments
- `kp` : next period capital
- `ksm` : KSModel instance
- `k` : current individual capital
- `K` : current aggregate capital
- `L` : current labor
- `zeps_i` : 
"""
function rhs_bellman(ksp::KSParameter,
            kp::AbstractFloat,value::Array{Float64,3},
            k::AbstractFloat,K::AbstractFloat,s_i::Integer)
    u,s_grid,beta,alpha,l_bar,delta, mu =
        ksp.u, ksp.s_grid, ksp.beta, ksp.alpha, ksp.l_bar, ksp.delta, ksp.mu
    z, eps = get_shock(s_grid,s_i)
    Kp,L = compute_Kp_L(K,s_i,kss.B,ksp) # Next period aggregate capital and current aggregate labor
    c = (r(alpha,z,K,L)+1-delta)*k+
          w(alpha,z,K,L)*(eps*l_bar+(1.0-eps)*mu)-kp # current consumption 
    expec = compute_expectation(kp,Kp,value,s_i,ksp)
    return u(c)+beta*expec
end

"""
Compute expectation of next period value

##### Arguments
- `kp` : next period individual capital
- `Kp` : next period aggregate capital
- `value` : given value 
- `s_i` : current shock state
- `ksp` : KSParameter instance
"""
function compute_expectation(
              kp::AbstractFloat,  # next period indicidual capital
              Kp::AbstractFloat,  # next period aggragte capital
              value::Array{Float64,3}, # next period value
              s_i::Integer, # index of current state,
              ksp::KSParameter 
              )
    k_grid, K_grid = ksp.k_grid, ksp.K_grid # unpack grid
    beta, P = ksp.beta, ksp.TransMat.P      # unpack parameters
    
    # compute expectations by summing up
    expec=0.0
    for s_n_i=1:ksp.s_size
        value_itp=interpolate((k_grid,K_grid),value[:,:,s_n_i],Gridded(Linear()))
        expec = expec + P[s_i,s_n_i]*value_itp[kp,Kp]
    end
    return expec
end
"""
Solve bellman equation for all states once

##### Arguments
- `k` : individual capital
- `K` : aggregate capital
- `s_i` : shock state index
- `ksp` : KSParameter
- `kss` : KSSolution
"""
function solve_bellman_once!(
            k_i::Integer,
            K_i::Integer,
            s_i::Integer,
            ksp::KSParameter,
            kss::KSSolution,
            )
    # obtain minimum and maximum of grid
    k_min, k_max = ksp.k_grid[1], ksp.k_grid[end]
    
    # unpack parameters
    alpha,delta,l_bar, mu = 
        ksp.alpha, ksp.delta, ksp.l_bar, ksp.mu
    
    # obtain state value
    k=ksp.k_grid[k_i]   # obtain individual capital value
    K=ksp.K_grid[K_i]   # obtain aggregate capital value
    z, eps = get_shock(ksp.s_grid,s_i) # obtain shock value
    
    Kp,L=compute_Kp_L(K,s_i,kss.B,ksp) # next aggregate capital and current aggregate labor
    # if kp>k_c_pos, consumption is negative 
    k_c_pos=(r(alpha,z,K,L)+1-delta)*k+
        w(alpha,z,K,L)*(eps*l_bar+(1-eps)*mu)
    obj(kp)=-rhs_bellman(ksp,kp,kss.value,k,K,s_i) # objective function
    res=optimize(obj,k_min,min(k_c_pos,k_max)) # maximize value
    # obtain result
    kss.k_opt[k_i,K_i,s_i]=Optim.minimizer(res) 
    kss.value[k_i,K_i,s_i]=-Optim.minimum(res)
    return nothing
end

"""
Solve bellman equation for all states until convergence

##### Arguments
- `ksp` : KSParameter
- `kss` : KSSolution
- `tol` : tolerance of value function difference
- `max_iter` : maximum number of iteration
"""
function solve_bellman!(
            ksp::KSParameter,
            kss::KSSolution;
            tol::AbstractFloat=1e-8,
            max_iter::Integer=100,
            Howard::Bool=false,
            Howard_n_iter::Integer=20,
            print_skip::Integer=10)
    counter_VFI=0  # counter
    while true
        counter_VFI += 1
        value_old=copy(kss.value) # guessed value
        # maximize value for all state
        [solve_bellman_once!(k_i,K_i,s_i,ksp,kss)
          for k_i in 1:ksp.k_size, K_i in 1:ksp.K_size, s_i in 1:ksp.s_size]
        # Howard's policy iteration
        !Howard || iterate_policy!(ksp,kss,n_iter=Howard_n_iter)
        # difference of guessed and new value
        dif=maxabs(value_old-kss.value)
        # print covergence process
        !(counter_VFI % print_skip ==0) ||
            println("VFI iteration $counter_VFI : dif = $dif")
        # if difference is sufficiently small
        if dif<tol
            println(" ** VFI converged successfully!! dif = $dif")
            break
        elseif counter_VFI >= max_iter
            println("VFI reached its maximum iteration : $max_iter")
            break
        end
    end
end

"""
Iterate the value function fixing the policy function

##### Arguments
- `ksp` : KSParameter instance
- `kss` : KSSolution instance
- `n_iter` : number of iterations
"""
function iterate_policy!(ksp::KSParameter,
            kss::KSSolution;n_iter::Int64=20)
    value=similar(kss.value)
    for i=1:n_iter
        # update value using policy
        value .= 
            [rhs_bellman(ksp,
                kss.k_opt[k_i,K_i,s_i],kss.value,
                ksp.k_grid[k_i],ksp.K_grid[K_i],s_i)
                for k_i in 1:ksp.k_size,
                    K_i in 1:ksp.K_size,
                    s_i in 1:ksp.s_size]
        kss.value.=copy(value)
    end
    return nothing
end
Out[18]:
iterate_policy!

Implementation

Let's skip the following steps in this section to save computational time.

  • consturuction of ksp instance since it is same
  • consturuction of kss instance to use the previous result as initial guess of the solution
  • draws of the shocks to use same ones

However, instead of constructing kss again, obtain value from the policy function derived by Euler method:

In [19]:
iterate_policy!(ksp,kss,n_iter=30);

Now, the model is solved by VFI in the next cell

In [20]:
# find ALM coefficient
@time K_ts = find_ALM_coef!(ksp,kss,z_shocks,eps_shocks,
               tol_ump=1e-8,max_iter_ump=10000,
               Howard_on=true,Howard_n_iter=20,
               tol_B=1e-8, max_iter_B=50,update_B=0.3,
               T_discard=1000,print_skip_VFI=10,
               method=:VFI, update_k=0.7);
 --- Iteration over ALM coefficient: 1 ---
VFI iteration 10 : dif = 4.364130337367499
VFI iteration 20 : dif = 0.5155086466810985
VFI iteration 30 : dif = 0.06101406344856741
VFI iteration 40 : dif = 0.007215234574232454
VFI iteration 50 : dif = 0.000852061281875649
VFI iteration 60 : dif = 0.00010046707080846318
VFI iteration 70 : dif = 1.1826375214241125e-5
VFI iteration 80 : dif = 1.3896425912207633e-6
VFI iteration 90 : dif = 1.6298264426950482e-7
VFI iteration 100 : dif = 1.9079607227467932e-8
 ** VFI converged successfully!! dif = 8.085862646112218e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.02551266404864358 and B = [0.130832,0.964894,0.106545,0.970074]
 --- Iteration over ALM coefficient: 2 ---
VFI iteration 10 : dif = 0.04110720950990299
VFI iteration 20 : dif = 0.004799456575142358
VFI iteration 30 : dif = 0.0005606688813486471
VFI iteration 40 : dif = 6.551400201715296e-5
VFI iteration 50 : dif = 7.655140279894113e-6
VFI iteration 60 : dif = 8.94252366379078e-7
VFI iteration 70 : dif = 1.0442380471431534e-7
VFI iteration 80 : dif = 1.2187229003757238e-8
 ** VFI converged successfully!! dif = 9.830898761720164e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.010571478106761914 and B = [0.138905,0.962971,0.113833,0.968363]
 --- Iteration over ALM coefficient: 3 ---
VFI iteration 10 : dif = 0.0004685646893847206
VFI iteration 20 : dif = 5.4314097894803126e-5
VFI iteration 30 : dif = 6.4477072214685904e-6
VFI iteration 40 : dif = 7.622459747835819e-7
VFI iteration 50 : dif = 8.993168876259006e-8
VFI iteration 60 : dif = 1.0597830168990185e-8
 ** VFI converged successfully!! dif = 8.557094588468317e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.006449099771311431 and B = [0.140108,0.962642,0.114784,0.968106]
 --- Iteration over ALM coefficient: 4 ---
VFI iteration 10 : dif = 0.0004316023300816596
VFI iteration 20 : dif = 5.057148842979586e-5
VFI iteration 30 : dif = 6.002880127198296e-6
VFI iteration 40 : dif = 7.094776037774864e-7
VFI iteration 50 : dif = 8.366657766600838e-8
VFI iteration 60 : dif = 1.0639979564075475e-8
 ** VFI converged successfully!! dif = 7.953985914355144e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0037177793034470913 and B = [0.140774,0.962464,0.11558,0.967893]
 --- Iteration over ALM coefficient: 5 ---
VFI iteration 10 : dif = 0.0001699352112325414
VFI iteration 20 : dif = 1.9677074277524298e-5
VFI iteration 30 : dif = 2.2995970141437283e-6
VFI iteration 40 : dif = 2.691489555672888e-7
VFI iteration 50 : dif = 3.161997597089794e-8
 ** VFI converged successfully!! dif = 8.755819180805702e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0021012252914709983 and B = [0.141146,0.962364,0.116081,0.967759]
 --- Iteration over ALM coefficient: 6 ---
VFI iteration 10 : dif = 0.00016965897827958543
VFI iteration 20 : dif = 1.848068291110394e-5
VFI iteration 30 : dif = 2.05708909106761e-6
VFI iteration 40 : dif = 2.325381274204119e-7
VFI iteration 50 : dif = 2.647169594638399e-8
 ** VFI converged successfully!! dif = 8.945391982706496e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0012679245607524897 and B = [0.141284,0.962326,0.116284,0.967703]
 --- Iteration over ALM coefficient: 7 ---
VFI iteration 10 : dif = 4.9276457787073014e-5
VFI iteration 20 : dif = 4.933374896154419e-6
VFI iteration 30 : dif = 5.170394956621749e-7
VFI iteration 40 : dif = 5.5685518418613356e-8
 ** VFI converged successfully!! dif = 9.463576589041622e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0007726858753848709 and B = [0.141347,0.962309,0.116399,0.967672]
 --- Iteration over ALM coefficient: 8 ---
VFI iteration 10 : dif = 3.44231706321807e-5
VFI iteration 20 : dif = 3.568709871615283e-6
VFI iteration 30 : dif = 3.8495443277497543e-7
VFI iteration 40 : dif = 4.244049023327534e-8
 ** VFI converged successfully!! dif = 9.12990572032868e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.00048209680162884794 and B = [0.141369,0.962302,0.116458,0.967656]
 --- Iteration over ALM coefficient: 9 ---
VFI iteration 10 : dif = 1.758010006369659e-5
VFI iteration 20 : dif = 1.7869241446533124e-6
VFI iteration 30 : dif = 1.9031887177334283e-7
VFI iteration 40 : dif = 2.07742800739652e-8
 ** VFI converged successfully!! dif = 8.598362910561264e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.0003055378304484685 and B = [0.141377,0.9623,0.11649,0.967647]
 --- Iteration over ALM coefficient: 10 ---
VFI iteration 10 : dif = 1.023216435669383e-5
VFI iteration 20 : dif = 1.03666968698235e-6
VFI iteration 30 : dif = 1.102799274121935e-7
VFI iteration 40 : dif = 1.2026134754705708e-8
 ** VFI converged successfully!! dif = 9.64348600973608e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.00019619175347125595 and B = [0.141377,0.9623,0.116508,0.967642]
 --- Iteration over ALM coefficient: 11 ---
VFI iteration 10 : dif = 6.056440042812028e-6
VFI iteration 20 : dif = 6.10131223766075e-7
VFI iteration 30 : dif = 6.479382363977493e-8
 ** VFI converged successfully!! dif = 8.801350759313209e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 0.00012730360336368762 and B = [0.141376,0.962301,0.116518,0.96764]
 --- Iteration over ALM coefficient: 12 ---
VFI iteration 10 : dif = 3.566279190181376e-6
VFI iteration 20 : dif = 3.556169758667238e-7
VFI iteration 30 : dif = 3.754672661671066e-8
 ** VFI converged successfully!! dif = 9.877851425699191e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 8.313458897782344e-5 and B = [0.141374,0.962301,0.116524,0.967638]
 --- Iteration over ALM coefficient: 13 ---
VFI iteration 10 : dif = 2.2697702775076323e-6
VFI iteration 20 : dif = 2.2698196744386223e-7
VFI iteration 30 : dif = 2.4034420675889123e-8
 ** VFI converged successfully!! dif = 9.869268069451209e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 5.463237199027082e-5 and B = [0.141373,0.962302,0.116527,0.967637]
 --- Iteration over ALM coefficient: 14 ---
VFI iteration 10 : dif = 1.2852221971115796e-6
VFI iteration 20 : dif = 1.2533530480141053e-7
VFI iteration 30 : dif = 1.3040221347182523e-8
 ** VFI converged successfully!! dif = 8.32721980259521e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.63293125280606e-5 and B = [0.141371,0.962302,0.116529,0.967636]
 --- Iteration over ALM coefficient: 15 ---
VFI iteration 10 : dif = 9.561501315147325e-7
VFI iteration 20 : dif = 9.57708721216477e-8
VFI iteration 30 : dif = 1.017394879454514e-8
 ** VFI converged successfully!! dif = 8.143899776769103e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.4483813326064974e-5 and B = [0.14137,0.962302,0.116531,0.967636]
 --- Iteration over ALM coefficient: 16 ---
VFI iteration 10 : dif = 5.508533149622963e-7
VFI iteration 20 : dif = 5.388682211560081e-8
 ** VFI converged successfully!! dif = 8.823121788736898e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.636877019500771e-5 and B = [0.14137,0.962302,0.116532,0.967636]
 --- Iteration over ALM coefficient: 17 ---
VFI iteration 10 : dif = 4.021055701741716e-7
VFI iteration 20 : dif = 4.010468046544702e-8
 ** VFI converged successfully!! dif = 8.30476665214519e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.0984150278420257e-5 and B = [0.141369,0.962302,0.116532,0.967636]
 --- Iteration over ALM coefficient: 18 ---
VFI iteration 10 : dif = 2.3281569383470924e-7
VFI iteration 20 : dif = 2.2501069452118827e-8
 ** VFI converged successfully!! dif = 9.03918362382683e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.2048826815218625e-6 and B = [0.141369,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 19 ---
VFI iteration 10 : dif = 1.230218344971945e-7
VFI iteration 20 : dif = 1.0416044915473321e-8
 ** VFI converged successfully!! dif = 8.250196970038814e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.988788248755371e-6 and B = [0.141369,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 20 ---
VFI iteration 10 : dif = 1.0984092568833148e-7
VFI iteration 20 : dif = 1.089637180484715e-8
 ** VFI converged successfully!! dif = 8.687948138685897e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.290747322109988e-6 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 21 ---
VFI iteration 10 : dif = 8.936746098697768e-8
VFI iteration 20 : dif = 9.178563686873531e-9
 ** VFI converged successfully!! dif = 9.178563686873531e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.158295437842961e-6 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 22 ---
VFI iteration 10 : dif = 3.8132498048071284e-8
 ** VFI converged successfully!! dif = 9.028951808431884e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.5509704746619057e-6 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 23 ---
VFI iteration 10 : dif = 1.3471407100951183e-8
 ** VFI converged successfully!! dif = 8.787935712462058e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.0163098228543888e-6 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 24 ---
VFI iteration 10 : dif = 4.011184273622348e-8
 ** VFI converged successfully!! dif = 8.319318567373557e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 6.248587137436257e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 25 ---
VFI iteration 10 : dif = 2.6097382033185568e-8
 ** VFI converged successfully!! dif = 8.564882136852248e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.880477443818364e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 26 ---
VFI iteration 10 : dif = 2.0421111912583e-8
 ** VFI converged successfully!! dif = 8.35751734484802e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.843728103413268e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 27 ---
VFI iteration 10 : dif = 3.812928639490565e-8
 ** VFI converged successfully!! dif = 8.57767190609593e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 2.919172302912054e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 28 ---
VFI iteration 10 : dif = 1.3114174635120435e-8
 ** VFI converged successfully!! dif = 8.181245902960654e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.6344598458006843e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 29 ---
VFI iteration 10 : dif = 5.4228507906373125e-8
 ** VFI converged successfully!! dif = 9.760071861819597e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 9.524729271959131e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 30 ---
VFI iteration 10 : dif = 3.021312977580237e-8
 ** VFI converged successfully!! dif = 8.285155672638211e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.177423369530977e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 31 ---
VFI iteration 10 : dif = 2.392872033851745e-8
 ** VFI converged successfully!! dif = 8.194376732717501e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.0872591149624355e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 32 ---
VFI iteration 10 : dif = 5.144397619005758e-8
 ** VFI converged successfully!! dif = 9.24671894608764e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.570228442410171e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 33 ---
VFI iteration 10 : dif = 2.3067457277647918e-8
 ** VFI converged successfully!! dif = 9.729632211019634e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.874108731709129e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 34 ---
VFI iteration 10 : dif = 7.674799462620285e-8
VFI iteration 20 : dif = 8.989559319161344e-9
 ** VFI converged successfully!! dif = 8.989559319161344e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.674765916561462e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 35 ---
VFI iteration 10 : dif = 3.5518326058081584e-8
 ** VFI converged successfully!! dif = 9.848633908404736e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.1407666089535695e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 36 ---
VFI iteration 10 : dif = 5.78716878862906e-8
 ** VFI converged successfully!! dif = 8.39969516164274e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 3.834874243158204e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 37 ---
VFI iteration 10 : dif = 8.088653657978284e-8
VFI iteration 20 : dif = 9.574307568982476e-9
 ** VFI converged successfully!! dif = 9.574307568982476e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.068997078165992e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 38 ---
VFI iteration 10 : dif = 6.087861947889905e-8
 ** VFI converged successfully!! dif = 8.969834652816644e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.422016042227028e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 39 ---
 ** VFI converged successfully!! dif = 9.783235555005376e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.2206419956750647e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 40 ---
VFI iteration 10 : dif = 7.822649195077247e-8
VFI iteration 20 : dif = 9.270820555684622e-9
 ** VFI converged successfully!! dif = 9.270820555684622e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.8181412958506726e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 41 ---
VFI iteration 10 : dif = 1.760491841196199e-8
 ** VFI converged successfully!! dif = 9.28258714338881e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 5.6862735645091256e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 42 ---
VFI iteration 10 : dif = 5.217941634327872e-8
 ** VFI converged successfully!! dif = 9.463349215366179e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 9.779224406647469e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 43 ---
VFI iteration 10 : dif = 2.28354224418581e-8
 ** VFI converged successfully!! dif = 9.649710364101338e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 4.013947332848211e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 44 ---
 ** VFI converged successfully!! dif = 7.641517640877282e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 8.247222312018909e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 45 ---
VFI iteration 10 : dif = 5.345765430320171e-8
 ** VFI converged successfully!! dif = 9.702773695607902e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 5.5271162438530475e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 46 ---
VFI iteration 10 : dif = 6.705204214085825e-8
 ** VFI converged successfully!! dif = 9.735970252222614e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 5.700652201678924e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 47 ---
VFI iteration 10 : dif = 2.999195203301497e-8
 ** VFI converged successfully!! dif = 8.3674080997298e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.5016738530437834e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 48 ---
VFI iteration 10 : dif = 1.08337872006814e-8
 ** VFI converged successfully!! dif = 8.718018307263264e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 1.779748049768326e-7 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 49 ---
VFI iteration 10 : dif = 1.6313094874931267e-8
 ** VFI converged successfully!! dif = 8.512699878338026e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 6.603349556044691e-8 and B = [0.141368,0.962303,0.116533,0.967635]
 --- Iteration over ALM coefficient: 50 ---
VFI iteration 10 : dif = 7.444867833328317e-8
VFI iteration 20 : dif = 8.854442512529204e-9
 ** VFI converged successfully!! dif = 8.854442512529204e-9
simulating aggregate path ... please wait ... 
difference of ALM coefficient is 7.825109485382065e-8 and B = [0.141368,0.962303,0.116533,0.967635]
----------------------------------------------------------------
Iteration over ALM coefficient reached its maximum (50)
----------------------------------------------------------------
26895.406260 seconds (27.08 G allocations: 12.690 TB, 25.15% gc time)

The following same exercises show that the main result is same as before

In [21]:
plot_ALM(ksp.z_grid,z_shocks,
         kss.B,K_ts)
Out[21]:
In [22]:
#kss.B  # Regression coefficient
println("Approximated aggregate capital law of motion")
println("log(K_{t+1})=$(kss.B[1])+$(kss.B[2])log(K_{t}) in good time (R2 = $(kss.R2[1]))")
println("log(K_{t+1})=$(kss.B[3])+$(kss.B[4])log(K_{t}) in bad time (R2 = $(kss.R2[2]))")
Approximated aggregate capital law of motion
log(K_{t+1})=0.14136817405539825+0.9623027152453734log(K_{t}) in good time (R2 = 0.9999992137407112)
log(K_{t+1})=0.11653324573318116+0.9676352636916012log(K_{t}) in bad time (R2 = 0.9999982189177621)
In [23]:
save("result_VFI.jld","kss",kss)
In [24]:
# Compute mean of capital implied by regression
mc=MarkovChain(ksp.TransMat.Pz)
sd=stationary_distributions(mc)[1]
logKg=kss.B[1]/(1-kss.B[2])
logKb=kss.B[3]/(1-kss.B[4])
meanK_reg=exp(sd[1]*logKg+sd[2]*logKb)
meanK_sim=mean(K_ts[1001:end])
println("mean of capital implied by regression is $meanK_reg")
println("mean of capital implied by simulation is $meanK_sim")
mean of capital implied by regression is 39.46272307237643
mean of capital implied by simulation is 39.70193508699299

Figures in Krusell-Smith

Figure 1

In [25]:
plot_Fig1(ksp,kss,K_ts)
Out[25]:

Figure 2

In [26]:
plot_Fig2(ksp,kss,40)
Out[26]:

Again, the figures are succesfully replicated.