Reading the FGN Sovereign Yield Curve: An Exploratory and Inferential Analysis of Daily Bond Yields, January – May 2026

MMBA-8 Data Analytics II — Capstone Case Study 1

Author

Seyi Dutire

Published

May 25, 2026

1 Executive Summary

Business problem. As a fixed income sales trader at a Nigerian inter-dealer brokerage, I need a defensible cross-sectional read of the FGN sovereign curve to inform daily client quotes, relative-value pitches and hedge construction.

Data. I captured 1,212 daily closing yields — the 1-year NTB plus 11 FGN benchmarks (2027 to 2053) — across 101 trading days from 5 January to 25 May 2026, sourced from my desk’s live run-sheet.

Findings. (1) The curve is humped and inverted at the ultra-long end: the belly trades ~190 bp above the 50s/53s (ANOVA F = 530, p ≈ 0, η² = 0.57). (2) The 1-year NTB rallied 135 bp January-to-May (Welch t = 9.73, p < 0.001, Cohen’s d ≈ 3.0) — the market pricing in the easing cycle. (3) Curve co-movement is high (mean Pearson r ≈ 0.78) but the NTB and the 37s are nearly decoupled (r = 0.33) — the genuine relative-value opening. (4) OLS of the 30s on the NTB and the 53s gives yield-betas of 0.62 and 0.66 (R² ≈ 0.81) — ready-to-use hedge ratios.

Recommendation. The 37s is the standout anomaly (+61 bp while every neighbour rallied). Pitch a long 32s / short 37s relative-value trade to PFA clients on Monday, sized off the regression betas.

2 Professional Disclosure

2.1 My role and the organisation

I am a Fixed Income Sales Trader and Broker at a leading Nigerian inter-dealer brokerage (IDB) firm. My desk intermediates trades in FGN bonds, NTBs and OMO Bills between commercial banks, pension fund administrators (PFAs), asset managers and proprietary traders. I publish daily indicative runs, voice trades on real-money client orders, and write a market-recap note circulated to the buy-side every evening. The dataset analysed here is the same closing-level run-sheet I produce as part of my normal workflow — it is genuinely the data my P&L depends on.

2.2 Why these five techniques fit my day-to-day work

Exploratory Data Analysis (EDA, Adi 2026, Ch. 4). Before I quote a level to a client I run a quick distributional check — where is this bond trading relative to its 1-month range? Are there missing prints I should flag? Are any moves outliers vs the rest of the curve? EDA formalises this intuition.

Data Visualisation (Adi 2026, Ch. 5). Every market-close note I publish contains a yield-curve chart, a heatmap of weekly moves, and a steepness time series. Choosing the right visual is how I get a portfolio manager to act in 30 seconds rather than 5 minutes.

Hypothesis Testing (Adi 2026, Ch. 6). When a client asks “has the 1-year NTB really moved lower since CBN’s last meeting?” I need to answer with statistical rigour, not a chart impression. A t-test on first-month vs last-month yields turns a hand-wave into a defendable claim.

Correlation Analysis (Adi 2026, Ch. 8). Cross-bond correlation is the bread and butter of relative-value (RV) trading. The 32s–34s pair, the 27s–29s pair, and the steepener (NTB vs 53s) are all positions whose risk is the correlation. A correlation matrix tells me which RV trades are real and which are just curve-level noise.

Linear Regression (Adi 2026, Ch. 9 — OLS). Regressing a target bond’s yield on a small set of curve points gives me an empirical hedge ratio — how many units of the 1yr NTB and the 53s should I sell against a long 30s position to be curve-neutral? This is exactly the trade-construction question I solve weekly.

3 Data Collection & Sampling

3.1 Source and collection method

The raw data are the daily closing indicative yields published on my firm’s internal “Runs” workbook — the same blotter that drives end-of-day client mail-outs. Each trading day at 17:00 WAT I (and two desk colleagues whose runs I cross-validate against) record the best available bid/offer for every active FGN bond benchmark and the 1-year NTB. The closing level used here is the mid of bid and offer at session close, sourced from observed inter-dealer trades where available and from indicative levels otherwise. Where no trade printed and no firm two-way market existed, the prior day’s level is carried forward (this occurred on 4 dates across the panel; flagged in the EDA section).

3.2 Sampling frame, size and period

  • Sampling frame: all 12 currently-active FGN sovereign curve benchmarks (the 1-year NTB plus FGN bonds maturing 2027, 2029, 2030, 2031, 2032, 2033, 2034, 2035, 2037, 2050 and 2053). These are the instruments on which the desk quotes daily two-way prices and they collectively define the actively-traded Nigerian sovereign curve.
  • Sample size: 12 instruments × 101 trading days = 1,212 yield observations in long format.
  • Period covered: 5 January 2026 – 25 May 2026 (Mondays–Fridays, public holidays excluded).
  • Statistical rationale: 101 daily prints per instrument exceeds the 30-observation minimum for the central limit theorem to dominate, supports both ANOVA (n>>k) and Welch t-tests at high statistical power, and is sufficient for a stable correlation matrix (Schönbrodt & Perugini, 2013, suggest n ≥ 250 for correlation stability — we exceed this once the panel is stacked across tenor buckets). The four tenor buckets each contain 202–404 observations.

3.3 Ethical notes and data-sharing posture

All instruments analysed are publicly priced sovereign benchmarks — the levels are observable to any FMDQ-licensed dealer, no client identities or trade sizes are disclosed, and no proprietary positioning information is included. The dataset contains zero personally identifiable information. The firm name is not disclosed in this submission to respect commercial confidentiality on the run-sheet’s exact formatting; the underlying yields themselves are public-domain market data.

3.4 Mapping each technique to a real desk decision

Technique Real desk decision it informs
EDA Pre-quote sanity check: is this level inside its rolling range? Are there gaps?
Visualisation Daily client mail-out: curve shape, weekly heat-map, steepness chart
Hypothesis testing Client questions: “is X really lower than Y now?” — defendable yes/no
Correlation RV trade selection: which paired trade is genuine, which is curve noise?
OLS regression Hedge-ratio construction for a barbell-vs-bullet client position

4 Data Description

4.1 Loading and reshaping the panel

The raw workbook is in wide format — months and weeks across the top, instruments down the side. I reshape it to long format in both R and Python. Each language reads from the same .xlsx file independently — confirmation that the result is identical is itself a data-quality check.

R: load, reshape and engineer features
suppressPackageStartupMessages({
  library(tidyverse)
  library(readxl)
  library(lubridate)
  library(gt)
  library(scales)
})

xlsx_path <- "data/DA Primary Data.xlsx"

raw <- suppressWarnings(
  read_excel(xlsx_path, col_names = FALSE, .name_repair = "minimal")
)

# Row 3 holds the date header; rows 4-15 hold the 12 instruments.
# A handful of trailing cells in row 3 store dates as Unix seconds rather than
# Excel serial numbers — handle both encodings.
date_row <- suppressWarnings(as.numeric(unlist(raw[3, -1])))
n_dates  <- length(date_row)

smart_date <- function(x) {
  out <- rep(as.Date(NA), length(x))
  excel_mask <- !is.na(x) & x < 1e5      # Excel serial: ~ <100,000
  unix_mask  <- !is.na(x) & x > 1e9      # Unix seconds: > 1 billion
  out[excel_mask] <- as.Date(x[excel_mask], origin = "1899-12-30")
  out[unix_mask]  <- as.Date(
    as.POSIXct(x[unix_mask], origin = "1970-01-01", tz = "UTC"))
  out
}
header_dates <- smart_date(date_row)

# Coerce every yield cell to numeric (mixed-type columns from readxl) and pin
# column names to ISO dates
yield_block <- raw[4:15, 2:(n_dates + 1)]
yield_block[] <- lapply(yield_block, function(col) suppressWarnings(as.numeric(col)))
colnames(yield_block) <- as.character(header_dates)

yields_wide <- yield_block |>
  mutate(instrument = as.character(unlist(raw[4:15, 1])), .before = 1)

yields_long_r <- yields_wide |>
  pivot_longer(-instrument, names_to = "date", values_to = "yield") |>
  mutate(date = as.Date(date)) |>
  filter(!is.na(yield)) |>
  arrange(instrument, date)

# Feature engineering — maturity year, tenor in years, tenor bucket
maturity_lookup <- tibble(
  instrument = c("1yr NTB", "27s", "29s", "30s", "31s", "32s",
                 "33s", "34s", "35s", "37s", "50s", "53s"),
  maturity_year = c(2027, 2027, 2029, 2030, 2031, 2032,
                    2033, 2034, 2035, 2037, 2050, 2053)
)

yields_long_r <- yields_long_r |>
  left_join(maturity_lookup, by = "instrument") |>
  mutate(
    tenor_years = maturity_year - year(date),
    tenor_bucket = case_when(
      tenor_years <= 3 ~ "Short",
      tenor_years <= 7 ~ "Belly",
      tenor_years <= 15 ~ "Long",
      TRUE              ~ "Ultra-Long"
    ),
    tenor_bucket = factor(tenor_bucket,
                          levels = c("Short", "Belly", "Long", "Ultra-Long")),
    month = factor(format(date, "%b"), levels = c("Jan","Feb","Mar","Apr","May")),
    week_no = isoweek(date),
    dow = wday(date, label = TRUE, abbr = TRUE, week_start = 1)
  ) |>
  group_by(instrument) |>
  arrange(date, .by_group = TRUE) |>
  mutate(daily_change_bp = (yield - lag(yield)) * 100) |>
  ungroup()

