``````

In [1]:

%matplotlib inline

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

In [2]:

import numpy as np
import seaborn as sn

import models
import orthogonal_collocation

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

In [15]:

# define some variables
t, X, i, pi = sym.symbols('t, X, i, pi')

# define some parameters
r = sym.symbols('r')

# assumes logarithmic utility from consumption
X_dot = i - pi - r

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

In [16]:

# define some parameters
epsilon, theta, rho, psi = sym.symbols('epsilon, theta, rho, psi')

# assumes convex adjustment costs a la Rotemberg
pi_dot = rho * pi - ((epsilon - 1) / theta) * (X**(1 + psi) - 1)

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

In [17]:

# define some parameters
istar, phi_pi, phi_X = sym.symbols('istar, phi_pi, phi_X')

# define some Taylor-esque rule for setting the nominal interest rate
taylor_rule = istar + phi_pi * pi + phi_X * sym.log(X)

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

In [14]:

rhs = {X: X_dot.subs({i: taylor_rule}), pi: pi_dot}

# define some boundary conditions
X0 = 5.0
pi0 = 1.0
bcs = {'lower': [X - X0, pi - pi0], 'upper': None}

basic_NK = models.BoundaryValueProblem(dependent_vars=[X, pi],
independent_var=t,
rhs=rhs,
boundary_conditions=bcs)

params = {'rho': 0.03, 'r': 0.04, 'epsilon': 1.04, 'theta': 0.5, 'psi': 0.5,
'istar': 1.0, 'phi_pi': 0.5, 'phi_X': 1.05}

basic_NK_solver = orthogonal_collocation.OrthogonalCollocationSolver(basic_NK, params)

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

In [17]:

# specify an initial guess
def kstar(g, n, alpha, delta, rho, theta, sigma):
"""Steady-state level of capital (per unit effective labor)."""
gamma = (sigma - 1) / sigma
if gamma == 0:
return (alpha / (delta + rho + theta * g))**(1 / (1 - alpha))
else:
return (1 - alpha)**(1 / gamma) * ((alpha / (delta + rho + theta * g))**(gamma / (gamma - 1)) - alpha)**(-1 / gamma)

def cstar(g, n, alpha, delta, rho, theta, sigma):
"""Steady-state level of consumption (per unit effective labor)."""
return kstar(g, n, alpha, delta, rho, theta, sigma)**alpha - (g + n + delta) * kstar(g, n, alpha, delta, rho, theta, sigma)

domain = [0, 200]
ts = np.linspace(domain[0], domain[1], 1000)
ks = kstar(**ramsey_params) - (kstar(**ramsey_params) - k0) * np.exp(-ts)
initial_capital_poly = np.polynomial.Chebyshev.fit(ts, ks, 50, domain)

cs = np.log(initial_capital_poly(ks))  # guess that consumption is concave function of capital
initial_consumption_poly = np.polynomial.Chebyshev.fit(ts, cs, 50, domain)
initial_ramsey_coefs = {k: initial_capital_poly.coef, c: initial_consumption_poly.coef}

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

In [18]:

ramsey_solution = ramsey_solver.solve(kind="Chebyshev",
coefs_dict=initial_ramsey_coefs,
domain=domain)

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

In [19]:

ramsey_solution.number_knots = 1000
ramsey_solution.solution.plot(subplots=True, style=['r', 'b'])

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

Out[19]:

array([<matplotlib.axes._subplots.AxesSubplot object at 0x1172bf790>,
<matplotlib.axes._subplots.AxesSubplot object at 0x1192651d0>], dtype=object)

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

In [20]:

ramsey_solution.residuals.plot(subplots=True, style=['r', 'b'])

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

Out[20]:

array([<matplotlib.axes._subplots.AxesSubplot object at 0x1193876d0>,
<matplotlib.axes._subplots.AxesSubplot object at 0x1194e4ed0>], dtype=object)

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

In [ ]:

``````