This notebook is an implementation of Jón DanÃelsson's Financial Risk Forecasting (Wiley, 2011) in R 4.0, with annotations and introductory examples. The introductory examples (Appendix) are similar to Appendix B in the original book.
Bullet point numbers correspond to the R/MATLAB Listing numbers in the original book, referred to henceforth as FRF.
More details can be found at the book website: https://www.financialriskforecasting.com/
Last updated: July 2020
Copyright 2011- 2020 Jón DanÃelsson. This code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. The GNU General Public License is available at: https://www.gnu.org/licenses/.
Chapter 0. Appendix - Introduction
Chapter 1. Financial Markets, Prices and Risk
Chapter 2. Univariate Volatility Modeling
Chapter 3. Multivariate Volatility Models
Chapter 4. Risk Measures
Chapter 5. Implementing Risk Forecasts
Chapter 6. Analytical Value-at-Risk for Options and Bonds
Chapter 7. Simulation Methods for VaR for Options and Bonds
Chapter 8. Backtesting and Stress Testing
Chapter 9. Extreme Value Theory
Created in R 4.0 (July 2020)
# Entering and Printing Data in R
# Listing R.1
# Last updated June 2018
#
#
x = 10 # assign x the value 10
print(x) # print x
# Vectors, Matrices and Sequences in R
# Listing R.2
# Last updated June 2018
#
#
y = c(1,3,5,7,9) # create vector using c()
print(y)
print(y[3]) # calling 3rd element (R indices start at 1)
print(dim(y)) # gives NULL since y is a vector, not a matrix
print(length(y)) # as expected, y has length 5
v = matrix(nrow=2,ncol=3) # fill a 2 x 3 matrix with NaN values (default)
print(dim(v)) # as expected, v is size (2,3)
w = matrix(c(1,2,3),nrow=6,ncol=3) # repeats matrix twice by rows, thrice by columns
print(w)
s = 1:10 # s is a list of integers from 1 to 10 inclusive
print(s)
# Basic Summary Statistics in R
# Listing R.3
# Last updated June 2018
#
#
y=matrix(c(3.1,4.15,9))
sum(y) # sum of all elements of y
prod(y) # product of all elements of y
max(y) # maximum value of y
min(y) # minimum value of y
range(y) # min, max value of y
mean(y) # arithmetic mean
median(y) # median
var(y) # variance
cov(y) # covar matrix = variance for single vector
cor(y) # corr matrix = [1] for single vector
sort(y) # sorting in ascending order
log(y) # natural log
# Calculating Moments in R
# Listing R.4
# Last updated June 2018
#
#
library(moments)
mean(y) # mean
var(y) # variance
sd(y) # unbiased standard deviation, by default
skewness(y) # skewness
kurtosis(y) # kurtosis
# Basic Matrix Operations in R
# Listing R.5
# Last updated June 2018
#
#
z = matrix(c(1,2,3,4),2,2) # z is a 2 x 2 matrix
x = matrix(c(1,2),1,2) # x is a 1 x 2 matrix
## Note: z * x is undefined since the two matrices are not conformable
z %*% t(x) # this evaluates to a 2 x 1 matrix
rbind(z,x) # "stacking" z and x vertically
cbind(z,t(x)) # "stacking z and x' horizontally
## Note: dimensions must match along the combining axis
# Statistical Distributions in R
# Listing R.6
# Last updated June 2018
#
#
q = seq(from = -3, to = 3, length = 7) # specify a set of values
p = seq(from = 0.1, to = 0.9, length = 9) # specify a set of probabilities
qnorm(p, mean = 0, sd = 1) # element-wise inverse Normal quantile
pt(q, df = 4) # element-wise cdf under Student-t(4)
dchisq(q, df = 2) # element-wise pdf under Chisq(2)
## Similar syntax for other distributions
## q for quantile, p for cdf, d for pdf
## followed by the abbreviation of the distribution
## One can also obtain pseudorandom samples from distributions
x = rt(100, df = 5) # Sampling 100 times from TDist with 5 df
y = rnorm(50, mean = 0, sd = 1) # Sampling 50 times from a standard normal
## Given data, we obtain MLE estimates of distribution parameters with package MASS:
library(MASS)
res = fitdistr(x, densfun = "normal") # Fitting x to normal dist
print(res)
# Statistical Tests in R
# Listing R.7
# Last updated June 2018
#
#
library(tseries)
x = rt(500, df = 5) # Create hypothetical dataset x
jarque.bera.test(x) # Jarque-Bera test for normality
Box.test(x, lag = 20, type = c("Ljung-Box")) # Ljung-Box test for serial correlation
# Time Series in R
# Listing R.8
# Last updated June 2018
#
#
x = rt(60, df = 5) # Create hypothetical dataset x
par(mfrow=c(1,2), pty='s')
acf(x,20) # autocorrelation for lags 1:20
pacf(x,20) # partial autocorrelation for lags 1:20
# Loops and Functions in R
# Listing R.9
# Last updated June 2018
#
#
## For loops
for (i in 3:7) # iterates through [3,4,5,6,7]
print(i^2)
## If-else loops
X = 10
if (X %% 3 == 0) {
print("X is a multiple of 3")
} else {
print("X is not a multiple of 3")
}
## Functions (example: a simple excess kurtosis function)
excess_kurtosis = function(x, excess = 3){ # note: excess optional, default=3
m4 = mean((x-mean(x))^4)
excess_kurt = m4/(sd(x)^4) - excess
excess_kurt
}
x = rt(60, df = 5) # Create hypothetical dataset x
excess_kurtosis(x)
# Basic Graphs in R
# Listing R.10
# Last updated June 2018
#
#
y = rnorm(50, mean = 0, sd = 1)
par(mfrow=c(2,2)) # sets up space for subplots
barplot(y) # bar plot
plot(y,type='l') # line plot
hist(y) # histogram
plot(y) # scatter plot
# Miscellaneous Useful Functions in R
# Listing R.11
# Last updated June 2018
#
#
## Convert objects from one type to another with as.integer() etc
## To check type, use typeof(object)
x = 8.0
print(typeof(x))
x = as.integer(x)
print(typeof(x))
# Download S&P500 data in R
# Listing 1.1
# Last updated August 2019
#
#
price = read.csv('index.csv')
y=diff(log(price$Index)) # calculate returns
plot(y) # plot returns
# Sample statistics in R
# Listing 1.3
# Last updated July 2020
#
#
library(moments)
library(tseries)
mean(y)
sd(y)
min(y)
max(y)
skewness(y)
kurtosis(y)
jarque.bera.test(y)
# ACF plots and the Ljung-Box test in R
# Listing 1.5
# Last updated July 2020
#
#
library(MASS)
library(stats)
par(mfrow=c(1,2), pty="s")
q = acf(y,20)
q1 = acf(y^2,20)
Box.test(y, lag = 20, type = c("Ljung-Box"))
Box.test(y^2, lag = 20, type = c("Ljung-Box"))
# QQ plots in R
# Listing 1.7
# Last updated June 2018
#
#
library(car)
par(mfrow=c(1,2), pty="s")
qqPlot(y)
qqPlot(y,distribution="t",df=5)
# Download stock prices in R
# Listing 1.9
# Last updated June 2018
#
#
p = read.csv('stocks.csv')
y=apply(log(p),2,diff)
print(cor(y)) # correlation matrix
# ARCH and GARCH estimation in R
# Listing 2.1
# Last updated July 2020
#
#
library(rugarch)
p = read.csv('index.csv')
## We multiply returns by 100 and de-mean them
y=diff(log(p$Index))*100
y=y-mean(y)
## GARCH(1,1)
spec1 = ugarchspec(variance.model = list( garchOrder = c(1, 1)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE))
res1 = ugarchfit(spec = spec1, data = y)
## ARCH(1)
spec2 = ugarchspec(variance.model = list( garchOrder = c(1, 0)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE))
res2 = ugarchfit(spec = spec2, data = y)
##Â tGARCH(1,1)
spec3 = ugarchspec(variance.model = list( garchOrder = c(1, 1)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE),
distribution.model = "std")
res3 = ugarchfit(spec = spec3, data = y)
## plot(res) shows various graphical analysis, works in command line
# Advanced ARCH and GARCH estimation in R
# Listing 2.3
# Last updated July 2020
#
#
## Normal APARCH(1,1)
spec4 = ugarchspec(variance.model = list(model="apARCH", garchOrder = c(1, 1)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE))
res4 = ugarchfit(spec = spec4, data = y)
## show(res4)
## Normal APARCH(1,1) with fixed delta
spec5 = ugarchspec(variance.model = list(model="apARCH", garchOrder = c(1, 1)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE), fixed.pars=list(delta=2))
res5 = ugarchfit(spec = spec5, data = y)
show(res5)
# Download stock prices in R
# Listing 3.1
# Last updated August 2019
#
#
p = read.csv('stocks.csv')
y=apply(log(p),2,diff) # calculate returns
y = y[,1:2] # consider first two stocks
y[,1] = y[,1]-mean(y[,1]) # subtract mean
y[,2] = y[,2]-mean(y[,2])
TT = dim(y)[1]
# EWMA in R
# Listing 3.3
# Last updated August 2019
#
#
## create a matrix to hold covariance matrix for each t
EWMA = matrix(nrow=TT,ncol=3)
lambda = 0.94
S = cov(y) # initial (t=1) covar matrix
EWMA[1,] = c(S)[c(1,4,2)] # extract var and covar
for (i in 2:dim(y)[1]){
S = lambda*S+(1-lambda)* y[i-1,] %*% t(y[i-1,])
EWMA[i,] = c(S)[c(1,4,2)]
}
EWMArho = EWMA[,3]/sqrt(EWMA[,1]*EWMA[,2]) # calculate correlations
print(head(EWMArho))
print(tail(EWMArho))
# GOGARCH in R
# Listing 3.5
# Last updated August 2019
#
#
library(rmgarch)
spec = gogarchspec(mean.model = list(armaOrder = c(0, 0),
include.mean =FALSE),
variance.model = list(model = "sGARCH",
garchOrder = c(1,1)) ,
distribution.model = "mvnorm"
)
fit = gogarchfit(spec = spec, data = y)
show(fit)
# DCC in R
# Listing 3.7
# Last updated August 2019
#
#
xspec = ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE))
uspec = multispec(replicate(2, xspec))
spec = dccspec(uspec = uspec, dccOrder = c(1, 1), distribution = 'mvnorm')
res = dccfit(spec, data = y)
H=res@mfit$H
DCCrho=vector(length=dim(y)[1])
for(i in 1:dim(y)[1]){
DCCrho[i] = H[1,2,i]/sqrt(H[1,1,i]*H[2,2,i])
}
# Sample statistics in R
# Listing 3.9
# Last updated August 2019
#
#
matplot(cbind(EWMArho,DCCrho),type='l',las=1,lty=1,col=2:3,ylab="")
mtext("Correlations",side=2,line=0.3,at=1,las=1,cex=0.8)
legend("bottomright",c("EWMA","DCC"),lty=1,col=2:3,bty="n",cex=0.7)
# ES in R
# Listing 4.1
# Last updated August 2016
#
#
p = c(0.5,0.1,0.05,0.025,0.01,0.001)
VaR = -qnorm(p)
ES = dnorm(qnorm(p))/p
cat("Probabilities:", paste0(p*100,"%"), "\n",
"VaR:", VaR, "\n",
"ES:", ES)
# Download stock prices in R
# Listing 5.1
# Last updated July 2020
#
#
p = read.csv('stocks.csv')
## Calculate returns. Note that first column is dates
y = apply(log(p),2,diff)
## Specify portfolio value and VaR probability
portfolio_value = 1000
p = 0.01
# Univariate HS in R
# Listing 5.3
# Last updated July 2020
#
#
y1 = y[,1] # select one asset
ys = sort(y1) # sort returns
op = ceiling(length(y1)*p) # p percent smallest, rounded up
VaR1 = -ys[op]*portfolio_value
print(VaR1)
# Multivariate HS in R
# Listing 5.5
# Last updated July 2020
#
#
w = matrix(c(0.3,0.2,0.5)) # vector of portfolio weights
## The number of columns of the left matrix must be the same as the number of rows of the right matrix
yp = y %*% w # obtain portfolio returns
yps = sort(yp)
VaR2 = -yps[op]*portfolio_value
print(VaR2)
# Univariate ES in R
# Listing 5.7
# Last updated August 2019
#
#
ES1 = -mean(ys[1:op])*portfolio_value
print(ES1)
# Normal VaR in R
# Listing 5.9
# Last updated August 2019
#
#
sigma = sd(y1) # estimate volatility
VaR3 = -sigma * qnorm(p) * portfolio_value
print(VaR3)
# Portfolio normal VaR in R
# Listing 5.11
# Last updated August 2019
#
#
sigma = sqrt(t(w) %*% cov(y) %*% w)[1] # portfolio volatility
## Note: the trailing [1] is to convert a single element matrix to float
VaR4 = -sigma * qnorm(p)*portfolio_value
print(VaR4)
# Student-t VaR in R
# Listing 5.13
# Last updated July 2020
#
#
library(QRM)
scy1 = (y1)*100 # scale the returns
res = fit.st(scy1)
sigma1 = unname(res$par.ests['sigma']/100) # rescale the volatility
nu = unname(res$par.ests['nu'])
## Note: We are removing the names of the fit.st output using unname()
VaR5 = - sigma1 * qt(df=nu,p=p) * portfolio_value
print(VaR5)
# Normal ES in R
# Listing 5.15
# Last updated June August 2019
#
#
sigma = sd(y1)
ES2 = sigma*dnorm(qnorm(p))/p * portfolio_value
print(ES2)
# Direct integration ES in R
# Listing 5.17
# Last updated July 2020
#
#
VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$value/p*portfolio_value
print(ES)
# MA normal VaR in R
# Listing 5.19
# Last updated June August 2019
#
#
WE=20
for (t in seq(length(y1)-5,length(y1))){
t1=t-WE+1
window= y1[t1:t] # estimation window
sigma=sd(window)
VaR6 = -sigma * qnorm(p) * portfolio_value
print(VaR6)
}
# EWMA VaR in R
# Listing 5.21
# Last updated July 2020
#
#
lambda = 0.94
s11 = var(y1) # initial variance, using unconditional
for (t in 2:length(y1)){
s11 = lambda * s11 + (1-lambda) * y1[t-1]^2
}
VaR7 = -qnorm(p) * sqrt(s11) * portfolio_value
print(VaR7)
# Three-asset EWMA VaR in R
# Listing 5.23
# Last updated August 2019
#
#
s = cov(y) # initial covariance
for (t in 2:dim(y)[1]){
s = lambda*s + (1-lambda)*y[t-1,] %*% t(y[t-1,])
}
sigma = sqrt(t(w) %*% s %*% w)[1] # portfolio vol
## Note: trailing[1] is to convert single element matrix to float
VaR8 = -sigma * qnorm(p) * portfolio_value
print(VaR8)
# Univariate GARCH in R
# Listing 5.25
# Last updated July 2020
#
#
library(rugarch)
spec = ugarchspec(variance.model = list( garchOrder = c(1, 1)),
mean.model = list( armaOrder = c(0,0),include.mean = FALSE))
res = ugarchfit(spec = spec, data = y1)
omega = res@fit$coef['omega']
alpha = res@fit$coef['alpha1']
beta = res@fit$coef['beta1']
sigma2 = omega + alpha * tail(y1,1)^2 + beta * tail(res@fit$var,1)
VaR9 = -sqrt(sigma2) * qnorm(p) * portfolio_value
names(VaR9)="VaR"
print(VaR9)
# Black-Scholes function in R
# Listing 6.1
# Last updated 2011
#
#
bs = function(X, P, r, sigma, T){
d1 = (log(P/X) + (r + 0.5*sigma^2)*(T))/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
Call = P*pnorm(d1,mean=0,sd=1)-X*exp(-r*(T))*pnorm(d2,mean=0,sd=1)
Put = X*exp(-r*(T))*pnorm(-d2,mean=0,sd=1)-P*pnorm(-d1,mean=0,sd=1)
Delta.Call = pnorm(d1, mean = 0, sd = 1)
Delta.Put = Delta.Call - 1
Gamma = dnorm(d1, mean = 0, sd = 1)/(P*sigma*sqrt(T))
return(list(Call=Call,Put=Put,Delta.Call=Delta.Call,Delta.Put=Delta.Put,Gamma=Gamma))
}
# Black-Scholes in R
# Listing 6.3
# Last updated July 2020
#
#
f = bs(X = 90, P = 100, r = 0.05, sigma = 0.2, T = 0.5)
print(f)
# Transformation in R
# Listing 7.1
# Last updated 2011
#
#
x = seq(-3, 3, by = 0.1)
plot(x, pnorm(x), type = "l")
# Various RNs in R
# Listing 7.3
# Last updated 2011
#
#
set.seed(12) # set seed
S = 10
runif(S)
rnorm(S)
rt(S,4)
# Price bond in R
# Listing 7.5
# Last updated July 2020
#
#
yield = c(5.00, 5.69, 6.09, 6.38, 6.61,
6.79, 6.94, 7.07, 7.19, 7.30) # yield curve
T = length(yield) # number of time periods
r = 0.07 # initial yield rate
Par = 10 # par value
coupon = r * Par # coupon payments
cc = rep(coupon, T) # vector of cash flows
cc[T] = cc[T] + Par # add par to cash flows
P = sum(cc/((1+yield/100)^(1:T))) # calculate price
print(P)
# Simulate yields in R
# Listing 7.7
# Last updated August 2019
#
#
set.seed(12) # set seed
sigma = 1.5 # daily yield volatiltiy
S = 8 # number of simulations
r = rnorm(S, 0, sigma) # generate random numbers
ysim = matrix(nrow=length(yield),ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
matplot(ysim,type='l')
# Simulate bond prices in R
# Listing 7.9
# Last updated November 2020
#
#
SP = rep(NA, length = S)
for (i in 1:S){ # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(1:T)))
}
SP = SP-(mean(SP) - P) # correct for mean
par(mfrow=c(1,2), pty="s")
barplot(SP)
hist(SP,probability=TRUE)
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))
S = 50000
r = rnorm(S, 0, sigma) # generate random numbers
ysim = matrix(nrow=length(yield),ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
SP = rep(NA, length = S)
for (i in 1:S){ # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(1:T)))
}
SP = SP-(mean(SP) - P) # correct for mean
par(mfrow=c(1,2), pty="s")
barplot(SP)
hist(SP,probability=TRUE)
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))
