I want to use scipy.optimize.root to solve a static equilibrium model. Specifically, I need find the root of a three-dimensional function system. The code is given as below:
# Import package
import numpy as np
from scipy.optimize import root,fsolve
# Parameters
Kbar = 10 # Total capital
Lbar = 20 # Total labour
α = 0.3
β = np.array([0.3,0.6])
p = np.array([1.0,0.1])
# The first element of the price array is p1=1 (constant), the second is p2 which needs to be optimized.
# Defining the objective function system
def markets(x):
# prices (excluded p1=1), x is a 3-dimentional array
p[1] = x[0] # p2 is the first element of input x
w = x[1] # wage w is the second element of input x
r = x[2] # interest rate r is the third elemendt of input x
# total income
Ybar = w * Lbar + r * Kbar
# get market equation
sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0])
sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1])
sol[2] = β[0]*α*Ybar/w + β[1]*(1-α)*Ybar/w - Lbar
return sol
# Initial guess x0 is a 3*1 array
x0 = np.zeros(3) + 5
# Find market equilibrium
res = root(markets, x0)
print(res)
However, the results report an error:
fjac: array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
fun: array([nan, nan, nan])
message: 'The iteration is not making good progress, as measured by the \n improvement from the last ten iterations.'
nfev: 17
qtf: array([ 0.89142371, 0.09796604, -4.69999992])
r: array([nan, nan, nan, nan, nan, nan])
status: 5
success: False
x: array([5., 5., 5.])
<ipython-input-37-15479e9f467a>:24: RuntimeWarning: invalid value encountered in double_scalars
sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0])
<ipython-input-37-15479e9f467a>:25: RuntimeWarning: invalid value encountered in double_scalars
sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1])
I searched for some information about the "RuntimeWarning: invalid value encountered in double_scalars" and know that this error would show when dividing by zero. However, I did not find any possible bug in my function.
Can anyone help me out of this? Many thanks.
First, this isn't a minimal working example (MWE) since you miss defining sol inside your markets function. So let's assume
# Defining the objective function system
def markets(x):
# prices (excluded p1=1), x is a 3-dimentional array
p[1] = x[0] # p2 is the first element of input x
w = x[1] # wage w is the second element of input x
r = x[2] # interest rate r is the third elemendt of input x
# total income
Ybar = w * Lbar + r * Kbar
# get market equation
sol = np.zeros(3)
sol[0] = 1/p[0]-(β[0]/w)**β[0]*((1-β[0])/r)**(1-β[0])
sol[1] = 1/p[1]-(β[1]/w)**β[1]*((1-β[1])/r)**(1-β[1])
sol[2] = β[0]*α*Ybar/w + β[1]*(1-α)*Ybar/w - Lbar
return sol
Then, as you already mentioned, the error is due to division by zero in the term 1/p[1]. Therefore, you need to find a root x > 0, i.e. we seek for a point x > 0 such that markets(x) = 0. However, the scipy.optimize.root method doesn't support simple variable bounds.
By writing the root problem as equivalent constrained minimization problem:
min ||markets(x)|| s.t. x >= eps
where eps > 0, you can solve the problem with the help of scipy.optimize.minimize:
from scipy.optimize import minimize
# Variable bounds
eps = 1.0e-4
bnds = [(eps, None) for _ in range(3)]
# initial point
x0 = np.zeros(3) + 5
# Solve the minimization problem
res = minimize(lambda x: np.linalg.norm(markets(x)), x0=x0, bounds=bnds)