Optimal Portfolio Allocation

An investment universe of the following risky assets with a dependence structure (correlation) is given:

\begin{equation*} \begin{matrix} \textbf{Asset} & \boldsymbol\mu & \boldsymbol\sigma & \boldsymbol\omega\\ A & 0.02 & 0.05 & \omega_1\\ B & 0.07 & 0.12 & \omega_2\\ C & 0.15 & 0.17 & \omega_3\\ D & 0.20 & 0.25 & \omega_4 \end{matrix} \qquad \qquad R = \begin{pmatrix} 1 & 0.3 & 0.3 & 0.3\\ 0.3 & 1 & 0.6 & 0.6\\ 0.3 & 0.6 & 1 & 0.6\\ 0.3 & 0.6 & 0.6 & 1 \end{pmatrix} \end{equation*}

Question 1. Consider minimum variance portfolio with a target return m.

\begin{equation*} \underset{\boldsymbol\omega}{\arg\min} \frac{1}{2}\omega'\Sigma\omega \qquad s.t. \quad \omega'1 = 1 \quad \mu_\Pi = \omega' \mu = m \end{equation*}

• Formulate the Lagrangian and give its partial derivatives.
• Write down the analytical solution for optimal allocations $ \omega^* $ (derivation not required).
• Inverse optimisation: generate above 700 random allocation sets (vectors) 4 × 1, those will not be optimal allocations.
  • Standardise each set to satisfy $\omega'1 = 1$, in fact you can generate 3 allocations and compute the 4th.
  • For each vector of allocations compute $\mu_\Pi = \omega'\mu$ and $\sigma_\Pi = \sqrt{\omega'\Sigma\omega}$.
  • Plot the cloud of points of $\mu_\Pi$ vertically on $\sigma_\Pi$ horizontally. Explain this plot.

Hint: Since we treat this as an inverse optimisation, there is no computation of $\omega^*$ from the ready formula.

Answers and code

