Introduction¶
This notebook is an implementation of Jón DanÃelsson's Financial Risk Forecasting (Wiley, 2011) in R 3.6, 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: August 2019
Copyright 2011- 2019 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/.
Appendix: An Introduction to R¶
Created in R 3.6 (August 2019)
- R.1: Entering and Printing Data
- R.2: Vectors, Matrices and Sequences
- R.3: Importing Data (to be updated)
- R.4: Basic Summary Statistics
- R.5: Calculating Moments
- R.6: Basic Matrix Operations
- R.7: Statistical Distributions
- R.8: Statistical Tests
- R.9: Time Series
- R.10: Loops and Functions
- R.11: Basic Graphs
- R.12: Miscellaneous Useful Functions
In [1]:
# Entering and Printing Data in R
# Listing R.1
# Last updated June 2018
#
#
x = 10 # assign x the value 10
print(x) # print x
In [2]:
# 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)
In [3]:
# Importing Data in R
# Listing R.3
# Last updated June 2018
#
#
## There are many data sources for financial data, for instance
## Yahoo Finance, AlphaVantage and Quandl. However, some of the
## free data sources have numerous issues with accuracy and
## handling of missing data, so only CSV importing is shown here.
##
## For csv data, one can use read.csv to read it
##
## Example:
## data = read.csv('Ch1aprices.csv', header=TRUE, sep=',')
## one can use the zoo() function from the package zoo
## to turn the data into a timeseries (see Listing 1.1/1.2)
In [4]:
# Basic Summary Statistics in R
# Listing R.4
# 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
In [5]:
# Calculating Moments in R
# Listing R.5
# 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
In [6]:
# Basic Matrix Operations in R
# Listing R.6
# 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
In [7]:
# Statistical Distributions in R
# Listing R.7
# 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)
In [8]:
# Statistical Tests in R
# Listing R.8
# 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
In [9]:
# Time Series in R
# Listing R.9
# 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
In [10]:
# Loops and Functions in R
# Listing R.10
# 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)
In [11]:
# Basic Graphs in R
# Listing R.11
# 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
In [12]:
# Miscellaneous Useful Functions in R
# Listing R.12
# 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))
Chapter 1: Financial Markets, Prices and Risk¶
- 1.1: Loading hypothetical stock prices, converting to returns, plotting returns
- 1.3: Summary statistics for returns timeseries
- 1.5: Autocorrelation function (ACF) plots, Ljung-Box test
- 1.7: Quantile-Quantile (QQ) plots
- 1.9: Correlation matrix between different stocks
In [1]:
# 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
In [4]:
# Sample statistics in R
# Listing 1.3
# Last updated August 2019
#
#
library(moments)
library(tseries)
mean(y)
sd(y)
min(y)
max(y)
skewness(y)
kurtosis(y)
acf(y,1)
acf(y^2,1)
jarque.bera.test(y)
Box.test(y, lag = 20, type = c("Ljung-Box"))
Box.test(y^2, lag = 20, type = c("Ljung-Box"))
In [2]:
# ACF plots and the Ljung-Box test in R
# Listing 1.5
# Last updated August 2019
#
#
library(MASS)
library(stats)
par(mfrow=c(1,2), pty="s")
q = acf(y,20)
q1 = acf(y^2,20)
In [16]:
# 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)
In [4]:
# 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
Chapter 2: Univariate Volatility Modelling¶
- 2.1: GARCH and t-GARCH estimation
- 2.3: APARCH estimation
In [9]:
# ARCH and GARCH estimation in R
# Listing 2.1
# Last updated August 2019
#
#
library(rugarch)
p = read.csv('index.csv')
y=diff(log(p$Index))*100
y=y-mean(y)
## We multiply returns by 100 and de-mean them
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)
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)
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
In [15]:
# Advanced ARCH and GARCH estimation in R
# Listing 2.3
# Last updated August 2019
#
#
library(rugarch)
## 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)
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)
Chapter 3: Multivariate Volatility Models¶
- 3.1: Loading hypothetical stock prices
- 3.3: EWMA estimation
- 3.5: OGARCH estimation
- 3.7: DCC estimation
- 3.9: Comparison of EWMA, OGARCH, DCC
In [17]:
# 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]
In [19]:
# 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))
In [32]:
# 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)
In [26]:
# DCC in R
# Listing 3.7
# Last updated August 2019
#
#
library(rmgarch)
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])
}
In [33]:
# 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)
Chapter 4: Risk Measures¶
- 4.1: Expected Shortfall (ES) estimation under normality assumption
In [25]:
# 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
print(ES)
Chapter 5: Implementing Risk Forecasts¶
- 5.1: Loading hypothetical stock prices, converting to returns
- 5.3: Univariate HS Value at Risk (VaR)
- 5.5: Multivariate HS VaR
- 5.7: Univariate ES VaR
- 5.9: Normal VaR
- 5.11: Portfolio Normal VaR
- 5.13: Student-t VaR
- 5.15: Normal ES VaR
- 5.17: Direct Integration Normal ES VaR
- 5.19: MA Normal VaR
- 5.21: EWMA VaR
- 5.23: Two-asset EWMA VaR
- 5.25: GARCH(1,1) VaR
In [87]:
# Download stock prices in R
# Listing 5.1
# Last updated August 2019
#
#
p = read.csv('stocks.csv')
y=apply(log(p),2,diff) # calculate returns. note first column is dates
portfolio_value = 1000
p = 0.01 # probability
# we want the number p*dim(y)[1] to be an integer, so
# remove the first observations to make it so
print(dim(y))
t1=dim(y)[1]-100*floor(dim(y)[1]*p)+1
y=y[t1:dim(y)[1],]
print(dim(y))
y1=y[,1]
In [88]:
# Univariate HS in R
# Listing 5.3
# Last updated August 2016
#
#
ys = sort(y1) # sort returns
op = length(y1)*p # p percent smallest
VaR1 = -ys[op]*portfolio_value
print(VaR1)
In [89]:
# Multivariate HS in R
# Listing 5.5
# Last updated August 2019
#
#
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)
In [90]:
# Univariate ES in R
# Listing 5.7
# Last updated August 2019
#
#
ES1 = -mean(ys[1:op])*portfolio_value
print(ES1)
In [91]:
# Normal VaR in R
# Listing 5.9
# Last updated August 2019
#
#
sigma = sd(y1) # estimate volatility
VaR3 = -sigma * qnorm(p) * portfolio_value
print(VaR3)
In [92]:
# 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)
In [93]:
# Student-t VaR in R
# Listing 5.13
# Last updated August August 2019
#
#
library(QRM)
scy1=(y1)*100 # scale the returns
res=fit.st(scy1)
sigma1=res$par.ests[3]/100 # rescale the volatility
nu=res$par.ests[1]
VaR5 = - sigma1 * qt(df=nu,p=p) * portfolio_value
print(VaR5)
In [74]:
# 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)
In [94]:
# Direct integration ES in R
# Listing 5.17
# Last updated August 2019
#
#
VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$portfolio_value/p*portfolio_value
print(ES)
In [81]:
# 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)
}
In [95]:
# EWMA VaR in R
# Listing 5.21
# Last updated August 2019
#
#
lambda = 0.94;
s11 = var(y1[1:30]); # initial variance
for (t in 2:length(y1)){
s11 = lambda * s11 + (1-lambda) * y1[t-1]^2
}
VaR7 = -qnorm(p) * sqrt(s11) * portfolio_value
print(VaR7)
In [96]:
# 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: [1] is to convert single element matrix to float
VaR8 = -sigma * qnorm(p) * portfolio_value
print(VaR8)
In [126]:
# Univariate GARCH in R
# Listing 5.25
# Last updated August 2019
#
#
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[1]
alpha = res@fit$coef[2]
beta = res@fit$coef[3]
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)
Chapter 6: Analytical Value-at-Risk for Options and Bonds¶
- 6.1: Black-Scholes function definition
- 6.3: Black-Scholes option price calculation example
In [104]:
# 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))
}
In [40]:
# Black-Scholes in R
# Listing 6.3
# Last updated August 2016
#
#
f=bs(90,100,0.05,0.2,0.5)
print(f)
Chapter 7: Simulation Methods for VaR for Options and Bonds¶
- 7.1: Plotting normal distribution transformation
- 7.3: Random number generation from Uniform(0,1), Normal(0,1)
- 7.5: Bond pricing using yield curve
- 7.7: Yield curve simulations
- 7.9: Bond price simulations
- 7.11: Black-Scholes analytical pricing of call
- 7.13: Black-Scholes Monte Carlo simulation pricing of call
- 7.15: Option density plots
- 7.17: VaR simulation of portfolio with only underlying
- 7.19: VaR simulation of portfolio with only call
- 7.21: VaR simulation of portfolio with call, put and underlying
- 7.23: Simulated two-asset returns
- 7.25: Two-asset portfolio VaR
- 7.27: Two-asset portfolio VaR with a call
In [41]:
# Transformation in R
# Listing 7.1
# Last updated 2011
#
#
x=seq(-3,3,by=0.1)
plot(x,pnorm(x),type="l")
In [42]:
# Various RNs in R
# Listing 7.3
# Last updated 2011
#
#
set.seed(12) # set seed
S=10
runif(S)
rnorm(S)
rt(S,4)
In [98]:
# Price bond in R
# Listing 7.5
# Last updated August 2019
#
#
yield=c(5.00, 5.69, 6.09, 6.38, 6.61,
6.79, 6.94, 7.07, 7.19, 7.30) # yield curve
r=0.07 # initial yield rate
Par=10 # par value
coupon=r*Par # coupon payments
cc=1:10*0+coupon # vector of cash flows
cc[10]=cc[10]+Par # add par to cash flows
P=sum(cc/((1+yield/100)^(1:length(yield)))) # calculate price
print(P)
In [106]:
# 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')