Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Transition Dynamics and MIT shock

The goal of this note is to show how to solve the transition dynamics of a standard incomplete market model. We use the Endogenous Grid Method (EGM) to solve the dynamic programming problem and non-stochastic simulation to solve for the distribution of agents in the economy.

After we solve for the Impulse Response Function of a TFP shock, we simulate the model using the method of Boppart-Krusell-Mitman (2018, JEDC). Finally, we show how to compute the transition path between two different steady states by simulating an exogenous change in the labor tax.

The model is pretty standard. Instead of solving for the Stationary Equilibrium, we must solve for the Sequential Equilibrium. That is, the equilibrium in the asset market (and any other) must hold for all periods $t$.

The production function is Cobb-Douglas. There is now time-varying total factor productivity $Z_t$. TFP follows an AR(1) process.

$$Y_t = Z_t K_t^\alpha L^{1-\alpha},$$

where $\log Z_t = \rho_z \log Z_{t-1} + \sigma_z \varepsilon_t$.

The consumption-savings problem is summarized by the following sequential Dynamic Programming problem:

$$V_t (a, s) = \max_{a' \geq -\phi} \{u((1+r_t )a + w_t \exp{s} -a' ) + \beta\mathbb{E}[V_{t+1} (a', s')|s] \},$$

where $s$ follows an AR(1) process:

$$s_t = \rho s_{t-1} + \sigma \varepsilon_t,$$

where $\varepsilon \sim N(0, 1)$.

The main difference is that the Value Function is indexed by the time index $t$. The time index summarizes all relevant time varying information for the household, such as prices.

Since we have to solve for the equilibrium in all periods, we have to aggregate the wealth of all households using the distribution $\lambda(a,s)$ of households and get the aggregate capital.

The distribution will be time-varying and have to satisfy the time-varying law-of-motion:

$$\lambda_t (\mathcal{A} \times \mathcal{S}) = \int_{A \times S} Q_t((a, s), \mathcal{A} \times \mathcal{S})) d \lambda_t.$$

Intuitively, since the transition function $Q_t$ depends on the household policy functions (that are time-varying), the distribution would be evolving during the transition.

Finally, the interest rate, $r_t$, must clear the asset market in all periods:

$$\int_{A\times S} a d\lambda_t(a,s) = K_t.$$

Numerical implementation

First, we define packages, auxiliary functions, and all functions used in solving the model in the steady state.
md"""
## Transition Dynamics and MIT shock

The goal of this note is to show how to solve the transition dynamics of a standard incomplete market model. We use the Endogenous Grid Method (EGM) to solve the dynamic programming problem and non-stochastic simulation to solve for the distribution of agents in the economy.

After we solve for the Impulse Response Function of a TFP shock, we simulate the model using the method of Boppart-Krusell-Mitman (2018, JEDC). Finally, we show how to compute the transition path between two different steady states by simulating an exogenous change in the labor tax.

The model is pretty standard. Instead of solving for the Stationary Equilibrium, we must solve for the Sequential Equilibrium. That is, the equilibrium in the asset market (and any other) must hold for all periods $t$.

The production function is Cobb-Douglas. There is now time-varying total factor productivity $Z_t$. TFP follows an AR(1) process.

$Y_t = Z_t K_t^\alpha L^{1-\alpha},$

where $\log Z_t = \rho_z \log Z_{t-1} + \sigma_z \varepsilon_t$.

The consumption-savings problem is summarized by the following sequential Dynamic Programming problem:

$V_t (a, s) = \max_{a' \geq -\phi} \{u((1+r_t )a + w_t \exp{s} -a' ) + \beta\mathbb{E}[V_{t+1} (a', s')|s] \},$

where $s$ follows an AR(1) process:

$s_t = \rho s_{t-1} + \sigma \varepsilon_t,$

where $\varepsilon \sim N(0, 1)$.

The main difference is that the Value Function is indexed by the time index $t$. The time index summarizes all relevant time varying information for the household, such as prices.

Since we have to solve for the equilibrium in all periods, we have to aggregate the wealth of all households using the distribution $\lambda(a,s)$ of households and get the aggregate capital.

The distribution will be time-varying and have to satisfy the time-varying law-of-motion:

$\lambda_t (\mathcal{A} \times \mathcal{S}) = \int_{A \times S} Q_t((a, s), \mathcal{A} \times \mathcal{S})) d \lambda_t.$

Intuitively, since the transition function $Q_t$ depends on the household policy functions (that are time-varying), the distribution would be evolving during the transition.

Finally, the interest rate, $r_t$, must clear the asset market in all periods:
1.9 ms
## ======================================================================== ##
# This block import the packages we use
## ======================================================================== ##

begin
using Parameters # enable @unpack
using Random, Distributions # enable cdf of Normal
using Roots # root finding routines: brent
using PlutoUI # printing stuff in the Pluto notebook
using Plots # plots: comment out if you wont use it, since it is heavy.
end
7.9 s
coordGridIntp! (generic function with 1 method)
## ======================================================================== ##
# This Defines some auxiliary functions
## ======================================================================== ##

begin


function rouwenhorst(N, ρ, σ; μ = 0.0)
q = (ρ+1.0)/2
nu = ((N-1.0)/(1.0-ρ^2))^(1/2)*σ
s = collect(range(μ/(1.0-ρ)-nu, μ/(1.0-ρ)+nu, length = N))

P = [q 1-q; 1-q q]
for i in 2:N-1
P = q*[P zeros(i,1); zeros(1,i+1)] + (1-q)*[zeros(i,1) P; zeros(1,i+1)] + (1-q)*[zeros(1,i+1); P zeros(i,1)] + q*[zeros(1,i+1); zeros(i,1) P]
P[2:i,:]=P[2:i,:]/2
end
return s, P
end


"""
linearGridIntp!(xgrid, ygrid, xtarget, ytarget)

