Plotting sensitivity of variable to parameter

Specific Domains Optimization (Mathematical)
4

Here’s how I’d write your model:

julia> using JuMP

julia> import Distributions

julia> import Ipopt

julia> import QuadGK

julia> function main(;
           num_firms = 40,
           θ = 2.44,
           μᴹᶜ = 1.62,
           σᴹᶜ = 0.73,
           μᵣᴹᶜ = 0.0,
           αᵘ = 1.14,
           u₀ = 1.17,
       )
           num_firms_half = div(num_firms, 2)
           χ_dist = Distributions.Normal(0, 1)
           incumbent_advantage = [fill(μᵣᴹᶜ, num_firms_half); zeros(num_firms_half)]
           function expected_MC(j)
               expected_cost, _ = QuadGK.quadgk(-4, 4; rtol = 1e-4) do χ
                   mcj = μᴹᶜ + incumbent_advantage[j]
                   mci = 4.0 * Distributions.cdf(Distributions.Normal(0, 1), χ)
                   return (mcj + mci) * Distributions.pdf(χ_dist, χ)
               end
               return expected_cost
           end
           function market_share(prices, firm_idx)
               share, _ = QuadGK.quadgk(-4, 4; rtol = 1e-4) do χ
                   utilities = map(1:num_firms) do k
                       return -αᵘ * prices[k] + ifelse(k <= num_firms_half, θ, 0.0)
                   end
                   exp_uj = exp(utilities[firm_idx])
                   exp_sum = sum(exp, utilities) + exp(u₀)
                   choice_prob = exp_uj / exp_sum        
                   return choice_prob * Distributions.pdf(χ_dist, χ)
               end
               return share
           end
           market_share_jump(j, x...) = market_share(collect(x), round(Int, j))
           expected_MCs = expected_MC.(1:num_firms)
           model = Model(Ipopt.Optimizer)
           set_silent(model)
           @variable(
               model,
               expected_MCs[j] <= p[j in 1:num_firms] <= 50,
               start = 1.5 * expected_MCs[j],
           )
           @operator(model, op_market_share, 1 + num_firms, market_share_jump)
           @expression(
               model,
               profit[j in 1:num_firms], 
               (p[j] - expected_MCs[j]) * op_market_share(j, p...),
           )
           @objective(model, Max, sum(profit))
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return value.(p)
       end
main (generic function with 1 method)

julia> solutions = [x => main(; μᴹᶜ = x) for x in range(1.0, 2.0; length = 11)];

julia> map(solutions) do (x, y)
           return x => extrema(y)
       end
11-element Vector{Pair{Float64, Tuple{Float64, Float64}}}:
 1.0 => (4.35244388352189, 4.352443948268603)
 1.1 => (4.418130951995438, 4.418134349252176)
 1.2 => (4.485512143991204, 4.485512166093286)
 1.3 => (4.554574506199629, 4.554575449707015)
 1.4 => (4.625314534369957, 4.625314701417606)
 1.5 => (4.697712923216194, 4.697713471679339)
 1.6 => (4.771746944243104, 4.771752086201073)
 1.7 => (4.847386725269032, 4.847387100544168)
 1.8 => (4.9245956732298195, 4.924596132256306)
 1.9 => (5.003329940153706, 5.0033315602162824)
 2.0 => (5.08354131761681, 5.083541948581439)

I assume that I got things correct, but I may have made a typo somewhere, etc.

Note that @NL is now legacy syntax. There’s a new nonlinear interface: Nonlinear Modeling · JuMP