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

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import sys
import warnings
warnings.filterwarnings("ignore")

#
np.random.seed(42)
sns.set_theme(style="whitegrid")
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['axes.titlesize'] = 14
plt.rcParams['axes.labelsize'] = 10

# Assets
ASSETS = {
    "SPY": "Equity (S&P 500)",
    "TLT": "Fixed Income (20+Y Treasury)",
    "SHY": "Cash-like (1-3Y Treasury)",
    "GLD": "Precious Metals (Gold)",
    "DBO": "Crude Oil"
}
TICKERS = list(ASSETS.keys())

# Dates
START_DATE = "2010-01-01"
END_DATE = "2022-12-30"

# Data Loading
print(f"Downloading data for {TICKERS}...")
raw_data = yf.download(TICKERS, start=START_DATE, end=END_DATE, auto_adjust=False, progress=False)

try:
    if 'Adj Close' in raw_data.columns:
        prices = raw_data['Adj Close']
    elif 'Close' in raw_data.columns:
        print("Note: Using 'Close' column.")
        prices = raw_data['Close']
    else:
        prices = raw_data.xs('Adj Close', level=0, axis=1)

    prices = prices.ffill().dropna()

    if prices.empty: raise ValueError("Data Empty.")
    print(f"Data Aligned & Loaded. Shape: {prices.shape}")

except Exception as e:
    sys.exit(f"CRITICAL DATA ERROR: {e}")

# Log Returns (Standard Differentiation, d=1)
log_returns = np.log(prices / prices.shift(1)).dropna()

# Squared Returns (Volatility)
squared_returns = log_returns ** 2

# Rolling Correlation
rolling_corr_spy_tlt = log_returns['SPY'].rolling(window=60).corr(log_returns['TLT']).dropna()

# Fractional Differentiation (Preserving Memory)
def get_weights_ffd(d, thres=1e-4):
    w, k = [1.], 1
    while True:
        w_k = -w[-1] / k * (d - k + 1)
        if abs(w_k) < thres: break
        w.append(w_k)
        k += 1
    return np.array(w[::-1]).reshape(-1, 1)

def frac_diff_ffd(series, d, thres=1e-4):
    w = get_weights_ffd(d, thres)
    width = len(w) - 1
    if width >= series.shape[0]: return pd.DataFrame(np.nan, index=series.index, columns=series.columns)
    df = {}
    series_val = series.fillna(method='ffill').dropna()
    for name in series_val.columns:
        series_f = series_val[[name]].fillna(method='ffill').dropna()
        df_ = pd.Series(0.0, index=series_f.index)
        for i in range(width, series_f.shape[0]):
            val = series_f.iloc[i-width:i+1].values
            if val.shape[0] == w.shape[0]:
                df_.iloc[i] = np.dot(w.T, val)[0,0]
            else:
                df_.iloc[i] = np.nan
        df[name] = df_
    return pd.DataFrame(df).iloc[width:]


frac_diffs = pd.DataFrame(index=prices.index)
d_values = {}
for ticker in TICKERS:
    best_d = 1.0
    for d in np.linspace(0.1, 1.0, 10):
        try:
            temp = frac_diff_ffd(prices[[ticker]], d=d, thres=1e-4).dropna()
            if len(temp) > 50 and adfuller(temp)[1] < 0.05:
                best_d = d
                break
        except: continue
    d_values[ticker] = best_d
    frac_diffs = frac_diffs.join(frac_diff_ffd(prices[[ticker]], d=best_d, thres=1e-4), how='outer')

#

# Price Dynamics
fig1, axes = plt.subplots(1, 2, figsize=(18, 7))

# Raw Prices
for col in prices.columns:
    axes[0].plot(prices.index, prices[col], label=col, linewidth=1.5)
axes[0].set_title("Figure 1: Raw Prices (Before Normalization)", fontweight='bold')
axes[0].set_ylabel("Price ($)")
axes[0].legend(loc="upper left")
axes[0].grid(True, alpha=0.3)

# Normalized Prices
norm_prices = prices / prices.iloc[0]
for col in norm_prices.columns:
    axes[1].plot(norm_prices.index, norm_prices[col], label=col, linewidth=1.5)