Form the Langrange function with two multipliers $\lambda$ and $\gamma$ for a given vector of weights $w$ \begin{equation*} L(w, \lambda, \gamma) = \frac{1}{2}w'\Sigma w + \lambda(m - w'\mu) + \gamma(1 - w'1) \end{equation*}

where $(m - w'\mu)$ is the return constraint and $(1 - w'1)$ is the budget constraint.

First order condition

The optimal portfolio is the point where the gradient of $L(w, \lambda, \gamma) = 0$.

Therefore $\frac{\partial L}{\partial w} = 0$.

\begin{equation*} \frac{\partial L}{\partial w} = \frac{1}{2} 2 \Sigma w + \lambda (-\mu) + \gamma(-1) = \Sigma w - \lambda \mu - \gamma 1 = 0 \end{equation*}

Second order condition

Given this is a minimization problem, we seek the Hessian to be greater than 1

\begin{equation*} \frac{\partial^2 L}{\partial w^2} = \Sigma \end{equation*}

As the covariance matrix is positive, the 2nd order condition is satisfied.

Therefore the optimal allocation $w^*$ can be derived

\begin{equation*} \Sigma w - \lambda \mu - \gamma 1 = 0 \end{equation*}

\begin{equation*} w^* = \Sigma^{-1} (\lambda \mu + \gamma 1) \end{equation*}

Substituting the above optimal weight $w^*$ into the constraint

\begin{equation*} \mu'w = \mu' \Sigma^{-1} (\lambda \mu + \gamma 1) = \mu' \Sigma^{-1} \lambda \mu + \mu' \Sigma^{-1}\gamma 1 = m \end{equation*} and \begin{equation*} 1'w = 1' \Sigma^{-1} (\lambda \mu + \gamma 1) = 1' \Sigma^{-1} \lambda \mu + 1' \Sigma^{-1}\gamma 1 = 1 \end{equation*}

Setting scalars for ease of presentation and computation \begin{equation*} A = 1\Sigma^{-1} 1' \end{equation*}

\begin{equation*} B = \mu' \Sigma^{-1} 1 = 1' \Sigma^{-1} \mu \end{equation*}

\begin{equation*} C = \mu' \Sigma^{-1} \mu \end{equation*}

The constraint functions then become

\begin{equation*} C \lambda + B \gamma = m \end{equation*}

\begin{equation*} B \lambda + A \gamma = 1 \end{equation*}

This is a system of 2 equations with 2 unknowns

\begin{equation*} \lambda = \frac{m - B \gamma}{C} \end{equation*}

therefore

\begin{equation*} \frac{Bm - B^2 \gamma}{C} + A \gamma = 1 \end{equation*}

\begin{equation*} Bm - B^2 \gamma + A C \gamma = C \end{equation*}

\begin{equation*} A C \gamma - B^2 \gamma = C - Bm \end{equation*}

\begin{equation*} \gamma = \frac{C - Bm}{AC - B^2} \end{equation*}

in the same manner

\begin{equation*} \lambda = \frac{Am - B}{AC - B^2} \end{equation*}

Finally, substituting back into optimal weight $w^*$

\begin{equation*} w^* = \Sigma^{-1} (\lambda \mu + \gamma 1) = \frac{1}{AC - B^2} \Sigma^{-1}[(A \mu - B1)m + (C1 - B \mu)] \end{equation*}

# Import the relevant packages
# Import numpy
import plotly.express as px
import numpy as np
from numpy import *
from numpy.linalg import multi_dot

# Import pandas
import pandas as pd

# Import plotly express
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 650, 450
import plotly.io as pio
pio.renderers.default='notebook'

# Set the matrices for mu and sigma
mu = np.array([0.02, 0.07, 0.15, 0.2]).reshape((4, 1))
sigma = np.array([0.05, 0.12, 0.17, 0.25]).reshape((4, 1))
R = np.array([[1, 0.3, 0.3, 0.3], [0.3, 1, 0.6, 0.6], [
             0.3, 0.6, 1, 0.6], [0.3, 0.6, 0.6, 1]]).reshape(4, 4)

# Compute the covariance matrix usings SRS where S = the diagonal standard deviation matrix
S = np.diag(sigma.flatten())
covariance_matrix = multi_dot([S, R, S])

# Set values for N and initiate an empty list for the
N = 10000
weights_arr = []

# Generate N random weights for the 4 assets
for i in range(N):
    weights = random.random(3)
    weights = np.append(weights, 1-sum(weights))
    # Standardise these weights
    weights /= sum(weights)
    weights_arr.append(weights)

# Add them to a dataframe
df = pd.DataFrame(weights_arr, columns=['w_1', 'w_2', 'w_3', 'w_4'])

df.head()

	w_1	w_2	w_3	w_4
0	0.261793	0.783503	0.306526	-0.351821
1	0.741104	0.086845	0.053350	0.118701
2	0.265768	0.614823	0.080235	0.039174
3	0.429468	0.579392	0.338292	-0.347152
4	0.159600	0.401197	0.992521	-0.553318

df.tail()

	w_1	w_2	w_3	w_4
9995	0.360908	0.620995	0.416362	-0.398265
9996	0.982305	0.434214	0.582299	-0.998818
9997	0.247895	0.809135	0.752616	-0.809646
9998	0.184516	0.410816	0.341352	0.063315
9999	0.539174	0.479185	0.677368	-0.695727

# loop through weights, extract matrix, calculate portfolio mean and sd, append back to df
for index, row in df.iterrows():
    # convert to a matrix of weights
    w = np.array(row).reshape(1, 4)
    # Compute mean for the portfolio, add it to the dataframe
    mu_pi = multi_dot([w, mu])
    df.loc[index, 'mu_pi'] = mu_pi.flatten()
    # Compute the sd for the portfolio, add it to the dataframe
    sigma_pi = np.sqrt(multi_dot([w, covariance_matrix, w.T]))
    df.loc[index, 'sigma_pi'] = sigma_pi.flatten()

df.head()

	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
0	0.261793	0.783503	0.306526	-0.351821	0.035696	0.101859
1	0.741104	0.086845	0.053350	0.118701	0.052644	0.065918
2	0.265768	0.614823	0.080235	0.039174	0.068223	0.094403
3	0.429468	0.579392	0.338292	-0.347152	0.030460	0.090670
4	0.159600	0.401197	0.992521	-0.553318	0.069490	0.155274

df.tail()

	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
9995	0.360908	0.620995	0.416362	-0.398265	0.033489	0.101324
9996	0.982305	0.434214	0.582299	-0.998818	-0.062378	0.187849
9997	0.247895	0.809135	0.752616	-0.809646	0.012561	0.165337
9998	0.184516	0.410816	0.341352	0.063315	0.096313	0.110845
9999	0.539174	0.479185	0.677368	-0.695727	0.006786	0.139408

# Plotting the results
fig = px.scatter(df, x='sigma_pi', y='mu_pi', color='mu_pi', color_continuous_scale='Viridis',
                 labels={'sigma_pi': 'Expected Volatility',
                         'mu_pi': 'Expected Return'},
                 title='Efficient Frontier highlighting maximum and minimum return')
idx_max_mu = df['mu_pi'].idxmax()
idx_min_sigma = df['sigma_pi'].idxmin()

fig.add_scatter(x=[df['sigma_pi'].loc[idx_max_mu]], y=[df['mu_pi'].loc[idx_max_mu]], marker=dict(color='green', size=10, symbol='cross'),
                name='Maximum Return').update(layout_showlegend=False)
fig.add_scatter(x=[df['sigma_pi'].loc[idx_min_sigma]], y=[df['mu_pi'].loc[idx_min_sigma]], marker=dict(color='red', size=10, symbol='cross'),
                name='Minimum Volatility').update(layout_showlegend=False)
fig.show()

The above plot builds up an example through simulated data showing the efficient frontier of a portfolio of 4 risky assets. The maximum return and minimum variance are marked with green adn red crosses respectively.

Optimal Portfolio Allocation

An investment universe of the following risky assets with a dependence structure (correlation) is given:

\begin{equation*} \begin{matrix} \textbf{Asset} & \boldsymbol\mu & \boldsymbol\sigma & \boldsymbol\omega\\ A & 0.02 & 0.05 & \omega_1\\ B & 0.07 & 0.12 & \omega_2\\ C & 0.15 & 0.17 & \omega_3\\ D & 0.20 & 0.25 & \omega_4 \end{matrix} \qquad \qquad R = \begin{pmatrix} 1 & 0.3 & 0.3 & 0.3\\ 0.3 & 1 & 0.6 & 0.6\\ 0.3 & 0.6 & 1 & 0.6\\ 0.3 & 0.6 & 0.6 & 1 \end{pmatrix} \end{equation*}

Question 2. Consider optimisation for a tangency portfolio (maximum Sharpe Ratio).

• Formulate optimisation expression.
• Formulate Lagrangian function and give its partial derivatives only.
• For the range of tangency portfolios given by $r_f = 50bps, 100bps, 150bps, 175bps$ optimal compute allocations (ready formula) and $\sigma_\Pi$ . Present results in a table.
  • Plot the efficient frontier in the presence of a risk-free asset for $r_f = 100bps, 175bps$

Answers and code

The tangency portfolio is invested wholly in risky assets and that maximises the sharpe ratio.

\begin{equation*} \underset{\boldsymbol\omega}{\arg\max} \quad \frac{\mu - r_f}{\sigma} = \frac{\omega' \mu - r_f}{(\omega' \Sigma \omega)^{\frac{1}{2}}} \qquad s.t. \quad \omega'1 = 1 \end{equation*}

