44"""
55
66import numpy as np
7- import quantecon as qe
87import time
98
10- from scipy .spatial import ConvexHull
11- from scipy .optimize import linprog , minimize , minimize_scalar
12- from scipy .interpolate import UnivariateSpline
9+ from scipy .optimize import linprog , minimize
1310import numpy .polynomial .chebyshev as cheb
1411
1512
@@ -30,7 +27,7 @@ def __init__(self, β, mbar, h_min, h_max, n_h, n_m, N_g):
3027 self .N_a = self .n_h * self .n_m
3128
3229 # Utility and production functions
33- uc = lambda c : np .log (c )
30+ uc = lambda c : np .log (np . maximum ( c , 1e-10 )) # Clip to 1e-10 to avoid log(0) or log(-ve )
3431 uc_p = lambda c : 1 / c
3532 v = lambda m : 1 / 500 * (mbar * m - 0.5 * m ** 2 )** 0.5
3633 v_p = lambda m : 0.5 / 500 * (mbar * m - 0.5 * m ** 2 )** (- 0.5 ) * (mbar - m )
@@ -306,7 +303,7 @@ def solve_bellman(self, θ_min, θ_max, order, disp=False, tol=1e-7, maxiters=10
306303 mbar = self .mbar
307304
308305 # Utility and production functions
309- uc = lambda c : np .log (c )
306+ uc = lambda c : np .log (np . maximum ( c , 1e-10 )) # Clip to 1e-10 to avoid log(0) or log(-ve )
310307 uc_p = lambda c : 1 / c
311308 v = lambda m : 1 / 500 * (mbar * m - 0.5 * m ** 2 )** 0.5
312309 v_p = lambda m : 0.5 / 500 * (mbar * m - 0.5 * m ** 2 )** (- 0.5 ) * (mbar - m )
@@ -343,13 +340,13 @@ def p_fun(x):
343340 scale = - 1 + 2 * (x [2 ] - θ_min )/ (θ_max - θ_min )
344341 p_fun = - (u (x [0 ], x [1 ]) \
345342 + self .β * np .dot (cheb .chebvander (scale , order - 1 ), c ))
346- return p_fun
343+ return p_fun . item ()
347344
348345 def p_fun2 (x ):
349346 scale = - 1 + 2 * (x [1 ] - θ_min )/ (θ_max - θ_min )
350347 p_fun = - (u (x [0 ],mbar ) \
351348 + self .β * np .dot (cheb .chebvander (scale , order - 1 ), c ))
352- return p_fun
349+ return p_fun . item ()
353350
354351 cons1 = ({'type' : 'eq' , 'fun' : lambda x : uc_p (f (x [0 ], x [1 ])) * x [1 ]
355352 * (x [0 ] - 1 ) + v_p (x [1 ]) * x [1 ] + self .β * x [2 ] - θ },
@@ -416,8 +413,9 @@ def p_fun2(x):
416413 h_grid = np .zeros (100 )
417414 for i in range (100 ):
418415 θ = θ_grid_fine [i ]
416+ x0 = np .clip (lb1 + (ub1 - lb1 )/ 2 , lb1 , ub1 )
419417 res = minimize (p_fun ,
420- lb1 + ( ub1 - lb1 ) / 2 ,
418+ x0 ,
421419 method = 'SLSQP' ,
422420 bounds = bnds1 ,
423421 constraints = cons1 ,
@@ -428,8 +426,9 @@ def p_fun2(x):
428426 θ_prime_grid [i ] = res .x [2 ]
429427 h_grid [i ] = res .x [0 ]
430428 m_grid [i ] = res .x [1 ]
429+ x0 = np .clip (lb2 + (ub2 - lb2 )/ 2 , lb2 , ub2 )
431430 res = minimize (p_fun2 ,
432- lb2 + ( ub2 - lb2 ) / 2 ,
431+ x0 ,
433432 method = 'SLSQP' ,
434433 bounds = bnds2 ,
435434 constraints = cons2 ,
@@ -441,7 +440,8 @@ def p_fun2(x):
441440 h_grid [i ] = res .x [0 ]
442441 m_grid [i ] = self .mbar
443442 scale = - 1 + 2 * (θ - θ_min )/ (θ_max - θ_min )
444- resid_grid [i ] = np .dot (cheb .chebvander (scale , order - 1 ), c ) - p
443+ resid_grid_val = np .dot (cheb .chebvander (scale , order - 1 ), c ) - p
444+ resid_grid [i ] = resid_grid_val .item ()
445445
446446 self .resid_grid = resid_grid
447447 self .θ_grid_fine = θ_grid_fine
@@ -465,13 +465,14 @@ def ValFun(x):
465465 res = minimize (ValFun ,
466466 (θ_min + θ_max )/ 2 ,
467467 bounds = [(θ_min , θ_max )])
468- θ_series [0 ] = res .x
468+ θ_series [0 ] = res .x . item ()
469469
470470 # Simulate
471471 for i in range (30 ):
472472 θ = θ_series [i ]
473+ x0 = np .clip (lb1 + (ub1 - lb1 )/ 2 , lb1 , ub1 )
473474 res = minimize (p_fun ,
474- lb1 + ( ub1 - lb1 ) / 2 ,
475+ x0 ,
475476 method = 'SLSQP' ,
476477 bounds = bnds1 ,
477478 constraints = cons1 ,
@@ -481,8 +482,9 @@ def ValFun(x):
481482 h_series [i ] = res .x [0 ]
482483 m_series [i ] = res .x [1 ]
483484 θ_series [i + 1 ] = res .x [2 ]
485+ x0 = np .clip (lb2 + (ub2 - lb2 )/ 2 , lb2 , ub2 )
484486 res2 = minimize (p_fun2 ,
485- lb2 + ( ub2 - lb2 ) / 2 ,
487+ x0 ,
486488 method = 'SLSQP' ,
487489 bounds = bnds2 ,
488490 constraints = cons2 ,
0 commit comments