quant-econ Solutions: Default Risk and Income Fluctuations

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)









    



Finished iteration 25 with dist of 0.3424484168091517
Finished iteration 50 with dist of 0.0982039407428843
Finished iteration 75 with dist of 0.029158662291511206
Finished iteration 100 with dist of 0.008729266837647742
Finished iteration 125 with dist of 0.002618400938121823
Finished iteration 150 with dist of 0.0007857709211798181
Finished iteration 175 with dist of 0.0002358324600919559
Finished iteration 200 with dist of 7.078195654486308e-5
Finished iteration 225 with dist of 2.1244388765495614e-5
Finished iteration 250 with dist of 6.376267929653068e-6
Finished iteration 275 with dist of 1.9137668516577833e-6
Finished iteration 300 with dist of 5.743961750681592e-7
Finished iteration 325 with dist of 1.7239873884022927e-7
Finished iteration 350 with dist of 5.174360495630026e-8
Finished iteration 375 with dist of 1.5530293495658043e-8

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]: