Solutions for http://quant-econ.net/jl/arellano.html
In [3]:
using QuantEcon, QuantEcon.Models
using Gadfly, Compose, ColorTypes, DataFrames
Compute the value function, policy and equilibrium prices
In [4]:
ae = ArellanoEconomy(β=.953, # time discount rate
γ=2., # risk aversion
r=0.017, # international interest rate
ρ=.945, # persistence in output
η=0.025, # st dev of output shock
θ=0.282, # prob of regaining access
ny=21, # number of points in y grid
nB=251) # number of points in B grid
# now solve the model on the grid.
vfi!(ae)
Compute the bond price schedule as seen in figure 3 of Arellano (2008)
In [6]:
# Create "Y High" and "Y Low" values as 5% devs from mean
high, low = mean(ae.ygrid)*1.05, mean(ae.ygrid)*.95
iy_high, iy_low = map(x->searchsortedfirst(ae.ygrid, x), (high, low))
# Extract a suitable plot grid
x = Float64[]
q_low = Float64[]
q_high = Float64[]
for i=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=repeat(x, outer=[2]), y=[q_low; q_high],
color=repeat([:Low, :High], inner=[length(x)]),
Guide.title("Bond price schedule q(y, B')"),
Guide.xlabel("B'"), Guide.ylabel("q"),
Guide.colorkey("y"), Geom.line)
Out[6]:
Draw a plot of the value functions
In [7]:
plot(x=repeat(ae.Bgrid, outer=[2]),
y=vec(ae.vf[:, [iy_low, iy_high]]),
color=repeat([:Low, :High], inner=[length(ae.Bgrid)]),
Guide.title("Value functions"),
Guide.xlabel("B"), Guide.ylabel("V(y,B)"),
Guide.colorkey("y"), Geom.line)
Out[7]:
Draw a heat map for default probability
In [8]:
plot(x_min=repeat(ae.Bgrid[1:end-1], inner=[ae.ny-1]),
x_max=repeat(ae.Bgrid[2:end], inner=[ae.ny-1]),
y_min=repeat(ae.ygrid[1:end-1], outer=[ae.nB-1]),
y_max=repeat(ae.ygrid[2:end], outer=[ae.nB-1]),
color=clamp(vec(ae.defprob[1:end-1, 1:end-1]'), 0, 1),
Geom.rectbin,
Guide.xlabel("B'"), Guide.ylabel("y"),
Guide.title("Probability of default"), Geom.rectbin,
Scale.y_continuous(minvalue=0.8, maxvalue=1.2),
Scale.x_continuous(minvalue=minimum(ae.Bgrid), maxvalue=0.0),
Scale.color_continuous(minvalue=0, maxvalue=1))
Out[8]:
Plot a time series of major variables simulated from the model.
In [31]:
# simulate
T = 250
y_vec, B_vec, q_vec, default_vec = simulate(ae, T)
# find starting and ending periods of recessions (if any)
if any(default_vec)
defs = find(default_vec)
def_breaks = diff(defs) .> 1
def_start = defs[[true; def_breaks]]
def_end = defs[[def_breaks; true]]
else
starts, ends = Int[], Int[]
end
# construct boxes that shade periods of default
def_box = Guide.annotation(compose(context(),
[rectangle(i[1], 0h, i[2]-i[1], 1h)
for i=zip(def_start, def_end)]...,
fill(RGBA(0.5, 0.5, 0.5, 0.2))))
# xy labels are common for all plots
xy_lab = [Guide.xlabel("time"), Guide.ylabel("")]
# now iterate over three variables and put them into an array
plots = Gadfly.Plot[]
for (vec, name) in [(y_vec, "Output"), (B_vec, "Foreign assets"), (q_vec, "Bond price")]
push!(plots,
plot(x=1:T, y=vec, Geom.line, def_box, Guide.title(name), xy_lab...))
end
# set final plot height and vertically stack the above three plots
set_default_plot_size(6inch, 8inch)
vstack(plots...)
Out[31]: