-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathInflation_with_maturity.jl
209 lines (171 loc) · 7.19 KB
/
Inflation_with_maturity.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
using LinearAlgebra, Statistics
using LaTeXStrings, QuantEcon, DataFrames, Plots, Random
## Based on QuantEcon lecture. Extended for long term debt, exogenous inflation and maturity choice.
function ArellanoEconomy(; beta = 0.953,
gamma = 2.0,
r = 0.017,
lambda = 0.9, z = 0.01,
rho = 0.945,
eta = 0.025,
theta = 0.282,
ny = 21,
nB = 251, nm = 10)
# create grids
Bgrid = collect(range(-0.4, 0.4, length = nB))
mc = tauchen(ny, rho, eta)
Pi = mc.p
ygrid = exp.(mc.state_values)
ydefgrid = min.(0.969 * mean(ygrid), ygrid)
infl_grid = LinRange(1.2, 1, ny) # exogenous inflation grid DECREASES with the state.
# infl_grid = LinRange(1, 1.2, ny) # inflation increases with the state
Pi_m = (1/nm) .* ones(nm, nm)
m_min = -0.01
m_max = 0.01
m_grid = LinRange(m_min, m_max, nm)
# define value functions
# notice ordered different than Python to take
# advantage of column major layout of Julia
vf = zeros(nB, ny, nm)
vd = zeros(1, ny, nm)
vc = zeros(nB, ny, nm)
policy = zeros(nB, ny, nm)
q = ones(nB, ny) .* (1 / (1 + r))
defprob = zeros(nB, ny, nm)
return (; beta, gamma, r, lambda, z, rho, eta, theta, ny,
nB, ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc,
policy, q, defprob, m_grid, nm, infl_grid)
end
u(ae, c) = c^(1 - ae.gamma) / (1 - ae.gamma)
function one_step_update!(ae,
EV,
EVd,
EVc)
# unpack stuff
(; beta, gamma, r, lambda, z, rho, eta, theta, ny, nB) = ae
(; ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob, m_grid, nm, infl_grid) = ae
zero_ind = searchsortedfirst(Bgrid, 0.0)
for im in 1:nm
m = ae.m_grid[im]
m_min = ae.m_grid[1]
for iy in 1:ny
y = ae.ygrid[iy]
ydef = ae.ydefgrid[iy]
infl = ae.infl_grid[iy]
# value of being in default with income y
defval = u(ae, ydef + m_min) +
beta * (theta * EVc[zero_ind, iy] + (1 - theta) * EVd[1, iy])
ae.vd[1, iy, im] = defval
for ib in 1:nB
B = ae.Bgrid[ib]
current_max = -1e14
pol_ind = 0
for ib_next in 1:nB
c = max(y + m - ae.q[ib_next, iy] * (Bgrid[ib_next] - (1-lambda)*B/infl) + (lambda + (1-lambda)*z)*B/infl, 1e-14)
mm = u(ae, c) + beta * EV[ib_next, iy]
if mm > current_max
current_max = mm
pol_ind = ib_next
end
end
# update value and policy functions
ae.vc[ib, iy, im] = current_max
ae.policy[ib, iy, im] = Int(pol_ind)
ae.vf[ib, iy, im] = defval > current_max ? defval : current_max
end
end
end
end
function compute_prices!(ae)
# unpack parameters
(; beta, gamma, r, lambda, z, rho, eta, theta, ny, nB) = ae
(; ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob, m_grid, nm, infl_grid) = ae
# create default values with a matching size
vd_compat = repeat(ae.vd, nB, 1, 1)
default_states = vd_compat .> ae.vc
# update default probabilities and prices
q_temp = zeros(nB, ny, nm)
for im in 1:nm
ae.defprob[:,:,im] .= default_states[:,:,im] * ae.Pi'
policies = Int.(ae.policy[:, :, im])
# Initialize an array to store the diagonal values
q_new_d = [ae.q[policies[ib, iy], iy] for ib in 1:nB, iy in 1:ny]
q_new_d = reshape(q_new_d, nB, ny) # Reshape the array to nB x ny
infl_mat = infl_grid .* ones(1, nB) # ny x nB matrix
infl_mat = infl_mat' # nB x ny matrix
q_temp[:,:,im] .= ( (1 ./ infl_mat) .* (1 .- default_states[:,:,im]) ## Double check that my infl_mat elementwise multiplication makes sense
.* (lambda .+ (1 .- lambda) .* (z .+ q_new_d) ) * ae.Pi' ) ## as an extension of equation (6) in Chatterjee and Eyigungor (2012)
end
# average over m
ae.q .= mean(q_temp, dims = 3) / (1 + r)
return
end
function vfi!(ae; tol = 1e-5, maxit = 10000)
# unpack stuff
(; beta, gamma, r, lambda, z, rho, eta, theta, ny, nB) = ae
(; ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob, m_grid, nm, infl_grid) = ae
Pit = Pi'
# Iteration stuff
it = 0
dist = 10.0
# allocate memory for update
V_upd = similar(ae.vf)
q_upd = similar(ae.q)
while dist > tol && it < maxit
it += 1
# compute expectations for this iterations
# (we need Pi' because of order value function dimensions)
copyto!(V_upd, ae.vf)
copyto!(q_upd, ae.q)
EV = dropdims(mean(ae.vf, dims = 3), dims = 3) * Pit
EVd = dropdims(mean(ae.vd, dims = 3), dims = 3) * Pit
EVc = dropdims(mean(ae.vc, dims = 3), dims = 3) * Pit
# update value function
one_step_update!(ae, EV, EVd, EVc)
# update prices
compute_prices!(ae)
dist = maximum(abs(x - y) for (x, y) in zip(V_upd, ae.vf))
dist2 = maximum(abs(x-y) for (x,y) in zip(q_upd, ae.q))
if it % 25 == 0
println("Finished iteration $(it) with dist of $(dist)")
end
end
end
ae = ArellanoEconomy(beta = 0.953, # time discount rate
gamma = 2.0, # risk aversion
r = 0.017, # international interest rate
lambda = 0.9, z = 0.01,
rho = 0.945, # persistence in output
eta = 0.025, # st dev of output shock
theta = 0.282, # prob of regaining access
ny = 21, # number of points in y grid
nB = 251,
nm = 5) # number of points in m grid
# now solve the model on the grid.
vfi!(ae)
# create "Y High" and "Y Low" values as 5% devs from mean
high, low = 1.05 * mean(ae.ygrid), 0.95 * mean(ae.ygrid)
iy_high, iy_low = map(x -> searchsortedfirst(ae.ygrid, x), (high, low))
# extract a suitable plot grid
x = zeros(0)
q_low = zeros(0)
q_high = zeros(0)
for i in 1:(ae.nB)
b = ae.Bgrid[i]
if -0.35 <= b <= 0 # to match fig 3 of Arellano
push!(x, b)
push!(q_low, ae.q[i, iy_low])
push!(q_high, ae.q[i, iy_high])
end
end
# generate plot
plot(x, q_low, label = "Low")
plot!(x, q_high, label = "High")
plot!(title = L"Bond price schedule $q(y, B^\prime)$",
xlabel = L"B^\prime", ylabel = L"q", legend_title = L"y",
legend = :topleft)
heatmap(ae.Bgrid[1:(end - 1)],
ae.ygrid[2:end],
reshape(clamp.(vec( mean(ae.defprob[1:(end - 1), 1:(end - 1),:], dims = 3) ), 0, 1), 250,
20)')
plot!(xlabel = L"B^\prime", ylabel = L"y", title = "Probability of default",
legend = :topleft)