<jupyter_start><jupyter_text>The basic model: amplification and persistance<jupyter_code># 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
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]
params = {'a': 1.01, 'm': 10.0, 'alpha': 0.33, 'R': 1.5, 'K_bar': 100, 'K0': 95}
Kstar
B0
Bstar
# 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])
nodes = pycollocation.PolynomialSolver.collocation_nodes(5, domain, "Chebyshev")
problem = pycollocation.TwoPointBVP(bcs_lower, bcs_upper, 1, 4, rhs, params)
solution = pycollocation.PolynomialSolver.solve({'kind': "Chebyshev"},
                                                initial_coefs,
                                                domain,
                                                nodes,
                                                problem)
pycollocation.PolynomialSolver._array_to_list(initial_coefs, 4)
initial_capital_poly.coef
solution.result
K_hat, B_hat, q_hat, u_hat = solution.functions
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()
K_resids, B_resids, q_resids, u_resids = solution.residuals(pts)
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()
basic_model_solver.solve(kind="Chebyshev",
                         coefs_dict=initial_coefs,
                         domain=domain)

basic_model_solver.result["success"]
basic_model_viz = pycollocation.Visualizer(basic_model_solver)
basic_model_viz.interpolation_knots = np.linspace(domain[0], domain[1], 1000)
basic_model_viz.solution.plot(subplots=True, style=['r', 'b'])
plt.show()<jupyter_output><empty_output><jupyter_text>Solution is not as accurate as I would like...<jupyter_code>basic_model_viz.residuals.plot(subplots=True, style=['r', 'b'])
plt.show()<jupyter_output><empty_output><jupyter_text>...actually, when using noramlized residuals eveything looks great!<jupyter_code>basic_model_viz.normalized_residuals.plot(logy=True, sharey=True)
plt.show()
assets = basic_model_viz.solution[['q', 'K']].prod(axis=1)
liabilities = basic_model_viz.solution.B
equity = assets - liabilities
leverage = assets / equity
leverage.plot()
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)

initial_condition = np.array([Kstar(a, m, alpha, R, K_bar), Bstar(a, m, alpha, R, K_bar), (R / (R - 1)) * a, a])
initial_condition
credit_cycles(0, initial_condition)
jacobian(0, initial_condition)
from scipy import linalg
from IPython.html.widgets import fixed, interact, FloatSliderWidget
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
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))
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<jupyter_output><empty_output><jupyter_text>Full model<jupyter_code>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}
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]
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}
# 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)
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}
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)
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))
full_model_solver.solve(kind="Chebyshev",
                        coefs_dict=initial_coefs,
                        domain=domain)

full_model_solver.result["success"]
full_model_viz = pycollocation.Visualizer(full_model_solver)
full_model_viz.interpolation_knots = np.linspace(domain[0], domain[1], 1000)
full_model_viz.solution.plot(subplots=True)
full_model_viz.normalized_residuals.plot(subplots=True)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1822775986756121, 'fraction_numerical': 0.03406813627254509, 'mean_word_length': 2.660076530612245, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '511934', 'license': 'MIT', 'n_tokens_mistral': 5010, 'n_tokens_neox': 4701, 'n_words': 1182, 'path': 'examples/credit-cycles.ipynb', 'repo_name': 'davidrpugh/bvp-solver', 'size': '11477'}