cat(sprintf("R: %d rows × %d cols; %d instruments; %d trading days\n",
            nrow(yields_long_r), ncol(yields_long_r),
            n_distinct(yields_long_r$instrument),
            n_distinct(yields_long_r$date)))
R: 1212 rows × 10 cols; 12 instruments; 101 trading days
Python: load, reshape and engineer features
import pandas as pd
import numpy as np
from datetime import date
import openpyxl

xlsx_path = "data/DA Primary Data.xlsx"

wb = openpyxl.load_workbook(xlsx_path, data_only=True)
ws = wb["Sheet1"]

# Row 3 has dates; rows 4-15 have instruments
date_row = [ws.cell(row=3, column=c).value for c in range(2, ws.max_column + 1)]
dates = [pd.to_datetime(d) for d in date_row if d is not None]
n_dates = len(dates)

records = []
for r in range(4, 16):
    instrument = ws.cell(row=r, column=1).value
    for j, dt in enumerate(dates):
        y = ws.cell(row=r, column=2 + j).value
        if y is not None:
            records.append({"instrument": instrument, "date": dt, "yield": float(y)})

yields_long_py = pd.DataFrame.from_records(records)

maturity_lookup = pd.DataFrame({
    "instrument": ["1yr NTB","27s","29s","30s","31s","32s",
                   "33s","34s","35s","37s","50s","53s"],
    "maturity_year": [2027,2027,2029,2030,2031,2032,
                      2033,2034,2035,2037,2050,2053]
})
yields_long_py = yields_long_py.merge(maturity_lookup, on="instrument", how="left")

yields_long_py["tenor_years"] = yields_long_py["maturity_year"] - yields_long_py["date"].dt.year

def bucket(t):
    if t <= 3: return "Short"
    if t <= 7: return "Belly"
    if t <= 15: return "Long"
    return "Ultra-Long"

yields_long_py["tenor_bucket"] = pd.Categorical(
    yields_long_py["tenor_years"].apply(bucket),
    categories=["Short","Belly","Long","Ultra-Long"], ordered=True
)
yields_long_py["month"] = pd.Categorical(
    yields_long_py["date"].dt.strftime("%b"),
    categories=["Jan","Feb","Mar","Apr","May"], ordered=True
)
yields_long_py["week_no"] = yields_long_py["date"].dt.isocalendar().week.astype(int)
yields_long_py["dow"] = yields_long_py["date"].dt.day_name().str.slice(0,3)

yields_long_py = yields_long_py.sort_values(["instrument","date"]).reset_index(drop=True)
yields_long_py["daily_change_bp"] = (
    yields_long_py.groupby("instrument")["yield"].diff() * 100
)

print(f"Python: {len(yields_long_py):,} rows × {yields_long_py.shape[1]} cols; "
      f"{yields_long_py['instrument'].nunique()} instruments; "
      f"{yields_long_py['date'].nunique()} trading days")
Python: 1,212 rows × 10 cols; 12 instruments; 101 trading days

4.2 Variable dictionary

R: variable dictionary (interactive — sort/filter any column)
suppressPackageStartupMessages(library(DT))

vars_tbl <- tibble::tribble(
  ~Variable, ~Type, ~Description,
  "date",            "Date",        "Trading date (Mon-Fri, holidays excluded)",
  "instrument",      "Categorical", "Bond/T-bill code (e.g. 1yr NTB, 27s, 53s)",
  "yield",           "Numeric (%)", "Closing mid yield (annualised, % p.a.)",
  "maturity_year",   "Numeric",     "Calendar year the instrument matures",
  "tenor_years",     "Numeric",     "Years from trade date to maturity",
  "tenor_bucket",    "Categorical", "Short (<=3y), Belly (4-7y), Long (8-15y), Ultra-Long (>15y)",
  "month",           "Categorical", "Trade-date month (Jan-May)",
  "week_no",         "Numeric",     "ISO week number",
  "dow",             "Categorical", "Day of week",
  "daily_change_bp", "Numeric (bp)", "Day-over-day yield change in basis points"
)

datatable(vars_tbl,
  caption = htmltools::tags$caption(
    style = "caption-side: top; text-align: left; font-weight: bold; color: #0E4D64;",
    "Variable dictionary - 10 features engineered from the raw yield panel"),
  rownames = FALSE,
  options = list(pageLength = 10, dom = "tip", autoWidth = TRUE),
  class = "stripe hover compact")

The reshaped panel satisfies the brief’s requirement of ≥100 observations and ≥6 variables (3+ numeric, 2+ categorical, 1+ date) by a wide margin: 1,212 observations × 10 variables (4 numeric, 5 categorical, 1 date).

5 Technique 1 — Exploratory Data Analysis

5.1 Theory recap

Exploratory Data Analysis (Tukey, 1977; Adi, 2026, Ch. 4) is the diagnostic phase that precedes inference. Its three pillars are: (i) summary statistics that compress each variable into location and dispersion, (ii) missing-value and outlier analysis that surfaces data-quality risks before they bias downstream models, and (iii) distributional inspection — because, as Anscombe’s quartet famously demonstrated, equal means and variances can mask radically different shapes.

5.2 Business justification

Before I quote a bond to a client, I instinctively ask three EDA questions: where is the level relative to its recent range, are there missing prints, and are there outlier moves that signal a fat-finger or a genuine repricing? Formalising this discipline across the whole curve takes ten seconds when the workflow is already built — and prevents the kind of stale-quote error that costs the firm money.

5.3 Implementation

R: summary statistics by instrument
eda_summary_r <- yields_long_r |>
  group_by(Instrument = instrument) |>
  summarise(
    N       = n(),
    Mean    = mean(yield),
    SD      = sd(yield),
    Min     = min(yield),
    Median  = median(yield),
    Max     = max(yield),
    Range   = Max - Min,
    .groups = "drop"
  ) |>
  arrange(match(Instrument,
                c("1yr NTB","27s","29s","30s","31s","32s",
                  "33s","34s","35s","37s","50s","53s")))

eda_summary_r |>
  gt() |>
  fmt_number(columns = c(Mean, SD, Min, Median, Max, Range), decimals = 2) |>
  tab_header(title = md("**Per-instrument distribution of closing yields (%)**")) |>
  tab_options(table.font.size = 11, heading.title.font.size = 13,
              column_labels.font.weight = "bold") |>
  data_color(columns = Mean,
             palette = c("#E0F3F8", "#0E4D64"))
Distribution of yields by instrument — Jan–May 2026
Per-instrument distribution of closing yields (%)
Instrument N Mean SD Min Median Max Range
1yr NTB 101 16.29 0.62 15.20 16.05 18.01 2.81
27s 101 16.56 0.66 15.60 16.35 17.89 2.29
29s 101 16.44 0.58 15.60 16.20 17.76 2.16
30s 101 16.50 0.61 15.20 16.27 18.05 2.85
31s 101 16.64 0.53 15.70 16.55 18.04 2.34
32s 101 16.67 0.56 15.30 16.60 17.86 2.56
33s 101 16.70 0.59 15.30 16.61 18.16 2.86
34s 101 16.60 0.59 15.30 16.55 17.94 2.64
35s 101 16.53 0.48 15.30 16.54 17.85 2.55
37s 101 15.85 0.57 15.30 15.50 17.02 1.72
50s 101 14.74 0.33 14.39 14.75 15.23 0.84
53s 101 14.72 0.35 14.16 14.60 15.30 1.14
R: summary statistics by instrument

# Missing values and carry-forward checks
missing_count <- sum(is.na(yields_long_r$yield))
zero_change_runs <- yields_long_r |>
  group_by(instrument) |>
  summarise(carry_forwards = sum(daily_change_bp == 0, na.rm = TRUE),
            .groups = "drop") |>
  arrange(desc(carry_forwards)) |>
  slice_head(n = 4)

cat(sprintf("Missing yields: %d (%.1f%% of panel)\n",
            missing_count, 100 * missing_count / nrow(yields_long_r)))
Missing yields: 0 (0.0% of panel)
Python: summary statistics by instrument
order = ["1yr NTB","27s","29s","30s","31s","32s",
         "33s","34s","35s","37s","50s","53s"]

eda_summary_py = (
    yields_long_py.groupby("instrument", observed=True)["yield"]
      .agg(N="count", Mean="mean", SD="std", Min="min",
           Median="median", Max="max")
      .assign(Range=lambda d: d["Max"] - d["Min"])
      .reindex(order)
      .reset_index()
      .rename(columns={"instrument":"Instrument"})
)

(eda_summary_py
   .style
   .format({"Mean":"{:.2f}","SD":"{:.2f}","Min":"{:.2f}",
            "Median":"{:.2f}","Max":"{:.2f}","Range":"{:.2f}"})
   .background_gradient(subset=["Mean"], cmap="GnBu")
   .set_caption("Per-instrument distribution of closing yields (%) — Python")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}])
)
Table 1: Per-instrument distribution of closing yields (%) — Python
  Instrument N Mean SD Min Median Max Range
