Quantitative Risk Management in R

1 Exploring market risk-factor data

https://www.qrmtutorial.org/

1.1 Introduction

1.1.1 Exploring risk-factor time series: equity indexes

In this exercise, you will look at an equity index and plot it for a particular range of dates. The data used in this exercise and in the rest of the course are contained in the package qrmdata. You also need the package xts to manipulate time series.

When the qrmdata library is attached, as it will be throughout the course, you can load a dataset with the data() command. For example, the command data(“FTSE”) loads the UK FTSE (Financial Times Stock Exchange) index, which you can then refer to as object FTSE.

If you want to extract the data from a certain date range, for example from April 1 to June 30, 2000, you can create a new object using the command ftse00 <- FTSE[“2000-04-01/2000-06-30”].

library(zoo)

package 㤼㸱zoo㤼㸲 was built under R version 4.0.3
Attaching package: 㤼㸱zoo㤼㸲

The following object is masked from 㤼㸱package:timeSeries㤼㸲:

    time<-

The following objects are masked from 㤼㸱package:base㤼㸲:

    as.Date, as.Date.numeric

Warning messages:
1: package ‘timeSeries’ was built under R version 4.0.3 
2: package ‘timeDate’ was built under R version 4.0.3 

library(xts)

package 㤼㸱xts㤼㸲 was built under R version 4.0.3

library(qrmdata)

package 㤼㸱qrmdata㤼㸲 was built under R version 4.0.3

# Load DJ index
data("DJ")

# Show head() and tail() of DJ index
head(DJ)

              ^DJI
1985-01-29 1292.62
1985-01-30 1287.88
1985-01-31 1286.77
1985-02-01 1277.72
1985-02-04 1290.08
1985-02-05 1285.23

tail(DJ)

               ^DJI
2015-12-23 17602.61
2015-12-24 17552.17
2015-12-28 17528.27
2015-12-29 17720.98
2015-12-30 17603.87
2015-12-31 17425.03

# Plot DJ index
plot(DJ)

# Extract 2008-2009 and assign to dj0809
dj0809 <- DJ["2008-01-01/2009-12-31"]

# Plot dj0809
plot(dj0809)

Take a good look at how it behaves and note how far the index fell in the 2008-2009 financial crisis.

1.1.2 Exploring risk-factor time series: individual equities

For some risk management applications, it is sufficient to model equity risk by looking at indexes. If you want a more detailed model of the risk in a portfolio of equities, you can drill down to the level of individual share prices.

In the previous chapter, you used DJ[“2008/2009”] to extract the Dow Jones data from certain rows of an xts object by specifying a date range index. To also extract data from particular columns, you can add a column identifier, like a string name or numeric index, in the brackets following a comma. To select multiple columns, include these column identifiers in a vector. This [rows, columns] format is consistent with indexing most other two dimensional objects in R.

data[index, colname]
data[index, c(col1index, col2index)]

The qrmdata package also includes data for certain constituents, or the stocks or companies part of a larger index. The Dow Jones constituents data are contained in “DJ_const”. In this exercise, you will view the names of all its stocks, select the Apple and Goldman Sachs share prices, and plot them using the command plot.zoo() to display multiple time series.

# Load DJ constituents data
data("DJ_const")

# Apply names() and head() to DJ_const
names(DJ_const)

 [1] "AAPL" "AXP"  "BA"   "CAT"  "CSCO" "CVX"  "DD"   "DIS"  "GE"   "GS"   "HD"   "IBM"  "INTC" "JNJ" 
[15] "JPM"  "KO"   "MCD"  "MMM"  "MRK"  "MSFT" "NKE"  "PFE"  "PG"   "TRV"  "UNH"  "UTX"  "V"    "VZ"  
[29] "WMT"  "XOM" 

head(DJ_const)

           AAPL AXP       BA      CAT CSCO CVX       DD      DIS       GE GS HD      IBM INTC JNJ JPM
1962-01-02   NA  NA 0.212905 0.593184   NA  NA 1.227958 0.061014 0.145967 NA NA 2.346625   NA  NA  NA
1962-01-03   NA  NA 0.217163 0.598964   NA  NA 1.229230 0.061832 0.144501 NA NA 2.367136   NA  NA  NA
1962-01-04   NA  NA 0.215034 0.614372   NA  NA 1.220331 0.061832 0.142794 NA NA 2.343548   NA  NA  NA
1962-01-05   NA  NA 0.210773 0.620152   NA  NA 1.187280 0.062039 0.139132 NA NA 2.297395   NA  NA  NA
1962-01-08   NA  NA 0.211306 0.624001   NA  NA 1.169484 0.061832 0.138889 NA NA 2.254319   NA  NA  NA
1962-01-09   NA  NA 0.211839 0.629777   NA  NA 1.173297 0.063060 0.139620 NA NA 2.280983   NA  NA  NA
                 KO MCD MMM MRK MSFT NKE PFE PG TRV UNH UTX  V VZ WMT XOM
1962-01-02 0.031031  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA
1962-01-03 0.030339  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA
1962-01-04 0.030571  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA
1962-01-05 0.029879  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA
1962-01-08 0.029570  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA
1962-01-09 0.030110  NA  NA  NA   NA  NA  NA NA  NA  NA  NA NA NA  NA  NA

# Extract AAPL and GS in 2008-09 and assign to stocks
stocks <- DJ_const["2008/2009", c("AAPL", "GS")]

# Plot stocks with plot.zoo()
plot.zoo(stocks)

1.1.3 Exploring risk-factor data: exchange rates

For a portfolio with risk exposure in different countries, it is necessary to consider the risk coming from foreign exchange (FX) rates. The qrmdata package includes FX rate data for many currencies, ranging from Swiss Francs to Japanese Yen, with respect to the USD (United States dollar) and GBP (Great Britain pound).

In this exercise, you will look at the datasets “EUR_USD” and “GBP_USD”, which contain the Euro and British pound exchange rates against the US dollar. Then, you will merge these time series and plot them together for the period 2010-2015.

# Load exchange rate data
data("GBP_USD")
data("EUR_USD")

# Plot the two exchange rates
plot(GBP_USD)

plot(EUR_USD)

# Plot a USD_GBP exchange rate
plot(1/GBP_USD)

# Merge the two exchange rates GBP_USD and EUR_USD
fx <- merge(GBP_USD, EUR_USD, all = TRUE)

# Extract 2010-15 data from fx and assign to fx0015
fx0015 <- fx["2010/2015"]

# Plot the exchange rates in fx0015
plot.zoo(fx0015)

Note that merging the EUR_USD and GBP_USD data, in that order, would have produced a different object from fx.

1.2 Risk-factor returns

1.2.1 Exploring return series

To analyze risk, the key task is to model the fluctuations in prices and rates over different time periods; these fluctuations are known as returns. To calculate the log-returns of the FTSE stock index and assign to ftse_x, apply the log() and diff() functions in succession:

ftse_x <- diff(log(FTSE))

differencing in this way will always give a NA in the first position of the time series, which can then be removed with diff(log(FTSE))[-1]. However, you will not need to do this in the course unless it is specified in the instructions.

In this exercise, you will calculate and plot log-return series for the equity and FX risk factors that you have previously encountered.

# Compute the log-returns of dj0809 and assign to dj0809_x
dj0809_x <- diff(log(dj0809))

# Plot the log-returns
plot(dj0809_x)

# Compute the log-returns of djstocks and assign to djstocks_x
djstocks <- DJ_const["2008/2009", c("AAPL", "GS")]
djstocks_x <- diff(log(djstocks))

# Plot the two share returns
plot.zoo(djstocks_x)

# Compute the log-returns of GBP_USD and assign to erate_x
erate_x <- diff(log(GBP_USD))

# Plot the log-returns
plot(erate_x)

The return series often just look like noise with some periods of larger fluctuations. You’ll discover later that they typically have a lot of interesting structure.

1.2.2 Different ways of plotting risk-factor and return series

You already know that you can use plot.zoo() to plot multiple time series. For a four-dimensional time series data, the call plot.zoo(data) creates four separate plots by default, unless you include the parameter plot.type = “single” to plot all four series in one plot. You can also add even more parameters such as col to specify different colors and type = “h” to get vertical bars instead of joining points, which can sometimes be a better way of displaying returns.

plot.zoo(x, plot.type, col = 1, type = "l", ...)

In this exercise, you will explore the plot.zoo() function to plot equity risk-factor data and the corresponding returns in different ways.

djstocks <- DJ_const["2008/2009", c("AAPL", "AXP", "BA", "CAT")]
# Plot djstocks in four separate plots
plot.zoo(djstocks)

# Plot djstocks in one plot and add legend
plot.zoo(djstocks, plot.type = "single", col = c(1:4))
legend(julian(x = as.Date("2009-01-01")), y = 70, legend = names(DJ_const)[1:4], fill = 1:4)

# Compute log-returns and assign to djstocks_x
djstocks_x <- diff(log(djstocks))

# Plot djstocks_x in four separate plots
plot.zoo(djstocks_x)

# Plot djstocks_x with vertical bars
plot.zoo(djstocks_x, type = "h")

Note how in late 2008 there were large returns for all series. That was the height of the financial crisis.

1.3 Aggregating log-returns

1.3.1 Aggregating log-return series

In statistics, aggregate data are data combined from several measurements. You just learned that you can compute compute weekly, monthly and quarterly log-returns by summing daily log-returns with the corresponding apply.weekly(), apply.monthly() and apply.quarterly() functions.

For example, you can use the following code to form the quarterly returns for a univariate time series data and multivariate time series mv_data:

# apply.quarterly(x, FUN, ...)
data_q = apply.quarterly(data, sum)
mv_data_q = apply.quarterly(mv_data, colSums)

In this exercise, you will practice aggregating time series data using these functions and plotting the results. The data DJ and DJ_const are available in your workspace, as are the objects djx, which contains daily log-returns of the Dow Jones index from 2000-2015, and djreturns, which contains the daily log-returns for the first four DJ_const stocks from 2000-2015. Use plot for univariate time series and plot.zoo for multivariate time series.

dj <- DJ["2000/2015"]
djx <- diff(log(dj))
djstocks <- DJ_const["2000/2015",c("AAPL", "AXP", "BA", "CAT")]
djreturns <- diff(log(djstocks))

# Plot djx
plot(djx)

