|
1
|
|
|
from nose import tools |
|
2
|
|
|
|
|
3
|
|
|
import numpy as np |
|
4
|
|
|
from scipy import stats |
|
5
|
|
|
|
|
6
|
|
|
from . import models |
|
7
|
|
|
from .. import basis_functions |
|
8
|
|
|
from .. import solvers |
|
9
|
|
|
|
|
10
|
|
|
|
|
11
|
|
|
def analytic_solution(t, k0, alpha, delta, g, n, s, **params): |
|
12
|
|
|
"""Analytic solution for model with Cobb-Douglas production.""" |
|
13
|
|
|
lmbda = (g + n + delta) * (1 - alpha) |
|
14
|
|
|
ks = (((s / (g + n + delta)) * (1 - np.exp(-lmbda * t)) + |
|
15
|
|
|
k0**(1 - alpha) * np.exp(-lmbda * t))**(1 / (1 - alpha))) |
|
16
|
|
|
return ks |
|
17
|
|
|
|
|
18
|
|
|
|
|
19
|
|
|
def cobb_douglas_output(k, alpha, **params): |
|
20
|
|
|
"""Intensive output has Cobb-Douglas functional form.""" |
|
21
|
|
|
return k**alpha |
|
22
|
|
|
|
|
23
|
|
|
|
|
24
|
|
|
def equilibrium_capital(alpha, delta, g, n, s, **params): |
|
25
|
|
|
"""Equilibrium value of capital (per unit effective labor supply).""" |
|
26
|
|
|
return (s / (g + n + delta))**(1 / (1 - alpha)) |
|
27
|
|
|
|
|
28
|
|
|
|
|
29
|
|
|
def generate_random_params(scale, seed): |
|
30
|
|
|
np.random.seed(seed) |
|
31
|
|
|
|
|
32
|
|
|
# random g, n, delta such that sum of these params is positive |
|
33
|
|
|
g, n = stats.norm.rvs(0.05, scale, size=2) |
|
34
|
|
|
delta, = stats.lognorm.rvs(scale, loc=g + n, size=1) |
|
35
|
|
|
assert g + n + delta > 0 |
|
36
|
|
|
|
|
37
|
|
|
# s and alpha must be on (0, 1) (but lower values are more reasonable) |
|
38
|
|
|
s, alpha = stats.beta.rvs(a=1, b=3, size=2) |
|
39
|
|
|
|
|
40
|
|
|
# choose k0 so that it is not too far from equilibrium |
|
41
|
|
|
kstar = equilibrium_capital(alpha, delta, g, n, s) |
|
42
|
|
|
k0, = stats.uniform.rvs(0.5 * kstar, 1.5 * kstar, size=1) |
|
43
|
|
|
assert k0 > 0 |
|
44
|
|
|
|
|
45
|
|
|
params = {'g': g, 'n': n, 'delta': delta, 's': s, 'alpha': alpha, |
|
46
|
|
|
'k0': k0} |
|
47
|
|
|
|
|
48
|
|
|
return params |
|
49
|
|
|
|
|
50
|
|
|
|
|
51
|
|
|
def initial_mesh(t, T, num, problem): |
|
52
|
|
|
ts = np.linspace(t, T, num) |
|
53
|
|
|
kstar = equilibrium_capital(**problem.params) |
|
54
|
|
|
ks = kstar - (kstar - problem.params['k0']) * np.exp(-ts) |
|
55
|
|
|
return ts, ks |
|
56
|
|
|
|
|
57
|
|
|
|
|
58
|
|
|
random_seed = np.random.randint(2147483647) |
|
59
|
|
|
random_params = generate_random_params(0.1, random_seed) |
|
60
|
|
|
test_problem = models.SolowModel(cobb_douglas_output, equilibrium_capital, |
|
61
|
|
|
random_params) |
|
62
|
|
|
|
|
63
|
|
|
polynomial_basis = basis_functions.PolynomialBasis() |
|
64
|
|
|
solver = solvers.Solver(polynomial_basis) |
|
65
|
|
|
|
|
66
|
|
|
|
|
67
|
|
|
def _test_polynomial_collocation(basis_kwargs, boundary_points, num=1000): |
|
68
|
|
|
"""Helper function for testing various kinds of polynomial collocation.""" |
|
69
|
|
|
ts, ks = initial_mesh(*boundary_points, num=num, problem=test_problem) |
|
70
|
|
|
k_poly = polynomial_basis.fit(ts, ks, **basis_kwargs) |
|
71
|
|
|
initial_coefs = k_poly.coef |
|
72
|
|
|
nodes = polynomial_basis.roots(**basis_kwargs) |
|
73
|
|
|
|
|
74
|
|
|
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs, |
|
75
|
|
|
nodes, test_problem) |
|
76
|
|
|
|
|
77
|
|
|
# check that solver terminated successfully |
|
78
|
|
|
msg = "Solver failed!\nSeed: {}\nModel params: {}\n" |
|
79
|
|
|
tools.assert_true(solution.result.success, |
|
80
|
|
|
msg=msg.format(random_seed, test_problem.params)) |
|
81
|
|
|
|
|
82
|
|
|
# compute the residuals |
|
83
|
|
|
normed_residuals = solution.normalize_residuals(ts) |
|
84
|
|
|
|
|
85
|
|
|
# check that residuals are close to zero on average |
|
86
|
|
|
tools.assert_true(np.mean(normed_residuals) < 1e-6, |
|
87
|
|
|
msg=msg.format(random_seed, test_problem.params)) |
|
88
|
|
|
|
|
89
|
|
|
# check that the numerical and analytic solutions are close |
|
90
|
|
|
numeric_soln = solution.evaluate_solution(ts) |
|
91
|
|
|
analytic_soln = analytic_solution(ts, **test_problem.params) |
|
92
|
|
|
tools.assert_true(np.mean(numeric_soln - analytic_soln) < 1e-6) |
|
93
|
|
|
|
|
94
|
|
|
|
|
95
|
|
|
def test_chebyshev_collocation(): |
|
96
|
|
|
"""Test collocation solver using Chebyshev polynomials for basis.""" |
|
97
|
|
|
boundary_points = (0, 100) |
|
98
|
|
|
basis_kwargs = {'kind': 'Chebyshev', 'degree': 50, 'domain': boundary_points} |
|
99
|
|
|
_test_polynomial_collocation(basis_kwargs, boundary_points) |
|
100
|
|
|
|
|
101
|
|
|
|
|
102
|
|
|
def test_legendre_collocation(): |
|
103
|
|
|
"""Test collocation solver using Legendre polynomials for basis.""" |
|
104
|
|
|
boundary_points = (0, 100) |
|
105
|
|
|
basis_kwargs = {'kind': 'Legendre', 'degree': 50, 'domain': boundary_points} |
|
106
|
|
|
_test_polynomial_collocation(basis_kwargs, boundary_points) |
|
107
|
|
|
|
davidrpugh /
pyCollocation
