Thanks Mauro. Fair point about complex code, and yes "ns" was missing. I've tried to simplify code as follows below (still get same error). As far as I can see it does seem like NLsolve is evaluating outside the bounds so, unless anyone has other suggestions, I will file an issue. Here's the updated code
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 > ns=50 > prod=linspace(xmin,xmax,ns) > l1=0.7 > l2=0.3 > 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])))) > ps2= wbar + (kappa*(1-beta*(1-sigma*((1-x[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))) > K=((r/eta)^(1/(eta-1)))*M > > pd1=(1-eta)*(K^eta)*(M^(-eta))*prod1 > > pd2=(1-eta)*(K^eta)*(M^(-eta))*prod2 > > fvec[1]=pd1-ps1 > fvec[2]=pd2-ps2 > end > > mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3, 0.3]) > > > On Tuesday, June 30, 2015 at 1:32:03 PM UTC+1, Mauro wrote: > > Maybe you could use complex numbers, so the domain error isn't thrown, > then take the real part at the end (after checking that im==0)? > Alternatively, you could file an issue at NLsolve stating that the > objective function is evaluated outside the bounds. Maybe it is > possible improve the algorithm, maybe not. > > (Also, to get better feedback and be kinder to the reviewers, try to > make a smaller example which does only uses the minimal number of > dependencies. And make sure it works as, I think, your example does not > work, unless `ns` is something which is exported by one of the used > packages I don't have installed.) > > On Tue, 2015-06-30 at 07:01, David Zentler-Munro < > [email protected] <javascript:>> wrote: > > '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 > >
