Risk Management Project

Market Risk Modelling and Risk Minimisation Strategy

Prepared by: Hai Nam Nguyen

This workbook develops a market risk model for a portfolio of financial assets, estimates risk measures (VaR and ES) under different assumptions, and implements a simple risk minimisation strategy. The project is structured into five stages: data preparation, risk measure estimation, Monte Carlo simulation, stressed VaR/ES estimation, and risk minimisation strategy development.

Technical stack: Python, NumPy, pandas, Matplotlib, seaborn, pycop, SciPy

1. Data and Environment Setup

Import libraries and daily closing prices of APPL and IBM from 31 December 1999 to 15 August 2025

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

from scipy import stats
from scipy import optimize
from scipy.stats import norm

df = pd.read_csv('AT2StockPriceData.csv', parse_dates= ["Date"], dayfirst= True)
df.index = pd.to_datetime(df["Date"], format="%d/%m/%Y")
df = df.set_index("Date").sort_index(ascending=True)
df = df.dropna()
display(df)
if df.isnull().values.any() == True: print("Data contains missing values.")

	AAPL	IBM
Date
1999-12-31	0.7716	53.6409
2000-01-03	0.8402	57.6784
2000-01-04	0.7694	55.7194
2000-01-05	0.7806	57.6784
2000-01-06	0.7130	56.6840
...	...	...
2025-08-11	227.1800	236.3000
2025-08-12	229.6500	234.7700
2025-08-13	233.3300	240.0700
2025-08-14	232.7800	237.1100
2025-08-15	231.5900	239.7200

6445 rows × 2 columns

2. Assumptions and notations

This section enumerates the assumptions that apply throughout the project and the notations used in the subsequent sections.

• The daily observations in the provided data are equally spaced. Assume that there are 250 trading days in a year.
• The risk management time horizon has length $\Delta t = 10/250$ corresponding to 10 trading days expressed in years.
• Neither stock pays dividends nor has corporate actions during the period of interest.
• The continuously compounded risk-free interest rate, denoted by $r$, is constant at 3.5% per annum for any maturity.
• Suppose the Black-Scholes-Merton model assumptions hold. Let $S_1(t)$ and $S_2(t)$ denote the current price of the APPL and IMB stocks, respectively. We assume that $\{S_i(t)\}_{t \ge 0}$ satisfy the stochastic differential equation:

$$ d S_i(t) = r S_i(t) dt + \sigma_i S_i(t) dW_i(t), \quad t \ge 0, \quad S_i(0) \ge 0 \quad i=1,2 $$

under the risk-neutral property measure $\mathbb{Q}$, where $W_1$ and $W_2$ are (possibly correlated) standard $\mathbb{Q}$-Brownian motions, and $\sigma_1$ and $\sigma_2$ are the constant volatilities of the two stocks.

• Let $C_{i,t}$ denote the current price of the European call option on the Apple stock ($i=1$) and on the IBM stock ($i=2$). Denote by $\Delta_{i,t}$, $\Gamma_{i,t}$, and $\Theta_{i,t}$ the current value of the option delta, gamma, and theta, respectively. Assume that these quantities can be computed using the Black-Scholes-Merton equations.

• Let $Z_i(t) = \ln S_i(t)$ denote the log-price of the $i$-th stock, $i=1,2$. The log prices will be used as the risk factors of the portfolio. Using the time series notation $Z_{i,t} = Z_i(t)$ denotes the current log-price of the $i$-th stock, $Z_{i,t+1} = Z_i(t+\Delta t)$ denotes the log-price of the stock at time $t+\Delta t$ (i.e. after 10 trading days) Let $X_i = Z_{i,t+1} - Z_{i,t}$ denote the log-return of the $i$-th stock over the risk management horizon for $i=1,2$. Let $\mathbf{X} = (X_{1,t+1}, X_{2,t+1})^\top$ denote the vector of risk factor changes.

3. Project Stage 1: Risk Factor Mapping

Construct a time series of overlapping 10-day log-returns, and determine the Black-Scholes-Merton (BMS) quantities and delta gamma loss approximations coefficients.

Mean, Covariance (EWMA) and Annualised Volatility of Risk Factor Changes

We first construct a time series of overlapping 10-day log-returns for each stock, and report the mean log-return for each stock. Then, using all available log-returns, we calculate and report the covariance matrix of the log-returns using the exponentially weighted moving average (EWMA) method (with smoothing constant 0.04 as of today), which is then used to compute the annualised volatilities which we denote by $\hat{\sigma}_{i}$ of the stocks.

np.random.seed(1234)

# Risk management time horizon length
delta_t = 10 / 250         # 10 trading days, assuming 250 trading days in a year

# Compute the 10-day stock price log returns or the 10-day risk factor changes
rfch = np.log(df[['AAPL','IBM']]) - np.log(df[['AAPL','IBM']].shift(10))         # Shift by 10 since date is increasing order
rfch = rfch.dropna()
rfch.columns = ['AAPL','IBM']

# Current value of each risk factor
aapl_price = df.iloc[-1]['AAPL'] 
ibm_price = df.iloc[-1]['IBM']    

# Mean of 10-day risk factor changes
rfch_mean = rfch.mean(axis = 0)
print("Mean of risk factor changes:")
print(f'Apple stock: {rfch_mean["AAPL"]:.4f}')
print(f'IBM stock: {rfch_mean["IBM"]:.4f}')

# Function for EWMA covariance matrix estimator of the returns time series
def EWMA_cov(returns_time_series, time, smoothing_constant):

    # - Inputs
    # returns_time_series: time series as a data frame
    # time: time (index) at which the EWMA covariance matrix estimate is to be calculated (must be an integer from 1 to
    #       no. of rows of returns_time_series - 1)
    # smoothing_constant (theta): EWMA smoothing parameter

    data = returns_time_series
    theta = smoothing_constant

    # - Return error message if time is incorrect
    if (time < 1 or time > (data.shape[0])):
        print("Invalid `time` value")

    # - Perform recursion to update EWMA covariance matrix up until specified time
    else:

        # - Extract the time series of specified length from the input time series
        data_extract = data.head(time)

        # - Normalize the returns and convert to a numpy array
        normalized = (data_extract - data_extract.mean()).fillna(0).to_numpy()

        # - Calculate the weights
        weights = (1 - theta) ** np.arange(len(data_extract))[::-1]

        # - Calculate EWMA estimate of covariance matrix
        cov = (weights * normalized.T) @ normalized * theta

        return cov
    
# EWMA covariance matrix of 10-day risk factor changes
rfch_cov_ewma = EWMA_cov(rfch, len(rfch), 0.04)
print(f'The EWMA Covariance matrix is:\n{np.round(rfch_cov_ewma, 4)}')

# Annual volatilities of the stocks based on 10-day log returns
aapl_ann_vol = np.sqrt(rfch_cov_ewma[0,0] / delta_t)
ibm_ann_vol = np.sqrt(rfch_cov_ewma[1,1] / delta_t)

print(f'Annualized Volatility of AAPL stock: {aapl_ann_vol:.4f}')
print(f'Annualized Volatility of IBM stock: {ibm_ann_vol:.4f}')

# Annual volatilities of the stocks based on 10-day log returns
aapl_ann_vol = np.sqrt(rfch_cov_ewma[0,0] / delta_t)
ibm_ann_vol = np.sqrt(rfch_cov_ewma[1,1] / delta_t)

print(f'Annualized Volatility of AAPL stock: {aapl_ann_vol:.4f}')
print(f'Annualized Volatility of IBM stock: {ibm_ann_vol:.4f}')

Mean of risk factor changes:
Apple stock: 0.0089
IBM stock: 0.0022
The EWMA Covariance matrix is:
[[ 0.0031 -0.0006]
 [-0.0006  0.0052]]
Annualized Volatility of AAPL stock: 0.2801
Annualized Volatility of IBM stock: 0.3607
Annualized Volatility of AAPL stock: 0.2801
Annualized Volatility of IBM stock: 0.3607

Black-Scholes-Merton option price, delta, gamma and theta Using the given information on the call options, the current stock prices, and the estimated annualized volatilities, we can calculate the (Black-Scholes-Merton) option price, delta, gamma, and theta. The continuously compounded risk-free interest rate, denoted by $r$, is constant at 3.5% per annum for any maturity. This means that Rho interest-rate sensitivity is irrelevant here because $r$ is fixed and not shocked.

# Pricing function for European calls and puts using Black-Scholes-Merton Formulas
def BSprice(spot_price, strike_price, time_to_maturity, risk_free_rate, dividend_yield, volatility, option_type):

    # Use mathematical notation for function inputs
    s = spot_price
    K = strike_price
    tau = time_to_maturity # (T-t) in notes
    r = risk_free_rate
    q = dividend_yield
    sigma = volatility

    # Calculate d1 and d2
    d1 = (np.log(s / K) + (r - q + 0.5 * sigma ** 2) * tau) / (sigma * np.sqrt(tau))
    d2 = d1 - sigma * np.sqrt(tau)

    # Calculate option prices
    if option_type == 'call':
        price = (np.exp(-q * tau) * s * stats.norm.cdf(d1, 0.0, 1.0) \
                 - np.exp(-r * tau) * K * stats.norm.cdf(d2, 0.0, 1.0))

    if option_type == 'put':
        price = (np.exp(-r * tau) * K * stats.norm.cdf(-d2, 0.0, 1.0) \
                 - np.exp(-q * tau) * s * stats.norm.cdf(-d1, 0.0, 1.0))

    return price

# Function to calculate option delta under the Black-Scholes-Merton model
def BSdelta(spot_price, strike_price, time_to_maturity, risk_free_rate, dividend_yield, volatility, option_type):

    # Use mathematical notation for function inputs
    s = spot_price
    K = strike_price
    tau = time_to_maturity # (T-t) in notes
    r = risk_free_rate
    q = dividend_yield
    sigma = volatility

    # Calculate d1 and d2
    d1 = (np.log(s / K) + (r - q + 0.5 * sigma ** 2) * tau) / (sigma * np.sqrt(tau))

    # Calculate option delta
    if option_type == 'call':
        value = stats.norm.cdf(d1, 0.0, 1.0)

    if option_type == 'put':
        value = stats.norm.cdf(d1, 0.0, 1.0) - 1

    return value

# Function to calculate option gamma under the Black-Scholes-Merton model
def BSgamma(spot_price, strike_price, time_to_maturity, risk_free_rate, dividend_yield, volatility):

    # Use mathematical notation for function inputs
    s = spot_price
    K = strike_price
    tau = time_to_maturity # (T-t) in notes
    r = risk_free_rate
    q = dividend_yield
    sigma = volatility

    # Compute option gamma (same for calls and puts)
    d1 = (np.log(s / K) + (r - q + 0.5 * sigma ** 2) * tau) / (sigma * np.sqrt(tau))
    value = stats.norm.pdf(d1, 0.0, 1.0) / (s * sigma * np.sqrt(tau))

    return value

# Function to calculate option theta (time decay) under the Black-Scholes-Merton model
def BStheta_call(spot_price, strike_price, time_to_maturity, risk_free_rate, dividend_yield, volatility):

    # Use mathematical notation for function inputs
    s = spot_price
    K = strike_price
    tau = time_to_maturity # (T-t) in notes
    r = risk_free_rate
    q = dividend_yield
    sigma = volatility

    # Compute call option theta
    d1 = (np.log(s / K) + (r - q + 0.5 * sigma ** 2) * tau) / (sigma * np.sqrt(tau))
    d2 = d1 - sigma * np.sqrt(tau)
    value = -(s * norm.pdf(d1, 0.0, 1.0) * sigma * np.exp(-q * tau)) / (2 * np.sqrt(tau)) - r * K * np.exp(-r * tau) * norm.cdf(d2, 0.0, 1.0) + q * s * norm.cdf(d1, 0.0, 1.0) * np.exp(-q * tau)

    return value

# Option prices using standard BSM formulas
aapl_call_price = BSprice(aapl_price, aapl_price, 3/12, 0.035, 0, aapl_ann_vol, 'call')
ibm_call_price = BSprice(ibm_price, ibm_price, 6/12, 0.035, 0, ibm_ann_vol, 'call')
print(f"The at-the-money call option price on AAPL stock is: $ {aapl_call_price:.4f}")
print(f"The at-the-money call option price on IBM stock is: $ {ibm_call_price:.4f}")

# Option delta
aapl_call_delta = BSdelta(aapl_price, aapl_price, 3/12 , 0.035, 0, aapl_ann_vol, "call")
ibm_call_delta = BSdelta(ibm_price, ibm_price, 6/12 , 0.035, 0, ibm_ann_vol, "call")
print(f"The option delta of the at-the-money call option on AAPL stock is: {aapl_call_delta:.4f}")
print(f"The option delta of the at-the-money call option on IBM stock is: {ibm_call_delta:.4f}")

# Option gamma
aapl_call_gamma = BSgamma(aapl_price, aapl_price, 3/12 , 0.035, 0, aapl_ann_vol)
ibm_call_gamma = BSgamma(ibm_price, ibm_price, 6/12 , 0.035, 0, ibm_ann_vol)
print(f"The option gamma of the at-the-money call option on AAPL stock is: {aapl_call_gamma:.4f}")
print(f"The option gamma of the at-the-money call option on IBM stock is: {ibm_call_gamma:.4f}")

