'm trying to solve a large number (50) of non-linear simultaneous equations
in Julia. For the moment I'm just trying to make this work with 2 equations
to get the syntax right etc. However, I've tried a variety of
packages/tools - NLsolve, nsolve in SymPy and NLOpt in JuMP (where I ignore
the objective function and just enter equality constraints)- without much
luck. I guess I should probably focus on making it work in one. I'd
appreciate any advice on choice of packages and if possible code.
Here's how I tried to do it in NLsolve (using it in mcpsolve mode so I can
impose constraints on the variables I am solving for - x[1] and x[2] -
which are unemployment rates and so bounded between zero and 1) :
using Distributions
using Devectorize
using Distances
using StatsBase
using NumericExtensions
using NLsolve
beta = 0.95
xmin= 0.73;
xmax = xmin+1;
sigma = 0.023;
eta = 0.3;
delta = 0.01;
gamma=0.5
kappa = 1
psi=0.5
prod=linspace(xmin,xmax,ns)
l1=0.7
l2=0.3
assets = linspace(0,amaximum,na);
wbar=1
r=((1/beta)-1-1e-6 +delta)
## Test code
function f!(x, fvec)
ps1= wbar + ((kappa*(1-beta*(1-sigma*((1-x[1])/x[1]))))/
(beta*((x[1]/(sigma*(1-x[1])))^(gamma/(1-gamma)))*(1/(2-x[1]))))
ps2= wbar + ((kappa*(1-beta*(1-sigma*((1-x[2])/x[1]))))/
(beta*((x[2]/(sigma*(1-x[2])))^(gamma/(1-gamma)))*(1/(2-x[2]))))
prod1=prod[1]
prod2=prod[50]
y1=(1-x[1])*l1
y2=(1-x[2])*l2
M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(psi/(psi-1))
K=((r/eta)^(1/(eta-1)))*M
pd1=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+
((prod2*y2)^((psi- 1)/psi)))^(1/(psi-1)))*
((prod1*y1)^(-1/psi))*prod1
pd2=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+
((prod2*y2)^((psi-1)/psi)))^(1/(psi-1)))*
((prod2*y2)^(-1/psi))*prod2
fvec[1]=pd1-ps1
fvec[2]=pd2-ps2
end
mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3, 0.3])
However, I get this error
[image: error message]
As I understand it, this occurs because I am trying to take the root of a
negative number. However, the bounds on x[1] and x[2] should stop this
happening. Any advice very welcome. I appreciate the formulas are pretty
ugly so let me know if any further simplifications helpful (I have tried!)
David