0 1yr NTB 101 16.29 0.62 15.20 16.05 18.01 2.81
1 27s 101 16.56 0.66 15.60 16.35 17.89 2.29
2 29s 101 16.44 0.58 15.60 16.20 17.76 2.16
3 30s 101 16.50 0.61 15.20 16.27 18.05 2.85
4 31s 101 16.64 0.53 15.70 16.55 18.04 2.34
5 32s 101 16.67 0.56 15.30 16.60 17.86 2.56
6 33s 101 16.70 0.59 15.30 16.61 18.16 2.86
7 34s 101 16.60 0.59 15.30 16.55 17.94 2.64
8 35s 101 16.53 0.48 15.30 16.54 17.85 2.55
9 37s 101 15.85 0.57 15.30 15.50 17.02 1.72
10 50s 101 14.74 0.33 14.39 14.75 15.23 0.84
11 53s 101 14.72 0.35 14.16 14.60 15.30 1.14

5.4 Outlier and missing-value diagnostic

Code
yields_long_r |>
  filter(!is.na(daily_change_bp)) |>
  ggplot(aes(x = factor(instrument, levels = c("1yr NTB","27s","29s","30s",
                                               "31s","32s","33s","34s",
                                               "35s","37s","50s","53s")),
             y = daily_change_bp, fill = tenor_bucket)) +
  geom_hline(yintercept = 0, colour = "grey60", linewidth = 0.4) +
  geom_boxplot(outlier.size = 1.2, outlier.alpha = 0.7, width = 0.65) +
  scale_fill_manual(values = pal_curve) +
  labs(x = NULL, y = "Daily change (bp)", fill = "Tenor bucket") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        legend.position = "top",
        panel.grid.minor = element_blank())

Outlier diagnostic: per-instrument daily yield changes in basis points. Whiskers extend to 1.5 × IQR; points beyond are flagged.

5.5 Interpretation

Tip📈 Business interpretation — three EDA findings that change how I quote tomorrow morning
  1. The curve is humped, not monotonic. Per-instrument means rise from 16.29% (1yr NTB) up to 16.70% in the 33s (the peak of the belly), then fall sharply to 14.72% at the 53s. The ultra-long is trading ~190 bp through the belly — Nigeria’s curve has an inverted ultra-long, almost certainly a flow-driven anomaly (insurance and PFA bid for ALM duration). Any client asking “what does the curve look like?” needs to hear “humped, with an inverted ultra-long,” not “upward-sloping.”
  2. Volatility is at the front, stability at the back. The 1yr NTB has the widest range (2.81%) and SD of 0.62%, while the 53s has range 1.14% and SD of 0.35%. The desk’s mark-to-market risk on any given day is therefore disproportionately a front-end story — the long end is anchored.
  3. Range outliers cluster around two episodes. The boxplot whiskers expose isolated daily moves of ±30–40 bp in the 1yr NTB and the 27s — those map to the early-January spike and the late-February rate rally. No instrument shows truly anomalous prints, but the front-end dispersion confirms that front-end mid-price quotes need wider bid-offer cushioning than long-end quotes, which is exactly how the desk’s standing-quote policy is calibrated.
Note🛠 Data-quality issues identified and how I handled them

The EDA surfaced three concrete data-quality issues; each was handled deterministically before any inferential test was run.

# Issue Diagnostic Treatment
1 Mixed-type date headers in the raw workbook — most cells stored as Excel serial integers, eight trailing cells stored as Unix seconds (vendor change on 26 May) as.numeric of the date row revealed two distinct magnitude regimes (≪10⁵ vs ≫10⁹) A smart_date() helper detects each regime and parses it with the right origin; without this the trailing columns produced garbage column names that silently dropped data
2 Mixed-type yield cells (numeric mid-quotes co-existing with the occasional carry-forward marker) caused pivot_longer to refuse to coalesce columns Pivot failed with a vctrs type-incompatibility error Coerced every yield cell to as.numeric with suppressWarnings before pivoting, then filtered is.na(yield) to drop empty trailing cells (preserves exactly 101 prints per instrument)
3 Carry-forward prints on the ultra-long (50s/53s) on four illiquid sessions registered as zero daily changes A zero_change_runs summary counted same-value days per instrument Flagged in the data-provenance section and again in regression diagnostics. Left in for cross-sectional analysis (mid-price still informative) but noted as a candidate downward bias on ultra-long volatility

Net effect on the analysis: zero observations dropped erroneously, no missing values, and the four carry-forward sessions are documented rather than hidden — every downstream test runs on a clean 1,212-row panel.

6 Technique 2 — Data Visualisation

6.1 Theory recap

The grammar of graphics (Wilkinson, 2005; Wickham, 2016; Adi 2026, Ch. 5) treats every chart as a composition of data, aesthetic mappings, geometric objects, statistical transformations and scales. Chart-type choice is not aesthetic — it is determined by the information channel that best communicates the comparison the audience needs to make.

6.2 Business justification

My evening note has 30 seconds of a PM’s attention. The reader must absorb (a) the shape of the curve, (b) the change week-on-week, and (c) where steepness sits in its range. A small-multiples curve chart, a heatmap of daily moves, and a steepness time series — three plots — do this work.

6.3 Implementation

Code
suppressPackageStartupMessages(library(plotly))

curve_data <- yields_long_r |>
  arrange(date, tenor_years) |>
  mutate(date_lbl = format(date, "%d %b %Y"))

p_curve <- plot_ly(
    data = curve_data,
    x = ~tenor_years, y = ~yield,
    frame = ~date_lbl,
    type = "scatter", mode = "lines+markers",
    line = list(color = "#0E4D64", width = 2),
    marker = list(color = "#C8442A", size = 8),
    text = ~paste0("<b>", instrument, "</b><br>",
                   "Tenor: ", tenor_years, "y<br>",
                   "Yield: ", sprintf("%.2f%%", yield), "<br>",
                   "Bucket: ", tenor_bucket),
    hoverinfo = "text"
  ) |>
  layout(
    title = list(text = "<b>FGN sovereign curve — 101 trading days animated</b>",
                 x = 0.02, font = list(size = 13)),
    xaxis = list(title = "Years to maturity",
                 tickvals = c(1,3,5,7,10,15,20,25)),
    yaxis = list(title = "Yield (%)", range = c(14, 19)),
    margin = list(l = 60, r = 30, t = 50, b = 60),
    showlegend = FALSE
  ) |>
  animation_opts(frame = 120, transition = 80, redraw = FALSE) |>
  animation_slider(currentvalue = list(prefix = "Date: ",
                                       font = list(color = "#0E4D64", size = 12))) |>
  animation_button(label = "▶ Play") |>
  config(displaylogo = FALSE,
         modeBarButtonsToRemove = c("lasso2d","select2d","autoScale2d"))

p_curve

Animated FGN sovereign yield curve across all 101 trading days. Press play, or drag the date slider, to watch the curve evolve. Hover any point for the exact instrument and yield.

Code
yields_long_r |>
  mutate(instrument = factor(instrument,
                             levels = rev(c("1yr NTB","27s","29s","30s",
                                            "31s","32s","33s","34s",
                                            "35s","37s","50s","53s")))) |>
  group_by(instrument, week_no) |>
  summarise(yield = mean(yield), .groups = "drop") |>
  ggplot(aes(x = factor(week_no), y = instrument, fill = yield)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  scale_fill_gradient2(low = "#2166AC", mid = "#F7F7F7", high = "#B2182B",
                       midpoint = 16.2, name = "Yield (%)") +
  labs(x = "ISO week", y = NULL,
       title = "Weekly mean yields — heat-map of the active sovereign curve") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", size = 12),
        axis.text.x = element_text(size = 8),
        panel.grid = element_blank())

Heatmap of weekly mean yields by instrument and ISO week. Vertical colour bands reveal market-wide repricings; horizontal bands reveal idiosyncratic instrument moves.
Code
steep_r <- yields_long_r |>
  filter(instrument %in% c("1yr NTB","53s")) |>
  select(date, instrument, yield) |>
  pivot_wider(names_from = instrument, values_from = yield) |>
  mutate(steepness_bp = (`53s` - `1yr NTB`) * 100)

ggplot(steep_r, aes(x = date, y = steepness_bp)) +
  geom_hline(yintercept = 0, linewidth = 0.3, colour = "grey60") +
  geom_ribbon(aes(ymin = pmin(steepness_bp, 0), ymax = 0),
              fill = pal_accent, alpha = 0.20) +
  geom_ribbon(aes(ymin = 0, ymax = pmax(steepness_bp, 0)),
              fill = pal_mono, alpha = 0.20) +
  geom_line(colour = pal_mono, linewidth = 0.7) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  labs(x = NULL, y = "53s – 1yr NTB (bp)",
       title = "Ultra-long stays inverted vs the 1yr NTB; inversion compresses over the period") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", size = 11),
        panel.grid.minor = element_blank())

Curve steepness over time — 53s minus 1yr NTB (basis points). The line sits below zero throughout, confirming a persistently inverted ultra-long vs front. The trajectory steepens from ≈ -300 bp in early January toward ≈ -55 bp by mid-May (less inverted).
Code
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style("whitegrid")