axes[1].set_title("Figure 1.1: Wealth Dynamics (Normalized to 1.0)", fontweight='bold')
axes[1].set_ylabel("Cumulative Return")
axes[1].axvline(pd.Timestamp('2018-01-01'), color='r', linestyle='--', label='Test Start')
axes[1].legend(loc="upper left")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Regime Changes & Correlation Matrix
fig2, axes = plt.subplots(1, 2, figsize=(18, 6))

# Dynamic Rolling Correlation
axes[0].plot(rolling_corr_spy_tlt.index, rolling_corr_spy_tlt, color='darkblue', linewidth=1)
axes[0].axhline(0, color='black', linestyle='--', alpha=0.5)
axes[0].fill_between(rolling_corr_spy_tlt.index, rolling_corr_spy_tlt, 0, where=(rolling_corr_spy_tlt > 0), color='red', alpha=0.1)
axes[0].fill_between(rolling_corr_spy_tlt.index, rolling_corr_spy_tlt, 0, where=(rolling_corr_spy_tlt < 0), color='green', alpha=0.1)
axes[0].set_title("Figure 1.2: Dynamic: 60-Day Rolling Corr (SPY vs TLT)", fontweight='bold')

# Static Correlation Matrix
corr_matrix = log_returns.corr()
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
sns.heatmap(corr_matrix, mask=mask, annot=True, cmap='coolwarm', center=0, square=True, ax=axes[1])
axes[1].set_title("Figure 1.3: Correlation Matrix", fontweight='bold')

plt.tight_layout()
plt.show()

# Fat Tail Analysis
fig3, axes = plt.subplots(2, 3, figsize=(18, 10))
axes_flat = axes.flatten()

for i, ticker in enumerate(TICKERS):
    ax = axes_flat[i]
    sns.histplot(log_returns[ticker], stat="density", kde=False, ax=ax, color='skyblue', alpha=0.6)
    mu, std = stats.norm.fit(log_returns[ticker])
    x = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 100)
    ax.plot(x, stats.norm.pdf(x, mu, std), 'r--', linewidth=2, label='Normal')

    kurt = log_returns[ticker].kurtosis()
    ax.set_title(f"{ticker} (Kurtosis={kurt:.2f})")
    ax.legend()
    if kurt > 3.0: ax.text(0.05, 0.9, "FAT TAIL", transform=ax.transAxes, color='red', fontweight='bold')

# Combined KDE
ax_final = axes_flat[5]
for ticker in TICKERS:
    sns.kdeplot(log_returns[ticker], ax=ax_final, fill=True, alpha=0.1, linewidth=1.5, label=ticker)
ax_final.set_title("Figure 1.3: Combined Density (KDE)")
ax_final.legend()
plt.tight_layout()
plt.show()

#
fig4, axes = plt.subplots(5, 2, figsize=(14, 15)) # Changed from 5,4 to 5,2
cols = ["ACF (Standard d=1)", "PACF (Standard d=1)"]

for i, col_name in enumerate(cols):
    axes[0, i].set_title(col_name, fontweight='bold', fontsize=14)

for i, ticker in enumerate(TICKERS):
    # Data
    std_data = log_returns[ticker]

    # Col 0: ACF Standard
    plot_acf(std_data, ax=axes[i, 0], lags=20, alpha=0.05, title="", auto_ylims=True)
    axes[i, 0].set_ylabel(ticker, fontweight='bold', fontsize=12)

    # Col 1: PACF Standard
    plot_pacf(std_data, ax=axes[i, 1], lags=20, alpha=0.05, title="", auto_ylims=True)

plt.suptitle("Serial Correlation: Standard Log Returns", y=1.01, fontsize=18, fontweight='bold')
plt.tight_layout()
plt.show()
#

# FIGURE 5: Stationarity Series Plot
fig5, axes = plt.subplots(2, 3, figsize=(18, 10))
axes_flat = axes.flatten()
for i, ticker in enumerate(TICKERS):
    ax = axes_flat[i]
    series = frac_diffs[ticker].dropna()
    ax.plot(series.index, series.values, color='purple', linewidth=0.8)
    ax.set_title(f"{ticker} Frac Diff Series (d={d_values[ticker]:.2f})", fontweight='bold')
axes_flat[5].axis('off')
plt.tight_layout()
plt.show()

#
adf_res = []
for ticker in TICKERS:
    res = adfuller(log_returns[ticker])
    adf_res.append({'Asset': ticker, 'P-Value': res[1], 'Result': "STATIONARY" if res[1]<0.05 else "FAIL"})
print(pd.DataFrame(adf_res))