# Plot weekly log-returns of djx
plot(apply.weekly(djx, sum), type = "h")

# Plot monthly log-returns of djx
plot(apply.monthly(djx, sum), type = "h")

# Plot djreturns
plot.zoo(djreturns)

# Plot monthly log-returns of djreturns
plot.zoo(apply.monthly(djreturns, colSums), type = "h")

These aggregation functions are extremely useful as for analyzing risk over longer time horizons.

Data scientists often use the aggregations that you have learned so far in combination with summary statistics to extract even more insights from data. Functions that calculate summary statistics include mean(), median(), and var(). The object sp contains daily log-returns for the S&P 500 index for the period 1960-2015; it is loaded in your workspace. To three decimal places, what is the average quarterly log-return for the S&P 500 from 1990-2010?

data("SP500")
sp <- SP500["1990/2010"]
sp <- diff(log(sp))[-1]
mean(apply.quarterly(sp, sum))

[1] 0.01490178

1.4 Exploring other kinds of risk factors

1.4.1 Commodities data

The plotting function pairs() creates a pairwise scatterplot of the components of a multivariate time series with two or more dimensions. It is used on a zoo object rather than an xts object.

A roughly circular shape of a scatterplot indicates a low correlation between the log-returns of two different commodities. Generally speaking, low correlation is good in a portfolio as it implies that the assets are diversified. High correlation, on the other hand, represents a risk that must be properly modelled.

In this exercise, you will look at gold and oil prices over a 25 year period, calculate their daily and monthly log-returns, and plot them. The data gold and oil, containing the daily prices from 1990-2015 of gold and Brent crude oil.

data("OIL_Brent")
oil <- OIL_Brent["1990/2015"]
data("GOLD")
gold <- GOLD["1990/2015"]

# Plot gold and oil prices
plot(gold)

plot(oil)

# Calculate daily log-returns
goldx <- diff(log(gold))
oilx <- diff(log(oil))

# Calculate monthly log-returns
goldx_m <- apply.monthly(goldx, sum)
oilx_m <- apply.monthly(oilx, sum)

# Merge goldx_m and oilx_m into coms
coms <- merge(goldx_m, oilx_m)

# Plot coms with vertical bars
plot.zoo(coms, type = "h")

# Make a pairwise scatterplot of coms
pairs(as.zoo(coms))

As you can see, gold and oil are well diversified commodities.

1.5 Interest-rate data

The object zcb contains daily values of Canadian zero-coupon-bond yields, expressed as percentages, for the period 2006-2015. Yields are the key risk-factor when it comes to analysing the interest-rate risk in a portfolio of bonds or other fixed-income products.

It is not so clear what is the best way of calculating risk-factor changes for yields. It is possible to compute log-returns, provided yields are not negative, and it is also possible to calculate simple returns. To compute the simple returns of a series, use only diff() instead of diff() and log().

In this exercise, you will plot time series of yields for fixed times to maturity, and plot risk-factor changes for these yields. You will also plot the whole yield curve on particular dates. The zcb data has been loaded into your workspace. A vector yield_cols containing the names of the columns corresponding to maturities of 1, 5 and 10 years has been created. A numerical vector maturity containing all the maturities in years has also been created.

data("ZCB_CAD")
zcb <- ZCB_CAD["2006/2015"]

maturity <- c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 
3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 5.75, 6, 6.25, 
6.5, 6.75, 7, 7.25, 7.5, 7.75, 8, 8.25, 8.5, 8.75, 9, 9.25, 9.5, 
9.75, 10, 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25, 
12.5, 12.75, 13, 13.25, 13.5, 13.75, 14, 14.25, 14.5, 14.75, 
15, 15.25, 15.5, 15.75, 16, 16.25, 16.5, 16.75, 17, 17.25, 17.5, 
17.75, 18, 18.25, 18.5, 18.75, 19, 19.25, 19.5, 19.75, 20, 20.25, 
20.5, 20.75, 21, 21.25, 21.5, 21.75, 22, 22.25, 22.5, 22.75, 
23, 23.25, 23.5, 23.75, 24, 24.25, 24.5, 24.75, 25, 25.25, 25.5, 
25.75, 26, 26.25, 26.5, 26.75, 27, 27.25, 27.5, 27.75, 28, 28.25, 
28.5, 28.75, 29, 29.25, 29.5, 29.75, 30)

yield_cols <- c("1.00y", "5.00y", "10.00y")
# Compute log-returns as zcb_x and simple returns as zcb_x2
zcb_x <- diff(log(zcb))
zcb_x2 <- diff(zcb)

# Plot zcb_x for 1, 5 and 10-year maturities
plot.zoo(zcb_x[, yield_cols])

# Plot zcb_x2 for 1, 5 and 10-year maturities
plot.zoo(zcb_x2[, yield_cols])

# Plot the yield curve for the first day of zcb
plot(maturity, zcb[1,], ylim = range(zcb), type = "l", ylab = "yield (%)", col = "red")

# Add a line for the last day of zcb
lines(maturity, zcb[nrow(zcb), ])

Yields can seem a little tricky to work with at first but they are fundamental for analyzing bond portfolios and many other financial products that depend on interest rates.

2 Real world returns are riskier than normal

2.1 The normal distribution

2.1.1 Graphical methods for assessing normality

To create a histogram with 20 buckets that represents the probability density of the FTSE data, as well as how to add a normal distribution to the existing plot as a red line:

hist(ftse, nclass = 20, probability = TRUE)
lines(ftse, dnorm(ftse, mean = mu, sd = sigma), col = "red")

As you can see, dnorm(x, mean, sd) calculates the probability density function (PDF) of the data x with the calculated sample mean and standard deviation; this is known as the method-of-moments.

Finally, to calculate an estimate of the density of data x, use density(x). This creates a so-called kernel-density estimate (KDE) using a non-parametric method that makes no assumptions about the underlying distribution.

The various plots suggest that the data are heavier tailed than normal, although you will learn about better graphical and numerical tests in future exercises.

In this exercise, you will fit a normal distribution to the log-returns of the Dow Jones index for 2008-2009 and compare the data with the fitted distribution using a histogram and a density plot. The object djx containing Dow Jones data is loaded into your workspace.

djx <- DJ["2008/2009"]
djx <- diff(log(djx))
djx <- as.numeric(djx)[-1]
djx <- sort(djx)

# Calculate average and standard deviation of djx
mu <- mean(djx)
sigma <- sd(djx)

# Plot histogram of djx
hist(djx, nclass = 20, probability = TRUE)

# Add the normal density as a red line to histogram
lines(djx, dnorm(djx, mean = mu, sd = sigma), col = "red")

# Plot non-parametric KDE of djx
plot(density(djx))

# Add the normal density as red line to KDE
lines(djx, dnorm(djx, mean = mu, sd = sigma), col = "red")

mu <- mean(djx)
sigma <- sd(djx)
te <- dnorm(djx, mean = mu, sd = sigma)
ord <- sort(djx)
head(ord)

[1] -0.08200514 -0.08014005 -0.07615925 -0.07234497 -0.05863401 -0.05725062

tail(ord)

[1] 0.05633869 0.06337643 0.06458521 0.06611575 0.10325919 0.10508346

te <- dnorm(ord, mean = mu, sd = sigma)
head(te)

[1] 0.004951580 0.007206443 0.015594272 0.031482814 0.291488669 0.355490098

The data don’t look very normal. Compare in particular the center and the tails of the histogram and density plot with the red normal curve.

2.2 Testing for normality

2.2.1 Q-Q plots for assessing normality

The quantile-quantile plot (Q-Q plot) is a better graphical method for revealing non-normality. In general, a Q-Q plot compares the quantiles of the data with the quantiles of a reference distribution; if the data are from a distribution of the same type (up to scaling and location), a reasonably straight line should be observed. You should know that the degrees of freedom (df) refer to the number of values or observations that can affect the system you are working with.

To generate 1000 normal data points with the rnorm() function, as well as how to use qqnorm() to create the Q-Q plot, and qqline() to add a straight line for reference:

data <- rnorm(1000, mean = 3, sd = 2)
qqnorm(data)
qqline(data)

You will create a Q-Q plot of the Dow Jones log-returns in djx against the normal reference distribution, which you will add as a visual guide. You will then compare the plot with simulated datasets from normal, Student t and uniform distributions generated with the rnorm(), rt() and runif() functions. Y

If the data are from a normal distribution the dots should be close to the red line (although there may be some deviation at the very end).

# Make a Q-Q plot of djx and add a red line
qqnorm(djx)
qqline(djx, col = "red")

# Calculate the length of djx as n
n <- length(djx)

# Generate n standard normal variables, make a Q-Q plot, add a red line
x1 <- rnorm(n)
qqnorm(x1)
qqline(x1, col = "red")

# Generate n Student t variables, make a Q-Q plot, add a red line
x2 <- rt(n, df = 4)
qqnorm(x2)
qqline(x2, col = "red")

# Generate n standard uniform variables, make a Q-Q plot, add red line
x3 <- runif(n)
qqnorm(x3)
qqline(x3, col = "red")

The Q-Q plot is a very effective tool that is widely used in applied statistical and econometric work.

2.3 Skewness, kurtosis and the Jarque-Bera test

2.3.1 Numerical tests of normality

The moments package contains functions for computing the kurtosis and skewness of data and well as for implementing the Jarque-Bera test, which is a test of normality based on these higher-order moments. In one command, it compares the skewness and kurtosis of the data with the theoretical values for the normal distribution, which are 0 and 3, respectively.

jarque.test(x)
skewness(x, na.rm = FALSE)
kurtosis(x, na.rm = FALSE)

In this exercise, you will calculate the skewness and kurtosis for the djx, the Dow Jones index from 2008-2011, and apply the Jarque-Bera test of normality. You will then apply the same methods to djreturns, which contains 29 of the Dow Jones stocks for the same period.

Recall that you can use apply(X, MARGIN, FUN, …) to apply functions over array margins. The MARGIN parameter is a vector indicating where the function will be applied; in this instance, you will use 2 to specify that the function FUN should be applied to the columns in matrix X.

dj <- DJ["2008/2011"]
djx <- diff(log(dj))
djx <- as.numeric(djx)[-1]
djx <- sort(djx)

djstocks <- DJ_const["2008/2011"]
djreturns <- diff(log(djstocks))[-1]

library(moments)

package 㤼㸱moments㤼㸲 was built under R version 4.0.3
Attaching package: 㤼㸱moments㤼㸲

The following objects are masked from 㤼㸱package:timeDate㤼㸲:

    kurtosis, skewness

# Calculate skewness and kurtosis of djx
skewness(djx)

[1] 0.001913285

kurtosis(djx)

[1] 9.058875

# Carry out a Jarque-Bera test for djx
jarque.test(djx)

    Jarque-Bera Normality Test

data:  djx
JB = 1541.8, p-value < 2.2e-16
alternative hypothesis: greater

# Calculate skewness and kurtosis of djreturns 
s <- apply(djreturns, 2, skewness)
k <- apply(djreturns, 2, kurtosis)

# Plot k against s and add text labels to identify stocks
plot(s, k, type = "n")
text(s, k, names(s), cex = 0.6)

# Carry out Jarque-Bera tests for each constituent in djreturns
apply(djreturns, 2, jarque.test)

$AAPL

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1819.5, p-value < 2.2e-16
alternative hypothesis: greater


$AXP

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1066.1, p-value < 2.2e-16
alternative hypothesis: greater


$BA

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 347.08, p-value < 2.2e-16
alternative hypothesis: greater


$CAT

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 304.64, p-value < 2.2e-16
alternative hypothesis: greater


$CSCO

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2758.9, p-value < 2.2e-16
alternative hypothesis: greater


$CVX

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5268.2, p-value < 2.2e-16
alternative hypothesis: greater


$DD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 459.61, p-value < 2.2e-16
alternative hypothesis: greater


$DIS

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1053.4, p-value < 2.2e-16
alternative hypothesis: greater


$GE

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1176.9, p-value < 2.2e-16
alternative hypothesis: greater


$GS

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 3513.6, p-value < 2.2e-16
alternative hypothesis: greater


$HD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 470.66, p-value < 2.2e-16
alternative hypothesis: greater


$IBM

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 660.46, p-value < 2.2e-16
alternative hypothesis: greater


$INTC

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 510.51, p-value < 2.2e-16
alternative hypothesis: greater


$JNJ

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5987.3, p-value < 2.2e-16
alternative hypothesis: greater


$JPM

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1906.3, p-value < 2.2e-16
alternative hypothesis: greater


$KO

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 4310.6, p-value < 2.2e-16
alternative hypothesis: greater


$MCD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 809.57, p-value < 2.2e-16
alternative hypothesis: greater


$MMM

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 379.55, p-value < 2.2e-16
alternative hypothesis: greater


$MRK

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2647.7, p-value < 2.2e-16
alternative hypothesis: greater


$MSFT

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2280.7, p-value < 2.2e-16
alternative hypothesis: greater


$NKE

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 944.42, p-value < 2.2e-16
alternative hypothesis: greater


$PFE

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 717.42, p-value < 2.2e-16
alternative hypothesis: greater


$PG

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1690.9, p-value < 2.2e-16
alternative hypothesis: greater


$TRV

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 8039.2, p-value < 2.2e-16
alternative hypothesis: greater


$UNH

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 10514, p-value < 2.2e-16
alternative hypothesis: greater


$UTX

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1044, p-value < 2.2e-16
alternative hypothesis: greater


$V

    Jarque-Bera Normality Test

data:  newX[, i]
JB = NA, p-value = NA
alternative hypothesis: greater


$VZ

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1969.1, p-value < 2.2e-16
alternative hypothesis: greater


$WMT

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2350.8, p-value < 2.2e-16
alternative hypothesis: greater


$XOM

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5812.4, p-value < 2.2e-16
alternative hypothesis: greater

The return distributions of the Dow Jones stocks all have high kurtosis and some of them are quite skewed.

2.3.2 Testing normality for longer time horizons

As returns are added together over longer time periods, a central limit effect takes place and returns tend to become more normal.

In this exercise, you will use aggregation functions that you learned in the first chapter to aggregate the data in djx_d, containing the daily log-returns for 29 of the Dow Jones stocks for the period 2000-2015.

djstocks <- DJ_const["2000/2015"]
djreturns <- diff(log(djstocks))[-1]
djx_d <- djreturns
# Calculate weekly and monthly log-returns from djx_d
djx_w <- apply.weekly(djx_d, colSums)
djx_m <- apply.monthly(djx_d, colSums)

# Calculate the p-value for each series in djx_d
apply(djx_d, 2, function(v){jarque.test(v)$p.value})

AAPL  AXP   BA  CAT CSCO  CVX   DD  DIS   GE   GS   HD  IBM INTC  JNJ  JPM   KO  MCD  MMM  MRK MSFT  NKE 
   0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
 PFE   PG  TRV  UNH  UTX    V   VZ  WMT  XOM 
   0    0    0    0    0   NA    0    0    0 

# Calculate the p-value for each series in djx_w
apply(djx_w, 2, function(v){jarque.test(v)$p.value})

AAPL  AXP   BA  CAT CSCO  CVX   DD  DIS   GE   GS   HD  IBM INTC  JNJ  JPM   KO  MCD  MMM  MRK MSFT  NKE 
   0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
 PFE   PG  TRV  UNH  UTX    V   VZ  WMT  XOM 
   0    0    0    0    0   NA    0    0    0 

# Calculate the p-value for each series in djx_m
apply(djx_m, 2, function(v){jarque.test(v)$p.value})

        AAPL          AXP           BA          CAT         CSCO          CVX           DD          DIS 
0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.110223e-16 1.677566e-01 5.059611e-03 4.043754e-11 
          GE           GS           HD          IBM         INTC          JNJ          JPM           KO 
3.086420e-13 3.259324e-02 5.493080e-02 0.000000e+00 0.000000e+00 5.472560e-08 1.982609e-06 5.641173e-07 
         MCD          MMM          MRK         MSFT          NKE          PFE           PG          TRV 
0.000000e+00 9.341348e-02 1.524618e-07 0.000000e+00 0.000000e+00 3.076613e-01 0.000000e+00 0.000000e+00 
         UNH          UTX            V           VZ          WMT          XOM 
0.000000e+00 0.000000e+00           NA 0.000000e+00 7.281104e-04 1.100037e-02 

Although the p-values get larger, all monthly returns other than for Chevron (CVX), 3M (MMM), and Pfizer (PFE) still fail the normality test.

2.3.3 Overlapping returns

When you aggregate series by summing daily log-returns into longer intervals, you analyze a smaller amount of observations. To preserve the quantity of data, you can calculate overlapping returns with the rollapplyr() function; this also creates strong correlations between observations.

There are 5 trading days in the average calendar week. By computing the 5-day moving sums of the log-returns of daily index data, you obtain approximate overlapping weekly returns ending on each calendar week. Similarly, calculating 21-day moving sums gives approximate overlapping monthly returns, and calculating 63-day moving sums gives approximate overlapping quarterly returns. Let’s look at an example with the Dow Jones daily return data in djx. Because 5 values are used to calculate each moving sum, the first 4 values in the result are NA. In this instance, we will use indexing to remove them:

dj <- DJ["2008/2011"]
djx <- diff(log(dj))[-1]
djx5 <- rollapplyr(djx, width = 5, FUN = sum)
head(djx5)

                  ^DJI
2008-01-03          NA
2008-01-04          NA
2008-01-07          NA
2008-01-08          NA
2008-01-09 -0.02394677
2008-01-10 -0.01571869

djx5 <- djx5[-(1:4)]

n this exercise, you will calculate moving sums of different intervals from djx, which is loaded in your workspace. You will then find the skewness and kurtosis of the resulting data and conduct the Jarque-Bera test just as you have in previous exercises. Do the overlapping returns appear more normal?

# Calculate a 21-day moving sum of djx
djx21 <- rollapplyr(djx, width = 21, FUN = sum)[-(1:20)]

# Calculate a 63-day moving sum of djx
djx63 <- rollapplyr(djx, width = 63, FUN = sum)[-(1:62)]

# Merge the three series and plot
djx2 <- merge(djx, djx21, djx63, all = FALSE)
plot.zoo(djx2)

# Compute the skewness and kurtosis for each series in djx2
apply(djx2, 2, skewness)

      X.DJI     X.DJI.1     X.DJI.2 
-0.01276806 -1.16018080 -0.87509434 

apply(djx2, 2, kurtosis)

   X.DJI  X.DJI.1  X.DJI.2 
9.264838 5.845321 3.845090 

# Conduct the Jarque-Bera test to each series in djx2
apply(djx2, 2, jarque.test)

$X.DJI

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 1547.1, p-value < 2.2e-16
alternative hypothesis: greater


$X.DJI.1

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 531.33, p-value < 2.2e-16
alternative hypothesis: greater


$X.DJI.2

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 148.89, p-value < 2.2e-16
alternative hypothesis: greater

These overlapping returns are highly correlated and even more difficult to interpret.

2.4 The Student t distribution

2.4.1 Fitting t distribution to data

A Student t distribution is generally a much better fit to daily, weekly, and monthly returns than a normal distribution.

You can create one by using the fit.st() function in the QRM package. The resulting fitted model has a parameter estimates component par.ests which can be assigned to a list tpars in order to store its values of nu, mu, and sigma for later use:

tfit <- fit.st(ftse)
tpars <- tfit$par.ests
tpars

You will fit a Student t distribution to the daily log-returns of the Dow Jones index from 2008-2011 contained in djx. Then, you will plot a histogram of the data and superimpose a red line to the plot showing the fitted t density.

library(QRM)

djx <- as.numeric(djx)
djx <- sort(djx)
# Fit a Student t distribution to djx
tfit <- fit.st(djx)

# Define tpars, nu, mu, and sigma
tpars <- tfit$par.ests
nu <- tpars[1]
mu <- tpars[2]
sigma <- tpars[3]

# Plot a histogram of djx
hist(djx, nclass = 20, probability = TRUE, ylim = range(0, 40))

# Compute the fitted t density at the values djx
yvals <- dt((djx - mu)/sigma, df = nu)/sigma

# Superimpose a red line to show the fitted t density
lines(djx, yvals, col = "red")

The fitted Student t distribution looks a lot better than the normal did.

2.4.2 Testing FX returns for normality

So far, the exercises in this chapter have examined the normality of equity index returns and individual equity returns.

