Final armington case.

Author: SebKrantz

src/main/optimal_network.jl

@@ -81,8 +81,13 @@ function optimal_network(param, graph; I0=nothing, Il=nothing, Iu=nothing, verbo
             recover_allocation = recover_allocation_mobility_cgc
         end
     elseif param.mobility == 0.5 && param.cong
-        model = model_partial_mobility_cgc(optimizer, auxdata)
-        recover_allocation = recover_allocation_partial_mobility_cgc
+        if all(sum(param.Zjn .> 0, dims = 2) .<= 1) # Armington case
+            model = model_partial_mobility_cgc_armington(optimizer, auxdata)
+            recover_allocation = recover_allocation_partial_mobility_cgc_armington
+        else
+            model = model_partial_mobility_cgc(optimizer, auxdata)
+            recover_allocation = recover_allocation_partial_mobility_cgc
+        end
     elseif param.mobility == 0 && param.cong
         if all(sum(param.Zjn .> 0, dims = 2) .<= 1) # Armington case
             model = model_fixed_cgc_armington(optimizer, auxdata)

src/models/model_partial_mobility_cgc_armington.jl

@@ -0,0 +1,97 @@
+# The primal case, with mobility within regions and cross-good congestion, and an Armington (1969) world where each location produces only one good
+
+function model_partial_mobility_cgc_armington(optimizer, auxdata)
+
+    # Extract parameters
+    param = auxdata.param
+    graph = auxdata.graph
+    region = graph.region
+    if length(region) != graph.J
+        error("length(region) = $(length(region)) does not match number of nodes = $(graph.J)")
+    end
+    kappa_ex_init = auxdata.kappa_ex
+    A = auxdata.edges.A
+    Apos = auxdata.edges.Apos
+    Aneg = auxdata.edges.Aneg
+    m = param.m # Vector of weights on each goods flow for aggregate congestion term
+    psigma = (param.sigma - 1) / param.sigma
+    beta_nu = (param.beta + 1) / param.nu
+    Lr = param.Lr
+    if length(param.omegar) != param.nregions
+        error("length(param.omegar) = $(length(param.omegar)) does not match number of regions = $(param.nregions)")
+    end
+    if length(Lr) != param.nregions
+        error("Populations Lr = $(length(Lr)) does not match number of regions = $(param.nregions)")
+    end
+
+    # Model
+    model = Model(optimizer)
+    set_string_names_on_creation(model, false)
+
+    # Variable declarations
+    @variable(model, ur[1:param.nregions], container=Array, start = 0.0)                          # Utility per capita in each region
+    @variable(model, Djn[1:graph.J, 1:param.N] >= 1e-8, container=Array, start = 1e-6)            # Consumption per good pre-transport cost (Dj)
+    @variable(model, cj[1:graph.J] >= 1e-8, container=Array, start = 1e-6)                        # Overall consumption bundle, including transport costs
+    @variable(model, Qin_direct[1:graph.ndeg, 1:param.N] >= 1e-8, container=Array, start = 0.0)   # Direct aggregate flow
+    @variable(model, Qin_indirect[1:graph.ndeg, 1:param.N] >= 1e-8, container=Array, start = 0.0) # Indirect aggregate flow
+    @variable(model, Lj[1:graph.J] >= 1e-8, container=Array)                        # Overall labour
+    # Calculate start values for Lj
+    pop_start = (Lr ./ gsum(ones(graph.J), param.nregions, region))[region]
+    set_start_value.(Lj, pop_start)
+
+    # Parameters: to be updated between solves
+    @variable(model, kappa_ex[i = 1:graph.ndeg] in Parameter(i), container=Array)
+    set_parameter_value.(kappa_ex, kappa_ex_init)
+
+    # Objective
+    @expression(model, U, sum(param.omegar .* Lr .* ur)) # Overall Utility
+    @objective(model, Max, U)
+
+    # Utility constraint (Lj * ur <= ... )
+    @constraint(model, Lj .* ur[region] - (cj .* Lj ./ param.alpha) .^ param.alpha .* (param.Hj ./ (1 - param.alpha)) .^ (1 - param.alpha) .<= -1e-8)
+
+    # Create the matrix B_direct (resp. B_indirect) of transport cost along the direction of the edge (resp. in edge opposite direction)
+    B_direct = @expression(model, ((Qin_direct .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
+    B_indirect = @expression(model, ((Qin_indirect .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
+    # Final good constraints
+    @expression(model, Dj, sum(Djn .^ psigma, dims=2) .^ (1 / psigma))
+    @constraint(model, cj .* Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
+
+    # Balanced flow constraints: same as full mobility
+    @expression(model, Yjn, param.Zjn .* (Lj .^ param.a))
+    @constraint(model, Pjn, Djn + A * Qin_direct - A * Qin_indirect - Yjn .<= -1e-8)
+
+    # Labor resource constraints (within each region)
+    @constraint(model, -1e-8 .<= gsum(Lj, param.nregions, region) .- Lr .<= 1e-8)
+    
+    return model
+end
+
+function recover_allocation_partial_mobility_cgc_armington(model, auxdata)
+    param = auxdata.param
+    graph = auxdata.graph
+    model_dict = model.obj_dict
+    results = Dict()
+
+    results[:welfare] = value(model_dict[:U])
+    results[:reg_pc_welfare] = value.(model_dict[:ur])
+    results[:Yjn] = value.(model_dict[:Yjn])
+    results[:Yj] = dropdims(sum(results[:Yjn], dims=2), dims = 2) 
+    results[:Lj] = value.(model_dict[:Lj])
+    results[:Ljn] = (param.Zjn .> 0) .* results[:Lj]
+    results[:Djn] = value.(model_dict[:Djn]) # Consumption per good pre-transport cost
+    results[:Dj] = dropdims(value.(model_dict[:Dj]), dims=2)
+    results[:cj] = value.(model_dict[:cj])
+    results[:Cj] = results[:cj] .* results[:Lj]
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:uj] = param.u.(results[:cj], results[:hj])
+    # Prices
+    results[:Pjn] = shadow_price.(model_dict[:Pjn])
+    results[:PCj] = dropdims(sum(results[:Pjn] .^ (1-param.sigma), dims=2), dims=2) .^ (1/(1-param.sigma))    
+    # Network flows
+    Qin_direct = value.(model_dict[:Qin_direct])
+    Qin_indirect = value.(model_dict[:Qin_indirect])
+    results[:Qin] = ifelse.(Qin_direct .> Qin_indirect, Qin_direct, -Qin_indirect)
+    results[:Qjkn] = gen_network_flows(results[:Qin], graph, param.N)
+    return results
+end