import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from polars_bloomberg import BQuery
from tulip.data.dataquery import get_us_otr_yields
from tulip.analysis.return_functions import calculate_excess_returns
import os
from dotenv import load_dotenv

load_dotenv()
has_terminal = os.environ.get("LOCAL_BLOOMBERG_LICENSE") == "true"

Sizing a Steepener Trade¶

Using a 5s30s Steepener Trade as Example

Calculating it with PCA on OTR Yields (Fixed-Income Native approach)¶

Define the front leg (10Y) and back leg (30Y) futures contracts for the steepener trade.

FRONT_LEG = "FVZ5 Comdty"
BACK_LEG = "WNZ5 Comdty"
FRONT_LEG_GENERIC = "FV1 Comdty"
BACK_LEG_GENERIC = "WN1 Comdty"

HALF_LIFE = 75

Fetch cheapest-to-deliver (CTD) bond information including conversion factors, durations, and contract values.

ctd_identifier = "FUT_CTD_TICKER"
with BQuery() as bq:
    ctd_data = bq.bdp(
        [FRONT_LEG, BACK_LEG],
        [
            ctd_identifier,
            "FUT_CNVS_FACTOR",
            "FUT_CNV_RISK_FRSK",
            "CONTRACT_VALUE",
            "PX_LAST",
        ],
    )
ctd_identifiers_df = [
    f"{x} Govt" for x in ctd_data.get_column(ctd_identifier).to_list()
]
cf = ctd_data.to_pandas().set_index("FUT_CTD_TICKER")[
    "FUT_CNVS_FACTOR"
]  # Conversion Factor
cf.index = ctd_identifiers_df

mod_dur = ctd_data.to_pandas().set_index("security")[
    "FUT_CNV_RISK_FRSK"
]  # Risk to determine hedge ratio
contract_value = ctd_data.to_pandas().set_index("security")[
    "CONTRACT_VALUE"
]  # Contract Value
ctd_px = ctd_data.to_pandas().set_index("security")["PX_LAST"]
#

Set up utility functions for PCA analysis and load US Treasury par curve data across all key rate durations.

# ---------- Utils ----------
def ew_weights(dates, half_life_days=HALF_LIFE):
    dates = pd.to_datetime(dates)
    age = (dates.max() - dates).days.values
    lam = np.log(2) / half_life_days
    w = np.exp(-lam * age)
    return w / w.sum()


def pca_from_changes(X_bp, K_keep=3, weights=None):
    X = X_bp.to_numpy(dtype=float)
    if weights is None:
        Xc = X - X.mean(axis=0, keepdims=True)
        S = (Xc.T @ Xc) / X.shape[0]
    else:
        w = np.asarray(weights, float)
        w = w / w.sum()
        # weighted mean
        mu = (w[:, None] * X).sum(axis=0, keepdims=True)
        Xc = X - mu
        # avoid np.diag(w): O(T^2) memory
        S = Xc.T @ (Xc * w[:, None])  # bp^2; weights already normalized
    eigvals, eigvecs = np.linalg.eigh(S)
    idx = eigvals.argsort()[::-1]
    lam = eigvals[idx][:K_keep]
    B = eigvecs[:, idx][:, :K_keep]
    if B[:, 0].sum() < 0:
        B[:, 0] *= -1  # stable PC1 sign
    Omega = np.diag(lam)
    return S, B, Omega, eigvals[idx]


def normalize_tenor_cols(cols):
    out = []
    for c in cols:
        if isinstance(c, (int, float)):
            out.append(float(c))
            continue
        s = str(c).strip().upper()
        if s.endswith("MO"):
            out.append(float(s[:-2]) / 12.0)
        elif s.endswith("M"):
            out.append(float(s[:-1]) / 12.0)
        elif s.endswith("YR"):
            out.append(float(s[:-2]))
        elif s.endswith("Y"):
            out.append(float(s[:-1]))
        else:
            out.append(float(s))
    return out


