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