Apols, below meant your suggestion to start at x0 = [.2, .05]. And yep I
still don't get convergence when starting there. Thanks for suggestion re
variable transformation.
On Wednesday, July 1, 2015 at 3:25:13 PM UTC+1, j verzani wrote:
>
> I'm not sure what starting points you are using, as there is a lot of
> "below" down there. The vpasolve solution is specifying a box to search in,
> unlike the nsolve method which just gets a starting point (matrix vs.
> vector). I looked at the mpmath.findroot functionality (where nsolve ends
> up) and doesn't look like there is a method for constrained problems. You
> might get away with a transformation of variables to re-express your
> domain. (Something like (pi/2+atan(x))/pi works well for the toy problem.)
>
> On Wednesday, July 1, 2015 at 9:54:07 AM UTC-4, David Zentler-Munro wrote:
>>
>> Hm, I've updated all packages and am still not getting convergence even
>> when starting with the guess shown below (or my matlab answer). Any ideas?
>> In the full code I will have 50 equations so it will not be practical to
>> graph solution for the starting guess.
>>
>> Finally, should I submit a pull request or issue to suggest incorporating
>> constraints into nsolve e.g. so I could specify that x lies between zero
>> and one?
>>
>> Thanks again for your help with this.
>>
>> David
>>
>> On Wednesday, July 1, 2015 at 2:19:41 PM UTC+1, j verzani wrote:
>>>
>>> No, that error is just the algorithm not converging. I would guess
>>> somewhere between your SymPy version and mine the `nsolve` routine changed
>>> a bit. However, with the answer given away, you can reverse engineer the
>>> convergence you want by starting nearby. Replace the `nsolve` call with
>>> this (so that it will work without updating to master):
>>>
>>> x0 = [.2, .05]
>>>
>>> sol=SymPy.sympy_meth(:nsolve, [eqn1, eqn2], [x1,x2], tuple(x0...))
>>>
>>>
>>>
>>> One way to get a decent initial guess would be to plot the functions.
>>> With PyPlot installed, the `plot_implicit` function shows roughly the
>>> location of your answer:
>>>
>>>
>>> plot_implicit(eqn1 * eqn2, (x1, 0, 1), (x2,0,1))
>>>
>>>
>>> The zeros of eqn1 and eqn2 seem to cross at (1,1) and your answer.
>>>
>>> On Wednesday, July 1, 2015 at 6:12:37 AM UTC-4, David Zentler-Munro
>>> wrote:
>>>>
>>>> Thanks loads J, this looks really helpful. You're right that the answer
>>>> looks out: x1 and x2 are unemployment rates so should be between zero and
>>>> one. The answer from matlab is (0.119, 0.055). Is there a way of putting
>>>> bounds on variables in nsolve?
>>>>
>>>> Also when I tried the code you suggest I got the error: " PyError
>>>> (:PyObject_Call) <type 'exceptions.ValueError'>ValueError('Could not find
>>>> root within given tolerance. (0.0742305 > 2.1684e-19)\nTry another
>>>> starting
>>>> point or tweak arguments.',)" Do I need to wait until you have updated
>>>> SymPy? (No worries if so!)
>>>>
>>>> Thanks,
>>>>
>>>> David
>>>>
>>>>
>>>>
>>>> On Tuesday, June 30, 2015 at 7:31:45 PM UTC+1, j verzani wrote:
>>>>>
>>>>> Here is how it can be done with SymPy (though I just had to check in a
>>>>> fix as there was an error with nsolve, so this will only work for now
>>>>> with
>>>>> master).
>>>>>
>>>>> Change the lines at the bottom to read:
>>>>>
>>>>> using SymPy
>>>>> x1,x2 = [symbols("x$i", real=true, positive=true) for i in 1:2]
>>>>>
>>>>> ps1= wbar +
>>>>> ((kappa*(1-beta*(1-sigma*((1-x1)/x1))))/(beta*((x1/(sigma*(1-x1)))^(gamma/(1-gamma)))*(1/(2-x1))));
>>>>> ps2= wbar +
>>>>> ((kappa*(1-beta*(1-sigma*((1-x2)/x2))))/(beta*((x2/(sigma*(1-x2)))^(gamma/(1-gamma)))*(1/(2-x2))));
>>>>>
>>>>> prod1=prod[1];
>>>>> prod2=prod[50];
>>>>> y1=(1-x1)*l1;
>>>>> y2=(1-x2)*l2;
>>>>>
>>>>> M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(psi/(psi-1));
>>>>>
>>>>> Mprime=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(1/(psi-1));
>>>>> K=((r/eta)^(1/(eta-1)))*M;
>>>>>
>>>>> pd1=(1-eta)*(K^eta)*(M^(-eta))*Mprime*((prod1*y1)^(-1/psi))*prod1;
>>>>> pd2=(1-eta)*(K^eta)*(M^(-eta))*Mprime*((prod2*y2)^(-1/psi))*prod2;
>>>>>
>>>>> eqn1=pd1-ps1;
>>>>> eqn2=pd2-ps2;
>>>>>
>>>>> sol=nsolve([eqn1, eqn2], [x1, x2],[1/2, 1/2])
>>>>>
>>>>>
>>>>> With this we get an answer, though I think it has escaped and isn't
>>>>> what you are looking for. (Perhaps you have a better starting point...)
>>>>>
>>>>>
>>>>> julia> sol
>>>>>
>>>>> 2x1 Array{Float64,2}:
>>>>>
>>>>> 214.956
>>>>>
>>>>> 300.851
>>>>>
>>>>>
>>>>> julia> subs(eqn1, (x1, sol[1]), (x2, sol[2]))
>>>>>
>>>>> -1.11022302462516e-16
>>>>>
>>>>>
>>>>> julia> subs(eqn2, (x1, sol[1]), (x2, sol[2]))
>>>>>
>>>>> -2.77555756156289e-17
>>>>>
>>>>> On Tuesday, June 30, 2015 at 1:20:41 PM UTC-4, David Zentler-Munro
>>>>> wrote:
>>>>>>
>>>>>> Thanks Tim. I know this is a little sacrilegious, but I decided to go
>>>>>> back to Matlab to see if I could solve my problem there. Eventually I
>>>>>> managed with the following code (note the code above was a simplified
>>>>>> version of my problem, the code below is the full version):
>>>>>>
>>>>>> beta = 0.95;
>>>>>> % Discount Factor(Krusell, Mukoyama, Sahin)
>>>>>> lambda0 = .90;
>>>>>> % Exog contact Rate @inbounds for unemployed to job (DZM)
>>>>>> lambda1 = 0.05;
>>>>>> % Exog job to job contact (DZM)
>>>>>> tao = .5;
>>>>>> % Proportion of job contacts resulting in move (JM Robin)
>>>>>> mu = 2;
>>>>>> % First shape parameter of beta productivity distribution (JM
>>>>>> Robin)
>>>>>> rho = 5.56;
>>>>>> % Second shape parameter of beta productivity distribution (JM
>>>>>> Robin)
>>>>>> xmin= 0.73;
>>>>>> % Lower bound @inbounds for beta human capital distribution (JM
>>>>>> Robin)
>>>>>> xmax = xmin+1;
>>>>>> % Upper bound @inbounds for beta human capital distribution (JM
>>>>>> Robin)
>>>>>> z0 = 0.0;
>>>>>> % parameter @inbounds for unemployed income (JM Robin)
>>>>>> nu = 0.64;
>>>>>> % parameter @inbounds for unemployed income (DZM)
>>>>>> sigma = 0.023;
>>>>>> % Job Destruction Rate (DZM)
>>>>>> alpha = 2;
>>>>>> % Coef of risk aversion in utility function (Krusell, Mukoyama,
>>>>>> Sahin)
>>>>>> TFP = 1;
>>>>>> % TFP (Krusell, Mukoyama, Sahin)
>>>>>> eta = 0.3;
>>>>>> % Index of labour @inbounds for CD production function (Krusell,
>>>>>> Mukoyama, Sahin)
>>>>>> delta = 0.01;
>>>>>> % Capital Depreciation Rate (Krusell, Mukoyama, Sahin)
>>>>>> tol1=1e-5;
>>>>>> % Tolerance parameter @inbounds for value function iteration
>>>>>> tol2=1e-5;
>>>>>> % Tolerance parameter @inbounds for eveything else
>>>>>> na=500;
>>>>>> % Number of points on household assets grid
>>>>>> ns=50;
>>>>>> % Number of points on human capital grid
>>>>>> amaximum=500;
>>>>>> % Max point on assets grid
>>>>>> maxiter1=10000;
>>>>>> maxiter=20000;
>>>>>> % Max number of iterations
>>>>>> bins=25;
>>>>>> % Number of bins @inbounds for wage and income distributions
>>>>>> zeta=0.97;
>>>>>> % Damping parameter (=weight put on new guess when updating in
>>>>>> algorithms)
>>>>>> zeta1=0.6;
>>>>>> zeta2=0.1;
>>>>>> zeta3=0.01;
>>>>>> mwruns=15;
>>>>>> gamma=0.5;
>>>>>> % Matching Function Parameter
>>>>>> kappa = 1;
>>>>>> % Vacancy Cost
>>>>>> psi=0.5;
>>>>>> prod=linspace(xmin,xmax,ns);
>>>>>> % Human Capital Grid (Values)
>>>>>> l1=0.7;
>>>>>> l2=0.3;
>>>>>> wbar=0.2;
>>>>>> r=((1/beta)-1-1e-6 +delta);
>>>>>>
>>>>>> syms x1 x2
>>>>>>
>>>>>> ps1= wbar +
>>>>>> ((kappa*(1-beta*(1-sigma*((1-x1)/x1))))/(beta*((x1/(sigma*(1-x1)))^(gamma/(1-gamma)))*(1/(2-x1))));
>>>>>> ps2= wbar +
>>>>>> ((kappa*(1-beta*(1-sigma*((1-x2)/x2))))/(beta*((x2/(sigma*(1-x2)))^(gamma/(1-gamma)))*(1/(2-x2))));
>>>>>>
>>>>>> prod1=prod(1);
>>>>>> prod2=prod(50);
>>>>>> y1=(1-x1)*l1;
>>>>>> y2=(1-x2)*l2;
>>>>>>
>>>>>> M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(psi/(psi-1));
>>>>>>
>>>>>> Mprime=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(1/(psi-1));
>>>>>> K=((r/eta)^(1/(eta-1)))*M;
>>>>>>
>>>>>> pd1=(1-eta)*(K^eta)*(M^(-eta))*Mprime*((prod1*y1)^(-1/psi))*prod1;
>>>>>> pd2=(1-eta)*(K^eta)*(M^(-eta))*Mprime*((prod2*y2)^(-1/psi))*prod2;
>>>>>>
>>>>>> eqn1=pd1-ps1;
>>>>>> eqn2=pd2-ps2;
>>>>>>
>>>>>> sol=vpasolve([0==eqn1, 0==eqn2], [x1 x2],[0 1;0 1]);
>>>>>> [sol.x1 sol.x2]
>>>>>>
>>>>>> If anyone has an idea how I can replicate this e.g. the vpasolve
>>>>>> call, in Julia I would be very grateful.
>>>>>>
>>>>>> Best,
>>>>>>
>>>>>> David
>>>>>>
>>>>>> On Tuesday, June 30, 2015 at 2:27:23 PM UTC+1, Tim Holy wrote:
>>>>>>>
>>>>>>> You could have your function return NaN when one or more arguments
>>>>>>> are out-of-
>>>>>>> bounds.
>>>>>>>
>>>>>>> If you want to find a solution very near such a boundary, your
>>>>>>> performance will
>>>>>>> presumably stink. I'm not sure what the state of the art here is
>>>>>>> (for
>>>>>>> minimization, I'd suggest an interior-point method).
>>>>>>>
>>>>>>> Best,
>>>>>>> --Tim
>>>>>>>
>>>>>>> On Tuesday, June 30, 2015 06:18:09 AM David Zentler-Munro wrote:
>>>>>>> > 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
>>>>>>>
>>>>>>>
