%load_ext autoreload
autoreload 2
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
import seaborn as sn
import pycollocation
# 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
98.97176708265565
B0
142.5
Bstar
199.9229695069644
# 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)
[array([ 98.25192341, 1.29704214, -0.91771876, 0.5470026 , -0.25019508, 0.10814735]), array([ 189.51562177, 18.75236135, -13.26818405, 7.90844809, -3.61726759, 1.56357157]), array([ 1.10152628, -0.16491558, 0.11668559, -0.06954998, 0.03181166, -0.01375066]), array([ 3.00254212, 0.04947467, -0.03500568, 0.02086499, -0.0095435 , 0.0041252 ])]
initial_capital_poly.coef
array([ 98.25192341, 1.29704214, -0.91771876, 0.5470026 , -0.25019508, 0.10814735])
solution.result
status: 5 success: False qtf: array([ 118.11507587, -73.62805396, 34.42542314, 91.23327592, -21.58362966, -46.31673836, -312.73178766, 49.73195376, -63.00744997, 17.69122132, -29.50229648, 6.74219956, -130.3121271 , 76.17524295, -0.89808815, 1.0618052 , 4.798498 , 8.7009872 , 8.70516571, 13.22207324, 3.44698481, -3.86957298, 10.27331385, 4.00065662]) nfev: 59 r: array([ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]) fun: array([ 9.91315830e+01, 1.27135105e+02, 1.31413999e+02, 1.32020647e+02, 1.32290429e+02, 1.05707121e+02, 1.30513462e+02, 1.32691054e+02, 1.32959922e+02, 1.35856846e+02, 3.31235599e+00, 3.91420017e+00, 4.02867648e+00, 4.04545691e+00, 4.02224183e+00, 2.61978349e+00, 2.24958520e+00, 2.05253856e+00, 2.01347433e+00, 2.01616608e+00, 1.31817469e-01, 9.31581637e-01, -2.02819269e+00, 2.02245781e+00]) x: array([ 9.82519234e+01, 1.29704214e+00, -9.17718765e-01, 5.47002603e-01, -2.50195078e-01, 1.08147352e-01, 1.89515622e+02, 1.87523613e+01, -1.32681840e+01, 7.90844809e+00, -3.61726759e+00, 1.56357157e+00, 1.10152628e+00, -1.64915583e-01, 1.16685588e-01, -6.95499784e-02, 3.18116626e-02, -1.37506585e-02, 3.00254212e+00, 4.94746748e-02, -3.50056763e-02, 2.08649935e-02, -9.54349878e-03, 4.12519754e-03]) message: 'The iteration is not making good progress, as measured by the \n improvement from the last ten iterations.' fjac: array([[ -7.23449284e-02, -1.30682715e-01, -1.39759317e-01, -1.41050339e-01, -1.41202692e-01, 3.95426739e-01, 3.01616125e-01, 2.86438732e-01, 2.84266727e-01, 2.84010180e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.73561927e-02, 1.47345086e-01, 2.56011589e-01, 2.81930285e-01, 2.85267356e-01, -4.22620469e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ -2.66774214e-02, -2.23107591e-02, -1.19285758e-01, -2.23129477e-01, -2.87427339e-01, -4.32154370e-01, -1.92558236e-01, 2.77589130e-02, 2.36578697e-01, 3.65436465e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.98970708e-02, -9.40681674e-02, 2.48102042e-02, 2.34633815e-01, 3.67053739e-01, 4.87750569e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ -2.97178802e-01, -2.85558434e-01, -1.72785009e-01, 1.52402057e-01, 4.59196777e-01, -2.14202659e-01, 1.87158167e-01, 3.72466946e-01, 1.41007900e-01, -2.11049032e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -1.48188023e-02, 9.14298210e-02, 3.32901526e-01, 1.39848677e-01, -2.11983431e-01, 3.13422790e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 5.28146766e-01, 1.95267774e-01, -2.43643372e-01, -1.07104666e-01, 5.24213406e-01, 1.67066461e-02, -3.38591405e-01, -1.94745874e-02, 2.91408159e-01, -2.19860019e-02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.15586606e-03, -1.65408163e-01, -1.74058898e-02, 2.89012865e-01, -2.20825123e-02, -1.57839943e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 5.52677579e-01, -2.21286529e-01, -1.69796027e-01, 3.48154980e-01, -4.44166173e-01, -1.07787892e-01, -1.27528422e-01, 3.43715089e-01, -1.21986510e-01, -1.06626072e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -7.45696437e-03, -6.22993180e-02, 3.07202752e-01, -1.20984683e-01, -1.07098758e-01, -3.95499428e-02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ -4.08984368e-01, 4.90842459e-01, -4.45438456e-01, 3.35336614e-01, -2.94714927e-01, 1.46249558e-01, -2.21291912e-01, -5.99810589e-03, 1.67251343e-01, -1.67153214e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.01176778e-02, -1.08105235e-01, -5.36096591e-03, 1.65879718e-01, -1.67895451e-01, -1.58579220e-02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ -1.70265717e-01, -3.19460977e-01, -2.73805037e-01, -3.10048709e-01, -7.61346773e-02, -1.85920395e-01, -2.55209155e-01, -2.65093370e-01, -2.67416326e-01, -3.32673475e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.70557503e-02, 8.65885515e-02, 1.49584906e-01, 1.63683377e-01, 1.00225145e-01, -2.93057307e-01, -4.32456489e-01, 0.00000000e+00, 0.00000000e+00], [ 2.47590809e-01, 2.74408382e-01, 2.45881859e-02, -2.11383274e-01, -9.35247906e-02, 3.25113573e-01, 3.49629266e-01, 9.61799178e-02, -1.05265609e-01, -3.79273507e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.47760040e-02, -5.19934032e-02, 6.89406236e-03, 1.72431489e-01, 1.27599927e-01, 4.45254886e-01, -4.11842488e-01, 0.00000000e+00, 0.00000000e+00], [ 9.37327963e-02, -1.46606110e-01, -4.01193779e-01, -1.65580856e-02, 7.06490470e-02, 3.80745437e-01, 1.58479795e-02, -4.44632021e-01, -2.83474846e-01, 3.35637525e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -1.18334213e-02, 1.28269849e-01, 2.71991467e-01, -8.46985520e-03, -1.71263185e-01, 3.32432157e-01, 1.92577229e-01, 0.00000000e+00, 0.00000000e+00], [ -2.83414228e-02, 3.26163402e-01, -1.46957875e-01, -1.18446516e-01, 9.00695487e-02, -3.16900138e-01, 3.03232421e-01, 2.06809388e-01, -5.56230375e-01, 3.17559118e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -7.13296098e-03, -3.32661064e-01, 6.27590573e-02, 1.63904858e-01, -1.42600460e-01, -2.03965721e-01, -1.54617258e-02, 0.00000000e+00, 0.00000000e+00], [ 2.34564474e-02, -1.38943727e-01, -1.27915550e-01, 4.62546429e-02, -4.78450016e-02, -1.27137189e-01, 4.33101023e-01, -4.14795567e-01, 3.12863768e-01, -1.88015531e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -5.11635946e-02, -6.02898633e-01, 1.55122842e-01, -1.74006048e-01, 1.01433347e-01, -9.22605186e-02, 1.05250693e-01, 0.00000000e+00, 0.00000000e+00], [ 2.19519761e-01, 2.61633451e-01, -1.05877212e-01, 2.89874769e-02, -1.07016927e-01, -4.13756530e-01, 3.74783977e-01, -2.76762319e-01, 1.49644156e-01, -1.03647216e-01, 1.11022302e-16, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 6.04121077e-02, 6.20930332e-01, 1.31972814e-02, 1.15730685e-01, -1.08153678e-01, -1.01030969e-01, 1.17425727e-01, 0.00000000e+00, 0.00000000e+00], [ -2.39397034e-02, -1.16204377e-02, -1.01631154e-02, -2.98988479e-01, -8.70629369e-02, 4.89410703e-02, -7.07384754e-02, 1.18028512e-01, -2.14251036e-01, -4.29743435e-01, -9.92360879e-03, -9.92360879e-03, -9.92360879e-03, -9.92360879e-03, -9.92360879e-03, -2.96430220e-03, -3.07783560e-02, -4.18985610e-03, 1.98638147e-01, 1.67680810e-01, 1.45585961e-02, 7.62007751e-01, 1.98472181e-02, 0.00000000e+00], [ -2.58930933e-02, 1.68292599e-01, 3.62741618e-01, 4.68711891e-01, 1.51572398e-01, -2.12050715e-02, -9.69750249e-02, -2.20265277e-01, -1.95251997e-01, -2.18331324e-02, -1.45873117e-01, -1.21836806e-01, -8.29452362e-02, -4.40536423e-02, -2.00173794e-02, -9.23296414e-03, 1.19499959e-02, 4.29331555e-01, 2.63268276e-01, 4.47081542e-01, -1.79328112e-02, -1.67287295e-03, -1.93751070e-02, 0.00000000e+00], [ 5.37137391e-02, -2.59418901e-01, -3.35850517e-01, 4.20931617e-01, 1.17590691e-01, 1.65422264e-03, 1.76535094e-01, 2.33379245e-02, -1.87082809e-01, -1.28016433e-02, 1.14854101e-01, 3.92577084e-02, -3.73712344e-02, -5.75259730e-02, -4.17452344e-02, 3.35038875e-02, 3.20095051e-02, -6.03835028e-01, 1.53909149e-01, 3.81368362e-01, 1.88648133e-02, -3.02161020e-04, -4.12720037e-02, 0.00000000e+00], [ -3.29272081e-02, 2.73455770e-01, -3.31833205e-01, -1.00278309e-01, 1.98147384e-01, -4.18105071e-02, -4.13318757e-02, 1.62985929e-01, -4.29569771e-02, -4.10103400e-02, -1.24235383e-01, 1.23016014e-02, 4.25628092e-02, -3.46191422e-02, -7.19288349e-02, -1.79453150e-02, 1.52906886e-01, 1.51852840e-01, -6.70229735e-01, 4.56695234e-01, -1.12316617e-02, 8.40942075e-04, -5.21690839e-02, 0.00000000e+00], [ -9.52388969e-03, 5.06361338e-02, 5.63077523e-02, -2.75038299e-02, 2.29377868e-03, -2.22927176e-03, -3.01627462e-02, -1.16180410e-02, 7.26877680e-03, -4.68706154e-03, 7.00835097e-01, -1.55779542e-01, -5.39777423e-04, -8.47090663e-02, -5.26455924e-01, 1.54634859e-01, -1.07415713e-02, 1.12691533e-01, -2.49437772e-02, 6.06917592e-04, 4.31703604e-03, 5.70128458e-03, -3.93851734e-01, 0.00000000e+00], [ -7.47531264e-03, 5.27229490e-02, -2.62943060e-03, 1.38653385e-02, 3.25091872e-02, -6.53558011e-03, -1.74325321e-02, 1.92192393e-03, -1.73892260e-02, -8.58267023e-03, 5.47574061e-01, -2.42327311e-01, 1.98678694e-01, -2.44500792e-01, 6.84787257e-01, -1.86242104e-02, 1.94956118e-02, 6.70149459e-02, -6.16480684e-02, 8.20344211e-02, -3.92623880e-03, 3.95966079e-03, 2.29467993e-01, 0.00000000e+00], [ -3.08190762e-04, -1.09237189e-04, 7.63461611e-03, 2.47482619e-02, 1.50596252e-02, 4.48223885e-04, -1.44205517e-03, -6.43194087e-03, -1.88340432e-02, -1.03163355e-02, -1.70167697e-02, 6.77440774e-02, 2.29976100e-01, 5.23745904e-01, 2.76184012e-01, 4.97008734e-01, -3.50913420e-02, 1.95960332e-02, 1.37974090e-03, 4.81302440e-02, 2.49159052e-02, 1.33031637e-02, -3.37537680e-01, 4.73879042e-01], [ -5.29923115e-03, 2.98045194e-02, 5.75729863e-02, 4.42660428e-02, 2.17533492e-02, -2.26555002e-03, -1.88820351e-02, -2.80204253e-02, -2.11249678e-02, -4.95504063e-03, 1.37781563e-01, 6.37843429e-01, 3.02493165e-01, -3.06232622e-02, -1.08532422e-01, 4.38415361e-01, -3.12736361e-02, 8.94905097e-02, 1.29148511e-02, 5.99677300e-02, 1.92633017e-02, 1.93534689e-03, 3.53419003e-01, -3.62538419e-01], [ -6.69678975e-03, 3.99755074e-02, 2.03416154e-02, -2.02766580e-02, 6.74238335e-03, -3.10430317e-03, -1.88668987e-02, -8.73358081e-04, 4.79013105e-03, -2.10129790e-03, 1.92690456e-01, 5.30658249e-01, -6.73296231e-01, -1.85540440e-01, 1.47823187e-01, -2.80536678e-02, 7.90341297e-03, 6.77377449e-02, -3.68989997e-02, 1.09044061e-02, -3.86038025e-03, 3.19199653e-04, -2.04000084e-03, 4.02272262e-01], [ 3.44016404e-03, -2.97622175e-02, 2.24635192e-02, -1.65146416e-02, -2.54129774e-02, 5.00644272e-03, 5.27701594e-03, -6.95841347e-03, 1.50585787e-02, 5.18835064e-03, -2.66725599e-01, 7.85506046e-02, 3.74645672e-01, -7.59281144e-01, -3.06562708e-02, 1.41407356e-01, -2.59402851e-02, -1.87365649e-02, 4.65463165e-02, -6.48722269e-02, 8.55963111e-03, -8.96414257e-04, -2.53726266e-01, 3.32202295e-01], [ -3.06493631e-03, 2.23544608e-02, 3.40913110e-02, 5.61908371e-02, 1.22625693e-02, -3.38185295e-03, -1.03684830e-02, -2.30733532e-02, -2.38365678e-02, -3.84450350e-03, 1.76846786e-01, 2.73456889e-01, 4.50031196e-01, 1.97469281e-01, -1.93699069e-01, -6.21928043e-01, 4.85708635e-02, 2.71809688e-02, 4.50957143e-02, 4.27767672e-02, -3.32320185e-02, 4.17854289e-03, 1.86599830e-01, 4.23077296e-01], [ -2.44277376e-03, 1.09185342e-02, 2.91809620e-02, 1.32301721e-02, -2.94210441e-03, 3.94976745e-05, -9.24834627e-03, -1.21019043e-02, -1.00694510e-02, -8.74592901e-03, 1.30078168e-02, 3.55520016e-01, 9.86585043e-02, 3.63597046e-02, 3.13252697e-01, -3.43533510e-01, 2.23627076e-02, 2.13305680e-02, 3.56848323e-02, 4.55686761e-03, -1.80944132e-02, 1.63984065e-02, -6.70069847e-01, -4.39122776e-01]])
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"]
True
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()
Solution is not as accurate as I would like...
basic_model_viz.residuals.plot(subplots=True, style=['r', 'b'])
plt.show()
...actually, when using noramlized residuals eveything looks great!
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()
<matplotlib.axes._subplots.AxesSubplot at 0x10af3fb50>
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
array([K_bar - m*(alpha/(R*a))**(1/(-alpha + 1)), a*(K_bar - m*(alpha/(R*a))**(1/(-alpha + 1)))/(R - 1), R*a/(R - 1), a], dtype=object)
credit_cycles(0, initial_condition)
array([ -8.64019967e-12, 0.00000000e+00, 0.00000000e+00, -9.88098492e-14])
jacobian(0, initial_condition)
array([[ 1.01000000e+02, -9.61904762e-01, 9.52163239e+02, -9.07029478e-01], [ 1.05000000e+02, -1.00000000e+00, 9.89872674e+02, 0.00000000e+00], [ 0.00000000e+00, 0.00000000e+00, 1.00000000e-02, -1.01000000e+00], [ 2.15879931e+00, 0.00000000e+00, 0.00000000e+00, -1.00000000e+00]])
from scipy import linalg
from IPython.html.widgets import fixed, interact, FloatSliderWidget
:0: FutureWarning: IPython widgets are experimental and may change in the future.
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))
[ 4.96886453e+04+86405.36943471j 4.96886453e+04-86405.36943471j -9.93761281e+04 +0.j 4.01197590e-08 +0.j ]
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
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))
[-0.96239413+0.j 0.03119783+0.31974317j 0.03119783-0.31974317j 0.05000000+0.j ] [ 0.96239413 0.32126158 0.32126158 0.05 ]
<function __main__.eigenvalues>
full_model_solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain)
full_model_solver.result["success"]
False
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)
array([<matplotlib.axes._subplots.AxesSubplot object at 0x112f80b90>, <matplotlib.axes._subplots.AxesSubplot object at 0x112ff11d0>, <matplotlib.axes._subplots.AxesSubplot object at 0x113065e10>, <matplotlib.axes._subplots.AxesSubplot object at 0x1130bec50>], dtype=object)
full_model_viz.normalized_residuals.plot(subplots=True)
plt.show()