Prepared for WAMS 2017 by John Stachurski
Code for studying stability and computing solutions to the recursion that defines lifetime utility in Epstein-Zin preference models, with nonstationary consumption.
In [1]:
import QuantEcon: rouwenhorst, MarkovChain, LinInterp
In [2]:
"""
Function to compute spectral radius of a matrix.
"""
compute_spec_rad(Q::Matrix) = maximum(abs, eigvals(Q))
Out[2]:
Epstein-Zin utility specification.
In [3]:
"""
Struct to store parameters of Epstein-Zin model.
"""
mutable struct EpsteinZin{T <: AbstractFloat}
ψ::T # Elasticity of intertemporal substitution
γ::T # Risk aversion parameter
β::T # Time discount factor
ζ::T # Preference factor, current consumption
θ::T # Derived parameter
end
"""
Simple EpsteinZin constructor.
"""
function EpsteinZin(ψ, γ, β)
ζ = 1 - β
θ = (1 - γ) / (1 - 1/ψ)
return EpsteinZin(ψ, γ, β, ζ, θ)
end
"""
EpsteinZin constructor for BY model.
"""
function EpsteinZinBY(; ψ=1.5, γ=10.0, β=0.998)
return EpsteinZin(ψ, γ, β)
end
"""
EpsteinZin constructor for BKY model.
"""
function EpsteinZinBKY(; ψ=1.5, γ=10.0, β=0.9989)
return EpsteinZin(ψ, γ, β)
end
Out[3]:
Convenience functions for AR1 models
In [4]:
"""
Struct to store parameters of AR1 model
X' = ρ X + b + σ W
"""
mutable struct AR1{T <: AbstractFloat}
ρ::T # Correlation coefficient
b::T # Intercept
σ::T # Volatility
end
"""
A constructor for AR1 with default values from Pohl, Schmedders and Wilms.
"""
function AR1(; ρ=0.91, b=0.0, σ=0.0343)
return AR1(ρ, b, σ)
end
"""
Convert a Gaussian AR1 process to a Markov Chain via Rouwenhorst's method.
"""
function ar1_to_mc(ar1::AR1, M::Integer)
return rouwenhorst(M, ar1.ρ, ar1.σ, ar1.b)
end
Out[4]:
Discretizing the SV model using two iterations of Rouwenhorst.
The model is
$$ z' = ρ z + s_z σ e' $$$$ (σ^2)' = v σ^2 + d + s_σ w' $$where $\{e\}$ and $\{w\}$ are IID and $N(0, 1)$.
In [5]:
"""
Two convenience functions to switch between a single index m that
points to elements of a matrix A_{ij} and the multi-index (ij). The
rule taking us from ij to m is that we start at row 1 and keep counting,
starting at the first row and working down. Thus,
m = (i - 1) * J + j
"""
multi_from_single(m, J) = div(m - 1, J) + 1, rem(m - 1, J) + 1
single_from_multi(i, j, J) = (i - 1) * J + j
"""
Struct for parameters of the SV model
z' = ρ z + s_z σ e'
(σ^2)' = v σ^2 + d + s_σ w'
"""
mutable struct StochasticVolatility{T <: AbstractFloat}
ρ::T
s_z::T
v::T
d::T
s_σ::T
end
"""
Discretize the SV model defined above. Returns a (2, M) matrix
x_states, each element x of which is a pair (z, σ) stacked vertically,
and a transition matrix Q such that
Q[m, mp] = probability of transition x_states[m] -> x_states[mp]
The strategy is to
1. Discretize the σ process to produce state values σ_1, ..., σ_I
2. For each σ_i,
* discretize the z process to get z_{i1}, ... z_{iJ}
In each case, discretization uses Rouwenhorst's method
The final states are constructed as
x_m = (z_{ij}, σ_i), where m = (i - 1) * J + j.
Each x_m vector is stacked as a column of x_states. The transition
probability Q[m, n] from x_m to x_n is computed from the transition matrices
arising from the discretization of σ and z discussed above.
"""
function discretize_sv(sv::StochasticVolatility, I, J;
fail_with_neg_σ=false,
verbose=false)
# Unpack names
ρ, s_z, v, d, s_σ = sv.ρ, sv.s_z, sv.v, sv.d, sv.s_σ
# Discretize σ first
mc = rouwenhorst(I, v, s_σ, d)
sig_Q, sig2 = mc.p, collect(mc.state_values)
# This gives σ^2 values so now we take the square root
σ_states = similar(sig2)
if fail_with_neg_σ == true
@assert all(sig2 .>= 0) "Discretization failed: negative σ values."
else
for i in 1:I
σ_states[i] = sig2[i] < 0 ? 1e-8 : sqrt(sig2[i])
end
end
# Allocate memory
M = I * J
z_states = Array{Float64}(I, J)
q = Array{Float64}(I, J, J)
x_states = Array{Float64}(2, M)
Q = Array{Float64}(M, M)
# Discretize z at each σ_i and record state values for z in z_states.
# Also, record transition probability from z_states[i, j] to
# z_states[i, jp] when σ = σ_i. Store it as q[i, j, jp].
for (i, σ) in enumerate(σ_states)
mc_z = rouwenhorst(J, ρ, s_z * σ, 0.0)
for j in 1:J
z_states[i, j] = mc_z.state_values[j]
for jp in 1:J
q[i, j, jp] = mc_z.p[j, jp]
end
end
end
# Compute x_states and Q
for m in 1:M
i, j = multi_from_single(m, J)
x_states[:, m] = [z_states[i, j], σ_states[i]]
for mp in 1:M
ip, jp = multi_from_single(mp, J)
Q[m, mp] = sig_Q[i, ip] * q[i, j, jp]
end
end
if verbose == true
return x_states, Q, z_states, σ_states
else
return x_states, Q
end
end
Out[5]:
In [6]:
function StochasticVolatilityBY(; ρ=0.979,
s_z=0.044, # ϕ_x in PSW
v=0.987, # ν (nu) in PSW, vee here
σ_bar=0.0078, # same in PSW
s_σ=2.3e-6) # ϕ_σ in PSW
d = σ_bar^2 * (1 - v)
return StochasticVolatility(ρ, s_z, v, d, s_σ)
end
function StochasticVolatilityBKY(; ρ=0.975,
s_z=0.038, # ϕ_x in PSW
v=0.999, # ν (nu) in PSW, vee here
σ_bar=0.0072, # same in PSW
s_σ=2.8e-6) # ϕ_σ in PSW
d = σ_bar^2 * (1 - v)
return StochasticVolatility(ρ, s_z, v, d, s_σ)
end
Out[6]:
In [7]:
function compute_K_bansal_yaron(ez::EpsteinZin,
sv::StochasticVolatility;
μ=0.0015,
I=10, # discretization in σ
J=10) # discretization in z for each σ
# Unpack parameters, allocate memory
ψ, γ, β, ζ, θ = ez.ψ, ez.γ, ez.β, ez.ζ, ez.θ
M = I * J
K = Array{Float64}(M, M)
# Discretize SV process
x, Q = discretize_sv(sv, I, J)
for m in 1:M
for mp in 1:M
i, j = multi_from_single(m, J)
z, σ = x[1, m], x[2, m]
a = exp((1 - γ) * (μ + z) + (1 - γ)^2 * σ^2 / 2)
K[m, mp] = a * Q[m, mp]
end
end
return β^θ * K
end
Out[7]:
In [8]:
ez_by = EpsteinZinBY()
sv_by = StochasticVolatilityBY()
Out[8]:
In [9]:
K = compute_K_bansal_yaron(ez_by, sv_by);
In [10]:
compute_spec_rad(K)
Out[10]:
In [11]:
using PyPlot
plt = PyPlot
Out[11]:
In [12]:
J = 20 # grid size
R = Array{Float64}(J, J)
Out[12]:
In [13]:
x_vals = linspace(1.25, 2.25, J) # ψ
y_vals = linspace(0.0005, 0.01, J) # μ
Out[13]:
In [14]:
for (i, x) in enumerate(x_vals)
for (j, y) in enumerate(y_vals)
ez = EpsteinZinBY(ψ=x)
@assert ez.θ < 0 "Detected non-negative theta value"
K = compute_K_bansal_yaron(ez, sv_by, μ=y)
R[i, j] = compute_spec_rad(K)
end
end
In [15]:
fig, ax = plt.subplots(figsize=(10, 5.7))
#lvs = [0.0, 0.8, 1.0, 1.4, 1.8, 2.2, 4.4]
#cls = [cm.jet(i) for i in np.linspace(0.4, 1, len(lvs))]
cs1 = ax[:contourf](x_vals,
y_vals,
R',
alpha=0.6)
#levels=lvs,
ctr1 = ax[:contour](x_vals,
y_vals,
R',
levels=[1.0])
plt.clabel(ctr1, inline=1, fontsize=13)
plt.colorbar(cs1, ax=ax)
ax[:set_title]("Spectral radius")
ax[:set_xlabel]("ψ", fontsize=16)
ax[:set_ylabel]("μ", fontsize=16)
ax[:annotate]("Bansal-Yaron",
xy=(1.5 + 0.024, 0.0015 - 0.0001),
xycoords="data",
xytext=(5, 40),
textcoords="offset points",
fontsize=12,
arrowprops=Dict("arrowstyle" => "->"))
ax[:plot]([1.5], [0.0015], "ko", alpha=0.6)
plt.savefig("by.pdf")
plt.show()