Tracking international bond index - GDX input#

BondIndex.py BondIndex.gdx

"""
## 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
from gamspy import Equation
from gamspy import Model
from gamspy import ModelStatus
from gamspy import Parameter
from gamspy import Set
from gamspy import Sum
from gamspy import Variable


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )
    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)

    ResultHandle = open("BondIndex.csv", "w", encoding="UTF-8")
    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)}'
    )
    ResultHandle.close()

    import math

    assert math.isclose(BondIndex.objective_value, 23.206448, rel_tol=0.001)


if __name__ == "__main__":
    main()