In [1]:
%load_ext autoreload
In [2]:
autoreload 2
In [3]:
%matplotlib inline
In [4]:
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
import seaborn as sn
import pycollocation
In [7]:
# need to specify some production function for gatherers
def output(k, alpha):
return k**alpha
def marginal_product_capital(k, alpha, **params):
return alpha * k**(alpha - 1)
def K_dot(t, K, B, q, u, a, R):
return (1 / u) * ((a + q) * K - R * B) - K
def B_dot(t, K, B, q, u, R):
return (1 / R) * q * K - B
def q_dot(t, K, B, q, u, R):
return (R - 1) * q - R * u
def u_dot(t, K, B, q, u, m, K_bar, R, **params):
return (1 / R) * marginal_product_capital((1 / m) * (K_bar - K), **params) - u
def rhs(t, K, B, q, u, a, m, K_bar, R, **params):
out = [K_dot(t, K, B, q, u, a, R),
B_dot(t, K, B, q, u, R),
q_dot(t, K, B, q, u, R),
u_dot(t, K, B, q, u, m, K_bar, R, **params)]
return out
In [82]:
def steady_state_capital(a, m, K_bar, R, alpha, **params):
return K_bar - m * (alpha / (a * R))**(1 / (1 - alpha))
def steady_state_debt(a, m, K_bar, R, **params):
Kstar = steady_state_capital(a, m, K_bar, R, **params)
return (a / (R - 1)) * Kstar
def steady_state_land_price(a, R, **params):
return (R / (R - 1)) * a
def steady_state_user_cost(a, **params):
return a
def bcs_lower(t, K, B, q, u, K0, **params):
return [K - K0]
def bcs_upper(t, K, B, q, u, a, m, K_bar, R, **params):
Bstar = steady_state_debt(a, m, K_bar, R, **params)
qstar = steady_state_land_price(a, R)
ustar = steady_state_user_cost(a)
return [B - Bstar, q - qstar, u - ustar]
In [83]:
params = {'a': 1.01, 'm': 10.0, 'alpha': 0.33, 'R': 1.5, 'K_bar': 100, 'K0': 95}
In [84]:
Kstar
Out[84]:
In [85]:
B0
Out[85]:
In [86]:
Bstar
Out[86]:
In [93]:
# specify an initial guess
domain = [0, 10]
ts = np.linspace(domain[0], domain[1], 1000)
Kstar = steady_state_capital(**params)
Ks = Kstar - (Kstar - params['K0']) * np.exp(-ts)
initial_capital_poly = np.polynomial.Chebyshev.fit(ts, Ks, 5, domain)
# initial value of debt is some multiple of capital stock
B0 = 1.5 * params['K0']
Bstar = steady_state_debt(**params)
Bs = Bstar - (Bstar - B0) * np.exp(-ts)
initial_debt_poly = np.polynomial.Chebyshev.fit(ts, Bs, 5, domain)
# starting with K0 < Kstar, must be that u0 > ustar
ustar = steady_state_user_cost(**params)
u0 = 1.5 * ustar
us = ustar - (ustar - u0) * np.exp(-ts)
initial_user_cost_poly = np.polynomial.Chebyshev.fit(ts, us, 5, domain)
# starting with K0 < Kstar, must be that q0 > qstar
qstar = steady_state_land_price(**params)
q0 = 1.05 * qstar
qs = qstar + (qstar - q0) * np.exp(-ts)
initial_land_price_poly = np.polynomial.Chebyshev.fit(ts, qs, 5, domain)
initial_coefs = np.hstack([initial_capital_poly.coef, initial_debt_poly.coef,
initial_user_cost_poly.coef, initial_land_price_poly.coef])
In [94]:
nodes = pycollocation.PolynomialSolver.collocation_nodes(5, domain, "Chebyshev")
In [95]:
problem = pycollocation.TwoPointBVP(bcs_lower, bcs_upper, 1, 4, rhs, params)
In [96]:
solution = pycollocation.PolynomialSolver.solve({'kind': "Chebyshev"},
initial_coefs,
domain,
nodes,
problem)
In [99]:
pycollocation.PolynomialSolver._array_to_list(initial_coefs, 4)
Out[99]:
In [100]:
initial_capital_poly.coef
Out[100]:
In [97]:
solution.result
Out[97]:
In [72]:
K_hat, B_hat, q_hat, u_hat = solution.functions
In [73]:
pts = np.linspace(domain[0], domain[1], 1000)
fig, axes = plt.subplots(4, 1)
axes[0].plot(pts, K_hat(pts))
axes[1].plot(pts, B_hat(pts))
axes[2].plot(pts, q_hat(pts))
axes[3].plot(pts, q_hat(pts))
fig.tight_layout()
plt.show()
In [74]:
K_resids, B_resids, q_resids, u_resids = solution.residuals(pts)
In [75]:
pts = np.linspace(domain[0], domain[1], 1000)
fig, axes = plt.subplots(4, 1)
axes[0].plot(pts, K_resids)
axes[1].plot(pts, B_resids)
axes[2].plot(pts, q_resids)
axes[3].plot(pts, u_resids)
fig.tight_layout()
plt.show()
In [10]:
basic_model_solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain)
In [11]:
basic_model_solver.result["success"]
Out[11]:
In [12]:
basic_model_viz = pycollocation.Visualizer(basic_model_solver)
In [13]:
basic_model_viz.interpolation_knots = np.linspace(domain[0], domain[1], 1000)
basic_model_viz.solution.plot(subplots=True, style=['r', 'b'])
plt.show()
Solution is not as accurate as I would like...
In [14]:
basic_model_viz.residuals.plot(subplots=True, style=['r', 'b'])
plt.show()
...actually, when using noramlized residuals eveything looks great!
In [16]:
basic_model_viz.normalized_residuals.plot(logy=True, sharey=True)
plt.show()
In [17]:
assets = basic_model_viz.solution[['q', 'K']].prod(axis=1)
liabilities = basic_model_viz.solution.B
equity = assets - liabilities
leverage = assets / equity
In [18]:
leverage.plot()
Out[18]:
In [165]:
def credit_cycles(t, X, a, m, alpha, R, K_bar):
out = np.array([(1 / X[3]) * ((a + X[2]) * X[0] - R * X[1]) - X[0],
(1 / R) * X[2] * X[0] - X[1],
(R - 1) * X[2] - R * X[3],
(alpha / R) * ((1 / m) * (K_bar - X[0]))**(alpha - 1) - X[3]])
return out
def jacobian(t, X, a, m, alpha, R, K_bar):
out = np.array([[((a + X[2]) / X[3]) - 1.0, -R / X[3], X[0] / X[3], -X[3]**(-2)],
[(1 / R) * X[2], -1.0, (1 / R) * X[0], 0.0],
[0.0, 0.0, R - 1, -R],
[-(1 / m) * (alpha - 1) * (alpha / R) * ((1 / m) * (K_bar - X[0]))**(alpha - 2), 0.0, 0.0, -1.0]])
return out
def Kstar(a, m, alpha, R, K_bar):
return K_bar - m * (alpha / (a * R))**(1 / (1 - alpha))
def Bstar(a, m, alpha, R, K_bar):
return (a / (R - 1)) * Kstar(a, m, alpha, R, K_bar)
In [166]:
initial_condition = np.array([Kstar(a, m, alpha, R, K_bar), Bstar(a, m, alpha, R, K_bar), (R / (R - 1)) * a, a])
In [167]:
initial_condition
Out[167]:
In [901]:
credit_cycles(0, initial_condition)
Out[901]:
In [868]:
jacobian(0, initial_condition)
Out[868]:
In [98]:
from scipy import linalg
In [99]:
from IPython.html.widgets import fixed, interact, FloatSliderWidget
In [170]:
def eigenvalues(a=1.0, m=1.0, alpha=0.33, R=1.05, K_bar=10.0):
steady_state = np.array([Kstar(a, m, alpha, R, K_bar),
Bstar(a, m, alpha, R, K_bar),
(R / (R - 1)) * a,
a])
vals, vecs = linalg.eig(jacobian(0, steady_state, a, m, alpha, R, K_bar))
print vals
In [172]:
interact(eigenvalues, a=(0.0, 1e3, 1e0), m=(0.0, 1e2, 1e-1), R=(0.0, 1e2, 1e-2), K_bar=(0.0, 1e4, 1e1))
In [1120]:
params = 2.0, 0.5, 0.33, 1.01, 500.0
problem = ivp.IVP(credit_cycles, jacobian)
problem.f_params = params
problem.jac_params = params
In [6]:
lamda, pi, phi, = sym.symbols('lamda, pi, phi')
# full model from Kiyotaki and Moore "credit-cycles" paper
K_dot = (pi / (phi + u)) * ((a + q + lamda * phi) * K - R * B) - pi * lamda * K
B_dot = (R - 1) * B + (phi * (1 - lamda) - a) * K
q_dot = (R - 1) * q - R * u
u_dot = (1 / R) * mpk.subs({k: (1 / m) * (K_bar - K)}) - u
rhs = {'K': K_dot, 'B': B_dot, 'q': q_dot, 'u': u_dot}
In [85]:
bcs = {}
ustar = ((pi * a - (1 - lamda) * (1 - R + pi * R) * phi) /
(lamda * pi + (1 - lamda) * (1 - R + pi * R)))
qstar = (R / (R - 1)) * ustar
Kstar = K_bar - m * (alpha / (ustar * R))**(1 / (1 - alpha))
Bstar = ((a - (1 - lamda) * phi) / (R - 1)) * Kstar
# initial conditions for K and B are given
K0 = 75
bcs['lower'] = [K - K0]
# boundary conditions on B, q and u can be written in terms of steady state values
bcs['upper'] = [B - Bstar, q - qstar, u - ustar]
In [116]:
params = {'a': 1.05, 'pi': 0.05, 'phi': 20.0, 'lamda': 0.975,'m': 1.0, 'alpha': 0.16,
'R': 1.01, 'K_bar': 100}
In [117]:
# set up the model and solver
full_model = pycollocation.SymbolicBoundaryValueProblem(dependent_vars=['K', 'B', 'q', 'u'],
independent_var='t',
rhs=rhs,
boundary_conditions=bcs,
params=params)
full_model_solver = pycollocation.OrthogonalPolynomialSolver(full_model)
In [125]:
def Kstar(a, phi, R, alpha, pi, m, lamda, K_bar):
return K_bar - m * (alpha / (ustar(a, phi, R, alpha, pi, m, lamda, K_bar) * R))**(1 / (1 - alpha))
def Bstar(a, phi, R, alpha, pi, m, lamda, K_bar):
return ((a - (1 - lamda) * phi) / (R - 1)) * Kstar(a, phi, R, alpha, pi, m, lamda, K_bar)
def qstar(a, phi, R, alpha, pi, m, lamda, K_bar):
return (R / (R - 1)) * ustar(a, phi, R, alpha, pi, m, lamda, K_bar)
def ustar(a, phi, R, alpha, pi, m, lamda, K_bar):
u = ((pi * a - (1 - lamda) * (1 - R + pi * R) * phi) /
(lamda * pi + (1 - lamda) * (1 - R + pi * R)))
return u
# specify an initial guess
domain = [0, 25]
ts = np.linspace(domain[0], domain[1], 1000)
Ks = Kstar(**params) - (Kstar(**params) - K0) * np.exp(-ts) * np.cos(2.0 * np.pi * ts)
initial_capital_poly = np.polynomial.Chebyshev.fit(ts, Ks, 25, domain)
# initial value of debt is some multiple of capital stock
B0 = 1.5 * K0
Bs = Bstar(**params) - (Bstar(**params) - B0) * np.exp(-ts) #* np.cos(2.0 * np.pi * ts)
initial_debt_poly = np.polynomial.Chebyshev.fit(ts, Bs, 25, domain)
# starting with K0 > Kstar, must be that u0 > ustar
us = ustar(**params) - (ustar(**params) - 1.5 * ustar(**params)) * np.exp(-ts) #* np.cos(2.0 * np.pi * ts)
initial_user_cost_poly = np.polynomial.Chebyshev.fit(ts, us, 25, domain)
# starting with K0 > Kstar, must be that q0 > qstar
qs = qstar(**params) - (qstar(**params) - 1.5 * qstar(**params)) * np.exp(-ts) #* np.cos(2.0 * np.pi * ts)
initial_land_price_poly = np.polynomial.Chebyshev.fit(ts, qs, 25, domain)
initial_coefs = {'K': initial_capital_poly.coef, 'B': initial_debt_poly.coef,
'u': initial_user_cost_poly.coef, 'q': initial_land_price_poly.coef}
In [119]:
def jacobian(t, X, a, phi, R, alpha, pi, m, lamda, K_bar):
out = np.array([[(pi / (phi + X[3])) * (a + X[2] + lamda * phi) - pi * lamda, -(pi / (phi + X[3])) * R, (pi / (phi + X[3])) * X[0], -(pi / (phi + X[3])**2) * ((a + X[2] + lamda * phi) * X[0] - R * X[1])],
[(R - 1) * X[1] + (phi * (1 - lamda) - a), (R - 1), 0.0, 0.0],
[0.0, 0.0, R - 1, -R],
[-(1 / m) * (alpha - 1) * (alpha / R) * ((1 / m) * (K_bar - X[0]))**(alpha - 2), 0.0, 0.0, -1.0]])
return out
def eigenvalues(a=1.0, phi=20.0, pi=0.05, lamda=0.975, m=1.0, alpha=0.33, R=1.05, K_bar=10.0):
steady_state = np.array([Kstar(a, phi, R, alpha, pi, m, lamda, K_bar),
Bstar(a, phi, R, alpha, pi, m, lamda, K_bar),
qstar(a, phi, R, alpha, pi, m, lamda, K_bar),
ustar(a, phi, R, alpha, pi, m, lamda, K_bar)])
vals, vecs = linalg.eig(jacobian(0, steady_state, a, phi, R, alpha, pi, m, lamda, K_bar))
print vals
print np.absolute(vals)
In [120]:
interact(eigenvalues, a=(1.0, 2.0, 1e-2), alpha=(1e-2, 1-1e-2, 1e-2), m=(0.0, 1e2, 1e-1), R=(0.0, 1e2, 1e-2), K_bar=(0.0, 1e4, 1e1))
Out[120]:
In [ ]:
In [126]:
full_model_solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain)
In [127]:
full_model_solver.result["success"]
Out[127]:
In [123]:
full_model_viz = pycollocation.Visualizer(full_model_solver)
In [124]:
full_model_viz.interpolation_knots = np.linspace(domain[0], domain[1], 1000)
full_model_viz.solution.plot(subplots=True)
Out[124]:
In [52]:
full_model_viz.normalized_residuals.plot(subplots=True)
plt.show()
In [ ]: