STATS 509: Homework 4

Question 1

Suppose \(X\) is Pareto with parameters \(a=20\) and \(\mu=1,\) and \(Z\) is random variable independent of \(X\) satisfying that \(P(Z=-1)=P(Z=1)=\frac{1}{2}\)

1. Compute Pearson and Spearman correlation coefficients for \(X\) and \(Y=Z X\).
2. Compute Pearson and Spearman correlation coefficients for \(X\) and \(Y=X^{2}\).
3. Compute Pearson and Spearman correlation coefficients for \(X\) and \(Y=X^{9}\).

\(\mathbf{Solution.}\qquad\) It will be useful to compute a couple moments for \(X\) before we proceed to compute the correlations. First we are given that the distribution of \(X\) is given by, \[ X\sim F(x)=\begin{cases} 1-\left(\frac{1}{x}\right)^{20} & \text{if }x\geq1\\ 0 & \text{if }x<1 \end{cases} \]

Then the \(k\)-th moment of \(X\) is given by, \[\begin{align*} \mathbb{E}\left[X^{k}\right] & =20\int_{1}^{\infty}x^{k}x^{-21}dx\\ & =20\int_{1}^{\infty}x^{k-21}dx\\ & =\frac{20}{k-20}\left[\left.x^{k-20}\right|_{1}^{\infty}=-1\right]\\ & =\frac{-20}{k-20}\\ & =\frac{20}{20-k}\quad\text{for }k\in\left\{ 1,\ldots,19\right\} \end{align*}\]

Thus it is simple to compute \(\mathrm{Var}\left(X\right)=\mathbb{\mathbb{E}}\left[X^{2}\right]-\mathbb{E}\left[X\right]^{2}=\frac{20}{18}-\left(\frac{20}{19}\right)^{2}\).

For \(Z\), we have that \[\begin{align*} \mathbb{E}\left[Z\right] & =\frac{1}{2}\left(1\right)+\frac{1}{2}\left(-1\right)=0\\ \\ \mathbb{E}\left[Z^{2}\right] & =\frac{1}{2}\left(1\right)^{2}+\frac{1}{2}\left(-1\right)^{2}=1 \end{align*}\]

Thus we can see that \(\mathrm{Var}\left(Z\right)=1\). Now we can proceed and compute the following correlations.

(a) Solution.

Pearson

\[ \rho:=\mathrm{Corr}\left(X,Y\right):=\frac{\mathrm{Cov}(X,ZX)}{\mathrm{SD}(X)\mathrm{SD}(ZX)} \] We compute each term separately in order to calculate the full correlation. First, \[\begin{align*} \mathrm{Cov}(X,ZX) & =\mathbb{E}\left[ZX^{2}\right]-\mathbb{E}\left[ZX\right]\mathbb{E}\left[X\right]\\ & =\mathbb{E}\left[Z\right]\mathbb{E}\left[X^{2}\right]-\mathbb{E}\left[Z\right]\mathbb{E}\left[X\right]^{2}\\ & =\mathbb{E}\left[Z\right]\left(\mathrm{Var}\left(X\right)\right)=0 \end{align*}\] The \(\mathrm{Var}\left(X\right)\) is a finite number. We check that \(\mathrm{Var}\left(ZX\right)\) is also finite. \[\begin{align*} \mathrm{Var}\left(ZX\right) & =\mathbb{E}\left[Z^{2}X^{2}\right]-\mathbb{E}\left[ZX\right]^{2}\\ & =\mathbb{E}\left[Z^{2}\right]\mathbb{E}\left[X^{2}\right]\\ & =\frac{20}{18} \end{align*}\] Thus \(\rho=0\).

Spearman

First we compute the distribution of \(Y=ZX\), \[\begin{align*} F_{Y}(y)=\mathbb{P}(Y\leq y) & =\mathbb{P}(ZX\leq y)\\ & =\mathbb{P}(ZX\leq y\mid Z=1)\mathbb{P}(Z=1)+\mathbb{P}(ZX\leq y\mid Z=-1)\mathbb{P}(Z=-1)\\ & =\frac{1}{2}\mathbb{P}(X\leq y)+\frac{1}{2}\mathbb{P}(-X\leq y)\\ F_{Y}(y) & =\frac{1}{2}F_{X}(y)+\frac{1}{2}\left(1-F_{X}(-y)\right) \end{align*}\] Note that \(\mathbb{E}\left[F_{X}(x)\right]=\mathbb{E}\left[F_{Y}(y)\right]=\frac{1}{2}\) since \(F_{X}(x),F_{Y}(y)\sim Unif\left(0,1\right)\). Next compute the Spearman correlation, \[\begin{align*} \rho_{S}(X,Y) & =\rho\left(F_{X}(X),F_{Y}(Y)\right)\\ & =\rho\left(F_{X}(X),\frac{1}{2}F_{X}(Y)+\frac{1}{2}\left(1-F_{X}(-Y)\right)\right) \end{align*}\]

For simplicity define \(F(x):=F_{X}(x)\). Let’s compute the covariance term of the correlation, \[\begin{align} & \mathrm{Cov}\left(F(X),\frac{1}{2}F(Y)+\frac{1}{2}(1-F(-Y))\right)\label{eq:1}\\ & =\frac{1}{2}\mathrm{Cov}(F(X),F(Y))-\frac{1}{2}\mathrm{Cov}(F(X),F(-Y))\nonumber \end{align}\] We will compute each of the covariance terms separately,\ \[\begin{align*} & \mathrm{Cov}(F(X),F(Y))\\ & =\mathbb{E}[F(X)F(Y)]-\mathbb{E}[F(X)]\mathbb{E}[F(Y)]\\ & =\mathbb{E}[F(X)F(ZX)]-\left(\frac{1}{2}\right)^{2}\\ & =\mathbb{E}[F(X)F(X)\mid Z=1]\mathbb{P}(Z=1)+\mathbb{E}[F(X)F(-X)\mid Z=-1]\mathbb{P}(Z=1)-\left(\frac{1}{2}\right)^{2}\\ & =\frac{1}{2}\mathbb{E}\left[F^{2}(X)\right]+\frac{1}{2}\mathbb{E}[F(X)F(-X)]-\left(\frac{1}{2}\right)^{2} \end{align*}\]

Similarly, \[\begin{align*} & \mathrm{Cov}(F(X),F(-Y))\\ & =\mathbb{E}[F(X)F(-Y)]-\mathbb{E}[F(X)]\mathbb{E}[F(-Y)]\\ & =\mathbb{E}[F(X)F(-ZX)]-\left(\frac{1}{2}\right)^{2}\\ & =\mathbb{E}[F(X)F(X)\mid Z=-1]\mathbb{P}(Z=-1)+\mathbb{E}[F(X)F(-X)\mid Z=1]\mathbb{P}(Z=1)-\left(\frac{1}{2}\right)^{2}\\ & =\frac{1}{2}\mathbb{E}\left[F^{2}(X)\right]+\frac{1}{2}\mathbb{E}[F(X)F(-X)]-\left(\frac{1}{2}\right)^{2} \end{align*}\]

Since both covariance terms are equal, then equation (\ref{eq:1}) is equal to \(0\). Thus, \(\rho_{S}(X,Y)=0\).

   

(b) Solution.

Pearson

\[ \rho:=\mathrm{Corr}\left(X,Y\right):=\frac{\mathrm{Cov}(X,X^{2})}{\mathrm{SD}(X)\mathrm{SD}(X^{2})} \] We compute each term separately in order to calculate the full correlation. First, \[\begin{align*} \mathrm{Cov}(X,X^{2}) & =\mathbb{E}\left[X^{3}\right]-\mathbb{E}\left[X\right]\mathbb{E}\left[X^{2}\right]\\ & =\frac{20}{17}-\frac{20}{19}\left(\frac{20}{18}\right) \end{align*}\]

Computation of the terms in the denominator is straightforward, so we plug in directly. Finally we have, \[\begin{align*} \rho:=\mathrm{Corr}\left(X,Y\right):=\frac{\frac{20}{17}-\frac{20}{19}\left(\frac{20}{18}\right)}{\sqrt{\frac{20}{18}-\left(\frac{20}{19}\right)^{2}}\sqrt{\frac{20}{16}-\left(\frac{20}{18}\right)^{2}}} & \approx 0.998 \end{align*}\]

Spearman

Since \(X\) is supported on \(x\geq1\), then \(Y=X^{2}\) is a monotonic transformation of \(X\) on this domain. Thus we readily have, \[\begin{align*} \rho_{S}= & \rho\left(F_{X}\left(X\right),F_{Y}\left(Y\right)\right)=1 \end{align*}\]

   

(c) Solution.

Pearson

\[ \rho:=\mathrm{Corr}\left(X,Y\right):=\frac{\mathrm{Cov}(X,X^{9})}{\mathrm{SD}(X)\mathrm{SD}(X^{9})} \] We compute each term separately in order to calculate the full correlation. First, \[\begin{align*} \mathrm{Cov}(X,X^{9}) & =\mathbb{E}\left[X^{10}\right]-\mathbb{E}\left[X\right]\mathbb{E}\left[X^{9}\right]\\ & =\frac{20}{10}-\frac{20}{19}\left(\frac{20}{11}\right) \end{align*}\]

Computation of the terms in the denominator is straightforward, so we plug in directly. Finally we have, \[\begin{align*} \rho:=\mathrm{Corr}\left(X,Y\right):=\frac{\frac{20}{10}-\frac{20}{19}\left(\frac{20}{11}\right)}{\sqrt{\frac{20}{18}-\left(\frac{20}{19}\right)^{2}}\sqrt{\frac{20}{2}-\left(\frac{20}{11}\right)^{2}}} & \approx 0.6 \end{align*}\]

Spearman

Since \(X\) is supported on \(x\geq1\), then \(Y=X^{9}\) is a monotonic transformation of \(X\) on this domain. Thus we readily have, \[\begin{align*} \rho_{S}= & \rho\left(F_{X}\left(X\right),F_{Y}\left(Y\right)\right)=1 \end{align*}\]

Question 2

Suppose \(R_{1}, R_{2}\) are returns on 2 assets, and suppose that \(E\left(R_{1}\right)=.02, \quad E\left(R_{2}\right)=.03,\text{ } \operatorname{Var}\left(R_{1}\right)=(.025)^{2}, \operatorname{Var}\left(R_{2}\right)=(.04)^{2},\) and \(\operatorname{Corr}\left(R_{1}, R_{2}\right)=0.5\).

1. What are \(E\left(0.4 R_{1}+0.6 R_{2}\right)\) and \(\operatorname{Var}\left(0.4 R_{1}+0.6 R_{2}\right) ?\)
2. Assuming that \(0 \leq w \leq 1,\) for what value of \(w\) is the expected value of \(\left(w \cdot R_{1}+(1-w) \cdot R_{2}\right)\) maximized? For what value is the variance of \(\left(w \cdot R_{1}+(1-w) \cdot R_{2}\right)\) minimized?
3. Assuming a portfolio of a \(\$ 1\) million, and a multi-variate normal distribution for \(R_{1}\) and \(R_{2},\) find the value \(w\) that minimizes the VaR at \(q=.005\) and why might that be useful? Also, report the associated expected shortfall with this portfolio. Hints: Note that the random variable \(\left[w R_{1}+(1-w) R_{2}\right]\) is normally distributed for any \(w,\) and it will be easiest to use R-package - the solution can be approximate, i.e., accurate to .01.
4. Repeat the above if have multi-variate t-distribution for \(R_{1}\) and \(R_{2},\) with \(\nu=5\).

(a) \(\mathbf{Solution.}\qquad\) \[\begin{align*} \mathbb{E}\left[0.4R_{1}+0.6R_{2}\right] & =0.4\left(0.02\right)+0.6\left(0.03\right)= 0.026\\ \\ \operatorname{Var}\left(0.4R_{1}+0.6R_{2}\right) & =(0.16)\mathrm{Var}\left(R_{1}\right)+0.36\mathrm{Var}\left(R_{2}\right)+2(0.4)(0.6)\mathrm{SD}\left(R_{1}\right)\mathrm{SD}\left(R_{2}\right)\mathrm{Corr}\left(R_{1},R_{2}\right)\\ & =(0.16)(0.025)^{2}+0.36(0.04)^{2}+2(0.4)(0.6)(0.025)(0.04)0.5\\ & = 0.001 \end{align*}\]

   

(b) \(\mathbf{Solution.}\qquad\) \[\begin{align*} \mathbb{E}\left[wR_{1}+\left(1-w\right)R_{2}\right] & =w\mathbb{E}\left[\mathbb{R}_{1}\right]+(1-w)\mathbb{E}\left[R_{2}\right]\\ & =0.02w+0.03(1-w)\\ & =0.03-0.01w \end{align*}\]