# ---------- Inputs ----------
KRD_BUCKET_TO_YEARS = {
    "KEY_RATE_DUR_3MO": 0.25,
    "KEY_RATE_DUR_6MO": 0.5,
    "KEY_RATE_DUR_1YR": 1,
    "KEY_RATE_DUR_2YR": 2,
    "KEY_RATE_DUR_3YR": 3,
    "KEY_RATE_DUR_4YR": 4,
    "KEY_RATE_DUR_5YR": 5,
    "KEY_RATE_DUR_6YR": 6,
    "KEY_RATE_DUR_7YR": 7,
    "KEY_RATE_DUR_8YR": 8,
    "KEY_RATE_DUR_9YR": 9,
    "KEY_RATE_DUR_10YR": 10,
    "KEY_RATE_DUR_15YR": 15,
    "KEY_RATE_DUR_20YR": 20,
    "KEY_RATE_DUR_25YR": 25,
    "KEY_RATE_DUR_30YR": 30,
}

# Measure of price sensitivity of the security to rate shifts at the X year point on the curve.
# The curve is shifted by 1 basis point and the result is scaled to a 100 basis point curve shift.
# The curve used is the curve for the selected bond's currency.

krd_fields = list(KRD_BUCKET_TO_YEARS.keys())

# ---------- Par curve (decimals) ---------------
# ===============================================

Y_par = (
    get_us_otr_yields(start_date=pd.Timestamp("1989-01-01")).astype(float) / 100
)  # in percentage points
Y_par.columns = normalize_tenor_cols(Y_par.columns)
Y_par = Y_par.sort_index(axis=1).dropna()

Y_par = Y_par.reindex(list(KRD_BUCKET_TO_YEARS.values()), axis=1).interpolate(
    axis=1, method="cubicspline", limit_area="inside"
)

Calculate PCA-based risk model: compute key rate durations, run PCA on yield changes, and derive the DV01-neutral hedge ratio and expected volatility.

# ---------- Data on the futures ----------------
# ===============================================

ctd_identifier = "FUT_CTD_TICKER"
with BQuery() as bq:
    fut_df = (
        bq.bdp(
            [FRONT_LEG, BACK_LEG],
            [ctd_identifier, "FUT_CNVS_FACTOR", "CONTRACT_VALUE", ""],
        )
        .to_pandas()
        .set_index("security")
    )

    future_history_generic_futures = (
        bq.bdh(
            [FRONT_LEG_GENERIC, BACK_LEG_GENERIC],
            ["CONTRACT_VALUE"],
            start_date=pd.Timestamp("2020-01-01"),
            end_date=pd.Timestamp.today(),
        )
        .to_pandas()
        .set_index(["security", "date"])
    )


fut_to_ctd = dict(zip(fut_df.index, [f"{x} Govt" for x in fut_df[ctd_identifier]]))

fut_to_ctd

CF_by_leg = fut_df["FUT_CNVS_FACTOR"]  # index = futures legs
CONTRACT_VAL_by_leg = fut_df["CONTRACT_VALUE"]  # index = futures legs

# CTD prices
with BQuery() as bq:
    px_df = (
        bq.bdp(list(fut_to_ctd.values()), ["PX_LAST"])
        .to_pandas()
        .set_index("security")["PX_LAST"]
    )

# CTD key-rate durations (years per $100)
with BQuery() as bq:
    krd_ctd_raw = (
        bq.bdp(list(fut_to_ctd.values()), krd_fields).to_pandas().set_index("security")
    )

# Rename KRD columns â†’ numeric years, sorted
krd_ctd = krd_ctd_raw.copy()
krd_ctd.columns = [KRD_BUCKET_TO_YEARS[c] for c in krd_ctd.columns]
krd_ctd = krd_ctd.sort_index(axis=1)

# ---------- 1) PCA on par changes in bp ----------
X_bp = (Y_par.diff().dropna() * 1e4).astype(float)  # bp
# overlap tenors across curve and KRD buckets
tenors = sorted(set(krd_ctd.columns).intersection(X_bp.columns))
assert len(tenors) >= 3, f"Too few overlapping buckets: {tenors}"
X_bp = X_bp[tenors]