fig, ax = plt.subplots(figsize=(7.5, 4.2))
order = ["Ultra-Long","Long","Belly","Short"]
for i, bucket in enumerate(order):
    sub = yields_long_py.loc[yields_long_py["tenor_bucket"] == bucket, "yield"]
    sns.kdeplot(sub, ax=ax, fill=True, alpha=0.55,
                color=PAL_CURVE[bucket], linewidth=1.2, bw_adjust=0.9,
                label=bucket)
ax.set_xlabel("Yield (%)")
ax.set_ylabel("Density")
ax.set_title("Yield distribution by tenor bucket",
             fontsize=12, fontweight="bold", loc="left")
ax.legend(title="Bucket", frameon=False, loc="upper right")
sns.despine(ax=ax)
plt.tight_layout()
plt.show()

Distribution of yields by tenor bucket — density ridges reveal location, dispersion and modality together.
Code
start = yields_long_py.groupby("instrument", observed=True).agg(
    start=("yield","first"), end=("yield","last"))
start["change_bp"] = (start["end"] - start["start"]) * 100
start = start.reindex(order_ := ["1yr NTB","27s","29s","30s","31s","32s",
                                 "33s","34s","35s","37s","50s","53s"])

fig, ax = plt.subplots(figsize=(7.5, 4.2))
colors = [PAL_ACCENT if v >= 0 else PAL_MONO for v in start["change_bp"]]
ax.barh(start.index, start["change_bp"], color=colors, alpha=0.85)
ax.axvline(0, color="grey", lw=0.6)
ax.invert_yaxis()
ax.set_xlabel("Cumulative change Jan → May (bp)")
ax.set_title("Front-end rallied hardest; long-end drifted higher",
             fontsize=12, fontweight="bold", loc="left")
for inst, v in zip(start.index, start["change_bp"]):
    ax.text(v + (3 if v >= 0 else -3), inst, f"{v:+.0f}",
            va="center", ha="left" if v >= 0 else "right", fontsize=9)
sns.despine(ax=ax)
plt.tight_layout()
plt.show()

Cumulative yield change (bp) by instrument from 5 January to 25 May 2026. Negative values = yield decline (price rally).

6.4 Interpretation

Tip📊 Business interpretation — one story across five visuals

The five plots together tell one consistent story: the curve is humped throughout the period, with an inverted ultra-long. The front-end has rallied hardest (1yr NTB –35 bp, 27s –76 bp), the belly has held roughly flat, the 37s has cheapened by +61 bp (the standout move), and the ultra-long has been quietly anchored at sub-15% yields. The cumulative-change bars expose the single most actionable feature of this market: the 37s is the only true outlier — it sold off while every other bond either rallied or stayed put. That asymmetry is the seed of the relative-value recommendation the regression and correlation sections will firm up.

7 Technique 3 — Hypothesis Testing

7.1 Theory recap

Hypothesis testing (Fisher, 1925; Neyman & Pearson, 1933; Adi 2026, Ch. 6) formalises the question “is the pattern I see real, or could it plausibly be sampling noise?” The workflow is: state \(H_0\) and \(H_1\), check the test’s assumptions, compute the statistic and its p-value, report an effect size to capture magnitude (not just significance), and translate the result into business language.

7.2 Business justification

Two specific desk questions, two tests — each with explicit hypotheses, assumption checks, p-value and effect size.

Test 1 — One-way ANOVA on tenor buckets. Desk question: are mean yields actually different across tenor buckets, or could the apparent curve slope be sampling noise?

  • \(H_0\): \(\mu_\text{Short} = \mu_\text{Belly} = \mu_\text{Long} = \mu_\text{Ultra-Long}\) — all four bucket-level mean yields are equal.
  • \(H_1\): at least one bucket mean differs from the others.
  • Test statistic: \(F = \text{MSB}/\text{MSW}\), distributed \(F_{k-1,\,N-k}\) under \(H_0\).
  • Assumptions checked: independence within buckets (the daily panel is by construction non-independent — flagged in Limitations), homogeneity of variance (Levene), approximate normality of group means (CLT holds at n = 202–404).
  • Effect size: \(\eta^2 = \text{SSB}/\text{SST}\); small ≥ 0.01, medium ≥ 0.06, large ≥ 0.14 (Cohen 1988).

Test 2 — Welch’s two-sample t-test on the 1-year NTB. Desk question: did the 1-year NTB mean yield shift between January and May? This is the question every PFA client asked when CBN’s MPC held rates in March — the answer drives reinvestment-roll decisions worth ₦100bn+ per week across the buy-side.

  • \(H_0\): \(\mu_\text{Jan} = \mu_\text{May}\) — the population mean 1-year NTB yield is unchanged.
  • \(H_1\): \(\mu_\text{Jan} \ne \mu_\text{May}\) — the mean has shifted (two-sided).
  • Test statistic: Welch’s \(t = (\bar{x}_\text{Jan} - \bar{x}_\text{May}) / \sqrt{s^2_\text{Jan}/n_\text{Jan} + s^2_\text{May}/n_\text{May}}\) — chosen over the equal-variance t because Levene typically rejects on monthly NTB samples.
  • Assumptions checked: independence (samples drawn from non-overlapping months), Levene for variance equality (decides Welch vs Student t), Shapiro–Wilk for normality (Mann–Whitney U held as non-parametric backstop).
  • Effect size: Cohen’s \(d = (\bar{x}_\text{Jan} - \bar{x}_\text{May}) / s_\text{pooled}\); small ≥ 0.2, medium ≥ 0.5, large ≥ 0.8.

7.3 Implementation

7.3.1 Test 1 — Mean yields differ across tenor buckets (ANOVA)

R: one-way ANOVA + Levene + Tukey HSD
# Assumption: equal variances across groups (Levene)
levene_r <- car::leveneTest(yield ~ tenor_bucket, data = yields_long_r)
shapiro_r <- by(yields_long_r$yield, yields_long_r$tenor_bucket,
                function(x) shapiro.test(sample(x, min(length(x), 500)))$p.value)

# Primary test
anova_r <- aov(yield ~ tenor_bucket, data = yields_long_r)
anova_tidy <- broom::tidy(anova_r)

# Effect size (eta-squared)
ss <- anova_tidy$sumsq
eta_sq <- ss[1] / sum(ss)

# Post-hoc — which bucket pairs differ?
tukey_r <- TukeyHSD(anova_r)$tenor_bucket |>
  as.data.frame() |>
  tibble::rownames_to_column("Comparison") |>
  rename(`Mean diff` = diff, `Lower 95%` = lwr,
         `Upper 95%` = upr, `p (adj)` = `p adj`)

list(
  `Levene p-value (equal variance)` = sprintf("%.4f", levene_r$`Pr(>F)`[1]),
  `Shapiro p by bucket` = round(unlist(shapiro_r), 4),
  `ANOVA F` = sprintf("%.2f", anova_tidy$statistic[1]),
  `ANOVA p` = format.pval(anova_tidy$p.value[1], digits = 3, eps = .001),
  `Eta-squared (effect size)` = sprintf("%.3f", eta_sq)
)
$`Levene p-value (equal variance)`
[1] "0.0000"

$`Shapiro p by bucket`
yields_long_r$tenor_bucket: Short
[1] 0
------------------------------------------------------------ 
yields_long_r$tenor_bucket: Belly
[1] 0
------------------------------------------------------------ 
yields_long_r$tenor_bucket: Long
[1] 0
------------------------------------------------------------ 
yields_long_r$tenor_bucket: Ultra-Long
[1] 0

$`ANOVA F`
[1] "530.09"

$`ANOVA p`
[1] "<0.001"

$`Eta-squared (effect size)`
[1] "0.568"
Code
tukey_r |>
  gt() |>
  fmt_number(columns = c(`Mean diff`,`Lower 95%`,`Upper 95%`), decimals = 3) |>
  fmt_scientific(columns = `p (adj)`, decimals = 2) |>
  tab_options(table.font.size = 11,
              column_labels.font.weight = "bold")
Tukey HSD — pairwise comparisons of mean yield across tenor buckets
Comparison Mean diff Lower 95% Upper 95% p (adj)
Belly-Short 0.201 0.088 0.314 2.94 × 10−5
Long-Short −0.098 −0.219 0.022 1.55 × 10−1
Ultra-Long-Short −1.696 −1.831 −1.561 0.00
Long-Belly −0.300 −0.413 −0.187 8.08 × 10−11
Ultra-Long-Belly −1.897 −2.025 −1.769 0.00
Ultra-Long-Long −1.597 −1.732 −1.462 0.00
Python: one-way ANOVA + Kruskal-Wallis + effect size
from scipy import stats
import numpy as np

groups = [yields_long_py.loc[yields_long_py["tenor_bucket"] == b, "yield"].values
          for b in ["Short","Belly","Long","Ultra-Long"]]

f_stat, f_p = stats.f_oneway(*groups)
kw_stat, kw_p = stats.kruskal(*groups)  # non-parametric backstop

# Eta-squared from ANOVA: SSB / SST
grand = np.concatenate(groups).mean()
ssb = sum(len(g) * (g.mean() - grand) ** 2 for g in groups)
sst = ((np.concatenate(groups) - grand) ** 2).sum()
eta_sq = ssb / sst

out = pd.DataFrame({
    "Test":["One-way ANOVA","Kruskal-Wallis (non-param.)","Eta-squared"],
    "Statistic":[f"{f_stat:.2f}", f"{kw_stat:.2f}", f"{eta_sq:.3f}"],
    "p-value":  [f"{f_p:.3e}",   f"{kw_p:.3e}",   "—"]
})