Returns a linear interpolation of the vector xtarget on grid (xgrid, ygrid).
Assume xgrid and xtarget is sorted.
# Arguments
- `xgrid::Array{AbstractFloat,1}` : xgrid.
- `ygrid::Array{AbstractFloat,1}` : ygrid.
- `xtarget::Array{AbstractFloat,1}` : grid to interpolate.
# Output (inplace)
- `ytarget::Array{AbstractFloat,1}`
"""
function linearGridIntp!(xgrid::Array{Float64,1}, ygrid::Array{Float64,1}, xtarget::Array{Float64,1}, ytarget::Array{Float64,1})

# check if is sorted
19.6 ms

Here we have the functions to solve for the stationary equilibrium. They are pretty much the same functions used in the previous notebook. The only difference is that we include fiscal policy parameters, such as tax and transfers, and a government budget constraint equilibrium condition (check the excess demand function!).

md"""
Here we have the functions to solve for the stationary equilibrium. They are pretty much the same functions used in the previous notebook. The only difference is that we include fiscal policy parameters, such as tax and transfers, and a government budget constraint equilibrium condition (check the excess demand function!).
"""

220 μs
ModelSolutionSS (generic function with 1 method)
## ======================================================================== ##
# Define all functions used to solve for the steady state.
## ======================================================================== ##

begin
function solveHHproblemEGM(param, r, w, τl, τk, Transf)
@unpack nA, nS, gA, gS, transS, β, γ, δ = param

# Define parameters
gAA = repeat(gA, 1, nS) # nA x nS matrix of asset grid.
gSS = repeat(gS', nA, 1) # nA x nS matrix of labor grid. Note the transposed
gYY = (1-τl)*w.*gSS .+ (1.0+r*(1-τk)).*gAA .+ Transf

maxiter = 2000
tol = 1E-13

# Pre-allocate value and policy funcctions
cp = gYY - repeat(gA, 1, nS) # initial guess is policy on-grid
#cpNew = copy(cp)
ap = copy(cp)
endogAA = copy(cp)


# function that runs one iteration
function endoGridIter(cp)
expect = β*(1.0+r*((1-τk)))* cp.^-γ*transS' # RHS of EE
c = expect.^(-1.0/γ) # use EE to get contemporaneous C
endogAA = (c + gAA - w.*gSS*(1-τl) .- Transf)/(1+r*(1-τk)) # note that the gAA is the policy

apNew = zeros(nA, nS)
for is = 1:nS
apNew[:, is] = linearGridIntp!(endogAA[:, is], gA, gA, apNew[:, is])

# Binding constraint
for ia = 1:nA
19.2 ms
Parameters

In the function below, we define the parameters. Discretization is standard. We simulate the model for 150 periods. Note that we are already simulating the response of $Z_t$ to the shock. We also define some policy parameters in the initial and final steady-state.

md"""
##### Parameters

In the function below, we define the parameters. Discretization is standard. We simulate the model for 150 periods. Note that we are already simulating the response of $Z_t$ to the shock. We also define some policy parameters in the initial and final steady-state.

"""
302 μs
setPar (generic function with 1 method)
function setPar(;
nA = 300, # Asset grid size
nS = 7, # Labor endowment grid size
α = 0.33,
δ = 0.05,
β = 0.96,
γ = 2.0, # CRRA parameter
ϕ = 0.0, # borrowing constraint
ρ = 0.9, # autocorrelation
σ = 0.03, # std deviation
Zbar = 0.0, # SS aggregate shock (in logs)
ρz = 0.95, #
nT = 150, # periods

# TAXES/GOV. PARAMETERS
sizeshock = 0.01, # size of shock (1%)
τl0 = 0.0, # labor tax ss 0
τl1 = 0.0, # labor tax ss 1
τk0 = 0.0, # capital tax ss 0
τk1 = 0.0, # capital tax ss 1
G0 = 0.0, # GOV EXPENDITURE 0
G1 = 0.0, # GOV EXPENDITURE 1

amax = 250.0, # maximum grid point
grid_growth = 0.025 # asset grid curvature (growth rate between points)

)


# ===== Define the labor grid =====
gS, transS = rouwenhorst(nS, ρ, σ)
gS = exp.(gS)

# compute invariant distribution of labor process (by iteration)
invS = ones(nS) / nS # initial guess
tol=1E-11; maxit=10^4
10.2 msparam
param = setPar()
79.5 μs
Dynamic Programming Problem

We now write a function that solves for the policy functions for all $t$. The function takes as inputs the parameters, an interest rate and wage sequence, a sequence of possible time-varying policy parameters, and the solution of the dynamic programming of the final steady-state (in this case, the policy function).

The policy function for consumption and savings, $g_{c, t}(a, s)$ and $g_{a, t}(a, s)$, are arrays of dimension: $n_A \times n_S \times n_T$ (basically a matrix for each $t$). We set the policy functions of the final steady-state at the last period.

The Endogenous Grid Method works exactly as in the case of the stationary model. Nevertheless, one has to be careful to use the correct prices, $r_t$ and $w_t$. We start from the last period and loop backward.
md"""
##### Dynamic Programming Problem