w = ew_weights(X_bp.index, HALF_LIFE)
S, B_arr, Omega, eigvals_desc = pca_from_changes(X_bp, K_keep=3, weights=w)
B = pd.DataFrame(
    B_arr, index=tenors, columns=[f"PC{k + 1}" for k in range(Omega.shape[0])]
)


def durations_to_ctd_kr01_per_contract(dur_row: pd.Series, contract_value):
    # KRD * 1000 / 100 = KRD * 10 for 100k contracts
    return dur_row.astype(float) * float(contract_value) / 10000.0


rows = []
for fut_leg, ctd in fut_to_ctd.items():
    dur = krd_ctd.loc[ctd].reindex(tenors).fillna(0.0)
    # dur units: $ per $100 face per 100bp (raw Bloomberg KRDs)
    kr01_ctd = dur / 100.0  # CTD KRDs already per contract, just convert to per bp
    # kr01_ctd units: $ per $100 face per 1bp (converted from 100bp to 1bp)
    # True DV01
    rows.append(kr01_ctd.rename(fut_leg))

KR01_CTD = pd.DataFrame(rows)  # index = futures legs, cols = tenors
# KR01_CTD units: $ per $100 face per 1bp
KR01_FUT = KR01_CTD.div(CF_by_leg, axis=0)  # divide by CF to get futures KR01
# KR01_FUT units: $ per $100 face per 1bp for futures (CTD risk â†’ futures risk)
KR01_FUT = KR01_FUT.mul(CONTRACT_VAL_by_leg, axis=0) / 100
# KR01_FUT final units: $ per contract per 1bp
# (scaled from $100 face to full contract size, e.g., $100,000)


# sanity: PV01 tie-out
pv01 = KR01_FUT.sum(axis=1)
assert np.all(pv01.values > 0), "PV01 must be positive per leg"
assert np.all(np.linalg.eigvalsh(S) >= -1e-10)

# ---------- 3) Instrument â†’ factor and covariance ----------
K_mat = KR01_FUT[tenors].to_numpy()  # [2 x J], units: $/bp per contract
L = K_mat @ B.to_numpy()  # [2 x K], units: $/bp per factor per contract
Sigma_inst = L @ Omega @ L.T  # [2 x 2], units: $Â² per contractÂ² (variance)

# ---------- 4) DV01-neutral hedge and sigma ----------
h = pv01.loc[FRONT_LEG] / pv01.loc[BACK_LEG]
x_dn = np.array([1.0, -float(h)])
sigma_daily = float(np.sqrt(x_dn @ Sigma_inst @ x_dn))
sigma_annual = sigma_daily * np.sqrt(252)

print("Tenors used:", tenors)
print("PV01 per contract ($/bp):")
print(pv01.round(2).to_string())
print("DV01-neutral h =", round(float(h), 6))
print("Daily sigma ($):", round(sigma_daily, 2))
print("Annualized sigma ($):", round(sigma_annual, 2))

# Optional diagnostics
explained = float(np.sum(np.diag(Omega)) / np.sum(np.linalg.eigvalsh(S)))
print("Variance explained by kept PCs:", round(explained, 4))

Y_hist = Y_par
# Daily yield changes in bp
dY_bp = Y_hist.diff().dropna() * 1e4

## Using Actuals:
# Use fixed KR01s from a recent date (or average over recent period)
# These are the same loadings from your PCA calculation
K_front = KR01_FUT.loc[FRONT_LEG, :].values  # $/bp per tenor
K_back = KR01_FUT.loc[BACK_LEG, :].values

# Fixed DV01-neutral hedge ratio
h_fixed = K_front.sum() / K_back.sum()

# Calculate daily PnL: sum(KR01_i * dY_i) for each leg
pnl_front = (dY_bp * K_front).sum(axis=1)  # $ from front leg
pnl_back = (dY_bp * K_back).sum(axis=1)  # $ from back leg

# DV01-neutral portfolio PnL
pnl_steepener = pnl_front - h_fixed * pnl_back

Checks¶

Check of Durations and Key Rate Durations

kr01_sum = krd_ctd_raw.sum(axis=1).div(CF_by_leg.values)
kr01_sum.index = [FRONT_LEG, BACK_LEG]  # Sum across tenors
print(f"Front leg - Sum of KR01s: ${kr01_sum[FRONT_LEG]:.2f}/bp")
print(f"Back leg - Sum of KR01s: ${kr01_sum[BACK_LEG]:.2f}/bp")

# Check 2: Compare to Bloomberg's total risk field
with BQuery() as bq:
    bbg_dv01 = (
        bq.bdp([FRONT_LEG, BACK_LEG], ["CONVENTIONAL_CTD_FORWARD_FRSK"])
        .to_pandas()
        .set_index("security")
    ).squeeze()

print(f"Front leg - Sum of KR01s: ${bbg_dv01[FRONT_LEG]:.2f}/bp")
print(f"Back leg - Sum of KR01s: ${bbg_dv01[BACK_LEG]:.2f}/bp")

Backtests¶

The PCA-predicted volatility by computing realized P&L using historical yield changes multiplied by today’s KR01 sensitivities¶

This isolates yield curve risk only (excluding futures basis risk and CTD switches) and validates whether our PCA decomposition accurately captures the covariance structure of yield changes - if the ratio is close to 1.0, it confirms the PCA model is well-calibrated.

Y_hist = Y_par
# Daily yield changes in bp
dY_bp = Y_hist.diff().dropna() * 1e4

# Use fixed KR01s from a recent date (or average over recent period)
# These are the same loadings from your PCA calculation
K_front = KR01_FUT.loc[FRONT_LEG, tenors].values  # $/bp per tenor
K_back = KR01_FUT.loc[BACK_LEG, tenors].values

# Fixed DV01-neutral hedge ratio
h_fixed = K_front.sum() / K_back.sum()

# Calculate daily PnL: sum(KR01_i * dY_i) for each leg
pnl_front = (dY_bp * K_front).sum(axis=1)  # $ from front leg
pnl_back = (dY_bp * K_back).sum(axis=1)  # $ from back leg

# DV01-neutral portfolio PnL
pnl_steepener = pnl_front - h_fixed * pnl_back

# Realized vol
sigma_realized_daily = pnl_steepener.std()
sigma_realized_annual = sigma_realized_daily * np.sqrt(252)

# Compare
print(f"PCA-predicted daily vol: ${sigma_daily:.2f}")
print(f"Realized daily vol (yields only): ${sigma_realized_daily:.2f}")
print(f"Ratio (realized/predicted): {sigma_realized_daily / sigma_daily:.3f}")

# Optional: rolling window comparison
rolling_vol = pnl_steepener.rolling(63).std() * np.sqrt(252)  # 3M rolling

## Using Generics
KR01_FUT_GEN = KR01_FUT.copy()
KR01_FUT_GEN.index = [FRONT_LEG_GENERIC, BACK_LEG_GENERIC]
K_front_gen = KR01_FUT_GEN.loc[FRONT_LEG_GENERIC, :].values  # $/bp per tenor
K_back_gen = KR01_FUT_GEN.loc[BACK_LEG_GENERIC, :].values
pnl_front_gen = (dY_bp * K_front_gen).sum(axis=1)  # $ from front leg
pnl_back_gen = (dY_bp * K_back_gen).sum(axis=1)
pnl_steepener_gen = pnl_front_gen - h_fixed * pnl_back_gen
sigma_realized_gen_daily = pnl_steepener_gen.std()
sigma_realized_gen_annual = sigma_realized_gen_daily * np.sqrt(252)

pnl_steepener_gen.std()
# pnl_steepener_gen.loc['20':].mean()

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

pnl = pnl_steepener_gen.dropna()

