HumphreyYang
commented
3 months ago Since Tom is actively working on this lecture, I suggest we use sequence unpacking to resolve this issue:
class ChangLQ:
"""
Class to solve LQ Chang model
"""
def __init__(self, α, α0, α1, α2, c, T=1000, θ_n=200):
# Record parameters
self.α, self.α0, self.α1 = α, α0, α1
self.α2, self.c, self.T, self.θ_n = α2, c, T, θ_n
# Create β using "Poor Man's Friedman Rule"
self.β = np.exp(-α1 / (α * α2))
# Solve the Ramsey Problem #
# LQ Matrices
R = -np.array([[α0, -α1 * α / 2],
[-α1 * α/2, -α2 * α**2 / 2]])
Q = -np.array([[-c / 2]])
A = np.array([[1, 0], [0, (1 + α) / α]])
B = np.array([[0], [-1 / α]])
# Solve LQ Problem (Subproblem 1)
lq = LQ(Q, R, A, B, beta=self.β)
self.P, self.F, self.d = lq.stationary_values()
# Solve Subproblem 2
self.θ_R = -self.P[0, 1] / self.P[1, 1]
# Find bliss level of θ
self.θ_B = np.log(self.β)
# Solve the Markov Perfect Equilibrium
self.μ_MPE = -α1 / ((1 + α) / α * c + α / (1 + α)
* α2 + α**2 / (1 + α) * α2)
self.θ_MPE = self.μ_MPE
self.μ_check = -α * α1 / (α2 * α**2 + c)
# Calculate value under MPE and Check economy
self.J_MPE = (α0 + α1 * (-α * self.μ_MPE) - α2 / 2
* (-α * self.μ_MPE)**2 - c/2 * self.μ_MPE**2) / (1 - self.β)
self.J_check = (α0 + α1 * (-α * self.μ_check) - α2/2
* (-α * self.μ_check)**2 - c / 2 * self.μ_check**2) \
/ (1 - self.β)
# Simulate Ramsey plan for large number of periods
θ_series = np.vstack((np.ones((1, T)), np.zeros((1, T))))
μ_series = np.zeros(T)
J_series = np.zeros(T)
θ_series[1, 0] = self.θ_R
[μ_series[0]] = -self.F.dot(θ_series[:, 0])
J_series[0] = -θ_series[:, 0] @ self.P @ θ_series[:, 0].T
for i in range(1, T):
θ_series[:, i] = (A - B @ self.F) @ θ_series[:, i-1]
[μ_series[i]] = -self.F @ θ_series[:, i]
J_series[i] = -θ_series[:, i] @ self.P @ θ_series[:, i].T
self.J_series = J_series
self.μ_series = μ_series
self.θ_series = θ_series
# Find the range of θ in Ramsey plan
θ_LB = min(θ_series[1, :])
θ_LB = min(θ_LB, self.θ_B)
θ_UB = max(θ_series[1, :])
θ_UB = max(θ_UB, self.θ_MPE)
θ_range = θ_UB - θ_LB
self.θ_LB = θ_LB - 0.05 * θ_range
self.θ_UB = θ_UB + 0.05 * θ_range
self.θ_range = θ_range
# Find value function and policy functions over range of θ
θ_space = np.linspace(self.θ_LB, self.θ_UB, 200)
J_space = np.zeros(200)
check_space = np.zeros(200)
μ_space = np.zeros(200)
θ_prime = np.zeros(200)
for i in range(200):
J_space[i] = - np.array((1, θ_space[i])) \
@ self.P @ np.array((1, θ_space[i])).T
[μ_space[i]] = - self.F @ np.array((1, θ_space[i]))
x_prime = (A - B @ self.F) @ np.array((1, θ_space[i]))
θ_prime[i] = x_prime[1]
check_space[i] = (α0 + α1 * (-α * θ_space[i]) -
α2/2 * (-α * θ_space[i])**2 - c/2 * θ_space[i]**2) / (1 - self.β)
J_LB = min(J_space)
J_UB = max(J_space)
J_range = J_UB - J_LB
self.J_LB = J_LB - 0.05 * J_range
self.J_UB = J_UB + 0.05 * J_range
self.J_range = J_range
self.J_space = J_space
self.θ_space = θ_space
self.μ_space = μ_space
self.θ_prime = θ_prime
self.check_space = check_space(i.e.
[μ_space[i]] = - self.F @ np.array((1, θ_space[i])))
to fix cases where single-element list is nested in the array.
DeprecationWarning is raised by NumPy 1.25:
in
calvo.md.Related issue: https://github.com/numpy/numpy/pull/10615 https://numpy.org/doc/stable/release/1.25.0-notes.html