To reinforce these ideas, you will apply similar ideas to exchange-rate log-returns. The dataset fx_d contains daily log-returns of the EUR/USD, GBP/USD and JPY/USD exchange rates for the period 2001-2015, and the dataset fx_m contains the corresponding monthly log-returns. Which of the monthly log-return series appears the most normal?

data("EUR_USD")
data("GBP_USD")
data("JPY_USD")
fx1 <- EUR_USD["2001/2015"]
fx2 <- GBP_USD["2001/2015"]
fx3 <- JPY_USD["2001/2015"]
fx <- merge(fx1, fx2, fx3)
fx_d <- apply(log(fx), 2, diff)
fx_m <- apply.monthly(fx_d, colSums)

# Plot the daily log-return series in fx_d
plot.zoo(fx_d, type = "h")

# Apply the Jarque-Bera test to each of the series in fx_d
apply(fx_d, 2, jarque.test)

$EUR.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2791.1, p-value < 2.2e-16
alternative hypothesis: greater


$GBP.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 10962, p-value < 2.2e-16
alternative hypothesis: greater


$JPY.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5039.7, p-value < 2.2e-16
alternative hypothesis: greater

# Plot the monthly log-return series in fx_m
plot.zoo(fx_m, type = "h")

# Apply the Jarque-Bera test to each of the series in fx_m
apply(fx_m, 2, jarque.test)

$EUR.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 22.717, p-value = 1.167e-05
alternative hypothesis: greater


$GBP.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 40.814, p-value = 1.372e-09
alternative hypothesis: greater


$JPY.USD

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 3.6561, p-value = 0.1607
alternative hypothesis: greater

# Fit a Student t distribution to each of the series in fx_m
apply(fx_m, 2, function(v){fit.st(v)$par.ests})

         EUR.USD      GBP.USD      JPY.USD
nu    4.91304346 6.0621630186 1.840870e+01
mu    0.00194676 0.0008555348 1.118179e-04
sigma 0.02376514 0.0204503515 2.553159e-02

JPY/USD exchange rate log-returns appear to be the most normal.

2.4.3 Testing interest-rate returns for normality

The object zcbx_m contains monthly log-return series for the 1-year, 5-year and 10-year Canadian zero-coupon bond yields. The object zcbx2_m contains the corresponding simple returns. Both are multivariate; they are loaded into your workspace.

In this exercise, you will plot these interest rate return series and then examine their normality with Q-Q plots and Jarque-Bera tests.

The log-returns show clearer evidence of non-normality than the simple returns in this case.

data("ZCB_CAD")
zcb <- ZCB_CAD[,c("1.00y", "5.00y", "10.00y")]
zcbx <- diff(log(zcb))
zcbx2 <- diff(zcb)
zcbx_m <- apply.monthly(zcbx, colSums)["2006/2015"]
zcbx2_m<- apply.monthly(zcbx2, colSums)["2006/2015"]

# Plot the interest-rate return series zcbx_m and zcbx2_m
plot.zoo(zcbx_m, type = "h")

plot.zoo(zcbx2_m, type = "h")

# Make Q-Q plots of the 3rd component series of zcbx_m and zcbx2_m
qqnorm(zcbx_m[,3])

qqnorm(zcbx2_m[,3])

# Compute the kurtosis of each series in zcbx_m and zcbx2_m
apply(zcbx_m, 2, kurtosis)

    1.00y     5.00y    10.00y 
 8.597937 13.899820  7.990000 

apply(zcbx2_m, 2, kurtosis)

   1.00y    5.00y   10.00y 
6.105796 3.737785 3.510665 

# Conduct the Jarque-Bera test on each series in zcbx_m and zcbx2_m
apply(zcbx_m, 2, jarque.test)

$`1.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 164.94, p-value < 2.2e-16
alternative hypothesis: greater


$`5.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 645.07, p-value < 2.2e-16
alternative hypothesis: greater


$`10.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 142.19, p-value < 2.2e-16
alternative hypothesis: greater

apply(zcbx2_m, 2, jarque.test)

$`1.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 77.413, p-value < 2.2e-16
alternative hypothesis: greater


