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

# Set a random seed for reproducibility
np.random.seed(42)

# Generate mock daily returns data
# Let's simulate 252 daily returns (approx. one trading year)
num_days = 252
mean_return = 0.0005  # 0.05% daily average return
volatility = 0.01    # 1% daily standard deviation

# Generate normally distributed daily returns
daily_returns = np.random.normal(mean_return, volatility, num_days)

# Create a Pandas Series for easier handling
returns_series = pd.Series(daily_returns, name='Daily Returns')

print("--- Sample Daily Returns ---")
print(returns_series.head())


def calculate_variance(data):
    """
    Calculates the sample variance of a given dataset.

    Args:
        data (np.ndarray or pd.Series): The input data (e.g., asset returns).

    Returns:
        float: The sample variance.

    Formula:
        variance = sum((x_i - mean)^2) / (N - 1)
    """
    mean = np.mean(data)
    variance = np.sum((data - mean)**2) / (len(data) - 1)
    return variance


def calculate_std_dev(data):
    """
    Calculates the sample standard deviation of a given dataset.

    Args:
        data (np.ndarray or pd.Series): The input data (e.g., asset returns).

    Returns:
        float: The sample standard deviation.

    Formula:
        std_dev = sqrt(variance)
    """
    return np.sqrt(calculate_variance(data))

# Calculate and display variance and standard deviation
calculated_var = calculate_variance(returns_series)
calculated_std = calculate_std_dev(returns_series)

print(f"\nCalculated Variance: {calculated_var:.6f}")
print(f"Calculated Standard Deviation: {calculated_std:.6f}")

# Verify with NumPy's built-in functions (ddof=1 for sample standard deviation)
np_var = np.var(returns_series, ddof=1)
np_std = np.std(returns_series, ddof=1)

print(f"NumPy Variance (ddof=1): {np_var:.6f}")
print(f"NumPy Standard Deviation (ddof=1): {np_std:.6f}")

# Interpretation of Results:
# The variance and standard deviation calculated from our mock daily returns provide an estimate
# of the volatility of these returns. A daily standard deviation of approximately 0.01 (1%)
# suggests that, on an average day, the returns deviate from the mean return by about 1%.

from scipy.stats import norm

# Define the confidence level for VaR and CVaR
confidence_level = 0.99  # 99% confidence


def calculate_historical_var(returns, confidence_level=0.99):
    """
    Calculates the Historical Value at Risk (VaR).

    Args:
        returns (pd.Series or np.ndarray): Historical daily returns.
        confidence_level (float): The confidence level (e.g., 0.99 for 99%).

    Returns:
        float: The Historical VaR (as a positive loss).

    Explanation:
        Historical VaR is calculated by sorting the historical returns
        and finding the return at the specified percentile.
        For a 99% confidence level, we look for the 1st percentile (1 - 0.99).
    """
    if not isinstance(returns, pd.Series):
        returns = pd.Series(returns)

    # Sort the returns in ascending order
    sorted_returns = returns.sort_values(ascending=True)

    # Calculate the index corresponding to the percentile
    # For 99% confidence, we want the 1st percentile (1 - 0.99)
    var_index = int(len(sorted_returns) * (1 - confidence_level))

    # The VaR is the return at this index. We multiply by -1 to represent it as a positive loss.
    historical_var = -sorted_returns.iloc[var_index]

    return historical_var


def calculate_parametric_var(mean_return, std_dev, confidence_level=0.99):
    """
    Calculates the Parametric (Variance-Covariance) Value at Risk (VaR).
    Assumes returns are normally distributed.

    Args:
        mean_return (float): The mean of the returns.
        std_dev (float): The standard deviation of the returns.
        confidence_level (float): The confidence level (e.g., 0.99 for 99%).

    Returns:
        float: The Parametric VaR (as a positive loss).

    Formula:
        VaR = - (mean_return + Z-score * std_dev)
        Where Z-score is the inverse CDF of the normal distribution at (1 - confidence_level).
    """
    # Z-score for the desired confidence level (e.g., for 99% VaR, we need the Z-score for 1% tail)
    z_score = norm.ppf(1 - confidence_level)

    # Calculate VaR. We multiply by -1 to represent it as a positive loss.
    parametric_var = -(mean_return + z_score * std_dev)

    return parametric_var

# --- Calculate VaR using both methods ---

# Historical VaR
historical_var_value = calculate_historical_var(returns_series, confidence_level)
print(f"Historical VaR ({confidence_level*100:.0f}%): {historical_var_value:.4f}")

# Parametric VaR
mean_returns = np.mean(returns_series)
std_dev_returns = np.std(returns_series, ddof=1) # Use sample std dev

parametric_var_value = calculate_parametric_var(mean_returns, std_dev_returns, confidence_level)
print(f"Parametric VaR ({confidence_level*100:.0f}%): {parametric_var_value:.4f}")

# Interpretation of Results:
# A Historical VaR of, for example, 0.0210 (2.10%) at 99% confidence means that, based on past data,
# there's a 1% chance of losing more than 2.10% of the portfolio value in a single day.
# The Parametric VaR makes an assumption about the distribution (normal in this case).
# If the actual returns are not perfectly normal, the parametric VaR might differ from the historical VaR.
# Both values give us an estimate of the maximum expected loss at a given confidence level.

