Solve system of nonlinear equations

New to Julia
1

H everyone

Apologies in advance for what may seem like a basic question, but I’m fairly new to julia. I’m looking to solve a system of nonlinear equations – currently I have 47 equations for 47 variables – and am looking for a solution more effective than MATLAB’s fsolve (unable to solve my problem) as I’m looking to find a steady state.

I think NLsolve would be most suitable for my needs, but I’m a bit unsure on some of the notation and syntax of julia.

Would the following format be appropriate for my problem:

using NLsolve

function f!(F, x)
F[1] = …
F[2] = …
…
F[47] = …
end

nlsolve(f!, initial_x, autodiff = :forward)

Where x is basically a vector of my endogenous variables, and initial_x is a vector of starting guesses for the vector of endogenous variables.

Any help or advice would be much appreciated.

2

Yes that looks right.

3

It is very unlikely that just switching to Julia will help you solve your problem, as the same general algorithms are used by pretty much all simple multivariate solvers. Usually the best approach is to look at the structure of your problem, and see if you can solve for blocks in variables (ie a block-triangular structure) independently of others, then break it up into simple ones. Good initial guesses also help a lot, which also need understanding of the problem.

Failing that, homotopy continuation methods may work, but ATM I am only aware of a Julia package for polynomials.

4

Yup, no general purpose continuation methods in Julia.

I think your (OP) best bet is to maybe show us the problem you’re trying to solve, because I think it’s very likely that NLsolve will fail it fsolve fails.

1 Like
5

I have yet to code up the problem in julia, but I’ve just copied over the problem I’m trying to solve over from Atom. It’s basically a simple macroeconomic general equilibrium programming problem. Here’s the following system of equations:

            alpha = 1 - 0.623  
            n = 1.01^0.25 - 1  
            phiI = 132 
            kappa = 12 
            phiK = 1 
            delta = 0.025  
            omega = 1.023^(-0.25) 
            gamma = 1.1^(-0.25)  
            upsilon = 0.4  
            rho = -3  
            beta = .999 
            xi = 0.6  
            g = 0.17  
            b = 4  
            varrho = 0.4  
            phiPI = 2  
            phiR = 0.5  
            Rbar = (1/beta) - 1 
            phiNX = 1  
            ystar = 1 
            Rstar = (1/0.99 - 1) 
            phiS = 0.000742  
            fstarbar = -0.5  

function f!(F,x)
F[1]=x[2]^(alpha)*x[3]^(1-alpha) - x[1] 

F[2]=(x[5]/(1-alpha))^(1-alpha)*(x[6]/alpha)^alpha - x[4] 
     
F[3)=(1-alpha)*x[4]*x[1]/x[3] - x[5] 

F[4]=alpha*x[4]*x[1]/x[2] - x[6] 

F[5]=(1-x[4]-(phiI/2)*(x[7]-1)^2)*x[1] - x[10] 

F[6]=kappa*(x[4]-1)*x[1] + (1+n)*(x[7]/x[6])*phiI*(x[7]-1)*x[7]*x[1] - phiI*(x[7]-1)*x[7]*x[1]
     
F[7]=x[13]*x[12] - (1-phiK*((x[12]/x[12])-1)^2)*x[12] - x[11] 
     
F[8]=((1-delta)*x[2]+(1-(phiK*(x[12]/x[12])-1)^2)*x[12]) - x[2] 
     
F[9]=(x[7]/x[14])*(x[13]*(1-delta)+x[6]) - x[13] 
     
F[10]=(1-omega)+gamma*x[15] - (1+n)*x[15]
     
F[11]=(x[17]*x[18])^(-1/rho)*x[19]*(((1-upsilon)/upsilon)*(1/(xi*x[5])))^(1-upsilon) - x[16] 
     
F[12]=(x[17]*x[18])*(((x[14]*x[20])/x[7])+x[21]+x[22]-xi*x[25]) - x[19]
     
F[13]=x[23] + ((gamma*x[7])/x[14])*x[22] - x[22] 
     
F[14]=varrho*(x[5] - x[34]) - x[23] 
     
F[15]=x[15] - ((1-upsilon)/upsilon)*(1/(xi*x[5]))*x[19] - x[24] 
     
F[16]=xi*x[5]*x[24] + ((gamma*x[7])/x[14])*x[21] - x[21] 
     
F[17]=xi*x[25] + ((gamma*x[7])/x[14])*x[26] - x[26] 
     
F[18]=((x[14]*x[20])/x[7]) - x[19] + (1+omega)*(((x[14]*x[29])/x[7])+ x[5] - x[28]) - x[20] 
     
F[19]=1 - ((x[17]*x[18])/(x[17]*x[18]))*gamma*x[14]^(rho/(1-rho))*beta^(rho/(1-rho))*(1/x[7])^((upsilon*rho)/(1-rho))*(x[5]/x[5])^(((1-upsilon)*rho)/(1-rho)) - (x[17]*x[18])
     