stats_df = log_returns.describe().T
stats_df['Kurtosis'] = log_returns.kurtosis()
stats_df['Loss Function'] = np.where(stats_df['Kurtosis'] > 3.0, "Huber Loss", "MSE")
print(stats_df[['Kurtosis', 'Loss Function']])

print("\n--- 3. OPTIMAL WINDOW SIZE (Volatility Memory) ---")
window_suggestions = {}
for ticker in TICKERS:
    pacf_vals, confint = pacf(squared_returns[ticker], nlags=40, alpha=0.05)
    sig_lags = np.where(np.abs(pacf_vals) > 0.1)[0]
    sig_lags = sig_lags[sig_lags > 0]
    window_suggestions[ticker] = sig_lags[-1] if len(sig_lags) > 0 else 1
print(f"Recommended Window: {max(window_suggestions.values())} days")

Downloading data for ['SPY', 'TLT', 'SHY', 'GLD', 'DBO']...
Data Aligned & Loaded. Shape: (3271, 5)

  Asset       P-Value      Result
0   SPY  4.904129e-23  STATIONARY
1   TLT  8.364964e-19  STATIONARY
2   SHY  6.711575e-22  STATIONARY
3   GLD  0.000000e+00  STATIONARY
4   DBO  0.000000e+00  STATIONARY
         Kurtosis Loss Function
Ticker                         
DBO      5.988149    Huber Loss
GLD      5.155523    Huber Loss
SHY      8.747785    Huber Loss
SPY     11.555626    Huber Loss
TLT      4.147508    Huber Loss

--- 3. OPTIMAL WINDOW SIZE (Volatility Memory) ---
Recommended Window: 30 days

# ============================================================================
# STEP 2: LSTM
# ============================================================================

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math

# Ensure reproducibility
tf.random.set_seed(42)
np.random.seed(42)

# Global Plot Settings
plt.rcParams['figure.figsize'] = (16, 6)
plt.rcParams['font.family'] = 'sans-serif'


# Helper: RSI Calculation
def calculate_rsi(series, period=14):
    delta = series.diff()
    gain = (delta.where(delta > 0, 0)).rolling(window=period).mean()
    loss = (-delta.where(delta < 0, 0)).rolling(window=period).mean()
    rs = gain / loss
    return 100 - (100 / (1 + rs))

# Frac Diff Prices (Trend Preservation)
feat_prices = frac_diffs.copy()
feat_prices.columns = [f"{col}_frac" for col in feat_prices.columns]

#Volatility (Log-Scaled to normalize spikes)
feat_vol = np.log1p(log_returns.rolling(window=5).std().fillna(0) * 100)
feat_vol.columns = [f"{col}_vol" for col in feat_vol.columns]

# RSI (Momentum Indicator)
feat_rsi = pd.DataFrame()
for ticker in TICKERS:
    feat_rsi[f"{ticker}_rsi"] = calculate_rsi(prices[ticker], period=14)
feat_rsi = feat_rsi.fillna(50)

# Combine Features
features = pd.concat([feat_prices, feat_vol, feat_rsi], axis=1).dropna()

# TARGET: 25-Day Cumulative Log Return
horizon = 25
targets = pd.DataFrame()
for ticker in TICKERS:
    targets[ticker] = log_returns[ticker].rolling(window=horizon).sum().shift(-horizon)

targets = targets.dropna()

# Align Data
common_idx = features.index.intersection(targets.index)
features = features.loc[common_idx]
targets = targets.loc[common_idx]

# Split Train/Test (Cutoff: Jan 1, 2018)
test_start = pd.Timestamp('2018-01-01')
train_idx = features.index[features.index < test_start]
test_idx = features.index[features.index >= test_start]

# Validation Split (Last 20% of Train)
val_cut = int(0.8 * len(train_idx))
train_idx_final = train_idx[:val_cut]
val_idx = train_idx[val_cut:]

# ---------------------------
# 2. MODEL DEFINITION (WITH BATCH NORM)
# ---------------------------

def create_advanced_lstm(input_shape):
    model = keras.Sequential([
        # Layer 1: LSTM
        layers.LSTM(64, return_sequences=True, input_shape=input_shape),
        layers.BatchNormalization(), # Helps stabilize learning
        layers.Dropout(0.2),

        # Layer 2: LSTM
        layers.LSTM(32, return_sequences=False),
        layers.BatchNormalization(),
        layers.Dropout(0.2),

        # Dense Layers
        layers.Dense(16, activation='relu'),
        layers.Dense(1) # Linear output
    ])

    model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001), loss='mse')
    return model

