|
1
|
|
|
""" |
|
2
|
|
|
Test case using the Ramsey-Cass-Koopmans model of optimal savings. |
|
3
|
|
|
|
|
4
|
|
|
Model assumes Cobb-Douglas production technology and Constant Relative Risk |
|
5
|
|
|
Aversion (CRRA) preferences. I then generate random parameter values |
|
6
|
|
|
consistent with a constant consumption-output ratio (i.e., constant savings |
|
7
|
|
|
rate). |
|
8
|
|
|
|
|
9
|
|
|
""" |
|
10
|
|
|
from nose import tools |
|
11
|
|
|
|
|
12
|
|
|
import numpy as np |
|
13
|
|
|
from scipy import stats |
|
14
|
|
|
|
|
15
|
|
|
from . import models |
|
16
|
|
|
from .. import basis_functions |
|
17
|
|
|
from .. import solvers |
|
18
|
|
|
|
|
19
|
|
|
|
|
20
|
|
|
def analytic_solution(t, A0, alpha, delta, g, K0, n, N0, theta, **params): |
|
21
|
|
|
"""Analytic solution for model with Cobb-Douglas production.""" |
|
22
|
|
|
k0 = K0 / (A0 * N0) |
|
23
|
|
|
lmbda = (g + n + delta) * (1 - alpha) |
|
24
|
|
|
ks = ((k0**(1 - alpha) * np.exp(-lmbda * t) + |
|
25
|
|
|
(1 / (theta * (g + n + delta))) * (1 - np.exp(-lmbda * t)))**(1 / (1 - alpha))) |
|
26
|
|
|
ys = cobb_douglas_output(ks, alpha, **params) |
|
27
|
|
|
cs = ((theta - 1) / theta) * ys |
|
28
|
|
|
|
|
29
|
|
|
return [ks, cs] |
|
30
|
|
|
|
|
31
|
|
|
|
|
32
|
|
|
def cobb_douglas_output(k, alpha, **params): |
|
33
|
|
|
"""Cobb-Douglas output (per unit effective labor supply).""" |
|
34
|
|
|
return k**alpha |
|
35
|
|
|
|
|
36
|
|
|
|
|
37
|
|
|
def cobb_douglas_mpk(k, alpha, **params): |
|
38
|
|
|
"""Marginal product of capital with Cobb-Douglas production.""" |
|
39
|
|
|
return alpha * k**(alpha - 1) |
|
40
|
|
|
|
|
41
|
|
|
|
|
42
|
|
|
def equilibrium_capital(alpha, delta, g, n, rho, theta, **params): |
|
43
|
|
|
"""Steady state value for capital stock (per unit effective labor).""" |
|
44
|
|
|
return (alpha / (delta + rho + theta * g))**(1 / (1 - alpha)) |
|
45
|
|
|
|
|
46
|
|
|
|
|
47
|
|
|
def generate_random_params(scale, seed): |
|
48
|
|
|
np.random.seed(seed) |
|
49
|
|
|
|
|
50
|
|
|
# random g, n, delta such that sum of these params is positive |
|
51
|
|
|
g, n = stats.norm.rvs(0.05, scale, size=2) |
|
52
|
|
|
delta, = stats.lognorm.rvs(scale, loc=g + n, size=1) |
|
53
|
|
|
assert g + n + delta > 0 |
|
54
|
|
|
|
|
55
|
|
|
# choose alpha consistent with theta > 1 |
|
56
|
|
|
upper_bound = 1 if (g + delta) < 0 else min(1, (g + delta) / (g + n + delta)) |
|
57
|
|
|
alpha, = stats.uniform.rvs(scale=upper_bound, size=1) |
|
58
|
|
|
assert 0 < alpha < upper_bound |
|
59
|
|
|
|
|
60
|
|
|
lower_bound = delta / (alpha * (g + n + delta) - g) |
|
61
|
|
|
theta, = stats.lognorm.rvs(scale, loc=lower_bound, size=1) |
|
62
|
|
|
rho = alpha * theta * (g + n + delta) - (delta * theta * g) |
|
63
|
|
|
assert rho > 0 |
|
64
|
|
|
|
|
65
|
|
|
# normalize A0 and N0 |
|
66
|
|
|
A0 = 1.0 |
|
67
|
|
|
N0 = A0 |
|
68
|
|
|
|
|
69
|
|
|
# choose K0 so that it is not too far from balanced growth path |
|
70
|
|
|
kstar = equilibrium_capital(g, n, alpha, delta, rho, theta) |
|
71
|
|
|
K0, = stats.uniform.rvs(0.5 * kstar, 1.5 * kstar, size=1) |
|
72
|
|
|
assert K0 > 0 |
|
73
|
|
|
|
|
74
|
|
|
params = {'g': g, 'n': n, 'delta': delta, 'rho': rho, 'alpha': alpha, |
|
75
|
|
|
'theta': theta, 'K0': K0, 'A0': A0, 'N0': N0} |
|
76
|
|
|
|
|
77
|
|
|
return params |
|
78
|
|
|
|
|
79
|
|
|
|
|
80
|
|
|
def initial_mesh(t, T, num, problem): |
|
81
|
|
|
# compute equilibrium values |
|
82
|
|
|
cstar = problem.equilibrium_consumption(**problem.params) |
|
83
|
|
|
kstar = problem.equilibrium_capital(**problem.params) |
|
84
|
|
|
ystar = problem.intensive_output(kstar, **problem.params) |
|
85
|
|
|
|
|
86
|
|
|
# create the mesh for capital |
|
87
|
|
|
ts = np.linspace(t, T, num) |
|
88
|
|
|
k0 = problem.params['K0'] / (problem.params['A0'] * problem.params['N0']) |
|
89
|
|
|
ks = kstar - (kstar - k0) * np.exp(-ts) |
|
90
|
|
|
|
|
91
|
|
|
# create the mesh for consumption |
|
92
|
|
|
s = 1 - (cstar / ystar) |
|
93
|
|
|
y0 = cobb_douglas_output(k0, **problem.params) |
|
94
|
|
|
c0 = (1 - s) * y0 |
|
95
|
|
|
cs = cstar - (cstar - c0) * np.exp(-ts) |
|
96
|
|
|
|
|
97
|
|
|
return ts, ks, cs |
|
98
|
|
|
|
|
99
|
|
|
|
|
100
|
|
|
def pratt_arrow_risk_aversion(t, c, theta, **params): |
|
101
|
|
|
"""Assume constant relative risk aversion""" |
|
102
|
|
|
return theta / c |
|
103
|
|
|
|
|
104
|
|
|
|
|
105
|
|
|
random_seed = np.random.randint(2147483647) |
|
106
|
|
|
random_params = generate_random_params(0.1, random_seed) |
|
107
|
|
|
test_problem = models.RamseyCassKoopmansModel(pratt_arrow_risk_aversion, |
|
108
|
|
|
cobb_douglas_output, |
|
109
|
|
|
equilibrium_capital, |
|
110
|
|
|
cobb_douglas_mpk, |
|
111
|
|
|
random_params) |
|
112
|
|
|
|
|
113
|
|
|
polynomial_basis = basis_functions.PolynomialBasis() |
|
114
|
|
|
solver = solvers.Solver(polynomial_basis) |
|
115
|
|
|
|
|
116
|
|
|
|
|
117
|
|
|
def _test_polynomial_collocation(basis_kwargs, boundary_points, num=1000): |
|
118
|
|
|
"""Helper function for testing various kinds of polynomial collocation.""" |
|
119
|
|
|
ts, ks, cs = initial_mesh(*boundary_points, num=num, problem=test_problem) |
|
120
|
|
|
k_poly = polynomial_basis.fit(ts, ks, **basis_kwargs) |
|
121
|
|
|
c_poly = polynomial_basis.fit(ts, cs, **basis_kwargs) |
|
122
|
|
|
initial_coefs = np.hstack([k_poly.coef, c_poly.coef]) |
|
123
|
|
|
nodes = polynomial_basis.roots(**basis_kwargs) |
|
124
|
|
|
|
|
125
|
|
|
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs, |
|
126
|
|
|
nodes, test_problem) |
|
127
|
|
|
|
|
128
|
|
|
# check that solver terminated successfully |
|
129
|
|
|
msg = "Solver failed!\nSeed: {}\nModel params: {}\n".format(random_seed, test_problem.params) |
|
130
|
|
|
tools.assert_true(solution.result.success, msg=msg) |
|
131
|
|
|
|
|
132
|
|
|
# compute the residuals |
|
133
|
|
|
normed_residuals = solution.normalize_residuals(ts) |
|
134
|
|
|
|
|
135
|
|
|
# check that residuals are close to zero on average |
|
136
|
|
|
for normed_residual in normed_residuals: |
|
137
|
|
|
tools.assert_true(np.mean(normed_residual) < 1e-6, msg=msg) |
|
138
|
|
|
|
|
139
|
|
|
# check that the numerical and analytic solutions are close |
|
140
|
|
|
numeric_solns = solution.evaluate_solution(ts) |
|
141
|
|
|
analytic_solns = analytic_solution(ts, **test_problem.params) |
|
142
|
|
|
|
|
143
|
|
|
for numeric_soln, analytic_soln in zip(numeric_solns, analytic_solns): |
|
144
|
|
|
error = numeric_soln - analytic_soln |
|
145
|
|
|
tools.assert_true(np.mean(error) < 1e-6, msg=msg) |
|
146
|
|
|
|
|
147
|
|
|
|
|
148
|
|
|
def test_chebyshev_collocation(): |
|
149
|
|
|
"""Test collocation solver using Chebyshev polynomials for basis.""" |
|
150
|
|
|
boundary_points = (0, 100) |
|
151
|
|
|
basis_kwargs = {'kind': 'Chebyshev', 'degree': 50, 'domain': boundary_points} |
|
152
|
|
|
_test_polynomial_collocation(basis_kwargs, boundary_points) |
|
153
|
|
|
|
|
154
|
|
|
|
|
155
|
|
|
def test_legendre_collocation(): |
|
156
|
|
|
"""Test collocation solver using Legendre polynomials for basis.""" |
|
157
|
|
|
boundary_points = (0, 100) |
|
158
|
|
|
basis_kwargs = {'kind': 'Legendre', 'degree': 50, 'domain': boundary_points} |
|
159
|
|
|
_test_polynomial_collocation(basis_kwargs, boundary_points) |
|
160
|
|
|
|
davidrpugh /
pyCollocation