# Call option theta of Apple and IBM
aapl_call_theta = BStheta_call(aapl_price, aapl_price, 3/12, 0.035, 0 , aapl_ann_vol)
ibm_call_theta = BStheta_call(ibm_price, ibm_price, 6/12, 0.035, 0 , ibm_ann_vol)
print(f"The option theta of the call option of AAPL stock is: {aapl_call_theta:.4f}")
print(f"The option theta of the call option of IBM stock is: {ibm_call_theta:.4f}")
print("Long calls usually have negative theta as they lose value as time passes")

The at-the-money call option price on AAPL stock is: $ 13.9077
The at-the-money call option price on IBM stock is: $ 26.2520
The option delta of the at-the-money call option on AAPL stock is: 0.5527
The option delta of the at-the-money call option on IBM stock is: 0.5778
The option gamma of the at-the-money call option on AAPL stock is: 0.0122
The option gamma of the at-the-money call option on IBM stock is: 0.0064
The option theta of the call option of AAPL stock is: -29.6484
The option theta of the call option of IBM stock is: -27.8571
Long calls usually have negative theta as they lose value as time passes

Delta-Gamma Approximation for the portfolio

• Portfolio Value:

We have the function of the portfolio value as follows, which includes going long on Apple and IBM call options and shorting Apple and IBM stocks:

$$ V(t, Z_{1,t}, Z_{2,t}) = C_1(t,Z_{1,t}) + C_2(t,Z_{2,t}) + \eta_1 e^{Z_{1,t}} + \eta_2 e^{Z_{2,t}} $$

Given that $\eta_1 = \eta_2 = 1$, the function then becomes:

$$ V(t, Z_{1,t}, Z_{2,t}) = C_1(t,Z_{1,t}) + C_2(t,Z_{2,t}) + e^{Z_{1,t}} + e^{Z_{2,t}} $$

Taylor Expansion of Portfolio Change $$ V(t+\Delta t, Z_{1,t}+X_{1,t+1}, Z_{2,t}+X_{2,t+1}) - V(t,Z_{1,t},Z_{2,t}) \approx \frac{\partial V}{\partial t}\Delta t + \frac{\partial V}{\partial Z_1} X_1 + \frac{\partial V}{\partial Z_2} X_2 + \tfrac{1}{2}\Bigg( \frac{\partial^2 V}{\partial Z_1^2} X_1^2 + \frac{\partial^2 V}{\partial Z_2^2} X_2^2 + 2\frac{\partial^2 V}{\partial Z_1 \partial Z_2} X_1X_2 \Bigg) $$

• Theta or Time sensitivity:

$$ \frac{\partial V}{\partial t} = \frac{\partial{}}{\partial t} (C_1(t,Z_{1,t})+C_2(t,Z_{2,t})) = \Theta_{1,\text{call}} + \Theta_{2,\text{call}} $$

• Delta sensitivities:

$$ \frac{\partial V}{\partial {Z_{1,t}}} = \frac{\partial{}}{\partial Z_{1,t}} (C_1(t,Z_{1,t}) - e^{Z_{1,t}}) = \frac{\partial C_1}{\partial S_1}\frac{\partial S_1}{\partial Z_1} - e^{Z_{1,t}} = e^{Z_{1,t}}\Delta_{1,\text{call}} - e^{Z_{1,t}} = e^{Z_{1,t}}(\Delta_{1,\text{call}}-1) $$

Similarly,

$$ \frac{\partial V}{\partial {Z_{2,t}}} = \frac{\partial C_2}{\partial S_2}\frac{\partial S_2}{\partial Z_2} - e^{Z_{2,t}} = e^{Z_{2,t}}(\Delta_{2,\text{call}}-1) $$

• Gamma sensitivities:

$$ \frac{\partial^2 V}{\partial {Z^2_{1,t}}} = \frac{\partial{}}{\partial {Z_{1,t}}} (e^{Z_{1,t}} \Delta_{1,call}) - \frac{\partial {e^{Z_{1,t}}}}{\partial{Z_{1,t}}} = e^{Z_{1,t}} \frac{\partial {\Delta_{1,call}}}{\partial{Z_{1,t}}} + \Delta_{1,call} \frac{\partial {e^{Z_{1,t}}}}{\partial{Z_{1,t}}} - e^{Z_{1,t}} = \Gamma_{1,\text{call}} e^{2Z_{1,t}} + \Delta_{1,\text{call}} e^{Z_{1,t}} - e^{Z_{1,t}}$$

$$\quad$$

With $e^{Z_{1,t}} \frac{\partial {\Delta_{1,call}}}{\partial{Z_{1,t}}} = \frac{\partial {\Delta_{1,call}}}{\partial{S_1}} \frac{\partial{S_1}} {\partial {Z_{1,t}}} = \Gamma_{1,\text{call}} e^{Z_{1,t}} e^{Z_{1,t}}$

Similarly,

$$ \frac{\partial^2 V}{\partial {Z^2_{2,t}}} = \Gamma_{2,\text{call}} e^{2Z_{2,t}} + \Delta_{2,\text{call}} e^{Z_{2,t}} - e^{Z_{2,t}} $$

The delta approximation for the portfolio loss over $\Delta t$ is: $$ \boxed{ L_{t+1}^{\Delta} = -(\Theta_{1,\text{call}}+\Theta_{2,\text{call}}) \Delta t - \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1,\text{call}} -1) \\ e^{Z_{2,t}}(\Delta_{2,\text{call}} -1) \end{bmatrix}^\top \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} } $$

The delta-gamma approximation for the portfolio loss over $\Delta t$ is:

$$ L_{t+1}^{\Delta,\Gamma} = L_{t+1}^{\Delta} - \tfrac{1}{2} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix}^\top \begin{bmatrix} \frac{\partial^2 V}{\partial Z_1^2} & 0 \\ 0 & \frac{\partial^2 V}{\partial Z_2^2} \end{bmatrix} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} $$

$$ \boxed{ L_{t+1}^{\Delta,\Gamma} = -(\Theta_{1,\text{call}}+\Theta_{2,\text{call}}) \Delta t - \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1,\text{call}} -1) \\ e^{Z_{2,t}}(\Delta_{2,\text{call}} -1) \end{bmatrix}^\top \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} - \tfrac{1}{2} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix}^\top \begin{bmatrix} \frac{\partial^2 V}{\partial Z_1^2} & 0 \\ 0 & \frac{\partial^2 V}{\partial Z_2^2} \end{bmatrix} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} } $$

Or we can express in the form of:

$$ \boxed{ L^{\Delta, \Gamma}_{t+1} = - \mathbf{c}_t - \mathbf{b}^\top_t X_{t+1} - \tfrac{1}{2} X^\top_{t+1} \mathbf{B}_t X_{t+1} }$$

Where, $$ \boxed{ {\mathbf{c}}_t = (\Theta_{1,\text{call}}+\Theta_{2,\text{call}}) \Delta t }$$

$$ \boxed{ {\mathbf{b}}_t = \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1,\text{call}} -1) \\ e^{Z_{2,t}}(\Delta_{2,\text{call}} -1) \end{bmatrix} } $$

$$ \boxed{ \mathbf{B}_t = \begin{bmatrix} \Gamma_{1,\text{call}} e^{2Z_{1,t}} + \Delta_{1,\text{call}} e^{Z_{1,t}} - e^{Z_{1,t}} & 0 \\ 0 & \Gamma_{2,\text{call}} e^{2Z_{2,t}} + \Delta_{2,\text{call}} e^{Z_{2,t}} - e^{Z_{2,t}} \end{bmatrix} } $$

# Quantities Apple stock and IBM stock
eta1 = 1
eta2 = 1

# Time decay (theta) term of the call options
time_decay = (aapl_call_theta + ibm_call_theta) * delta_t           # per 10-trading day time horizon
print(f"Time decay (theta) term c_t of the call options: {time_decay:.4f}")

# Call option exposure of AAPL stock
aapl_call_exposure = np.array([aapl_price * aapl_call_delta - eta1 * aapl_price, 0])

# Call option exposure of ibm stock
ibm_call_exposure = np.array([0, ibm_price * ibm_call_delta - eta2 * ibm_price])

# Portfolio exposure b_t
portfolio_exposure = aapl_call_exposure + ibm_call_exposure
print(f"Portfolio exposure b_t: {np.round(portfolio_exposure, 4)}")

# APPL call option exposure (second-order)
aapl_call_exposure_SO = np.zeros([rfch.shape[1], rfch.shape[1]])
aapl_call_exposure_SO[0,0] = aapl_call_gamma * aapl_price **2  + aapl_call_delta * aapl_price - eta1*aapl_price

# ibm call option exposure (second-order)
ibm_call_exposure_SO = np.zeros([rfch.shape[1], rfch.shape[1]])
ibm_call_exposure_SO[1,1] = ibm_call_gamma * ibm_price **2  + ibm_call_delta * ibm_price - eta2*ibm_price

# Portfolio exposure (second-order)
portfolio_exposure_SO = aapl_call_exposure_SO + ibm_call_exposure_SO
print(f"Portfolio exposure (second-order) B_t is:")
print(f"\n{np.round(portfolio_exposure_SO, 4)}")

Time decay (theta) term c_t of the call options: -2.3002
Portfolio exposure b_t: [-103.5886 -101.2216]
Portfolio exposure (second-order) B_t is:

[[550.2792   0.    ]
 [  0.     266.5728]]

3. Project Stage 2: Market Risk Modelling

Estimate the 10-day 99% Value-at-Risk (VaR) and Expected Shortfall (ES) of the portfolio using different methods and model assumptions.

Analytical Approach

We assume a linear loss operator $L_{t+1} = -c_t - \mathbf{b}_t^\top \mathbf{X}_{t+1}$, where $\mathbf{X}_{t+1} \sim N_d(\boldsymbol{\mu}_{t+1}, \boldsymbol{\Sigma}_{t+1})$ and $c_t \in \mathbb{R}$ and $\mathbf{b}_t \in \mathbb{R}^d$ are constants dependent only on information available at time $t$.

At 99% confidence level, where $\alpha \in (0,1)$, the VaR and ES are calculated as follows: $${\text{VaR}}_{\alpha}(L_{t+1}) = \mu_L + \sigma_L \Phi^{-1}(\alpha), \quad{\text{ES}}_{\alpha}(L_{t+1}) = \mu_L + \sigma_L \frac {{\phi}({\Phi}^{-1}({\alpha}))} {1-{\alpha}}$$ where ${\Phi}$ and ${\phi}$ are the cdf and pdf of standard normal distribution $N(0,1)$.

# Functions for VaR and ES under the variance-covariance approach
def VaR_ES_vcv(constant_loss, risk_exposure, risk_factor_mean, risk_factor_cov_mat, conf_level):

    # Shorthand notation for inputs
    c = constant_loss
    b = risk_exposure
    mu = risk_factor_mean
    Sigma = risk_factor_cov_mat
    alpha = conf_level

    # Compute mean and variance of loss
    loss_mean = -c - np.transpose(b) @ mu
    loss_var = np.transpose(b) @ Sigma @ b

    # Compute VaR
    VaR = loss_mean + np.sqrt(loss_var) * stats.norm.ppf(alpha, loc = 0.0, scale = 1.0)
    VaR = max(VaR, 0)

    # Compute ES
    ES = loss_mean + np.sqrt(loss_var) * stats.norm.pdf(stats.norm.ppf(alpha, loc = 0.0, scale = 1.0), loc = 0.0, scale = 1.0) / (1 - alpha)
    ES = max(ES, 0)

    # Return output
    return np.array([VaR, ES])

# Confidence level
conf_level = 0.99

# Undiversified VaR and ES
# Apple call option risk measures
aapl_call_risk_measures = VaR_ES_vcv(aapl_call_theta * delta_t, aapl_call_exposure, rfch_mean, rfch_cov_ewma, conf_level)
print(f'VaR and ES for Apple call option: \n{np.round(aapl_call_risk_measures, 4)}')

# IBM call option risk measures
ibm_call_risk_measures = VaR_ES_vcv(ibm_call_theta * delta_t, ibm_call_exposure, rfch_mean, rfch_cov_ewma, conf_level)
print(f'VaR and ES for IBM call option: \n{np.round(ibm_call_risk_measures, 4)}')

# Undiversified portfolio risk measures
portfolio_risk_measures = aapl_call_risk_measures + ibm_call_risk_measures
print(f'VaR and ES for portfolio (undiversified): \n{np.round(portfolio_risk_measures, 4)}')

# Diversified portfolio risk measures
portfolio_risk_measures_d = VaR_ES_vcv(time_decay, portfolio_exposure, rfch_mean, rfch_cov_ewma, conf_level)
print(f'VaR and ES for portfolio (diversified):\n{np.round(portfolio_risk_measures_d, 4)}')

VaR and ES for Apple call option: 
[15.605  17.5716]
VaR and ES for IBM call option: 
[18.3294 20.804 ]
VaR and ES for portfolio (undiversified): 
[33.9344 38.3756]
VaR and ES for portfolio (diversified):
[23.6119 26.5496]

# Compare undiversified and diversified VaR and ES
print(f'Difference in VaR: {(portfolio_risk_measures[0] - portfolio_risk_measures_d[0]):.4f}')
print(f'Difference in ES: {(portfolio_risk_measures[1] - portfolio_risk_measures_d[1]):.4f}')

Difference in VaR: 10.3225
Difference in ES: 11.8261

We could cleary that by having a diversified portfolio of Apple and IBM call options, the risk measures $\text{VaR}$ and $\text{ES}$ are significantly lower than without diversification.

Historical Simulation Approach

Delta-Gamma Loss Approximation: Using the standard historical simulation approach and the delta-gamma loss approximation, we can calculate the 10-day 99% VaR and ES of the portfolio. From these calculations, can then calculate the losses of the portfolio and plot the histogram of the losses, which shows a right-skewed distribution with a long tail, indicating that there are some extreme losses that can occur with low probability.

