!pip install numpy pandas matplotlib seaborn scipy yfinance --quiet

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats, optimize
import warnings

warnings.filterwarnings('ignore')
%matplotlib inline
plt.rcParams['figure.figsize'] = (14, 5)
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False
print('Imports ready.')

# ── Data Source Toggle ─────────────────────────────
USE_LIVE_DATA = False
TICKERS       = ['BTC-USD','ETH-USD','SOL-USD','BNB-USD','AVAX-USD',
                  'MATIC-USD','LINK-USD','DOT-USD']
START_DATE    = '2020-01-01'
END_DATE      = '2024-12-31'
# ──────────────────────────────────────────────────
ASSETS               = ['BTC','ETH','SOL','BNB','AVAX','MATIC','LINK','DOT']
MAX_HORIZON          = 30   # test up to 30-day holding periods
TRANSACTION_COST_BPS = 20   # round-trip cost in basis points
TOP_N                = 3    # assets per leg
FACTORS              = {
    'MOM_7d':  {'lookback': 7,  'type': 'momentum'},
    'MOM_21d': {'lookback': 21, 'type': 'momentum'},
    'LOW_VOL': {'lookback': 21, 'type': 'low_vol'},
}
print('Config ready.')

rng = np.random.default_rng(42)

if USE_LIVE_DATA:
    import yfinance as yf
    prices = yf.download(TICKERS, start=START_DATE, end=END_DATE)['Close']
    prices.columns = ASSETS
    prices.dropna(inplace=True)
    returns = prices.pct_change().dropna()
    print(f'Live data: {len(returns)} days')
else:
    n    = 1200
    vols = np.array([0.65,0.75,1.20,0.70,1.10,1.40,1.15,0.90]) / np.sqrt(252)
    mu   = np.array([0.30,0.25,0.50,0.20,0.45,0.60,0.40,0.35]) / 252
    ar   = np.array([0.12,0.09,0.07,0.10,0.06,0.08,0.11,0.05])
    data = np.zeros((n, 8))
    for t in range(1, n):
        data[t] = ar * data[t-1] + rng.standard_normal(8) * vols + mu
    idx = pd.date_range('2020-01-01', periods=n, freq='B')
    returns = pd.DataFrame(data, columns=ASSETS, index=idx)
    print(f'Synthetic AR data: {n} days')

def build_signal(returns, lookback, factor_type):
    """
    Build factor signal from returns.

    Parameters
    ----------
    returns     : pd.DataFrame  Daily returns.
    lookback    : int           Signal lookback window.
    factor_type : str           'momentum' or 'low_vol'.

    Returns
    -------
    pd.DataFrame  Signal values (higher = prefer to go long).
    """
    if factor_type == 'momentum':
        return returns.rolling(lookback).mean()
    elif factor_type == 'low_vol':
        return -returns.rolling(lookback).std()  # negative vol → low vol → long
    return returns.rolling(lookback).mean()


def compute_factor_decay(returns, lookback, factor_type, max_h, top_n):
    """
    Compute the mean IC of a factor at each forward horizon from 1 to max_h.

    Parameters
    ----------
    returns     : pd.DataFrame  Daily returns.
    lookback    : int           Signal lookback window.
    factor_type : str           Factor type ('momentum' or 'low_vol').
    max_h       : int           Maximum forward horizon.
    top_n       : int           Assets per leg (not used directly but stored for metadata).

    Returns
    -------
    list  Mean IC at each horizon 1..max_h.
    """
    signal  = build_signal(returns, lookback, factor_type)
    ic_by_h = []
    for h in range(1, max_h + 1):
        fwd = returns.rolling(h).sum().shift(-h)
        daily_ics = []
        for t_idx in range(lookback, len(returns) - h):
            s = signal.iloc[t_idx]
            f = fwd.iloc[t_idx]
            valid = s.notna() & f.notna()
            if valid.sum() >= 4:
                ic, _ = stats.spearmanr(s[valid], f[valid])
                daily_ics.append(ic)
        ic_by_h.append(np.nanmean(daily_ics) if daily_ics else 0.0)
    return ic_by_h


decay_results = {}
for fname, fparams in FACTORS.items():
    print(f'Computing decay for {fname}...')
    decay_results[fname] = compute_factor_decay(
        returns, fparams['lookback'], fparams['type'], MAX_HORIZON, TOP_N)

decay_df = pd.DataFrame(decay_results, index=range(1, MAX_HORIZON + 1))
decay_df.index.name = 'horizon'
print('\n', decay_df.round(4).to_string())

tc = TRANSACTION_COST_BPS / 10000  # convert bps to decimal
horizons = np.arange(1, MAX_HORIZON + 1)

print(f'Transaction cost: {TRANSACTION_COST_BPS} bps = {tc:.4f}')
for fname in FACTORS:
    cumIC = np.cumsum(decay_results[fname])
    net   = cumIC - tc  # subtract fixed one-way trading cost
    opt_h = horizons[np.argmax(net)]
    print(f'{fname}: optimal hold = {opt_h} days (net IC at opt_h = {max(net):.4f})')

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
fig.suptitle('Factor Decay Curves', fontsize=13, fontweight='bold')

colors = ['#e53935','#1976d2','#43a047']
ax1 = axes[0]
for (fname, ics), color in zip(decay_results.items(), colors):
    ax1.plot(horizons, ics, lw=2, color=color, label=fname)
ax1.axhline(0, color='black', lw=0.8, ls='--')
ax1.set_xlabel('Forward Horizon (days)'); ax1.set_ylabel('IC (Spearman)')
ax1.legend(fontsize=9); ax1.set_title('Factor IC Decay by Horizon')

ax2 = axes[1]
for (fname, ics), color in zip(decay_results.items(), colors):
    cumIC = np.cumsum(ics)
    net   = cumIC - tc
    ax2.plot(horizons, net, lw=2, color=color, label=fname)
    opt_h = horizons[np.argmax(net)]
    ax2.axvline(opt_h, color=color, lw=0.8, ls=':', alpha=0.7)
ax2.axhline(0, color='black', lw=0.8, ls='--', label=f'Zero net (TC={TRANSACTION_COST_BPS}bps)')
ax2.set_xlabel('Forward Horizon (days)'); ax2.set_ylabel('Cumulative IC Net of TC')
ax2.legend(fontsize=8); ax2.set_title('Net Factor Return vs Holding Period')

plt.tight_layout(); plt.show()

decay_df.to_csv('factor_decay_curve.csv')
print('Saved: factor_decay_curve.csv')