So the expected value is maximized when \(w=0\). Next we minimize the variance of the portfolio, \[\begin{align*} & \operatorname{Var}\left(wR_{1}+\left(1-w\right)R_{2}\right)\\ & =w^{2}\mathrm{Var}\left(R_{1}\right)+\left(1-w\right)^{2}\mathrm{Var}\left(R_{2}\right)+2w\left(1-w\right)\mathrm{SD}\left(R_{1}\right)\mathrm{SD}\left(R_{2}\right)\mathrm{Corr}\left(R_{1},R_{2}\right) \end{align*}\]

Taking derivative with respect to \(w\), \[\begin{align*} 2w\mathrm{Var}\left(R_{1}\right)-2(1-\omega)\mathrm{Var}\left(R_{2}\right)+2\sqrt{\mathrm{Var}\left(R_{1}\right)\mathrm{Var}(2)}\mathrm{Corr}\left(R_{1},R_{2}\right) & =0\\ -4w\sqrt{\mathrm{Var}\left(R_{1}\right)\mathrm{Var}(R_{2})}\mathrm{Corr}\left(R_{1},R_{2}\right)\\ 2w\left(\mathrm{Var}\left(R_{1}\right)-2\sqrt{\mathrm{Var}\left(R_{1}\right)\mathrm{Var}(R_{2})}\mathrm{Corr}\left(R_{1},R_{2}\right)+\mathrm{Var}\left(R_{2}\right)\right) & =0\\ -2\mathrm{Var}\left(R_{2}\right)+2\sqrt{\mathrm{Var}\left(R_{1}\right)\mathrm{Var}(R_{2})}\mathrm{Corr}\left(R_{1},R_{2}\right) \end{align*}\]

Substituting the known values, \[\begin{align*} 0.001w & = 0.001 \end{align*}\]

So the variance is minimized when \(w=1\).

   

(c) \(\mathbf{Solution.}\qquad\) We iterate through a sequence of values of \(w\in[0,1]\). Finding \(w\) such that VaR is minimized for \(q=0.005\) allows us minimize a large potential loss in our portfolio from a rare event. In the code below we set up our paramers and compute the VaR.

# QUESTION 2C -----------------------------------------------------------------
set.seed(509)
w = seq(0,1,.01)
n = length(w)
ES = rep(0, n)
draws = 2 * 1e7
capital = 1e6

# Set up parameters
sigma1 = 0.025; sigma2 = 0.04
mu1 = .02; mu2 = .03
rho = 0.5
alpha = 0.005

# Construct mu vector and covariance matrix
cov_12 = sigma1 * sigma2 * rho
mu_R = matrix(c(mu1, mu2), nrow = 2)
Sigma_R = matrix(c(sigma1^2, cov_12, 
                 cov_12, sigma2^2), nrow = 2, ncol = 2)

# Compute Vector of Averages and SDs
W = cbind(w, 1 - w)
mu = W %*% mu_R
sigma = sqrt(diag(W %*% Sigma_R %*% t(W)))

# Vector of VaRs
VaR = -qnorm(alpha, mean = mu, sd = sigma)

# Compute ES
i_min = which.min(VaR)
returns = rnorm(draws, mean = mu[i_min], sd = sigma[i_min])
ES = -mean(returns[returns < -VaR[i_min]])

# Results
w_star = w[i_min]
r_VaR_min = VaR[i_min]
VaR_min = r_VaR_min * capital 
shortfall = ES * capital

We find that the optimal \(w\) that minimizes VaR is \(w^*=\) 0.82. The corresponding relative VaR is \(\tilde{\mathrm{VaR}}_{0.005}=\) 0.042, so \(\mathrm{VaR}_{0.005}=\) $42321.6. Finally, we find that shortfall is about $50256.58.

(d) \(\mathbf{Solution.}\qquad\) We perform same analysis as part (c) using t-distribution with \(\nu=5\),

# QUESTION 2D -----------------------------------------------------------------
nu = 5
VaR_t = -(sigma * sqrt((nu - 2) / nu) * qt(alpha, nu) + mu)

i_min = which.min(VaR_t)
returns = sigma[i_min] * sqrt((nu - 2) / nu) * rt(draws, nu) + mu[i_min]
ES_t = -mean(returns[returns < -VaR_t[i_min]])

# Results
w_star = w[i_min]
r_VaRt_min = VaR_t[i_min]
VaRt_min = r_VaRt_min * capital 
shortfall = ES_t * capital