$`5.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5.7791, p-value = 0.0556
alternative hypothesis: greater


$`10.00y`

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 4.0851, p-value = 0.1297
alternative hypothesis: greater

Note how the simple monthly returns for the 5 and 10 year yields did not fail the normality test.

2.4.4 Testing gold price returns for normality

The object goldx_q contains quarterly log-returns of the gold price from the beginning of 1990 to the end of 2015.

Test the data for normality using the Jarque-Bera test, then fit a Student t distribution and find the estimated degree of freedom ν^ to the nearest integer.

data("GOLD")
gold <- GOLD["1990/2015"]
goldx <- diff(log(gold))[-1]
goldx_q <- as.numeric(apply.quarterly(goldx, sum)[,1])
jarque.test(goldx_q)

    Jarque-Bera Normality Test

data:  goldx_q
JB = 35.01, p-value = 2.498e-08
alternative hypothesis: greater

tfit <- fit.st(goldx_q)

tpars <- tfit$par.ests

tpars[1]

      nu 
9.917565 

The data fail a test of normality and \(\hat{\nu}=10\)

3 Real world returns are volatile and correlated

3.1 Characteristics of volatile return series

In this exercise you will plot the Dow Jones log-returns for 2008-2011 alongside independent and identically distributed (iid) normal data and iid Student t data.

The xts object djx contains the Dow Jones index and the objects npars and tpars contain the parameter estimates that are obtained when a normal distribution and a t distribution are fitted to djx. you will generate a normal sample from the fitted model by generating standard normal data, scaling them with the second component npars[2] and shifting them with the first component npars[1]. For the Student t sample you will do something similar but note that this time the first component tpars[1] contains the degree of freedom parameter and tpars[2] and tpars[3] contain the location and scale parameters.

After making the plots, you should compare the behavior of the real returns with that of the iid returns, particularly around the 2008 financial crisis.

djx <- DJ["2007-12-31/2011-12-31"]
djx <- diff(log(djx))[-1]
npars <- c(mu = mean(djx), sigma = sd(djx))
tfit <- fit.st(djx)
tpars <- tfit$par.ests

set.seed(123)
# Compute the length n of djx 
n <- length(djx)

#  Generate a normal sample of size n with parameters given by npars
ndata <- rnorm(n)*npars[2] + npars[1]

# Generate a t-distributed sample of size n with paramaters given by tpars
tdata <- rt(n, df = tpars[1])*tpars[3] + tpars[2]

# Make ndata and tdata into xts objects
ndatax <- xts(ndata, time(djx))
tdatax <- xts(tdata, time(djx))

# Merge djx, ndatax, and tdatax and plot
alldata <- merge(djx, ndatax, tdatax)
plot.zoo(alldata, type = "h", ylim = range(alldata))

Note how the real returns look more volatile than the iid returns. There was clearly a period of large movements around the 2008 financial crisis followed by a quieter period in 2010.

3.2 Estimating serial correlations

3.2.1 Using acf plots to reveal volatility

This exercise follows on from the previous R exercise where we looked for visible signs of volatility in a financial time series. For the Dow Jones returns from 2008-11 in djx and the simulated normal and t-distributed data in ndata and tdata, respectively, you will calculate and plot the sample autocorrelation functions (acf) using the command acf().

While very little evidence of serial correlation is found in these plots, the picture changes dramatically when we look at absolute or squared return data. The real returns in the Dow Jones return series djx behave very differently to the simulated data. The serial correlation in absolute or squared returns is a consequence of volatility, which causes large returns to be followed by further large returns, although not necessarily of the same sign.

# Set up a plot region to show 3 plots at a time
par(mfrow = c(3, 1))

# Plot the acfs of djx, ndata and tdata
acf(djx)
acf(ndata)

acf(tdata)

# Plot the acfs of the absolute values
acf(abs(djx))
acf(abs(ndata))

acf(abs(tdata))

# Plot the acfs of the squares of the values
acf(djx^2)
acf(ndata^2)

acf(tdata^2)

It is very common to see strong serial correlation in absolute and squared return series.

3.3 The Ljung-Box test

3.3.1 Applying Ljung-Box tests to return data

This code applies the Ljung-Box test to the ftse data with a lag of 10:

Box.test(ftse, lag = 10, type = "Ljung")

You will carry out a Ljung-Box test for serial correlation on the time series djx which contains the Dow Jones daily index returns for 2008-2011, as well as on all the individual equity return series in djall which contains the Dow Jones data for 2006-2015. You will implement this test on both the raw return series and the absolute values of the series, which you can calculate with abs().

You should notice that while the hypothesis of no serial correlation is rejected for many of the raw return series, it is rejected overwhelmingly for all of the absolute value series.

djall <- DJ_const["2006-12-31/2015-12-31"]
djall <- diff(log(djall))[-1]

# Apply the Ljung-Box test to djx
Box.test(djx, lag = 10, type = "Ljung")

    Box-Ljung test

data:  djx
X-squared = 34.228, df = 10, p-value = 0.000169

# Apply the Ljung-Box test to absolute values of djx
Box.test(abs(djx), lag = 10, type = "Ljung")

    Box-Ljung test

data:  abs(djx)
X-squared = 1083.9, df = 10, p-value < 2.2e-16

# Apply the Ljung-Box test to all return series in djall
apply(djall, 2, Box.test, lag = 10, type = "Ljung")

$AAPL

    Box-Ljung test

data:  newX[, i]
X-squared = 15.063, df = 10, p-value = 0.1298


$AXP

    Box-Ljung test

data:  newX[, i]
X-squared = 45.453, df = 10, p-value = 1.802e-06


$BA

    Box-Ljung test

data:  newX[, i]
X-squared = 19.929, df = 10, p-value = 0.02993


$CAT

    Box-Ljung test

data:  newX[, i]
X-squared = 16.596, df = 10, p-value = 0.08379


$CSCO

    Box-Ljung test

data:  newX[, i]
X-squared = 15.694, df = 10, p-value = 0.1087


$CVX

    Box-Ljung test

data:  newX[, i]
X-squared = 63.987, df = 10, p-value = 6.331e-10


$DD

    Box-Ljung test

data:  newX[, i]
X-squared = 28.199, df = 10, p-value = 0.001677


$DIS

    Box-Ljung test

data:  newX[, i]
X-squared = 31.457, df = 10, p-value = 0.0004929


$GE

    Box-Ljung test

data:  newX[, i]
X-squared = 39.402, df = 10, p-value = 2.159e-05


$GS

    Box-Ljung test

data:  newX[, i]
X-squared = 35.037, df = 10, p-value = 0.0001231


$HD

    Box-Ljung test

data:  newX[, i]
X-squared = 12.51, df = 10, p-value = 0.2524


$IBM

    Box-Ljung test

data:  newX[, i]
X-squared = 27.601, df = 10, p-value = 0.002091


$INTC

    Box-Ljung test

data:  newX[, i]
X-squared = 56.863, df = 10, p-value = 1.414e-08


$JNJ

    Box-Ljung test

data:  newX[, i]
X-squared = 56.638, df = 10, p-value = 1.559e-08


$JPM

    Box-Ljung test

data:  newX[, i]
X-squared = 50.188, df = 10, p-value = 2.465e-07


$KO

    Box-Ljung test

data:  newX[, i]
X-squared = 37.219, df = 10, p-value = 5.186e-05


$MCD

    Box-Ljung test

data:  newX[, i]
X-squared = 40.098, df = 10, p-value = 1.628e-05


$MMM

    Box-Ljung test

data:  newX[, i]
X-squared = 28.234, df = 10, p-value = 0.001656


$MRK

    Box-Ljung test

data:  newX[, i]
X-squared = 19.687, df = 10, p-value = 0.03235


$MSFT

    Box-Ljung test

data:  newX[, i]
X-squared = 40.66, df = 10, p-value = 1.296e-05


$NKE

    Box-Ljung test

data:  newX[, i]
X-squared = 34.322, df = 10, p-value = 0.0001629


$PFE

    Box-Ljung test

data:  newX[, i]
X-squared = 49.19, df = 10, p-value = 3.759e-07


$PG

    Box-Ljung test

data:  newX[, i]
X-squared = 48.028, df = 10, p-value = 6.133e-07


$TRV

    Box-Ljung test

data:  newX[, i]
X-squared = 131.74, df = 10, p-value < 2.2e-16


$UNH

    Box-Ljung test

data:  newX[, i]
X-squared = 47.264, df = 10, p-value = 8.455e-07


$UTX

    Box-Ljung test

data:  newX[, i]
X-squared = 37.93, df = 10, p-value = 3.903e-05


$V

    Box-Ljung test

data:  newX[, i]
X-squared = 35.988, df = 10, p-value = 8.457e-05


$VZ

    Box-Ljung test

data:  newX[, i]
X-squared = 37.096, df = 10, p-value = 5.448e-05


$WMT

    Box-Ljung test

data:  newX[, i]
X-squared = 44.655, df = 10, p-value = 2.51e-06


$XOM

    Box-Ljung test

data:  newX[, i]
X-squared = 117.54, df = 10, p-value < 2.2e-16

# Apply the Ljung-Box test to absolute values of all returns in djall
apply(abs(djall), 2, Box.test, lag = 10, type = "Ljung")

$AAPL

    Box-Ljung test

data:  newX[, i]
X-squared = 910.73, df = 10, p-value < 2.2e-16


$AXP

    Box-Ljung test

data:  newX[, i]
X-squared = 3111.9, df = 10, p-value < 2.2e-16


$BA

    Box-Ljung test

data:  newX[, i]
X-squared = 1077, df = 10, p-value < 2.2e-16


$CAT

    Box-Ljung test

data:  newX[, i]
X-squared = 1567.6, df = 10, p-value < 2.2e-16


$CSCO

    Box-Ljung test

data:  newX[, i]
X-squared = 646.52, df = 10, p-value < 2.2e-16


$CVX

    Box-Ljung test

data:  newX[, i]
X-squared = 2203.8, df = 10, p-value < 2.2e-16


$DD

    Box-Ljung test

data:  newX[, i]
X-squared = 1432.9, df = 10, p-value < 2.2e-16


$DIS

    Box-Ljung test

data:  newX[, i]
X-squared = 1820.8, df = 10, p-value < 2.2e-16


$GE

    Box-Ljung test

data:  newX[, i]
X-squared = 2336.6, df = 10, p-value < 2.2e-16


$GS

    Box-Ljung test

data:  newX[, i]
X-squared = 2078.9, df = 10, p-value < 2.2e-16


$HD

    Box-Ljung test

data:  newX[, i]
X-squared = 1843.2, df = 10, p-value < 2.2e-16


$IBM

    Box-Ljung test

data:  newX[, i]
X-squared = 1025, df = 10, p-value < 2.2e-16


$INTC

    Box-Ljung test

data:  newX[, i]
X-squared = 1063.8, df = 10, p-value < 2.2e-16


$JNJ

    Box-Ljung test

data:  newX[, i]
X-squared = 1496.1, df = 10, p-value < 2.2e-16


$JPM

    Box-Ljung test

data:  newX[, i]
X-squared = 3038.8, df = 10, p-value < 2.2e-16


$KO

    Box-Ljung test

data:  newX[, i]
X-squared = 1526, df = 10, p-value < 2.2e-16


$MCD

    Box-Ljung test

data:  newX[, i]
X-squared = 1197.3, df = 10, p-value < 2.2e-16


$MMM

    Box-Ljung test

data:  newX[, i]
X-squared = 1272.8, df = 10, p-value < 2.2e-16


$MRK

    Box-Ljung test

data:  newX[, i]
X-squared = 1402.9, df = 10, p-value < 2.2e-16


$MSFT

    Box-Ljung test

data:  newX[, i]
X-squared = 824.11, df = 10, p-value < 2.2e-16


$NKE

    Box-Ljung test

data:  newX[, i]
X-squared = 1174.9, df = 10, p-value < 2.2e-16


$PFE

    Box-Ljung test

data:  newX[, i]
X-squared = 1187.9, df = 10, p-value < 2.2e-16


$PG

    Box-Ljung test

data:  newX[, i]
X-squared = 1115.4, df = 10, p-value < 2.2e-16


$TRV

    Box-Ljung test

data:  newX[, i]
X-squared = 2873.8, df = 10, p-value < 2.2e-16


$UNH

    Box-Ljung test

data:  newX[, i]
X-squared = 1407.5, df = 10, p-value < 2.2e-16


$UTX

    Box-Ljung test

data:  newX[, i]
X-squared = 1462.8, df = 10, p-value < 2.2e-16


$V

    Box-Ljung test

data:  newX[, i]
X-squared = 1165.2, df = 10, p-value < 2.2e-16


$VZ

    Box-Ljung test

data:  newX[, i]
X-squared = 1916.2, df = 10, p-value < 2.2e-16


$WMT

    Box-Ljung test

data:  newX[, i]
X-squared = 685.86, df = 10, p-value < 2.2e-16


$XOM

    Box-Ljung test

data:  newX[, i]
X-squared = 1918.2, df = 10, p-value < 2.2e-16

Note that the Ljung-Box test did not reject the null hypothesis when applied to the raw returns for AAPL, CAT, CSCO and HD.

3.3.2 Applying Ljung-Box tests to longer-interval returns

What happens when you apply the same analyses as in the previous exercise to the monthly returns rather than the daily returns? Does the amount of serial dependence in the data appear to increase or decrease?

# Create monthly log-returns from djx
djx_m <- apply.monthly(djx, sum)

# Apply Ljung-Box tests to raw and absolute values of djx_m
Box.test(djx_m, lag = 10, type = "Ljung")

    Box-Ljung test

data:  djx_m
X-squared = 21.099, df = 10, p-value = 0.02041

Box.test(abs(djx_m), lag = 10, type = "Ljung")

    Box-Ljung test

data:  abs(djx_m)
X-squared = 12.78, df = 10, p-value = 0.2362

# Create monthly log-returns from djall
djall_m <- apply.monthly(djall, colSums)

# Apply Ljung-Box tests to raw and absolute values of djall_m
apply(djall_m, 2, Box.test, lag = 10, type = "Ljung")

$AAPL

    Box-Ljung test

data:  newX[, i]
X-squared = 8.7875, df = 10, p-value = 0.5524


$AXP

    Box-Ljung test

data:  newX[, i]
X-squared = 24.336, df = 10, p-value = 0.006756


$BA

    Box-Ljung test

data:  newX[, i]
X-squared = 15.021, df = 10, p-value = 0.1313


$CAT

    Box-Ljung test

data:  newX[, i]
X-squared = 24.815, df = 10, p-value = 0.005707


$CSCO

    Box-Ljung test

data:  newX[, i]
X-squared = 11.565, df = 10, p-value = 0.3152


$CVX

    Box-Ljung test

data:  newX[, i]
X-squared = 8.576, df = 10, p-value = 0.5728


$DD

    Box-Ljung test

data:  newX[, i]
X-squared = 16.261, df = 10, p-value = 0.0924


$DIS

    Box-Ljung test

data:  newX[, i]
X-squared = 25.73, df = 10, p-value = 0.004122


$GE

    Box-Ljung test

data:  newX[, i]
X-squared = 13.919, df = 10, p-value = 0.1767


$GS

    Box-Ljung test

data:  newX[, i]
X-squared = 8.2406, df = 10, p-value = 0.6053


$HD

    Box-Ljung test

data:  newX[, i]
X-squared = 15.412, df = 10, p-value = 0.1178


$IBM

    Box-Ljung test

data:  newX[, i]
X-squared = 7.6343, df = 10, p-value = 0.6645


$INTC

    Box-Ljung test

data:  newX[, i]
X-squared = 7.7879, df = 10, p-value = 0.6495


$JNJ

    Box-Ljung test

data:  newX[, i]
X-squared = 12.386, df = 10, p-value = 0.2601


$JPM

    Box-Ljung test

data:  newX[, i]
X-squared = 8.1258, df = 10, p-value = 0.6165


$KO

    Box-Ljung test

data:  newX[, i]
X-squared = 23.391, df = 10, p-value = 0.00939


$MCD

    Box-Ljung test

data:  newX[, i]
X-squared = 15.718, df = 10, p-value = 0.108


$MMM

    Box-Ljung test

data:  newX[, i]
X-squared = 10.849, df = 10, p-value = 0.3694


$MRK

    Box-Ljung test

data:  newX[, i]
X-squared = 9.9017, df = 10, p-value = 0.4492


$MSFT

    Box-Ljung test

data:  newX[, i]
X-squared = 6.1226, df = 10, p-value = 0.8049


$NKE

    Box-Ljung test

data:  newX[, i]
X-squared = 6.1524, df = 10, p-value = 0.8023


$PFE

    Box-Ljung test

data:  newX[, i]
X-squared = 11.468, df = 10, p-value = 0.3222


$PG

    Box-Ljung test

data:  newX[, i]
X-squared = 6.427, df = 10, p-value = 0.7782


$TRV

    Box-Ljung test

data:  newX[, i]
X-squared = 11.718, df = 10, p-value = 0.3044


$UNH

    Box-Ljung test

data:  newX[, i]
X-squared = 19.614, df = 10, p-value = 0.03312


$UTX

    Box-Ljung test

data:  newX[, i]
X-squared = 20.689, df = 10, p-value = 0.02337


$V

    Box-Ljung test

data:  newX[, i]
X-squared = 10.395, df = 10, p-value = 0.4066


$VZ

    Box-Ljung test

data:  newX[, i]
X-squared = 15.69, df = 10, p-value = 0.1089


$WMT

    Box-Ljung test

data:  newX[, i]
X-squared = 6.0142, df = 10, p-value = 0.8141


$XOM

    Box-Ljung test

data:  newX[, i]
X-squared = 8.3262, df = 10, p-value = 0.597

apply(abs(djall_m), 2, Box.test, lag = 10, type = "Ljung")

$AAPL

    Box-Ljung test

data:  newX[, i]
X-squared = 25.433, df = 10, p-value = 0.004582


$AXP

    Box-Ljung test

data:  newX[, i]
X-squared = 63.524, df = 10, p-value = 7.759e-10


$BA

    Box-Ljung test

data:  newX[, i]
X-squared = 30.285, df = 10, p-value = 0.0007693


$CAT

    Box-Ljung test

data:  newX[, i]
X-squared = 40.407, df = 10, p-value = 1.436e-05


$CSCO

    Box-Ljung test

data:  newX[, i]
X-squared = 14.867, df = 10, p-value = 0.137


$CVX

    Box-Ljung test

data:  newX[, i]
X-squared = 12.093, df = 10, p-value = 0.2789


$DD

    Box-Ljung test

data:  newX[, i]
X-squared = 21.997, df = 10, p-value = 0.01512


$DIS

    Box-Ljung test

data:  newX[, i]
X-squared = 19.021, df = 10, p-value = 0.03999


$GE

    Box-Ljung test

data:  newX[, i]
X-squared = 93.633, df = 10, p-value = 9.992e-16


$GS

    Box-Ljung test

data:  newX[, i]
X-squared = 27.807, df = 10, p-value = 0.001938


$HD

    Box-Ljung test

data:  newX[, i]
X-squared = 10.328, df = 10, p-value = 0.4122


$IBM

    Box-Ljung test

data:  newX[, i]
X-squared = 7.8389, df = 10, p-value = 0.6446


$INTC

    Box-Ljung test

data:  newX[, i]
X-squared = 4.981, df = 10, p-value = 0.8924


$JNJ

    Box-Ljung test

data:  newX[, i]
X-squared = 12.031, df = 10, p-value = 0.283


$JPM

    Box-Ljung test

data:  newX[, i]
X-squared = 50.321, df = 10, p-value = 2.329e-07


$KO

    Box-Ljung test

data:  newX[, i]
X-squared = 12.959, df = 10, p-value = 0.2259


$MCD

    Box-Ljung test

data:  newX[, i]
X-squared = 8.647, df = 10, p-value = 0.5659


$MMM

    Box-Ljung test

data:  newX[, i]
X-squared = 9.2386, df = 10, p-value = 0.5096


$MRK

    Box-Ljung test

data:  newX[, i]
X-squared = 22.105, df = 10, p-value = 0.01458


$MSFT

    Box-Ljung test

data:  newX[, i]
X-squared = 4.8293, df = 10, p-value = 0.9023


$NKE

    Box-Ljung test

data:  newX[, i]
X-squared = 6.8513, df = 10, p-value = 0.7394


$PFE

    Box-Ljung test

data:  newX[, i]
X-squared = 17.793, df = 10, p-value = 0.05855


$PG

    Box-Ljung test

data:  newX[, i]
X-squared = 19.134, df = 10, p-value = 0.03859


$TRV

    Box-Ljung test

data:  newX[, i]
X-squared = 6.0277, df = 10, p-value = 0.8129


$UNH

    Box-Ljung test

data:  newX[, i]
X-squared = 58.968, df = 10, p-value = 5.679e-09


$UTX

    Box-Ljung test

data:  newX[, i]
X-squared = 16.905, df = 10, p-value = 0.0765


$V

    Box-Ljung test

data:  newX[, i]
X-squared = 14.577, df = 10, p-value = 0.1483


$VZ

    Box-Ljung test

data:  newX[, i]
X-squared = 10.138, df = 10, p-value = 0.4285


$WMT

    Box-Ljung test

data:  newX[, i]
X-squared = 4.7833, df = 10, p-value = 0.9052


$XOM

    Box-Ljung test

data:  newX[, i]
X-squared = 10.339, df = 10, p-value = 0.4113

The amount of serial dependence appears to decrease when we move from daily to monthly returns.

3.4 Looking at the extremes in volatile return series

3.4.1 Extreme values in volatile time series

When you take a long series of iid data, such as several thousand observations, and select a small subset of the most extreme observations, like less than 100, then these extremes appear at random and the spaces or gaps between the extremes follow a distribution that is very close to exponential. When we carry out the same exercise for a volatile financial log-return series then the extremes appear in clusters during periods of high volatility. This is another feature of real log-return data that we need to take account of when building models.

You will investigate the irregular time series djx_extremes which contains the 100 most extreme negative log-returns of the Dow Jones index between 1985 and 2015. You will compare it with iid_extremes which contains the 100 most extreme values in an iid series of the same length as djx_extremes. To do this, you will use the object exp_quantiles, which contains 100 theoretical quantiles of the standard exponential distribution. These can be used to construct a Q-Q plot of each dataset against the exponential reference distribution.

set.seed(123)

#load data and get log returns
djreturns <- diff(log(DJ))[-1]

#get mean and sigma to generate normal distribution
npars <- c(mu = mean(djreturns), sigma = sd(djreturns))
n <- length(djreturns)

#calculate extreme losses
dj_losses <- -djreturns
djx_extremes <- dj_losses[dj_losses>0.0274]

#calculate extreme losses on iid dist
iid <- rnorm(n)*npars[2] + npars[1]
iid <- as.xts(iid, order.by = index(djreturns))
iid_losses <- -iid
iid_extremes <- iid_losses[iid_losses>0.02459]

#get 100 theroical quantiles of the standar exponential distribution
exp_quantiles <- qexp(ppoints(100))

# Partition plotting area into 3 pieces
par(mfrow = c(1, 3))

# Plot djx_extremes
plot(djx_extremes, type = "h")

# Compute the spaces between the times of the extremes
djx_spaces <- diff(time(djx_extremes))

# Make a histogram of these spaces
hist(as.numeric(djx_spaces))

# Make a Q-Q plot of djx_spaces against exp_quantiles
qqplot(exp_quantiles, djx_spaces)

# Carry out the previous 4 steps for iid_extremes
plot(iid_extremes, type = "h")
iid_spaces <- diff(time(iid_extremes))
hist(as.numeric(iid_spaces))

qqplot(exp_quantiles, iid_spaces)

Note how strongly clustered the extreme returns are compared with the extremes in the iid series.

3.4.2 Cross correlations between risk-factor return series

Many risk-factor returns are correlated with each other in the same time period. However, in the same way that there tends to be only weak serial correlation within series, there tends to be only weak cross correlation between series in different time periods.

The picture changes dramatically when we look at the absolute values, which are often strongly correlated both within and between series.

In this exercise, you will investigate cross correlations between the daily log-returns of the Dow Jones, FTSE100 and SMI equity indexes. When the function acf() is applied to a multivariate time series, we obtain a matrix of plots with the usual sample acf plots on the diagonal, and plots of the correlations between different series at different lags off the diagonal.

One thing to note here is that the US and European series are slightly out of step. The European markets tend to follow the US, so we see some evidence of cross correlation between US returns on one day and European returns on the next.

data("SMI")
data("FTSE")
smix <- diff(log(SMI))[-1]["2005/2015"]
ftsex <- diff(log(FTSE))[-1]["2005/2015"]
djx <- diff(log(DJ))[-1]["2005/2015"]
indexes <- merge(djx, ftsex, smix, all = TRUE)
indexes <- na.omit(indexes) 

# Make a time series plot of indexes with plot.zoo and a pairwise scatterplot with pairs
plot.zoo(indexes)

pairs(as.zoo(indexes))

# Calculate the sample correlation matrix of indexes
cor(indexes)

           X.DJI    X.FTSE    X.SSMI
X.DJI  1.0000000 0.5903534 0.5408565
X.FTSE 0.5903534 1.0000000 0.8226086
X.SSMI 0.5408565 0.8226086 1.0000000

# Plot the sample acfs and cross-correlation functions for the returns in indexes
acf(indexes)

# Plot the sample acfs and cross-correlations functions for the absolute values of indexes
acf(abs(indexes))

Do you see the evidence that the US market leads the European markets in the third of the four plots?

3.5 The stylized facts of return series

3.5.1 Volatility and correlation of FX returns

You will discover evidence of volatility and serial dependence in daily and weekly exchange-rate log-returns. The dataset fx contains daily log-returns for the “EUR_USD”, “GBP_USD”, “JPY_USD” and “CHF_USD” exchange rates while fx_w contains the corresponding weekly log-returns. Both are in your workspace.

Note that foreign exchange trading takes place every day of the week although the lower volume of trading at weekends leads to an unusual weekly correlation cycle which will be evident in one of the pictures you create.

data("CHF_USD")
fx1 <- diff(log(EUR_USD))[-1]["2011/2016"]
fx2 <- diff(log(GBP_USD))[-1]["2011/2016"]
fx3 <- diff(log(JPY_USD))[-1]["2011/2016"]
fx4 <- diff(log(CHF_USD))[-1]["2011/2016"]
fx <- merge(fx1, fx2, fx3, fx4, all = TRUE)
fx_w<- apply.weekly(fx, colSums)

# Plot fx and fx_w
plot.zoo(fx, type = "h")

plot.zoo(fx_w, type = "h")

# Make acf plots of fx and the absolute values of fx
acf(fx)

acf(abs(fx))

# Apply the Ljung-Box test to the components of fx and their absolute values
apply(fx, 2, Box.test, lag = 10, type = "Ljung")

$EUR.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 195.33, df = 10, p-value < 2.2e-16


$GBP.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 206.44, df = 10, p-value < 2.2e-16


$JPY.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 141.56, df = 10, p-value < 2.2e-16


$CHF.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 310.46, df = 10, p-value < 2.2e-16

apply(abs(fx), 2, Box.test, lag = 10, type = "Ljung")

$EUR.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 307.05, df = 10, p-value < 2.2e-16


$GBP.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 235.89, df = 10, p-value < 2.2e-16


$JPY.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 264.54, df = 10, p-value < 2.2e-16


$CHF.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 376.51, df = 10, p-value < 2.2e-16

# Make acf plots of fx_w and the absolute values of fx_w
acf(fx_w)

acf(abs(fx_w))

# Apply the Ljung-Box test to the components of fx_w and their absolute values
apply(fx_w, 2, Box.test, lag = 10, type = "Ljung")

$EUR.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 10.364, df = 10, p-value = 0.4091


$GBP.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 8.8896, df = 10, p-value = 0.5426


$JPY.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 7.7924, df = 10, p-value = 0.6491


$CHF.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 10.961, df = 10, p-value = 0.3605

apply(abs(fx_w), 2, Box.test, lag = 10, type = "Ljung")

$EUR.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 29.168, df = 10, p-value = 0.00117


$GBP.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 47.097, df = 10, p-value = 9.066e-07


$JPY.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 46.153, df = 10, p-value = 1.345e-06


$CHF.USD

    Box-Ljung test

data:  newX[, i]
X-squared = 34.449, df = 10, p-value = 0.000155

There is strong evidence of serial dependence in the absolute values of the weekly log-returns.

3.5.2 Volatility and correlation of interest-rate data

You will explore whether volatility and serial dependence is also a feature of daily and monthly interest-rate log-returns. The dataset zcb_x contains daily log-returns for the 1-year, 5-year and 10-year Canadian zero-coupon bond yields while zcbx_m contains the corresponding monthly log-returns.

# Make acf plots of zcb_x and the absolute values of zcb_x
acf(zcbx[-1])

acf(abs(zcbx[-1]))

# Apply the Ljung-Box test to the components of zcb_x and their absolute values
apply(zcbx, 2, Box.test, lag = 10, type = "Ljung")

$`1.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 34.333, df = 10, p-value = 0.0001622


