I have a system of non-linear equations. I do not have a good guess about initial values. And I want at least one set of all positive roots cause in economics, negative values to these variables do not make much sense.
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 15 21:48:56 2016
@author: Nick
"""
import scipy as sp
from scipy.optimize import root, fsolve
import numpy as np
#from scipy.optimize import *
el = 1.1
eg = el
ej = 10
om = 0.3
omg = 0.3
rhog = 0.8
xi = 0.9
mun = 2
pidss = 0.02
muc = 0.001
ec = 2.00 # sims obtains 2.47
beta = 0.998
h = 0.8
kappa = 4.00
n = 1/3.0
alpha = 1/3.0
delta = 0.025
egs = eg
oms = 0.2
omgs = oms
rhom = 0.7
psiygap = 1.000
psipi = 2.500
rhoicu = 0.800
taudss = 0.01 # steady state tax on domestic consumption (setting it as 0 would create algebraic difficulties)
taumss = 0.01 # steady state tax on imported consumption for domestic country
taukss = 0.01 # steady state tax on rental income from capital for domestic country block
taunss = 0.01 # steady state tax on labor for domestic country
tauydss = 0.05
gss = 0.23 # steady state government spending as a propostion of gdp for domestic country block
gsss = 0.23 # steady state government spending as a propostion of gdp for foreign country block
taudsss = 0.01
taumsss = 0.01
tauksss = 0.01
taunsss = 0.01
tauydsss = 0.01 # steady state tax rate on output for foreign country block
tauss = 1.0 # Steady state terms of trade
icu = ((1+pidss)/beta) - 1
mc = ((ej - 1)/ej)
r = (1/taukss) * ((1/beta) - (1-delta))
rs = (1-tauksss) * ((1/beta) - (1-delta))
KN = (mc*alpha/r)**(1/(1-alpha))
KNs = (mc*alpha/rs)**(1/(1-alpha))
psigma = (1-xi) * (1/(1-tauydss) - xi)**(-1)
psigmas = (1-xi) * (1/(1-tauydsss) - xi)**(-1)
w = (1-alpha) * mc * (KN)**(alpha)
z = np.zeros(16)
def fun(z):
Yd = z[0]
N = z[1]
X = z[2]
I = z[3]
Cd = z[4]
Cm = z[5]
Gd = z[6]
Gm = z[7]
Yds = z[8]
Ns = z[9]
Xs = z[10]
Is = z[11]
Cds = z[12]
Cms = z[13]
Gds = z[14]
Gms = z[15]
print (z)
f = np.zeros(16)
f[0] = N - ( (X - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/r)** (1-taunss) )
f[1] = Yd - ( Cd + Gd + I + ((1-n)/n) *(Cms + Gms) )
f[2] = Yd - ( (KN)**(alpha) * (psigma/(1-tauydss)**(ej)) )
f[3] = Cd - ( X * ((1-om)/(1+taudss)**(el)) *((1-om)*(1+taudss)**(1-el) + om * (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[4] = Gd - ( ((gss*Yd * (1-omg))/(1+taudss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[5] = I - ( delta* KN * N )
f[6] = Cm -( (X * (1-om)/(1+tauydss)**(el) ) *((1-om)*(1+taudss)**(1-el) + om* (1+taumss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[7] = Gm - ( ((gss*Yd * (omg))/(1+taumss)**(eg) ) *((1-omg)*(1+taudss)**(1-eg) + omg* (1+taumss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[8] = Ns - ( (Xs - muc)**(-ec) * ((1-alpha)/(mun)) * (mc)**(1/(1-alpha)) * (alpha/rs)** (1-taunsss) )
f[9] = Yds - ( Cds + Gds + Is + (n/(1-n)) *(Cm + Gm) )
f[10] = Yds - ( (KNs)**(alpha) * (psigmas/(1-tauydsss)**(ej) ) )
f[11] = Cds - ( Xs * ((1-oms)/(1+taudsss)**(el))* ((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[12] = Gds - ( ((gsss*Yds * (1-omgs))/(1+taudsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
f[13] = Is - ( delta* KNs * Ns )
f[14] = Cms -( (Xs * (1-oms)/(1+tauydsss)**(el) ) *((1-oms)*(1+taudsss)**(1-el) + oms* (1+taumsss)**(1-el) * tauss**(1-el))**(el/(1-el)) )
f[15] = Gms - ( ((gsss*Yds * (omgs))/(1+taumsss)**(eg) ) *((1-omgs)*(1+taudsss)**(1-eg) + omgs* (1+taumsss)**(1-eg) * tauss**(1-eg))**(eg/(1-eg) ) )
return f
z = sp.optimize.root(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100], method='lm')
#z = fsolve(fun, [0,0,0.0,0,1,1,1,1,1,1,1,1,1,1,1,1])
print(z)
The solution are the following roots
success: True
x: array([ 3.64725445e-01, 1.02848541e-06, -1.86761721e+02,
9.52089296e-10, -1.30733205e+02, -1.25265418e+02,
5.87207967e-02, 2.51660557e-02, 3.36422990e+00,
5.18324506e-04, 8.17060628e+01, 4.87111630e-04,
6.53648502e+01, 6.53648502e+01, 6.19018302e-01,
1.54754576e-01])
Given an initial estimate of the root, a numerical root-finding algorithm moves along a certain direction in the variables-space until it finds a root. Obviously, with this approach, there is no point in requiring that the returned root is bounded within a certain interval - it all depends on how good the initial estimate is (and the search method used by the algorithm).
Another approach that can return bounded roots is to pose the root-finding problem as an optimization (e.g., minimization) problem, since providing constraints in an optimization problem is meaningful. However, you have to provide a properly objective function whose minimum occurs at a root of the original function (there are many options for such a function, typically the choice is heuristic).
One such function is the sum-of-squares f[0]**2 + f[1]**2 + ... + f[15]**2. Clearly, this function has a minimum value of zero, which is achieved when each of the individual terms of the sum are zero, i.e., at their roots. You can use Scipy's least_squares for performing this minimization, which also allows for providing bounds for the optimization variables.
Without any bounds on the variables and using the same initial root estimate, least_squares returns the same solution as root:
from scipy.optimize import least_squares
z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100])
print(z_ls.x)
print(z_ls.cost)
[ 3.6473e-01 1.0285e-06 -1.8676e+02 9.5209e-10 -1.3073e+02 -1.2527e+02 5.8721e-02 2.5166e-02 3.3642e+00 5.1832e-04 8.1706e+01 4.8711e-04 6.5365e+01 6.5365e+01 6.1902e-01 1.5475e-01] 4.16527754459e-26
(Note that z_ls.cost is the sum-of-squares evaluated at this point, which is (within numerical precision) zero.)
Now, use least_squares constraining the estimate to be non-negative:
z_ls = least_squares(fun, [100,100,70,30,50,20,50,20,100,100,100,100,100,100,100,100],bounds = (0,np.inf))
print(z_ls.x)
print(z_ls.cost)
[ 5.9581e-01 2.1229e+01 4.2108e-02 1.1820e-37 2.0493e-33 1.1914e-33 5.8857e-37 9.3812e-37 3.4508e+00 2.5054e+00 1.1516e+00 2.2395e+00 8.0630e-01 4.2258e-01 5.1994e-01 1.3867e-37] 0.237262813475
The returned estimate has indeed non-negative elements. However, z_ls.cost is (significantly) greater than zero, indicating that this solution is not a root. This implies one of the two:
If you do not have any insights on the above, the only thing you can do is try different initialization values and hope that a desired root is returned (either directly by root or a minimization formulation as the above).