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