etf_simulation.py

from dataclasses import dataclass
from random import gauss
import math

type Eur = float
type TimeUnit = int
type PercentPoint = float


class Security:
    """Abstract common core of Bonds and ETFs."""

    def advance_time(self):
        """Advance maturity by 1 month."""
        raise NotImplementedError()

    def market_value(self, market_rate: PercentPoint) -> Eur:
        """Value on the market, function of the current "market rate".
        """
        raise NotImplementedError()

    def cash_flow(self) -> Eur:
        """Cash flow at the current maturity.

        Cash flows come at coupon time or at expiration, otherwise 0.
        """
        raise NotImplementedError()


_YEAR: TimeUnit = 12


def net_present_value(cash: Eur,
                      discount_rate: PercentPoint,
                      time: TimeUnit) -> Eur:
    """Return the present value of a future cash flow.
    """
    return cash / (1 + discount_rate / 100)**time


@dataclass
class Bond(Security):
    face_value: Eur
    expiration: TimeUnit
    coupon_rate: PercentPoint  # one _YEAR
    time_from_emission: TimeUnit = 0
    coupon_time: TimeUnit = _YEAR // 2

    def __post_init__(self):
        assert self.face_value > 0, \
            f"Bond must have a positive value, not {self.face_value}."
        assert self.expiration > 0, \
            f"Bond must expire in the future, not in {self.expiration}."
        assert self.coupon_rate >= 0, \
            f"Bond must have a non negative coupon rate, not {self.coupon_rate}."
        assert self.time_from_emission >= 0, \
            f"Bond must have a non negative time_from_emission, not {self.time_from_emission}."
        assert self.time_from_emission < self.expiration, \
            "Bond must have a time_from_emission less than expiration."
        assert self.coupon_time > 0 \
            and self.expiration % self.coupon_time == 0, \
            "Coupon time must be positive and a divisor of expiration."

        self.monthly_rate = self.coupon_rate / _YEAR
        self.coupon_value = self.face_value \
            * (self.monthly_rate / 100) * self.coupon_time

    def market_value(self, market_rate: PercentPoint) -> Eur:
        value: Eur = 0.0
        monthly_market_rate = market_rate / _YEAR
        for t in range(self.time_from_emission, self.expiration + 1):
            if t % self.coupon_time == 0:
                value += net_present_value(self.coupon_value,
                                           monthly_market_rate,
                                           t - self.time_from_emission)
            if t == self.expiration:
                value += net_present_value(self.face_value,
                                           monthly_market_rate,
                                           t - self.time_from_emission)
        return value

    def find_face_value(self,
                        market_value: Eur,
                        market_rate: PercentPoint) -> Eur:
        """Return the face value needed to get a specific market value
           in the last time unit, at the specific market rate.

        >>> import math
        >>> c_r = 20.0
        >>> m_r = 2.0
        >>> ex = 24
        >>> target = 1000
        >>> b = Bond(face_value=100, coupon_rate=c_r, expiration=ex)
        >>> f_v = b.find_face_value(target, m_r)
        >>> x = Bond(face_value=f_v, coupon_rate=c_r,
        ...          expiration=ex, time_from_emission=ex-1)
        >>> math.isclose(x.market_value(m_r), target, abs_tol=1e-3)
        True
        """
        def f(x: Eur):
            b = Bond(face_value=x,
                     expiration=self.expiration,
                     time_from_emission=self.expiration-1,
                     coupon_rate=self.coupon_rate)
            return b.market_value(market_rate)

        low = 0.001
        high = 10 * market_value
        while high - low > 1e-3:
            m = (low + high) / 2
            f_m = f(m)
            if f_m < market_value:
                low = m
            else:
                high = m
        return (low + high) / 2

    def cash_flow(self) -> Eur:
        cash = 0.0
        if self.is_expired():
            return 0.0
        if self.time_from_emission > 0 \
           and self.time_from_emission % self.coupon_time == 0:
            cash += self.coupon_value
        if self.time_from_emission == self.expiration:
            cash += self.face_value
        return cash

    def advance_time(self):
        self.time_from_emission += 1

    def is_expired(self) -> bool:
        """Return true when time_from_emission is greater than expiration.
        """
        return self.time_from_emission > self.expiration

    def duration(self, market_rate: PercentPoint) -> float:
        """Return Macaulay duration at current market rate.

        >>> import math
        >>> b = Bond(face_value=100, coupon_rate=20.0, expiration=24)
        >>> math.isclose(b.duration(4), 1.777, abs_tol=0.001)
        True
        """
        values: list[Eur] = []
        year_weights: list[float] = []
        for t in range(self.time_from_emission, self.expiration+1):
            if t > 0 and t % self.coupon_time == 0:
                values.append(net_present_value(self.coupon_value,
                                                market_rate
                                                / _YEAR
                                                * self.coupon_time,
                                                len(values)+1))
                year_weights.append(t / _YEAR)
        values.append(net_present_value(self.face_value,
                                        market_rate / _YEAR * self.coupon_time,
                                        len(values)))
        year_weights.append(self.expiration / _YEAR)
        v = sum(values)
        return sum(y * x/v for x, y in zip(values, year_weights))