def prepare_multivariate_sequences(data, target_series, window=60):
    X, y = [], []
    if len(data) <= window: return np.array(X), np.array(y)
    for i in range(window, len(data)):
        X.append(data[i-window:i])
        y.append(target_series[i])
    return np.array(X), np.array(y)

# ---------------------------
# 3. TRAINING
# ---------------------------

window_size = 60
models = {}
predictions = {}
metrics_list = []

print(f"\n{'='*80}")
print("STARTING TRAINING & VISUALIZATION LOOP")
print(f"{'='*80}\n")

for ticker in TICKERS:
    print(f"Processing Asset: {ticker}...")

    # --- A. Data Prep ---
    # Select specific features for this ticker
    ticker_cols = [c for c in features.columns if ticker in c]
    feature_data = features[ticker_cols].values
    target_data = targets[ticker].values

    # Scaling (Fit on Train only)
    scaler = MinMaxScaler(feature_range=(-1, 1))

    # Slicing
    X_train_raw = feature_data[features.index.isin(train_idx_final)]
    y_train_raw = target_data[features.index.isin(train_idx_final)]

    X_val_raw = feature_data[features.index.isin(val_idx)]
    y_val_raw = target_data[features.index.isin(val_idx)]

    X_test_raw = feature_data[features.index.isin(test_idx)]
    y_test_raw = target_data[features.index.isin(test_idx)]

    # Transform
    X_train_s = scaler.fit_transform(X_train_raw)
    X_val_s = scaler.transform(X_val_raw)
    X_test_s = scaler.transform(X_test_raw)

    # Sequence Creation
    train_X, train_y = prepare_multivariate_sequences(X_train_s, y_train_raw, window_size)
    val_X, val_y = prepare_multivariate_sequences(X_val_s, y_val_raw, window_size)
    test_X, test_y = prepare_multivariate_sequences(X_test_s, y_test_raw, window_size)

    # --- B. Training ---
    model = create_advanced_lstm((train_X.shape[1], train_X.shape[2]))

    history = model.fit(
        train_X, train_y,
        validation_data=(val_X, val_y),
        epochs=50,
        batch_size=32,
        verbose=0,
        callbacks=[
            keras.callbacks.EarlyStopping(monitor='val_loss', patience=8, restore_best_weights=True)
        ]
    )

    # Store Model & Predictions
    models[ticker] = model
    test_pred = model.predict(test_X, verbose=0).flatten()

    # --- ADDED: METRICS CALCULATION (For Questions 2b & 2c) ---
    train_pred = model.predict(train_X, verbose=0).flatten()
    val_pred = model.predict(val_X, verbose=0).flatten()

    # Calculate RMSE
    rmse_train = math.sqrt(mean_squared_error(train_y, train_pred))
    rmse_val = math.sqrt(mean_squared_error(val_y, val_pred))
    rmse_test = math.sqrt(mean_squared_error(y_test_raw[window_size:], test_pred))

    metrics_list.append({
        'Asset': ticker,
        'Train RMSE': rmse_train,
        'Val RMSE': rmse_val,
        'Test RMSE': rmse_test
    })
    # -----------------------------------------------------------

    pred_dates = features.index[features.index.isin(test_idx)][window_size:]

    # Save for Backtest Later
    predictions[ticker] = {
        'test_pred': test_pred,
        'test_dates': pred_dates
    }


    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 6))

    # Chart 1: Actual vs Predicted (Test Set)
    ax1.plot(pred_dates, y_test_raw[window_size:], label='Actual Return (25d)', color='grey', alpha=0.5, linewidth=1)
    ax1.plot(pred_dates, test_pred, label='LSTM Prediction', color='blue', linewidth=1.5)
    ax1.axhline(0, color='black', linestyle='--', linewidth=0.8)
    ax1.set_title(f"{ticker}: Out-of-Sample Performance (Test Set)", fontweight='bold')
    ax1.set_ylabel("Cumulative Log Return")
    ax1.legend(loc='upper left')
    ax1.grid(True, alpha=0.3)

    # Chart 2: Loss Function (Train vs Val)
    ax2.plot(history.history['loss'], label='Training Loss', color='navy')
    ax2.plot(history.history['val_loss'], label='Validation Loss', color='orange')
    ax2.set_title(f"{ticker}: Model Convergence (Loss Curve)", fontweight='bold')
    ax2.set_xlabel("Epochs")
    ax2.set_ylabel("MSE Loss")
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