We now write a function that solves for the policy functions for all $t$. The function takes as inputs the parameters, an interest rate and wage sequence, a sequence of possible time-varying policy parameters, and the solution of the dynamic programming of the final steady-state (in this case, the policy function).

The policy function for consumption and savings, $g_{c, t}(a, s)$ and $g_{a, t}(a, s)$, are arrays of dimension: $n_A \times n_S \times n_T$ (basically a matrix for each $t$). We set the policy functions of the final steady-state at the last period.

The Endogenous Grid Method works exactly as in the case of the stationary model. Nevertheless, one has to be careful to use the correct prices, $r_t$ and $w_t$. We start from the last period and loop backward.

"""
425 μs
RetrievePolicySequence (generic function with 1 method)
function RetrievePolicySequence(param, rseq, wseq, TransSeq, polSSend)
@unpack nA, nS, gA, gS, transS, β, γ, δ, nT , τkseq, τlseq = param

# Define parameters
gAA = repeat(gA, 1, nS)
gSS = repeat(gS', nA, 1) # transposed

# Pre-allocate value and policy funcctions
cpseq = zeros(nA, nS, nT)
apseq = zeros(nA, nS, nT)

apseq[:,:, nT] = polSSend.ap
cpseq[:,:, nT] = polSSend.cp


for it = nT-1:-1:2

expect = β*(1.0+rseq[it+1]*(1-τkseq[it+1]))* cpseq[:,:, it+1].^-γ*transS' # RHS of EE
c = expect.^(-1.0/γ)
endogAA = (c + gAA - wseq[it].*gSS*(1-τlseq[it]) .- TransSeq[it] )/(1+rseq[it]*(1-τkseq[it]) )

apNew = zeros(nA, nS)
for is = 1:nS
apNew[:, is] = linearGridIntp!(endogAA[:, is], gA, gA, apNew[:, is])

# Binding constraint
for ia = 1:nA
if apNew[ia, is] < gA[1] # binding constraint
apNew[ia, is] = gA[1]
else
break # exploit monotinicity of ia.
end
end
end

apseq[:,:, it] = apNew
4.7 ms
Distribution Sequence

The distribution sequence is an array of dimension $n_A \times n_S \times n_T$, basicallly a histogram for each $t$.

The function is similar to the one used for the stationary equilibrium. It takes as arguments parameters, the policy function sequence, and the distribution of the initial steady state.

Starting from the initial distribution, we loop forward and construct the next period distribution from the exogenous Markov chain and the policy function. In every loop, we first ``discretize'' the policy decision between two histogram bins (a la Young (2005)), then compute the mass of agents in every bin for the next period using the current distribution and the Markov-chain of the labor endowment process.
md"""
##### Distribution Sequence

The distribution sequence is an array of dimension $n_A \times n_S \times n_T$, basicallly a histogram for each $t$.

The function is similar to the one used for the stationary equilibrium. It takes as arguments parameters, the policy function sequence, and the distribution of the initial steady state.

Starting from the initial distribution, we loop forward and construct the next period distribution from the exogenous Markov chain and the policy function. In every loop, we first ``discretize'' the policy decision between two histogram bins (a la Young (2005)), then compute the mass of agents in every bin for the next period using the current distribution and the Markov-chain of the labor endowment process.


"""
393 μs
solveDistSeq (generic function with 1 method)
function solveDistSeq(param, apseq, dsnbeg)
@unpack nS, nA, gA, transS, nT = param

dsnseq = zeros(nA, nS, nT) # sequence of histograms
dsnseq[:, :, 1] = dsnbeg # initial steady-state
dsnseq[:, :, 2] = dsnbeg # shock hits at this period, but distribution depends of the policies from last period, so still the same.


for it = 2:nT-2

# === 1. RETRIEVE GRID AND WEIGHT OF INTERPOLATED POLICIES ============== #
ibelow = fill(0, nA, nS)
iweight = zeros(nA, nS)
for is = 1:nS
iweight[:, is], ibelow[:, is] = coordGridIntp!(gA, apseq[:, is, it], ibelow[:, is], iweight[:, is], robust = true)
end

# === 2. COMPUTE NEXT DISTRIBUTION GIVEN PREVIOUS DISTRIBUTION ========== #
dsnNew = zeros(nA, nS)
for ia = 1:nA, is = 1:nS
if dsnseq[ia, is, it] > 0.0
dsnNew[ibelow[ia, is], is] += iweight[ia, is] * dsnseq[ia, is, it]
dsnNew[ibelow[ia, is] + 1, is] += (1-iweight[ia, is]) * dsnseq[ia, is, it]
end
end
dsnseq[:, :, it+1] = dsnNew*transS # apply the markov-chain of labor process

end

return dsnseq
end

4.5 ms
Solving for the Equilibrium

Finally, we define a function that computes the excess demand for all periods $t$, and a function to solve for the equilibrium.

To compute the excess demand, we take as inputs the solution of the initial and final steady-state and a sequence of capital.

Using the sequence of capital and the time-varying $Z_t$, we use the first-order conditions of the firm problem to retrieve the sequence of prices, $r_t$ and $w_t$. Given the prices, we also check for the government budget constraint (if necessary) and get the lump-sum transfer that adjusts the budget. Now that we have the price and policy sequences, we solve the dynamic programming problem and simulate the distributions. Finally, we aggregate the distribution and compute a new sequence of capital.

If the sequence of capital is close to the guess, we stop. Otherwise, we update our guess using a convex combination between the new sequence and the guess.
md"""
##### Solving for the Equilibrium

