Solve a system of N equations with N unknowns using Julia

383 Views Asked by At

I have :

  • a set of N locations which can be workplace or residence
  • a vector of observed workers L_i, with i in N
  • a vector of observed residents R_n, with n in N
  • a matrix of distance observed between all pair residence n and workplace i
  • a shape parameter epsilon

Setting N=3, epsilon=5, and

d = [1 1.5 3 ; 1.5 1 1.5 ; 3 1.5 1] #distance matrix
L_i = [13  69  18] #vector of workers in each workplace
R_n = [27; 63; 10]

I want to find the vector of wages (size N) that solve this system of N equations, enter image description here

with l all the workplaces.

Do I need to implement an iterative algorithm on the vectors of workers and wages? Or is it possible to directly solve this system ?

I tried this,

w_i = [1 ; 1 ; 1]
er=1
n =1
while  er>1e-3
    L_i =   ( (w_i ./ d).^ϵ ) ./ sum( ( (w_i ./ d).^ϵ), dims=1) * R 
    er = maximum(abs.(L .- L_i))
    w_i = 0.7.*w_i + 0.3.*w_i.*((L .- L_i) ./ L_i)
    n = n+1
end
1

There are 1 best solutions below

6
Fabrizio Leone On BEST ANSWER

If L and R are given (i.e., do not depend on w_i), you should set up a non-linear search to get (a vector of) wages from that gravity equation (subject to normalising one w_i, of course).

Here's a minimal example. I hope it helps.

# Call Packages
using Random, NLsolve, LinearAlgebra

# Set seeds
Random.seed!(1704)

# Variables and parameters
N    = 10
R    = rand(N)
L    = rand(N) * 0.5
d    = ones(N, N) .+ Symmetric(rand(N, N)) / 10.0 
d[diagind(d)] .= 1.0
ε    = -3.0

# Define objective function
function obj_fun(x, d, R, L, ε) 
    
    # Find shares
    S_mat  = (x ./ d).^ε
    den    = sum(S_mat, dims = 1)
    s      = S_mat ./ den

    # Normalize last wage
    x[end] = 1.0

    # Define loss function
    loss   = L .- s * R

    # Return 
    return loss
    
end

# Run optimization
x₀   = ones(N) 
res  = nlsolve(x -> obj_fun(x, d, R, L, ε), x₀, show_trace = true)

# Equilibrium vector of wages
w    = res.zero