R

It is easy to replicate the calculations we did above by hand using R. First, we need to collect our stock data. We can do this using the getSymbols function. The input to the getSymbols function is the stock tickers from the New York Stock Exchange accessed through Yahoo Finance ®. It will pull all daily opening, closing, high, and low prices for the given stock tickers and save them in dataframes under the same names. The last function will extract only the last number of observations from the data sets. Here we are looking at the last 500 days. Every day you run this code, the last 500 days will change so the calculations could be different. Specifically we are looking only at the closing prices for each stock. The ROC function will give us the daily returns for each of these stocks now saved as the variables msft_r and appl_r.

tickers = c("AAPL", "MSFT")

getSymbols(tickers)

[1] "AAPL" "MSFT"

stocks <- cbind(last(AAPL[,4], '500 days'), last(MSFT[,4], '500 days'))

stocks$msft_r <- ROC(stocks$MSFT.Close)
stocks$aapl_r <- ROC(stocks$AAPL.Close)

head(stocks)

           AAPL.Close MSFT.Close       msft_r       aapl_r
2019-02-26    43.5825     112.36           NA           NA
2019-02-27    43.7175     112.17 -0.001692451  0.003092763
2019-02-28    43.2875     112.03 -0.001248876 -0.009884593
2019-03-01    43.7425     112.53  0.004453160  0.010456306
2019-03-04    43.9625     112.26 -0.002402217  0.005016874
2019-03-05    43.8825     111.70 -0.005000948 -0.001821436

Next, we just calculate the historical numbers we need for our calculations, namely the correlation between and the variance of the stocks. Instead of the correlation, we could also use the covariance directly. The cov, var, and cor functions give us the covariance, variances, and correlations respectively. From there we just recreate the equation written above in code. Lastly, we use the qnorm function to get the quantile from the normal distribution for our VaR calculation.

var.msft <- var(stocks$msft_r, na.rm=TRUE)
var.aapl <- var(stocks$aapl_r, na.rm=TRUE)
cov.m.a <- cov(stocks$msft_r, stocks$aapl_r, use="pairwise.complete.obs")
cor.m.a <- cor(stocks$msft_r, stocks$aapl_r, use="pairwise.complete.obs")

VaR.percentile = 0.01
AAPL.inv <- 100000
MSFT.inv <- 200000

var.port <- (MSFT.inv/(MSFT.inv+AAPL.inv))^2*var.msft + (AAPL.inv/(MSFT.inv+AAPL.inv))^2*var.aapl + 2*(AAPL.inv/(MSFT.inv+AAPL.inv))*(MSFT.inv/(MSFT.inv+AAPL.inv))*cov.m.a

VaR.DN.port <- (AAPL.inv+MSFT.inv)*qnorm(VaR.percentile)*sqrt(var.port)
dollar(VaR.DN.port)

[1] "-$14,729.81"

If the number above is different than the number we calculated by hand, then it is simply a result of a different time period and different data underlying it. The code is written to always use the most up to date data for each stock.

Another great aspect of the VaR under the normality assumption is the ease of calculating any \(n\)-day VaR from the 1-day VaR. This is defined simply as \(VaR_n = \sqrt(n) \times VaR_1\). Another great reason why the normality assumption in the Delta-Normal approach makes the calculation of VaR numbers rather easy.

SAS

We will use the functions in the R section of the code to get the necessary data in a dataset called stocks. First, we will use PROC CORR to get us the necessary information on the variances, covariances, and correlations. You could easily just write these numbers down, but we are going to use the out = option to store the needed information into a dataset called covar_results. We then use the DATA step to create MACRO variables that store the values of the covariance between and variances of the stocks.

proc corr data=SimRisk.stocks cov out=covar_results plots;
    var msft_r aapl_r;
run;

data _null_;
    set covar_results(where = (_type_ eq 'COV'));
    if upcase(_name_) eq 'MSFT_R' then do;
        call symput("var_msft",msft_r);
        call symput("covar",aapl_r);
    end;
    if upcase(_name_) eq 'AAPL_R' then do;
        call symput("var_aapl",aapl_r);
    end;
run;

data _null_;
    set covar_results(where = (_type_ eq 'CORR'));
    if upcase(_name_) eq 'MSFT_R' then do;
        call symput("corr",aapl_r);
    end;
run;

Next, we need to use PROC IML to do the calculations by hand. PROC IML is the interactive matrix language of SAS. It is similar in design to R. We start by defining some MACRO variables for the amount we invested in each stock as well as the desired VaR percentile. From there we invoke PROC IML and put in calculations as we would in R. We print out the results of these equations that we defined above.

%let msft_position = 200000;
%let aapl_position = 100000;
%let var_percentile = 0.01;

proc iml;
title 'VaR Results';

msft_p_weight = &msft_position/(&msft_position + &aapl_position);
aapl_p_weight = &aapl_position/(&msft_position + &aapl_position);

P_variance =  &var_msft*(msft_p_weight)**2 + &var_aapl*(aapl_p_weight)**2 
                   + 2*aapl_p_weight*msft_p_weight*&covar;
P_StdDev=sqrt(P_variance);

VaR_normal = (&msft_position + &aapl_position)*PROBIT(&var_percentile)*SQRT(P_variance);

print "Daily VaR (Percentile level: &var_percentile); Delta-Normal" VaR_normal[format=dollar15.2];

quit;

If the number above is different than the number we calculated by hand, then it is simply a result of a different time period and different data underlying it. The code is written to always use the most up to date data for each stock.

Another great aspect of the VaR under the normality assumption is the ease of calculating any \(n\)-day VaR from the 1-day VaR. This is defined simply as \(VaR_n = \sqrt(n) \times VaR_1\). Another great reason why the normality assumption in the Delta-Normal approach makes the calculation of VaR numbers rather easy.

R

Running this simulation is rather easy to do in R. First, we define our correlation structure (here R) and calculate the Cholesky decomposition of this matrix to multiply into our data. Next, we simulate 10,000 draws from the respective normal distributions for our Microsoft and Apple returns using the rnorm function. We multiply in the correlation structure and then simply calculate the return on the portfolio of these two sets of returns. From there it is a simple call of the quantile function to find the quantile of the distribution to get the VaR. A histogram of the distribution of changes in the portfolio is also provided.

n.simulations <- 10000
R <- matrix(data=cbind(1,cor.m.a, cor.m.a, 1), nrow=2)
U <- t(chol(R))

msft.r <- rnorm(n=n.simulations, mean=0, sd=sqrt(var.msft))
aapl.r <- rnorm(n=n.simulations, mean=0, sd=sqrt(var.aapl))
Both.r <- cbind(msft.r, aapl.r)
port.r <- U %*% t(Both.r)
port.r <- t(port.r)

value.change <- MSFT.inv*port.r[,1] + AAPL.inv*port.r[,2]
value <- (MSFT.inv + AAPL.inv) + value.change

VaR <- quantile(value.change, VaR.percentile, na.rm=TRUE)
VaR.label <- dollar(VaR)

hist(value.change/1000, breaks=50, main='1 Day Value Change Distribution', xlab='Value Change', col="lightblue")
breaks = c(-40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60)
axis(1, at = breaks, labels = paste("$", breaks, "K", sep = ""))
abline(v = VaR/1000, col="red", lwd=2)
mtext(paste("Value at Risk",VaR.label, sep=" = "), at=VaR/1000, col="red")

As you can see from above, the one day, 1% VaR for this portfolio (assuming normality of returns for each asset) is $14,476.19.

SAS

Running this simulation is a little complicated in SAS due to correlation structure we are imposing on the returns. For more detail on the code below and the approach to adding in correlation structures in simulations in SAS please see the code notes for Introducion to Simulation. First, we define our correlation structure using the DATA step. This will repeat the desired correlation matrix the number of times we are simulating. Next, we again use the DATA step to define starting points for each of the variables in the simulation - namely the returns for Microsoft and Apple. From there PROC MODEL simulates and correlates the necessary returns.

data Corr_Matrix;
    do i=1 to 10000;
        _type_ = "corr";
        _name_ = "msft_r";

        msft_r = 1.0; 
        aapl_r = &corr;

        output;
        _name_ = "aapl_r";

        msft_r = &corr; 
        aapl_r = 1.0;

        output;
    end;
run;

data Naive;
    do i=1 to 10000;
        msft_r=0;
        aapl_r=0;
        output;
    end;
run;

proc model noprint;

    msft_r = 0;
    errormodel msft_r ~Normal(&var_msft);

    aapl_r = 0;
    errormodel aapl_r ~Normal(&var_aapl);

    solve msft_r aapl_r/ random=1 sdata=Corr_Matrix
    data=Naive out=mc_stocks(keep=msft_r aapl_r i);
        by i;

run;
quit;

data mc_stocks;
    set mc_stocks;
    by i;
    if first.i then delete;
    rename i=simulation;
    value_change = &msft_position*msft_r + &aapl_position*aapl_r;
    value = &msft_position + &aapl_position + value_change;
    format value_change dollar15.2;
run;

%let var_clevel=%sysevalf(100*&var_percentile);
    
proc univariate data=mc_stocks;
    var value_change;
    format value_change dollar15.2;
    output out=percentiles pctlpts = &var_clevel pctlpre=P;
    histogram value_change / kernel normal;
  ods select Quantiles;
run;

As you can see from above, the one day, 1% VaR for this portfolio (assuming normality of returns for each asset) is $14,615.71.

R

Let’s first start with the direct normal distribution calculation approach. Under the Delta-Normal approach you can easily calculate the confidence interval for the value at risk based on the equation defined above. The code below just recreates that function using the qchisq and qnorm functions were necessary to get values from the \(\chi^2\) and normal distributions.

sigma.low <- sqrt(var.port*(length(stocks$AAPL.Close)-1)/qchisq((1-(VaR.percentile/2)),length(stocks$AAPL.Close)-1) )
sigma.up <- sqrt(var.port*(length(stocks$AAPL.Close)-1)/qchisq((VaR.percentile/2),length(stocks$AAPL.Close)-1) )

VaR.DN.port <- (AAPL.inv+MSFT.inv)*qnorm(VaR.percentile)*sqrt(var.port)
VaR.L <-  (AAPL.inv+MSFT.inv)*qnorm(VaR.percentile)*(sigma.low)
VaR.U <- (AAPL.inv+MSFT.inv)*qnorm(VaR.percentile)*(sigma.up)

print(paste("The VaR is ", dollar(VaR.DN.port), " with confidence interval (", dollar(VaR.L), " , ", dollar(VaR.U), ")", sep = "" ))

[1] "The VaR is -$14,729.81 with confidence interval (-$13,614.27 , -$16,029.05)"

To perform the bootstrap approach in R, we first need to define the number of bootstrap samples as well as the size of each of these samples. From there we loop through the simulation values we already obtained earlier called value.change using a for loop. Using the sample function we can sample values from our original simulation and calculate the VaR of this new sub-sample using the quantile function at our percentile. We repeat this process over and over until we form many different values for our VaR from which to build our confidence interval. From there we simply look at the values that define the middle 95% of all of our VaR values, again using the quantile function to define the upper and lower bounds of our confidence interval. Lastly, we plot the VaR as well as its confidence interval bounds on a histogram of the value.change distribution.

n.bootstraps <- 1000
sample.size <- 1000

VaR.boot <- rep(0,n.bootstraps)
for(i in 1:n.bootstraps){
  bootstrap.sample <- sample(value.change, size=sample.size)
  VaR.boot[i] <- quantile(bootstrap.sample, VaR.percentile, na.rm=TRUE)
}

VaR.boot.U <- quantile(VaR.boot, 0.975, na.rm=TRUE)
VaR.boot.L <- quantile(VaR.boot, 0.025, na.rm=TRUE)

print(paste("The VaR is ", dollar(mean(VaR.boot)), " with confidence interval (", dollar(VaR.boot.L), " , ", dollar(VaR.boot.U), ")", sep = "" ))

[1] "The VaR is -$14,033.25 with confidence interval (-$15,606.72 , -$12,821.69)"

hist(value.change/1000, breaks=50, main='1 Day Value Change Distribution', xlab='Value Change', col="lightblue")
breaks = c(-20, -10, 0, 10, 20)
axis(1, at = breaks, labels = paste("$", breaks, "K", sep = ""))
abline(v = VaR/1000, col="red", lwd=2)
mtext(paste("Value at Risk",VaR.label, sep=" = "), at = VaR/1000, col="red")
abline(v = VaR.boot.L/1000, col="blue", lwd=2, lty="dashed")
abline(v = VaR.boot.U/1000, col="blue", lwd=2, lty="dashed")

SAS

It is actually rather easy to calculate the confidence intervals for a percentile in SAS under the normality assumption. PROC UNIVARIATE and the cipctlnormal option calculates confidence intervals for percentiles of any variable. Here we define the value_change variable from our dataset to make this calculation.

proc univariate data=mc_stocks cipctlnormal;
    var value_change;
  ods select Quantiles;
run;

To perform the bootstrap approach in SAS, we first need to define the number of bootstrap samples as well as the size of each of these samples. From there we use PROC SURVEYSELECT to sample from our dataset many times. The method = srs option specifies that we want to use simple random sampling. We want sample sizes that are 10% of the size of the simulation which is defined in the samprate = option. Lastly, we repeat this sampling over and over again using the rep = option. Make sure to put noprint so you don’t see results from each of these 1000 samples!

Next, we just calculate percentiles using PROC UNIVARIATE for each of the samples (called replicates) using the by statement. We will output the percentiles from each of these calculations using the output statement with the out = option. Specifically, we want the 1% quantile (our VaR values) which we can define using the pctlpts option. We will label this quantile as P1. We define this name using the pctlpre option to define the prefix on the number of the quantile (here 1). Another PROC UNIVARIATE call will give us the percentiles of all of these VaR calculations from the previous PROC UNIVARIATE. Lastly, we print these results.

%let n_bootstraps=1000;
%let bootstrap_prop=0.1;
%let var_percentile = 0.01;
%let var_clevel=%sysevalf(100*&var_percentile);

proc surveyselect data=mc_stocks out=outboot seed=12345 method=srs 
                  samprate=&bootstrap_prop rep=&n_bootstraps noprint;
run; 

proc univariate data=outboot noprint;
    var value_change;
    output out=boot_var pctlpts = &var_clevel pctlpre=P;
    by replicate;
run;

proc univariate data=boot_var;
    var P1;
    output out=var_ci pctlpts = 2.5 97.5 pctlpre=P;
  ods select Quantiles;
run;

proc print data=var_ci;
run;

R

To get the estimated parameters from the Pareto distribution we can use the gpdFit function. The first input is the data itself (careful, the function requires positive values) and the second input input is the type = c("mle") option to specify using maximum likelihood estimation to estimate the parameters. From there we use the tailRisk function to estimate the VaR and ES from the output of the gpdFit function.

pareto.fit <-gpdFit(as.numeric(stocks$port_v[-1])*-1, type = c("mle"))
tailRisk(pareto.fit, prob = 0.99)

    Prob      VaR       ES
95% 0.99 19419.31 27720.33

dollar(tailRisk(pareto.fit, prob = 0.99)[,2]*-1)

[1] "-$19,419.31"

dollar(tailRisk(pareto.fit, prob = 0.99)[,3]*-1)

[1] "-$27,720.33"

We can see the values from the Pareto estimation above. To compare this to the normal distribution, the following plot shows the estimated distribution (in the histogram) of the Pareto distribution compared to the theoretical normal distribution (the density line) for the same data.

hist(-1*rgpd(n = 10000, mu = pareto.fit@parameter$u, beta = pareto.fit@fit$par.ests[2], xi = pareto.fit@fit$par.ests[1]), breaks = 50, main = "Comparison of Simulated tail to Normal", xlab = "Portfolio Value", freq = FALSE, yaxt = "n", ylab = "", col = "lightblue")
curve(dnorm(x, mean(stocks$port_v, na.rm = TRUE), sd(stocks$port_v, na.rm = TRUE)), add = TRUE, col = 'red', lwd = 2)

As we can see above, the data estimates a distribution with a much thicker tail than the normal distribution would otherwise account for. This aligns with the historical simulation approach that had a much riskier value for VaR and ES compared to the Delta-Normal approach. All of the approaches are compared in the table below:

Comparison of VaR and ES Across Approaches

SAS

In SAS, we first need to isolate the tail fo the data to fit the Pareto distribution with. SAS requires that the data values are all positive so we convert our negative values to positive using the DATA step and multiplying our values by -1. Next, we use PROC UNIVARIATE to calculate the percentiles of our data. We want to isolate the 5% tail to build our Pareto distribution off of. We specify the variable of interest, neg_port_v, in our VAR statement. We then use the OUTPUT statement to save the percentiles into a new dataset called percentiles by using the pctlpts = 95 and pctlpre = P options we have defined before. Next, we use another DATA step to create a MACRO variable for the value of the 5% quantile of interest (remember we are working with the 95% instead of 5% because we flipped all of our values when taking the negative). We then use a last DATA step to delete all the data values not in the tail of interest.

%let var_percentile=0.01;

data neg_stocks;
    set SimRisk.stocks;
    neg_port_v = -1*port_v;
run;

proc univariate data=neg_stocks;
    var neg_port_v;
    output out=percentiles pctlpts = 95 pctlpre=P;
  ods select Quantiles;
run;

proc print data=percentiles;
run;

data _null_;
    set percentiles;
    call symput("tail",P95);
run;

data stocks_tail;
    set neg_stocks;
    if neg_port_v <= &tail then delete;
run;