# ============================================================
# Empirical (historical) VaR / ES
# ============================================================
var95_raw = np.percentile(pnl, 5)  # 5th percentile
var99_raw = np.percentile(pnl, 1)  # 1st percentile

es95_raw = pnl[pnl <= var95_raw].mean()
es99_raw = pnl[pnl <= var99_raw].mean()

# Convert to positive "loss" numbers if you want
var95 = -var95_raw
var99 = -var99_raw
es95 = -es95_raw
es99 = -es99_raw

# ============================================================
# Normal-assumption VaR / ES (fit μ, σ to pnl)
# ============================================================
mu = pnl.mean()
sigma = pnl.std(ddof=1)

alpha95, alpha99 = 0.05, 0.01
z95, z99 = norm.ppf(alpha95), norm.ppf(alpha99)

# VaR under Normal(μ, σ²)
var95_norm_raw = mu + sigma * z95
var99_norm_raw = mu + sigma * z99

# ES under Normal(μ, σ²): ES_α = μ - σ * φ(zα) / α
es95_norm_raw = mu - sigma * norm.pdf(z95) / alpha95
es99_norm_raw = mu - sigma * norm.pdf(z99) / alpha99

var95_norm = -var95_norm_raw
var99_norm = -var99_norm_raw
es95_norm = -es95_norm_raw
es99_norm = -es99_norm_raw


print("Empirical Mean and Vol")
print(f"Mean : {mu:.4f}, Vol : {sigma:.4f}")
print(f"Last five years: {pnl_steepener_gen.loc['2020':].std():.4f}")

print("Empirical (historical)")
print(f"95% VaR: {var95:.4f}, 99% VaR: {var99:.4f}")
print(f"95% ES : {es95:.4f}, 99% ES : {es99:.4f}")

print("\nNormal-fit")
print(f"95% VaR: {var95_norm:.4f}, 99% VaR: {var99_norm:.4f}")
print(f"95% ES : {es95_norm:.4f}, 99% ES : {es99_norm:.4f}")

# ============================================================
# Plot histogram with both empirical and normal VaR lines
# ============================================================
plt.figure(figsize=(10, 5))
plt.hist(pnl, bins=50, alpha=0.7)

# empirical
plt.axvline(var95_raw, linestyle="--", label="Empirical 95% VaR")
plt.axvline(var99_raw, linestyle=":", label="Empirical 99% VaR")

# normal-assumption
plt.axvline(var95_norm_raw, linestyle="--", color="red", label="Normal 95% VaR")
plt.axvline(var99_norm_raw, linestyle=":", color="red", label="Normal 99% VaR")

plt.title("Daily P&L Distribution: Empirical vs Normal VaR")
plt.xlabel("Daily P&L")
plt.ylabel("Frequency")
plt.legend()
plt.tight_layout()
plt.show()

Another Check: The actual dollar movement of the PNL:

note: The PCA number should typically be lower than historical realized, as it excludes basis volatility and CTD switch risk.

with BQuery() as bq:
    future_history_actual_futures = (
        bq.bdh(
            [FRONT_LEG, BACK_LEG],
            ["CONTRACT_VALUE", "YLD_YTM_MID", "CONVENTIONAL_CTD_FORWARD_FRSK"],
            start_date=pd.Timestamp("2020-01-01"),
            end_date=pd.Timestamp.today(),
        )
        .to_pandas()
        .set_index(["security", "date"])
    )

future_history_actual_futures = future_history_actual_futures.unstack(
    "security"
).ffill()

# Assume the correct 'hedge' ratio is applied
print("\nAssume the correct 'hedge' ratio is applied")
print("=" * 40)
H = future_history_actual_futures.xs(
    "CONVENTIONAL_CTD_FORWARD_FRSK", level=0, axis=1
).dropna()
H = (H[FRONT_LEG] / H[BACK_LEG]).shift(1)
PNL = (
    future_history_actual_futures.xs("CONTRACT_VALUE", level=0, axis=1).dropna().diff()
)
PNL[BACK_LEG] = PNL[BACK_LEG] * H * -1
PNL = PNL.dropna().sum(axis=1)
print("Average Daily PNL in dollars:", round(PNL.pow(2).mean() ** 0.5, 2))
print(
    "Average EWM Daily PNL in dollars:",
    round(PNL.pow(2).ewm(halflife=HALF_LIFE).mean().iloc[-1] ** 0.5),
)