Finally, we define a function that computes the excess demand for all periods $t$, and a function to solve for the equilibrium.

To compute the excess demand, we take as inputs the solution of the initial and final steady-state and a sequence of capital.

Using the sequence of capital and the time-varying $Z_t$, we use the first-order conditions of the firm problem to retrieve the sequence of prices, $r_t$ and $w_t$. Given the prices, we also check for the government budget constraint (if necessary) and get the lump-sum transfer that adjusts the budget. Now that we have the price and policy sequences, we solve the dynamic programming problem and simulate the distributions. Finally, we aggregate the distribution and compute a new sequence of capital.

If the sequence of capital is close to the guess, we stop. Otherwise, we update our guess using a convex combination between the new sequence and the guess.


"""
806 μs
ExcessDemandTrans (generic function with 1 method)
function ExcessDemandTrans(param, solSS0, solSS1, Kseq)
@unpack gA, nS, Lbar, α, δ, nT, zShock, Gseq, τlseq, τkseq = param

# Price sequences
rseq = zeros(nT) # first and last are SS
wseq = zeros(nT)

rseq[1] = solSS0[4]; rseq[nT] = solSS1[4]
wseq[1] = solSS0[3]; wseq[nT] = solSS1[3]

for it = 2:nT-1
rseq[it] = zShock[it]*α*(Kseq[it]/Lbar)^(α-1) - δ
wseq[it] = zShock[it]*(1-α)*(Kseq[it]/Lbar)^α
end

# Policy Sequence:
Y = zShock.*Kseq.^α*Lbar^(1-α) # total GDP
aggG = Y.*Gseq # compute aggregate government expenditure
TransSeq = @. rseq*Kseq*τkseq + wseq*Lbar*τlseq - aggG #transfer adjust

# Retrieve Policy Sequence
apseq, cpseq = RetrievePolicySequence(param, rseq, wseq, TransSeq, solSS1[1])
apseq[:, :, 1] = solSS0[1].ap
cpseq[:, :, 1] = solSS0[1].cp

# Retrieve Invariant Distribution
dsnseq = solveDistSeq(param, apseq, solSS0[2])
dsnseq[:, :, nT] = solSS1[2]

# Aggregate Asset Supply
Kseq_next = zeros(nT); Kseq_next[1] = solSS0[6]; Kseq_next[2] = solSS0[6]; Kseq_next[nT] = solSS1[6]

for it = 3:nT-1
Kseq_next[it] = sum(dsnseq[:,:,it].*repeat(gA, 1, nS) ) # get aggregate asset supply
end

49.9 ms
solveModelTrans (generic function with 1 method)
# Solve for the transition
function solveModelTrans(param)

nT = param.nT
# ===== SOLVE FIRST STEADY STATE ======= #
timevarParam = [param.zShock[1], param.τlseq[1], param.τkseq[1], param.Gseq[1]]
solSS0 = ModelSolutionSS(param, timevarParam)

# ===== SOLVE FINAl STEADY STATE ======= #
timevarParam = [param.zShock[nT], param.τlseq[nT], param.τkseq[nT], param.Gseq[nT]]
solSS1 = ModelSolutionSS(param, timevarParam)

# ===== ITERATE TO FIND EQ ======= #
Kseq = zeros(nT);
# guess: linear increase/decrease between the two SS
Kseq[1] = solSS0[6]; Kseq[nT] = solSS1[6]
d = (Kseq[nT] - Kseq[1])/(nT-1)
for i = 2:nT-1
Kseq[i] = Kseq[1] + d*(i-1)
end

maxiter = 5000
tol = 1E-7
damp = 0.2 # damp parameter

for iter = 1:maxiter
Kseq_next, = ExcessDemandTrans(param, solSS0, solSS1, Kseq)

d = maximum(abs.(Kseq_next - Kseq)) # ===== Check if has converged =====
Kseq = damp*Kseq_next + (1-damp)*Kseq # update guess with dampen parameter

if d < tol
println("Tol. achieved: $d")
break
end

2.5 msres
res = solveModelTrans(setPar()) # zshock

❔❔
Interest Rate Guess: 6.7815769938924015
Tol. achieved: 9.947598300641403e-14
Max iterations achieved. Invariant distribution did not converge. Tol 8.909360204754566e-8

Excess Demand: 96.1397276103629

Interest Rate Guess: 16.2701772192957
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 3.439927856452307e-11

Excess Demand: -16.2701772192956

Interest Rate Guess: 14.896799958063223
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.247913025645094e-11

Excess Demand: -14.89679995806296

Interest Rate Guess: 10.229976972924012
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 8.384340444145266e-11

Excess Demand: -10.226492237115986

Interest Rate Guess: 7.948282740177204
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.427261749017646e-11

Excess Demand: -7.8608968681130325

Interest Rate Guess: 7.364929867034802
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.659049948762544e-11

Excess Demand: -7.061141391938606

Interest Rate Guess: 7.073253430463602
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.957469571109101e-11

Excess Demand: -6.278015754218829

Interest Rate Guess: 6.927415212178001
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.961940300451388e-11

Excess Demand: -5.111187462042691

Interest Rate Guess: 6.854496103035201
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99667484208322e-11

Excess Demand: -2.9154107967896126

Interest Rate Guess: 6.818036548463801
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.996253911431774e-11

Excess Demand: 1.4938346092149244

Interest Rate Guess: 6.830388898909332
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994309720096073e-11

Excess Demand: -0.7371822913077688

Interest Rate Guess: 6.826307381641426
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99492155706605e-11

Excess Demand: -0.135175145008688

Interest Rate Guess: 6.825466963017688
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99956792716028e-11

Excess Demand: 0.002618784776688976

Interest Rate Guess: 6.825482935241583
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11

Excess Demand: -4.899470898589442e-5
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11

Interest Rate Guess: 6.7815769938924015
Tol. achieved: 9.947598300641403e-14
Max iterations achieved. Invariant distribution did not converge. Tol 8.909360204754566e-8

Excess Demand: 96.1397276103629

Interest Rate Guess: 16.2701772192957
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 3.439927856452307e-11

Excess Demand: -16.2701772192956

Interest Rate Guess: 14.896799958063223
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.247913025645094e-11

Excess Demand: -14.89679995806296

Interest Rate Guess: 10.229976972924012
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 8.384340444145266e-11

Excess Demand: -10.226492237115986

Interest Rate Guess: 7.948282740177204
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.427261749017646e-11

Excess Demand: -7.8608968681130325

Interest Rate Guess: 7.364929867034802
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.659049948762544e-11

Excess Demand: -7.061141391938606

Interest Rate Guess: 7.073253430463602
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.957469571109101e-11

Excess Demand: -6.278015754218829

Interest Rate Guess: 6.927415212178001
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.961940300451388e-11

Excess Demand: -5.111187462042691

Interest Rate Guess: 6.854496103035201
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99667484208322e-11

Excess Demand: -2.9154107967896126

Interest Rate Guess: 6.818036548463801
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.996253911431774e-11

Excess Demand: 1.4938346092149244

Interest Rate Guess: 6.830388898909332
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994309720096073e-11

Excess Demand: -0.7371822913077688

Interest Rate Guess: 6.826307381641426
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99492155706605e-11

Excess Demand: -0.135175145008688

Interest Rate Guess: 6.825466963017688
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99956792716028e-11

Excess Demand: 0.002618784776688976

Interest Rate Guess: 6.825482935241583
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11

Excess Demand: -4.899470898589442e-5
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11
Iter: 1
Tol: 0.2970231411907909
Iter: 2
Tol: 0.6525127611008177
Iter: 3
Tol: 0.6343309020866901
Iter: 4
Tol: 0.4400920361989762
Iter: 5
Tol: 0.35663046906266427
Iter: 6
Tol: 0.28421519583596666
Iter: 7
Tol: 0.21642352497995265
Iter: 8
Tol: 0.17630290163553575
Iter: 9
Tol: 0.13433225259457515
Iter: 10
Tol: 0.10853632579252359
Iter: 11
Tol: 0.0835877722098477
Iter: 12
Tol: 0.06687098435772842
Iter: 13
Tol: 0.05192561907601778
Iter: 14
Tol: 0.04127964167367626
Iter: 15
Tol: 0.03220153269912451
Iter: 16
Tol: 0.025517893778433276
Iter: 17
Tol: 0.019948294147966728
Iter: 18
Tol: 0.015786867057212106
Iter: 19
Tol: 0.012350920571252644
Iter: 20
Tol: 0.009770260911588835
Iter: 21
Tol: 0.007645328601490142
Iter: 22
Tol: 0.006047450128869869
Iter: 23
Tol: 0.004732245112109901
Iter: 24
Tol: 0.003743219547346577
Iter: 25
Tol: 0.002929157643086988
Iter: 26
Tol: 0.0023168997037741335
Iter: 27
Tol: 0.0018131420116329622
Iter: 28
Tol: 0.0014340196966573515
Iter: 29
Tol: 0.0011223647351030763
Iter: 30
Tol: 0.0008875485376211856
Iter: 31
Tol: 0.0006947769240586155
Iter: 32
Tol: 0.0005493169301153955
Iter: 33
Tol: 0.00043009199047006064
Iter: 34
Tol: 0.00033997902367488564
Iter: 35
Tol: 0.0002662429503050845
Iter: 36
Tol: 0.0002104179897726155
Iter: 37
Tol: 0.0001648135376655091
Iter: 38
Tol: 0.00013023191114580612
Iter: 39
Tol: 0.00010202437888384708
Iter: 40
Tol: 8.060407525789515e-5
Iter: 41
Tol: 6.315540251122798e-5
Iter: 42
Tol: 4.988869361444159e-5
Iter: 43
Tol: 3.909419726610963e-5
Iter: 44
Tol: 3.087825869307892e-5
Iter: 45
Tol: 2.4199673888780637e-5
Iter: 46
Tol: 1.9112111561803147e-5
Iter: 47
Tol: 1.4979680252658056e-5
Iter: 48
Tol: 1.1829577477584508e-5
Iter: 49
Tol: 9.272395232962083e-6
Iter: 50
Tol: 7.322070658233315e-6
Iter: 51
Tol: 5.739555398065477e-6
Iter: 52
Tol: 4.532127663381402e-6
Iter: 53
Tol: 3.552728665923155e-6
Iter: 54
Tol: 2.805261560823169e-6
Iter: 55
Tol: 2.1990941050376023e-6
Iter: 56
Tol: 1.7363894855293438e-6
Iter: 57
Tol: 1.361206632743972e-6
Iter: 58
Tol: 1.074788771759927e-6
Iter: 59
Tol: 8.425644217879835e-7
Iter: 60
Tol: 6.652744781732167e-7
Iter: 61
Tol: 5.215322245177845e-7
Iter: 62
Tol: 4.117940077819071e-7
Iter: 63
Tol: 3.2281856121585406e-7
Iter: 64
Tol: 2.548945987612683e-7
Iter: 65
Tol: 1.9981840893024128e-7
Iter: 66
Tol: 1.5777658646243253e-7
Iter: 67
Tol: 1.2368366952841825e-7
Tol. achieved: 9.766202602179419e-8
25.6 s

After we have computed the equilibrium, we can compute some statistics and plot the impulse response function.
md"""