$`5.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 17.186, df = 10, p-value = 0.07036


$`10.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 26.048, df = 10, p-value = 0.003676

apply(abs(zcbx), 2, Box.test, lag = 10, type = "Ljung")

$`1.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 2988.9, df = 10, p-value < 2.2e-16


$`5.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 5839.9, df = 10, p-value < 2.2e-16


$`10.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 3281.4, df = 10, p-value < 2.2e-16

# Make acf plots of zcbx_m and the absolute values of zcbx_m
acf(zcbx_m)

acf(abs(zcbx_m))

# Apply the Ljung-Box test to the components of zcbx_m and their absolute values
apply(zcbx_m, 2, Box.test, lag = 10, type = "Ljung")

$`1.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 8.3928, df = 10, p-value = 0.5905


$`5.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 5.4943, df = 10, p-value = 0.8558


$`10.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 4.246, df = 10, p-value = 0.9356

apply(abs(zcbx_m), 2, Box.test, lag = 10, type = "Ljung")

$`1.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 41.993, df = 10, p-value = 7.52e-06


$`5.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 8.993, df = 10, p-value = 0.5328


$`10.00y`

    Box-Ljung test

data:  newX[, i]
X-squared = 11.912, df = 10, p-value = 0.291

The monthly log-returns for the 5 and 10 year yields don’t appear to show much serial dependence.

