diff --git a/source/_static/lecture_specific/amss/recursive_allocation.py b/source/_static/lecture_specific/amss/recursive_allocation.py index 43fde64..f5495f1 100644 --- a/source/_static/lecture_specific/amss/recursive_allocation.py +++ b/source/_static/lecture_specific/amss/recursive_allocation.py @@ -1,298 +1,212 @@ -import numpy as np -from scipy.optimize import fmin_slsqp -from scipy.optimize import root -from quantecon import MarkovChain - - -class RecursiveAllocationAMSS: - - def __init__(self, model, μgrid, tol_diff=1e-7, tol=1e-7): - - self.β, self.π, self.G = model.β, model.π, model.G - self.mc, self.S = MarkovChain(self.π), len(model.π) # Number of states - self.Θ, self.model, self.μgrid = model.Θ, model, μgrid - self.tol_diff, self.tol = tol_diff, tol - - # Find the first best allocation - self.solve_time1_bellman() - self.T.time_0 = True # Bellman equation now solves time 0 problem - - def solve_time1_bellman(self): - ''' - Solve the time 1 Bellman equation for calibration model and - initial grid μgrid0 - ''' - model, μgrid0 = self.model, self.μgrid - π = model.π - S = len(model.π) - - # First get initial fit from Lucas Stokey solution. - # Need to change things to be ex ante - pp = SequentialAllocation(model) - interp = interpolator_factory(2, None) - - def incomplete_allocation(μ_, s_): - c, n, x, V = pp.time1_value(μ_) - return c, n, π[s_] @ x, π[s_] @ V - cf, nf, xgrid, Vf, xprimef = [], [], [], [], [] - for s_ in range(S): - c, n, x, V = zip(*map(lambda μ: incomplete_allocation(μ, s_), μgrid0)) - c, n = np.vstack(c).T, np.vstack(n).T - x, V = np.hstack(x), np.hstack(V) - xprimes = np.vstack([x] * S) - cf.append(interp(x, c)) - nf.append(interp(x, n)) - Vf.append(interp(x, V)) - xgrid.append(x) - xprimef.append(interp(x, xprimes)) - cf, nf, xprimef = fun_vstack(cf), fun_vstack(nf), fun_vstack(xprimef) - Vf = fun_hstack(Vf) - policies = [cf, nf, xprimef] - - # Create xgrid - x = np.vstack(xgrid).T - xbar = [x.min(0).max(), x.max(0).min()] - xgrid = np.linspace(xbar[0], xbar[1], len(μgrid0)) - self.xgrid = xgrid - - # Now iterate on Bellman equation - T = BellmanEquation(model, xgrid, policies, tol=self.tol) - diff = 1 - while diff > self.tol_diff: - PF = T(Vf) - - Vfnew, policies = self.fit_policy_function(PF) - diff = np.abs((Vf(xgrid) - Vfnew(xgrid)) / Vf(xgrid)).max() - - print(diff) - Vf = Vfnew - - # Store value function policies and Bellman Equations - self.Vf = Vf - self.policies = policies - self.T = T - - def fit_policy_function(self, PF): - ''' - Fits the policy functions - ''' - S, xgrid = len(self.π), self.xgrid - interp = interpolator_factory(3, 0) - cf, nf, xprimef, Tf, Vf = [], [], [], [], [] - for s_ in range(S): - PFvec = np.vstack([PF(x, s_) for x in self.xgrid]).T - Vf.append(interp(xgrid, PFvec[0, :])) - cf.append(interp(xgrid, PFvec[1:1 + S])) - nf.append(interp(xgrid, PFvec[1 + S:1 + 2 * S])) - xprimef.append(interp(xgrid, PFvec[1 + 2 * S:1 + 3 * S])) - Tf.append(interp(xgrid, PFvec[1 + 3 * S:])) - policies = fun_vstack(cf), fun_vstack( - nf), fun_vstack(xprimef), fun_vstack(Tf) - Vf = fun_hstack(Vf) - return Vf, policies - - def Τ(self, c, n): - ''' - Computes Τ given c and n - ''' - model = self.model - Uc, Un = model.Uc(c, n), model.Un(c, n) - - return 1 + Un / (self.Θ * Uc) - - def time0_allocation(self, B_, s0): - ''' - Finds the optimal allocation given initial government debt B_ and - state s_0 - ''' - PF = self.T(self.Vf) - z0 = PF(B_, s0) - c0, n0, xprime0, T0 = z0[1:] - return c0, n0, xprime0, T0 - - def simulate(self, B_, s_0, T, sHist=None): - ''' - Simulates planners policies for T periods - ''' - model, π = self.model, self.π - Uc = model.Uc - cf, nf, xprimef, Tf = self.policies - - if sHist is None: - sHist = simulate_markov(π, s_0, T) - - cHist, nHist, Bhist, xHist, ΤHist, THist, μHist = np.zeros((7, T)) - # Time 0 - cHist[0], nHist[0], xHist[0], THist[0] = self.time0_allocation(B_, s_0) - ΤHist[0] = self.Τ(cHist[0], nHist[0])[s_0] - Bhist[0] = B_ - μHist[0] = self.Vf[s_0](xHist[0]) - - # Time 1 onward - for t in range(1, T): - s_, x, s = sHist[t - 1], xHist[t - 1], sHist[t] - c, n, xprime, T = cf[s_, :](x), nf[s_, :]( - x), xprimef[s_, :](x), Tf[s_, :](x) - - Τ = self.Τ(c, n)[s] - u_c = Uc(c, n) - Eu_c = π[s_, :] @ u_c - - μHist[t] = self.Vf[s](xprime[s]) - - cHist[t], nHist[t], Bhist[t], ΤHist[t] = c[s], n[s], x / Eu_c, Τ - xHist[t], THist[t] = xprime[s], T[s] - return np.array([cHist, nHist, Bhist, ΤHist, THist, μHist, sHist, xHist]) - - -class BellmanEquation: - ''' - Bellman equation for the continuation of the Lucas-Stokey Problem - ''' - - def __init__(self, model, xgrid, policies0, tol, maxiter=1000): - - self.β, self.π, self.G = model.β, model.π, model.G - self.S = len(model.π) # Number of states - self.Θ, self.model, self.tol = model.Θ, model, tol - self.maxiter = maxiter - - self.xbar = [min(xgrid), max(xgrid)] - self.time_0 = False - - self.z0 = {} - cf, nf, xprimef = policies0 - - for s_ in range(self.S): - for x in xgrid: - self.z0[x, s_] = np.hstack([cf[s_, :](x), - nf[s_, :](x), - xprimef[s_, :](x), - np.zeros(self.S)]) - - self.find_first_best() - - def find_first_best(self): - ''' - Find the first best allocation - ''' - model = self.model - S, Θ, Uc, Un, G = self.S, self.Θ, model.Uc, model.Un, self.G - - def res(z): - c = z[:S] - n = z[S:] - return np.hstack([Θ * Uc(c, n) + Un(c, n), Θ * n - c - G]) - - res = root(res, 0.5 * np.ones(2 * S)) - if not res.success: - raise Exception('Could not find first best') - - self.cFB = res.x[:S] - self.nFB = res.x[S:] - IFB = Uc(self.cFB, self.nFB) * self.cFB + \ - Un(self.cFB, self.nFB) * self.nFB - - self.xFB = np.linalg.solve(np.eye(S) - self.β * self.π, IFB) - - self.zFB = {} - for s in range(S): - self.zFB[s] = np.hstack( - [self.cFB[s], self.nFB[s], self.π[s] @ self.xFB, 0.]) - - def __call__(self, Vf): - ''' - Given continuation value function next period return value function this - period return T(V) and optimal policies - ''' - if not self.time_0: - def PF(x, s): return self.get_policies_time1(x, s, Vf) - else: - def PF(B_, s0): return self.get_policies_time0(B_, s0, Vf) - return PF - - def get_policies_time1(self, x, s_, Vf): - ''' - Finds the optimal policies - ''' - model, β, Θ, G, S, π = self.model, self.β, self.Θ, self.G, self.S, self.π - U, Uc, Un = model.U, model.Uc, model.Un - - def objf(z): - c, n, xprime = z[:S], z[S:2 * S], z[2 * S:3 * S] - - Vprime = np.empty(S) +class AMSS: + # WARNING: THE CODE IS EXTREMELY SENSITIVE TO CHOCIES OF PARAMETERS. + # DO NOT CHANGE THE PARAMETERS AND EXPECT IT TO WORK + + def __init__(self, pref, β, Π, g, x_grid, bounds_v): + self.β, self.Π, self.g = β, Π, g + self.x_grid = x_grid + self.n = x_grid[0][2] + self.S = len(Π) + self.bounds = bounds_v + self.pref = pref + + self.T_v, self.T_w = bellman_operator_factory(Π, β, x_grid, g, + bounds_v) + + self.V_solved = False + self.W_solved = False + + def compute_V(self, V, σ_v_star, tol_vfi, maxitr, print_itr): + + T_v = self.T_v + + self.success = False + + V_new = np.zeros_like(V) + + Δ = 1.0 + for itr in range(maxitr): + T_v(V, V_new, σ_v_star, self.pref) + + Δ = np.max(np.abs(V_new - V)) + + if Δ < tol_vfi: + self.V_solved = True + print('Successfully completed VFI after %i iterations' + % (itr+1)) + break + + if (itr + 1) % print_itr == 0: + print('Error at iteration %i : ' % (itr + 1), Δ) + + V[:] = V_new[:] + + self.V = V + self.σ_v_star = σ_v_star + + return V, σ_v_star + + def compute_W(self, b_0, W, σ_w_star): + T_w = self.T_w + V = self.V + + T_w(W, σ_w_star, V, b_0, self.pref) + + self.W = W + self.σ_w_star = σ_w_star + self.W_solved = True + print('Succesfully solved the time 0 problem.') + + return W, σ_w_star + + def solve(self, V, σ_v_star, b_0, W, σ_w_star, tol_vfi=1e-7, + maxitr=1000, print_itr=10): + print("===============") + print("Solve time 1 problem") + print("===============") + self.compute_V(V, σ_v_star, tol_vfi, maxitr, print_itr) + print("===============") + print("Solve time 0 problem") + print("===============") + self.compute_W(b_0, W, σ_w_star) + + def simulate(self, s_hist, b_0): + if not (self.V_solved and self.W_solved): + msg = "V and W need to be successfully computed before simulation." + raise ValueError(msg) + + pref = self.pref + x_grid, g, β, S = self.x_grid, self.g, self.β, self.S + σ_v_star, σ_w_star = self.σ_v_star, self.σ_w_star + + T = len(s_hist) + s_0 = s_hist[0] + + # Pre-allocate + n_hist = np.zeros(T) + x_hist = np.zeros(T) + c_hist = np.zeros(T) + τ_hist = np.zeros(T) + b_hist = np.zeros(T) + g_hist = np.zeros(T) + + # Compute t = 0 + l_0, T_0 = σ_w_star[s_0] + c_0 = (1 - l_0) - g[s_0] + x_0 = (-pref.Uc(c_0, l_0) * (c_0 - T_0 - b_0) + + pref.Ul(c_0, l_0) * (1 - l_0)) + + n_hist[0] = (1 - l_0) + x_hist[0] = x_0 + c_hist[0] = c_0 + τ_hist[0] = 1 - pref.Ul(c_0, l_0) / pref.Uc(c_0, l_0) + b_hist[0] = b_0 + g_hist[0] = g[s_0] + + # Compute t > 0 + for t in range(T - 1): + x_ = x_hist[t] + s_ = s_hist[t] + l = np.zeros(S) + T = np.zeros(S) for s in range(S): - Vprime[s] = Vf[s](xprime[s]) - - return -π[s_] @ (U(c, n) + β * Vprime) - - def objf_prime(x): - - epsilon = 1e-7 - x0 = np.asfarray(x) - f0 = np.atleast_1d(objf(x0)) - jac = np.zeros([len(x0), len(f0)]) - dx = np.zeros(len(x0)) - for i in range(len(x0)): - dx[i] = epsilon - jac[i] = (objf(x0+dx) - f0)/epsilon - dx[i] = 0.0 - - return jac.transpose() - - def cons(z): - c, n, xprime, T = z[:S], z[S:2 * S], z[2 * S:3 * S], z[3 * S:] - u_c = Uc(c, n) - Eu_c = π[s_] @ u_c - return np.hstack([ - x * u_c / Eu_c - u_c * (c - T) - Un(c, n) * n - β * xprime, - Θ * n - c - G]) - - if model.transfers: - bounds = [(0., 100)] * S + [(0., 100)] * S + \ - [self.xbar] * S + [(0., 100.)] * S - else: - bounds = [(0., 100)] * S + [(0., 100)] * S + \ - [self.xbar] * S + [(0., 0.)] * S - out, fx, _, imode, smode = fmin_slsqp(objf, self.z0[x, s_], - f_eqcons=cons, bounds=bounds, - fprime=objf_prime, full_output=True, - iprint=0, acc=self.tol, iter=self.maxiter) - - if imode > 0: - raise Exception(smode) - - self.z0[x, s_] = out - return np.hstack([-fx, out]) - - def get_policies_time0(self, B_, s0, Vf): - ''' - Finds the optimal policies - ''' - model, β, Θ, G = self.model, self.β, self.Θ, self.G - U, Uc, Un = model.U, model.Uc, model.Un - - def objf(z): - c, n, xprime = z[:-1] - - return -(U(c, n) + β * Vf[s0](xprime)) - - def cons(z): - c, n, xprime, T = z - return np.hstack([ - -Uc(c, n) * (c - B_ - T) - Un(c, n) * n - β * xprime, - (Θ * n - c - G)[s0]]) - - if model.transfers: - bounds = [(0., 100), (0., 100), self.xbar, (0., 100.)] - else: - bounds = [(0., 100), (0., 100), self.xbar, (0., 0.)] - out, fx, _, imode, smode = fmin_slsqp(objf, self.zFB[s0], f_eqcons=cons, - bounds=bounds, full_output=True, - iprint=0) - - if imode > 0: - raise Exception(smode) - - return np.hstack([-fx, out]) + x_arr = np.array([x_]) + l[s] = eval_linear(x_grid, σ_v_star[s_, :, s], x_arr) + T[s] = eval_linear(x_grid, σ_v_star[s_, :, S+s], x_arr) + + c = (1 - l) - g + u_c = pref.Uc(c, l) + Eu_c = Π[s_] @ u_c + + x = u_c * x_ / (β * Eu_c) - u_c * (c - T) + pref.Ul(c, l) * (1 - l) + + c_next = c[s_hist[t+1]] + l_next = l[s_hist[t+1]] + + x_hist[t+1] = x[s_hist[t+1]] + n_hist[t+1] = 1 - l_next + c_hist[t+1] = c_next + τ_hist[t+1] = 1 - pref.Ul(c_next, l_next) / pref.Uc(c_next, l_next) + b_hist[t+1] = x_ / (β * Eu_c) + g_hist[t+1] = g[s_hist[t+1]] + + return c_hist, n_hist, b_hist, τ_hist, g_hist, n_hist + + +def obj_factory(Π, β, x_grid, g): + S = len(Π) + + @njit + def obj_V(σ, state, V, pref): + # Unpack state + s_, x_ = state + + l = σ[:S] + T = σ[S:] + + c = (1 - l) - g + u_c = pref.Uc(c, l) + Eu_c = Π[s_] @ u_c + x = u_c * x_ / (β * Eu_c) - u_c * (c - T) + pref.Ul(c, l) * (1 - l) + + V_next = np.zeros(S) + + for s in range(S): + V_next[s] = eval_linear(x_grid, V[s], np.array([x[s]])) + + out = Π[s_] @ (pref.U(c, l) + β * V_next) + + return out + + @njit + def obj_W(σ, state, V, pref): + # Unpack state + s_, b_0 = state + l, T = σ + + c = (1 - l) - g[s_] + x = -pref.Uc(c, l) * (c - T - b_0) + pref.Ul(c, l) * (1 - l) + + V_next = eval_linear(x_grid, V[s_], np.array([x])) + + out = pref.U(c, l) + β * V_next + + return out + + return obj_V, obj_W + + +def bellman_operator_factory(Π, β, x_grid, g, bounds_v): + obj_V, obj_W = obj_factory(Π, β, x_grid, g) + n = x_grid[0][2] + S = len(Π) + x_nodes = nodes(x_grid) + + @njit(parallel=True) + def T_v(V, V_new, σ_star, pref): + for s_ in prange(S): + for x_i in prange(n): + state = (s_, x_nodes[x_i]) + x0 = σ_star[s_, x_i] + res = optimize.nelder_mead(obj_V, x0, bounds=bounds_v, + args=(state, V, pref)) + + if res.success: + V_new[s_, x_i] = res.fun + σ_star[s_, x_i] = res.x + else: + print("Optimization routine failed.") + + bounds_w = np.array([[-9.0, 1.0], [0., 10.]]) + + def T_w(W, σ_star, V, b_0, pref): + for s_ in prange(S): + state = (s_, b_0) + x0 = σ_star[s_] + res = optimize.nelder_mead(obj_W, x0, bounds=bounds_w, + args=(state, V, pref)) + + W[s_] = res.fun + σ_star[s_] = res.x + + return T_v, T_w diff --git a/source/_static/lecture_specific/amss2/crra_utility.py b/source/_static/lecture_specific/amss2/crra_utility.py new file mode 100644 index 0000000..700bbc9 --- /dev/null +++ b/source/_static/lecture_specific/amss2/crra_utility.py @@ -0,0 +1,38 @@ +import numpy as np + + +class CRRAutility: + + def __init__(self, + β=0.9, + σ=2, + γ=2, + π=0.5*np.ones((2, 2)), + G=np.array([0.1, 0.2]), + Θ=np.ones(2), + transfers=False): + + self.β, self.σ, self.γ = β, σ, γ + self.π, self.G, self.Θ, self.transfers = π, G, Θ, transfers + + # Utility function + def U(self, c, n): + σ = self.σ + if σ == 1.: + U = np.log(c) + else: + U = (c**(1 - σ) - 1) / (1 - σ) + return U - n**(1 + self.γ) / (1 + self.γ) + + # Derivatives of utility function + def Uc(self, c, n): + return c**(-self.σ) + + def Ucc(self, c, n): + return -self.σ * c**(-self.σ - 1) + + def Un(self, c, n): + return -n**self.γ + + def Unn(self, c, n): + return -self.γ * n**(self.γ - 1) \ No newline at end of file diff --git a/source/_static/lecture_specific/amss2/log_utility.py b/source/_static/lecture_specific/amss2/log_utility.py new file mode 100644 index 0000000..decb608 --- /dev/null +++ b/source/_static/lecture_specific/amss2/log_utility.py @@ -0,0 +1,31 @@ +import numpy as np + +class LogUtility: + + def __init__(self, + β=0.9, + ψ=0.69, + π=0.5*np.ones((2, 2)), + G=np.array([0.1, 0.2]), + Θ=np.ones(2), + transfers=False): + + self.β, self.ψ, self.π = β, ψ, π + self.G, self.Θ, self.transfers = G, Θ, transfers + + # Utility function + def U(self, c, n): + return np.log(c) + self.ψ * np.log(1 - n) + + # Derivatives of utility function + def Uc(self, c, n): + return 1 / c + + def Ucc(self, c, n): + return -c**(-2) + + def Un(self, c, n): + return -self.ψ / (1 - n) + + def Unn(self, c, n): + return -self.ψ / (1 - n)**2 \ No newline at end of file diff --git a/source/_static/lecture_specific/amss2/recursive_allocation.py b/source/_static/lecture_specific/amss2/recursive_allocation.py new file mode 100644 index 0000000..43fde64 --- /dev/null +++ b/source/_static/lecture_specific/amss2/recursive_allocation.py @@ -0,0 +1,298 @@ +import numpy as np +from scipy.optimize import fmin_slsqp +from scipy.optimize import root +from quantecon import MarkovChain + + +class RecursiveAllocationAMSS: + + def __init__(self, model, μgrid, tol_diff=1e-7, tol=1e-7): + + self.β, self.π, self.G = model.β, model.π, model.G + self.mc, self.S = MarkovChain(self.π), len(model.π) # Number of states + self.Θ, self.model, self.μgrid = model.Θ, model, μgrid + self.tol_diff, self.tol = tol_diff, tol + + # Find the first best allocation + self.solve_time1_bellman() + self.T.time_0 = True # Bellman equation now solves time 0 problem + + def solve_time1_bellman(self): + ''' + Solve the time 1 Bellman equation for calibration model and + initial grid μgrid0 + ''' + model, μgrid0 = self.model, self.μgrid + π = model.π + S = len(model.π) + + # First get initial fit from Lucas Stokey solution. + # Need to change things to be ex ante + pp = SequentialAllocation(model) + interp = interpolator_factory(2, None) + + def incomplete_allocation(μ_, s_): + c, n, x, V = pp.time1_value(μ_) + return c, n, π[s_] @ x, π[s_] @ V + cf, nf, xgrid, Vf, xprimef = [], [], [], [], [] + for s_ in range(S): + c, n, x, V = zip(*map(lambda μ: incomplete_allocation(μ, s_), μgrid0)) + c, n = np.vstack(c).T, np.vstack(n).T + x, V = np.hstack(x), np.hstack(V) + xprimes = np.vstack([x] * S) + cf.append(interp(x, c)) + nf.append(interp(x, n)) + Vf.append(interp(x, V)) + xgrid.append(x) + xprimef.append(interp(x, xprimes)) + cf, nf, xprimef = fun_vstack(cf), fun_vstack(nf), fun_vstack(xprimef) + Vf = fun_hstack(Vf) + policies = [cf, nf, xprimef] + + # Create xgrid + x = np.vstack(xgrid).T + xbar = [x.min(0).max(), x.max(0).min()] + xgrid = np.linspace(xbar[0], xbar[1], len(μgrid0)) + self.xgrid = xgrid + + # Now iterate on Bellman equation + T = BellmanEquation(model, xgrid, policies, tol=self.tol) + diff = 1 + while diff > self.tol_diff: + PF = T(Vf) + + Vfnew, policies = self.fit_policy_function(PF) + diff = np.abs((Vf(xgrid) - Vfnew(xgrid)) / Vf(xgrid)).max() + + print(diff) + Vf = Vfnew + + # Store value function policies and Bellman Equations + self.Vf = Vf + self.policies = policies + self.T = T + + def fit_policy_function(self, PF): + ''' + Fits the policy functions + ''' + S, xgrid = len(self.π), self.xgrid + interp = interpolator_factory(3, 0) + cf, nf, xprimef, Tf, Vf = [], [], [], [], [] + for s_ in range(S): + PFvec = np.vstack([PF(x, s_) for x in self.xgrid]).T + Vf.append(interp(xgrid, PFvec[0, :])) + cf.append(interp(xgrid, PFvec[1:1 + S])) + nf.append(interp(xgrid, PFvec[1 + S:1 + 2 * S])) + xprimef.append(interp(xgrid, PFvec[1 + 2 * S:1 + 3 * S])) + Tf.append(interp(xgrid, PFvec[1 + 3 * S:])) + policies = fun_vstack(cf), fun_vstack( + nf), fun_vstack(xprimef), fun_vstack(Tf) + Vf = fun_hstack(Vf) + return Vf, policies + + def Τ(self, c, n): + ''' + Computes Τ given c and n + ''' + model = self.model + Uc, Un = model.Uc(c, n), model.Un(c, n) + + return 1 + Un / (self.Θ * Uc) + + def time0_allocation(self, B_, s0): + ''' + Finds the optimal allocation given initial government debt B_ and + state s_0 + ''' + PF = self.T(self.Vf) + z0 = PF(B_, s0) + c0, n0, xprime0, T0 = z0[1:] + return c0, n0, xprime0, T0 + + def simulate(self, B_, s_0, T, sHist=None): + ''' + Simulates planners policies for T periods + ''' + model, π = self.model, self.π + Uc = model.Uc + cf, nf, xprimef, Tf = self.policies + + if sHist is None: + sHist = simulate_markov(π, s_0, T) + + cHist, nHist, Bhist, xHist, ΤHist, THist, μHist = np.zeros((7, T)) + # Time 0 + cHist[0], nHist[0], xHist[0], THist[0] = self.time0_allocation(B_, s_0) + ΤHist[0] = self.Τ(cHist[0], nHist[0])[s_0] + Bhist[0] = B_ + μHist[0] = self.Vf[s_0](xHist[0]) + + # Time 1 onward + for t in range(1, T): + s_, x, s = sHist[t - 1], xHist[t - 1], sHist[t] + c, n, xprime, T = cf[s_, :](x), nf[s_, :]( + x), xprimef[s_, :](x), Tf[s_, :](x) + + Τ = self.Τ(c, n)[s] + u_c = Uc(c, n) + Eu_c = π[s_, :] @ u_c + + μHist[t] = self.Vf[s](xprime[s]) + + cHist[t], nHist[t], Bhist[t], ΤHist[t] = c[s], n[s], x / Eu_c, Τ + xHist[t], THist[t] = xprime[s], T[s] + return np.array([cHist, nHist, Bhist, ΤHist, THist, μHist, sHist, xHist]) + + +class BellmanEquation: + ''' + Bellman equation for the continuation of the Lucas-Stokey Problem + ''' + + def __init__(self, model, xgrid, policies0, tol, maxiter=1000): + + self.β, self.π, self.G = model.β, model.π, model.G + self.S = len(model.π) # Number of states + self.Θ, self.model, self.tol = model.Θ, model, tol + self.maxiter = maxiter + + self.xbar = [min(xgrid), max(xgrid)] + self.time_0 = False + + self.z0 = {} + cf, nf, xprimef = policies0 + + for s_ in range(self.S): + for x in xgrid: + self.z0[x, s_] = np.hstack([cf[s_, :](x), + nf[s_, :](x), + xprimef[s_, :](x), + np.zeros(self.S)]) + + self.find_first_best() + + def find_first_best(self): + ''' + Find the first best allocation + ''' + model = self.model + S, Θ, Uc, Un, G = self.S, self.Θ, model.Uc, model.Un, self.G + + def res(z): + c = z[:S] + n = z[S:] + return np.hstack([Θ * Uc(c, n) + Un(c, n), Θ * n - c - G]) + + res = root(res, 0.5 * np.ones(2 * S)) + if not res.success: + raise Exception('Could not find first best') + + self.cFB = res.x[:S] + self.nFB = res.x[S:] + IFB = Uc(self.cFB, self.nFB) * self.cFB + \ + Un(self.cFB, self.nFB) * self.nFB + + self.xFB = np.linalg.solve(np.eye(S) - self.β * self.π, IFB) + + self.zFB = {} + for s in range(S): + self.zFB[s] = np.hstack( + [self.cFB[s], self.nFB[s], self.π[s] @ self.xFB, 0.]) + + def __call__(self, Vf): + ''' + Given continuation value function next period return value function this + period return T(V) and optimal policies + ''' + if not self.time_0: + def PF(x, s): return self.get_policies_time1(x, s, Vf) + else: + def PF(B_, s0): return self.get_policies_time0(B_, s0, Vf) + return PF + + def get_policies_time1(self, x, s_, Vf): + ''' + Finds the optimal policies + ''' + model, β, Θ, G, S, π = self.model, self.β, self.Θ, self.G, self.S, self.π + U, Uc, Un = model.U, model.Uc, model.Un + + def objf(z): + c, n, xprime = z[:S], z[S:2 * S], z[2 * S:3 * S] + + Vprime = np.empty(S) + for s in range(S): + Vprime[s] = Vf[s](xprime[s]) + + return -π[s_] @ (U(c, n) + β * Vprime) + + def objf_prime(x): + + epsilon = 1e-7 + x0 = np.asfarray(x) + f0 = np.atleast_1d(objf(x0)) + jac = np.zeros([len(x0), len(f0)]) + dx = np.zeros(len(x0)) + for i in range(len(x0)): + dx[i] = epsilon + jac[i] = (objf(x0+dx) - f0)/epsilon + dx[i] = 0.0 + + return jac.transpose() + + def cons(z): + c, n, xprime, T = z[:S], z[S:2 * S], z[2 * S:3 * S], z[3 * S:] + u_c = Uc(c, n) + Eu_c = π[s_] @ u_c + return np.hstack([ + x * u_c / Eu_c - u_c * (c - T) - Un(c, n) * n - β * xprime, + Θ * n - c - G]) + + if model.transfers: + bounds = [(0., 100)] * S + [(0., 100)] * S + \ + [self.xbar] * S + [(0., 100.)] * S + else: + bounds = [(0., 100)] * S + [(0., 100)] * S + \ + [self.xbar] * S + [(0., 0.)] * S + out, fx, _, imode, smode = fmin_slsqp(objf, self.z0[x, s_], + f_eqcons=cons, bounds=bounds, + fprime=objf_prime, full_output=True, + iprint=0, acc=self.tol, iter=self.maxiter) + + if imode > 0: + raise Exception(smode) + + self.z0[x, s_] = out + return np.hstack([-fx, out]) + + def get_policies_time0(self, B_, s0, Vf): + ''' + Finds the optimal policies + ''' + model, β, Θ, G = self.model, self.β, self.Θ, self.G + U, Uc, Un = model.U, model.Uc, model.Un + + def objf(z): + c, n, xprime = z[:-1] + + return -(U(c, n) + β * Vf[s0](xprime)) + + def cons(z): + c, n, xprime, T = z + return np.hstack([ + -Uc(c, n) * (c - B_ - T) - Un(c, n) * n - β * xprime, + (Θ * n - c - G)[s0]]) + + if model.transfers: + bounds = [(0., 100), (0., 100), self.xbar, (0., 100.)] + else: + bounds = [(0., 100), (0., 100), self.xbar, (0., 0.)] + out, fx, _, imode, smode = fmin_slsqp(objf, self.zFB[s0], f_eqcons=cons, + bounds=bounds, full_output=True, + iprint=0) + + if imode > 0: + raise Exception(smode) + + return np.hstack([-fx, out]) diff --git a/source/_static/lecture_specific/amss2/sequential_allocation.py b/source/_static/lecture_specific/amss2/sequential_allocation.py new file mode 100644 index 0000000..f8d5744 --- /dev/null +++ b/source/_static/lecture_specific/amss2/sequential_allocation.py @@ -0,0 +1,158 @@ +import numpy as np +from scipy.optimize import root +from quantecon import MarkovChain + + +class SequentialAllocation: + + ''' + Class that takes CESutility or BGPutility object as input returns + planner's allocation as a function of the multiplier on the + implementability constraint μ. + ''' + + def __init__(self, model): + + # Initialize from model object attributes + self.β, self.π, self.G = model.β, model.π, model.G + self.mc, self.Θ = MarkovChain(self.π), model.Θ + self.S = len(model.π) # Number of states + self.model = model + + # Find the first best allocation + self.find_first_best() + + def find_first_best(self): + ''' + Find the first best allocation + ''' + model = self.model + S, Θ, G = self.S, self.Θ, self.G + Uc, Un = model.Uc, model.Un + + def res(z): + c = z[:S] + n = z[S:] + return np.hstack([Θ * Uc(c, n) + Un(c, n), Θ * n - c - G]) + + res = root(res, 0.5 * np.ones(2 * S)) + + if not res.success: + raise Exception('Could not find first best') + + self.cFB = res.x[:S] + self.nFB = res.x[S:] + + # Multiplier on the resource constraint + self.ΞFB = Uc(self.cFB, self.nFB) + self.zFB = np.hstack([self.cFB, self.nFB, self.ΞFB]) + + def time1_allocation(self, μ): + ''' + Computes optimal allocation for time t >= 1 for a given μ + ''' + model = self.model + S, Θ, G = self.S, self.Θ, self.G + Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn + + def FOC(z): + c = z[:S] + n = z[S:2 * S] + Ξ = z[2 * S:] + # FOC of c + return np.hstack([Uc(c, n) - μ * (Ucc(c, n) * c + Uc(c, n)) - Ξ, + Un(c, n) - μ * (Unn(c, n) * n + Un(c, n)) \ + + Θ * Ξ, # FOC of n + Θ * n - c - G]) + + # Find the root of the first-order condition + res = root(FOC, self.zFB) + if not res.success: + raise Exception('Could not find LS allocation.') + z = res.x + c, n, Ξ = z[:S], z[S:2 * S], z[2 * S:] + + # Compute x + I = Uc(c, n) * c + Un(c, n) * n + x = np.linalg.solve(np.eye(S) - self.β * self.π, I) + + return c, n, x, Ξ + + def time0_allocation(self, B_, s_0): + ''' + Finds the optimal allocation given initial government debt B_ and + state s_0 + ''' + model, π, Θ, G, β = self.model, self.π, self.Θ, self.G, self.β + Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn + + # First order conditions of planner's problem + def FOC(z): + μ, c, n, Ξ = z + xprime = self.time1_allocation(μ)[2] + return np.hstack([Uc(c, n) * (c - B_) + Un(c, n) * n + β * π[s_0] + @ xprime, + Uc(c, n) - μ * (Ucc(c, n) + * (c - B_) + Uc(c, n)) - Ξ, + Un(c, n) - μ * (Unn(c, n) * n + + Un(c, n)) + Θ[s_0] * Ξ, + (Θ * n - c - G)[s_0]]) + + # Find root + res = root(FOC, np.array( + [0, self.cFB[s_0], self.nFB[s_0], self.ΞFB[s_0]])) + if not res.success: + raise Exception('Could not find time 0 LS allocation.') + + return res.x + + def time1_value(self, μ): + ''' + Find the value associated with multiplier μ + ''' + c, n, x, Ξ = self.time1_allocation(μ) + U = self.model.U(c, n) + V = np.linalg.solve(np.eye(self.S) - self.β * self.π, U) + return c, n, x, V + + def Τ(self, c, n): + ''' + Computes Τ given c, n + ''' + model = self.model + Uc, Un = model.Uc(c, n), model.Un(c, n) + + return 1 + Un / (self.Θ * Uc) + + def simulate(self, B_, s_0, T, sHist=None): + ''' + Simulates planners policies for T periods + ''' + model, π, β = self.model, self.π, self.β + Uc = model.Uc + + if sHist is None: + sHist = self.mc.simulate(T, s_0) + + cHist, nHist, Bhist, ΤHist, μHist = np.zeros((5, T)) + RHist = np.zeros(T - 1) + + # Time 0 + μ, cHist[0], nHist[0], _ = self.time0_allocation(B_, s_0) + ΤHist[0] = self.Τ(cHist[0], nHist[0])[s_0] + Bhist[0] = B_ + μHist[0] = μ + + # Time 1 onward + for t in range(1, T): + c, n, x, Ξ = self.time1_allocation(μ) + Τ = self.Τ(c, n) + u_c = Uc(c, n) + s = sHist[t] + Eu_c = π[sHist[t - 1]] @ u_c + cHist[t], nHist[t], Bhist[t], ΤHist[t] = c[s], n[s], x[s] / u_c[s], \ + Τ[s] + RHist[t - 1] = Uc(cHist[t - 1], nHist[t - 1]) / (β * Eu_c) + μHist[t] = μ + + return np.array([cHist, nHist, Bhist, ΤHist, sHist, μHist, RHist]) diff --git a/source/_static/lecture_specific/amss/utilities.py b/source/_static/lecture_specific/amss2/utilities.py similarity index 100% rename from source/_static/lecture_specific/amss/utilities.py rename to source/_static/lecture_specific/amss2/utilities.py diff --git a/source/rst/amss.rst b/source/rst/amss.rst index d7e804f..f282415 100644 --- a/source/rst/amss.rst +++ b/source/rst/amss.rst @@ -14,6 +14,7 @@ In addition to what's in Anaconda, this lecture will need the following librarie :class: hide-output !pip install --upgrade quantecon + !pip install interpolation Overview ======== @@ -22,13 +23,15 @@ Let's start with following imports: .. code-block:: ipython - import numpy as np - import matplotlib.pyplot as plt - %matplotlib inline - from scipy.optimize import root, fmin_slsqp - from scipy.interpolate import UnivariateSpline - from quantecon import MarkovChain + import numpy as np + import matplotlib.pyplot as plt + from scipy.optimize import root + from interpolation.splines import eval_linear, UCGrid, nodes + from quantecon import optimize, MarkovChain + from numba import njit, prange, float64 + from numba.experimental import jitclass + %matplotlib inline In :doc:`an earlier lecture`, we described a model of optimal taxation with state-contingent debt due to @@ -170,7 +173,7 @@ yields: b_t(s^{t-1}) = z_t(s^t) + \beta \sum_{s^{t+1}\vert s^t} \pi_{t+1}(s^{t+1} | s^t) { u_c(s^{t+1}) \over u_c(s^{t}) } \; b_{t+1}(s^t) -Components of :math:`z_t(s^t)` on the right side depend on :math:`s^t`, but the left side is required to depend only +Components of :math:`z_t(s^t)` on the right side depend on :math:`s^t`, but the left side is required to depend only on :math:`s^{t-1}` . **This is what it means for one-period government debt to be risk-free**. @@ -717,14 +720,6 @@ Examples We now turn to some examples. - -We will first build some useful functions for solving the model - -.. literalinclude:: /_static/lecture_specific/amss/utilities.py - - - - Anticipated One-Period War ---------------------------------- @@ -804,59 +799,83 @@ Paths with circles are histories in which there is peace, while those with triangle denote war. +.. code-block:: python3 + + # WARNING: DO NOT EXPECT THE CODE TO WORK IF YOU CHANGE PARAMETERS + σ = 2 + γ = 2 + β = 0.9 + Π = np.array([[0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0], + [0, 0, 0, 0.5, 0.5, 0], + [0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 1]]) + g = np.array([0.1, 0.1, 0.1, 0.2, 0.1, 0.1]) + + x_min = -1.5555 + x_max = 17.339 + x_num = 300 + + x_grid = UCGrid((x_min, x_max, x_num)) + + crra_pref = CRRAutility(β=β, σ=σ, γ=γ) + + S = len(Π) + bounds_v = np.vstack([np.hstack([np.ones(S) * -10., np.zeros(S)]), + np.hstack([np.ones(S) - g, np.ones(S) * 10])]).T + + amss_model = AMSS(crra_pref, β, Π, g, x_grid, bounds_v) + .. code-block:: python3 - # Initialize μgrid for value function iteration - μ_grid = np.linspace(-0.7, 0.01, 300) + # WARNING: DO NOT EXPECT THE CODE TO WORK IF YOU CHANGE PARAMETERS + V = np.zeros((len(Π), x_num)) + V[:] = -nodes(x_grid).T ** 2 + + σ_v_star = np.ones((S, x_num, S * 2)) + σ_v_star[:, :, :S] = 0.0 - time_example = CRRAutility() + W = np.empty(len(Π)) + b_0 = 1.0 + σ_w_star = np.ones((S, 2)) + σ_w_star[:, 0] = -0.05 - time_example.π = np.array([[0, 1, 0, 0, 0, 0], - [0, 0, 1, 0, 0, 0], - [0, 0, 0, 0.5, 0.5, 0], - [0, 0, 0, 0, 0, 1], - [0, 0, 0, 0, 0, 1], - [0, 0, 0, 0, 0, 1]]) +.. code-block:: python3 - time_example.G = np.array([0.1, 0.1, 0.1, 0.2, 0.1, 0.1]) - time_example.Θ = np.ones(6) # Θ can in principle be random + %%time - time_example.transfers = True # Government can use transfers - # Solve sequential problem - time_sequential = SequentialAllocation(time_example) - # Solve recursive problem - time_bellman = RecursiveAllocationAMSS(time_example, μ_grid) + amss_model.solve(V, σ_v_star, b_0, W, σ_w_star) - sHist_h = np.array([0, 1, 2, 3, 5, 5, 5]) - sHist_l = np.array([0, 1, 2, 4, 5, 5, 5]) - sim_seq_h = time_sequential.simulate(1, 0, 7, sHist_h) - sim_bel_h = time_bellman.simulate(1, 0, 7, sHist_h) - sim_seq_l = time_sequential.simulate(1, 0, 7, sHist_l) - sim_bel_l = time_bellman.simulate(1, 0, 7, sHist_l) +.. code-block:: python3 - # Government spending paths - sim_seq_l[4] = time_example.G[sHist_l] - sim_seq_h[4] = time_example.G[sHist_h] - sim_bel_l[4] = time_example.G[sHist_l] - sim_bel_h[4] = time_example.G[sHist_h] + # Solve the LS model + ls_model = SequentialLS(crra_pref, g=g, π=Π) - # Output paths - sim_seq_l[5] = time_example.Θ[sHist_l] * sim_seq_l[1] - sim_seq_h[5] = time_example.Θ[sHist_h] * sim_seq_h[1] - sim_bel_l[5] = time_example.Θ[sHist_l] * sim_bel_l[1] - sim_bel_h[5] = time_example.Θ[sHist_h] * sim_bel_h[1] +.. code-block:: python3 + + # WARNING: DO NOT EXPECT THE CODE TO WORK IF YOU CHANGE PARAMETERS + s_hist_h = np.array([0, 1, 2, 3, 5, 5, 5]) + s_hist_l = np.array([0, 1, 2, 4, 5, 5, 5]) + + sim_h_amss = amss_model.simulate(s_hist_h, b_0) + sim_l_amss = amss_model.simulate(s_hist_l, b_0) + + sim_h_ls = ls_model.simulate(b_0, 0, 7, s_hist_h) + sim_l_ls = ls_model.simulate(b_0, 0, 7, s_hist_l) fig, axes = plt.subplots(3, 2, figsize=(14, 10)) titles = ['Consumption', 'Labor Supply', 'Government Debt', 'Tax Rate', 'Government Spending', 'Output'] - for ax, title, sim_l, sim_h, bel_l, bel_h in zip(axes.flatten(), titles, - sim_seq_l, sim_seq_h, - sim_bel_l, sim_bel_h): - ax.plot(sim_l, '-ok', sim_h, '-^k', bel_l, '-or', bel_h, '-^r', alpha=0.7) + for ax, title, ls_l, ls_h, amss_l, amss_h in zip(axes.flatten(), titles, + sim_l_ls, sim_h_ls, + sim_l_amss, sim_h_amss): + ax.plot(ls_l, '-ok', ls_h, '-^k', amss_l, '-or', amss_h, '-^r', + alpha=0.7) ax.set(title=title) ax.grid() @@ -864,8 +883,6 @@ triangle denote war. plt.show() - - How a Ramsey planner responds to war depends on the structure of the asset market. If it is able to trade state-contingent debt, then at time :math:`t=2` @@ -904,7 +921,7 @@ Without state-contingent debt, the optimal tax rate is history dependent. * A war at time :math:`t=3` causes a permanent **increase** in the tax rate. -* Peace at time :math:`t=3` causes a permanent **reduction** in the tax rate. +* Peace at time :math:`t=3` causes a permanent **reduction** in the tax rate. Perpetual War Alert @@ -937,37 +954,65 @@ state-contingent debt (circles) and the economy with only a risk-free bond (triangles). +.. code-block:: python3 + + # WARNING: DO NOT EXPECT THE CODE TO WORK IF YOU CHANGE PARAMETERS + ψ = 0.69 + Π = 0.5 * np.ones((2, 2)) + β = 0.9 + g = np.array([0.1, 0.2]) + + x_min = -3.4107 + x_max = 3.709 + x_num = 300 + + x_grid = UCGrid((x_min, x_max, x_num)) + log_pref = LogUtility(β=β, ψ=ψ) + + S = len(Π) + bounds_v = np.vstack([np.zeros(2 * S), np.hstack([1 - g, np.ones(S)]) ]).T + + V = np.zeros((len(Π), x_num)) + V[:] = -(nodes(x_grid).T + x_max) ** 2 / 14 + + σ_v_star = 1 - np.ones((S, x_num, S * 2)) * 0.55 + + W = np.empty(len(Π)) + b_0 = 0.5 + σ_w_star = 1 - np.ones((S, 2)) * 0.55 + + amss_model = AMSS(log_pref, β, Π, g, x_grid, bounds_v) + .. code-block:: python3 - log_example = LogUtility() - log_example.transfers = True # Government can use transfers - log_sequential = SequentialAllocation(log_example) # Solve sequential problem - log_bellman = RecursiveAllocationAMSS(log_example, μ_grid) + %%time - T = 20 - sHist = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, - 0, 0, 0, 1, 1, 1, 1, 1, 1, 0]) + amss_model.solve(V, σ_v_star, b_0, W, σ_w_star, tol_vfi=3e-5, maxitr=3000, + print_itr=100) - # Simulate - sim_seq = log_sequential.simulate(0.5, 0, T, sHist) - sim_bel = log_bellman.simulate(0.5, 0, T, sHist) +.. code-block:: python3 - titles = ['Consumption', 'Labor Supply', 'Government Debt', - 'Tax Rate', 'Government Spending', 'Output'] + ls_model = SequentialLS(log_pref, g=g, π=Π) # Solve sequential problem + +.. code-block:: python3 - # Government spending paths - sim_seq[4] = log_example.G[sHist] - sim_bel[4] = log_example.G[sHist] + # WARNING: DO NOT EXPECT THE CODE TO WORK IF YOU CHANGE PARAMETERS + s_hist = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, + 0, 0, 0, 1, 1, 1, 1, 1, 1, 0]) - # Output paths - sim_seq[5] = log_example.Θ[sHist] * sim_seq[1] - sim_bel[5] = log_example.Θ[sHist] * sim_bel[1] + T = len(s_hist) + + sim_amss = amss_model.simulate(s_hist, b_0) + sim_ls = ls_model.simulate(0.5, 0, T, s_hist) + + titles = ['Consumption', 'Labor Supply', 'Government Debt', + 'Tax Rate', 'Government Spending', 'Output'] fig, axes = plt.subplots(3, 2, figsize=(14, 10)) - for ax, title, seq, bel in zip(axes.flatten(), titles, sim_seq, sim_bel): - ax.plot(seq, '-ok', bel, '-^b') + for ax, title, ls, amss in zip(axes.flatten(), titles, sim_ls, sim_amss): + ax.plot(ls, '-ok', amss, '-^b') ax.set(title=title) ax.grid() @@ -995,30 +1040,28 @@ This outcome reflects the presence of a force for **precautionary saving** that In :doc:`this subsequent lecture` and :doc:`this subsequent lecture`, some ultimate consequences of that force are explored. +.. code-block:: python3 + + T = 200 + s_0 = 0 + mc = MarkovChain(Π) + + s_hist_long = mc.simulate(T, init=s_0, random_state=5) .. code-block:: python3 - T = 200 # Set T to 200 periods - sim_seq_long = log_sequential.simulate(0.5, 0, T) - sHist_long = sim_seq_long[-3] - sim_bel_long = log_bellman.simulate(0.5, 0, T, sHist_long) + sim_amss = amss_model.simulate(s_hist_long, b_0) + sim_ls = ls_model.simulate(0.5, 0, T, s_hist_long) titles = ['Consumption', 'Labor Supply', 'Government Debt', 'Tax Rate', 'Government Spending', 'Output'] - # Government spending paths - sim_seq_long[4] = log_example.G[sHist_long] - sim_bel_long[4] = log_example.G[sHist_long] - - # Output paths - sim_seq_long[5] = log_example.Θ[sHist_long] * sim_seq_long[1] - sim_bel_long[5] = log_example.Θ[sHist_long] * sim_bel_long[1] fig, axes = plt.subplots(3, 2, figsize=(14, 10)) - for ax, title, seq, bel in zip(axes.flatten(), titles, sim_seq_long, \ - sim_bel_long): - ax.plot(seq, '-k', bel, '-.b', alpha=0.5) + for ax, title, ls, amss in zip(axes.flatten(), titles, sim_ls, \ + sim_amss): + ax.plot(ls, '-k', amss, '-.b', alpha=0.5) ax.set(title=title) ax.grid() diff --git a/source/rst/amss2.rst b/source/rst/amss2.rst index 87fbf86..fe6c353 100644 --- a/source/rst/amss2.rst +++ b/source/rst/amss2.rst @@ -263,7 +263,7 @@ In equation :eq:`amss2_TS_barg10`, it is understood that :math:`c` and :math:`g` The CRRA utility function is represented in the following class. -.. literalinclude:: /_static/lecture_specific/opt_tax_recur/crra_utility.py +.. literalinclude:: /_static/lecture_specific/amss2/crra_utility.py Example Economy @@ -295,14 +295,14 @@ The code is mostly taken or adapted from the earlier lectures :doc:`optimal tax :doc:`optimal taxation with state-contingent debt`. -.. literalinclude:: /_static/lecture_specific/opt_tax_recur/sequential_allocation.py +.. literalinclude:: /_static/lecture_specific/amss2/sequential_allocation.py :class: collapse -.. literalinclude:: /_static/lecture_specific/amss/recursive_allocation.py +.. literalinclude:: /_static/lecture_specific/amss2/recursive_allocation.py :class: collapse -.. literalinclude:: /_static/lecture_specific/amss/utilities.py +.. literalinclude:: /_static/lecture_specific/amss2/utilities.py :class: collapse @@ -433,7 +433,7 @@ To solve the equations for :math:`c_0, b_0`, we use SciPy's fsolve function .. code-block:: python3 - c0, b0 = fsolve(solve_cb, np.array([1., -1.], dtype='float64'), + c0, b0 = fsolve(solve_cb, np.array([1., -1.], dtype='float64'), args=(Φ_star, b[0], 1), xtol=1.0e-12) c0, b0 @@ -790,4 +790,4 @@ Now let's compute the implied meantime to get to within 0.01 of the limit The slow rate of convergence and the implied time of getting within one percent of the limiting value do a good job of approximating our long simulation above. -In :doc:`a subsequent lecture` we shall study an extension of the model in which the force highlighted in this lecture causes government debt to converge to a nontrivial distribution instead of the single debt level discovered here. +In :doc:`a subsequent lecture` we shall study an extension of the model in which the force highlighted in this lecture causes government debt to converge to a nontrivial distribution instead of the single debt level discovered here. diff --git a/source/rst/amss3.rst b/source/rst/amss3.rst index 8e2d7f1..2372770 100644 --- a/source/rst/amss3.rst +++ b/source/rst/amss3.rst @@ -113,7 +113,7 @@ We set preference parameters The following Python code sets up the economy -.. literalinclude:: /_static/lecture_specific/opt_tax_recur/crra_utility.py +.. literalinclude:: /_static/lecture_specific/amss2/crra_utility.py First and Second Moments @@ -150,13 +150,13 @@ We begin by showing the code that we used in earlier lectures on the AMSS model. Here it is -.. literalinclude:: /_static/lecture_specific/opt_tax_recur/sequential_allocation.py +.. literalinclude:: /_static/lecture_specific/amss2/sequential_allocation.py :class: collapse -.. literalinclude:: /_static/lecture_specific/amss/recursive_allocation.py +.. literalinclude:: /_static/lecture_specific/amss2/recursive_allocation.py :class: collapse -.. literalinclude:: /_static/lecture_specific/amss/utilities.py +.. literalinclude:: /_static/lecture_specific/amss2/utilities.py :class: collapse Next, we show the code that we use to generate a very long simulation starting from initial @@ -215,7 +215,7 @@ We obtain the following graph for the histogram of the last 100,000 observations .. figure:: /_static/lecture_specific/amss3/amss3_g3.png The black vertical line denotes the sample mean for the last 100,000 observations included in the histogram; the green vertical line denotes the -value of :math:`\frac{ {\mathcal B}^*}{E u_c}`, associated with a sample from our approximation to +value of :math:`\frac{ {\mathcal B}^*}{E u_c}`, associated with a sample from our approximation to the ergodic distribution where :math:`{\mathcal B}^*` is a regression coefficient to be described below; the red vertical line denotes an approximation by :cite:`BEGS1` to the mean of the ergodic distribution that can be computed **before** the ergodic distribution has been approximated, as described below. @@ -315,7 +315,7 @@ BEGS interpret random variations in the right side of :eq:`eq_fiscal_risk_1` as Asymptotic Mean ---------------- -BEGS give conditions under which the ergodic mean of :math:`{\mathcal B}_t` is approximated by +BEGS give conditions under which the ergodic mean of :math:`{\mathcal B}_t` is approximated by .. math:: :label: prelim_formula_1 @@ -339,7 +339,7 @@ of equation :eq:`eq_fiscal_risk_1`. Expressing formula :eq:`prelim_formula_1` in terms of our notation tells us that the ergodic mean of the par value :math:`b` of government debt in the -AMSS model should be approximately +AMSS model should be approximately .. math:: :label: key_formula_1 @@ -577,7 +577,7 @@ But because we are studying an IID case, :math:`\pi` has identical rows and we need only to compute objects for one row of :math:`\pi`. This explains why at some places below we set :math:`s=0` just to pick -off the first row of :math:`\pi`. +off the first row of :math:`\pi`. Running the code