After we have computed the equilibrium, we can compute some statistics and plot the impulse response function.
"""
203 μs
begin
@unpack zShock, nT, gA, nS = param
gT = collect(range(1,nT, length = nT))
# function that calculates gini (see wikipedia formula):
function calGini2(dist, values)
cumulative = cumsum(dist.*values)/sum(dist.*values) # cum distribution in %
B = sum(cumulative.*dist) # total area below lorenz curve
# with cns distribution this should be 0.5, but since it is a histogram it may deviate a bit
AreaBelow45 = sum(dist.*(cumsum(dist)./sum(dist))) # A + B
gini = (AreaBelow45 - B)/AreaBelow45
return gini
end
gini = zeros(nT,3)
share_const = zeros(nT)
for it = 1:nT # compute stats for all t
dsn = res[2][:,:,it] # distribution at t
asset_distribution = sum(dsn, dims = 2)[:]
gini[it,1] = calGini2(asset_distribution, gA) # gini everybody
gini[it,2] = calGini2(dsn[:,1]./sum(dsn[:,1]), gA) # gini income poor
gini[it,3] = calGini2(dsn[:,nS]./sum(dsn[:,nS]), gA) # gini income rich
share_const[it] = dsn[1]
end

# compute variables in % and in p.p. for impulse response function
zper=log.(zShock)*100
Kper = @. (log(res[1]) - log(res[1][1]))*100
rpp = @. (res[6] - res[6][1])*100
wper = @. (log(res[5]) - log(res[5][1]))*100

share_pp = (share_const .- share_const[1])*100
gini1_pp = ( gini[:,1] .- gini[1,1])*100
gini2_pp = ( gini[:,2] .- gini[1,2])*100
gini3_pp = ( gini[:,3] .- gini[1,3])*100
226 ms
begin # impulse response 1
p1 = plot(gT, zper, label="Z", lw = 3, ylabel = "%", xlabel = "t")
p2 = plot(gT, Kper, label="K", lw = 3, ylabel = "%", xlabel = "t")
p3 = plot(gT, rpp, label="r", lw = 3, ylabel = "p.p", xlabel = "t")
p4 = plot(gT, wper, label="w", lw = 3, ylabel = "%", xlabel = "t")
plot(p1, p2, p3, p4, layout = (2, 2))
end
4.9 ms
begin # impulse response 2
p5 = plot(gT, gini1_pp, title="Gini Wealth (All)", lw = 3, ylabel = "p.p", xlabel = "t")
p6 = plot(gT, gini2_pp, title="Gini Wealth (Smin)", lw = 3, ylabel = "p.p", xlabel = "t")
p7 = plot(gT, gini3_pp, title="Gini Wealth (Smax)", lw = 3, ylabel = "p.p", xlabel = "t")
p8 = plot(gT, share_pp, title="Share Constrained", lw = 3, ylabel = "p.p", xlabel = "t")
plot(p5, p6, p7, p8, layout = (2, 2), legends = false)
end
4.8 ms
Boppart-Krusell-Mitman

BKM have shown that we can use the IRF to simulate the model. We now use their procedure to simulate a sequence of capital. We simulate 500 periods.

The trick is to compute the capital at $t$ as linear function of past innovations:

$$K_t = \sum_{s=0}^\infty \varepsilon_{t-s}K_s ,$$

where $K_s$ is the response of capital today of shock that happened $s$ periods before. In practice, we do not actually use infinity, but a large enough number of periods such that the shock has faded out.

The beauty of this method is that, because of the linear assumption, $K_s$ is exactly what we have computed using the IRF.

Hence, to simulate a sequence of capital, we just need to simulate a sequence of innovations and take the linear product with the sequence of capital from the IRF.

md"""
##### Boppart-Krusell-Mitman