4 Estimating portfolio value-at-risk

4.1 Value-at-risk and expected shortfall

4.1.1 Computing VaR and ES for normal distribution

The standard function qnorm() calculates quantiles of a normal distribution from the probability p, the mean, and standard deviation, and thus can be used to calculate value-at-risk (VaR). The function ESnorm() from the QRM package calculates the expected shortfall (ES) for a normal distribution from the probability p, location parameter mu, and scale parameter sd: qnorm(p, mean = 0, sd = 1) ESnorm(p, mu = 0, sd = 1)

Common numeric values for p include 0.95 and 0.99 for confidence levels of 95% and 99%, respectively. you will compute and display VaR and ES for a normal distribution N(μ,σ2) with mean μ and standard deviation σ. In the process, you will use the new functions for sequence generation and adding straight lines to a plot. You can read about their arguments by typing in ?seq and ?abline to your console.

The variables mu and sigma contain the estimated mean and standard deviation of the Dow Jones index returns for 2008-2009 contained in djx.

djx <- diff(log(dj))[-1,]["2008/2009"]
sigma <- sd(djx)
mu <- mean(djx)

# Make a sequence of 100 x-values going from -4*sigma to 4*sigma
xvals <- seq(from = -4*sigma, to = 4*sigma, length.out = 100)

# Compute the density of a N(mu, sigma^2) distribution at xvals
ndens <- dnorm(xvals, mean = mu, sd = sigma)

# Plot ndens against xvals
plot(xvals, ndens, type = "l")

# Compute the 99% VaR and 99% ES of a N(mu, sigma^2) distribution
VaR99 <- qnorm(0.99, mean = mu, sd = sigma)
VaR99

[1] 0.04612495

ES99 <- ESnorm(0.99, mu = mu, sd = sigma)
ES99

[1] 0.05290841

# Draw vertical lines at VaR99 and ES99 in red and green
abline(v = VaR99, col = "red")

abline(v = ES99, col = "green")

– ES99 is only 14.7% bigger than VaR99. For heavy-tailed distributions, the difference can be much greater.

4.2 International equity portfolio

4.2.1 Examining risk factors for international equity portfolio

The UK investor in UK, US, and Swiss equities is exposed to 5 risk factors; the data is contained in riskfactors, a multivariate dataset.

In this exercise, you will recall some of the tests and techniques that you learned earlier for showing that these risk factors are heavier tailed than normal, highly volatile and subject to profound serial dependencies.

data("USD_GBP")
data("CHF_GBP")
X.FTSE <- FTSE["2000/2012"]
X.GSPC <- SP500["2000/2012"]
X.SSMI <- SMI["2000/2012"]
USD.GBP <- USD_GBP["2000/2012"]
CHF.GBP <- CHF_GBP["2000/2012"]
riskfactors <- merge(X.FTSE,  X.GSPC, X.SSMI, USD.GBP, CHF.GBP , all = TRUE)
riskfactors <- na.omit(riskfactors)

# Plot the risk-factor data
plot.zoo(riskfactors)

# Calculate the log-returns, assign to returns, and plot
returns <- diff(log(riskfactors))[-1, ]
plot.zoo(returns)

# Use apply() to carry out the Jarque-Bera test for all 5 series
apply(returns, 2, jarque.test)

$X.FTSE

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 4209.4, p-value < 2.2e-16
alternative hypothesis: greater


$X.GSPC

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 6961.6, p-value < 2.2e-16
alternative hypothesis: greater


$X.SSMI

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 5158.1, p-value < 2.2e-16
alternative hypothesis: greater


$USD.GBP

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 2586.4, p-value < 2.2e-16
alternative hypothesis: greater


$CHF.GBP

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 7759.2, p-value < 2.2e-16
alternative hypothesis: greater

# Make a Q-Q plot against normal for the 5th return series and add a reference line
qqnorm(returns[, 5])
qqline(returns[, 5])

# Make a picture of the sample acfs for returns and their absolute values
acf(returns)

acf(abs(returns))

All five risk factors are clearly non-normal and show strong serial and cross dependencies.

4.2.2 Historical simulation

Suppose that a UK investor has invested 30% of her wealth in the FTSE index, 40% in the S&P 500 index, and 30% in the SMI index.

For different vectors of log-returns for the 5 risk factors, the function lossop() computes the loss or gain incurred by the investor when her total wealth is 1. The function can also be applied to a 5-dimensional time series of log-returns to obtain a time series of historically-simulated losses and gains corresponding to each vector of log-returns in the time series.

The function lossop() is the so-called loss operator for the portfolio and has been specially written for this exercise. In general, for each new portfolio, a specific function has to be written to compute portfolio losses and gains.

In this exercise, you will form historically simulated losses and examine them. This is a necessary prelude to using these data to estimate VaR and ES.

lossop <- function (xseries, wts = c(0.3, 0.4, 0.3)){
    if (is.xts(xseries)) 
        x <- coredata(xseries)
    else if (is.matrix(xseries)) 
        x <- xseries
    else x <- matrix(xseries, nrow = 1)
    ll <- apply(x, 1, function(x, wts) {
        1 - (wts[1] * exp(x[1]) + wts[2] * exp(x[2] + x[4]) + 
            wts[3] * exp(x[3] + x[5]))
    }, wts = wts)
    if (is.xts(xseries)) 
        ll <- xts(ll, time(xseries))
    ll
}

# Calculate the loss from a log-return of -0.1 for all risk factors
lossop(rep(-0.1, 5))

[1] 0.1554372

# Apply lossop() to returns and plot hslosses
hslosses <- lossop(returns)
plot(hslosses)

# Form a Q-Q plot of hslosses against normal
qqnorm(hslosses)

# Plot the sample acf of hslosses and their absolute values
acf(hslosses)

acf(abs(hslosses))

Note how the features of the underlying risk-factor returns (heavy tails and serial dependence) are present in the historically simulated losses.

4.2.3 Estimating VaR and ES

Now you are ready to estimate VaR and ES for the international equity investor using the historically simulated losses and gains in hslosses.

You will do this by two methods. First, you will apply a simple non-parametric method using a sample quantile to estimate VaR and the average of values exceeding the sample quantile to estimate ES.

Then, you will compare these estimates with the values obtained when you assume that the hslosses have a normal distribution. Obviously, this is a very bad assumption and you should compare the two sets of estimates to see which are more conservative.