Now we can use PROC UNIVARIATE to fit our data to the Pareto distribution. To do this we use the HISTOGRAM statement on our variable neg_port_v with the pareto(theta=&tail percent=80) option. This will fit the parameters of the Pareto distribution to our dataset. We can output these parameters using the ODS OUTPUT statements. The output from PROC UNIVARIATE provides the parameter estimates for the Pareto distribution as well as the estimate of our quantile of interest, the VaR. We save this VaR as a MACRO variable and print it out using PROC IML.

proc univariate data=stocks_tail;
    histogram neg_port_v / pareto(theta=&tail percent=80);
    ods output ParameterEstimates=PE FitQuantiles=FQ;
    ods select histogram ParameterEstimates FitQuantiles;
run;

data _null_;
    set PE;
    if upcase(Symbol) eq 'THETA' then do;
        call symput("threshold",Estimate);
    end;
    if upcase(Symbol) eq 'SIGMA' then do;
        call symput("beta",Estimate);
    end;
    if upcase(Symbol) eq 'ALPHA' then do;
        call symput("xi",Estimate);
    end;
run;

data _null_;
    set FQ;
    call symput("VaR",EstQuantile);
run;

proc iml;
title 'EVT Results';

VaR_EVT = &VaR*-1;

print "Daily VaR (Percentile level: &var_percentile); EVT" VaR_EVT[format=dollar15.2];

quit;

Lastly, we use one final PROC IML call to calculate the ES using the equation defined above.

proc iml;
title 'EVT Results';

VaR_EVT = &VaR*-1;
ES_EVT = ((&VaR + &beta - &xi*&threshold)/(1-&xi))*-1;

print "Daily CVaR/ES (Percentile level: &var_percentile); EVT" ES_EVT[format=dollar15.2];

quit;

These values align with the historical simulation approach that had a much riskier value for VaR and ES compared to the Delta-Normal approach. All of the approaches are compared in the table below:

Comparison of VaR and ES Across Approaches

Delta-Gamma VaR

The Delta-Gamma approximation approach tries to improve on the Delta-Normal approach by taking into account higher derivatives than the first derivative only. In the expansion of the derivatives we can look at the first two:

\[ dV = \frac{\partial V}{\partial RF}\cdot dRF + \frac{1}{2}\cdot \frac{\partial^2 V}{\partial RF^2}\cdot dRF^2 + \cdots \]

\[ dV = \Delta \cdot dRF + \frac{1}{2} \cdot \Gamma \cdot dRF^2 \]

Essentially, this approach takes the variance component into account in building the Gamma portion of the equation. Delta-Normal typically underestimates VaR when assumptions do not hold. The Delta-Gamma approach tries to correct for this, but it is computationally more intensive.

Delta-Gamma Monte Carlo

The assumption of normality still might not be reliable even with the Delta-Gamma adjustment to Delta-Normal. The Delta-Gamma MC simulation approach simulates the risk factor before using Taylor-series expansion to create simulated movements in the portfolio’s value. Just like with regular simulation, the VaR is calculated using the simulated distribution of the portfolio.

For a comparison, if you have a large number of things to evaluate in a portfolio with few options, the Delta-Normal approach is fast and efficient enough. However, if you have only a few sources of risk and substantial complexity, the Delta-Gamma approach is better. If you have both a large number of risk factors and substantial complexity, simulation is probably the best approach.

Marginal VaR

The marginal VaR tries to answer the question of how much the VaR will change if we invest one more dollar in a position of the portfolio. In other words, the first derivative of VaR with respect to the weight of a certain position in the portfolio. Assuming normality, this calculation is quite easy:

\[ \Delta VaR = \alpha \times \frac{Cov(R_i, R_p)}{\sigma_p} \]

Incremental VaR

The incremental VaR tries to measure the change in VaR due to the addition of a new position in the portfolio. This is used whenever an institution wants to evaluate the effect of a proposed trade or change in the portfolio of interest. This is slightly different from the marginal VaR that only takes into account an additional dollar of one position instead of adding or subtracting an entire (or multiple) positions.

This is trickier to calculate as you now need to calculate the VaR under both portfolios and take their difference. There are Taylor series approximations to his, but they are not as accurate as calculating the full VaR for both portfolios and looking at the difference.

Component VaR

The component VaR decomposes the VaR into its basic components ina way that takes into account the diversification of the portfolio. It is a “partition” of the portfolio VaR that indicates the change of VaR if a given component was deleted. The component VaR is defined in terms of the marginal VaR:

\[ VaR_{Comp} = \Delta VaR \times (\$ \; value \; of \; component) \]

The sum of all component VaR’s is the total portfolio VaR.