The delta-gamma loss approximation is given by

$$\begin{align*} L_{\text{call}, t+1}^{\Delta, \Gamma} & = L_{\text{call}, t+1}^\Delta - \frac{1}{2} \mathbf{X}_{t+1}^\top \mathbf{B}_{\text{call}, t} \mathbf{X}_{t+1} \\ \end{align*}$$

$$ L_{t+1}^{\Delta,\Gamma} = -(\Theta_{1,\text{call}}+\Theta_{2,\text{call}}) \Delta t - \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1,\text{call}} -1) \\ e^{Z_{2,t}}(\Delta_{2,\text{call}} -1) \end{bmatrix}^\top \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} - \tfrac{1}{2} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix}^\top \begin{bmatrix} \frac{\partial^2 V}{\partial Z_1^2} & 0 \\ 0 & \frac{\partial^2 V}{\partial Z_2^2} \end{bmatrix} \begin{bmatrix} X_{1,t+1} \\ X_{2,t+1} \end{bmatrix} $$

# Define VaR and ES from empirical loss distribution
def VaR_ES_emp(loss_scenarios, conf_level, hs_weight_type, age_weight_decay = 0.94):
    
    # Inputs:
    # loss_scenarios - array of historically simulated loss amounts, arranged from least recent to most recent
    # conf_level - risk measure confidence level
    # hs_sim_type - "equal" or "age"
    #     "equal": traditional historical simulation, loss scenarios equally weighted and quantile is linearly interpolated
    #     "BRW": BRW approach with exponentially decaying weights
    #          : Requires "age_weight_decay" parameter
    
    if hs_weight_type == "equal":
        
        # Determine the VaR (empirical quantile, linearly interpolated)
        VaR = np.quantile(loss_scenarios, conf_level, method = "higher")     
        
        # Determine the ES
        ES = np.maximum(loss_scenarios - VaR, 0).mean() / (1 - conf_level) + VaR

    elif hs_weight_type == "BRW":

        # Calculate exponentially decaying weights
        weights = age_weight_decay ** np.arange(len(loss_scenarios))[::-1] * (1 - age_weight_decay) / (1 - age_weight_decay ** len(loss_scenarios))

        # Append weights to historical simulations and sort in increasing order of the loss
        loss_weighted = np.column_stack((loss_scenarios, weights))
        loss_weighted = loss_weighted[loss_weighted[:,0].argsort()]

        # Calculate cumulative weights and append to loss_weighted
        cumulative_weight = np.cumsum(loss_weighted[:,1])
        loss_weighted = np.column_stack((loss_weighted, cumulative_weight))
        
        # Determine the VaR
        loss_exceed = loss_weighted[(loss_weighted[:,2] >= conf_level)]
        VaR = min(loss_exceed[:,0])

        # Calculate the ES
        ES = sum((loss_exceed[:,0] - VaR) * loss_exceed[:,1]) / (1 - conf_level) + VaR

    # Return output
    return np.array([VaR, ES])