# Estimate the 99th sample percentile of the distribution of hslosses
quantile(hslosses, 0.99)

       99% 
0.03076173 

# Estimate the 99% ES
mean(hslosses[hslosses >= quantile(hslosses, 0.99)])

[1] 0.04184655

# Estimate the mean and standard deviation of hslosses
mu <- mean(hslosses)
sigma <- sd(hslosses)

# Compute the 99% quantile of a normal distribution
qnorm(0.99, mean = mu, sd = sigma)

[1] 0.02602973

# Compute the 99% ES of a normal distribution
ESnorm(0.99, mu = mu, sd = sigma)

[1] 0.02983979

This is what we expected. The estimates derived from a normal assumption are much less conservative than the estimates derived using the non-parametric method. That is not to say that the latter method is the best possible, but it is a reasonable first attempt to estimate the risk measures.

4.3 Option portfolio and Black Scholes

4.3.1 Compute Black-Scholes price of an option

The Black_Scholes() function in the package qrmtools can be used to price European call and put options using the standard Black-Scholes options pricing formula for a non-dividend-paying stock.

In this exercise you will price in succession: an out-of-the-money European call, an in-the-money European call, an in-the-money European put and an out-of-the-money European put. An option is in-the-money if immediate exercise would result in a positive payout and out-of-the-money if it would not.

The main point of the exercise is to understand the different risk factors that go into the price calculation: the current stock price, the current volatility and the current interest rate.

# Set the interest rate r to be 0.01, the volatility sigma to be 0.2 and the strike K to be 100
r <- 0.01
sigma <- 0.2
K <- 100


library(qrmtools)

package 㤼㸱qrmtools㤼㸲 was built under R version 4.0.3Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Attaching package: 㤼㸱qrmtools㤼㸲

The following objects are masked from 㤼㸱package:QRM㤼㸲:

    dGEV, dGPD, pGEV, pGPD, qGEV, qGPD, rGEV, rGPD

The following object is masked from 㤼㸱package:timeSeries㤼㸲:

    returns

# Look at the arguments of the Black_Scholes function
args(Black_Scholes)

function (t, S, r, sigma, K, T, type = c("call", "put")) 
NULL

# Price a European call option that matures in one year if the current stock price is 80
Black_Scholes(0, 80, r, sigma, K, 1, "call")

[1] 1.302245

# Price a European call option that matures in one year if the current stock price is 120
Black_Scholes(0, 120, r, sigma, K, 1, "call")

[1] 22.94188

# Price a European put option that matures in one year if the current stock price is 80
Black_Scholes(0, 80, r, sigma, K, 1,"put")

[1] 20.30723

# Price a European put option that matures in one year if the current stock price is 120
Black_Scholes(0, 120, r, sigma, K, 1,"put")

[1] 1.94686

Did you see how dramatically the option values changed as the stock price moved from out-of-the-money to in-the-money or vice versa?

4.3.2 Equity and implied volatility risk factors

To analyze the risk of a portfolio consisting of an option, it is necessary to consider changes in all three risk factors: stock price, volatility and interest rates. Here, you will focus on the first two of these risk factors and assume that interest rates do not change much over short time intervals. The daily risk-factor values for the period 1990-2010 are contained in riskfactors and the corresponding log-returns in returns; both multivariate datasets are loaded in your workspace.

Volatility is a new risk factor that hasn’t been considered so far in this course. It is represented by the VIX index which is constructed from the implied volatilities of a wide range of options on the S&P 500 index:

data("VIX")
X.GSPC <- SP500["1990/2010"]
X.VIX <- VIX["1990/2010"]
riskfactors <- merge(X.GSPC, X.VIX, all = TRUE)
riskfactors <- na.omit(riskfactors)
returns <- diff(log(riskfactors))[-1, ]

names(returns)

[1] "X.GSPC" "X.VIX" 

You will be able to verify whether the log-returns of volatility behave like other return data you have encountered, and to see how they vary with the log-returns of the S&P 500 index.

# Plot the risk factors and the log-returns
plot.zoo(riskfactors)

plot.zoo(returns)

# Make a scatterplot of the two return series
plot(as.matrix(returns))

# Apply the Jarque-Bera test to the returns and make a Q-Q plot of the volatility log-returns
apply(returns, 2, jarque.test)

$X.GSPC

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 17397, p-value < 2.2e-16
alternative hypothesis: greater


$X.VIX

    Jarque-Bera Normality Test

data:  newX[, i]
JB = 4412.6, p-value < 2.2e-16
alternative hypothesis: greater

qqnorm(returns[, 2])

# Create the sample acf of the returns and absolute returns
acf(returns)

acf(abs(returns))

# Calculate the correlation between the log-returns
cor(returns)

           X.GSPC      X.VIX
X.GSPC  1.0000000 -0.6978308
X.VIX  -0.6978308  1.0000000

It is clear that the log-returns of the VIX index show the same stylized facts as other returns that you have analyzed - non-normality, heavy tails, volatility, serial dependence in the absolute values but not the raw values. Moreover, they are negatively correlated with the log-returns of the SP500 index

4.4 Historical simulation for the option example

4.4.1 Historical simulation of losses for option portfolio

Suppose that an investor has invested one unit of wealth in a single European call option on the S&P 500 index. The function lossop() computes the loss or gain incurred by the investor over a one-day time horizon due to changes in the log stock price or changes in the log volatility. As before, this function has been written specially for the particular portfolio in this exercise: lossop(xseries, S, sigma) The first argument contains the log returns corresponding to the stock price and volatility risk factors, either in a series or in form c(stock_risk, volatility_risk), S is the current stock price, and sigma is the current volatility.

Changes in the interest rate over the time horizon will be neglected as being of lesser importance.

In this exercise, you will form the historically simulated losses for the option portfolio and examine their properties before estimating VaR and ES in the next exercise. The interest rate, strike price, and maturity have been set to r = 0.01, K = 100 and T = 1, respectively.

lossop <- function (xseries, r = 0.01, K = 100, T = 1, sigma = 0.2, S = 100){
    if (is.xts(xseries)) 
        x <- coredata(xseries)
    else if (is.matrix(xseries)) 
        x <- xseries
    else x <- matrix(xseries, nrow = 1)
    ll <- apply(x, 1, function(x, r, K, T, sigma, S) {
        deltat <- 1/250
        V_t0 <- Black_Scholes(0, S, r, sigma, K, T, "call")
        V_t1 = Black_Scholes(deltat, exp(log(S) + x[1]), r, exp(log(sigma) + 
            x[2]), K, T, "call")
        -(V_t1 - V_t0)/V_t0
    }, r = r, K = K, T = T, sigma = sigma, S = S)
    if (is.xts(xseries)) 
        ll <- xts(ll, time(xseries))
    ll
}

# Calculate the first loss
lossop(c(-0.1, -0.1), S = 80, sigma = 0.2)

[1] 0.8030928

# Calculate the second loss
lossop(c(-0.1, 0.1), S = 100, sigma = 0.2)

[1] 0.4380754

# Create and plot hslosses
hslosses <- lossop(returns, S = 100, sigma = 0.2)
plot(hslosses)

# Form a Q-Q plot of hslosses against normal
qqnorm(hslosses)

# Plot the sample acf of raw data and absolute values in hslosses
acf(hslosses)

acf(abs(hslosses))

4.4.2 Estimating VaR and ES for option portfolio

Now you are ready to estimate VaR and ES for the investor in the European call option using the historically simulated losses and gains in hslosses.

Once again, you will do this by two methods. First, you will apply a non-parametric method using a sample quantile to estimate VaR and calculate the average of values exceeding the same quantile to estimate ES.

Then, you will compare these estimates with the values obtained when you assume that the hslosses have a normal distribution. Like in the previous exercise, this is a bad assumption and you should compare the two sets of estimates to see which are more conservative.

# Estimate the 99.5% percentile of the distribution
quantile(hslosses, 0.995)

    99.5% 
0.1618376 

# Estimate the 99.5% ES
mean(hslosses[hslosses >= quantile(hslosses, 0.995)])

[1] 0.2211743

# Estimate the mean and standard deviation of hslosses
mu <- mean(hslosses)
sigma <- sd(hslosses)

# Compute the 99.5% quantile of a normal distribution
qnorm(0.995, mean = mu, sd = sigma)

[1] 0.1441452

# Compute the 99.5% ES of a normal distribution
ESnorm(0.995, mu = mu, sd = sigma)

[1] 0.1622111

It will be no surprise to you now that the normal distribution greatly underestimates these measures of risk.

4.4.3 Computing VaR for weekly losses

In this final exercise, you will test your understanding by computing an empirical estimate of VaR for weekly losses in the returns data. You will have to repeat the analysis of the previous exercise, but this time, you need to:

Find the weekly log-returns of returns using apply.weekly(). Use these weekly log-returns to simulate the losses of the two risk factors through lossop(). Note that the lossop() function has been adjusted in your workspace so that it correctly calculates the losses and gains of the option portfolio for a one-week time horizon. It still takes in arguments as follows:

lossop(xseries, S, sigma)

lossop <- function (xseries, r = 0.01, K = 100, T = 1, sigma = 0.2, S = 100){
    if (is.xts(xseries)) 
        x <- coredata(xseries)
    else if (is.matrix(xseries)) 
        x <- xseries
    else x <- matrix(xseries, nrow = 1)
    ll <- apply(x, 1, function(x, r, K, T, sigma, S) {
        deltat <- 5/250
        V_t0 <- Black_Scholes(0, S, r, sigma, K, T, "call")
        V_t1 = Black_Scholes(deltat, exp(log(S) + x[1]), r, exp(log(sigma) + 
            x[2]), K, T, "call")
        -(V_t1 - V_t0)/V_t0
    }, r = r, K = K, T = T, sigma = sigma, S = S)
    if (is.xts(xseries)) 
        ll <- xts(ll, time(xseries))
    ll
}

Your challenge is to compute the 99% VaR for weekly changes in value of the European call option in returns when the current stock price is S = 120 and the current volatility is sigma = 0.25.

return_w <- apply.weekly(returns, colSums)
hslosses <- lossop(return_w, S = 120, sigma = 0.25)
quantile(hslosses, 0.995)

    99.5% 
0.2402091