"""
## GAMSSOURCE: https://www.gams.com/latest/finlib_ml/libhtml/finlib_BondIndex.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## DATAFILES: BondIndex.gdx
Tracking international bond index - GDX input
* BondIndexGDX.gms: Tracking international bond index - GDX input.
* Consiglio, Nielsen and Zenios.
* PRACTICAL FINANCIAL OPTIMIZATION: A Library of GAMS Models, Section 8.2
* Last modified: Apr 2008.
"""
from __future__ import annotations
import os
from pathlib import Path
import numpy as np
from gamspy import (
Container,
Equation,
Model,
ModelStatus,
Parameter,
Set,
Sum,
Variable,
)
def main():
m = Container(
system_directory=os.getenv("SYSTEM_DIRECTORY", None),
)
m.read(
load_from=str(Path(__file__).parent.absolute()) + "/BondIndex.gdx",
symbol_names=[
"Bonds",
"Currencies",
"Scenarios",
"j",
"i",
"l",
"SS",
"JxI",
"ExchangeRates0",
"ExchangeRates1",
"Price0",
"Price1",
"InitialHoldings",
"Accruals0",
"Accruals1",
"Outstanding",
"ReinvestmentRate",
"IndexReturns",
"pr",
],
)
# SETS #
j, i, l = m.getSymbols(["j", "i", "l"])
JxI = m.getSymbols(["JxI"])[0]
i_recs = i.records.element_text.tolist()
# PARAMETERS #
(
ExchangeRates0,
ExchangeRates1,
Price0,
Price1,
InitialHoldings,
Accruals0,
Accruals1,
Outstanding,
ReinvestmentRate,
IndexReturns,
pr,
) = m.getSymbols(
[
"ExchangeRates0",
"ExchangeRates1",
"Price0",
"Price1",
"InitialHoldings",
"Accruals0",
"Accruals1",
"Outstanding",
"ReinvestmentRate",
"IndexReturns",
"pr",
]
)
# SCALARS #
TrnCstB = Parameter(
m,
name="TrnCstB",
records=0.0025,
description="Buying transaction costs",
)
TrnCstS = Parameter(
m,
name="TrnCstS",
records=0.0015,
description="Selling transaction costs",
)
CashInfusion = Parameter(
m,
name="CashInfusion",
records=100000,
description="Available budget infused",
)
EpsTolerance = Parameter(
m, name="EpsTolerance", records=0.10, description="Tolerance"
)
UpprBnd = Parameter(
m,
name="UpprBnd",
records=0.1,
description="Maximum holding percentage for each bond",
)
CHFtrade = Parameter(
m,
name="CHFtrade",
records=15000000,
description="Maximum trading value allowed for Swiss bonds (in CHF)",
)
HoldVal = Parameter(
m, name="HoldVal", description="Value of the initial holdings"
)
InitAccrCash = Parameter(
m,
name="InitAccrCash",
description="Accrued cash originated by the initial holdings",
)
InitVal = Parameter(
m, name="InitVal", description="Initial portfolio value"
)
# Calculate initial portfolio value
HoldVal[...] = Sum(
JxI[j, i], ExchangeRates0[j] * InitialHoldings[i] * Price0[i]
)
InitAccrCash[...] = Sum(
JxI[j, i], ExchangeRates0[j] * InitialHoldings[i] * Accruals0[i]
)
InitVal[...] = CashInfusion + InitAccrCash + HoldVal
# VARIABLES #
X0 = Variable(
m,
name="X0",
type="positive",
domain="*",
description="Face value bought today.",
)
Y0 = Variable(
m,
name="Y0",
type="positive",
domain="*",
description="Face value sold today.",
)
Z0 = Variable(
m,
name="Z0",
type="positive",
domain="*",
description="Face value hold today for the next period.",
)
Cash = Variable(
m,
name="Cash",
type="positive",
description=(
"Amount of cash resulting from trading (sell and buy) today."
),
)
FinalCash = Variable(
m,
name="FinalCash",
type="positive",
domain=l,
description="Amount of cash resulting from portfolio liquidation",
)
# Set the upper bound on the holdings
for jj, ii, _ in JxI.records.itertuples(index=False):
Z0.up[ii] = InitVal / ExchangeRates0[jj] / Price0[ii] * UpprBnd
# Set the limit on trading (sell or buy)
# CHF bonds for liquidity reasons
X0.up[i].where[JxI["CHF", i]] = CHFtrade / Price0[i]
Y0.up[i].where[JxI["CHF", i]] = CHFtrade / Price0[i]
# EQUATIONS #
CashInventoryCon = Equation(
m,
name="CashInventoryCon",
type="regular",
description="Cash balance equation today.",
)
FinalCashCon = Equation(
m,
name="FinalCashCon",
type="regular",
domain=l,
description="Cash balance equations at the end of first stage.",
)
InventoryCon = Equation(
m,
name="InventoryCon",
type="regular",
domain=i,
description="Constraints defining the asset inventory balance",
)
MADCon = Equation(
m,
name="MADCon",
type="regular",
domain=l,
description="MAD constraints",
)
Obj = 1000 * Sum(l, pr[l] * (FinalCash[l] / InitVal - 1))
CashInventoryCon[...] = (
CashInfusion
+ Sum(JxI[j, i], ExchangeRates0[j] * Y0[i] * Price0[i] * (1 - TrnCstS))
== Sum(
JxI[j, i], ExchangeRates0[j] * X0[i] * Price0[i] * (1 + TrnCstB)
)
+ Cash
)
FinalCashCon[l] = (
Sum(JxI[j, i], ExchangeRates1[j, l] * Accruals1[i, l] * Z0[i])
+ ReinvestmentRate[l] * Cash
+ Sum(
JxI[j, i],
ExchangeRates1[j, l]
* Z0[i]
* Outstanding[i, l]
* Price1[i, l]
* (1 - TrnCstS),
)
== FinalCash[l]
)
InventoryCon[i] = InitialHoldings[i] + X0[i] == Y0[i] + Z0[i]
MADCon[l] = (FinalCash[l] / InitVal - 1) - IndexReturns[l] >= -EpsTolerance
BondIndex = Model(
m,
name="BondIndex",
equations=m.getEquations(),
problem="LP",
sense="MAX",
objective=Obj,
)
# Find a feasible EpsTolerance
low = Parameter(m, name="low", description="Lower bisection value")
high = Parameter(m, name="high", description="Upper bisection value")
high[...] = -np.inf
low[...] = 0
EpsTolerance[...] = 0.01
while high.records.value[0] <= 0:
BondIndex.solve()
if BondIndex.status in [
ModelStatus.OptimalGlobal,
ModelStatus.OptimalLocal,
]:
high[...] = EpsTolerance
else:
EpsTolerance[...] = 2 * EpsTolerance
# Find a small feasible EpsTolerance via bisection
while True:
EpsTolerance[...] = (low.records.value[0] + high.records.value[0]) / 2
BondIndex.solve()
if BondIndex.status in [
ModelStatus.OptimalGlobal,
ModelStatus.OptimalLocal,
]:
high[...] = EpsTolerance
else:
low[...] = EpsTolerance
if (
(high.records.value[0] - low.records.value[0]) < 0.005
) and BondIndex.status in [
ModelStatus.OptimalGlobal,
ModelStatus.OptimalLocal,
]:
break
CurrentValue = Parameter(
m, name="CurrentValue", domain=i, description="Holdings in USD"
)
CurrentValue[i] = Sum(JxI[j, i], ExchangeRates0[j]) * Price0[i] * Z0.l[i]
ColHeaders = Set(
m,
name="ColHeaders",
records=["FaceValue", "USDValue", "Percent"],
description="Column headers",
)
SummaryReport = Parameter(
m, name="SummaryReport", domain=["*", ColHeaders]
)
SummaryReport[i, "FaceValue"] = Z0.l[i]
SummaryReport[i, "USDValue"] = CurrentValue[i]
SummaryReport[i, "Percent"] = CurrentValue[i] / InitVal * 100
SummaryReport["Total", ColHeaders] = Sum(i, SummaryReport[i, ColHeaders])
print("Summary Report: \n", SummaryReport.pivot().round(3))
print("\nEpsTolerance: \n", round(EpsTolerance.records.value[0], 3))
print(
"\nObjective Function Value: \n", round(BondIndex.objective_value, 3)
)
print("\nInitVal: \n", round(InitVal.records.value[0], 3))
print(CurrentValue.records)
with open("BondIndex.csv", "w", encoding="UTF-8") as ResultHandle:
ResultHandle.write(
f'"Objective Function", {round(BondIndex.objective_value,3)}\n'
)
ResultHandle.write(
f'"Final Epsilon", {round(EpsTolerance.records.value[0],3)}\n'
)
ResultHandle.write(
'"Initial Portfolio Value in USD",'
f" {round(InitVal.records.value[0],3)}\n"
)
ResultHandle.write("\n")
ResultHandle.write(
'"Bond number","CUSIP code","Holdings in unit of face value","Holdings'
' in USD","Percentage of the portfolio value"\n'
)
for ii in Z0.records.itertuples():
if ii.level == 0:
continue
cv_value = CurrentValue.records.loc[
CurrentValue.records["i"] == ii.uni, "value"
].values[0]
ResultHandle.write(
f'"{ii.uni}","{i_recs[ii.Index]}",{round(ii.level,3)},{round(cv_value,2)},{round(100*(cv_value/InitVal.records.value[0]),2)}\n'
)
ResultHandle.write(
'"Cash in US'
f' dollar",{Cash.records.level[0]},{((Cash.records.level[0]/InitVal.records.value[0])*100)}'
)
import math
assert math.isclose(BondIndex.objective_value, 23.206448, rel_tol=0.001)
if __name__ == "__main__":
main()
