``````

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

``````

## The basic model: amplification and persistance

``````

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]:

98.97176708265565

``````
``````

In [85]:

B0

``````
``````

Out[85]:

142.5

``````
``````

In [86]:

Bstar

``````
``````

Out[86]:

199.9229695069644

``````
``````

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]:

[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 ])]

``````
``````

In [100]:

initial_capital_poly.coef

``````
``````

Out[100]:

array([ 98.25192341,   1.29704214,  -0.91771876,   0.5470026 ,
-0.25019508,   0.10814735])

``````
``````

In [97]:

solution.result

``````
``````

Out[97]:

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]])

``````
``````

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]:

True

``````
``````

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]:

<matplotlib.axes._subplots.AxesSubplot at 0x10af3fb50>

``````
``````

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]:

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)

``````
``````

In [901]:

credit_cycles(0, initial_condition)

``````
``````

Out[901]:

array([ -8.64019967e-12,   0.00000000e+00,   0.00000000e+00,
-9.88098492e-14])

``````
``````

In [868]:

jacobian(0, initial_condition)

``````
``````

Out[868]:

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]])

``````
``````

In [98]:

from scipy import linalg

``````
``````

In [99]:

from IPython.html.widgets import fixed, interact, FloatSliderWidget

``````
``````

:0: FutureWarning: IPython widgets are experimental and may change in the future.

``````
``````

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))

``````
``````

[  4.96886453e+04+86405.36943471j   4.96886453e+04-86405.36943471j
-9.93761281e+04    +0.j           4.01197590e-08    +0.j        ]

``````
``````

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

``````

## Full model

``````

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))

``````
``````

[-0.96239413+0.j          0.03119783+0.31974317j  0.03119783-0.31974317j
0.05000000+0.j        ]
[ 0.96239413  0.32126158  0.32126158  0.05      ]

Out[120]:

<function __main__.eigenvalues>

``````
``````

In [ ]:

``````
``````

In [126]:

full_model_solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain)

``````
``````

In [127]:

full_model_solver.result["success"]

``````
``````

Out[127]:

False

``````
``````

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]:

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)

``````
``````

In [52]:

full_model_viz.normalized_residuals.plot(subplots=True)
plt.show()

``````
``````

``````
``````

In [ ]:

``````