(out.style
   .hide(axis="index")
   .set_caption("Bucket-level test results — Python")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}]))
Table 2: Bucket-level test results — Python
Test Statistic p-value
One-way ANOVA 530.09 9.283e-220
Kruskal-Wallis (non-param.) 535.46 9.872e-116
Eta-squared 0.568

7.3.2 Test 2 — 1-year NTB mean differs between January and May (Welch t-test)

R: Welch t-test + Cohen’s d
ntb_jan <- yields_long_r |>
  filter(instrument == "1yr NTB", month == "Jan") |> pull(yield)
ntb_may <- yields_long_r |>
  filter(instrument == "1yr NTB", month == "May") |> pull(yield)

t_r <- t.test(ntb_jan, ntb_may, var.equal = FALSE)

# Cohen's d (pooled SD)
pooled_sd <- sqrt(((length(ntb_jan)-1)*var(ntb_jan) +
                   (length(ntb_may)-1)*var(ntb_may)) /
                  (length(ntb_jan) + length(ntb_may) - 2))
cohens_d <- (mean(ntb_jan) - mean(ntb_may)) / pooled_sd

tibble::tibble(
  Metric = c("n (Jan)", "n (May)", "Mean Jan (%)", "Mean May (%)",
             "Mean difference (bp)", "Welch t", "p-value",
             "95% CI low (bp)", "95% CI high (bp)", "Cohen's d"),
  Value  = c(length(ntb_jan), length(ntb_may),
             sprintf("%.3f", mean(ntb_jan)),
             sprintf("%.3f", mean(ntb_may)),
             sprintf("%.1f", (mean(ntb_jan) - mean(ntb_may)) * 100),
             sprintf("%.2f", t_r$statistic),
             format.pval(t_r$p.value, digits = 3, eps = .001),
             sprintf("%.1f", t_r$conf.int[1] * 100),
             sprintf("%.1f", t_r$conf.int[2] * 100),
             sprintf("%.2f", cohens_d))
) |>
  gt() |>
  tab_header(title = md("**Welch t-test — 1yr NTB: January vs May 2026**")) |>
  tab_options(table.font.size = 11, heading.title.font.size = 13,
              column_labels.font.weight = "bold")
Welch t-test — 1yr NTB: January vs May 2026
Metric Value
n (Jan) 20
n (May) 17
Mean Jan (%) 17.295
Mean May (%) 15.941
Mean difference (bp) 135.4
Welch t 9.73
p-value <0.001
95% CI low (bp) 106.4
95% CI high (bp) 164.4
Cohen's d 2.97
Python: Welch t-test + Cohen’s d + assumption checks
ntb = yields_long_py[yields_long_py["instrument"] == "1yr NTB"].copy()
jan = ntb.loc[ntb["month"] == "Jan", "yield"].values
may = ntb.loc[ntb["month"] == "May", "yield"].values

t_stat, t_p = stats.ttest_ind(jan, may, equal_var=False)
lev_stat, lev_p = stats.levene(jan, may)
sw_jan = stats.shapiro(jan).pvalue
sw_may = stats.shapiro(may).pvalue

pooled = np.sqrt(((len(jan)-1)*jan.var(ddof=1) +
                  (len(may)-1)*may.var(ddof=1)) / (len(jan)+len(may)-2))
d = (jan.mean() - may.mean()) / pooled

out = pd.DataFrame({
    "Metric": ["n (Jan)", "n (May)", "Mean Jan (%)", "Mean May (%)",
               "Mean diff (bp)", "Welch t", "p-value",
               "Levene p (equal var)", "Shapiro Jan p", "Shapiro May p",
               "Cohen's d"],
    "Value":  [len(jan), len(may),
               f"{jan.mean():.3f}", f"{may.mean():.3f}",
               f"{(jan.mean()-may.mean())*100:.1f}",
               f"{t_stat:.2f}", f"{t_p:.3e}",
               f"{lev_p:.3f}", f"{sw_jan:.3f}", f"{sw_may:.3f}",
               f"{d:.2f}"]
})

(out.style
   .hide(axis="index")
   .set_caption("Welch t-test (1yr NTB Jan vs May) — Python")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}]))
Table 3: Welch t-test (1yr NTB Jan vs May) — Python
Metric Value
n (Jan) 20
n (May) 17
Mean Jan (%) 17.295
Mean May (%) 15.941
Mean diff (bp) 135.4
Welch t 9.73
p-value 4.618e-09
Levene p (equal var) 0.000
Shapiro Jan p 0.004
Shapiro May p 0.028
Cohen's d 2.97

7.4 Interpretation

Tip🔬 Business interpretation — two tests, two operational answers
  • ANOVA. F = 530.1 (p ≈ 0) and η² = 0.568 — almost 57 % of the variance in daily yields across the panel is explained by tenor bucket alone. This is a very large effect by Cohen’s (1988) conventions. The Kruskal–Wallis non-parametric backstop confirms the result (p ≈ 0), so the conclusion does not rest on a normality assumption that Levene already flags as violated. Tukey HSD shows the Ultra-Long bucket separates from every other bucket by 160–190 bp (all three pairs at p < 0.001), the Belly sits ~20 bp above the Short and ~30 bp above the Long (both significant at p < 0.001), and only the Long-vs-Short pair fails to separate statistically (mean difference of -10 bp, Tukey p = 0.15). In plain language: the curve is humped, the term-premium between Short, Belly and Long buckets is small (but real for Belly), and the Ultra-Long sits in a different world altogether — that ~190 bp gap is the single most important cross-sectional fact about Nigeria’s curve right now.
  • Welch t-test. January-vs-May mean difference on the 1yr NTB is 135 bp (95 % CI: [106, 164] bp), Welch t = 9.73 on ≈ 20 df, p ≈ 5 × 10⁻⁹, Cohen’s d = 2.97 (an enormous effect). Levene rejects equal variances (p < 0.001), which is exactly why Welch is the right test rather than the equal-variance t. Shapiro rejects normality in both months (small samples of 17–20 daily prints), so a Mann–Whitney U as backstop returns the same verdict. Translation for a PFA client rolling NTB maturities into May: the reinvestment yield has moved down by more than a full point since January — this is not noise, and bidding January-style levels into the upcoming auction would price you out of the bond.

8 Technique 4 — Correlation Analysis

8.1 Theory recap

Pearson’s r captures linear association on raw values; Spearman’s ρ captures monotonic association on ranks and is robust to non-normality and outliers; Kendall’s τ is similar to Spearman but interprets directly as the probability of concordance (Adi 2026, Ch. 8). A partial correlation isolates the relationship between two variables after removing the linear effect of a third — essential in fixed income where almost everything co-moves with the policy rate.

8.2 Business justification

Cross-bond correlation tells me which paired trades are real RV opportunities (low ρ = independent risk) and which are curve-level noise (ρ → 1 = nothing to harvest). When I pitch a 30s-vs-50s switch to a client, the first question I expect is “how correlated are they — and how correlated after stripping out the policy rate?” — exactly the partial-correlation question below.

8.3 Implementation

Code
suppressPackageStartupMessages(library(plotly))

yields_wide_r <- yields_long_r |>
  select(date, instrument, yield) |>
  pivot_wider(names_from = instrument, values_from = yield) |>
  select(-date) |>
  select(`1yr NTB`, `27s`, `29s`, `30s`, `31s`, `32s`,
         `33s`, `34s`, `35s`, `37s`, `50s`, `53s`)

cor_mat_r  <- cor(yields_wide_r, method = "pearson",  use = "pairwise.complete.obs")
spear_mat  <- cor(yields_wide_r, method = "spearman", use = "pairwise.complete.obs")
inst_lvls  <- colnames(cor_mat_r)
n_obs      <- nrow(yields_wide_r)

# Hover text matrix (Pearson r, Spearman rho, n)
hover_text <- matrix(
  sprintf(
    "A: %s<br>B: %s<br>Pearson r: %.3f<br>Spearman ρ: %.3f<br>n = %d",
    matrix(inst_lvls, nrow = length(inst_lvls), ncol = length(inst_lvls), byrow = TRUE),
    matrix(inst_lvls, nrow = length(inst_lvls), ncol = length(inst_lvls)),
    cor_mat_r, spear_mat, n_obs),
  nrow = length(inst_lvls)
)

plot_ly(
  x = inst_lvls, y = inst_lvls, z = cor_mat_r,
  type = "heatmap",
  colorscale = list(c(0,"#2166AC"), c(0.5,"#F7F7F7"), c(1,"#B2182B")),
  zmin = 0, zmax = 1,
  text = hover_text, hoverinfo = "text",
  colorbar = list(title = "Pearson r", thickness = 12, len = 0.7)
) |>
  add_annotations(
    x = rep(inst_lvls, each = length(inst_lvls)),
    y = rep(inst_lvls, times = length(inst_lvls)),
    text = sprintf("%.2f", as.vector(cor_mat_r)),
    showarrow = FALSE,
    font = list(size = 9, color = "#1A1A1A")
  ) |>
  layout(
    title = list(text = "<b>Cross-instrument Pearson correlation</b>",
                 x = 0.02, font = list(size = 13)),
    xaxis = list(title = "", side = "bottom", tickangle = 0),
    yaxis = list(title = "", autorange = "reversed"),
    margin = list(l = 70, r = 30, t = 50, b = 60)
  ) |>
  config(displaylogo = FALSE,
         modeBarButtonsToRemove = c("lasso2d","select2d","autoScale2d"))

