import sys import numpy as np from numpy import sqrt, eye, zeros, cumsum from numpy.random import randn import scipy.linalg import matplotlib.pyplot as plt from collections import namedtuple from quantecon import nullspace, mc_sample_path, var_quadratic_sum # == Set up a namedtuple to store data on the model economy == # Economy = namedtuple('economy', ('β', # Discount factor 'Sg', # Govt spending selector matrix 'Sd', # Exogenous endowment selector matrix 'Sb', # Utility parameter selector matrix 'Ss', # Coupon payments selector matrix 'discrete', # Discrete or continuous -- boolean 'proc')) # Stochastic process parameters # == Set up a namedtuple to store return values for compute_paths() == # Path = namedtuple('path', ('g', # Govt spending 'd', # Endowment 'b', # Utility shift parameter 's', # Coupon payment on existing debt 'c', # Consumption 'l', # Labor 'p', # Price 'τ', # Tax rate 'rvn', # Revenue 'B', # Govt debt 'R', # Risk free gross return 'π', # One-period risk-free interest rate 'Π', # Cumulative rate of return, adjusted 'ξ')) # Adjustment factor for Π def compute_paths(T, econ): """ Compute simulated time paths for exogenous and endogenous variables. Parameters =========== T: int Length of the simulation econ: a namedtuple of type 'Economy', containing β - Discount factor Sg - Govt spending selector matrix Sd - Exogenous endowment selector matrix Sb - Utility parameter selector matrix Ss - Coupon payments selector matrix discrete - Discrete exogenous process (True or False) proc - Stochastic process parameters Returns ======== path: a namedtuple of type 'Path', containing g - Govt spending d - Endowment b - Utility shift parameter s - Coupon payment on existing debt c - Consumption l - Labor p - Price τ - Tax rate rvn - Revenue B - Govt debt R - Risk free gross return π - One-period risk-free interest rate Π - Cumulative rate of return, adjusted ξ - Adjustment factor for Π The corresponding values are flat numpy ndarrays. """ # == Simplify names == # β, Sg, Sd, Sb, Ss = econ.β, econ.Sg, econ.Sd, econ.Sb, econ.Ss if econ.discrete: P, x_vals = econ.proc else: A, C = econ.proc # == Simulate the exogenous process x == # if econ.discrete: state = mc_sample_path(P, init=0, sample_size=T) x = x_vals[:, state] else: # == Generate an initial condition x0 satisfying x0 = A x0 == # nx, nx = A.shape x0 = nullspace((eye(nx) - A)) x0 = -x0 if (x0[nx-1] < 0) else x0 x0 = x0 / x0[nx-1] # == Generate a time series x of length T starting from x0 == # nx, nw = C.shape x = zeros((nx, T)) w = randn(nw, T) x[:, 0] = x0.T for t in range(1, T): x[:, t] = A @ x[:, t-1] + C @ w[:, t] # == Compute exogenous variable sequences == # g, d, b, s = ((S @ x).flatten() for S in (Sg, Sd, Sb, Ss)) # == Solve for Lagrange multiplier in the govt budget constraint == # # In fact we solve for ν = lambda / (1 + 2*lambda). Here ν is the # solution to a quadratic equation a(ν**2 - ν) + b = 0 where # a and b are expected discounted sums of quadratic forms of the state. Sm = Sb - Sd - Ss # == Compute a and b == # if econ.discrete: ns = P.shape[0] F = scipy.linalg.inv(eye(ns) - β * P) a0 = 0.5 * (F @ (x_vals.T @ Sm.T)**2)[0] H = ((Sb - Sd + Sg) @ x_vals) * ((Sg - Ss) @ x_vals) b0 = 0.5 * (F @ H.T)[0] a0, b0 = float(a0), float(b0) else: H = Sm.T @ Sm a0 = 0.5 * var_quadratic_sum(A, C, H, β, x0) H = (Sb - Sd + Sg).T @ (Sg + Ss) b0 = 0.5 * var_quadratic_sum(A, C, H, β, x0) # == Test that ν has a real solution before assigning == # warning_msg = """ Hint: you probably set government spending too {}. Elect a {} Congress and start over. """ disc = a0**2 - 4 * a0 * b0 if disc >= 0: ν = 0.5 * (a0 - sqrt(disc)) / a0 else: print("There is no Ramsey equilibrium for these parameters.") print(warning_msg.format('high', 'Republican')) sys.exit(0) # == Test that the Lagrange multiplier has the right sign == # if ν * (0.5 - ν) < 0: print("Negative multiplier on the government budget constraint.") print(warning_msg.format('low', 'Democratic')) sys.exit(0) # == Solve for the allocation given ν and x == # Sc = 0.5 * (Sb + Sd - Sg - ν * Sm) Sl = 0.5 * (Sb - Sd + Sg - ν * Sm) c = (Sc @ x).flatten() l = (Sl @ x).flatten() p = ((Sb - Sc) @ x).flatten() # Price without normalization τ = 1 - l / (b - c) rvn = l * τ # == Compute remaining variables == # if econ.discrete: H = ((Sb - Sc) @ x_vals) * ((Sl - Sg) @ x_vals) - (Sl @ x_vals)**2 temp = (F @ H.T).flatten() B = temp[state] / p H = (P[state, :] @ x_vals.T @ (Sb - Sc).T).flatten() R = p / (β * H) temp = ((P[state, :] @ x_vals.T @ (Sb - Sc).T)).flatten() ξ = p[1:] / temp[:T-1] else: H = Sl.T @ Sl - (Sb - Sc).T @ (Sl - Sg) L = np.empty(T) for t in range(T): L[t] = var_quadratic_sum(A, C, H, β, x[:, t]) B = L / p Rinv = (β * ((Sb - Sc) @ A @ x)).flatten() / p R = 1 / Rinv AF1 = (Sb - Sc) @ x[:, 1:] AF2 = (Sb - Sc) @ A @ x[:, :T-1] ξ = AF1 / AF2 ξ = ξ.flatten() π = B[1:] - R[:T-1] * B[:T-1] - rvn[:T-1] + g[:T-1] Π = cumsum(π * ξ) # == Prepare return values == # path = Path(g=g, d=d, b=b, s=s, c=c, l=l, p=p, τ=τ, rvn=rvn, B=B, R=R, π=π, Π=Π, ξ=ξ) return path def gen_fig_1(path): """ The parameter is the path namedtuple returned by compute_paths(). See the docstring of that function for details. """ T = len(path.c) # == Prepare axes == # num_rows, num_cols = 2, 2 fig, axes = plt.subplots(num_rows, num_cols, figsize=(14, 10)) plt.subplots_adjust(hspace=0.4) for i in range(num_rows): for j in range(num_cols): axes[i, j].grid() axes[i, j].set_xlabel('Time') bbox = (0., 1.02, 1., .102) legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'} p_args = {'lw': 2, 'alpha': 0.7} # == Plot consumption, govt expenditure and revenue == # ax = axes[0, 0] ax.plot(path.rvn, label=r'$\tau_t \ell_t$', **p_args) ax.plot(path.g, label='$g_t$', **p_args) ax.plot(path.c, label='$c_t$', **p_args) ax.legend(ncol=3, **legend_args) # == Plot govt expenditure and debt == # ax = axes[0, 1] ax.plot(list(range(1, T+1)), path.rvn, label=r'$\tau_t \ell_t$', **p_args) ax.plot(list(range(1, T+1)), path.g, label='$g_t$', **p_args) ax.plot(list(range(1, T)), path.B[1:T], label='$B_{t+1}$', **p_args) ax.legend(ncol=3, **legend_args) # == Plot risk free return == # ax = axes[1, 0] ax.plot(list(range(1, T+1)), path.R - 1, label='$R_t - 1$', **p_args) ax.legend(ncol=1, **legend_args) # == Plot revenue, expenditure and risk free rate == # ax = axes[1, 1] ax.plot(list(range(1, T+1)), path.rvn, label=r'$\tau_t \ell_t$', **p_args) ax.plot(list(range(1, T+1)), path.g, label='$g_t$', **p_args) axes[1, 1].plot(list(range(1, T)), path.π, label=r'$\pi_{t+1}$', **p_args) ax.legend(ncol=3, **legend_args) plt.show() def gen_fig_2(path): """ The parameter is the path namedtuple returned by compute_paths(). See the docstring of that function for details. """ T = len(path.c) # == Prepare axes == # num_rows, num_cols = 2, 1 fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 10)) plt.subplots_adjust(hspace=0.5) bbox = (0., 1.02, 1., .102) bbox = (0., 1.02, 1., .102) legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'} p_args = {'lw': 2, 'alpha': 0.7} # == Plot adjustment factor == # ax = axes[0] ax.plot(list(range(2, T+1)), path.ξ, label=r'$\xi_t$', **p_args) ax.grid() ax.set_xlabel('Time') ax.legend(ncol=1, **legend_args) # == Plot adjusted cumulative return == # ax = axes[1] ax.plot(list(range(2, T+1)), path.Π, label=r'$\Pi_t$', **p_args) ax.grid() ax.set_xlabel('Time') ax.legend(ncol=1, **legend_args) plt.show()