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