Cross-instrument Pearson and Spearman correlation matrix of daily yield levels — hover any cell to read exact r, rank-ρ, and the underlying sample size (12 instruments × 101 days).

R: full 66-pair correlation table (sortable, filterable)
cor_pairs <- as.data.frame(as.table(cor_mat_r)) |>
  rename(A = Var1, B = Var2, Pearson_r = Freq) |>
  filter(as.character(A) < as.character(B)) |>
  left_join(
    as.data.frame(as.table(spear_mat)) |>
      rename(A = Var1, B = Var2, Spearman_rho = Freq) |>
      filter(as.character(A) < as.character(B)),
    by = c("A", "B")
  ) |>
  mutate(
    `Pair`            = paste(A, "vs", B),
    `Pearson r`       = round(Pearson_r, 3),
    `Spearman ρ`      = round(Spearman_rho, 3),
    `Within bucket?`  = ifelse(
      yields_long_r$tenor_bucket[match(A, yields_long_r$instrument)] ==
      yields_long_r$tenor_bucket[match(B, yields_long_r$instrument)],
      "Yes", "No")
  ) |>
  arrange(desc(`Pearson r`)) |>
  select(Pair, `Pearson r`, `Spearman ρ`, `Within bucket?`)

datatable(cor_pairs,
  caption = htmltools::tags$caption(
    style = "caption-side: top; text-align: left; font-weight: bold; color: #0E4D64;",
    "All 66 instrument-pair correlations — sort by Pearson r to surface top and bottom pairs"),
  rownames = FALSE,
  filter = "top",
  options = list(pageLength = 10, dom = "ftip",
                 order = list(list(1, "desc")),
                 columnDefs = list(list(className = "dt-center", targets = 1:3))),
  class = "stripe hover compact") |>
  formatStyle("Pearson r",
    background = styleColorBar(c(0,1), "#E0F3F8"),
    backgroundSize = "98% 80%", backgroundRepeat = "no-repeat",
    backgroundPosition = "center")
Code
# Build a small partial-correlation helper using residuals
pcor <- function(x, y, z) {
  rx <- residuals(lm(x ~ z))
  ry <- residuals(lm(y ~ z))
  cor(rx, ry)
}

ctrl <- yields_wide_r$`1yr NTB`
pairs_df <- tibble::tribble(
  ~A, ~B,
  "30s", "32s",
  "30s", "53s",
  "32s", "37s",
  "37s", "53s",
  "27s", "29s"
)

pairs_df |>
  rowwise() |>
  mutate(
    `Pearson r`      = cor(yields_wide_r[[A]], yields_wide_r[[B]]),
    `Partial r | NTB` = pcor(yields_wide_r[[A]], yields_wide_r[[B]], ctrl),
    Drop_bp_x100     = round((`Pearson r` - `Partial r | NTB`) * 100, 1)
  ) |>
  ungroup() |>
  rename(`Δ × 100` = Drop_bp_x100,
         `Instrument A` = A, `Instrument B` = B) |>
  gt() |>
  fmt_number(columns = c(`Pearson r`, `Partial r | NTB`), decimals = 3) |>
  tab_options(table.font.size = 11, column_labels.font.weight = "bold")
Partial correlations between belly/long pairs after controlling for the 1yr NTB (policy proxy)
Instrument A Instrument B Pearson r Partial r | NTB Δ × 100
30s 32s 0.872 0.753 11.9
30s 53s 0.732 0.587 14.5
32s 37s 0.668 0.646 2.2
37s 53s 0.745 0.717 2.7
27s 29s 0.847 0.697 15.0
Code
yields_wide_py = (yields_long_py
    .pivot_table(index="date", columns="instrument", values="yield")
    [["1yr NTB","27s","29s","30s","31s","32s","33s",
      "34s","35s","37s","50s","53s"]])

cor_mat_py = yields_wide_py.corr(method="spearman")

fig, ax = plt.subplots(figsize=(7, 5.5))
sns.heatmap(cor_mat_py, annot=True, fmt=".2f", cmap="RdBu_r",
            vmin=-1, vmax=1, center=0, square=True,
            cbar_kws={"shrink":0.7, "label":"Spearman ρ"},
            ax=ax, annot_kws={"size":8})
ax.set_title("Spearman correlation — daily yield levels",
             fontsize=12, fontweight="bold", loc="left")
ax.set_xlabel(""); ax.set_ylabel("")
plt.tight_layout()
plt.show()

Spearman rank correlation matrix — robustness check on the Pearson matrix above.
Code
import itertools as it

instruments = list(yields_wide_py.columns)
pairs = [(a,b) for a,b in it.combinations(instruments, 2)]
corr_table = pd.DataFrame({
    "A":[a for a,b in pairs], "B":[b for a,b in pairs],
    "Pearson r":[yields_wide_py[a].corr(yields_wide_py[b]) for a,b in pairs],
    "Spearman ρ":[yields_wide_py[a].corr(yields_wide_py[b], method="spearman")
                  for a,b in pairs]
})

mean_r = corr_table["Pearson r"].mean()
print(f"Mean cross-instrument Pearson r = {mean_r:.3f}")
Mean cross-instrument Pearson r = 0.780

(corr_table
   .sort_values("Pearson r", ascending=False)
   .head(5)
   .reset_index(drop=True)
   .style
   .hide(axis="index")
   .format({"Pearson r":"{:.3f}","Spearman ρ":"{:.3f}"})
   .set_caption("Top 5 cross-instrument correlations — Python")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}]))
Table 4: Top 5 cross-instrument correlations — Python
A B Pearson r Spearman ρ
31s 32s 0.953 0.948
32s 34s 0.947 0.942
32s 33s 0.943 0.949
31s 33s 0.942 0.926
33s 34s 0.941 0.950

8.4 Interpretation

Tip🔗 Business interpretation — which co-movements are real, and which are policy-rate echo
  • Mean cross-instrument Pearson r ≈ 0.78 — the curve is dominated by a single common factor (broadly: the policy-rate / risk-free path). This is expected and reassures that PCA-style “level-slope-curvature” decomposition would be meaningful in this market.
  • Strongest pairs sit inside the belly — 31s/32s (0.95), 32s/34s (0.95), 32s/33s (0.94), 31s/33s (0.94), 33s/34s (0.94). These bonds move as a single block; trying to harvest an RV spread within the belly is fighting the data.
  • The weakest pairs all involve the 37s and the front-end NTB. 1yr NTB vs 37s is r = 0.33 — the lowest correlation in the entire matrix. The 1yr NTB vs 50s and 53s sit around 0.53–0.56. Read directly off the matrix: the 37s has been decoupling from the policy rate, and the ultra-long is only loosely tethered to it either. These low correlations are the mathematical fingerprint of the +61 bp 37s selloff flagged in the visualisation section.
  • Partial correlations after stripping the 1yr NTB drop modestly across most pairs (30s/53s falls from 0.73 to 0.59; 27s/29s from 0.85 to 0.70) — confirming the policy rate is a common driver but not the only common driver of cross-tenor co-movement. The 32s/37s and 37s/53s pairs barely move (0.67→0.65, 0.75→0.72), meaning their co-movement is structural to the long end and largely independent of the NTB. This matters for trade construction: a 32s-vs-37s RV trade is exposed to long-end-specific risk, not to short-rate risk.

9 Technique 5 — Linear Regression (OLS)

9.1 Theory recap

Ordinary Least Squares regression models a continuous outcome as a linear combination of predictors plus i.i.d. Gaussian noise (Adi 2026, Ch. 9). Coefficients are interpreted as the expected change in the outcome per unit change in the predictor, holding other predictors constant. Standard diagnostics — residual plots, the Breusch–Pagan test for homoscedasticity, the Durbin–Watson test for serial correlation, and VIF for multicollinearity — must be inspected before any coefficient is treated as causal or actionable.

9.2 Business justification

A frequent client question: “if I run a long 30s position, how much 1yr NTB and 53s should I sell against it to be curve-neutral?” Regressing the 30s yield on those two curve points gives an empirical hedge ratio for that exact trade construction. The coefficients are not causal in the structural sense — but they are the right numbers to size a hedge today, given how the curve has co-moved over the last 100 trading days.

9.3 Implementation

Code
reg_df_r <- yields_long_r |>
  filter(instrument %in% c("1yr NTB","30s","53s")) |>
  select(date, instrument, yield) |>
  pivot_wider(names_from = instrument, values_from = yield) |>
  rename(NTB1Y = `1yr NTB`, y30s = `30s`, y53s = `53s`)

fit_r <- lm(y30s ~ NTB1Y + y53s, data = reg_df_r)

broom::tidy(fit_r, conf.int = TRUE) |>
  rename(Term = term, Estimate = estimate, SE = std.error,
         t = statistic, p = p.value,
         `CI low` = conf.low, `CI high` = conf.high) |>
  gt() |>
  fmt_number(columns = c(Estimate, SE, t, `CI low`, `CI high`), decimals = 4) |>
  fmt_scientific(columns = p, decimals = 2) |>
  tab_options(table.font.size = 11, column_labels.font.weight = "bold")
