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