F[20]=(x[18])^(-1/rho)*x[28]*(((1-upsilon)/upsilon)*(1/x[5]))^(1-upsilon) - x[27]
     
F[21]=x[18]*(((x[14]*x[29])/x[7])+ x[30] + x[31] - x[34]) - x[28]
     
F[22]=((x[7]*(x[32]-omega))/(x[32]*x[14]))*x[22] + ((x[7]*omega)/(x[32]*x[14]))*x[31] - x[31] 
     
F[23]=1 - ((1-upsilon)/upsilon)*(1/x[5])*x[28] - x[33]
     
F[24]=x[5]*x[33] + omega*((x[30]*x[7])/(x[32]*x[14])) + (1-omega)*(x[17])^((rho-1)/rho)*(1/xi)^(1-upsilon)*((x[21]*x[7])/(x[32]*x[14])) - x[30] 
     
F[25]=x[34] + omega*((x[35]*x[7])/(x[32]*x[14])) + (1-omega)*(x[17])^((rho-1)/rho)*(1/xi)^(1-upsilon)*((x[26]*x[7])/(x[32]*x[14])) - x[35] 
     
F[26]=omega*((x[14]*x[29])/x[7] + x[5]*x[33] - x[28]) - x[29] 
     
F[27]=1 - beta^(1/(1-rho))*(1/x[7])^((rho*upsilon)/(1-rho))*x[14]^(rho/(1-rho))*x[32]^(rho/(1-rho))*(x[5]/x[5])^(((upsilon-1)*rho)/(1-rho)) - x[18] 
     
F[28]=omega + (1-omega)*x[17]^(-(1-rho)/rho)*(1/xi)^(1-upsilon) - x[32] 
     
F[29]=x[19] + x[28] - x[36] 
     
F[30]=x[20]/x[37] - x[38] 
     
F[31]=omega*(1 - x[17]*x[18])*x[38]*((x[14]*A)/x[7]) - (1+n)*(x[38] - (1-omega))*x[37] 
     
F[32]=x[16] + x[27] - x[39] 
     
F[33]=(x[14]/x[7])*x[40] + x[41] + x[42] - x[43] - (1-n)*x[40] 
     
F[34]=x[34] + xi*x[25] - x[43] 
     
F[35]=g - x[41] 
     
F[36]=b - (1+n)*x[40] 
     
F[37]=varrho*PSI*((1-alpha)*x[4]*x[1] - x[43]) - x[42] 
     
F[38]=phiR*Rbar + (1-phiR)*x[14] + phiPI*x[7] - x[14] 

F[39]=(1+n)*(x[9]+x[11])/x[9] - x[14]/x[7] 
     
F[40]=(1+n)*(x[8]+x[10])/x[8] - x[14]/x[7] 
     
F[41]=x[33] + xi*x[24] - x[3] 
     
F[42]=x[36] + x[12] + x[41] + x[45] - x[1]*(1 - (phiI/2)*(x[7]-1)^2 - phiK*((x[12]/x[12])-1)^2) 
     
F[43]=x[44]^(phiNX)*ystar - x[45] 
     
F[44]=x[2] + x[40] + x[8] + x[44]*x[47] - x[37] 
     
F[45]=(1 + Rstar + x[46])*(x[44]/x[44]) - 1 - x[14] 
     
F[46]=(1 + Rstar + x[46])*x[44]*x[47] + x[45] - x[44]*x[47] 
     
F[47]=phiS*(exp(x[47]-fstarbar)-1) - x[46]
6

I am a little bit confused, I see 23 things which look like parameters (alpha, n, phiI, …) and no variables (accesses to x). How do you get the 47 variables?

7

Variables are in all caps probably?

8

Oh that is it probably. So maybe one could formulate this as a polynomial system after some variable substitution and clearing of denominators.

1 Like
9

Usually, thinking about the mathematical structure of the problem first (before programming) has large and immediate payoffs. This looks like Euler and technology equation residuals from a DSGE model, so maybe you can derive the steady state capital, then consumption and labor, etc.

1 Like
10

Apologies, the code was essentially pulled out of a MATLAB script. I’ll edit my post now.

11

Correct, it’s a DSGE model.

Worst case scenario, I may have to manually log linearise the model and put into Dynare, setting initial values equal to zero.

12

As I said, it should be possible to simplify a lot by understanding the model. Brute-forcing your way to a steady state will be costly (it is actually the most costly part of solving large-ish DSGE models by perturbation methods; second is the linear part, the rest is pretty mechanical).

You can find examples on how to do this in most texts, eg the handbook chapter of Fernández-Villaverde et al. (Sorry if you already know this).

1 Like
13

One thing to consider is the range of the valid parameter space. If you have some variables, like variances, which have to be positive, solve for the log of that variable and exp it in the objective.

Related, if there are outputs like consumption and you have log utility, then the marginal utility 1/c only makes sense for positive consumption. You might try truncating consumption to 0 so you get an Inf (and fail on that particular guess), instead of crossing the discontinuity at 0, which might make the optimizer diverge.

1 Like