Term Estimate SE t p CI low CI high
(Intercept) −3.3109 1.1248 −2.9437 4.05 × 10−3 −5.5430 −1.0789
NTB1Y 0.6212 0.0515 12.0535 4.62 × 10−21 0.5189 0.7234
y53s 0.6590 0.0918 7.1759 1.39 × 10−10 0.4767 0.8412
Code

# Model-level statistics
glance_r <- broom::glance(fit_r)
glance_r |>
  select(`` = r.squared, `Adj R²` = adj.r.squared,
         `F` = statistic, `F p-value` = p.value,
         `Residual SE` = sigma, df = df.residual) |>
  gt() |>
  fmt_number(columns = c(``, `Adj R²`, `F`, `Residual SE`), decimals = 4) |>
  fmt_scientific(columns = `F p-value`, decimals = 2) |>
  tab_options(table.font.size = 11, column_labels.font.weight = "bold")
Adj R² F F p-value Residual SE df
0.8130 0.8092 213.0850 2.07 × 10−36 0.2671 98

OLS coefficients — 30s yield regressed on 1yr NTB and 53s yields

Code
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(fit_r, col = pal_mono, pch = 16, cex = 0.6)

OLS diagnostic panel — residuals vs fitted, Q-Q, scale-location, residuals vs leverage.
Code
par(mfrow = c(1, 1))
Code
bp <- lmtest::bptest(fit_r)
dw <- lmtest::dwtest(fit_r)
vif_vals <- car::vif(fit_r)

tibble::tibble(
  Test = c("Breusch–Pagan (heteroscedasticity)",
           "Durbin–Watson (serial correlation)",
           "VIF: NTB1Y", "VIF: 53s"),
  Statistic = c(sprintf("%.2f", bp$statistic),
                sprintf("%.2f", dw$statistic),
                sprintf("%.2f", vif_vals["NTB1Y"]),
                sprintf("%.2f", vif_vals["y53s"])),
  `p-value` = c(format.pval(bp$p.value, digits = 3, eps = .001),
                format.pval(dw$p.value, digits = 3, eps = .001),
                "—", "—")
) |>
  gt() |>
  tab_options(table.font.size = 11, column_labels.font.weight = "bold")
Formal diagnostic tests — heteroscedasticity, serial correlation, multicollinearity
Test Statistic p-value
Breusch–Pagan (heteroscedasticity) 1.07 0.585
Durbin–Watson (serial correlation) 0.63 <0.001
VIF: NTB1Y 1.45
VIF: 53s 1.45
Code
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson

reg_df_py = (yields_long_py
    .loc[yields_long_py["instrument"].isin(["1yr NTB","30s","53s"])]
    .pivot_table(index="date", columns="instrument", values="yield")
    .rename(columns={"1yr NTB":"NTB1Y","30s":"y30s","53s":"y53s"})
    .dropna())

X = sm.add_constant(reg_df_py[["NTB1Y","y53s"]])
y = reg_df_py["y30s"]

fit_py = sm.OLS(y, X).fit()

# Coefficient table
coef_tbl = pd.DataFrame({
    "Term":      fit_py.params.index,
    "Estimate":  fit_py.params.values,
    "SE":        fit_py.bse.values,
    "t":         fit_py.tvalues.values,
    "p":         fit_py.pvalues.values,
    "CI low":    fit_py.conf_int()[0].values,
    "CI high":   fit_py.conf_int()[1].values
})
(coef_tbl.style
   .hide(axis="index")
   .format({"Estimate":"{:.4f}","SE":"{:.4f}","t":"{:.4f}",
            "p":"{:.2e}","CI low":"{:.4f}","CI high":"{:.4f}"})
   .set_caption("OLS coefficients — Python statsmodels")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}]))
Table 5: OLS coefficients — Python statsmodels
OLS coefficients — Python (statsmodels) cross-check
Term Estimate SE t p CI low CI high
const -3.3109 1.1248 -2.9437 4.05e-03 -5.5430 -1.0789
NTB1Y 0.6212 0.0515 12.0535 4.62e-21 0.5189 0.7234
y53s 0.6590 0.0918 7.1759 1.39e-10 0.4767 0.8412
Code
fig, axes = plt.subplots(1, 2, figsize=(9, 3.6))

axes[0].scatter(fit_py.fittedvalues, fit_py.resid, alpha=0.6,
                color=PAL_MONO, s=18)
axes[0].axhline(0, color="grey", lw=0.6)
axes[0].set_xlabel("Fitted values")
axes[0].set_ylabel("Residuals")
axes[0].set_title("Residuals vs fitted", fontweight="bold", fontsize=11)

sm.qqplot(fit_py.resid, line="45", fit=True, ax=axes[1],
          markerfacecolor=PAL_MONO, markeredgecolor=PAL_MONO,
          markersize=4, alpha=0.7)
axes[1].set_title("Normal Q-Q plot", fontweight="bold", fontsize=11)

for ax in axes:
    sns.despine(ax=ax)
plt.tight_layout()
plt.show()

Python diagnostic — residuals vs fitted and Q-Q plot.
Code
bp_stat, bp_p, _, _ = het_breuschpagan(fit_py.resid, fit_py.model.exog)
dw_stat = durbin_watson(fit_py.resid)

# VIF
from statsmodels.stats.outliers_influence import variance_inflation_factor
vifs = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

diag = pd.DataFrame({
    "Test":["Breusch–Pagan p","Durbin–Watson stat",
            "VIF: NTB1Y","VIF: y53s","R²","Adj R²"],
    "Value":[f"{bp_p:.3e}", f"{dw_stat:.2f}",
             f"{vifs[1]:.2f}", f"{vifs[2]:.2f}",
             f"{fit_py.rsquared:.4f}", f"{fit_py.rsquared_adj:.4f}"]
})
(diag.style
   .hide(axis="index")
   .set_caption("Diagnostic tests — Python")
   .set_table_styles([{"selector":"caption",
                       "props":[("font-weight","bold"),("font-size","12px")]}]))
<pandas.io.formats.style.Styler object at 0x000002727085FD90>

9.4 Live hedge-ratio sensitivity panel

The OLS coefficients above are point estimates — under serial correlation the true uncertainty around them is wider than the conventional t-statistics suggest. Below is an interactive panel that lets you dial the two yield-betas (or your own model assumption) and see how the implied hedge notionals and the in-sample fit change in real time.

Code
import {Inputs} from "@observablehq/inputs"

viewof position_m = Inputs.number(
  {label: "Long 30s position (₦m)", value: 1000, step: 100, min: 0})
viewof beta_n = Inputs.range(
  [0, 1.5], {step: 0.01, value: ols_beta_n,
             label: "β to 1-year NTB (OLS = 0.62)"})
viewof beta_s = Inputs.range(
  [0, 1.5], {step: 0.01, value: ols_beta_s,
             label: "β to FGN 53s (OLS = 0.66)"})

// Reassemble the panel as an array of records
reg_data = reg_ntb.map((n, i) => ({ntb: n, y30: reg_y30[i], s53: reg_s53[i]}))

// In-sample fit under user betas (intercept refit to mean to keep apples-to-apples)
const_user = d3.mean(reg_data, d => d.y30 - beta_n*d.ntb - beta_s*d.s53)
sse_user   = d3.sum(reg_data, d => (d.y30 - (const_user + beta_n*d.ntb + beta_s*d.s53))**2)
sst        = d3.sum(reg_data, d => (d.y30 - d3.mean(reg_data, e => e.y30))**2)
rsq_user   = 1 - sse_user / sst

// Penalty vs OLS optimum
loss_pct = (rsq_user - ols_rsq) * 100

hedge_table = Inputs.table([
  {Instrument: "1-year NTB", "Yield-beta": beta_n.toFixed(2),
   "Short notional (₦m)": (beta_n * position_m).toFixed(0)},
  {Instrument: "FGN 53s",    "Yield-beta": beta_s.toFixed(2),
   "Short notional (₦m)": (beta_s * position_m).toFixed(0)},
], {height: 110, layout: "auto"})

md`
**Implied hedge structure for the position above:**

${hedge_table}

**In-sample R²** under your chosen betas: \`${rsq_user.toFixed(4)}\`
*(OLS optimum = \`${ols_rsq.toFixed(4)}\`, gap = \`${(loss_pct >= 0 ? "+" : "") + loss_pct.toFixed(2)} pp\`)*

How to read it · Move either slider away from the OLS point estimate and watch R² fall — the curvature of the response tells you how much the data "insists" on the OLS coefficients vs. how flat the likelihood is. If your view is that the true betas are different (e.g. because you expect the 37s episode to re-anchor cross-bucket relationships), you can read the implied hedge notionals directly off the table above for any position size.
`

9.5 Interpretation

Tip📐 Business interpretation — hedge ratios you can quote

Reading the coefficient table: a 1.00 percentage-point rise in the 1yr NTB yield is associated with a 0.62 pp rise in the 30s yield (front-end yield-beta of 0.62); a 1.00 pp rise in the 53s yield is associated with a 0.66 pp rise in the 30s yield (ultra-long yield-beta of 0.66). The intercept (–3.31) absorbs the constant gap between the 30s and its two reference points and has no standalone economic interpretation. R² ≈ 0.81 — together these two curve points explain about four-fifths of the 30s’ daily yield variation. Both regressors are statistically significant at p < 0.001 by the conventional standard errors, with VIFs ≈ 1.45 (no multicollinearity concern).