def calculate_cvar(returns, confidence_level=0.99):
    """
    Calculates the Conditional Value at Risk (CVaR) / Expected Shortfall.

    Args:
        returns (pd.Series or np.ndarray): Historical daily returns.
        confidence_level (float): The confidence level (e.g., 0.99 for 99%).

    Returns:
        float: The CVaR (as a positive loss).

    Explanation:
        CVaR is the average of the losses that are worse than the VaR.
        First, we calculate the VaR. Then, we filter the returns to only include
        those that are worse than the VaR (i.e., less than -VaR, since VaR is a positive loss).
        Finally, we take the average of these extreme losses.
    """
    if not isinstance(returns, pd.Series):
        returns = pd.Series(returns)

    # Calculate VaR at the given confidence level (as a negative return threshold)
    var_threshold = returns.quantile(1 - confidence_level)

    # Identify returns that are worse than the VaR threshold
    extreme_losses = returns[returns < var_threshold]

    if extreme_losses.empty:
        return 0.0 # No losses exceeded VaR

    # CVaR is the average of these extreme losses (multiplied by -1 for positive loss)
    cvar_value = -np.mean(extreme_losses)

    return cvar_value

# --- Calculate CVaR ---
cvar_value = calculate_cvar(returns_series, confidence_level)
print(f"CVaR ({confidence_level*100:.0f}%): {cvar_value:.4f}")

# Interpretation of Results:
# A CVaR of, for example, 0.0250 (2.50%) at 99% confidence means that if a loss does exceed
# the 99% VaR, the average loss experienced will be 2.50%. This value is typically higher
# than VaR, providing a more conservative estimate of tail risk.

plt.figure(figsize=(10, 6))
sns.histplot(returns_series, bins=30, kde=True, color='skyblue', stat='density')

# Convert VaR and CVaR to negative values for plotting on the returns axis
# Note: our calculate_historical_var and calculate_cvar return positive loss values
plot_historical_var = -historical_var_value
plot_cvar = -cvar_value

plt.axvline(x=plot_historical_var, color='red', linestyle='--', label=f'Historical VaR ({confidence_level*100:.0f}%)')
plt.axvline(x=plot_cvar, color='purple', linestyle=':', label=f'CVaR ({confidence_level*100:.0f}%)')

plt.title('Distribution of Daily Returns with VaR and CVaR', fontsize=14)
plt.xlabel('Daily Return', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()

# Interpretation of Visualization 1:
# The histogram illustrates the frequency of different daily return values. The peak around 0.0005
# confirms our simulated mean return. The red dashed line marks the Historical VaR threshold;
# returns to the left of this line represent the 1% (1 - confidence_level) worst-case scenarios.
# The purple dotted line indicates the CVaR, which is the average of all returns falling to the
# left of the VaR line. As expected, CVaR is further into the tail (more negative) than VaR,
# highlighting the average expected loss during extreme events.

# Calculate cumulative returns
cumulative_returns = (1 + returns_series).cumprod() - 1

# Assuming an initial portfolio value for visualization purposes
initial_portfolio_value = 1000
portfolio_value = initial_portfolio_value * (1 + cumulative_returns)

# Calculate daily VaR and CVaR as a percentage of current portfolio value
daily_var_abs = historical_var_value * portfolio_value.shift(1).fillna(initial_portfolio_value)
daily_cvar_abs = cvar_value * portfolio_value.shift(1).fillna(initial_portfolio_value)

plt.figure(figsize=(12, 7))

# Plot Portfolio Value
plt.plot(portfolio_value.index, portfolio_value, label='Portfolio Value', color='blue', linewidth=1.5)

# Plot VaR threshold (as a potential loss from previous day's close)
# We'll plot it as portfolio_value - daily_var_abs
plt.plot(portfolio_value.index, portfolio_value - daily_var_abs, label=f'VaR ({confidence_level*100:.0f}%) Lower Bound', color='red', linestyle='--', alpha=0.7)

# Plot CVaR threshold (as a potential loss from previous day's close)
plt.plot(portfolio_value.index, portfolio_value - daily_cvar_abs, label=f'CVaR ({confidence_level*100:.0f}%) Lower Bound', color='purple', linestyle=':', alpha=0.7)

plt.title('Simulated Portfolio Value with Daily VaR and CVaR Risk Bounds', fontsize=14)
plt.xlabel('Trading Day', fontsize=12)
plt.ylabel('Portfolio Value ($)', fontsize=12)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()

# Interpretation of Visualization 2:
# This chart shows how an initial investment would grow or shrink based on the simulated daily returns.
# The red dashed line represents the estimated lower bound for the portfolio value on a given day,
# as indicated by the 99% VaR. This means that, based on historical volatility, there's a 1% chance
# the portfolio could fall below this red line. The purple dotted line, representing CVaR, is a more
# severe lower bound, indicating the average level the portfolio could fall to if it breaches the VaR.
# This visualization helps in understanding the dynamic nature of risk bounds relative to portfolio performance.