Form the Langrange function with one multiplier $\lambda$ for a given vector of weights $w$ \begin{equation*} L(w, \lambda) = (w' \mu - r_f)(w' \Sigma w)^{-\frac{1}{2}} + \lambda(1 - w'1) \end{equation*}

Derivatives

Using a combination of the chain and product rules to obtain the first order derivative

$$ \frac{\partial L}{\partial w} = \mu (w' \Sigma w)^{-\frac{1}{2}} - \frac{1}{2}(w' \mu - r_f)(w' \Sigma w)^{-\frac{3}{2}}2 \Sigma w + \lambda(-1) $$

simplifying

$$ \frac{\partial L}{\partial w} = \mu (w' \Sigma w)^{-\frac{1}{2}} - (w' \mu - r_f)(w' \Sigma w)^{-\frac{3}{2}} \Sigma w - \lambda $$

and the second order derivative

$$ \frac{\partial^2 L}{\partial w^2} = -\frac{1}{2}\mu (w' \Sigma w)^{-\frac{5}{2}} 2\Sigma w - -\frac{3}{2} (w'\Sigma w)^{-\frac{5}{2}} 2 \Sigma w \Sigma w + (w' \Sigma w)^{-\frac{3}{2}} \Sigma $$

simplifying

$$ \frac{\partial^2 L}{\partial w^2} = -\mu (w' \Sigma w)^{-\frac{5}{2}} \Sigma w + 3 (w'\Sigma w)^{-\frac{5}{2}} \Sigma w \Sigma w - (w' \Sigma w)^{-\frac{3}{2}} \Sigma $$

The ready formula to calculate the optimal weights of the tangency portfolio is given by

$$ w_t = \frac{\Sigma^{-1} (\mu - r_f1)}{B - Ar_f} $$

where A and B are scalars as previously defined in question 1

\begin{equation*} A = 1' \Sigma^{-1} 1 \end{equation*}

\begin{equation*} B = \mu' \Sigma^{-1} 1 = 1' \Sigma^{-1} \mu \end{equation*}

# Import the relevant packages
# Import numpy
import plotly.io as pio
import plotly.express as px
import numpy as np
from numpy import *
from numpy.linalg import multi_dot, inv

# Import pandas
import pandas as pd

# Import plotly express
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 650, 450
pio.renderers.default = 'notebook'

# Set the risk free rates
r = [0.005, 0.01, 0.015, 0.0175]

# Set the matrices for mu and sigma and 1
mu = np.array([0.02, 0.07, 0.15, 0.2]).reshape((4, 1))
sigma = np.array([0.05, 0.12, 0.17, 0.25]).reshape((4, 1))
R = np.array([[1, 0.3, 0.3, 0.3], [0.3, 1, 0.6, 0.6], [
             0.3, 0.6, 1, 0.6], [0.3, 0.6, 0.6, 1]]).reshape(4, 4)
one = np.ones((4, 1))

# Compute the covariance matrix usings SRS where S = the diagonal standard deviation matrix
S = np.diag(sigma.flatten())
covariance_matrix = multi_dot([S, R, S])

# Computer the inverse covariance matrix
inverse_covariance_matrix = inv(covariance_matrix)
# Compute A and B
A = multi_dot((one.T, inverse_covariance_matrix, one)).flatten()[0]
B = multi_dot((one.T, inverse_covariance_matrix, mu)).flatten()[0]

# For each risk free rate, compute the optimal weights and add them to a dataframe
arr = []
for rate in r:
    risk_premium = mu - multi_dot((rate, one))
    denominator = (B - A*rate)
    weights = multi_dot([inverse_covariance_matrix, risk_premium]) / denominator
    # Add the risk free rate and the sigma_pi
    sigma_pi = np.sqrt(multi_dot((weights.T, covariance_matrix, weights)))
    mu_pi = multi_dot([weights.T, mu])
    row = np.insert(weights, 0, rate)
    row = np.append(row, mu_pi)
    row = np.append(row, sigma_pi)
    arr.append(row)

df_tangency = pd.DataFrame(
    arr, columns=['rf', 'w_1', 'w_2', 'w_3', 'w_4', 'mu_pi', 'sigma_pi'])
df_tangency

	rf	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
0	0.0050	0.016835	-0.229367	0.814340	0.398192	0.186070	0.196511
1	0.0100	-0.745937	-0.510569	1.490249	0.766257	0.326130	0.350665
2	0.0150	-8.644854	-3.422571	8.489651	4.577774	1.776525	1.972392
3	0.0175	8.103502	2.751851	-6.351431	-3.503922	-1.298799	1.473515

# Simulate random portfolios.
# Set values for N and initiate an empty list for the
N = 10000
weights_arr = []

# Generate N random weights for the 4 assets
for i in range(N):
    weights = random.uniform(-3, 3, 4)
    # weights = np.append(weights, 1-sum(weights))
    weights_arr.append(weights)

# Add them to a dataframe
df = pd.DataFrame(weights_arr, columns=['w_1', 'w_2', 'w_3', 'w_4'])

# loop through weights, extract matrix, calculate portfolio mean and sd, append back to df
for index, row in df.iterrows():
    # convert to a matrix of weights
    w = np.array(row).reshape(1, 4)
    # Compute mean for the portfolio, add it to the dataframe
    mu_pi = multi_dot([w, mu])
    df.loc[index, 'mu_pi'] = mu_pi.flatten()
    # Compute the sd for the portfolio, add it to the dataframe
    sigma_pi = np.sqrt(multi_dot([w, covariance_matrix, w.T]))
    df.loc[index, 'sigma_pi'] = sigma_pi.flatten()

df.head()

	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
0	0.249282	-2.659690	-0.830326	-1.052972	-0.516336	0.620728
1	-1.723734	-1.308708	-1.522946	-0.996548	-0.553836	0.607229
2	-1.115281	1.781101	1.480385	1.696261	0.663681	0.750717
3	1.079599	2.928616	-2.902509	0.406833	-0.127415	0.391805
4	1.279221	-0.105250	2.303874	0.235610	0.410920	0.445845

# Plotting the results
fig = px.scatter(df, x='sigma_pi', y='mu_pi',
                 labels={'sigma_pi': 'Expected Volatility',
                         'mu_pi': 'Expected Return'},
                 title='Efficient Frontier highlighting maximum and minimum return', opacity=0.5)

sigma_100 = df_tangency['sigma_pi'].iloc[1]
mu_100 = df_tangency['mu_pi'].iloc[1]
sigma_175 = df_tangency['sigma_pi'].iloc[3]
mu_175 = df_tangency['mu_pi'].iloc[3]

fig.add_shape(type="line",
              x0=0, y0=0.01, x1=sigma_100, y1=mu_100,
              line=dict(color="#32CD32", width=2), opacity=1
              )
fig.add_shape(type="line",
              x0=0, y0=0.0175, x1=sigma_175, y1=mu_175,
              line=dict(color="#0096FF", width=2), opacity=1,
              )
fig.add_scatter(x=sigma_100.flatten(), y=mu_100.flatten(), marker_symbol='star', mode='markers', marker_color='#32CD32',
                marker_size=10, name='TP 100bps', hovertemplate='Expected Return: %{y:.2f}\nVol: %{x:.2f}')
fig.add_scatter(x=sigma_175.flatten(), y=mu_175.flatten(), marker_symbol='star', mode='markers', marker_color='#0096FF',
                marker_size=10, name='TP 100bps', hovertemplate='Expected Return: %{y:.2f}\nVol: %{x:.2f}')
fig.show()

Products and Market Risk

Question 3. Implement the multi-step binomial method as described in Binomial Method lecture with the following variables and parameters: stock $S = 100$, interest rate $r = 0.05$ (continuously compounded) for a call option with strike $E = 100$, and maturity $T = 1$.

• Use any suitable parametrisation for up and down moves $uS$, $vS$.
• Compute the option value for a range of volatilities $[0.05, . . . , 0.80]$ and plot the result. Set trees to have a minimum four time steps.
• Now, compute and plot the value of one option, $\sigma_{imp} = 0.2$ as you increase the number of time steps $NTS = 4, 5, . . . , 50$.

Hint: This is a computational problem, best coded in Python. Plots must be presented.

Answers and code

The parameters for risk neutral binomial option pricing are defined below:

$$ u = 1 + \sigma \sqrt{\delta t} $$

$$ v = 1 - \sigma \sqrt{\delta t} $$

$$ p' = \frac{1}{2} + \frac{r \sqrt{\delta t}}{2 \sigma} $$

On each step of the tree of $N$ steps $S = u^{n_u}v^{n_v}S$ where $n_u$ is the number of upward steps and $n_v$ the number of downwards steps where $n_u + n_v = N$ at the end of the tree.

The payoff at each time step is

$$Call = \max [S - E, 0]$$ $$Put = \max [E - S, 0]$$

The value of the option at each time step is

$$ V = \frac{1}{e^{r \delta t}} p'V^{+} + (1 - p')V^{-} $$

where $V^{+}$ is the calculated option one step up in the future and $V^{-}$ is the calculated option one step down in the future.

Finally, delta hedging is calculated at each step to be

$$ \Delta = \frac{V^{+} - V^{-}}{(u - v)S} $$

# Import module
import numpy as np
import pandas as pd
# Import plotly express
import plotly.express as px
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 800, 600
import plotly.io as pio
pio.renderers.default='notebook'

# Credit to Kannan Singaravelu for the below function from the Binomial Models python lab
def binomial_option(spot, strike, rate, sigma, time, steps, output=0, is_call=True):
    """
    binomial_option(spot, strike, rate, sigma, time, steps, output=0)

    Function to calculate binomial option pricing for european call option

    Params
    ------
    spot       -int or float    - spot price
    strike     -int or float    - strike price
    rate       -float           - interest rate
    time       -int or float    - expiration time
    steps      -int             - number of time steps
    output     -int             - [0: price, 1: payoff, 2: option value, 3: option delta]

    Returns
    --------
    out: ndarray
    An array object of price, payoff, option value and delta as specified by the output flag

    """

    # define parameters
    ts = time / steps
    u = 1 + sigma*np.sqrt(ts)
    v = 1 - sigma*np.sqrt(ts)
    p = 0.5 + rate * np.sqrt(ts) / (2*sigma)
    df = 1/(np.exp(rate*ts))
    call_put_flag = 1 if is_call else -1

    # initialize the arrays
    px = np.zeros((steps+1, steps+1))
    cp = np.zeros((steps+1, steps+1))
    V = np.zeros((steps+1, steps+1))
    d = np.zeros((steps+1, steps+1))

    # binomial loop : forward loop
    for j in range(steps+1):
        for i in range(j+1):
            px[i, j] = spot * np.power(v, i) * np.power(u, j-i)
            cp[i, j] = np.maximum((px[i, j] - strike)*call_put_flag, 0)

    # reverse loop
    for j in range(steps+1, 0, -1):
        for i in range(j):
            if (j == steps+1):
                V[i, j-1] = cp[i, j-1]
                d[i, j-1] = 0
            else:
                V[i, j-1] = df*(p*V[i, j]+(1-p)*V[i+1, j])
                d[i, j-1] = (V[i, j]-V[i+1, j])/(px[i, j]-px[i+1, j])

    results = np.around(px, 2), np.around(cp, 2), np.around(V, 2), np.around(d, 4)

    return results[output]

spot = 100
strike = 100
rate = 0.05
T = 1
steps = 4
sigma_arr = [0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8]
call_price_arr = []
put_price_arr = []

# loop through sigmas and calculate the option price at T = 0
for sigma in sigma_arr:
    call_price_arr.append(binomial_option(spot, strike, rate, sigma, T, steps, 2, True)[0][0])
    put_price_arr.append((binomial_option(spot, strike, rate, sigma, T, steps, 2, False)[0][0]))
    
data = {'Vol %': sigma_arr, 'Call Price': call_price_arr, 'Put Price': put_price_arr}
df = pd.DataFrame(data)
df['Vol %'] = df['Vol %']*100
df

	Vol %	Call Price	Put Price
0	5.0	5.13	0.28
1	10.0	6.60	1.76
2	15.0	8.39	3.54
3	20.0	10.28	5.44
4	25.0	12.24	7.40
5	30.0	14.25	9.40
6	35.0	16.29	11.44
7	40.0	18.36	13.51
8	45.0	20.46	15.61
9	50.0	22.57	17.73
10	55.0	24.71	19.86
11	60.0	26.85	22.00
12	65.0	29.00	24.15
13	70.0	31.15	26.30
14	75.0	33.30	28.45
15	80.0	35.44	30.59

# Plot the results
fig = px.line(df, x='Vol %', y=['Call Price', 'Put Price'],
                 title='Binomial option value with increasing sigma', labels={'x': 'Vol', 'y': 'Price'}, markers=True)
fig.update_layout(xaxis_title='Volatility %',
                  yaxis_title='Option Value', 
                  legend_title_text='Option type')
fig.show()

The above clearly demonstrates that an option's value increases with increasing volatility.

# Now calculate the option price with increasing time steps
spot = 100
strike = 100
rate = 0.05
T = 1
# Initialise array between 4 and 50
steps_arr = [n for n in range(4,51)]
sigma = 0.2
call_price_arr = []
put_price_arr = []

# loop through sigmas and calculate the option price at T = 0
for N in steps_arr:
    call_price_arr.append(binomial_option(spot, strike, rate, sigma, T, N, 2, True)[0][0])
    put_price_arr.append(binomial_option(spot, strike, rate, sigma, T, N, 2, False)[0][0])
    

data = {'Steps': steps_arr, 'Call Price': call_price_arr, 'Put Price': put_price_arr}
df = pd.DataFrame(data)
df.head()

	Steps	Call Price	Put Price
0	4	10.28	5.44
1	5	10.73	5.87
2	6	10.38	5.52
3	7	10.64	5.78
4	8	10.42	5.56

df.tail()

	Steps	Call Price	Put Price
42	46	10.48	5.60
43	47	10.45	5.57
44	48	10.48	5.60
45	49	10.45	5.57
46	50	10.48	5.60

# Plot the results
fig = px.line(df, x='Steps', y='Call Price',
                 title='Binomial call option value with increasing time steps', labels={'x': 'Steps', 'y': 'Price'}, markers=True)
fig.update_layout(xaxis_title="Time steps",
                  yaxis_title="Option Value")
fig.show()

# Plot the results
fig = px.line(df, x='Steps', y='Put Price',
              title='Binomial put option value with increasing time steps', labels={'x': 'Steps', 'y': 'Price'}, markers=True)
fig.update_layout(xaxis_title="Time steps",
                  yaxis_title="Option Value")
fig.update_traces(line_color='blue')
fig.show()

This plot shows that the option value stabilises at between 22 and 26 time steps (slightly less for a put) and then shows more variability at the extreme ends of the range (with more variability in the price at less time steps).

Products and Market Risk

Question 4. Use the ready formula for Expected Shortfall in order to compute the standardised value of Expected Shortfall for $N (0, 1)$.

• Compute for the following range of percentiles $[99.95; 99.75; 99.5; 99.25; 99; 98.5; 98; 97.5]$ Provide a table.
• The formula to use, and $1 − c$ refers to $1 − 99.95$ and so on,

$$ ES_c(X) = \mu − \sigma \frac{\phi (\Phi^{-1}(1 - c))}{1−c} $$

Hint: For the derivation and explanation of the formula please refer to VaR and ES lecture solutions.

Answers and code

# Import modules
from scipy.stats import norm
import pandas as pd

# Set the percentiles, mu and sd
percentiles = [0.9995, 0.9975, 0.995, 0.9925, 0.99, 0.985, 0.98, 0.975]
mu = 0
sd = 1

# Define a function to compute the expected shortfall

def expected_shortfall(mu, sd, percentile):
    # SciPy ppf function is Percent point function (inverse of cdf — percentiles).
    big_phi_inverse = norm.ppf(1-percentile)
    # Take the pdf of big phi inverse
    phi = norm.pdf(big_phi_inverse)
    return mu - (sd * (phi/(1-percentile)))

# Compute the ES for each percentile
es_arr = []

for percentile in percentiles:
    es = expected_shortfall(mu, sd, percentile)
    es_arr.append([percentile, es])

df = pd.DataFrame(es_arr, columns=['Percentile', 'Expected shortfall'])

df

	Percentile	Expected shortfall
0	0.9995	-3.554381
1	0.9975	-3.104357
2	0.9950	-2.891949
3	0.9925	-2.761240
4	0.9900	-2.665214
5	0.9850	-2.524695
6	0.9800	-2.420907
7	0.9750	-2.337803

Products and Market Risk

Question 5. For this questions, use S&P500 index data provided in order to implement backtesting for $99%/10day$ Value at Risk and report the following:

(a) The count and percentage of VaR breaches.

(b) The count and percentage of consecutive VaR breaches. Example: 1, 1, 1 indicates two consecutive occurrences.

(c) Provide a plot which identifies the breaches.

$$ VaR_{10D,t} = Factor × \sigma_t × \sqrt{10} $$

• Compute the rolling standard deviation $\sigma_t$ from 21 daily returns.
• Timescale of the standard deviation is ‘daily’ regardless of how many returns are in the sample. We project from 1-day to 10-day using the additivity of variance $\sigma_{10D} = \sqrt{\sigma_t × 10}$.
• VaR is fixed at time t and compared to the return realised from $t$ to $t + 10$. A breach occurs when that forward realised 10-day return $ln(\frac{S_{t+10}}{S_t})$ is below the $VaR_t$ quantity.

$r_{10D,t+10} < VaR_{10D,t} \quad $ means breach, given both numbers are negative.

Answers and code

# Import modules
import pandas as pd
import numpy as np
from scipy.stats import norm
import plotly.express as px

# Import plotly express
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 800, 600
import plotly.io as pio
pio.renderers.default='notebook'

# Set probability value
c = 0.01
# Compute VaR factor
factor = norm.ppf(1 - c)
# import the data
df = pd.read_csv('Data_SP500.csv', index_col='Date')
# Compute the daily returns
df['Log Daily Return'] = np.log1p(df.pct_change(1))
# Compute the 10 day return
df['Log 10D Return'] = np.log1p(df['SP500'].pct_change(10))
# Compute rolling 21 day standard deviation for 1D and 10D
df['21D Rolling 1D SD'] = df['Log Daily Return'].rolling(21).std()
df['21D Rolling 10D SD'] = df['21D Rolling 1D SD'] * np.sqrt(10)
# Compute VaR for each row
df['VaR (Return)'] = -factor * df['21D Rolling 10D SD']
# Add the t-10 VaR to the row for comparison
df['VaR to Compare'] = df['VaR (Return)'].shift(10)
# Compute the breaches, 1 for yes and 0 for no
df['Breach'] = df.apply(lambda row: 1 if (row['Log 10D Return'] < row['VaR to Compare']) else 0, axis=1)
# Drop rows that have NaN values
df.dropna(inplace=True)
# Compute consecutive breaches
df['Consecutive breach'] = df['Breach'] + df['Breach'].shift(1)
df['Consecutive breach'] = df['Consecutive breach'].apply(lambda row: 1 if (row > 1) else 0)

df.head()

	SP500	Log Daily Return	Log 10D Return	21D Rolling 1D SD	21D Rolling 10D SD	VaR (Return)	VaR to Compare	Breach	Consecutive breach
Date
07/03/2013	1544.260010	0.001815	0.027468	0.007087	0.022410	-0.052134	-0.043908	0	0
08/03/2013	1551.180054	0.004471	0.023205	0.007125	0.022532	-0.052416	-0.045926	0	0
11/03/2013	1556.219971	0.003244	0.044928	0.007103	0.022462	-0.052255	-0.055272	0	0
12/03/2013	1552.479980	-0.002406	0.036431	0.007083	0.022399	-0.052107	-0.055465	0	0
13/03/2013	1554.520020	0.001313	0.025098	0.007073	0.022366	-0.052031	-0.059057	0	0

# Count the number and % of breaches
breaches = df['Breach'].sum()
percent_breaches = breaches / df['Breach'].count()
consecutive_breaches = df['Consecutive breach'].sum()
percent_consecutive_breaches = consecutive_breaches / df['Consecutive breach'].count()

print('Number of breaches: ', breaches)
print('% of breaches: ', f'{round(percent_breaches*100,3)}%')
print('Number of consecutive breaches: ', consecutive_breaches)
print('% of consecutive breaches: ', f'{round(percent_consecutive_breaches*100,3)}%')

Number of breaches:  25
% of breaches:  2.051%
Number of consecutive breaches:  14
% of consecutive breaches:  1.148%

fig = px.bar(df, x=df.index, y='Log 10D Return', title='Backtesting VaR')
fig.update_traces(marker_color='blue',
                  marker_line_color='blue',
                  selector=dict(type="bar"))
fig.add_scatter(x=df.index, y=df['VaR to Compare'], name='VaR - 10D', mode='lines', opacity=1, marker_color='red')
fig.add_scatter(x=df[df['Breach'] > 0].index, y=df[df['Breach'] > 0]['Log 10D Return'], name='Breach', mode='markers', opacity=1)
fig.show()

Products and Market Risk

Question 6. Re-implement backtesting using assumptions above (as necessary) and the EWMA variance forecast equation below. This is also known as RiskMetrics approach.

The tutor will not reconfirm how to compute this model, but you can use the variance of your computed log-returns for the entire dataset to initialise the scheme.

$$ \sigma_{t+1 | t}^2 = \lambda \sigma_{t | t−1}^2 + (1 − \lambda) r^2_t $$

with $\lambda = 0.72$ value set to minimise out of sample forecasting error, and $r_t$ refers to a return.

Provide the same deliverables (a), (b), and (c) as in the previous Question.

Answers and code

# Import modules
import pandas as pd
import numpy as np
from scipy.stats import norm
import plotly.express as px

# Import plotly express
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 800, 600
import plotly.io as pio
pio.renderers.default='notebook'

# Set probability value, assumed to be the same as the prior question
c = 0.01
# Set lambda
l = 0.72
# Compute VaR factor
factor = norm.ppf(1 - c)
# import the data
df = pd.read_csv('Data_SP500.csv', index_col='Date')
# Compute the daily returns
df['Log Daily Return'] = np.log1p(df.pct_change(1))
# Compute the 10 day return
df['Log 10D Return'] = np.log1p(df['SP500'].pct_change(10))
# Compute initial variance
initial_var = df['Log Daily Return'].var()
df['EWMA var'] = initial_var

# Compute EWMA variance
for index, row in df.iterrows():
    date_range = list(df.index)
    # Skip the first iteration as this is the initial variance
    if index == date_range[0]:
        continue
    # Find the index of the prior day
    prior_day_idx = date_range.index(index) - 1
    # Find the prior date
    prior_day = date_range[prior_day_idx]
    # This represents variance at t-1
    prior_var = df.loc[prior_day]['EWMA var']
    r_squared = row[1]**2
    row['EWMA var'] = l*prior_var + (1-l)*r_squared

# Take square root to compute standard deviations
df['1D SD'] = np.sqrt(df['EWMA var'])
df['10D SD'] = df['1D SD'] * np.sqrt(10)

# Compute VaR for each row
df['VaR (Return)'] = -factor * df['10D SD']
# Add the t-10 VaR to the row for comparison
df['VaR to Compare'] = df['VaR (Return)'].shift(10)
# Compute the breaches, 1 for yes and 0 for no
df['Breach'] = df.apply(lambda row: 1 if (
    row['Log 10D Return'] < row['VaR to Compare']) else 0, axis=1)

# Drop rows that have NaN values
df.dropna(inplace=True)

# Compute consecutive breaches
df['Consecutive breach'] = df['Breach'] + df['Breach'].shift(1)
df['Consecutive breach'] = df['Consecutive breach'].apply(
    lambda row: 1 if (row > 1) else 0)

df.head()

	SP500	Log Daily Return	Log 10D Return	EWMA var	1D SD	10D SD	VaR (Return)	VaR to Compare	Breach	Consecutive breach
Date
05/02/2013	1511.290039	0.010363	0.012471	0.000078	0.008837	0.027944	-0.065008	-0.055164	0	0
06/02/2013	1512.119995	0.000549	0.011513	0.000056	0.007504	0.023729	-0.055202	-0.047174	0	0
07/02/2013	1509.390015	-0.001807	0.009700	0.000041	0.006439	0.020361	-0.047366	-0.040029	0	0
08/02/2013	1517.930054	0.005642	0.009911	0.000039	0.006226	0.019688	-0.045801	-0.040007	0	0
11/02/2013	1517.010010	-0.000606	0.011156	0.000028	0.005293	0.016736	-0.038935	-0.034704	0	0

# Count the number and % of breaches
breaches = df['Breach'].sum()
percent_breaches = breaches / df['Breach'].count()
consecutive_breaches = df['Consecutive breach'].sum()
percent_consecutive_breaches = consecutive_breaches / df['Consecutive breach'].count()

print('Number of breaches: ', breaches)
print('% of breaches: ', f'{round(percent_breaches*100,3)}%')
print('Number of consecutive breaches: ', consecutive_breaches)
print('% of consecutive breaches: ', f'{round(percent_consecutive_breaches*100,3)}%')

Number of breaches:  32
% of breaches:  2.581%
Number of consecutive breaches:  17
% of consecutive breaches:  1.371%

fig = px.bar(df, x=df.index, y='Log 10D Return', title='Backtesting VaR')
fig.update_traces(marker_color='blue',
                  marker_line_color='blue',
                  selector=dict(type="bar"))
fig.add_scatter(x=df.index, y=df['VaR to Compare'], name='VaR - 10D', mode='lines', opacity=1, marker_color='red')
fig.add_scatter(x=df[df['Breach'] > 0].index, y=df[df['Breach'] > 0]['Log 10D Return'], name='Breach', mode='markers', opacity=1)
fig.show()