# --- ADDED: PRINT METRICS REPORT (Satisfies 2b & 2c) ---
print("\n" + "-"*50)
print("MODEL PERFORMANCE REPORT (RMSE)")
print("="*50)
print(pd.DataFrame(metrics_list).round(5))
print("-" * 50)


# ---------------------------
# 4. TRADING STRATEGY (BACKTEST)
# ---------------------------
# (Keep the same logic as before for the portfolio construction)

print(f"\n{'='*80}")
print("CALCULATING STRATEGY RETURNS...")
print(f"{'='*80}\n")

# Reconstruct Predictions DataFrame
pred_df = pd.DataFrame()
for ticker in TICKERS:
    p_series = pd.Series(predictions[ticker]['test_pred'], index=predictions[ticker]['test_dates'])
    pred_df[ticker] = p_series
pred_df = pred_df.dropna()

portfolio_log_returns = []
portfolio_dates = []
rebalance_freq = 25

for i in range(0, len(pred_df), rebalance_freq):
    current_date = pred_df.index[i]
    if i + rebalance_freq >= len(pred_df): break

    # Rank
    ranks = pred_df.iloc[i].rank(ascending=False)
    longs = ranks[ranks <= 2].index.tolist()
    shorts = ranks[ranks > 3].index.tolist()

    # Real Returns
    start_loc = log_returns.index.get_loc(current_date)
    period_rets = log_returns.iloc[start_loc : start_loc+rebalance_freq]

    # Calculate Strategy Return
    l_ret = period_rets[longs].mean(axis=1) if longs else 0
    s_ret = period_rets[shorts].mean(axis=1) if shorts else 0
    strat_ret = 0.5 * l_ret - 0.5 * s_ret # Long - Short

    portfolio_log_returns.extend(strat_ret.values)
    portfolio_dates.extend(strat_ret.index)

# Final Comparison Plot
port_series = pd.Series(portfolio_log_returns, index=portfolio_dates)
avail_assets = pred_df.columns.tolist()
bench_ret = log_returns.loc[port_series.index, avail_assets].mean(axis=1)

# Cumulative (Simple)
cum_strat = np.exp(port_series.cumsum()) - 1
cum_bench = np.exp(bench_ret.cumsum()) - 1

plt.figure(figsize=(14, 7))
plt.plot(cum_strat.index, cum_strat, label='LSTM Long/Short Strategy', color='blue', linewidth=2)
plt.plot(cum_bench.index, cum_bench, label='Equal Weight Benchmark', color='grey', linestyle='--', linewidth=2)
plt.title("Final Backtest: Strategy vs Benchmark", fontweight='bold', fontsize=14)
plt.ylabel("Cumulative Return")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# --- ADDED: DETAILED METRICS FOR QUESTION 2e ---
total_strat = cum_strat.iloc[-1]
total_bench = cum_bench.iloc[-1]
ann_vol_strat = port_series.std() * np.sqrt(252)
sharpe_strat = (port_series.mean() * 252) / (port_series.std() * np.sqrt(252))

print("\n" + "-"*50)
print("FINAL STRATEGY vs BENCHMARK REPORT (Step 2e)")
print("-"*50)
print(f"Strategy Total Return:   {total_strat:.2%}")
print(f"Benchmark Total Return:  {total_bench:.2%}")
print(f"Strategy Annual Vol:     {ann_vol_strat:.2%}")
print(f"Strategy Sharpe Ratio:   {sharpe_strat:.2f}")
print("="*50)

================================================================================
STARTING TRAINING & VISUALIZATION LOOP
================================================================================

Processing Asset: SPY...

Processing Asset: TLT...

Processing Asset: SHY...

Processing Asset: GLD...

Processing Asset: DBO...

--------------------------------------------------
MODEL PERFORMANCE REPORT (RMSE)
==================================================
  Asset  Train RMSE  Val RMSE  Test RMSE
0   SPY     0.07878   0.03553    0.15012
1   TLT     0.08032   0.08631    0.12396
2   SHY     0.02135   0.01969    0.06880
3   GLD     0.06397   0.05418    0.04844
4   DBO     0.10213   0.06058    0.11714
--------------------------------------------------

