# Copyright (C) 2023 Luigi Ballabio
#
# This file is released under the terms of the 3-Clause BSD License
# (see https://opensource.org/license/bsd-3-clause/)
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER “AS IS” AND ANY
# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
# OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


# Required modules

import QuantLib as ql
import csv
from matplotlib import pyplot as plt
from matplotlib.dates import YearLocator
from matplotlib.ticker import PercentFormatter


def new_plot():
    f = plt.figure(figsize=(8, 5))
    ax = f.add_subplot(1, 1, 1)

    ax.xaxis.grid(True, "major", color="lightgray")
    ax.yaxis.grid(True, "major", color="lightgray")
    ax.xaxis.set_major_locator(YearLocator(10))

    return f, ax


# Global evaluation date

today = ql.Date(21, ql.September, 2021)
ql.Settings.instance().evaluationDate = today


# Read bond data from a CSV file


def from_iso(date):
    return ql.Date(date, "%Y-%m-%d")


def apply_types(row):
    return [
        from_iso(row[0]),
        from_iso(row[1]),
        float(row[2]),
        float(row[3]),
    ]


with open("2023-09.csv") as f:
    reader = csv.reader(f)
    next(reader)  # skip header
    data = [apply_types(row) for row in reader]
data.sort(key=lambda r: r[1])


# Create bonds and helpers

helpers = []
bonds = []
for start, maturity, coupon, price in data:
    schedule = ql.Schedule(
        start,
        maturity,
        ql.Period(1, ql.Years),
        ql.TARGET(),
        ql.ModifiedFollowing,
        ql.ModifiedFollowing,
        ql.DateGeneration.Backward,
        False,
    )
    bond = ql.FixedRateBond(
        3,
        100.0,
        schedule,
        [coupon / 100.0],
        ql.Actual360(),
        ql.ModifiedFollowing,
    )
    bonds.append(bond)
    helpers.append(ql.BondHelper(ql.QuoteHandle(ql.SimpleQuote(price)), bond))

discount_curve = ql.RelinkableYieldTermStructureHandle()
bond_engine = ql.DiscountingBondEngine(discount_curve)
for b in bonds:
    b.setPricingEngine(bond_engine)

# Fit to a few curve models

methods = {
    "Nelson/Siegel": ql.NelsonSiegelFitting(),
    "Exp. splines": ql.ExponentialSplinesFitting(True),
    "B splines": ql.CubicBSplinesFitting(
        [
            -30.0,
            -20.0,
            0.0,
            5.0,
            10.0,
            15.0,
            20.0,
            25.0,
            30.0,
            40.0,
            50.0,
        ],
        True,
    ),
    "Svensson": ql.SvenssonFitting(),
}

tolerance = 1e-8
max_iterations = 5000
day_count = ql.Actual360()

curves = {
    tag: ql.FittedBondDiscountCurve(
        today,
        helpers,
        day_count,
        methods[tag],
        tolerance,
        max_iterations,
    )
    for tag in methods
}


# Plot the curves

f, ax = new_plot()
ax.yaxis.set_major_formatter(PercentFormatter(1.0))
styles = iter(["-", "--", ":", "-."])

dates = [today + ql.Period(i, ql.Months) for i in range(12 * 30 + 1)]

for tag in curves:
    rates = [curves[tag].zeroRate(d, day_count, ql.Continuous).rate() for d in dates]
    ax.plot_date(
        [d.to_date() for d in dates],
        rates,
        next(styles),
        label=tag,
    )
ax.legend(loc="best")
f.savefig("2023-09.1.pdf")

# Plot the prices

quoted_prices = [row[-1] for row in data]


def prices(tag):
    discount_curve.linkTo(curves[tag])
    return [b.cleanPrice() for b in bonds]


def errors(tag):
    return [q - p for p, q in zip(prices(tag), quoted_prices)]


f, ax = new_plot()
maturities = [r[1].to_date() for r in data]
markers = iter([".", "+", "x", "1"])
for tag in curves:
    ax.plot_date(
        maturities,
        errors(tag),
        next(markers),
        label=tag,
    )
ax.legend(loc="best")
f.savefig("2023-09.2.pdf")


f, ax = new_plot()
ps = prices("Svensson")
qs = quoted_prices
ax.plot_date(maturities, qs, "P", label="quoted")
ax.plot_date(maturities, ps, "o", label="Svensson")
ax.legend(loc="best")
for m, p, q in zip(maturities, ps, qs):
    ax.plot_date([m, m], [p, q], "-", color="grey")
f.savefig("2023-09.3.pdf")


# Another point of view

quoted_yields = [
    b.bondYield(
        ql.BondPrice(p, ql.BondPrice.Clean),
        day_count,
        ql.Compounded,
        ql.Annual,
    )
    for b, p in zip(bonds, quoted_prices)
]


def yields(tag):
    discount_curve.linkTo(curves[tag])
    return [b.bondYield(day_count, ql.Compounded, ql.Annual) for b in bonds]


f, ax = new_plot()
ax.yaxis.set_major_formatter(PercentFormatter(1.0))
ys = yields("Svensson")
qys = quoted_yields
ax.plot_date(maturities, qys, ".", label="quoted")
ax.plot_date(maturities, ys, "x", label="Svensson")
ax.legend(loc="best")
f.savefig("2023-09.4.pdf")


# Remove outliers

ys = yields("Svensson")
qys = quoted_yields

filtered_helpers = [h for h, y1, y2 in zip(helpers, ys, qys) if abs(y1 - y2) < 0.005]

curves["Svensson (new)"] = ql.FittedBondDiscountCurve(
    today,
    filtered_helpers,
    day_count,
    ql.SvenssonFitting(),
    tolerance,
    max_iterations,
)

# Improved results

f, ax = new_plot()
ax.yaxis.set_major_formatter(PercentFormatter(1.0))
ys = yields("Svensson")
ys2 = yields("Svensson (new)")
qys = quoted_yields
ax.plot_date(maturities, qys, ".", label="quoted")
ax.plot_date(maturities, ys2, "x", label="Svensson (new)")
ax.legend(loc="best")
f.savefig("2023-09.5.pdf")


f, ax = new_plot()
ps = prices("Svensson (new)")
qs = quoted_prices
ax.plot_date(maturities, qs, "P", label="quoted")
ax.plot_date(maturities, ps, "o", label="Svensson (new)")
ax.legend(loc="best")
for m, p, q in zip(maturities, ps, qs):
    ax.plot_date([m, m], [p, q], "-", color="grey")
f.savefig("2023-09.6.pdf")