print("\n")

print("Assume a fixed hedge ratio of 0.40")
print("=" * 40)
PNL = (
    future_history_actual_futures.xs("CONTRACT_VALUE", level=0, axis=1).dropna().diff()
)
PNL[BACK_LEG] = PNL[BACK_LEG] * 0.40 * -1
PNL = PNL.dropna().sum(axis=1)
print("Average Daily PNL in dollars:", round(PNL.pow(2).mean() ** 0.5, 2))
print(
    "Average EWM Daily PNL in dollars:",
    round(PNL.pow(2).ewm(halflife=HALF_LIFE).mean().iloc[-1] ** 0.5),
)

Checks: Using the actual PNL from the generics

with BQuery() as bq:
    future_history_generic_futures = (
        bq.bdh(
            [FRONT_LEG_GENERIC, BACK_LEG_GENERIC],
            ["CONTRACT_VALUE"],
            start_date=pd.Timestamp("2020-01-01"),
            end_date=pd.Timestamp.today(),
        )
        .to_pandas()
        .set_index(["security", "date"])
    )

future_history_generic_futures = future_history_generic_futures.unstack(
    "security"
).ffill()

Using Contract Value (Issues with rolls)

# Assume a
PNL = (
    future_history_generic_futures.xs("CONTRACT_VALUE", level=0, axis=1).dropna().diff()
)
PNL[BACK_LEG_GENERIC] = PNL[BACK_LEG_GENERIC] * 0.40 * -1
PNL = PNL.dropna().sum(axis=1)
print("Average Daily PNL in dollars:", round(PNL.pow(2).mean() ** 0.5, 2))
print(
    "Average EWM Daily PNL in dollars:",
    round(PNL.pow(2).ewm(halflife=HALF_LIFE).mean().iloc[-1] ** 0.5),
)

With Roll Adjustment to Prices (should account for the roll)

Note maybe we should use this logic for calculate_excess

with BQuery() as bq:
    future_history_generic_futures = (
        bq.bdh(
            [FRONT_LEG_GENERIC, BACK_LEG_GENERIC],
            ["PX_LAST", "CONTRACT_VALUE"],
            start_date=pd.Timestamp("2020-01-01"),
            end_date=pd.Timestamp.today(),
        )
        .to_pandas()
        .set_index(["security", "date"])
    )

    cont_size = (
        bq.bdp(
            [FRONT_LEG_GENERIC, BACK_LEG_GENERIC],
            ["FUT_CONT_SIZE"],
        )
        .to_pandas()
        .set_index("security")["FUT_CONT_SIZE"]
    )

future_history_generic_futures = future_history_generic_futures.unstack(
    "security"
).ffill()

roll_adj = (
    future_history_generic_futures.xs("PX_LAST", level=0, axis=1)
    .mul(cont_size)
    .div(100)
)
PNL = roll_adj.dropna().diff()
PNL[BACK_LEG_GENERIC] = PNL[BACK_LEG_GENERIC] * 0.40 * -1
PNL = PNL.dropna().sum(axis=1)
print("Average Daily PNL in dollars:", round(PNL.pow(2).mean() ** 0.5, 2))
print(
    "Average EWM Daily PNL in dollars:",
    round(PNL.pow(2).ewm(halflife=HALF_LIFE).mean().iloc[-1] ** 0.5),
)

Brian’s check (a bit above)