================================================================================
CALCULATING STRATEGY RETURNS...
================================================================================

--------------------------------------------------
FINAL STRATEGY vs BENCHMARK REPORT (Step 2e)
--------------------------------------------------
Strategy Total Return:   -4.96%
Benchmark Total Return:  19.38%
Strategy Annual Vol:     9.16%
Strategy Sharpe Ratio:   -0.12
==================================================

# ============================================================================
# STEP 3: MULTI-OUTPUT LSTM
# ============================================================================

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math

# Ensure reproducibility
tf.random.set_seed(42)
np.random.seed(42)

# Global Plot Settings
plt.rcParams['figure.figsize'] = (16, 6)
plt.rcParams['font.family'] = 'sans-serif'
plt.style.use('seaborn-v0_8-whitegrid')

# ---------------------------
# 1. DATA PREPARATION (GLOBAL SCALING)
# ---------------------------
print("Preparing High-Performance Data Structures...")

# We use the 'features' and 'targets' defined in Step 2.
# CRITICAL FIX: Scale Features and Targets INDEPENDENTLY for better convergence.

# Split Indices
test_start = pd.Timestamp('2018-01-01')
train_idx = features.index[features.index < test_start]
test_idx = features.index[features.index >= test_start]

# Validation Split
val_cut = int(0.8 * len(train_idx))
train_idx_final = train_idx[:val_cut]
val_idx = train_idx[val_cut:]

# Global Scalers
scaler_X = MinMaxScaler((-1, 1))
scaler_y = MinMaxScaler((-1, 1))

# Fit Scalers on TRAIN set only
X_train_raw = features.loc[train_idx_final].values
y_train_raw = targets.loc[train_idx_final].values

scaler_X.fit(X_train_raw)
scaler_y.fit(y_train_raw)

# Transform All Data
X_scaled = scaler_X.transform(features.values)
y_scaled = scaler_y.transform(targets.values)

def prepare_multi_output_sequences(X_data, y_data, window=60):
    X, y = [], []
    if len(X_data) <= window: return np.array(X), np.array(y)
    for i in range(window, len(X_data)):
        X.append(X_data[i-window:i])
        y.append(y_data[i])
    return np.array(X), np.array(y)

# Create Sequences
window_size = 60

# Train
mask_train = features.index.isin(train_idx_final)
X_tr, y_tr = prepare_multi_output_sequences(X_scaled[mask_train], y_scaled[mask_train], window_size)

# Val
mask_val = features.index.isin(val_idx)
X_val, y_val = prepare_multi_output_sequences(X_scaled[mask_val], y_scaled[mask_val], window_size)

# Test
mask_test = features.index.isin(test_idx)
X_te, y_te = prepare_multi_output_sequences(X_scaled[mask_test], y_scaled[mask_test], window_size)

print(f"Input Shape: {X_tr.shape} (Samples, Window, Features)")
print(f"Output Shape: {y_tr.shape} (Samples, 5 Assets)")

# ---------------------------
# 2. ENCODER-DECODER ARCHITECTURE (The "Brain" Upgrade)
# ---------------------------
# This architecture forces the model to compress market info into a latent state
# before predicting, which captures cross-asset correlations better.

def create_encoder_decoder_model(input_shape, num_outputs=5):
    inputs = layers.Input(shape=input_shape)

    # --- ENCODER (Read the Market) ---
    # Process the sequence
    x = layers.LSTM(64, return_sequences=True)(inputs)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(0.3)(x)

    # Compress to Context Vector (The "Gist" of the market)
    x = layers.LSTM(32, return_sequences=False)(x)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(0.3)(x)

    # --- BOTTLENECK (Latent Space) ---
    x = layers.RepeatVector(num_outputs)(x) # Prepare for decoding per asset

    # --- DECODER (Predict the Future) ---
    x = layers.LSTM(32, return_sequences=True)(x)
    x = layers.BatchNormalization()(x)
    x = layers.Dropout(0.3)(x)

    # Flatten for final dense layer
    x = layers.Flatten()(x)

    # Interpretation
    x = layers.Dense(64, activation='relu')(x)

    # Output Layer (5 neurons, one for each asset)
    outputs = layers.Dense(num_outputs)(x)

    model = keras.Model(inputs=inputs, outputs=outputs)

    # We use a slightly lower learning rate for stability in multi-output
    model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.0005), loss='mse')
    return model

