Parts Supply Problem w/ 10 Types w/ Random p(i) (PS10_S_MN)#

ps10_s_mn.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_ps10_s_mn.html
## LICENSETYPE: Demo
## MODELTYPE: NLP
## KEYWORDS: nonlinear programming, contract theory, principal-agent problem, adverse selection, parts supply problem


Parts Supply Problem w/ 10 Types w/ Random p(i) (PS10_S_MN)

Hideo Hashimoto, Kojun Hamada, and Nobuhiro Hosoe, "A Numerical Approach
to the Contract Theory: the Case of Adverse Selection", GRIPS Discussion
Paper 11-27, National Graduate Institute for Policy Studies, Tokyo, Japan,
March 2012.
"""

from __future__ import annotations

import os
import time

from gamspy import Alias
from gamspy import Card
from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Number
from gamspy import Options
from gamspy import Ord
from gamspy import Parameter
from gamspy import Problem
from gamspy import Sense
from gamspy import Set
from gamspy import Sum
from gamspy import Variable
from gamspy.math import Round
from gamspy.math import uniform


def main():
    start = time.time()

    # Decreased no. of draw from 1001 to 11 for convenience
    # Otherwise, it takes a lot of time
    NUM_DRAWS = 11

    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )

    # Sets
    i = Set(
        m,
        name="i",
        records=[str(i) for i in range(10)],
        description="type of supplier",
    )

    t = Set(
        m,
        name="t",
        records=[str(i) for i in range(1, NUM_DRAWS)],
        description="no. of Monte-Carlo draws",
    )

    j = Alias(m, name="j", alias_with=i)

    # Parameters
    theta = Parameter(m, name="theta", domain=i, description="efficiency")
    pt = Parameter(
        m, name="pt", domain=[i, t], description="probability of type"
    )
    p = Parameter(m, name="p", domain=i, description="probability of type")

    theta[i] = Ord(i) / Card(i)

    # Generating probability
    for tt, _ in t.records.itertuples(index=False):
        pt[i, tt] = uniform(0, 1)
    pt[i, t] = pt[i, t] / Sum(j, pt[j, t])

    # Parameters
    F = Parameter(
        m,
        name="F",
        domain=[i, t],
        description="cumulative probability (Itho p. 42)",
    )
    noMHRC0 = Parameter(
        m,
        name="noMHRC0",
        domain=[i, t],
        description="no MHRC combination between i and i-1",
    )
    # (MHRC: monotone hazard rate condition)
    noMHRC = Parameter(
        m, name="noMHRC", domain=t, description=">=1: no MHRC case"
    )

    F[i, t] = Sum(j.where[Ord(j) <= Ord(i)], pt[j, t])
    noMHRC0[i, t].where[Ord(i) < Card(i)] = Number(1).where[
        F[i, t] / pt[i.lead(1), t] < F[i.lag(1), t] / pt[i, t]
    ]
    noMHRC[t].where[Sum(i, noMHRC0[i, t]) >= 1] = 1

    ru = Parameter(m, name="ru", records=0, description="reservation utility")

    # Definition of Primal/Dual Variables
    x = Variable(m, name="x", type="positive", domain=i, description="quality")
    b = Variable(
        m, name="b", type="positive", domain=i, description="maker's revenue"
    )
    w = Variable(m, name="w", type="positive", domain=i, description="price")

    # Equations
    rev = Equation(
        m,
        name="rev",
        domain=i,
        description="maker's revenue function",
    )
    pc = Equation(
        m,
        name="pc",
        domain=i,
        description="participation constraint",
    )
    licd = Equation(
        m,
        name="licd",
        domain=i,
        description="incentive compatibility constraint",
    )
    licu = Equation(
        m,
        name="licu",
        domain=i,
        description="incentive compatibility constraint",
    )
    ic = Equation(
        m,
        name="ic",
        domain=[i, j],
        description="global incentive compatibility constraint",
    )
    mn = Equation(
        m,
        name="mn",
        domain=i,
        description="monotonicity constraint",
    )

    obj = Sum(i, p[i] * (b[i] - w[i]))  # Objective Function; maker's utility

    rev[i] = b[i] == x[i] ** (0.5)

    pc[i] = w[i] - theta[i] * x[i] >= ru

    licd[i] = w[i] - theta[i] * x[i] >= w[i.lead(1)] - theta[i] * x[i.lead(1)]

    licu[i] = w[i] - theta[i] * x[i] >= w[i.lag(1)] - theta[i] * x[i.lag(1)]

    ic[i, j] = w[i] - theta[i] * x[i] >= w[j] - theta[i] * x[j]

    mn[i] = x[i] >= x[i.lead(1)]

    # Setting Lower Bounds on Variables to Avoid Division by Zero
    x.lo[i] = 0.0001

    # Models
    SB_lic = Model(
        m,
        name="SB_lic",
        equations=[rev, pc, licd],
        problem=Problem.NLP,
        sense=Sense.MAX,
        objective=obj,
    )
    SB_lic2 = Model(
        m,
        name="SB_lic2",
        equations=[rev, pc, licd, mn],
        problem=Problem.NLP,
        sense=Sense.MAX,
        objective=obj,
    )

    # Parameters
    Util_lic = Parameter(
        m, name="Util_lic", domain=t, description="util solved w/o MN"
    )
    Util_lic2 = Parameter(
        m, name="Util_lic2", domain=t, description="util solved w/ MN"
    )
    Util_gap = Parameter(
        m,
        name="Util_gap",
        domain=t,
        description="gap between these two util",
    )
    x_lic = Parameter(
        m, name="x_lic", domain=[i, t], description="x solved in w/o MN"
    )
    x_lic2 = Parameter(
        m, name="x_lic2", domain=[i, t], description="x solved in w/ MN"
    )
    MN_lic = Parameter(
        m,
        name="MN_lic",
        domain=t,
        description="monotonicity of x solved w/o MN",
    )
    MN_lic2 = Parameter(
        m,
        name="MN_lic2",
        domain=t,
        description="monotonicity of x solved w/ MN",
    )

    for tt, _ in t.records.itertuples(index=False):
        p[i] = pt[i, tt]

        #  Solving the model w/o MN
        SB_lic.solve(options=Options(solver_link_type=5))

        Util_lic[tt] = SB_lic.objective_value
        x_lic[i, tt] = x.l[i]
        MN_lic[tt] = Sum(
            i, Number(1).where[Round(x.l[i], 10) < Round(x.l[i.lead(1)], 10)]
        )

        #  Solving the model w/ MN
        SB_lic2.solve(options=Options(solver_link_type=5))

        Util_lic2[tt] = SB_lic2.objective_value
        x_lic2[i, tt] = x.l[i]
        MN_lic2[tt] = Sum(
            i, Number(1).where[Round(x.l[i], 10) < Round(x.l[i.lead(1)], 10)]
        )

    Util_gap[t] = Number(1).where[
        Round(Util_lic[t], 10) != Round(Util_lic2[t], 10)
    ]

    # Computing probability that MHRC and MN holds.
    p_noMHRC = Parameter(
        m, name="p_noMHRC", description="no MHRC case          [%]"
    )
    p_noMN_lic = Parameter(
        m, name="p_noMN_lic", description="no MN case            [%]"
    )
    p_Util_gap = Parameter(
        m, name="p_Util_gap", description="no util-equality case [%]"
    )

    p_noMHRC[...] = Sum(t.where[noMHRC[t] > 0], 1) / Card(t) * 100
    p_noMN_lic[...] = Sum(t.where[MN_lic[t] > 0], 1) / Card(t) * 100
    p_Util_gap[...] = Sum(t.where[Util_gap[t] > 0], 1) / Card(t) * 100

    print(f"no MHRC case: {p_noMHRC.records.value[0]}%")
    print(f"no MN case: {p_noMN_lic.records.value[0]}%")
    print(f"no util-equality case: {p_Util_gap.records.value[0]}%")
    print(f"Time Elapsed: {round(time.time()-start, 2)}s")


if __name__ == "__main__":
    main()