part4.py

import numpy as np
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')


def compute_residuals(model, X_scaled, y_true):
    y_pred = model.predict(X_scaled)
    residuals = y_true - y_pred
    return residuals


def get_current_state_from_real_data(part1_results, lag_y=5, lag_sigma=5):
    prices = part1_results['prices']
    sigma = part1_results['sigma_extended']
    
    max_lag = max(lag_y, lag_sigma)
    
    state_prices = prices[-lag_y:]
    state_sigma = sigma[-lag_sigma:]
    
    return state_prices, state_sigma


def simulate_paths(model, scaler, part1_results, lag_y=5, lag_sigma=5,
                   steps_ahead=3, n_sim=10000, residuals=None, verbose=True):
    if verbose:
        print(f"\nSimulating {n_sim} paths for {steps_ahead} steps ahead...")
    
    state_prices_0, state_sigma_0 = get_current_state_from_real_data(
        part1_results, lag_y=lag_y, lag_sigma=lag_sigma
    )
    
    sims = np.zeros((n_sim, steps_ahead))
    
    for sim_idx in range(n_sim):
        prices_state = state_prices_0.copy()
        sigma_state = state_sigma_0.copy()
        
        for h in range(steps_ahead):
            features = np.concatenate([prices_state, sigma_state]).reshape(1, -1)
            features_scaled = scaler.transform(features)
            mu_hat = model.predict(features_scaled)[0]

            if residuals is not None and len(residuals) > 0:
                eps = np.random.choice(residuals)
                y_next = mu_hat + eps
            else:
                y_next = mu_hat
            
            sims[sim_idx, h] = y_next
            
            prices_state = np.concatenate([prices_state[1:], [y_next]])
            sigma_state = np.concatenate([sigma_state[1:], [sigma_state[-1]]])
        
        if verbose and (sim_idx + 1) % max(1, n_sim // 10) == 0:
            print(f"  Simulated {sim_idx + 1} / {n_sim} paths")
    
    if verbose:
        print("Simulation complete.")
    
    return sims


def fit_kde_for_horizons(sims, bandwidth='scott', verbose=True):
    n_sim, steps_ahead = sims.shape
    kde_funcs = []
    
    if verbose:
        print(f"Using bandwidth method: {bandwidth}")
    
    for h in range(steps_ahead):
        samples_h = sims[:, h]
        kde = gaussian_kde(samples_h, bw_method=bandwidth)
        kde_funcs.append(kde)
        
        if verbose:
            print(f"\nHorizon {h+1}:")
            print(f"  Mean: {samples_h.mean():.4f}")
            print(f"  Std:  {samples_h.std():.4f}")
            for q in [1, 5, 25, 50, 75, 95, 99]:
                print(f"  {q:2d}th pct: {np.percentile(samples_h, q):.4f}")
    
    if verbose:
        print("\nKDE estimation complete.")
    
    return kde_funcs


def plot_kde_for_horizons(kde_funcs, sims, last_price=None, n_points=500):
    n_sim, steps_ahead = sims.shape
    fig, axes = plt.subplots(steps_ahead, 1, figsize=(8, 3*steps_ahead))
    if steps_ahead == 1:
        axes = [axes]

    for h, ax in enumerate(axes):
        samples_h = sims[:, h]
        kde = kde_funcs[h]

        x_min, x_max = samples_h.min(), samples_h.max()
        padding = (x_max - x_min) * 0.2 if x_max > x_min else 0.1
        x_lo = x_min - padding
        x_hi = x_max + padding

        xs = np.linspace(x_lo, x_hi, n_points)
        ys = kde(xs)

        ax.hist(samples_h, bins=40, density=True, alpha=0.3, edgecolor='k', label='MC samples')
        ax.plot(xs, ys, color='blue', label=f"KDE (Horizon {h+1})")
        ax.fill_between(xs, ys, alpha=0.2, color='blue')

        mean_sim = samples_h.mean()
        ax.axvline(mean_sim, color='red', linestyle='--', label=f'Mean sim = {mean_sim:.3f}')
        if last_price is not None:
            ax.axvline(last_price, color='green', linestyle=':', label=f'Last price = {last_price:.3f}')

        ax.set_xlim(x_lo, x_hi)
        ax.set_title(f"KDE - Horizon {h+1}")
        ax.set_xlabel("Price")
        ax.set_ylabel("Density")
        ax.grid(alpha=0.2)
        ax.legend()

    plt.tight_layout()
    plt.show()


def main4(part1_results, part2_results, part3_results,
          lag_y=5, lag_sigma=5, steps_ahead=3, n_sim=10000,
          bandwidth='scott',
          verbose=True, visualize=True):
    model = part3_results['model']
    scaler = part3_results['scaler']
    X_train_scaled = part3_results['X_train']
    y_train_split = part3_results['y_train']
    
    last_price = part1_results['prices'][-1]
    
    if verbose:
        print("\nComputing residuals from training data...")
    residuals = compute_residuals(model, X_train_scaled, y_train_split)
    
    sims = simulate_paths(
        model, scaler, part1_results,
        lag_y=lag_y, lag_sigma=lag_sigma,
        steps_ahead=steps_ahead, n_sim=n_sim,
        residuals=residuals, verbose=verbose
    )
    
    kde_funcs = fit_kde_for_horizons(sims, bandwidth=bandwidth, verbose=verbose)
    
    terminal = sims[:, -1]
    ret = terminal / last_price - 1.0
    if verbose:
        print("\n[Diagnostics] Terminal price and return (Horizon):")
        print(f"  Last real price: {last_price:.4f}")
        print(f"  Mean terminal price: {terminal.mean():.4f}")
        print(f"  Mean return: {ret.mean():.4%}")
        for q in [1, 5, 25, 50, 75, 95, 99]:
            print(f"  {q:2d}th pct price:  {np.percentile(terminal, q):.4f}")
            print(f"  {q:2d}th pct return: {np.percentile(ret, q):.4%}")
    
    if visualize:
        plot_kde_for_horizons(kde_funcs, sims, last_price=last_price, n_points=500)
    
    results = {
        'sims': sims,
        'kde_funcs': kde_funcs,
        'residuals': residuals
    }
    
    return results


if __name__ == "__main__":
    np.random.seed(42)
    prices = np.linspace(100, 110, 50)
    sigma = np.linspace(0.01, 0.02, 50)
    part1_mock = {
        'prices': prices,
        'sigma_extended': sigma
    }
    
    from xgboost import XGBRegressor
    from sklearn.preprocessing import StandardScaler
    
    X = np.random.randn(1000, 10)
    y = X[:, 0] * 2 + np.random.randn(1000) * 0.1
    scaler = StandardScaler().fit(X)
    X_scaled = scaler.transform(X)
    model = XGBRegressor(n_estimators=10).fit(X_scaled, y)
    
    part3_mock = {
        'model': model,
        'scaler': scaler,
        'X_train': X_scaled,
        'y_train': y
    }
    
    part2_mock = {}
    
    results4 = main4(part1_mock, part2_mock, part3_mock,
                     lag_y=5, lag_sigma=5,
                     steps_ahead=3, n_sim=1000,
                     bandwidth='scott',
                     verbose=True, visualize=True)
    
    print("\nPart 4 test complete.")