# ---------------------------
# 3. TRAINING & EVALUATION
# ---------------------------

print(f"\n{'='*80}")
print("TRAINING ENCODER-DECODER MULTI-OUTPUT MODEL")
print(f"{'='*80}\n")

model_mo = create_encoder_decoder_model((X_tr.shape[1], X_tr.shape[2]), num_outputs=5)

# Train with Early Stopping
history_mo = model_mo.fit(
    X_tr, y_tr,
    validation_data=(X_val, y_val),
    epochs=60, # More epochs because learning rate is lower
    batch_size=32,
    verbose=0,
    callbacks=[
        keras.callbacks.EarlyStopping(monitor='val_loss', patience=12, restore_best_weights=True)
    ]
)

# Generate Predictions (Scaled)
train_pred_scaled = model_mo.predict(X_tr, verbose=0)
val_pred_scaled = model_mo.predict(X_val, verbose=0)
test_pred_scaled = model_mo.predict(X_te, verbose=0)

# Inverse Transform (Back to Real Returns)
train_pred = scaler_y.inverse_transform(train_pred_scaled)
val_pred = scaler_y.inverse_transform(val_pred_scaled)
test_pred = scaler_y.inverse_transform(test_pred_scaled)

# Get Actuals (Inverse Scaled)
y_tr_real = scaler_y.inverse_transform(y_tr)
y_val_real = scaler_y.inverse_transform(y_val)
y_te_real = scaler_y.inverse_transform(y_te)

# ---------------------------
# 4. DASHBOARD & METRICS
# ---------------------------

mo_metrics = []
pred_dates = features.index[mask_test][window_size:]

# Store predictions for strategy
mo_predictions = pd.DataFrame(test_pred, index=pred_dates, columns=TICKERS)

for i, ticker in enumerate(TICKERS):
    # Metrics
    rmse_tr = math.sqrt(mean_squared_error(y_tr_real[:, i], train_pred[:, i]))
    rmse_val = math.sqrt(mean_squared_error(y_val_real[:, i], val_pred[:, i]))
    rmse_te = math.sqrt(mean_squared_error(y_te_real[:, i], test_pred[:, i]))
    mae_te = mean_absolute_error(y_te_real[:, i], test_pred[:, i])

    mo_metrics.append({
        'Asset': ticker,
        'Train RMSE': rmse_tr,
        'Val RMSE': rmse_val,
        'Test RMSE': rmse_te,
        'Test MAE': mae_te
    })

    # --- VISUALIZATION DASHBOARD ---
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 6))

    # Plot 1: Fit
    ax1.plot(pred_dates, y_te_real[:, i], label='Actual', color='grey', alpha=0.6)
    ax1.plot(pred_dates, test_pred[:, i], label='Multi-Output Prediction', color='purple', linewidth=1.5)
    ax1.axhline(0, color='black', linestyle='--', linewidth=0.8)
    ax1.set_title(f"{ticker}: Multi-Output Model Fit", fontweight='bold')
    ax1.set_ylabel("Cumulative Log Return")
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # Plot 2: Loss (Global)
    ax2.plot(history_mo.history['loss'], label='Global Train Loss', color='navy')
    ax2.plot(history_mo.history['val_loss'], label='Global Val Loss', color='orange')
    ax2.set_title(f"Global Model Convergence (Loss Curve)", fontweight='bold')
    ax2.set_xlabel("Epochs")
    ax2.set_ylabel("MSE Loss")
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

print("\n" + "="*50)
print("TABLE 3: MULTI-OUTPUT MODEL METRICS")
print("="*50)
print(pd.DataFrame(mo_metrics).round(5))
print("-" * 50)

# ---------------------------
# 5. STRATEGY & COMPARISON (3c & 3d)
# ---------------------------

print(f"\n{'='*80}")
print("MULTI-OUTPUT STRATEGY BACKTEST (3c & 3d)")
print(f"{'='*80}\n")

mo_portfolio_returns = []
mo_portfolio_dates = []
rebalance_freq = 25

# Reuse indices from pred_dates logic
aligned_preds = mo_predictions.dropna()