BKM have shown that we can use the IRF to simulate the model. We now use their procedure to simulate a sequence of capital. We simulate 500 periods.

The trick is to compute the capital at $t$ as linear function of past innovations:

$K_t = \sum_{s=0}^\infty \varepsilon_{t-s}K_s ,$

where $K_s$ is the response of capital today of shock that happened $s$ periods before. In practice, we do not actually use infinity, but a large enough number of periods such that the shock has faded out.

The beauty of this method is that, because of the linear assumption, $K_s$ is exactly what we have computed using the IRF.

Hence, to simulate a sequence of capital, we just need to simulate a sequence of innovations and take the linear product with the sequence of capital from the IRF.


"""
450 μs
begin
@unpack Zbar, ρz = param
simul = 500
Random.seed!(123)
Kirf = res[1][nT:-1:2] # ignore first period since it was SS and reverse the array
Kss = res[1][1]
nShocks = length(Kirf) # period of shocks (nT - 1)
Ksimul = zeros(simul)
#shocks = rand(Uniform(-1, 1), simul+nShocks)
shocks = rand(Normal(0, 1), simul+nShocks)
log_z_Simul = zeros(simul+nShocks)
log_z_Simul[1] = Zbar*ρz + shocks[1] # first period z is in the SS
for it = 2:(simul+nShocks)
log_z_Simul[it] = log_z_Simul[it-1]*ρz + shocks[it] # AR(1)
end
for it = 1:simul
# simulations are in deviations of SS
# recall that Kirf is reversed
Ksimul[it] = sum(shocks[it:it+nShocks-1].*(Kirf./Kss .-1))*100

end

gSimul = collect(range(1,simul,step=1))
pK = plot(gSimul, Ksimul, lw = 3, ylabel = "% deviations of SS", xlabel = "t", label = "K")
pZ = plot(gSimul, log_z_Simul[nShocks+1:nShocks+simul], lw = 3, ylabel = "% deviations of SS", xlabel = "t", label = "log Z")
plot(pK, pZ, layout = (2, 1))
end
4.6 ms
Transition between different Steady-states

As a final experiment, we solve for the deterministic transition between different steady-states. We assume that at $t=2$, the government decides to change the labor tax rate from 0 to 20%. The economy slowly converges to the new steady-state.

md"""
##### Transition between different Steady-states