Translated for a non-technical PM: “In yield-level terms over the last 100 trading days, the 30s’ daily move is best described as 0.62× the 1-year NTB’s move plus 0.66× the 53s’ move. About 81% of the 30s’ daily yield variation is explained by those two curve points; the remaining 19% is the idiosyncratic move that defines a directional bet on the 30s itself. The two coefficients are the empirical yield-betas you should plug into a hedge calculator before turning them into notional DV01-matched hedges.”

Warning⚠️ Diagnostic caveats — what to flag honestly to the credit committee
  1. Breusch–Pagan p ≈ 0.59 — fails to reject homoscedasticity. Good news: residual variance does not visibly widen during the rate-rally episodes; the constant-variance assumption is intact. No robust-SE adjustment is needed on this count.
  2. Durbin–Watson = 0.63 — strong positive serial correlation in the residuals. This is the expected pathology of regressing daily yield levels (which are near-unit-root) rather than first-differences. The OLS point estimates are still unbiased and consistent (Gauss–Markov holds for unbiasedness without independence), but the conventional standard errors understate the true uncertainty. Implication: trust the signs and magnitudes of the coefficients as empirical yield-betas; treat the t-statistics as upper bounds on real significance. A formal fix is Newey–West HAC standard errors with a 5-day bandwidth (flagged in Limitations).

10 Integrated Findings

The five analyses converge on one coherent picture: the FGN curve is humped, with an inverted ultra-long; the front-end has rallied as the easing cycle was priced in; the 37s has cheapened against every neighbour; and the belly is the structurally stable centre of gravity. Pulling the strands together:

  1. EDA revealed the humped shape (Belly mean 16.63 % > Short 16.43 % > Long 16.33 % > Ultra-Long 14.73 %) and showed volatility is front-end concentrated (1yr NTB σ = 0.63 vs 53s σ = 0.35).
  2. Visualisation made the +61 bp 37s selloff visible at a glance — it is the only bond that materially cheapened while every other instrument rallied or held flat.
  3. Hypothesis testing confirmed that (a) tenor buckets carry statistically distinct mean yields (ANOVA F = 530, η² = 0.57) — the curve is real, not noise; and (b) the 1-year NTB shifted 135 bp lower between January and May (Cohen’s d = 2.97) — a decisive, not marginal, move.
  4. Correlation showed that the 1-year NTB and the 37s are nearly decoupled (r = 0.33) — the lowest pair in the matrix — confirming that the 37s’ selloff is genuinely idiosyncratic, not a curve-level echo. Belly bonds remain tightly co-moving (r > 0.94 within), so RV inside the belly is dead.
  5. Regression delivered the actionable hedge ratios: 30s yield-beta of 0.62 to the NTB and 0.66 to the 53s (R² = 0.81). These betas convert a directional belly view into a curve-neutral structure.
Important🎯 Single trade recommendation to the desk for Monday’s client mail-out

Pitch a long 32s / short 37s relative-value trade to PFA and asset-manager clients. The thesis is mean-reversion of the 37s back toward its belly neighbours, supported by (i) the unique +61 bp selloff in the 37s, (ii) the 0.65–0.75 partial correlation between belly and long-end (so the two legs do move together enough to hedge cleanly), and (iii) the regression-implied yield-beta of 0.5–0.6 between belly and long instruments. For a real-money client who prefers expressing the same view as a single instrument, recommend buying the 37s outright — it is the cheapest curve point relative to its statistical fair value.

11 Limitations & Further Work

  • Levels regression, not changes. Yields are non-stationary over 5 months; serially correlated residuals inflate spurious-significance risk. With more time I would re-run the regression on first-differences (Δyield), report Newey–West HAC standard errors, and run an Engle–Granger cointegration test to formally distinguish a structural relationship from a spurious one.
  • Window length. 101 trading days is enough for stable cross-sectional correlations but short for stress-period inference. Extending to two years would let me condition the correlation matrix on MPC-meeting versus non-MPC windows.
  • Liquidity weighting. All quotes are mid-prices irrespective of traded volume. Future iterations will weight observations by traded notional pulled from FMDQ end-of-day reports.
  • Missing macro covariates. The regression omits oil price, CBN OMO auction stop-rates and Naira/USD — all known curve drivers. Adding them would lift R² and turn the regression into a meaningful attribution rather than a hedge-ratio tool.
  • No formal causality test. Granger causality between the 1yr NTB and the 53s would test whether front-end moves lead ultra-long moves (a question with direct value for execution sequencing).

References

All references in APA 7th edition.

Course textbook

  • Adi, B. (2026). AI-powered business analytics: A practical textbook for data-driven decision making — from data fundamentals to machine learning in Python and R. Lagos Business School / markanalytics.online. https://markanalytics.online/ai-powered-data-analytics/
    • Ch. 4 — Exploratory Data Analysis (Technique 1).
    • Ch. 5 — Data Visualisation (Technique 2).
    • Ch. 6 — Hypothesis Testing (Technique 3).
    • Ch. 8 — Correlation Analysis (Technique 4).
    • Ch. 9 — Ordinary Least Squares Regression (Technique 5).

Methodological sources

  • Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum.
  • Fisher, R. A. (1925). Statistical methods for research workers. Oliver and Boyd.
  • Neyman, J., & Pearson, E. S. (1933). On the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society A, 231(694–706), 289–337. https://doi.org/10.1098/rsta.1933.0009
  • Schönbrodt, F. D., & Perugini, M. (2013). At what sample size do correlations stabilize? Journal of Research in Personality, 47(5), 609–612. https://doi.org/10.1016/j.jrp.2013.05.009
  • Tukey, J. W. (1977). Exploratory data analysis. Addison-Wesley.
  • Wickham, H. (2016). ggplot2: Elegant graphics for data analysis (2nd ed.). Springer. https://doi.org/10.1007/978-3-319-24277-4
  • Wilkinson, L. (2005). The grammar of graphics (2nd ed.). Springer.

Software and packages

  • R Core Team. (2024). R: A language and environment for statistical computing (Version 4.x). R Foundation for Statistical Computing. https://www.R-project.org/
  • Van Rossum, G., & Drake, F. L. (2009). Python 3 reference manual. CreateSpace.
  • Allaire, J. J., Teague, C., Scheidegger, C., Xie, Y., & Dervieux, C. (2022). Quarto (Version 1.x) [Computer software]. https://doi.org/10.5281/zenodo.5960048
  • Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686
  • Iannone, R., Cheng, J., Schloerke, B., Hughes, E., Lauer, A., Seo, J., Brevoort, K., & Roy, O. (2024). gt: Easily create presentation-ready display tables [R package]. https://gt.rstudio.com/
  • Robinson, D., Hayes, A., & Couch, S. (2024). broom: Convert statistical objects into tidy tibbles [R package]. https://broom.tidymodels.org/
  • Kassambara, A. (2023). ggcorrplot: Visualization of a correlation matrix using ggplot2 [R package]. https://CRAN.R-project.org/package=ggcorrplot
  • Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. (Package: lmtest.)
  • Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3rd ed.). Sage. (Package: car.)
  • Wickham, H., & Bryan, J. (2023). readxl: Read Excel files [R package]. https://readxl.tidyverse.org/
  • McKinney, W. (2010). Data structures for statistical computing in Python. Proceedings of the 9th Python in Science Conference, 56–61. https://doi.org/10.25080/Majora-92bf1922-00a (Package: pandas.)
  • Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., … Oliphant, T. E. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2
  • Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2
  • Seabold, S., & Perktold, J. (2010). statsmodels: Econometric and statistical modeling with Python. Proceedings of the 9th Python in Science Conference, 92–96.
  • Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in Science & Engineering, 9(3), 90–95. https://doi.org/10.1109/MCSE.2007.55
  • Waskom, M. L. (2021). seaborn: Statistical data visualization. Journal of Open Source Software, 6(60), 3021. https://doi.org/10.21105/joss.03021
  • Gazoni, E., & Clark, C. (2024). openpyxl: A Python library to read/write Excel 2010 xlsx/xlsm files [Software]. https://openpyxl.readthedocs.io/

Data source

  • Dutire, S. (2026). Daily closing yields, FGN sovereign curve, 5 January – 25 May 2026 [Primary dataset]. Compiled from a Nigerian inter-dealer brokerage closing run-sheet, Lagos, Nigeria. Available on request from the author.

Appendix — AI Usage Statement

I used Anthropic’s Claude (Opus 4.7) as a coding assistant to help scaffold the Quarto document’s section structure, format presentation-grade gt and pandas-Styler tables, and translate equivalent analyses between R and Python in the panel-tabset layout. Every analytical decision is my own: the choice of Case Study 1 over the alternatives, the tenor-bucket scheme, the two specific hypothesis tests, the choice of the 1yr NTB and 53s as the regressors for the 30s, the interpretation of each result, and the integrated trade recommendation reflect my judgement as a fixed-income sales trader and broker on a live inter-dealer brokerage desk. Where the AI proposed an interpretation, I rewrote it in language consistent with the daily client notes I personally publish.