Non negative roots for a system of non linear equations (python)

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])

=================

  

 

Is there any point in the code where you tell the solver that for your solution, you nee the constraint zi≥0zi≥0z_i\geq 0?
– Roland
Oct 20 at 23:14

  

 

No, Is there a way I can include non-negativity constraints without messing up with the system consistency?
– Nck
Oct 20 at 23:15

  

 

You could try to minimize the function with scipy.optimize.minimize to see if you can get solutions of the system where the variables are nonegative, but the solution is at most zero.
– Roland
Oct 20 at 23:18

  

 

i have a vector of functions, a system. I am not sure if minimize can do this work. Though my familiarity with python is at best amature.
– Nck
Oct 20 at 23:23

  

 

@Nck, you can minimize the RSS
– uranix
2 days ago

=================

=================