Immunization models#

Immunization.py

"""
## GAMSSOURCE: https://www.gams.com/latest/finlib_ml/libhtml/finlib_Immunization.html
## LICENSETYPE: Demo
## MODELTYPE: LP


Immunization models

Immunization.gms: Immunization models.
Consiglio, Nielsen and Zenios.
PRACTICAL FINANCIAL OPTIMIZATION: A Library of GAMS Models, Section 4.4
Last modified: Apr 2008.
"""

from __future__ import annotations

import os

import numpy as np
import pandas as pd

import gamspy.math as gams_math
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 Ord
from gamspy import Parameter
from gamspy import Sense
from gamspy import Set
from gamspy import Sum
from gamspy import Variable
from gamspy.math import sqr


def main():
    # Define container
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )

    # Bond data. Prices, coupons and maturities from the Danish market
    bond_data_recs = pd.DataFrame(
        np.array([
            [112.35, 2006, 8],
            [105.33, 2003, 8],
            [111.25, 2007, 7],
            [107.30, 2004, 7],
            [107.62, 2011, 6],
            [106.68, 2009, 6],
            [101.93, 2002, 6],
            [101.30, 2005, 5],
            [101.61, 2003, 5],
            [100.06, 2002, 4],
        ]),
        columns=["Price", "Maturity", "Coupon"],
        index=[
            "DS-8-06",
            "DS-8-03",
            "DS-7-07",
            "DS-7-04",
            "DS-6-11",
            "DS-6-09",
            "DS-6-02",
            "DS-5-05",
            "DS-5-03",
            "DS-4-02",
        ],
    )

    bond_data_recs = bond_data_recs.reset_index().melt(
        id_vars="index", var_name="Category", value_name="Value"
    )

    # SETS #
    Time = Set(
        m,
        name="Time",
        records=[str(i) for i in range(2001, 2012)],
        description="Time periods",
    )

    Bonds = Set(
        m,
        name="Bonds",
        records=[
            "DS-8-06",
            "DS-8-03",
            "DS-7-07",
            "DS-7-04",
            "DS-6-11",
            "DS-6-09",
            "DS-6-02",
            "DS-5-05",
            "DS-5-03",
            "DS-4-02",
        ],
        description="Bonds universe",
    )

    # ALIASES #
    t = Alias(m, name="t", alias_with=Time)
    i = Alias(m, name="i", alias_with=Bonds)

    # ALIAS (Time, t, t1, t2)

    # SCALARS #
    Now = Parameter(m, name="Now", description="Current year")
    Horizon = Parameter(m, name="Horizon", description="End of the Horizon")

    Now[...] = 2001
    Horizon[...] = Card(t) - 1

    # PARAMETER #
    tau = Parameter(m, name="tau", domain=t, description="Time in years")

    # Note: time starts from 0
    tau[t] = Ord(t) - 1

    Coupon = Parameter(m, name="Coupon", domain=i, description="Coupons")
    Maturity = Parameter(
        m, name="Maturity", domain=i, description="Maturities"
    )
    F = Parameter(m, name="F", domain=[t, i], description="Cashflows")
    BondData = Parameter(
        m,
        name="BondData",
        domain=[i, "*"],
        records=bond_data_recs,
        description="Bonds data",
    )
    Liability = Parameter(
        m, name="Liability", domain=t, description="Stream of liabilities"
    )

    # Copy/transform data. Note division by 100 to get unit data, and
    # subtraction of "Now" from Maturity date (so consistent with tau):

    Coupon[i] = BondData[i, "Coupon"] / 100
    Maturity[i] = BondData[i, "Maturity"] - Now

    # Calculate the ex-coupon cashflow of Bond i in year t:

    F[t, i] = (
        Number(1).where[tau[t] == Maturity[i]]
        + Coupon[i].where[(tau[t] <= Maturity[i]) & (tau[t] > 0)]
    )

    Liability.setRecords(
        np.array([
            0,
            80000,
            100000,
            110000,
            120000,
            140000,
            120000,
            90000,
            50000,
            75000,
            150000,
        ])
    )

    r = Parameter(
        m,
        name="r",
        domain=t,
        records=np.array([
            0,
            0.0422,
            0.0440,
            0.0450,
            0.0466,
            0.0480,
            0.0482,
            0.0485,
            0.0488,
            0.0491,
            0.0493,
        ]),
        description="spot rates",
    )
    y = Parameter(
        m,
        name="y",
        domain=i,
        records=np.array([
            0.0501,
            0.0500,
            0.0469,
            0.0426,
            0.0489,
            0.0485,
            0.0392,
            0.0453,
            0.0406,
            0.0386,
        ]),
        description="yield rates",
    )

    # The following are the Present value, Fischer-Weil duration (D^FW)
    # and Convexity (Q_i), for both the bonds and the liabilities:
    # Present value, Fisher & Weil duration, and convexity for
    # the bonds.
    PV = Parameter(
        m, name="PV", domain=i, description="Present value of assets"
    )
    Dur = Parameter(m, name="Dur", domain=i, description="Duration of assets")
    Conv = Parameter(
        m, name="Conv", domain=i, description="Convexity of assets"
    )

    # Present value, Fisher & Weil duration, and convexity for
    # the liability.
    PV_Liab = Parameter(
        m, name="PV_Liab", description="Present value of liability"
    )
    Dur_Liab = Parameter(
        m, name="Dur_Liab", description="Duration of liability"
    )
    Conv_Liab = Parameter(
        m, name="Conv_Liab", description="Convexity of liability"
    )

    PV[i] = Sum(t, F[t, i] * gams_math.exp(-r[t] * tau[t]))

    Dur[i] = (1.0 / PV[i]) * Sum(
        t, tau[t] * F[t, i] * gams_math.exp(-r[t] * tau[t])
    )

    Conv[i] = (1.0 / PV[i]) * Sum(
        t, sqr(tau[t]) * F[t, i] * gams_math.exp(-r[t] * tau[t])
    )

    print("PV: \n", PV.records)
    print("Dur: \n", Dur.records)
    print("Conv: \n", Conv.records)

    # Calculate the corresponding amounts for Liabilities. Use its PV as its
    # "price".

    PV_Liab[...] = Sum(t, Liability[t] * gams_math.exp(-r[t] * tau[t]))

    Dur_Liab[...] = (1.0 / PV_Liab) * Sum(
        t, tau[t] * Liability[t] * gams_math.exp(-r[t] * tau[t])
    )

    Conv_Liab[...] = (1.0 / PV_Liab) * Sum(
        t, sqr(tau[t]) * Liability[t] * gams_math.exp(-r[t] * tau[t])
    )

    print("PV_Liab: ", PV_Liab.records.value[0])
    print("Dur_Liab: ", Dur_Liab.records.value[0])
    print("Conv_Liab: ", Conv_Liab.records.value[0])

    # Build a sequence of increasingly sophisticated immunuzation models.

    # VARIABLES #
    x = Variable(
        m,
        name="x",
        type="positive",
        domain=i,
        description="Holdings of bonds (amount of face value)",
    )

    # EQUATIONS #
    PresentValueMatch = Equation(
        m,
        name="PresentValueMatch",
        description=(
            "Equation matching the present value of asset and liability"
        ),
    )
    DurationMatch = Equation(
        m,
        name="DurationMatch",
        description="Equation matching the duration of asset and liability",
    )
    ConvexityMatch = Equation(
        m,
        name="ConvexityMatch",
        description="Equation matching the convexity of asset and liability",
    )

    z = Sum(i, Dur[i] * PV[i] * y[i] * x[i]) / (PV_Liab * Dur_Liab)

    PresentValueMatch[...] = Sum(i, PV[i] * x[i]) == PV_Liab

    DurationMatch[...] = Sum(i, Dur[i] * PV[i] * x[i]) == PV_Liab * Dur_Liab

    ConvexityMatch[...] = Sum(i, Conv[i] * PV[i] * x[i]) >= PV_Liab * Conv_Liab

    ImmunizationOne = Model(
        m,
        name="ImmunizationOne",
        equations=[PresentValueMatch, DurationMatch],
        problem="LP",
        sense=Sense.MAX,
        objective=z,
    )
    ImmunizationOne.solve()

    Convexity = Parameter(m, name="Convexity")
    Convexity[...] = (1.0 / PV_Liab) * Sum(i, Conv[i] * PV[i] * x.l[i])
    x_results = []

    x_results.append(x.records.level.tolist())
    print("Convexity: ", Convexity.records.value[0])
    print("Conv_Liab: ", Conv_Liab.records.value[0])

    ImmunizationTwo = Model(
        m,
        name="ImmunizationTwo",
        equations=[PresentValueMatch, DurationMatch, ConvexityMatch],
        problem="LP",
        sense=Sense.MAX,
        objective=z,
    )
    ImmunizationTwo.solve()

    DurationMatch.l[...] = DurationMatch.l / PV_Liab

    ConvexityMatch.l[...] = ConvexityMatch.l / PV_Liab

    x_results.append(x.records.level.tolist())
    print("PresentValueMatch: ", PresentValueMatch.records.level[0])
    print("DurationMatch: ", DurationMatch.records.level[0])
    print("ConvexityMatch: ", ConvexityMatch.records.level[0])

    z = (1.0 / PV_Liab) * Sum(i, Conv[i] * PV[i] * x[i])

    ImmunizationThree = Model(
        m,
        name="ImmunizationThree",
        equations=[PresentValueMatch, DurationMatch],
        problem="LP",
        sense=Sense.MIN,
        objective=z,
    )
    ImmunizationThree.solve()

    x_results.append(x.records.level.tolist())
    x_results = pd.DataFrame(
        np.array(x_results).T,
        columns=["Model1", "Model2", "Model3"],
        index=Bonds.records.uni,
    )
    print("Objective Function Value: ", ImmunizationThree.objective_value)

    print(x_results)


if __name__ == "__main__":
    main()