We find that the optimal \(w\) that minimizes VaR is \(w^*=\) 0.83. The corresponding relative VaR is \(\tilde{\mathrm{VaR}}_{0.005}=\) 0.056, so \(\mathrm{VaR}_{0.005}=\) $55937.5. Finally, we find that shortfall is about $79310.98.

Question 3

The data for this problem is in the Data folder on Canvas.

1. Generate pairwise scatter diagram on the 4 assets in Problem 1 on page 177 of Ruppert / Matteson via the R-commands given there, and summarize your findings/observations from these plots.

   library(MASS) # needed for cov.trob
   library(mnormt) # needed for dmt
   df = seq(2.5, 8, 0.01)
   n = length(df)
   loglik_profile = rep(0, n)
   for (i in 1:n) {
     fit = cov.trob(Berndt, nu = df[i])
     mu = as.vector(fit$center)
     sigma = matrix(fit$cov, nrow = 4)
     loglik_profile[i] = sum(log(dmt(Berndt, mean = fit$center,
                                     S= fit$cov, df = df[i])))
     }

2. Do Problem 2 on page 178 of Ruppert/Matteson. Using the results produced by the code above, find the MLE of \(\nu\) and a \(90 \%\) profile likelihood confidence interval for \(\nu .\) Include your R code with your work. Also, plot the profile log-likelihood and indicate the MLE and the confidence interval on the plot.
3. Generate \(\mathrm{QQ}\) plots for the estimated distributions in (b) relative to the asset returns and summarize what they show.

4. Based on (b) and using \(q=.005,\) derive an estimate of the VaR and expected shortfall of a portfolio of 10 million dollars evenly split between the 4 assets.

(a) \(\mathbf{Solution.}\qquad\) In figure 1 we show pair plot of the assets. We see that CITCRP is highly correlated with CONTIL and DATGEN. On the other hand, CONED has the lowest correlation with CONTIL and DATGEN.

# QUESTION 3A -----------------------------------------------------------------
berndtInvest = read.csv("berndtInvest.csv")
Berndt = berndtInvest[, 2:5]

# Generate Pair Plot
ggpairs(Berndt, axisLabels = "none",
        upper = list(continuous = wrap("points", 
                                       colour = "steelblue"), combo = "dot"),
        lower = list(continuous = "cor", combo = "dot"),
        diag = list(continuous = function(data, mapping, ...) {
          ggally_barDiag(data = data, mapping = mapping, 
                         alpha = .9, fill = "steelblue", bins = 20)
          }))

Figure 1: Pair Plot of Assets

   

(b) \(\mathbf{Solution.}\qquad\) In the R code below we take the maximum of the profile log likelihood to find \(\hat{\nu}_{MLE}\). To get a 90% confidence interval, we intersect the horizontal line in figure 2 which is \(\ell_{max}(\hat{\nu}_{MLE}) - \chi^2_{0.9,1}\) with the profile log likelihood.

# QUESTION 3B -----------------------------------------------------------------
library(MASS) # needed for cov.trob
library(mnormt) # needed for dmt
df = seq(2.5, 8, 0.01)
n = length(df)
loglik_profile = rep(0, n)
for (i in 1:n) {
  fit = cov.trob(Berndt, nu = df[i])
  mu = as.vector(fit$center)
  sigma = matrix(fit$cov, nrow = 4)
  loglik_profile[i] = sum(log(dmt(Berndt, mean = fit$center,
                                  S = fit$cov, df = df[i])))
}

# Find nu that maximizes loglik profile
dt = data.frame("nu" = df, "loglik" = loglik_profile)

i_max = which.max(loglik_profile)

nu_MLE = df[i_max]


# Split profile curve into two halves to find intersection
left = loglik_profile[1:i_max]
right = loglik_profile[i_max + 1:length(loglik_profile)]

# Build CI
h = max(loglik_profile) - 0.5*qchisq(0.9, df = 1)
ci_left = df[which.min(abs(left - h))]
ci_right = df[i_max + which.min(abs(right - h))]

# Label MLE
mle_point = dt[i_max,]

