# 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
## 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]

# 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

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

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()
```