Market tâtonnement in two-score admissions model

TwoScoreAdmissionsEquilibrium.gif

TwoScoreAdmissionsEquilibrium.jl

#=  School-choice markets are like any economic market in that they are characterized
    by prices (score cutoffs), supply constraints (school capacity), and consumer demand
    (student interest). But a unique challenge in school-choice is how to model the
    preferences of schools over students and vice-versa. In this example, I use
    the multinomial logit (MNL) choice model for student preferences and assume
    that schools score students using on a convex combination of orthogonal
    components of their characteristics. I visualize the polyhedra corresponding
    to possible sets of admitted schools and the price path.  =#

using Combinatorics
using Plots, Plots.PlotMeasures
using Polyhedra
using LinearAlgebra

# Instance data
m = 2                       # Number of schools
qualities = [1.1, 0.4]      # Preferability of each school
blends = [.2 .8;
          .6 .4]            # School assessment functions
capacities = [.4, .2]       # Capacity of each school

# Demand function (defined below)
demand(cut) = demands_MNL_ttests(qualities, blends, cut)

init_cutoffs = [0.1, 0.1]   # Initial cutoffs
T = 76                      # Number of iterations to do

"""
    demands_MNL_ttests(qualities, blends, cutoffs)

Return demand for each school given a set of cutoffs and ignoring capacity, using
multinomial logit choice model with one student profile and `t` test scores,
which represent orthogonal dimensions of student quality. Each school uses
a convex combination of the tests scores to evaluate students; the rows of `blends`
give these combinations. `qualities` are the MNL preferability parameters.

Simplification of `demands_pMNL_ttests()` from `DeferredAcceptance` package.
"""
function demands_MNL_ttests(qualities   ::AbstractArray{<:AbstractFloat, 1},
                            blends      ::AbstractArray{<:AbstractFloat, 2},
                            cutoffs     ::AbstractArray{<:AbstractFloat, 1};
                            )::AbstractArray{<:AbstractFloat, 1}

    (m, t) = size(blends)
    demands = zeros(m)
    γ = exp.(qualities)
    bounds = vcat([HalfSpace(-col, 0.) for col in eachcol(I(t))],
                  [HalfSpace(1col, 1.) for col in eachcol(I(t))])

    for C♯ in powerset(1:m) # Possible choice sets
        if !isempty(C♯)
            # Probability of having this choice set is the volume of
            # this m-dimensional polyhedron.
            hspaces = [c in C♯ ? HalfSpace(-blends[c, :], -cutoffs[c]) :
                          HalfSpace( blends[c, :],  cutoffs[c])
                       for c in 1:m]
            poly = polyhedron(hrep(vcat(bounds, hspaces)))
            vol = volume(poly)

            if vol > 0
                mult = vol / sum(γ[C♯])

                for c in C♯
                    demands[c] += mult * γ[c]
                end
            end
        end
    end

    return demands
end

"""
Tatonnement procedure for school-choice market. If a school gets too many applicants,
it raises its cutoff; if it gets too few, it lowers it.

Simplification of `nonatomic_tatonnement()` from `DeferredAcceptance` package,
and modified to return the full price and demand sequence.
"""
function tatonnement(init_cutoffs::AbstractArray{<:AbstractFloat, 1};
                     α           ::Union{AbstractFloat, AbstractArray{<:AbstractFloat, 1}}=1.,
                     verbose     ::Bool=false,
                     maxit       ::Int=500)

    UB = max.(0., 1 .- capacities)
    cut = zeros(Float64, m, maxit)
    dem = zeros(Float64, m, maxit)
    cut[:, 1] = init_cutoffs

    for nit in 1:maxit-1
        verbose ? println("Round $nit") : nothing
        dem[:, nit] = demand(cut[:, nit])
        excess_demand = dem[:, nit] - capacities

        cut[:, nit + 1] = max.(0, min.(UB, cut[:, nit] + α .* excess_demand))
    end

    return cut, dem
end



# Get the time series
cut_container, dem_container = tatonnement(init_cutoffs, α=0.08, verbose=false, maxit=T)

function p_closure()

    P = plot(xlabel="Antarctic Institute of Technology",
             ylabel="Lunar University of Foreign Studies",
             title="Cutoff and demand paths",
             xlim=(0,1), ylim=(0,1),
             legend=:topright)

    # Gray arrows indicating excess demand
    verts = [(-1, -0.5), (1, 0), (-1, 0.5)]
    triangle(ang) = Shape([(x*cos(ang) - y*sin(ang), x*sin(ang) + y*cos(ang))
                           for (x, y) in verts])

    meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)))
    stp = 0.04
    x, y = meshgrid(0.0+stp:stp:1.0-stp, 0.0+stp:stp:1.0-stp)

    len = size(x)[1]
    for i in 1:len
          u, v = demand([x[i], y[i]]) - capacities
          ang = atan(v, u)
          scatter!(P, [x[i]], [y[i]],
                   marker=triangle(ang),
                   color=:black,
                   label=nothing,
                   msw=0,
                   ms=25*sqrt(u^2 + v^2),
                   ma=0.1)
    end

    return P
end

function q_closure()

    Q = plot(xlim=(0,1),
        ylim=(0,1),
        xlabel="Language score",
        ylabel="Math score",
        title="Dynamic cutoffs in two-score admissions model:\nAdmissions polyhedra",
        legend=:topright)
end

function anim_doer()
    P = p_closure()
    Q = q_closure()

    anim = @animate for it = vcat(1:T-1, repeat([T-1], 25))
        p = deepcopy(P)
        q = deepcopy(Q)

        plot!(p, cut_container[1, 1:it],
                 cut_container[2, 1:it],
                 label="Cutoff",
                 shape=:circle,
                 color=:goldenrod,
                 msw=0,
                 ms=2.5,
                 lc=:black,
                 la=0.6)

        plot!(p,
              dem_container[1, 1:it],
              dem_container[2, 1:it],
              label="Demand",
              shape=:circle,
              color=:crimson,
              msw=0,
              ms=2.5,
              lc=:black,
              la=0.6)

        plot!(p, [capacities[1], capacities[1]], [0, 1], st=:line, lc=:black, ls=:dash, la=0.4,
              label="Capacity")
        plot!(p, [0, 1], [capacities[2], capacities[2]], st=:line, lc=:black, ls=:dash, la=0.4,
              label=nothing)

        annotate!(q, 0.95, 0.8, text("t = $it", :white, :right, 10))

        for (c, blend) in enumerate(eachrow(blends))
            plot!(q,
                  [0, 1],
                  [cut_container[c, it] / blend[2], (cut_container[c, it] - blend[1]) / blend[2]],
                  fill=(1, 0.5, :auto),
                  color=(:olivedrab, :dodgerblue)[c],
                  label=("Antarctic IT", "Lunar UFS")[c])
        end

        plot(q, p,
             layout=grid(2, 1),
             size=(500, 1000),
             titlefontsize=12,
             leftmargin=20px)
    end

    return anim
end

gif(anim_doer(), fps=30, "TwoScoreAdmissionsEquilibrium.gif")