# Construct loss scenarios
loss_scenarios = - time_decay - np.matmul(rfch.to_numpy(), portfolio_exposure) - 0.5 * np.diag(np.matmul(rfch.to_numpy(), np.matmul(portfolio_exposure_SO, rfch.to_numpy().T)))

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(loss_scenarios, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

Full Valuation Approach: Using the full valuation approach, we can calculate the 10-day 99% VaR and ES of the portfolio. The full valuation approach involves revaluing the portfolio at each historical scenario, which captures the non-linearities in the portfolio's value changes more accurately than the delta-gamma approximation. The resulting loss distribution from the full valuation approach is typically more accurate and may show a different shape compared to the delta-gamma approximation, especially in cases where there are significant non-linearities in the portfolio's value changes.

# Risk factor scenarios
aapl_price_scenarios = aapl_price * np.exp(rfch['AAPL'])
ibm_price_scenarios = ibm_price * np.exp(rfch['IBM'])

# Call option price scenarios
aapl_call_scenarios = BSprice(aapl_price_scenarios, aapl_price, 3/12 - delta_t, 0.035 , 0, aapl_ann_vol, 'call')
ibm_call_scenarios = BSprice(ibm_price_scenarios, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
portfolio_value = aapl_call_price + ibm_call_price - aapl_price - ibm_price
portfolio_scenarios = aapl_call_scenarios + ibm_call_scenarios - aapl_price_scenarios - ibm_price_scenarios
loss_scenarios_fv = -(portfolio_scenarios - portfolio_value)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(loss_scenarios_fv, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

Standard Historical Simulation using Delta-Gamma Approximation method: Another attempt is to use the standard historical simulation approach with the delta-gamma approximation method, which involves calculating the losses using the delta-gamma approximation for each historical scenario. This method is computationally less intensive than the full valuation approach while still capturing some of the non-linearities in the portfolio's value changes. The resulting loss distribution may be similar to that obtained from the full valuation approach, but it may not capture all the nuances of the portfolio's risk profile as accurately as the full valuation method.

rm_hs_dg = VaR_ES_emp(loss_scenarios, conf_level = conf_level, hs_weight_type = "equal")  # loss scenarios sorted from least recent to most recent

print(f'Portfolio VaR: {rm_hs_dg[0]:.4f}')
print(f'Portfolio ES: {rm_hs_dg[1]:.4f}')

Portfolio VaR: 20.0579
Portfolio ES: 23.0180

Standard Historical Simulation using the Full Valuation method: This method involves calculating the losses using the full valuation approach for each historical scenario, which captures all the non-linearities in the portfolio's value changes. This method is computationally intensive but provides the most accurate representation of the portfolio's risk profile. The resulting loss distribution may show significant differences compared to the delta-gamma approximation, especially in cases where there are large non-linearities in the portfolio's value changes.

rm_hs_fv = VaR_ES_emp(loss_scenarios_fv, conf_level = conf_level, hs_weight_type = "equal")

print(f'Portfolio VaR: {rm_hs_fv[0]:.4f}')
print(f'Portfolio ES: {rm_hs_fv[1]:.4f}')

Portfolio VaR: 19.9890
Portfolio ES: 23.1564

Comparison of VaR and ES Estimates from Historical Simulation and Delta-Gamma Approximation

# Comparison of the results of the standard historical simulations
if rm_hs_dg[0] - rm_hs_fv[0] > 0:
    print(f"The Delta-Gamma Approximation method gives a higher VaR of {(rm_hs_dg[0] - rm_hs_fv[0]):.4f} than the Full Valuation Approach")
else:
    print(f"The Full Valuation Approach gives a higher VaR of {(rm_hs_fv[0] - rm_hs_dg[0]):.4f} than the Delta-Gamma Approximation method")

if rm_hs_dg[1] - rm_hs_fv[1] > 0:
    print(f"The Delta-Gamma Approximation method gives a higher ES of {(rm_hs_dg[1] - rm_hs_fv[1]):.4f} than the Full Valuation Approach")
else:
    print(f"The Full Valuation Approach gives a higher ES of {(rm_hs_fv[1] - rm_hs_dg[1]):.4f} than the Delta-Gamma Approximation method")

The Delta-Gamma Approximation method gives a higher VaR of 0.0689 than the Full Valuation Approach
The Full Valuation Approach gives a higher ES of 0.1384 than the Delta-Gamma Approximation method

Comments on the results of Stage 2:

• The Delta-Gamma Loss Approximation method gives a higher VaR given that it freezes the Greeks (theta, delta, gamma) while in the Full Valuation method, the Greeks change with the stock value and with the passage of time. Moreoever, for big moves in losses, the Delta-Gamma Approximation approach can exaggerate large-down moves, so the 99% confidence level may land a little higher than the valuation's.

• The Full Valuation approach gives a higher ES as this approach enforces the option's true convex pricing and time decay, which can make the far tail a bit heavier even if the VaR threshold itself is slightly lower. Furthermore, given the 99% confidence level, a couple of extreme full-valuation repriced outcomes can lift ES more than they can move VaR.

The age-weighted historical simulation approach.

We can see that as we introduce weights which exponentiallly decay the older the observation is, the VaR and ES increase relative to the equal-weights case. The difference (albeit small) tells us that the most recent period has heavier tails than the long-run average as giving more weights to recent data lifts the tail risk measures.

# Calculate VaR and ES using age-weighted historical simulation with decay parameter 0.94
rm_hs_aw = VaR_ES_emp(loss_scenarios, conf_level = conf_level, hs_weight_type = "BRW", age_weight_decay = 0.94)

# Calculate VaR and ES using age-weighted historical simulation with the full valuation approach where decay parameter is 0.94
rm_hs_aw_fv = VaR_ES_emp(loss_scenarios_fv, conf_level = conf_level, hs_weight_type = "BRW", age_weight_decay = 0.94)

# Gather the risk measures from the standard and age-weighted historical simulations
shs_aw = pd.DataFrame({
    'Method' : ['Standard HS', 'Age-weighted HS'],
    'VaR': [f"{rm_hs_aw[0]:.4f}", f"{rm_hs_aw_fv[0]:.4f}"],
    'ES': [f"{rm_hs_aw[1]:.4f}", f"{rm_hs_aw_fv[1]:.4f}"]
})
print(shs_aw)

            Method      VaR       ES
0      Standard HS  11.9942  13.2211
1  Age-weighted HS  12.0524  13.2682

4. Project Stage 3: Monte Carlo Simulation

Run Monte Carlo simulations to estimate the 10-day 99% VaR and ES of the portfolio under different model assumptions, including the Black-Scholes-Merton model and a t-copula model with t-distributed marginals.

Log-likelihood function for the Student's $t$ copula

In the two-dimensional case, the copula density of the Student's $t$ copula is given by

$$c_{\nu, \mathbf{P}}(u_1, u_2) = \frac{\mathcal{f}_{\nu,\mathbf{P}}(F^{-1}_\nu(u_1), F^{-1}_\nu(u_2))}{f_\nu(F^{-1}_\nu(u_1)) f_\nu(F^{-1}_\nu(u_2))}$$

where:

• $f_{\nu, \mathbf{P}}$ is the joint density function of a two-dimensional Student's $t$ distribution with zero mean, correlation matrix $\mathbf{P}$, and $\nu$ degrees of freedom,
• $f_\nu$ is the density function of a univariate Student's $t$ distribution with $\nu$ degrees of freedom (loc $ = 0$, scale $ = 1$)
• $F_\nu^{-1}$ is the inverse of the df corresponding to $f_\nu$

# Construct the negative log-likelihood function for the Student's t copula
def neg_loglik_t(pars, u, v):
    # Parameters for the multivariate t distribution
    nu = pars[0]
    rho = pars[1]
    cov_mat = [[1, rho], [rho, 1]]

    # Compute the inverse of the univariate cdf at the pseudo-observations
    F_inv1 = stats.t.ppf(u, df = nu)
    F_inv2 = stats.t.ppf(v, df = nu)

    # Calculate the log-likelihood
    numerator = stats.multivariate_t.pdf(np.array([F_inv1, F_inv2]).T, loc = None, shape = cov_mat, df = nu, allow_singular = False)
    denominator = stats.t.pdf(F_inv1, df = nu) * stats.t.pdf(F_inv2, df = nu)

    ratio = numerator / denominator
    ratio = ratio[~np.isnan(ratio)]    # In case the process above generates NaNs

    value = -np.sum(np.log(ratio))
    return value

# Form pseudo-observations for the log-returns using the empirical marginal distributions
rfch_rank = stats.rankdata(rfch[['AAPL', 'IBM']], axis = 0)
aapl_ps_obs_em = rfch_rank[:,0] / (len(rfch_rank) + 1)
ibm_ps_obs_em = rfch_rank[:,1] / (len(rfch_rank) + 1)

# Find MLE estimates of t copula
init_pars = [2, 0.1]
bnds = ((0, None), (-0.999999, 0.999999))       
tcop_emp_pars = optimize.minimize(neg_loglik_t, init_pars, args = (aapl_ps_obs_em, ibm_ps_obs_em), bounds = bnds)

# Extract outputs
tcop_emp_df = tcop_emp_pars.x[0]
tcop_emp_rho = tcop_emp_pars.x[1]
tcop_emp_cov = [[1, tcop_emp_rho], [tcop_emp_rho, 1]]

print(f'Under the empirical marginal distributions, t copula df: {tcop_emp_df:.4f}')
print(f'Under the empirical marginal distributions, t copula rho: {tcop_emp_rho:.4f}')
print(f'Under the empirical marginal distributions, t copula covariance matrix: \n{np.round(tcop_emp_cov, 4)}')

Under the empirical marginal distributions, t copula df: 7.8926
Under the empirical marginal distributions, t copula rho: 0.3478
Under the empirical marginal distributions, t copula covariance matrix: 
[[1.     0.3478]
 [0.3478 1.    ]]

Tail dependence

# Calculate the coeffcients of upper and lower tail dependence of the log returns
# Empirical copula - AAPL and IBM
em_cop = pycop.empirical(rfch[['AAPL','IBM']].T.values)
utdc_aapl_ibm = em_cop.UTDC(0.99)
ltdc_aapl_ibm = em_cop.LTDC(0.01)

# Non-parametric tail dependence coefficient estimates
print(f'Coefficient of lower tail dependence of log returns: {ltdc_aapl_ibm:.4f}')
print(f'Coefficient of upper tail dependence of log returns: {utdc_aapl_ibm:.4f}')

Coefficient of lower tail dependence of log returns: 0.1399
Coefficient of upper tail dependence of log returns: 0.1663

Since $C$ is the bivariate Student's $t$ copula $C^{t}_{\nu,\rho}$, for some $\rho \in (-1,1]$, then:

$$ \lambda \;=\; 2\, t_{\nu+1}\!\left(-\sqrt{\frac{(\nu+1)(1-\rho)}{1+\rho}}\right) $$

where $t_\nu$ is the density of a univariate $t$ distribution with $\nu$ degrees of freedom. Thus, if $\rho>-1$, $X_1$ and $X_2$ are asymptotically dependent in both the upper and lower tails.

# Calculate the coeffcients of upper and lower tail dependence implied by the bivariate t copula
def tail_dep_t(nu, rho):
    lambda_U = 2 * stats.t.cdf(-np.sqrt((nu + 1) * (1 - rho) / (1 + rho)), df = nu + 1)
    lambda_L = lambda_U
    return lambda_U, lambda_L

tail_dep_tcop = tail_dep_t(tcop_emp_df, tcop_emp_rho)
print(f'Coefficient of upper tail dependence (t copula): {tail_dep_tcop[0]:.4f}')
print(f'Coefficient of lower tail dependence (t copula): {tail_dep_tcop[1]:.4f}')

Coefficient of upper tail dependence (t copula): 0.0682
Coefficient of lower tail dependence (t copula): 0.0682

# Pairwise scatterplot
ps_obs_df = pd.DataFrame(np.column_stack((aapl_ps_obs_em, ibm_ps_obs_em)))
ps_obs_df.columns = ["AAPL", "IBM"]
sns.pairplot(ps_obs_df, plot_kws = dict(alpha = 0.1 , edgecolor = 'k'))

<seaborn.axisgrid.PairGrid at 0x233e3c11310>

# Calculate difference in tail dependence coefficients are consistent with the historical data
print(f'Difference in lower tail dependence coefficients: {(ltdc_aapl_ibm - tail_dep_tcop[1]):.4f}')
print(f'Difference in upper tail dependence coefficients: {(utdc_aapl_ibm - tail_dep_tcop[0]):.4f}')

Difference in lower tail dependence coefficients: 0.0716
Difference in upper tail dependence coefficients: 0.0980

While the $t$-copula forces symmetric upper/lower tail dependence, the data show stronger and asymmetric tail co-movement (upper TD > lower TD). We can conclude that the fitted $t$-copula is underestimating tail dependence when compared to the historical data.

Monte Carlo VaR/ES with the Student’s $t$ copula

# Number of simulations
N_sim = 10000

# Apple log-return mean and standard deviation
aapl_rfch_mean = rfch_mean['AAPL']
aapl_rfch_sd = np.sqrt(rfch_cov_ewma[0,0])

# IBM log-return mean and standard deviation
ibm_rfch_mean = rfch_mean['IBM']
ibm_rfch_sd = np.sqrt(rfch_cov_ewma[1,1])

# Simulate from the t copula
mult_t_sim = stats.multivariate_t.rvs(loc = None, shape = tcop_emp_cov, df = tcop_emp_df, size = N_sim)

# Evaluate the univariate t cdf at the simulated values
tcop_sim = stats.t.cdf(mult_t_sim, df = tcop_emp_df)

# Apply quantile transform to generate simulated log-returns
aapl_rfch_tcop_sim = stats.norm.ppf(tcop_sim[:,0], loc = aapl_rfch_mean, scale = aapl_rfch_sd)
ibm_rfch_tcop_sim = stats.norm.ppf(tcop_sim[:,1], loc = ibm_rfch_mean, scale = ibm_rfch_sd)

# Compile into a single array
rfch_sim_tcop = np.array([aapl_rfch_tcop_sim, ibm_rfch_tcop_sim]).T

Loss Scenarios and Risk Measures under the Delta-Gamma Loss Approximation using the $t$ Copula

# Loss scenarios - for-loop approach
loss_scenarios_mc_dg_tcop = np.zeros(len(rfch_sim_tcop))

for i in range(len(rfch_sim_tcop)):
        loss_scenarios_mc_dg_tcop[i] = - time_decay - np.matmul(rfch_sim_tcop[i,:], portfolio_exposure) - 0.5 * np.matmul(rfch_sim_tcop[i,:], np.matmul(portfolio_exposure_SO, rfch_sim_tcop[i,:].T))
print(loss_scenarios_mc_dg_tcop)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(loss_scenarios_mc_dg_tcop, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
rm_mc_dg_tcop = VaR_ES_emp(loss_scenarios_mc_dg_tcop, conf_level = conf_level, hs_weight_type = "equal")

print('Risk measures using Monte Carlo simulation with t copula and Delta-Gamma approximation:')
print(f'Portfolio VaR: {rm_mc_dg_tcop[0]:.4f}')
print(f'Portfolio ES: {rm_mc_dg_tcop[1]:.4f}')

[ -4.60431667 -15.91770899  -8.0179341  ...   0.86718708  -2.30733299
   9.98013658]

Risk measures using Monte Carlo simulation with t copula and Delta-Gamma approximation:
Portfolio VaR: 22.5522
Portfolio ES: 24.5085

Risk Measures under the Full Valuation Approach using the $t$ Copula

# Risk factor scenarios
aapl_price_scenarios_mc_fv_tcop = aapl_price * np.exp(aapl_rfch_tcop_sim)
ibm_price_scenarios_mc_fv_tcop = ibm_price * np.exp(ibm_rfch_tcop_sim)

# Option price scenarios
aapl_call_scenarios_mc_fv_tcop = BSprice(aapl_price_scenarios_mc_fv_tcop, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
ibm_call_scenarios_mc_fv_tcop = BSprice(ibm_price_scenarios_mc_fv_tcop, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
portfolio_value_mc_fv_tcop = aapl_call_price + ibm_call_price - aapl_price - ibm_price
portfolio_scenarios_mc_fv_tcop = aapl_call_scenarios_mc_fv_tcop + ibm_call_scenarios_mc_fv_tcop - aapl_price_scenarios_mc_fv_tcop - ibm_price_scenarios_mc_fv_tcop
loss_scenarios_mc_fv_tcop = -(portfolio_scenarios_mc_fv_tcop - portfolio_value_mc_fv_tcop)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(loss_scenarios_mc_fv_tcop, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
rm_mc_fv_tcop = VaR_ES_emp(loss_scenarios_mc_fv_tcop, conf_level = conf_level, hs_weight_type = "equal")
print("Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:")
print(f'Portfolio VaR: {rm_mc_fv_tcop[0]:.4f}')
print(f'Portfolio ES: {rm_mc_fv_tcop[1]:.4f}')

Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:
Portfolio VaR: 22.3543
Portfolio ES: 24.3266

The two methods are close, but the delta–gamma approximation approach gives slightly more conservative results (higher VaR/ES) than the full valuation approach. The time decay effect may contribute to the difference.

Monte Carlo VaR/ES with the bivariate Gaussian distribution

Risk Measures under the Delta-Gamma Loss Approximation using bivariate Gaussian distribution

# Parameters of joint Gaussian distribution for log-returns
log_ret_mean_vec = np.array([aapl_rfch_mean, ibm_rfch_mean])
log_ret_cov_mat = rfch_cov_ewma     # Instead of rfch_cov

# Simulate log-returns from multivariate Gaussian distribution
log_ret_sim_Gauss = stats.multivariate_normal.rvs(mean = log_ret_mean_vec, cov = log_ret_cov_mat, size = N_sim)

# Calculate loss scenarios using the delta-gamma loss approximation
portfolio_loss_mc_dg_Gauss = - time_decay - log_ret_sim_Gauss @ portfolio_exposure - 0.5 * np.diag(log_ret_sim_Gauss @ portfolio_exposure_SO @ log_ret_sim_Gauss.T)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(portfolio_loss_mc_dg_Gauss, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Determine the VaR and ES using empirical distribution (equiv to standard historical simulation)
portfolio_risk_measures_mc_dg_Gauss = VaR_ES_emp(portfolio_loss_mc_dg_Gauss, conf_level = conf_level, hs_weight_type = "equal")
print("Risk measures using Monte Carlo simulation with Gaussian copula and Delta-Gamma approximation:")
print(f'Portfolio VaR: {portfolio_risk_measures_mc_dg_Gauss[0]:.4f}')
print(f'Portfolio ES: {portfolio_risk_measures_mc_dg_Gauss[1]:.4f}')

Risk measures using Monte Carlo simulation with Gaussian copula and Delta-Gamma approximation:
Portfolio VaR: 19.0655
Portfolio ES: 20.9675

Risk Measures Using the Full Valuation Approach

Using bivariate Gaussian distribution

# Risk factor scenarios
aapl_price_scenarios_mc_fv_Gauss = aapl_price * np.exp(log_ret_sim_Gauss[:,0])
ibm_price_scenarios_mc_fv_Gauss = ibm_price * np.exp(log_ret_sim_Gauss[:,1])

# Option price scenarios
aapl_call_scenarios_mv_fv_Gauss = BSprice(aapl_price_scenarios_mc_fv_Gauss, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
ibm_call_scenarios_mv_fv_Gauss = BSprice(ibm_price_scenarios_mc_fv_Gauss, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
portfolio_value_mc_fv_Gauss = aapl_call_price + ibm_call_price - aapl_price - ibm_price
portfolio_scenarios_mc_fv_Gauss = aapl_call_scenarios_mv_fv_Gauss + ibm_call_scenarios_mv_fv_Gauss - aapl_price_scenarios_mc_fv_Gauss - ibm_price_scenarios_mc_fv_Gauss
loss_scenarios_mc_fv_Gauss = -(portfolio_scenarios_mc_fv_Gauss - portfolio_value_mc_fv_Gauss)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(loss_scenarios_mc_fv_Gauss, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
rm_mc_fv_Gauss = VaR_ES_emp(loss_scenarios_mc_fv_Gauss, conf_level = conf_level, hs_weight_type = "equal")
print("Risk measures using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:")
print(f'Portfolio VaR: {rm_mc_fv_Gauss[0]:.4f}')
print(f'Portfolio ES: {rm_mc_fv_Gauss[1]:.4f}')

Risk measures using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:
Portfolio VaR: 18.9505
Portfolio ES: 20.8049

# Summarise the results of all the approaches in a DataFrame
results_df = pd.DataFrame({
    'Method': ['Analytical', 'Analytical', 'D-G', 'FV', 'D-G', 'FV', 'D-G', 'FV'],
    'Approach': ['Var-Cov (Undiversified)', 'Var-Cov (Diversified)', 'SHS', 'SHS', 'A-WHS', 'A-WHS', 'MC (t copula)', 'MC (t copula)'],
    'VaR': [portfolio_risk_measures[0], portfolio_risk_measures_d[0], rm_hs_dg[0], rm_hs_fv[0], rm_hs_aw[0], rm_hs_aw_fv[0], rm_mc_dg_tcop[0], rm_mc_fv_tcop[0]],
    'ES': [portfolio_risk_measures[1], portfolio_risk_measures_d[1], rm_hs_dg[1], rm_hs_fv[1], rm_hs_aw[1], rm_hs_aw_fv[1], rm_mc_dg_tcop[1], rm_mc_fv_tcop[1]]
})
results_df['VaR'] = results_df['VaR'].map('{:.4f}'.format)
results_df['ES'] = results_df['ES'].map('{:.4f}'.format)
display(results_df)

	Method	Approach	VaR	ES
0	Analytical	Var-Cov (Undiversified)	33.9344	38.3756
1	Analytical	Var-Cov (Diversified)	23.6119	26.5496
2	D-G	SHS	20.0579	23.0180
3	FV	SHS	19.9890	23.1564
4	D-G	A-WHS	11.9942	13.2211
5	FV	A-WHS	12.0524	13.2682
6	D-G	MC (t copula)	22.5522	24.5085
7	FV	MC (t copula)	22.3543	24.3266

Undiversified vs diversified (Var–Cov): The undiversified figures are essentially the sums of standalone risks (ignores correlation), so they are much higher. The diversified Var–Cov uses the covariance matrix and thus lowering the results.

Delta–Gamma vs Full re-valuation: Across the Historical simulation and Monte Carlo simulation methods, the Delta-Gamma Approximation apporoach gives slightly higher results than the Full Valuation approach. It should be expected as the quadratic approximation can overstate losses for large moves and freezes Greeks, while the other approach reprices the options (honors the option floor and changing Greeks).

Standard HS vs Age-Weighted HS: Age-weighted is much lower here. This tells us that the most recent window is calmer than the long-run sample as giving more weight to recent observations shrinks tail estimates.

MC $t$-copula vs HS: The Monte Carlo simulation with a t-copula lands above the Historical simulation. The reason is that the $t$-copula injects positive tail dependence ($\lambda >0$), making joint extremes more frequent than HS’s finite historical tail.

5. Project Stage 4: Stressed VaR and ES

Estimate the 10-day 99% Stressed VaR and Stressed ES of the portfolio using different approaches for different stressed periods

Using Analytical approach to calculate the stressed VaR and ES for different stressed periods

# Extract financial stress period data from Stress Period 1: 03-12-2007 to 30-06-2009
sp1_rfch = rfch.loc['2007-12-03':'2009-06-30']

# Mean and covariance matrix of stressed risk factor changes from the sample
sp1_rfch_mean = sp1_rfch.mean()
sp1_rfch_cov_sample = np.cov(sp1_rfch, rowvar = False)

# Marginal distribution parameters of Apple log-returns
sp1_aapl_mean = sp1_rfch_mean['AAPL']
sp1_aapl_sd = np.sqrt(sp1_rfch_cov_sample[0,0])

# Marginal distribution parameters of IBM log-returns
sp1_ibm_mean = sp1_rfch_mean['IBM']
sp1_ibm_sd = np.sqrt(sp1_rfch_cov_sample[1,1])

# Portfolio-level stressed VaR and ES
sp1_rm_analytical = VaR_ES_vcv(time_decay, portfolio_exposure, sp1_rfch_mean, sp1_rfch_cov_sample, conf_level)
print('For Stress Period 1 and using Analytical Approach')
print(f'Portfolio-level stressed VaR and ES (diversified) using analytical approach: \n{np.round(sp1_rm_analytical, 4)}')

For Stress Period 1 and using Analytical Approach
Portfolio-level stressed VaR and ES (diversified) using analytical approach: 
[37.4988 42.6964]

# Extract financial stress period data from Stress Period 2: 30-01-2020 to 01-02-2021
sp2_rfch = rfch.loc['2020-01-30':'2021-02-01']

# Mean and covariance matrix of stressed risk factor changes from the sample
sp2_rfch_mean = sp2_rfch.mean()
sp2_rfch_cov_sample = np.cov(sp2_rfch, rowvar = False)

# Marginal distribution parameters of Apple log-returns
sp2_aapl_mean = sp2_rfch_mean['AAPL']
sp2_aapl_sd = np.sqrt(sp2_rfch_cov_sample[0,0])

# Marginal distribution parameters of IBM log-returns
sp2_ibm_mean = sp2_rfch_mean['IBM']
sp2_ibm_sd = np.sqrt(sp2_rfch_cov_sample[1,1])

# Portfolio-level stressed VaR and ES
sp2_rm_analytical = VaR_ES_vcv(time_decay, portfolio_exposure, sp2_rfch_mean, sp2_rfch_cov_sample, conf_level)
print('For Stress Period 2 and using Analytical Approach')
print(f'Portfolio-level stressed VaR: {sp2_rm_analytical[0]:.4f}')
print(f'Portfolio-level stressed ES: {sp2_rm_analytical[1]:.4f}')

For Stress Period 2 and using Analytical Approach
Portfolio-level stressed VaR: 36.4899
Portfolio-level stressed ES: 41.1845

Using Standard historical simulation approach with full valuation to calculate the stressed VaR and ES for different stressed periods

Stress Period 1: 03-12-2007 to 30-06-2009

# Risk factor scenarios
sp1_aapl_price_scenarios = aapl_price * np.exp(sp1_rfch['AAPL'])
sp1_ibm_price_scenarios = ibm_price * np.exp(sp1_rfch['IBM'])

# Option price scenarios
sp1_aapl_call_scenarios = BSprice(sp1_aapl_price_scenarios, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp1_ibm_call_scenarios = BSprice(sp1_ibm_price_scenarios, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp1_portfolio_value = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp1_portfolio_scenarios = sp1_aapl_call_scenarios + sp1_ibm_call_scenarios - sp1_aapl_price_scenarios - sp1_ibm_price_scenarios
sp1_loss_scenarios = -(sp1_portfolio_scenarios - sp1_portfolio_value)

# Calculate risk measures empirical loss distribution
sp1_hs_fv = VaR_ES_emp(sp1_loss_scenarios, conf_level = conf_level, hs_weight_type = "equal")

print("For Stress Period 1 using full valuation approach:")
print(f"Portfolio-level stressed VaR (diversified) : {sp1_hs_fv[0]:.4f}")
print(f"Portfolio-level stressed ES (diversified) : {sp1_hs_fv[1]:.4f}")

For Stress Period 1 using full valuation approach:
Portfolio-level stressed VaR (diversified) : 22.4713
Portfolio-level stressed ES (diversified) : 23.6675

Stress Period 2: 30-01-2020 to 01-02-2021

# Risk factor scenarios
sp2_aapl_price_scenarios = aapl_price * np.exp(sp2_rfch['AAPL'])
sp2_ibm_price_scenarios = ibm_price * np.exp(sp2_rfch['IBM'])

# Option price scenarios
sp2_aapl_call_scenarios = BSprice(sp2_aapl_price_scenarios, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp2_ibm_call_scenarios = BSprice(sp2_ibm_price_scenarios, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp2_portfolio_value = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp2_portfolio_scenarios = sp2_aapl_call_scenarios + sp2_ibm_call_scenarios - sp2_aapl_price_scenarios - sp2_ibm_price_scenarios
sp2_loss_scenarios = -(sp2_portfolio_scenarios - sp2_portfolio_value)

# Calculate risk measures empirical loss distribution
sp2_hs_fv = VaR_ES_emp(sp2_loss_scenarios, conf_level = conf_level, hs_weight_type = "equal")

print("For Stress Period 2 using full valuation approach:")
print(f"Portfolio-level stressed VaR (diversified) : {sp2_hs_fv[0]:.4f}")
print(f"Portfolio-level stressed ES (diversified) : {sp2_hs_fv[1]:.4f}")

For Stress Period 2 using full valuation approach:
Portfolio-level stressed VaR (diversified) : 20.6504
Portfolio-level stressed ES (diversified) : 23.2192

Monte Carlo simulation approach with full valuation under the Student's $t$ copula

Stress Period 1: 03-12-2007 to 30-06-2009

# Form pseudo-observations for the log-returns using the empirical marginal distributions
sp1_rfch_rank = stats.rankdata(sp1_rfch[['AAPL', 'IBM']], axis = 0)
sp1_aapl_ps_obs_em = sp1_rfch_rank[:,0] / (len(sp1_rfch_rank) + 1)
sp1_ibm_ps_obs_em = sp1_rfch_rank[:,1] / (len(sp1_rfch_rank) + 1)

# Function of MLE estimates of t copula
def MLE_tcop(ps_obs1, ps_obs2):
    # Initial parameter values and bounds
    mle_init_pars = [2, 0.1]
    mle_bnds = ((0, None), (-0.999999, 0.999999))
    mle_tcop_emp_pars = optimize.minimize(neg_loglik_t, mle_init_pars, args = (ps_obs1, ps_obs2), bounds = mle_bnds)

    # Extract outputs
    mle_tcop_emp_df = mle_tcop_emp_pars.x[0]
    mle_tcop_emp_rho = mle_tcop_emp_pars.x[1]
    mle_tcop_emp_cov = [[1, mle_tcop_emp_rho], [mle_tcop_emp_rho, 1]]
    return mle_tcop_emp_df, mle_tcop_emp_rho, mle_tcop_emp_cov

# MLE estimates of t copula parameters for Stress Period 1
sp1_MLE_tcop = MLE_tcop(sp1_aapl_ps_obs_em, sp1_ibm_ps_obs_em)
print('Under the empirical marginal distributions,')
print(f't copula df: {sp1_MLE_tcop[0]:.4f}')
print(f't copula rho: {sp1_MLE_tcop[1]:.4f}')
print(f't copula covariance matrix: \n{np.round(sp1_MLE_tcop[2], 4)}')

# Simulate from the t copula
sp1_mult_t_sim = stats.multivariate_t.rvs(loc = None, shape = sp1_MLE_tcop[2], df = sp1_MLE_tcop[0], size = N_sim)

# Evaluate the univariate t cdf at the simulated values
sp1_tcop_sim = stats.t.cdf(sp1_mult_t_sim, df = sp1_MLE_tcop[0])

# Apply quantile transform to generate simulated log-returns
sp1_aapl_rfch_tcop_sim = stats.norm.ppf(sp1_tcop_sim[:,0], loc = sp1_aapl_mean, scale = sp1_aapl_sd)
sp1_ibm_rfch_tcop_sim = stats.norm.ppf(sp1_tcop_sim[:,1], loc = sp1_ibm_mean, scale = sp1_ibm_sd)

# Compile into a single array
sp1_rfch_sim_tcop = np.array([sp1_aapl_rfch_tcop_sim, sp1_ibm_rfch_tcop_sim]).T

# Risk factor scenarios
sp1_aapl_price_scenarios_mc_fv_tcop = aapl_price * np.exp(sp1_aapl_rfch_tcop_sim)
sp1_ibm_price_scenarios_mc_fv_tcop = ibm_price * np.exp(sp1_ibm_rfch_tcop_sim)

# Option price scenarios
sp1_aapl_call_scenarios_mc_fv_tcop = BSprice(sp1_aapl_price_scenarios_mc_fv_tcop, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp1_ibm_call_scenarios_mc_fv_tcop = BSprice(sp1_ibm_price_scenarios_mc_fv_tcop, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp1_portfolio_value_mc_fv_tcop = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp1_portfolio_scenarios_mc_fv_tcop = sp1_aapl_call_scenarios_mc_fv_tcop + sp1_ibm_call_scenarios_mc_fv_tcop - sp1_aapl_price_scenarios_mc_fv_tcop - sp1_ibm_price_scenarios_mc_fv_tcop
sp1_loss_scenarios_mc_fv_tcop = -(sp1_portfolio_scenarios_mc_fv_tcop - sp1_portfolio_value_mc_fv_tcop)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(sp1_loss_scenarios_mc_fv_tcop, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios of Stress Period 1 (GFC)')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
sp1_rm_mc_fv_tcop = VaR_ES_emp(sp1_loss_scenarios_mc_fv_tcop, conf_level = conf_level, hs_weight_type = "equal")
print("Stress Period 1 - Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:")
print(f'Portfolio VaR: {sp1_rm_mc_fv_tcop[0]:.4f}')
print(f'Portfolio ES: {sp1_rm_mc_fv_tcop[1]:.4f}')

Under the empirical marginal distributions,
t copula df: 2.8737
t copula rho: 0.5100
t copula covariance matrix: 
[[1.   0.51]
 [0.51 1.  ]]

Stress Period 1 - Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:
Portfolio VaR: 23.4262
Portfolio ES: 25.0646

Stress Period 2: 30-01-2020 to 01-02-2021

# Form pseudo-observations for the log-returns using the empirical marginal distributions
sp2_rfch_rank = stats.rankdata(sp2_rfch[['AAPL', 'IBM']], axis = 0)
sp2_aapl_ps_obs_em = sp2_rfch_rank[:,0] / (len(sp2_rfch_rank) + 1)
sp2_ibm_ps_obs_em = sp2_rfch_rank[:,1] / (len(sp2_rfch_rank) + 1)

# MLE estimates of t copula parameters for Stress Period 2
sp2_MLE_tcop = MLE_tcop(sp2_aapl_ps_obs_em, sp2_ibm_ps_obs_em)      # the function is defined above for stress period 1
print('Under the empirical marginal distributions,')
print(f't copula df: {sp2_MLE_tcop[0]:.4f}')
print(f't copula rho: {sp2_MLE_tcop[1]:.4f}')
print(f't copula covariance matrix: \n{np.round(sp2_MLE_tcop[2], 4)}')

# Simulate from the t copula
sp2_mult_t_sim = stats.multivariate_t.rvs(loc = None, shape = sp2_MLE_tcop[2], df = sp2_MLE_tcop[0], size = N_sim)

# Evaluate the univariate t cdf at the simulated values
sp2_tcop_sim = stats.t.cdf(sp2_mult_t_sim, df = sp2_MLE_tcop[0])

# Apply quantile transform to generate simulated log-returns
sp2_aapl_rfch_tcop_sim = stats.norm.ppf(sp2_tcop_sim[:,0], loc = sp2_aapl_mean, scale = sp2_aapl_sd)
sp2_ibm_rfch_tcop_sim = stats.norm.ppf(sp2_tcop_sim[:,1], loc = sp2_ibm_mean, scale = sp2_ibm_sd)

# Compile into a single array
sp2_rfch_sim_tcop = np.array([sp2_aapl_rfch_tcop_sim, sp2_ibm_rfch_tcop_sim]).T

# Risk factor scenarios
sp2_aapl_price_scenarios_mc_fv_tcop = aapl_price * np.exp(sp2_aapl_rfch_tcop_sim)
sp2_ibm_price_scenarios_mc_fv_tcop = ibm_price * np.exp(sp2_ibm_rfch_tcop_sim)

# Option price scenarios
sp2_aapl_call_scenarios_mc_fv_tcop = BSprice(sp2_aapl_price_scenarios_mc_fv_tcop, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp2_ibm_call_scenarios_mc_fv_tcop = BSprice(sp2_ibm_price_scenarios_mc_fv_tcop, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp2_portfolio_value_mc_fv_tcop = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp2_portfolio_scenarios_mc_fv_tcop = sp2_aapl_call_scenarios_mc_fv_tcop + sp2_ibm_call_scenarios_mc_fv_tcop - sp2_aapl_price_scenarios_mc_fv_tcop - sp2_ibm_price_scenarios_mc_fv_tcop
sp2_loss_scenarios_mc_fv_tcop = -(sp2_portfolio_scenarios_mc_fv_tcop - sp2_portfolio_value_mc_fv_tcop)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(sp2_loss_scenarios_mc_fv_tcop, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios of Stress Period 2 (Covid-19)')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
sp2_rm_mc_fv_tcop = VaR_ES_emp(sp2_loss_scenarios_mc_fv_tcop, conf_level = conf_level, hs_weight_type = "equal")
print("Stress Period 2 - Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:")
print(f'Portfolio VaR: {sp2_rm_mc_fv_tcop[0]:.4f}')
print(f'Portfolio ES: {sp2_rm_mc_fv_tcop[1]:.4f}')

Under the empirical marginal distributions,
t copula df: 32.5851
t copula rho: 0.3542
t copula covariance matrix: 
[[1.     0.3542]
 [0.3542 1.    ]]

Stress Period 2 - Risk measures using Monte Carlo simulation with t copula and Full Valuation Approach:
Portfolio VaR: 24.0218
Portfolio ES: 25.7071

Monte Carlo simulation approach with full valuation under the bivariate Gaussian distribution

Stress Period 1: 03-12-2007 to 30-06-2009

# Parameters of joint Gaussian distribution for log-returns
sp1_rfch_mean_vec = np.array([sp1_rfch_mean['AAPL'], sp1_rfch_mean['IBM']])

# The covariance matrix is already defined as sp1_rfch_cov_sample

# Simulate log-returns from multivariate Gaussian distribution
sp1_rfch_sim_Gauss = stats.multivariate_normal.rvs(mean = sp1_rfch_mean_vec, cov = sp1_rfch_cov_sample, size = N_sim)

# Risk factor scenarios
sp1_aapl_price_scenarios_mc_fv_Gauss = aapl_price * np.exp(sp1_rfch_sim_Gauss[:,0])
sp1_ibm_price_scenarios_mc_fv_Gauss = ibm_price * np.exp(sp1_rfch_sim_Gauss[:,1])

# Option price scenarios
sp1_aapl_call_scenarios_mv_fv_Gauss = BSprice(sp1_aapl_price_scenarios_mc_fv_Gauss, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp1_ibm_call_scenarios_mv_fv_Gauss = BSprice(sp1_ibm_price_scenarios_mc_fv_Gauss, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp1_portfolio_value_mc_fv_Gauss = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp1_portfolio_scenarios_mc_fv_Gauss = sp1_aapl_call_scenarios_mv_fv_Gauss + sp1_ibm_call_scenarios_mv_fv_Gauss - sp1_aapl_price_scenarios_mc_fv_Gauss - sp1_ibm_price_scenarios_mc_fv_Gauss
sp1_loss_scenarios_mc_fv_Gauss = -(sp1_portfolio_scenarios_mc_fv_Gauss - sp1_portfolio_value_mc_fv_Gauss)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(sp1_loss_scenarios_mc_fv_Gauss, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios of Stress Period 1 (GFC)')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
sp1_rm_mc_fv_Gauss = VaR_ES_emp(sp1_loss_scenarios_mc_fv_Gauss, conf_level, hs_weight_type = "equal")
print("Risk measures of Stressed Period 1 using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:")
print(f'Portfolio VaR: {sp1_rm_mc_fv_Gauss[0]:.4f}')
print(f'Portfolio ES: {sp1_rm_mc_fv_Gauss[1]:.4f}')

Risk measures of Stressed Period 1 using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:
Portfolio VaR: 22.7180
Portfolio ES: 24.2098

Stress Period 2: 30-01-2020 to 01-02-2021

# Parameters of joint Gaussian distribution for log-returns
sp2_rfch_mean_vec = np.array([sp2_rfch_mean['AAPL'], sp2_rfch_mean['IBM']])

# The covariance matrix is already defined as sp2_rfch_cov_sample

# Simulate log-returns from multivariate Gaussian distribution
sp2_rfch_sim_Gauss = stats.multivariate_normal.rvs(mean = sp2_rfch_mean_vec, cov = sp2_rfch_cov_sample, size = N_sim)

# Risk factor scenarios
sp2_aapl_price_scenarios_mc_fv_Gauss = aapl_price * np.exp(sp2_rfch_sim_Gauss[:,0])
sp2_ibm_price_scenarios_mc_fv_Gauss = ibm_price * np.exp(sp2_rfch_sim_Gauss[:,1])

# Option price scenarios
sp2_aapl_call_scenarios_mv_fv_Gauss = BSprice(sp2_aapl_price_scenarios_mc_fv_Gauss, aapl_price, 3/12 - delta_t, 0.035, 0, aapl_ann_vol, 'call')
sp2_ibm_call_scenarios_mv_fv_Gauss = BSprice(sp2_ibm_price_scenarios_mc_fv_Gauss, ibm_price, 6/12 - delta_t, 0.035, 0, ibm_ann_vol, 'call')

# Portfolio loss scenarios
sp2_portfolio_value_mc_fv_Gauss = aapl_call_price + ibm_call_price - aapl_price - ibm_price
sp2_portfolio_scenarios_mc_fv_Gauss = sp2_aapl_call_scenarios_mv_fv_Gauss + sp2_ibm_call_scenarios_mv_fv_Gauss - sp2_aapl_price_scenarios_mc_fv_Gauss - sp2_ibm_price_scenarios_mc_fv_Gauss
sp2_loss_scenarios_mc_fv_Gauss = -(sp2_portfolio_scenarios_mc_fv_Gauss - sp2_portfolio_value_mc_fv_Gauss)

# Histogram of loss scenarios
plt.figure(figsize=(8, 5))
plt.hist(sp2_loss_scenarios_mc_fv_Gauss, bins=100, color='skyblue', edgecolor='black')
plt.title('Histogram of Loss Scenarios of Stress Period 2 (Covid-19)')
plt.xlabel('Loss')
plt.ylabel('Frequency')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# Calculate risk measures empirical loss distribution
sp2_rm_mc_fv_Gauss = VaR_ES_emp(sp2_loss_scenarios_mc_fv_Gauss, conf_level, hs_weight_type = "equal")
print("Risk measures of Stressed Period 2 using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:")
print(f'Portfolio VaR: {sp2_rm_mc_fv_Gauss[0]:.4f}')
print(f'Portfolio ES: {sp2_rm_mc_fv_Gauss[1]:.4f}')

Risk measures of Stressed Period 2 using Monte Carlo simulation with Gaussian copula and Full Valuation Approach:
Portfolio VaR: 24.8027
Portfolio ES: 26.4915

Summary of Stage 4 results in relation to Stressed VaR and ES

stress_period_results = pd.DataFrame({
    'Method': ['Analytical', 'SH Simulation', 'MC Simulation (t copula)', 'MC Simulation (Gaussian)'],
    'Stress Period 1 VaR': [sp1_rm_analytical[0], sp1_hs_fv[0], sp1_rm_mc_fv_tcop[0], sp1_rm_mc_fv_Gauss[0]],
    'Stress Period 1 ES': [sp1_rm_analytical[1], sp1_hs_fv[1], sp1_rm_mc_fv_tcop[1], sp1_rm_mc_fv_Gauss[1]],
    'Stress Period 2 VaR': [sp2_rm_analytical[0], sp2_hs_fv[0], sp2_rm_mc_fv_tcop[0], sp2_rm_mc_fv_Gauss[0]],
    'Stress Period 2 ES': [sp2_rm_analytical[1], sp2_hs_fv[1], sp2_rm_mc_fv_tcop[1], sp2_rm_mc_fv_Gauss[1]]
})
print("Table 1: Stressed 10-day 99% VaR and ES using log-returns during Stress Periods 1 and 2")
display(stress_period_results)

Table 1: Stressed 10-day 99% VaR and ES using log-returns during Stress Periods 1 and 2

	Method	Stress Period 1 VaR	Stress Period 1 ES	Stress Period 2 VaR	Stress Period 2 ES
0	Analytical	37.498796	42.696374	36.489895	41.184499
1	SH Simulation	22.471338	23.667521	20.650395	23.219213
2	MC Simulation (t copula)	23.426151	25.064626	24.021811	25.707079
3	MC Simulation (Gaussian)	22.717951	24.209754	24.802650	26.491505

# Calculate and print the tail dependence coefficients and Pearson correlation for Stress Period 1
sp1_check_nu, sp1_check_rho = sp1_MLE_tcop[0], sp1_MLE_tcop[1]
lambda_u = 2 * stats.t.cdf(-np.sqrt((sp1_check_nu+1)*(1-sp1_check_rho)/(1+sp1_check_rho)), df=sp1_check_nu+1)  # t-copula upper tail dep
print(f"[SP1] t-copula: nu={sp1_check_nu:.4f}, rho={sp1_check_rho:.4f}, lambda_u={lambda_u:.4f}")

# Pearson correlation for Stress Period 1
pearson_rho = sp1_rfch[['AAPL','IBM']].corr().iloc[0,1]
print(f"[SP1] Gaussian Pearson rho={pearson_rho:.4f}")

# Calculate and print the tail dependence coefficients and Pearson correlation for Stress Period 2
sp2_check_nu, sp2_check_rho = sp2_MLE_tcop[0], sp2_MLE_tcop[1]
lambda_u = 2 * stats.t.cdf(-np.sqrt((sp2_check_nu+1)*(1-sp2_check_rho)/(1+sp2_check_rho)), df=sp2_check_nu+1)  # t-copula upper tail dep
print(f"[SP2] t-copula: nu={sp2_check_nu:.4f}, rho={sp2_check_rho:.4f}, lambda_u={lambda_u:.4f}")

# Pearson correlation for Stress Period 2
pearson_rho = sp2_rfch[['AAPL','IBM']].corr().iloc[0,1]
print(f"[SP2] Gaussian Pearson rho={pearson_rho:.4f}")

[SP1] t-copula: nu=2.8737, rho=0.5100, lambda_u=0.3268
[SP1] Gaussian Pearson rho=0.5298
[SP2] t-copula: nu=32.5851, rho=0.3542, lambda_u=0.0003
[SP2] Gaussian Pearson rho=0.5284

# Compare the risk measures obtained using the stress period data with those obtained using the full sample data in Section 4.2
print("Comparison of Risk Measures between Full Sample and Stress Periods:")
comparison_df = pd.DataFrame({
    'Method': ['Analytical', 'SH Simulation', 'MC Simulation (t copula)', 'MC Simulation (Gaussian)'],
    'Full Sample VaR': [portfolio_risk_measures_d[0], rm_hs_fv[0], rm_mc_fv_tcop[0], rm_mc_fv_Gauss[0]],
    'Full Sample ES': [portfolio_risk_measures_d[1], rm_hs_fv[1], rm_mc_fv_tcop[1], rm_mc_fv_Gauss[1]],
    'Stress Period 1 VaR': [sp1_rm_analytical[0], sp1_hs_fv[0], sp1_rm_mc_fv_tcop[0], sp1_rm_mc_fv_Gauss[0]],
    'Stress Period 1 ES': [sp1_rm_analytical[1], sp1_hs_fv[1], sp1_rm_mc_fv_tcop[1], sp1_rm_mc_fv_Gauss[1]],
    'Stress Period 2 VaR': [sp2_rm_analytical[0], sp2_hs_fv[0], sp2_rm_mc_fv_tcop[0], sp2_rm_mc_fv_Gauss[0]],
    'Stress Period 2 ES': [sp2_rm_analytical[1], sp2_hs_fv[1], sp2_rm_mc_fv_tcop[1], sp2_rm_mc_fv_Gauss[1]]
})
display(comparison_df)

Comparison of Risk Measures between Full Sample and Stress Periods:

	Method	Full Sample VaR	Full Sample ES	Stress Period 1 VaR	Stress Period 1 ES	Stress Period 2 VaR	Stress Period 2 ES
0	Analytical	23.611937	26.549553	37.498796	42.696374	36.489895	41.184499
1	SH Simulation	19.989002	23.156404	22.471338	23.667521	20.650395	23.219213
2	MC Simulation (t copula)	22.354283	24.326579	23.426151	25.064626	24.021811	25.707079
3	MC Simulation (Gaussian)	18.950455	20.804926	22.717951	24.209754	24.802650	26.491505

Analytical: This method consistently delivers the largest VaR/ES in our tables as it linearly approximates the losses and assumes Normally distributed risk-factor changes. In crisis windows with very high volatility, the Normal quantile multiplier plus large variances make the analytical measure conservative relative to HS/MC. This does not exclude the possibility that the method overstates the measures. Thus, using the risk measures calculated from this method is not appropriate as stressed risk measure estimates.

Age-weighted HS in stress windows: We do not use age-weighted HS inside the already-restricted stress periods. These windows are short and already focus on stressed days; applying an additional decay would make the result noisy without clear benefit.

Standard Historical Simulation: We can observe that this method, most of the time, gives the lowest VaR/ES. Despite its simplicity, its risk measures are not recommended to be used as stressed risk measure estimates as it ignores parametric distribution.

Monte Carlo (full valuation) — Gaussian vs $t$-copula:

• In the full sample and in Stress Period 1 (GFC), the $t$-copula produces higher VaR/ES than the Gaussian copula. This is expected as there exists positive tail dependence given that $\nu$ > 0 but not large and $\lambda$ > 0. Hence, it is appropriate to use the VaR/ES from the MC $t$-copula for Stressed Period 1.

• In Stress Period 2 (COVID-19) it can be observed that the $t$-copula VaR/ES is below the Gaussian’s. The fitted t-copula there has large $\nu$ (near-Gaussian) and little or no tail dependence ($\lambda$ = 0 ). In that case, the dependence looks essentially Gaussian, and if the Gaussian fit also has a higher correlation, its VaR/ES can exceed the $t$-copula’s. It is also supported by Cristiano Villa, Francisco J. Rubio (2018) that as $\nu$ tends to $+\infty$, the multivariate $t$ distribution converges to the multivariate Normal distribution. A large estimated $\nu$ therefore explains why COVID-period results look Gaussian. Here, in this stressed period, it is appropriate to use the VaR/ES from the MC Gaussian as it gives higher VaR/ES than the $t$-copula.

6. Project Stage 5: A Simple Risk Minimisation Strategy

Determine the number of shares that are required to minimise the risk of the portfolio.

Stock position and Delta Hedging

The Portfolio value is: $V(t, Z_{1,t}, Z_{2,t}) = C_1 (t,Z_{1,t}) + C_2 (t,Z_{2,t}) - \eta_1 e^{Z_{1,t}} - \eta_2 e^{Z_{2,t}}$ , where the number of shares $\eta_1 > 0$ and $\eta_2 > 0$

The delta gamma loss approximation is defined as follows: $$ L^{\Delta,\Gamma}_{t+1} = - (\Theta_{1,\text{call}} +\Theta_{2,\text{call}}) \Delta t - \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1, \text{call}} - \eta_1) \\ e^{Z_{2,t}}(\Delta_{2, \text{call}} - \eta_2) \\ \end{bmatrix}^\top \begin{bmatrix} X_{1, t+1} \\ X_{2, t+1} \end{bmatrix} $$ $$- \tfrac{1}{2} \begin{bmatrix} X_{1, t+1} \\ X_{2, t+1} \end{bmatrix}^\top \begin{bmatrix} e^{Z_{1,t}}(e^{Z_{1,t}} \Gamma_{1,\text{call}} + \Delta_{1,\text{call}} - \eta_1) & 0 \\ 0 & e^{Z_{2,t}}(e^{Z_{2,t}} \Gamma_{2,\text{call}} + \Delta_{2,\text{call}} - \eta_2) \end{bmatrix} \begin{bmatrix} X_{1, t+1} \\ X_{2, t+1} \end{bmatrix} $$

Where the delta loss approximation is defined as:

$$L^{\Delta}_{t+1} = -\mathbf{c}_t - \mathbf{b}^\top_t X_{i, t+1} $$

$$L^{\Delta}_{t+1} = - (\Theta_{1,\text{call}} +\Theta_{2,\text{call}}) \Delta t - \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1, \text{call}} - \eta_1) \\ e^{Z_{2,t}}(\Delta_{2, \text{call}} - \eta_2) \\ \end{bmatrix}^\top \begin{bmatrix} X_{1, t+1} \\ X_{2, t+1} \end{bmatrix} $$

Suppose $X_{t+1}\sim N_d(\mu_{t+1}, \Sigma_{t+1}) $ , we then have $L^{\Delta}_{t+1}\sim N({\mu}_L, {\sigma}^2_{L})$ where

$$\mu_L := E[L_{T+1}] = -\mathbf{c}_t - \mathbf{b}^\top_t\mu_{t+1}$$ $$\sigma^2_L := Var(L_{t+1}) = \mathbf{b}^\top_t \Sigma_{t+1} \mathbf{b}_t$$

To minimise the portfolio loss or the risk of the portfolio, we need minimise the variance $\sigma^2_L$. Since $\Sigma$ is possitive definite, the quadratic form is minised at:

$$\boxed{\mathbf{b}_t = 0}$$

where $\mathbf{b}_t = \begin{bmatrix} \frac{\partial V}{\partial Z_1} \\ \frac{\partial V}{\partial Z_2} \end{bmatrix} = \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1, \text{call}} - \eta_1) \\ e^{Z_{2,t}}(\Delta_{2, \text{call}} - \eta_2) \\ \end{bmatrix}$

The stock position must be chosen so that: $$ \frac{\partial V}{\partial Z_1} = e^{Z_{1,t}}\Delta_{1,\text{call}} - e^{Z_{1,t}}\eta_1 = 0 \implies \eta_1^\star = \Delta_{1,\text{call}} $$

$$ \frac{\partial V}{\partial Z_2} = e^{Z_{2,t}}\Delta_{2,\text{call}} - e^{Z_{2,t}}\eta_2 = 0 \implies \eta_2^\star = \Delta_{2,\text{call}} $$

This is essentially the same as Delta hedging (Alòs, E., & Merino, R., 2023) where having a hedge position would reduce the uncertainty of the portfolio loss compared to an unhedged position (Chiarella, C., He, X.-Z., & Sklibosios Nikitopoulos, C., 2015)

Optimal stock position using Analytical approach

At 99% confidence level, where $\alpha \in (0,1)$, the VaR and ES are calculated as follows: $${\text{VaR}}_{\alpha}(L_{t+1}) = \mu_L + \sigma_L \Phi^{-1}(\alpha), \quad{\text{ES}}_{\alpha}(L_{t+1}) = \mu_L + \sigma_L \frac {{\phi}({\Phi}^{-1}({\alpha}))} {1-{\alpha}}$$ where ${\Phi}$ and ${\phi}$ are the cdf and pdf of standard normal distribution $N(0,1)$

# Delta hedge short positions in the underlying stocks
eta1_star = aapl_call_delta
eta2_star = ibm_call_delta
print(f"By Delta hedging, short Apple stock by {eta1_star:.4f} shares")      
print(f"By Delta hedging, short IBM stock by {eta2_star:.4f} shares")

# Construct a function which takes input the number of shares and returns the 10-day 99% VaR and ES under analytical approach
def minimise_var_es_analytical(eta1, eta2):
    # time-decay term
    c_t = (aapl_call_theta + ibm_call_theta) * delta_t

    # call exposures
    b_t = np.array([
        aapl_price*aapl_call_delta - eta1*aapl_price,
        ibm_price*ibm_call_delta - eta2*ibm_price
    ], dtype=float)

    # mean and standard deviation of the linear loss distribution
    mu_L    = -c_t - b_t @ rfch_mean
    sigma_L = float(np.sqrt(b_t @ rfch_cov_ewma @ b_t))

    return mu_L + norm.ppf(0.99) * sigma_L           # 99% confidence level

# Minimise VaR over eta1, eta2 ≥ 0
def etas_min_hedge_risk(eta1_start=1.0, eta2_start=1.0):        #eta1_start and eta2_start: initial guess
    obj = lambda e: minimise_var_es_analytical(e[0], e[1])
    res = optimize.minimize(obj, x0=[eta1_start, eta2_start], bounds=[(0,None),(0,None)], method="L-BFGS-B")
    return res.x, res.fun

eta_minimise, var_minimise = etas_min_hedge_risk(aapl_call_delta, ibm_call_delta)  # start near deltas
print(f"By Analytical Approach, short {eta_minimise[0]:.4f} Apple shares")
print(f"By Analytical Approach, short {eta_minimise[1]:.4f} IBM shares")
print(f"The minimised VaR is expected to be: {var_minimise:.4f}")

# Compare with the results in Task 1
if eta_minimise[0] - eta1_star == 0:
    print("The optimal position for Apple stock is the same as the Delta hedge")
else:
    print(f"The optimal position for Apple stock is different from the Delta hedge by {(abs(eta_minimise[0] - eta1_star)):.4f} shares")
if eta_minimise[1] - eta2_star == 0:
    print("The optimal position for IBM stock is the same as the Delta hedge")
else:
    print(f"The optimal position for IBM stock is different from the Delta hedge by {(abs(eta_minimise[1] - eta2_star)):.4f} shares")

By Delta hedging, short Apple stock by 0.5527 shares
By Delta hedging, short IBM stock by 0.5778 shares
By Analytical Approach, short 0.5527 Apple shares
By Analytical Approach, short 0.5778 IBM shares
The minimised VaR is expected to be: 2.3002
The optimal position for Apple stock is the same as the Delta hedge
The optimal position for IBM stock is the same as the Delta hedge

Optimal stock position under the standard historical simulation method with a full valuation approach

We would expect the optimal stock position to shift (down) due to gamma and theta effects and historical tail dependence. The optimal stock position is not exactly the same as the delta of the call options, which is consistent with the fact that we are using a full valuation approach and historical simulation method, which captures non-linearities and tail dependence in the data.

# Function for the loss scenarios under full valuation approach
def minimise_loss_scenarios_fv(eta1, eta2, conf_level):

    # portfolio loss scenarios
    minimise_portfolio_value = aapl_call_price + ibm_call_price - eta1*aapl_price - eta2*ibm_price
    minimise_portfolio_scenarios = aapl_call_scenarios + ibm_call_scenarios - eta1*aapl_price_scenarios - eta2*ibm_price_scenarios
    minimise_loss_scenarios_fv = -(minimise_portfolio_scenarios - minimise_portfolio_value)

    # risk measure VaR using standard historical simulation
    var_es_hs_fv = VaR_ES_emp(minimise_loss_scenarios_fv, conf_level=conf_level, hs_weight_type="equal")
    return var_es_hs_fv[0]           # return VaR only
    
# Minimise VaR over eta1, eta2 ≥ 0
def etas_min_hedge_risk_hs_fv(eta1_start=1.0, eta2_start=1.0, conf_level=0.99):        #eta1_start and eta2_start: initial guess
    obj = lambda e: minimise_loss_scenarios_fv(e[0], e[1], conf_level)
    hs_fv = optimize.minimize(obj, x0=[eta1_start, eta2_start], bounds=[(0,None),(0,None)], method="L-BFGS-B")
    return hs_fv.x, hs_fv.fun

eta_minimise_hs_fv, var_minimise_hs_fv = etas_min_hedge_risk_hs_fv(aapl_call_delta, ibm_call_delta, conf_level)  # start near deltas
print("The optimal stock position under HS and full valuation approach should be:")
print(f"Short {eta_minimise_hs_fv[0]:.4f} Apple shares")
print(f"Short {eta_minimise_hs_fv[1]:.4f} IBM shares")
print(f"The minimised VaR is expected to be: {var_minimise_hs_fv:.4f}")

# Apple stock
if eta_minimise_hs_fv[0] == eta_minimise[0]:
    print("The optimal position for Apple stock and IBM stock under HS and FV is the same as the Analytical approach")
else:
    print(f"The optimal position for Apple stock and IBM stock under HS and FV is different from the Analytical approach by {abs(eta_minimise_hs_fv[0] - eta_minimise[0]):.4f} shares")

# IBM stock
if eta_minimise_hs_fv[1] == eta_minimise[1]:
    print("The optimal position for IBM stock under HS and FV is the same as the Analytical approach")
else:
    print(f"The optimal position for IBM stock under HS and FV is different from the Analytical approach by {abs(eta_minimise_hs_fv[1] - eta_minimise[1]):.4f} shares")

# VaR results
if var_minimise_hs_fv == var_minimise:
    print("The minimised VaR under HS and FV is the same as the Analytical approach")
else:
    print(f"The minimised VaR under HS and FV is different from the Analytical approach by {abs(var_minimise_hs_fv - var_minimise):.4f}")

The optimal stock position under HS and full valuation approach should be:
Short 0.5480 Apple shares
Short 0.5726 IBM shares
The minimised VaR is expected to be: 2.3563
The optimal position for Apple stock and IBM stock under HS and FV is different from the Analytical approach by 0.0047 shares
The optimal position for IBM stock under HS and FV is different from the Analytical approach by 0.0051 shares
The minimised VaR under HS and FV is different from the Analytical approach by 0.0561

Optimal stock position under Monte Carlo simulation method with a full valuation approach assuming that the joint distribution of the log returns follow the $t$ copula

# Function of the loss scenarios under the Monte Carlo simulation with t copula and full valuation approach
def minimise_loss_scenarios_mc_fv_tcop(eta1, eta2, conf_level):

    # portfolio loss scenarios
    minimise_portfolio_value_mc_fv_tcop = aapl_call_price + ibm_call_price - eta1*aapl_price - eta2*ibm_price
    minimise_portfolio_scenarios_mc_fv_tcop = aapl_call_scenarios_mc_fv_tcop + ibm_call_scenarios_mc_fv_tcop - eta1*aapl_price_scenarios_mc_fv_tcop - eta2*ibm_price_scenarios_mc_fv_tcop
    minimise_loss_scenarios_mc_fv_tcop = -(minimise_portfolio_scenarios_mc_fv_tcop - minimise_portfolio_value_mc_fv_tcop)

    # risk measure VaR using standard historical simulation
    var_es_mc = VaR_ES_emp(minimise_loss_scenarios_mc_fv_tcop, conf_level=conf_level, hs_weight_type="equal")
    return var_es_mc[0]           # return VaR only

# Minimise VaR over eta1, eta2 ≥ 0
def etas_min_hedge_risk_mc_fv_tcop(eta1_start=1.0, eta2_start=1.0, conf_level=0.99):        #eta1_start and eta2_start: initial guess
    obj = lambda e: minimise_loss_scenarios_mc_fv_tcop(e[0], e[1], conf_level)
    mc_fv_tcop = optimize.minimize(obj, x0=[eta1_start, eta2_start], bounds=[(0,None),(0,None)], method="L-BFGS-B")
    return mc_fv_tcop.x, mc_fv_tcop.fun

eta_minimise_mc_fv_tcop, var_minimise_mc_fv_tcop = etas_min_hedge_risk_mc_fv_tcop(aapl_call_delta, ibm_call_delta, conf_level)  # start near deltas
print("The optimal stock position under Monte Carlo simulation with t copula and full valuation approach should be:")
print(f'Short {eta_minimise_mc_fv_tcop[0]:.4f} Apple shares')
print(f'Short {eta_minimise_mc_fv_tcop[1]:.4f} IBM shares')
print(f"The minimised VaR is expected to be: {var_minimise_mc_fv_tcop:.4f}")

The optimal stock position under Monte Carlo simulation with t copula and full valuation approach should be:
Short 0.5545 Apple shares
Short 0.5745 IBM shares
The minimised VaR is expected to be: 2.3509

hedging_results = pd.DataFrame({
    'Method': ['Analytical', 'Standard HS', 'MCS (t copula)'],
    'Optimal Apple': [eta_minimise[0], eta_minimise_hs_fv[0], eta_minimise_mc_fv_tcop[0]],
    'Optimal IBM': [eta_minimise[1], eta_minimise_hs_fv[1], eta_minimise_mc_fv_tcop[1]],
    'Minimised VaR': [var_minimise, var_minimise_hs_fv, var_minimise_mc_fv_tcop]
}).to_markdown(index=False)
print(hedging_results)

| Method         |   Optimal Apple |   Optimal IBM |   Minimised VaR |
|:---------------|----------------:|--------------:|----------------:|
| Analytical     |        0.552707 |      0.577751 |         2.30022 |
| Standard HS    |        0.547998 |      0.572612 |         2.35631 |
| MCS (t copula) |        0.55447  |      0.574479 |         2.35087 |

Under the Analytical approach, the optimal $\eta_1$ and $\eta_2$ (numbers of Apple stocks and IBM stocks) are the same as the numbers of the Delta hedging. This approach does not consider the effects of time (theta) and gamma.​ ​

Under the Historical simulation - Full valuation approach, the optimal $\eta_1$ and $\eta_2$ have shifted (down) due to gamma and theta effects plus historical tail dependence.

Under the Monte Carlo - Full valuation approach - $t$ copula joint distribution, the optimal $\eta_1$ and $\eta_2$

We can also comment on the riskiness of the portfolio over the risk management time horizon under the optimal number of shares as even though the linear exposure $\mathbf{b}_t$ is zero, there are still risks introduced over the time horizon $\Delta t$ such as:

• Gamma: large moves in and tail dependence under $t$-copula create losses through the call exposures (second order) $\frac{1}{2} \mathbf{X}_{t+1}^\top \mathbf{B}_{\text{call}, t} \mathbf{X}_{t+1}$ term.
• Theta (time decay): shifts the mean of the loss by $\mathbf{c}_t$ value as the more time passes, the lower the option value becomes.
• Dependence risk: exhibits on upper and lower tail dependence.
• Imperfect hedge: even if all the assumptions of the Black/Scholes holds, in practice, hedges are still imperfect. This results in an uncertain loss (Chiarella, C., He, X.-Z., & Sklibosios Nikitopoulos, C. (2015)).

7. Project Summary and Conclusions

Conclusions on risk measures calculated for Stage 2 and Stage 3

# Display all risk measures in a tabular format
all_risk_measures = pd.DataFrame({
    'Model': ['Analytical', 'Analytical', 'Historical Simulation', 'Historical Simulation', 'Age-Weighted', 'Age-Weighted', 'Monte Carlo Simulation', 'Monte Carlo Simulation','Monte Carlo Simulation', 'Monte Carlo Simulation'],
    'Approach': ['Undiversified','Diversified', 'Delta-Gamma', 'Full Valuation', 'Delta-Gamma','Full Valuation', 'Delta-Gamma (t copula)' ,'Full Valuation (t copula)','Delta-Gamma (Gaussian)' ,'Full Valuation FV (Gaussian)'],
    'VaR': [portfolio_risk_measures[0] ,portfolio_risk_measures_d[0], rm_hs_dg[0], rm_hs_fv[0], rm_hs_aw[0], rm_hs_aw_fv[0], rm_mc_dg_tcop[0], rm_mc_fv_tcop[0], portfolio_risk_measures_mc_dg_Gauss[0] ,rm_mc_fv_Gauss[0]],
    'ES': [portfolio_risk_measures[1], portfolio_risk_measures_d[1], rm_hs_dg[1], rm_hs_fv[1], rm_hs_aw[1], rm_hs_aw_fv[1], rm_mc_dg_tcop[1], rm_mc_fv_tcop[1], portfolio_risk_measures_mc_dg_Gauss[1], rm_mc_fv_Gauss[1]]
})
display(all_risk_measures)

	Model	Approach	VaR	ES
0	Analytical	Undiversified	33.934413	38.375647
1	Analytical	Diversified	23.611937	26.549553
2	Historical Simulation	Delta-Gamma	20.057940	23.018045
3	Historical Simulation	Full Valuation	19.989002	23.156404
4	Age-Weighted	Delta-Gamma	11.994233	13.221086
5	Age-Weighted	Full Valuation	12.052425	13.268244
6	Monte Carlo Simulation	Delta-Gamma (t copula)	22.552221	24.508520
7	Monte Carlo Simulation	Full Valuation (t copula)	22.354283	24.326579
8	Monte Carlo Simulation	Delta-Gamma (Gaussian)	19.065532	20.967462
9	Monte Carlo Simulation	Full Valuation FV (Gaussian)	18.950455	20.804926

The most appropriate method to use is the Monte Carlo (full re-valuation) with a fitted $t$ copula for dependence.

What we have done in Stage 2 and Stage 3 is to estimate the 10-day 99% VaR and ES of the portfolio using different methods and model assumptions, which are:

• Analytical – Assumes linearity of losses and risk factors follow the Normal distribution.
• Delta–Gamma Approximation – Adds curvature; still assumes Normal distribution.
• Historical Simulation (HS) – Replays past market moves with equal weighting.
• Age-weighted HS – Similar to HS, but gives more weight to recent data.
• MC full valuation – $t$ copula – Simulates and fully reprices with fat-tailed dependence.
• MC full valuation – Gaussian copula – Simulates and fully reprice with Gaussian dependence.

The data are telling us that in the portfolio with going long call options and shorting stocks, the loss is non-linear (gamma/theta effects). Moreover, using pseudo-observations has showed that there is presence of positive tail dependence, and that's why MC ($t$ copula) comes out higher than MC (Gaussian), as the $t$ copula captures both upper and lower tails whereas the Gaussian copula does not. This results in higher loss estimates under the $t$ copula, as it accounts for both upper and lower tail dependence - something the Gaussian copula misses.

Moreover, the MC ($t$ copula) framework allows for practice stress testing of parameters such as:

• Degrees of freedom ($\nu$)
• Correlation ($\rho$)
• Volatility
• Maturities

Strengths and trade-offs of using MC ($t$ copula):

Strengths of MC with $t$ copula:

• Handles non-linear portfolios effectively.
• Captures joint tail behavior better than Gaussian models.
• Transparent: losses come from real repricing, not approximations (in the case of Analytical, Delta-Gamma, HS and Age-weighted HS).

Limitations of MC with $t$ copula:

• Requires a parametric fit ($\rho$, $ν$).
• Less efficient in low dimensional setting.
• In the case of doing a full valuation of large portfolios, it is often computationally expensive.
• If the data have asymmetric tails where upper tail dependence is different from lower tail dependence, the $t$ copula can't capture them.
• Simulation error – use enough paths and set a seed.

To sum up, the MC full re-valuation with $t$ copula is the best balance of realism and control as it handles option curvature and the heavy and joint tails observed in the data.

Conclusions of the stress-testing approaches in Stage 4

Strengths:

• Parameters are re-estimated on the GFC (SP1) and COVID-19 (SP2) windows, so today’s portfolio is tested against actual stressed volatility and dependence.
• Monte Carlo simulation using the full valuation approach captures Gamma and Theta sensitivities, instead of relying only on linear loss (delta loss approximation).
• The three methods Analytical, Historical Simulation and Monte Carlo Simulation provide results to cross-checks and compare which assumptions drive VaR/ES.
• When comparing the Gaussian copula and $t$ copula, we can understand the impact of tail dependence ($\lambda >0$ when $\nu$ is finite).

Limitations:

• The results from the testing windows depend on the chosen dates and length and different slices of a crisis can move VaR/ES.
• The Analytical approach assumes linearity of losses and that the risk factor changes follow the Normal distribution, as it seems to over-estimates the VaR (this method results in the higest VaR).
• The Historical Simulation method assume the distribution of losses/risk factor changes is stationery over the selected stressed windows, explaning the under-estimated VaR (this method generates the lowest VaR)
• The MC $t$-copula is symmetric and can not capture upper/lower tail asymmetry.
• The MC Gaussian copula may miss very heavy univariate tails.
• Parameter uncertainty in $\rho$ and $\nu$ between the two data periods.

Comments on the results:

• For the Stressed Period 1 (GFC): Fitted t-copula shows low $\nu$ ~3 and high tail dependence $\lambda$ approx 0.3. The MC t-copula gives the most appropriate VaR compared to MC Gaussian, although the Analytical approach is largest overall ()

• For the Stressed Period 2 (COVID-19): Fitted t-copula has large $\nu$ (~33) and $\lambda$ approx 0 (near-Gaussian dependence). The MC Gaussian seems to give the most appropriate VaR compared to MC $t$ copula and the Analytical approach remains high.

• The Historical Simulation tends to be lowest in both windows.

• After optimal hedging, risk is mainly from gamma/theta and joint extremes; delta exposure is largely neutralized.

Recommendation:

• Use MC full valuation as the primary stressed measure. Under:
  • SP 1 (GFC): the prefered approach is using the $t$ copula as it captures tail dependence.
  • SP 2 (COVID-19): we should take the higher of t-copula and Gaussian as dependence is near-Gaussian, explained by the large $\nu$.
• Report the worst of each stressed period as the final stressed VaR, with Analytical/HS shown as context.
• Note parameter uncertainty especially $\nu$ and $\rho$.

Conclusions of the simple risk minimisation strategy in Stage 5

What the strategy is: We choose short stock positions ($\eta_1$, $\eta_2$ >0) to minimise 10-day risk under the delta (linear) loss model with log‐return risk factors: $$ L^{\Delta}_{t+1} = -\mathbf{c}_t - \mathbf{b}^\top_t X_{i, t+1} $$

$$ \mathbf{c}_t = (\Theta_{1,\text{call}} +\Theta_{2,\text{call}}) \Delta t$$

$$ \mathbf{b}_t= \begin{bmatrix} e^{Z_{1,t}}(\Delta_{1, \text{call}} - \eta_1) \\ e^{Z_{2,t}}(\Delta_{2, \text{call}} - \eta_2) \\ \end{bmatrix} $$

Because VaR/ES under this model are monotone in $\sigma_L = \sqrt{Var(L_{t+1})} = \sqrt{\mathbf{b}^\top_t \Sigma \mathbf{b}_t}$, the minimum is at $\mathbf{b}_t=0$, i.e.

$$ \boxed{\eta_i^\star=\Delta_i \quad \text{(delta–neutral in log price)}} $$

Strengths of the simple risk minimisation strategy:

• The stretegy is simple and fast to implement as it provides a closed-form solution with no simulation needed.
• It is intuitive as going short on the delta removes the first-order exposure to log-return shocks.
• It only requires parameters such as current prices, deltas and a (EWMA) covariance estimate.

Limitations of the simple risk minimisation strategy:

• The method only removes the delta linear risk and leave portfolio subject to residual risks from gamma (convexity, second-order derivative), theta (time decay) and potentially higher order Greeks.
• The delta–normal model assumes risk factors to follow the elliptical/Normal distribution and does not consider the heavy tails or tail dependence (e.g., $t$ copula), which can leave more tail risk than the model suggests.
• Imperfect hedge always exists and as such errors in deltas, vol, and the covariance matrix can underhedge or overhedge (Chiarella, C., He, X.-Z., & Sklibosios Nikitopoulos, C., 2015).
• Other risks are not considered such as liquidity risk (mentioned in Choy, S. K., & Wei, J. (2020) that a negative liquidity risk premium exists in delta-hedge option returns) and default risk (high default risk is related to low average returns on options, Vasquez, A., & Xiao, X. (2024))

How to generalize to more complex portfolios:

• Multi-asset, multiple hedging instruments (delta variance minimisation).
  Let $d$ be the vector of portfolio log-deltas to risk factors, and let $A$ be the matrix of factor loadings of the hedging instruments (each column is a hedge’s delta to each factor). Minimise $\eta^\star=(A^\top\Sigma A)^{-1}A^\top\Sigma d$ by using the Least Squares method (Jorion, 2007).

• Include curvature (delta–gamma hedging): Choose hedge options so that both delta and gamma are near zero (solve a small linear system), or minimise the Variance of ($\mathbf{b}^\top X+\tfrac12 X^\top {\mathbf{B}} X$) using scenarios (Dowd, 2005).

• Practical constraints: Add transaction-cost penalties and consider other risks (e.g. liquidity risk) to keep the solution tradeable.

REFERENCES

• Alòs, E., Merino, R., & Rolloos, F. (2023). Introduction to Financial Derivatives with Python (First edition, Vol. 1). CRC Press. https://doi.org/10.1201/9781003266730.

• Chiarella, C., He, X.-Z., & Sklibosios Nikitopoulos, C. (2015). Derivative Security Pricing : Techniques, Methods and Applications (1st ed. 2015.). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-662-45906-5.

• Choy, S. K., & Wei, J. (2020). Liquidity risk and expected option returns. Journal of Banking & Finance, 111, Article 105700. https://doi.org/10.1016/j.jbankfin.2019.105700.

• Dowd, K. (2005). Measuring market risk (2nd ed.). Wiley.

• Hull, J. C. (2022). Options, futures, and other derivatives (Eleventh edition, Global Edition.). Pearson.

• Jorion, P. (2007). Value at Risk: The new benchmark for managing financial risk (3rd ed.). McGraw–Hill.

• OpenAI. (2025). ChatGPT (October 25 version) response to “Explain differences between delta–gamma and full re-valuation VaR and ES” [Large language model].https://chat.openai.com/.

• Vasquez, A., & Xiao, X. (2024). Default Risk and Option Returns. Management Science, 70(4), 2144–2167. https://doi.org/10.1287/mnsc.2023.4796.

• Villa, C., & Rubio, F. J. (2018). Objective priors for the number of degrees of freedom of a multivariate t distribution and the t-copula. Computational Statistics & Data Analysis, 124, 197–219. https://doi.org/10.1016/j.csda.2018.03.010.