From 7bd97a065b6ae8c2a4945e4bb381790417df4959 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sat, 4 Jul 2026 22:58:20 +0900 Subject: [PATCH 1/2] TST: Keep next states within the interpolation domain in test and benchmark models Co-Authored-By: Claude Fable 5 --- benchmark/benchmarks.jl | 23 +++++++++++++++-------- test/test_cdp_multidim.jl | 22 ++++++++++++++++------ 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index fda40b1..c8ed5d8 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -27,10 +27,14 @@ using QuantEcon: qnwlogn, qnwnorm #= Model definitions =# # 1-D stochastic optimal growth -function growth_model_1d(basis) +# The next state is kept within the interpolation domain [s_min, s_max]: +# feasible actions can otherwise map s' outside (even negative), where the +# fitted value function is extrapolated -- catastrophically so for the +# Chebyshev basis, making VFI diverge. +function growth_model_1d(basis, s_min, s_max) alpha = 0.65 f(s, x) = log(x) - g(s, x, e) = e * s^alpha - x + g(s, x, e) = clamp(e * s^alpha - x, s_min, s_max) shocks, weights = qnwlogn(9, 0.0, 0.01) x_lb(s) = 1e-8 x_ub(s) = s @@ -38,8 +42,9 @@ function growth_model_1d(basis) end # 2-D stochastic optimal growth with leisure (Santos, 1999, Sec. 7.3; -# same model as in test/test_cdp_multidim.jl) -function growth_model_2d(basis) +# same model as in test/test_cdp_multidim.jl). As in the 1-D model, the +# capital stock is kept within the interpolation domain [k_min, k_max]. +function growth_model_2d(basis, k_min, k_max) beta = 0.95 lambda = 1 / 3 A = 10.0 @@ -68,7 +73,7 @@ function growth_model_2d(basis) function g(s, x, e) k, logz = s z = exp(logz) - kp = kprime_from_x(k, z, x) + kp = clamp(kprime_from_x(k, z, x), k_min, k_max) logzp = rho * logz + e return (kp, logzp) end @@ -87,12 +92,14 @@ logz_min, logz_max = -0.32, 0.32 nk, nlogz = 43, 3 cases = [ - ("1d_cheb", growth_model_1d(Basis(ChebParams(50, s_min_1d, s_max_1d)))), + ("1d_cheb", growth_model_1d( + Basis(ChebParams(50, s_min_1d, s_max_1d)), s_min_1d, s_max_1d)), ("1d_spline", growth_model_1d( - Basis(SplineParams(99, s_min_1d, s_max_1d, 3)))), + Basis(SplineParams(99, s_min_1d, s_max_1d, 3)), s_min_1d, s_max_1d)), ("2d_spline", growth_model_2d( Basis(SplineParams(nk - 1, k_min, k_max, 2), - SplineParams(nlogz - 1, logz_min, logz_max, 2)))), + SplineParams(nlogz - 1, logz_min, logz_max, 2)), + k_min, k_max)), ] eval_grids = Dict( diff --git a/test/test_cdp_multidim.jl b/test/test_cdp_multidim.jl index ba6b825..b8c5375 100644 --- a/test/test_cdp_multidim.jl +++ b/test/test_cdp_multidim.jl @@ -59,10 +59,17 @@ end # Transition function + # During the inner maximization, candidate actions can map the capital + # stock outside the interpolation domain [k_min, k_max] (e.g. x -> 0 + # gives k' ~ z*A*k^alpha, up to ~30), where the fitted value function is + # extrapolated. For global polynomial (Chebyshev) bases such + # extrapolation is catastrophically large, so keep the next state within + # the domain. The true optimal policy maps the domain well inside + # [k_min, k_max], so the clamp does not affect the analytical checks. function g(s, x, e) k, logz = s z = exp(logz) - kp = kprime_from_x(k, z, x) + kp = clamp(kprime_from_x(k, z, x), k_min, k_max) logzp = rho * logz + e return (kp, logzp) end @@ -109,13 +116,16 @@ ), ( "Spline", - # Spline degree: quadratic over k, linear over logz - # Santos (1999) footnote 17: shape-preserving spline over k, - # cubic spline over z + # Spline degree: quadratic over k and over logz (with only 3 nodes + # over logz, cubic is infeasible). Santos (1999) footnote 17 uses + # shape-preserving spline over k and cubic spline over z, so the + # Table 20 policy error (1.92e-4) need not be attained here; the + # collocation solution's policy error is ~1.1e-2, hence the safety + # factor 100. Basis(SplineParams(breaks_k, k_min, k_max, dk), SplineParams(breaks_z, logz_min, logz_max, dz)), - 1.92e-2, # policy_tol: Santos (1999) Table 20 - 9.58e-1, # value_tol: Santos (1999) Table 20 + 1.92e-4 * 100, # policy_tol: Santos (1999) Table 20 with safety factor 100 + 9.58e-1, # value_tol: Santos (1999) Table 20 ), ( "Chebyshev", From 1ac757e570ca73d79397835668d7e8b7eceebab5 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Sat, 4 Jul 2026 23:18:04 +0900 Subject: [PATCH 2/2] Address Copilot review: restrict actions instead of clamping in 1-D benchmark model Tighten x_lb/x_ub using the quadrature shock extrema so that the next state stays within the interpolation domain for all shocks, keeping the transition law smooth (no kink from clamping). The 2-D Santos model keeps the clamp: its domain-respecting action bound has no closed form. Co-Authored-By: Claude Fable 5 --- benchmark/benchmarks.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index c8ed5d8..de049ad 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -27,23 +27,28 @@ using QuantEcon: qnwlogn, qnwnorm #= Model definitions =# # 1-D stochastic optimal growth -# The next state is kept within the interpolation domain [s_min, s_max]: -# feasible actions can otherwise map s' outside (even negative), where the -# fitted value function is extrapolated -- catastrophically so for the -# Chebyshev basis, making VFI diverge. +# The action bounds are tightened so that s' stays within the interpolation +# domain [s_min, s_max] for all quadrature shock nodes: feasible actions can +# otherwise map s' outside (even negative), where the fitted value function +# is extrapolated -- catastrophically so for the Chebyshev basis, making VFI +# diverge. Restricting the action set (rather than clamping s' inside `g`) +# keeps the transition law smooth. function growth_model_1d(basis, s_min, s_max) alpha = 0.65 f(s, x) = log(x) - g(s, x, e) = clamp(e * s^alpha - x, s_min, s_max) + g(s, x, e) = e * s^alpha - x shocks, weights = qnwlogn(9, 0.0, 0.01) - x_lb(s) = 1e-8 - x_ub(s) = s + e_min, e_max = extrema(shocks) + x_lb(s) = max(1e-8, e_max * s^alpha - s_max) + x_ub(s) = min(s, e_min * s^alpha - s_min) return ContinuousDP(f, g, 0.95, shocks, weights, x_lb, x_ub, basis) end # 2-D stochastic optimal growth with leisure (Santos, 1999, Sec. 7.3; -# same model as in test/test_cdp_multidim.jl). As in the 1-D model, the -# capital stock is kept within the interpolation domain [k_min, k_max]. +# same model as in test/test_cdp_multidim.jl). The capital stock is kept +# within the interpolation domain [k_min, k_max] by clamping in `g`: unlike +# in the 1-D model, the action bound ensuring k' <= k_max has no closed +# form (k' is transcendental in x). function growth_model_2d(basis, k_min, k_max) beta = 0.95 lambda = 1 / 3