def simulate_single_bond(b: Bond, market: list[PercentPoint]):
    total_cf = 0.0
    print(b)
    for t, mr in enumerate(market):
        total_cf += b.cash_flow()
        print(f"Month {t:2d}, the market rate is {mr:1.2f}%. "
              + f"Market value {b.market_value(mr):.2f}, "
              + f"CF: {b.cash_flow():.2f}")
        b.advance_time()
    print(f"Total cash flow: {total_cf:.2f}")


class ETF(Security):
    def __init__(self,
                 roll_cycle: TimeUnit,
                 market_rate: PercentPoint,
                 capital: Eur):
        assert roll_cycle > 0, \
            f"Roll cycle must be positive, not {roll_cycle}."
        assert market_rate >= 0, \
            f"Current market rate cannot be negative, not {roll_cycle}."
        assert capital > 0, \
            "The capital to be invested must be positive."
        self.roll_cycle = roll_cycle
        self.bb = []
        for m in range(self.roll_cycle - 1):
            self.bb.append(Bond(face_value=capital / self.roll_cycle,
                                expiration=self.roll_cycle,
                                coupon_rate=market_rate,
                                time_from_emission=m))
        value = self.market_value(market_rate)
        if value < capital:
            last = self.bb[-1].find_face_value(capital-value, market_rate)
            self.bb.append(Bond(face_value=last,
                                expiration=self.roll_cycle,
                                coupon_rate=market_rate,
                                time_from_emission=self.roll_cycle-1))

    def market_value(self, market_rate: PercentPoint) -> Eur:
        return math.fsum(b.market_value(market_rate) for b in self.bb)

    def advance_time(self):
        for b in self.bb:
            b.advance_time()

    def replace_expired(self, market_rate: PercentPoint):
        """Buy a new Bond for each one is_expired today.

        The new Bond has the same face value
        and will expire at the end of my roll cycle.
        """
        new_bb = []
        for b in self.bb:
            if b.is_expired():
                new_bb.append(Bond(face_value=b.face_value,
                                   expiration=self.roll_cycle,
                                   coupon_rate=market_rate))
            else:
                new_bb.append(b)
        self.bb = new_bb

    def cash_flow(self) -> Eur:
        cf = 0.0
        for b in self.bb:
            cf += b.cash_flow()
            if b.time_from_emission == b.expiration:
                cf -= b.face_value
        return cf

    def average_duration(self, market_rate: PercentPoint) -> float:
        return sum(b.duration(market_rate) for b in self.bb) / len(self.bb)


def simulate_etf(e: ETF, market: list[PercentPoint]):
    total_cf = 0.0
    print(f'ETF composition ({e.market_value(market[0]):.2f}EUR):\n'
          + '\n'.join([str(b) for b in e.bb]))
    for t, mr in enumerate(market):
        total_cf += e.cash_flow()
        print(f"Month {t:2d}, the market rate is {mr:1.2f}%. "
              + f"Market value {e.market_value(mr):.2f}, "
              + f"Average duration {e.average_duration(mr):.2f} years, "
              + f"Cash flow: {e.cash_flow():.2f}EUR")
        e.advance_time()
        e.replace_expired(mr)
    print(f"Total cash flow: {total_cf:.2f}EUR")
    print(f"Total value: {total_cf+e.market_value(market[-1]):.2f}EUR")


if __name__ == "__main__":
    # Market rate starts at 2%,
    # then each month it changes with a delta with mean 0, stdev 0.25
    market = [2.0]
    for i in range(1, 5*_YEAR + 1):
        market.append(max(market[i-1] - gauss(0, .25), 0.))

    roll_cycle = 3*_YEAR

    b = Bond(face_value=1000*roll_cycle,
             expiration=roll_cycle,
             coupon_rate=market[0])
    simulate_single_bond(b, market)

    etf = ETF(3*_YEAR, market[0], 36000)
    simulate_etf(etf, market)