As a final experiment, we solve for the deterministic transition between different steady-states. We assume that at $t=2$, the government decides to change the labor tax rate from 0 to 20%. The economy slowly converges to the new steady-state.
"""
249 μs
begin
param2 = setPar(sizeshock = 0.0, τl1 = 0.2)
res2 = solveModelTrans(param2) # Increase in labor tax // ps. must change initial guess to experiment with capital tax
end


❔❔
Interest Rate Guess: 6.7815769938924015
Tol. achieved: 9.947598300641403e-14
Max iterations achieved. Invariant distribution did not converge. Tol 8.909360204754566e-8

Excess Demand: 96.1397276103629

Interest Rate Guess: 16.2701772192957
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 3.439927856452307e-11

Excess Demand: -16.2701772192956

Interest Rate Guess: 14.896799958063223
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.247913025645094e-11

Excess Demand: -14.89679995806296

Interest Rate Guess: 10.229976972924012
Tol. achieved: 5.684341886080802e-14
Tol. achieved: 8.384340444145266e-11

Excess Demand: -10.226492237115986

Interest Rate Guess: 7.948282740177204
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.427261749017646e-11

Excess Demand: -7.8608968681130325

Interest Rate Guess: 7.364929867034802
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.659049948762544e-11

Excess Demand: -7.061141391938606

Interest Rate Guess: 7.073253430463602
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.957469571109101e-11

Excess Demand: -6.278015754218829

Interest Rate Guess: 6.927415212178001
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.961940300451388e-11

Excess Demand: -5.111187462042691

Interest Rate Guess: 6.854496103035201
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99667484208322e-11

Excess Demand: -2.9154107967896126

Interest Rate Guess: 6.818036548463801
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.996253911431774e-11

Excess Demand: 1.4938346092149244

Interest Rate Guess: 6.830388898909332
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994309720096073e-11

Excess Demand: -0.7371822913077688

Interest Rate Guess: 6.826307381641426
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99492155706605e-11

Excess Demand: -0.135175145008688

Interest Rate Guess: 6.825466963017688
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.99956792716028e-11

Excess Demand: 0.002618784776688976

Interest Rate Guess: 6.825482935241583
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11

Excess Demand: -4.899470898589442e-5
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.994291939180444e-11

Interest Rate Guess: 6.7815769938924015
Tol. achieved: 9.947598300641403e-14
Max iterations achieved. Invariant distribution did not converge. Tol 4.708622348015973e-8

Excess Demand: 79.46027589335081

Interest Rate Guess: 16.2701772192957
Tol. achieved: 2.842170943040401e-14
Tol. achieved: 2.3517838718123407e-11

Excess Demand: -16.27017721929559

Interest Rate Guess: 14.657511633677938
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 2.912802464505315e-11

Excess Demand: -14.657511633677823

Interest Rate Guess: 10.11033281073137
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.315739846194049e-11

Excess Demand: -10.109642530341421

Interest Rate Guess: 7.918371699629043
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.114914378827166e-11

Excess Demand: -7.878072168077321

Interest Rate Guess: 7.349974346760722
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.893347252543094e-11

Excess Demand: -7.182628188101754

Interest Rate Guess: 7.065775670326562
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.935967326679673e-11

Excess Demand: -6.588100728586724

Interest Rate Guess: 6.923676332109482
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.977827591933774e-11

Excess Demand: -5.790377801771859

Interest Rate Guess: 6.852626663000941
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.991876076531625e-11

Excess Demand: -4.353908828712924

Interest Rate Guess: 6.817101828446671
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.998236266683946e-11

Excess Demand: -1.4492111226370321

Interest Rate Guess: 6.800265445528016
Tol. achieved: 8.526512829121202e-14
Max iterations achieved. Invariant distribution did not converge. Tol 2.118685086312938e-10

Excess Demand: 3.8375202926622904

Interest Rate Guess: 6.808683636987343
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.996429552183717e-11

Excess Demand: 0.37917572231460905

Interest Rate Guess: 6.810429423901084
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.997582796350546e-11

Excess Demand: -0.08799884576473271

Interest Rate Guess: 6.810100580586706
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.998254047599575e-11

Excess Demand: -0.004365157259480945

Interest Rate Guess: 6.810083612357462
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.998973090480368e-11

Excess Demand: 2.689066196381873e-6
Tol. achieved: 8.526512829121202e-14
Tol. achieved: 9.998973090480368e-11
Iter: 1
Tol: 0.0982518416599083
Iter: 2
Tol: 0.1949876578293459
Iter: 3
Tol: 0.1758245031766421
Iter: 4
Tol: 0.13912320804336886
Iter: 5
Tol: 0.11991650148743016
Iter: 6
Tol: 0.10022398243463826
Iter: 7
Tol: 0.08414141823639198
Iter: 8
Tol: 0.07095048325955045
Iter: 9
Tol: 0.05947602654959905
Iter: 10
Tol: 0.05008456313488452
Iter: 11
Tol: 0.04206765649741229
Iter: 12
Tol: 0.03536180017791235
Iter: 13
Tol: 0.029741017434847805
Iter: 14
Tol: 0.024978364208545045
Iter: 15
Tol: 0.021018372564653554
Iter: 16
Tol: 0.01764869457450491
Iter: 17
Tol: 0.01485116114857199
Iter: 18
Tol: 0.012471241869954852
Iter: 19
Tol: 0.010492823592223921
Iter: 20
Tol: 0.008812908119161023
Iter: 21
Tol: 0.0074134577394655565
Iter: 22
Tol: 0.006227694679542317
Iter: 23
Tol: 0.005237885543085241
Iter: 24
Tol: 0.004400763878907199
Iter: 25
Tol: 0.003700828622651642
Iter: 26
Tol: 0.0031097214794595374
Iter: 27
Tol: 0.002614862679585883
Iter: 28
Tol: 0.002197399126973032
Iter: 29
Tol: 0.001847584273305003
Iter: 30
Tol: 0.0015527157692316607
Iter: 31
Tol: 0.001305460260104141
Iter: 32
Tol: 0.001097164496895786
Iter: 33
Tol: 0.0009224140220682742
Iter: 34
Tol: 0.0007752633604427572
Iter: 35
Tol: 0.000651763641491776
Iter: 36
Tol: 0.0005478039466666473
Iter: 37
Tol: 0.000460527712354164
Iter: 38
Tol: 0.0003870792411078128
Iter: 39
Tol: 0.00032540363824029583
Iter: 40
Tol: 0.0002735103257904825
Iter: 41
Tol: 0.0002299268942929089
Iter: 42
Tol: 0.00019326219433768443
Iter: 43
Tol: 0.00016246421125121202
Iter: 44
Tol: 0.00013655876382490106
Iter: 45
Tol: 0.00011479582269924293
Iter: 46
Tol: 9.649211960294224e-5
Iter: 47
Tol: 8.111381371112003e-5
Iter: 48
Tol: 6.818106718675665e-5
Iter: 49
Tol: 5.7314413309583756e-5
Iter: 50
Tol: 4.817652649791171e-5
Iter: 51
Tol: 4.0497957405349894e-5
Iter: 52
Tol: 3.404136441620409e-5
Iter: 53
Tol: 2.8615581983970628e-5
Iter: 54
Tol: 2.4053498520970606e-5
Iter: 55
Tol: 2.0219583456437817e-5
Iter: 56
Tol: 1.6996104608146823e-5
Iter: 57
Tol: 1.428702989691999e-5
Iter: 58
Tol: 1.200937495582366e-5
Iter: 59
Tol: 1.0095127455755915e-5
Iter: 60
Tol: 8.48577129097805e-6
Iter: 61
Tol: 7.133156682215258e-6
Iter: 62
Tol: 5.996007152120342e-6
Iter: 63
Tol: 5.040246420229266e-6
Iter: 64
Tol: 4.236750655195465e-6
Iter: 65
Tol: 3.5614091347113686e-6
Iter: 66
Tol: 2.9936676639863435e-6
Iter: 67
Tol: 2.5164712837977277e-6
Iter: 68
Tol: 2.115311096950734e-6
Iter: 69
Tol: 1.7781246572923237e-6
Iter: 70
Tol: 1.4946686706096557e-6
Iter: 71
Tol: 1.2564132951808915e-6
Iter: 72
Tol: 1.056125734422153e-6
Iter: 73
Tol: 8.877749690938685e-7
Iter: 74
Tol: 7.46253378203221e-7
Iter: 75
Tol: 6.272968811060764e-7
Iter: 76
Tol: 5.272985674764641e-7
Iter: 77
Tol: 4.432443398982855e-7
Iter: 78
Tol: 3.725863946257846e-7
Iter: 79
Tol: 3.1319393212214663e-7
Iter: 80
Tol: 2.632675162317355e-7
Iter: 81
Tol: 2.2130108412454774e-7
Iter: 82
Tol: 1.8602344642459911e-7
Iter: 83
Tol: 1.5637021011372099e-7
Iter: 84
Tol: 1.3144321986402474e-7
Iter: 85
Tol: 1.1049013526331919e-7
Tol. achieved: 9.287688929049409e-8
35.4 s
begin
# compute variables in % and in p.p. for impulse response function
τlseq =param2.τlseq*100
Kper2 = @. (log(res2[1]) - log(res2[1][1]))*100
rpp2 = @. (res2[6] - res2[6][1])*100
wper2 = @. (log(res2[5]) - log(res2[5][1]))*100

p21 = plot(gT, τlseq, title="Labor Tax", lw = 3, ylabel = "%", xlabel = "t")
p22 = plot(gT, Kper2, title="K", lw = 3, ylabel = "%", xlabel = "t")
p23 = plot(gT, rpp2, title="r", lw = 3, ylabel = "p.p", xlabel = "t")
p24 = plot(gT, wper2, title="w", lw = 3, ylabel = "%", xlabel = "t")
plot(p21, p22, p23, p24, layout = (2, 2), legends = false)

end
4.7 ms