for i in range(0, len(aligned_preds), rebalance_freq):
    current_date = aligned_preds.index[i]
    if i + rebalance_freq >= len(aligned_preds): break

    # 1. Rank based on Multi-Output Predictions
    row = aligned_preds.iloc[i]
    ranks = row.rank(ascending=False)

    longs = ranks[ranks <= 2].index.tolist()
    shorts = ranks[ranks > 3].index.tolist()

    # 2. Real Returns
    start_loc = log_returns.index.get_loc(current_date)
    period_rets = log_returns.iloc[start_loc : start_loc+rebalance_freq]

    l_ret = period_rets[longs].mean(axis=1) if longs else 0
    s_ret = period_rets[shorts].mean(axis=1) if shorts else 0
    strat_ret = 0.5 * l_ret - 0.5 * s_ret

    mo_portfolio_returns.extend(strat_ret.values)
    mo_portfolio_dates.extend(strat_ret.index)

# Create Series
mo_series = pd.Series(mo_portfolio_returns, index=mo_portfolio_dates)
bench_series_mo = log_returns.loc[mo_series.index, TICKERS].mean(axis=1)

# Cumulative Returns
cum_mo_strat = np.exp(mo_series.cumsum()) - 1
cum_bench_mo = np.exp(bench_series_mo.cumsum()) - 1
# Note: 'cum_strat' from Step 2 should be available in memory.
# If running fresh, ensure Step 2 variable 'cum_strat' is saved or recalculated.
cum_step2_strat = cum_strat # Re-using variable from previous cell

# Comparison Plot
plt.figure(figsize=(16, 8))
plt.plot(cum_mo_strat.index, cum_mo_strat, label='Step 3: Multi-Output Strategy (Enc-Dec)', color='purple', linewidth=2.5)
plt.plot(cum_step2_strat.index, cum_step2_strat, label='Step 2: Single-Output Strategy', color='blue', linewidth=1.5, alpha=0.7)
plt.plot(cum_bench_mo.index, cum_bench_mo, label='Benchmark (Equal Weight)', color='grey', linestyle='--', linewidth=2)

plt.title("Step 3d: Model Comparison (Multi-Output vs Single-Output vs Benchmark)", fontweight='bold', fontsize=16)
plt.ylabel("Cumulative Return")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Final Metrics
total_mo = cum_mo_strat.iloc[-1]
total_s2 = cum_step2_strat.iloc[-1]
total_bn = cum_bench_mo.iloc[-1]
vol_mo = mo_series.std() * np.sqrt(252)
sharpe_mo = (mo_series.mean() * 252) / (mo_series.std() * np.sqrt(252))

print("\n" + "="*60)
print("TABLE 4: FINAL COMPARATIVE REPORT (Step 3d)")
print("="*60)
print(f"Multi-Output Strategy Return:  {total_mo:.2%}")
print(f"Single-Output Strategy Return: {total_s2:.2%}")
print(f"Benchmark Return:              {total_bn:.2%}")
print("-" * 60)
print(f"MO Strategy Volatility:        {vol_mo:.2%}")
print(f"MO Strategy Sharpe Ratio:      {sharpe_mo:.2f}")
print("="*60)

Preparing High-Performance Data Structures...
Input Shape: (1240, 60, 15) (Samples, Window, Features)
Output Shape: (1240, 5) (Samples, 5 Assets)

================================================================================
TRAINING ENCODER-DECODER MULTI-OUTPUT MODEL
================================================================================

==================================================
TABLE 3: MULTI-OUTPUT MODEL METRICS
==================================================
  Asset  Train RMSE  Val RMSE  Test RMSE  Test MAE
0   SPY     0.03148   0.02847    0.06085   0.04211
1   TLT     0.03818   0.02528    0.05296   0.04027
2   SHY     0.00205   0.00238    0.00547   0.00380
3   GLD     0.05357   0.03726    0.04692   0.03783
4   DBO     0.10478   0.11528    0.14396   0.11855
--------------------------------------------------

================================================================================
MULTI-OUTPUT STRATEGY BACKTEST (3c & 3d)
================================================================================

============================================================
TABLE 4: FINAL COMPARATIVE REPORT (Step 3d)
============================================================
Multi-Output Strategy Return:  -12.36%
Single-Output Strategy Return: 22.32%
Benchmark Return:              19.38%
------------------------------------------------------------
MO Strategy Volatility:        10.74%
MO Strategy Sharpe Ratio:      -0.27
============================================================

Introduction¶

Step 1¶

Results¶

Step 2¶

Results¶

Step 3¶

Results¶

Step 4¶

Conclusion¶

References¶