with BQuery() as bq:
    history = (
        bq.bdh(
            ["GT10 Govt", "GT30 Govt"],  # , "GT20 Govt"
            ["YLD_YTM_MID"],  # "PX_LAST",
            start_date=pd.Timestamp("2020-01-01"),
            end_date=pd.Timestamp.today(),
        )
        .to_pandas()
        .set_index(["security", "date"])
        .squeeze()
        .unstack(level=0)
    )


spread_chg = history.loc[:, "GT10 Govt"] - history.loc[:, "GT30 Govt"]
spread_chg_bps = spread_chg * 100
spread_chg_bps.pow(2).ewm(halflife=HALF_LIFE).mean().iloc[-1] ** 0.5
# moves 24 bps daily
print(
    "Avg Daily Move",
    ctd_data.to_pandas()["FUT_CNV_RISK_FRSK"][0] * spread_chg_bps.std(),
)

How the risk system would account for this position:¶

from tulip.risk.models.cov_estimators import *

kate_fast = (15, 15 * 2)
kate_slow = (126, 126 * 2)

excess_return = calculate_excess_returns([FRONT_LEG_GENERIC, BACK_LEG_GENERIC])
CONTRACT_VAL_front, CONTRACT_VAL_back = CONTRACT_VAL_by_leg.to_list()
# Use GROSS notional as imaginary NAV
NAV = abs(CONTRACT_VAL_front) + abs(h * CONTRACT_VAL_back)
# Weights (now both will be ~0.5 in magnitude)
weights = np.array(
    [
        CONTRACT_VAL_front / NAV,  # ~0.48
        -h * CONTRACT_VAL_back / NAV,
    ]
)  # ~-0.52

# Fast
kate_ewm_fast_cov = ewm_covariance(
    ex_ret=excess_return,
    hl_vola=kate_fast[0],
    hl_corr=kate_fast[1],
    give_last_date=True,
)

cov_ann_fast = (kate_ewm_fast_cov["cov"]) / 252  # It comes already annualized

# Checks
print(f"Covariance diagonal (variance) front leg: {cov_ann_fast.iloc[0, 0]:.6f}")
print(f"Daily std dev front leg: {np.sqrt(cov_ann_fast.iloc[0, 0]):.6f}")
print(f"Annualized vol: {np.sqrt(cov_ann_fast.iloc[0, 0] * 252):.4f}")

print(f"Weight front: {weights[0]:.4f}")
print(f"Weight back: {weights[1]:.4f}")
print(f"Sum of abs weights: {abs(weights).sum():.4f}")  # = 1.0

# Calculate risk
portfolio_vol = np.sqrt(weights @ cov_ann_fast @ weights.T)
dollar_vol = portfolio_vol * NAV
print(f"Portfolio daily vol ($): {dollar_vol:.2f}")

# Slow
kate_ewm_slow_cov = ewm_covariance(
    ex_ret=excess_return,
    hl_vola=kate_slow[0],
    hl_corr=kate_slow[1],
    give_last_date=True,
)

cov_ann_slow = (kate_ewm_slow_cov["cov"]) / 252  # It comes already annualized

# Checks
print(f"Covariance diagonal (variance) front leg: {cov_ann_slow.iloc[0, 0]:.6f}")
print(f"Daily std dev front leg: {np.sqrt(cov_ann_slow.iloc[0, 0]):.6f}")
print(f"Annualized vol: {np.sqrt(cov_ann_slow.iloc[0, 0] * 252):.4f}")

print(f"Weight front: {weights[0]:.4f}")
print(f"Weight back: {weights[1]:.4f}")
print(f"Sum of abs weights: {abs(weights).sum():.4f}")  # = 1.0

# Calculate risk
portfolio_vol_slow = np.sqrt(weights @ cov_ann_slow @ weights.T)
dollar_vol_slow = portfolio_vol_slow * NAV
print(f"Portfolio daily vol ($): {dollar_vol_slow:.2f}")

midpoint_daily_vol = (dollar_vol_slow + dollar_vol) / 2
print(f"Midpoint daily vol ($): {midpoint_daily_vol:.2f}")

from tulip_mania.notebook_related import notebook_updated

notebook_updated()