# Create Plot
ggplot(dt, aes(x = nu, y = loglik)) + 
  geom_line(colour = "steelblue") + 
  geom_hline(yintercept = h, linetype = "dashed", colour = "darkred") + 
  geom_point(data=mle_point, colour="darkred") +  # this adds a red point
  geom_text(data=mle_point, label="MLE", vjust = 2) +
  geom_ribbon(data = subset(dt, ci_left <= nu & nu <= ci_right),
              aes(ymin=min(loglik_profile), ymax=loglik), alpha = 0.5, 
              fill = "steelblue") + 
  ggtitle(TeX("Profile Log-likelihood for $\\nu$")) + 
  xlab(TeX("$\\nu$")) + 
  ylab('log likelihood') + 
  theme(plot.title = element_text(hjust = 0.5, size = 10))

Figure 2: A profile likelihood confidence interval for \(\nu\).

From figure 2, we determine that \(\hat{\nu}_{MLE} =\) 4.55 and \(\ell_{max}(\hat{\nu}_{MLE}) - \chi^2_{0.9,1} =\) 537.827. Our 90% confidence interval on \(\nu\) is [3.22,6.59]. In the next problems we will use \(\nu =\) 5 when we fit t-distributions.

   

(c) \(\mathbf{Solution.}\qquad\) In the QQ plots displayed in figure 3, we see that the tails CONTIL are more heavy tailed than the tails of our fitted t-distribution. On the other hand, we are slightly overfitting the tails of CITCRP. Overall the tails for CONED and DATGEN are fitted well.

# QUESTION 3C -----------------------------------------------------------------
# Function for generating t-plot of scaled t distribution
qqt_scale = function(x, nu, asset){
  xsort = sort(x)
  n = length(xsort)
  q = (1 / (n + 1)) * c(1:n)
  qq = qt(q, nu)
  
  # Construct dataframe to plot in ggplot
  df = data.frame("quantiles" = qq, "data" = xsort)
  
  # Compute theoretical quantiles
  xx = quantile(qq, c(0.25, 0.75))
  yy = quantile(xsort, c(0.25, 0.75))
  slope = (yy[2] - yy[1]) / (xx[2] - xx[1])
  inter = yy[1] - slope*xx[1]
  
  g = ggplot(df, aes(x = quantiles, y = data)) +
    geom_point(colour = "darkred", shape = 1) + 
    ggtitle(paste0("t-Plot of ", asset)) + 
    theme(text = element_text(size = 7),
          plot.title = element_text(hjust = 0.5)) + 
    geom_abline(intercept = inter, slope = slope, colour = "steelblue") + 
    xlab('Theoretical Quantiles') + 
    ylab('Sample Quantiles')
  
  return(g)
}

# Construct qq plots for each asset
nu_MLE = round(nu_MLE)
plot_list = list()
i = 1

for (col in colnames(Berndt)) {
  plot_list[[i]] = local({
    p1 = qqt_scale(as.matrix(Berndt[col]), nu_MLE, col)
  })
  i = i + 1
}

# plot all subplots
n <- length(plot_list)
nCol <- ceiling(sqrt(n))
grid.arrange(grobs = plot_list, ncol = nCol,
             top = textGrob("QQ-Plot of Assets with t-distribution"),
                            gp = gpar(fontsize = 13, font = 8))

Figure 3: QQ plots of assets using fitted t-distribution.

   

(d) \(\mathbf{Solution.}\qquad\) Below we estimate the VaR of our portfolio returns when \(q=0.005\) using our fitted t-distribution and compute the expected shortfall with a portfolio of 10 million dollars.

# QUESTION 3D -----------------------------------------------------------------
capital = 1e7
w = matrix(rep(1 / 4, 4), ncol = 1)
mu = as.numeric(colMeans(Berndt) %*% w)
sigma = as.numeric(t(w) %*% cov(Berndt) %*% w)

# Compute relative VaR
nu = nu_MLE
r_VaR_t = as.numeric(-(sigma * sqrt((nu - 2) / nu) * qt(alpha, nu) + mu))

# Compute Expected Shortfall
returns = sigma * sqrt((nu - 2) / nu) * rt(draws, nu) + mu
ES_t = -mean(returns[returns < -r_VaR_t])

# VaR and Shortfall
VaR_t = capital * r_VaR_t
shortfall_t = capital * ES_t

The corresponding relative VaR is \(\tilde{\mathrm{VaR}}_{0.005}=\) 0.008, so \(\mathrm{VaR}_{0.005}=\) $76559.81. Finally, we find that shortfall is about $127423.5.