# Introduction¶

This notebook is an implementation of Jón Daníelsson's Financial Risk Forecasting (Wiley, 2011) in R 3.4.3, 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: June 2018

Copyright 2011, 2016, 2018 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.4.3 (June 2018)

• 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

[1] 10

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)

[1] 1 3 5 7 9
[1] 5
NULL
[1] 5
[1] 2 3
[,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    1    1    1
[5,]    2    2    2
[6,]    3    3    3
[1]  1  2  3  4  5  6  7  8  9 10

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.
##
##
## Example:
## 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

16.25
115.785
9
3.1
1. 3.1
2. 9
5.41666666666667
4.15
 9.90583
 9.90583
 1
1. 3.1
2. 4.15
3. 9
 1.1314 1.42311 2.19722
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

5.41666666666667
 9.90583
3.14735338551827
0.61960293299089
1.5
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

 7 10
 1 3 2 4 1 2
 1 3 1 2 4 2
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)

1. -1.2815515655446
2. -0.841621233572914
3. -0.524400512708041
4. -0.2533471031358
5. 0
6. 0.2533471031358
7. 0.524400512708041
8. 0.841621233572914
9. 1.2815515655446
1. 0.0199709840358594
2. 0.0580582617584078
3. 0.186950483150029
4. 0.5
5. 0.813049516849971
6. 0.941941738241592
7. 0.980029015964141
1. 0
2. 0
3. 0
4. 0.5
5. 0.303265329856317
6. 0.183939720585721
7. 0.111565080074215
      mean          sd
0.02829141   1.42600592
(0.14260059) (0.10083385)

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

	Jarque Bera Test

data:  x
X-squared = 48.511, df = 2, p-value = 2.924e-11

	Box-Ljung test

data:  x
X-squared = 21.823, df = 20, p-value = 0.3502

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)

[1] 9
[1] 16
[1] 25
[1] 36
[1] 49
[1] "X is not a multiple of 3"

2.9748560343394
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))

[1] "double"
[1] "integer"


# Chapter 1: Financial Markets, Prices and Risk¶

• 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 [13]:
# Download S&P500 data in R
# Listing 1.1
# Last updated June 2018
#
#

library(tseries)
library(zoo)

y=diff(log(price))  # calculate returns
plot(y)             # plot returns
y=coredata(y)       # get core data

In [14]:
# Sample statistics in R
# Listing 1.3
# Last updated June 2018
#
#

library(moments)

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"))

0.000258159900825998
0.0100057286437631
-0.101955486273023
0.106735895029117
Index: 0.152633269896337
Index: 16.981171461092
	Jarque Bera Test

data:  y
X-squared = 46251, df = 2, p-value < 2.2e-16

	Box-Ljung test

data:  y
X-squared = 93.488, df = 20, p-value = 1.809e-11

	Box-Ljung test

data:  y^2
X-squared = 2543, df = 20, p-value < 2.2e-16

In [15]:
# ACF plots and the Ljung-Box test in R
# Listing 1.5
# Last updated June 2018
#
#

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 [17]:
# Download stock prices in R
# Listing 1.9
# Last updated June 2018
#
#

y=diff(log(p))
print(cor(y)) # correlation matrix

          A         B         C
A 1.0000000 0.2296842 0.2126192
B 0.2296842 1.0000000 0.1450511
C 0.2126192 0.1450511 1.0000000


# Chapter 2: Univariate Volatility Modelling¶

• 2.1: GARCH and t-GARCH estimation
• 2.3: APARCH estimation
In [18]:
# ARCH and GARCH estimation in R
# Listing 2.1
# Last updated June 2018
#
#

library(tseries)

y=diff(log(p))*100
y=y-mean(y)

## We multiply returns by 100 and de-mean them

library(fGarch)

garchFit(~ garch(1,0), data = y,include.mean=FALSE)
garchFit(~ garch(4,0), data = y,include.mean=FALSE)
garchFit(~ garch(4,1), data = y,include.mean=FALSE)
garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="std",trace=F)
res=garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="sstd",trace=F)
## plot(res)

## plot(res) shows various graphical analysis, works in command line

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(1, 0)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               1 0
Max GARCH Order:           1
Maximum Order:             1
Conditional Dist:          norm
h.start:                   2
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              1.000573

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.977953e-17 1.977953e-17    0.0    FALSE
omega   1.000000e-06 1.000000e+02    0.1     TRUE
alpha1  1.000000e-08 1.000000e+00    0.1     TRUE
gamma1 -1.000000e+00 1.000000e+00    0.1    FALSE
delta   0.000000e+00 2.000000e+00    2.0    FALSE
skew    1.000000e-01 1.000000e+01    1.0    FALSE
shape   1.000000e+00 1.000000e+01    4.0    FALSE
Index List of Parameters to be Optimized:
omega alpha1
2      3
Persistence:                  0.1

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     16797.442: 0.100000 0.100000
1:     7857.9063:  1.06968 0.344367
2:     7851.0951:  1.07156 0.309265
3:     7691.1536: 0.801715 0.230117
4:     7671.9202: 0.710551 0.267167
5:     7668.4564: 0.695436 0.364405
6:     7667.5808: 0.673414 0.349908
7:     7667.3737: 0.684009 0.340375
8:     7667.3737: 0.683887 0.340279
9:     7667.3737: 0.683883 0.340296

Final Estimate of the Negative LLH:
LLH:  7670.624    norm LLH:  1.351414
omega    alpha1
0.6846670 0.3402963

R-optimhess Difference Approximated Hessian Matrix:
omega    alpha1
omega  -4440.117 -1189.811
alpha1 -1189.811 -1746.219
attr(,"time")
Time difference of 0.2310231 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.349035 secs

Title:
GARCH Modelling

Call:
garchFit(formula = ~garch(1, 0), data = y, include.mean = FALSE)

Mean and Variance Equation:
data ~ garch(1, 0)
<environment: 0x00000000256dcba8>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega   alpha1
0.68467  0.34030

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega    0.68467     0.01660    41.25   <2e-16 ***
alpha1   0.34030     0.02647    12.86   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-7670.624    normalized:  -1.351414

Description:
Thu Jul 05 09:39:31 2018 by user: fanj5

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(4, 0)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               4 0
Max GARCH Order:           4
Maximum Order:             4
Conditional Dist:          norm
h.start:                   5
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              1.000573

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.977953e-17 1.977953e-17  0.000    FALSE
omega   1.000000e-06 1.000000e+02  0.100     TRUE
alpha1  1.000000e-08 1.000000e+00  0.025     TRUE
alpha2  1.000000e-08 1.000000e+00  0.025     TRUE
alpha3  1.000000e-08 1.000000e+00  0.025     TRUE
alpha4  1.000000e-08 1.000000e+00  0.025     TRUE
gamma1 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma2 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma3 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma4 -1.000000e+00 1.000000e+00  0.100    FALSE
delta   0.000000e+00 2.000000e+00  2.000    FALSE
skew    1.000000e-01 1.000000e+01  1.000    FALSE
shape   1.000000e+00 1.000000e+01  4.000    FALSE
Index List of Parameters to be Optimized:
omega alpha1 alpha2 alpha3 alpha4
2      3      4      5      6
Persistence:                  0.1

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     13453.426: 0.100000 0.0250000 0.0250000 0.0250000 0.0250000
1:     7947.5194: 0.766879 0.416384 0.393311 0.394380 0.385552
2:     7836.0369: 0.394424 0.659988 0.529821 0.556091 0.470103
3:     7570.1274: 0.188761 0.569629 0.434980 0.466180 0.374129
4:     7400.2949: 0.254371 0.431193 0.292256 0.342890 0.241043
5:     7324.0323: 0.286334 0.338362 0.205786 0.267016 0.168650
6:     7293.0743: 0.408281 0.253747 0.153533 0.213959 0.145472
7:     7282.1500: 0.359125 0.219652 0.132971 0.186340 0.125381
8:     7278.2992: 0.389663 0.201469 0.171289 0.146613 0.154680
9:     7276.3805: 0.383261 0.185352 0.157128 0.168508 0.138276
10:     7276.2321: 0.388222 0.184365 0.157080 0.170252 0.141319
11:     7276.1961: 0.386767 0.179634 0.153546 0.169314 0.141475
12:     7276.1886: 0.390209 0.179084 0.155771 0.166470 0.145040
13:     7276.1567: 0.389449 0.179508 0.153142 0.171417 0.142641
14:     7276.1546: 0.389085 0.179029 0.153070 0.171071 0.142566
15:     7276.1538: 0.389487 0.178699 0.153390 0.170928 0.142877
16:     7276.1534: 0.389286 0.178496 0.153393 0.170794 0.142757
17:     7276.1533: 0.389383 0.178459 0.153452 0.170757 0.142829
18:     7276.1532: 0.389359 0.178400 0.153455 0.170692 0.142833
19:     7276.1532: 0.389404 0.178377 0.153488 0.170661 0.142865
20:     7276.1532: 0.389390 0.178355 0.153488 0.170633 0.142861
21:     7276.1532: 0.389412 0.178344 0.153506 0.170614 0.142873
22:     7276.1532: 0.389406 0.178337 0.153506 0.170604 0.142871
23:     7276.1532: 0.389416 0.178327 0.153519 0.170584 0.142876
24:     7276.1532: 0.389414 0.178327 0.153519 0.170581 0.142872
25:     7276.1532: 0.389415 0.178325 0.153520 0.170578 0.142876
26:     7276.1532: 0.389415 0.178325 0.153520 0.170578 0.142875
27:     7276.1532: 0.389415 0.178325 0.153520 0.170578 0.142874
28:     7276.1532: 0.389415 0.178325 0.153520 0.170578 0.142874
29:     7276.1532: 0.389415 0.178325 0.153520 0.170578 0.142874

Final Estimate of the Negative LLH:
LLH:  7279.404    norm LLH:  1.282488
omega    alpha1    alpha2    alpha3    alpha4
0.3898611 0.1783252 0.1535197 0.1705782 0.1428745

R-optimhess Difference Approximated Hessian Matrix:
omega     alpha1     alpha2     alpha3     alpha4
omega  -7504.084 -1937.3228 -2076.9591 -2177.6312 -2195.9507
alpha1 -1937.323 -3294.6402 -1007.5122  -836.7555  -920.6524
alpha2 -2076.959 -1007.5122 -3512.2848 -1050.3154  -872.9675
alpha3 -2177.631  -836.7555 -1050.3154 -2965.3649  -887.0863
alpha4 -2195.951  -920.6524  -872.9675  -887.0863 -4012.0553
attr(,"time")
Time difference of 0.1500149 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.885088 secs

Title:
GARCH Modelling

Call:
garchFit(formula = ~garch(4, 0), data = y, include.mean = FALSE)

Mean and Variance Equation:
data ~ garch(4, 0)
<environment: 0x0000000026460d38>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega   alpha1   alpha2   alpha3   alpha4
0.38986  0.17833  0.15352  0.17058  0.14287

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega    0.38986     0.01467   26.569  < 2e-16 ***
alpha1   0.17833     0.01933    9.223  < 2e-16 ***
alpha2   0.15352     0.01897    8.094 6.66e-16 ***
alpha3   0.17058     0.02113    8.073 6.66e-16 ***
alpha4   0.14287     0.01742    8.201 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-7279.404    normalized:  -1.282488

Description:
Thu Jul 05 09:39:32 2018 by user: fanj5

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(4, 1)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               4 1
Max GARCH Order:           4
Maximum Order:             4
Conditional Dist:          norm
h.start:                   5
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              1.000573

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.977953e-17 1.977953e-17  0.000    FALSE
omega   1.000000e-06 1.000000e+02  0.100     TRUE
alpha1  1.000000e-08 1.000000e+00  0.025     TRUE
alpha2  1.000000e-08 1.000000e+00  0.025     TRUE
alpha3  1.000000e-08 1.000000e+00  0.025     TRUE
alpha4  1.000000e-08 1.000000e+00  0.025     TRUE
gamma1 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma2 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma3 -1.000000e+00 1.000000e+00  0.100    FALSE
gamma4 -1.000000e+00 1.000000e+00  0.100    FALSE
beta1   1.000000e-08 1.000000e+00  0.800     TRUE
delta   0.000000e+00 2.000000e+00  2.000    FALSE
skew    1.000000e-01 1.000000e+01  1.000    FALSE
shape   1.000000e+00 1.000000e+01  4.000    FALSE
Index List of Parameters to be Optimized:
omega alpha1 alpha2 alpha3 alpha4  beta1
2      3      4      5      6     11
Persistence:                  0.9

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     7165.7198: 0.100000 0.0250000 0.0250000 0.0250000 0.0250000 0.800000
1:     7109.1967: 0.0825475 0.0255746 0.0243028 0.0235361 0.0242867 0.792114
2:     7071.4966: 0.0730530 0.0354922 0.0321951 0.0303386 0.0325046 0.794490
3:     7063.6076: 0.0547706 0.0349453 0.0304894 0.0277453 0.0309942 0.789613
4:     7032.8817: 0.0516722 0.0423519 0.0362660 0.0317462 0.0382105 0.803908
5:     7018.0850: 0.0377531 0.0426118 0.0348197 0.0275753 0.0381511 0.816436
6:     7006.5201: 0.0337169 0.0442224 0.0342497 0.0225642 0.0380474 0.834489
7:     6995.0451: 0.0270589 0.0439473 0.0318539 0.0139028 0.0346248 0.849765
8:     6972.3749: 0.0164921 0.0422046 0.0268601 1.00000e-08 0.0177477 0.901740
9:     6969.2620: 0.0152882 0.0416307 0.0262682 1.00000e-08 0.0170805 0.901041
10:     6969.0439: 0.0154206 0.0413251 0.0255587 1.00000e-08 0.0155547 0.901387
11:     6968.3478: 0.0157871 0.0431875 0.0255366 1.00000e-08 0.0131081 0.903018
12:     6967.8004: 0.0151576 0.0445899 0.0256731 1.00000e-08 0.0102252 0.904266
13:     6966.9421: 0.0145706 0.0456234 0.0267802 1.00000e-08 0.00485967 0.908808
14:     6966.9319: 0.0141104 0.0435070 0.0262002 1.00000e-08 0.00473697 0.909987
15:     6966.7066: 0.0144670 0.0428438 0.0265538 1.00000e-08 0.00465038 0.910940
16:     6966.6159: 0.0141635 0.0434128 0.0273319 1.00000e-08 0.00390984 0.910744
17:     6966.5821: 0.0140831 0.0437446 0.0274338 1.00000e-08 0.00317827 0.911716
18:     6966.4603: 0.0139729 0.0444497 0.0267042 1.00000e-08 0.00243714 0.911847
19:     6966.4082: 0.0141390 0.0450441 0.0265411 1.00000e-08 0.00142341 0.912263
20:     6966.3775: 0.0138485 0.0441489 0.0271645 1.00000e-08 0.00120085 0.912796
21:     6966.3621: 0.0138036 0.0433420 0.0278375 1.00000e-08 0.00107025 0.913492
22:     6966.3113: 0.0136578 0.0444053 0.0273375 1.00000e-08 0.000621025 0.913427
23:     6966.2952: 0.0137060 0.0440261 0.0283940 1.00000e-08 3.40461e-05 0.913461
24:     6966.2707: 0.0134085 0.0446216 0.0265024 1.00000e-08 1.00000e-08 0.914921
25:     6966.2695: 0.0133741 0.0446273 0.0265026 1.00000e-08 1.00000e-08 0.914907
26:     6966.2686: 0.0133759 0.0446575 0.0265240 1.00000e-08 1.00000e-08 0.914914
27:     6966.2678: 0.0133375 0.0447097 0.0264937 1.00000e-08 1.00000e-08 0.914891
28:     6966.2661: 0.0133601 0.0448333 0.0264112 1.00000e-08 1.00000e-08 0.914888
29:     6966.2643: 0.0133450 0.0450560 0.0262138 1.00000e-08 1.00000e-08 0.914847
30:     6966.2628: 0.0133996 0.0455006 0.0258141 1.00000e-08 1.00000e-08 0.914811
31:     6966.2613: 0.0133567 0.0459534 0.0254299 1.00000e-08 1.00000e-08 0.914725
32:     6966.2603: 0.0133379 0.0458781 0.0256410 1.00000e-08 1.00000e-08 0.914734
33:     6966.2597: 0.0133971 0.0458157 0.0257720 1.00000e-08 1.00000e-08 0.914572
34:     6966.2596: 0.0133898 0.0458145 0.0257706 1.00000e-08 1.00000e-08 0.914569
35:     6966.2596: 0.0133873 0.0458203 0.0257758 1.00000e-08 1.00000e-08 0.914570
36:     6966.2595: 0.0133845 0.0458156 0.0257887 1.00000e-08 1.00000e-08 0.914561
37:     6966.2595: 0.0133889 0.0458367 0.0257954 1.00000e-08 1.00000e-08 0.914536
38:     6966.2595: 0.0133880 0.0458629 0.0257741 1.00000e-08 1.00000e-08 0.914535
39:     6966.2595: 0.0133877 0.0458512 0.0257966 1.00000e-08 1.00000e-08 0.914529
40:     6966.2595: 0.0133891 0.0458635 0.0257932 1.00000e-08 1.00000e-08 0.914517
41:     6966.2595: 0.0133891 0.0458610 0.0257938 1.00000e-08 1.00000e-08 0.914519

Final Estimate of the Negative LLH:
LLH:  6969.51    norm LLH:  1.227891
omega     alpha1     alpha2     alpha3     alpha4      beta1
0.01340442 0.04586097 0.02579384 0.00000001 0.00000001 0.91451944

R-optimhess Difference Approximated Hessian Matrix:
omega    alpha1    alpha2    alpha3    alpha4     beta1
omega  -1547018.7 -529529.8 -537948.8 -538943.1 -549006.4 -755728.9
alpha1  -529529.8 -299566.9 -296707.6 -295346.1 -300555.7 -351582.5
alpha2  -537948.8 -296707.6 -301397.7 -299505.1 -304762.0 -357220.9
alpha3  -538943.1 -295346.1 -299505.1 -306415.6 -312032.4 -363256.3
alpha4  -549006.4 -300555.7 -304762.0 -312032.4 -321365.9 -372529.4
beta1   -755728.9 -351582.5 -357220.9 -363256.3 -372529.4 -476914.3
attr(,"time")
Time difference of 0.2420168 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.833076 secs

Title:
GARCH Modelling

Call:
garchFit(formula = ~garch(4, 1), data = y, include.mean = FALSE)

Mean and Variance Equation:
data ~ garch(4, 1)
<environment: 0x0000000024f721b8>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega      alpha1      alpha2      alpha3      alpha4       beta1
0.01340442  0.04586097  0.02579384  0.00000001  0.00000001  0.91451944

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega  1.340e-02   1.860e-03    7.208 5.68e-13 ***
alpha1 4.586e-02   1.160e-02    3.955 7.66e-05 ***
alpha2 2.579e-02   1.532e-02    1.683   0.0923 .
alpha3 1.000e-08   2.043e-02    0.000   1.0000
alpha4 1.000e-08   1.735e-02    0.000   1.0000
beta1  9.145e-01   6.801e-03  134.472  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6969.51    normalized:  -1.227891

Description:
Thu Jul 05 09:39:33 2018 by user: fanj5

Title:
GARCH Modelling

Call:
garchFit(formula = ~garch(1, 1), data = y, cond.dist = "std",
include.mean = FALSE, trace = F)

Mean and Variance Equation:
data ~ garch(1, 1)
<environment: 0x0000000026709118>
[data = y]

Conditional Distribution:
std

Coefficient(s):
omega    alpha1     beta1     shape
0.017738  0.080087  0.903693  3.915327

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega   0.017738    0.003114    5.697 1.22e-08 ***
alpha1  0.080087    0.009024    8.875  < 2e-16 ***
beta1   0.903693    0.009440   95.730  < 2e-16 ***
shape   3.915327    0.215269   18.188  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6553.022    normalized:  -1.154514

Description:
Thu Jul 05 09:39:33 2018 by user: fanj5

In [19]:
# Advanced ARCH and GARCH estimation in R
# Listing 2.3
# Last updated June 2018
#
#

## normal APARCH(1,1)
print(garchFit(~ aparch(1,1),data=y,include.mean=FALSE,trace=F))
## fixing delta at 2 (or to any value)
print(garchFit(~ aparch(1,1),data=y,include.mean=FALSE,trace=F,include.delta=F,delta=2))
## Student-t conditional distribution
print(garchFit(~ aparch(1,1),data=y,include.mean=FALSE,cond.dist="std",trace=F))
## normal APARCH(2,2)
print(garchFit(~ aparch(2,2),data=y,include.mean=FALSE,trace=F))

Title:
GARCH Modelling

Call:
garchFit(formula = ~aparch(1, 1), data = y, include.mean = FALSE,
trace = F)

Mean and Variance Equation:
data ~ aparch(1, 1)
<environment: 0x00000000256b9040>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega     alpha1     gamma1      beta1      delta
0.0130295  0.0680740  0.0093648  0.9181769  2.0000000

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega   0.013029    0.001691    7.707 1.29e-14 ***
alpha1  0.068074    0.007012    9.708  < 2e-16 ***
gamma1  0.009365    0.031768    0.295    0.768
beta1   0.918177    0.005419  169.436  < 2e-16 ***
delta   2.000000    0.174164   11.483  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6970.276    normalized:  -1.228026

Description:
Thu Jul 05 09:39:36 2018 by user: fanj5

Title:
GARCH Modelling

Call:
garchFit(formula = ~aparch(1, 1), data = y, delta = 2, include.mean = FALSE,
include.delta = F, trace = F)

Mean and Variance Equation:
data ~ aparch(1, 1)
<environment: 0x0000000028066d38>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega     alpha1     gamma1      beta1
0.0130295  0.0680740  0.0093648  0.9181769

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega   0.013029    0.001682    7.748 9.33e-15 ***
alpha1  0.068074    0.005103   13.340  < 2e-16 ***
gamma1  0.009365    0.031726    0.295    0.768
beta1   0.918177    0.005302  173.182  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6970.276    normalized:  -1.228026

Description:
Thu Jul 05 09:39:36 2018 by user: fanj5

Title:
GARCH Modelling

Call:
garchFit(formula = ~aparch(1, 1), data = y, cond.dist = "std",
include.mean = FALSE, trace = F)

Mean and Variance Equation:
data ~ aparch(1, 1)
<environment: 0x0000000027fcd1e0>
[data = y]

Conditional Distribution:
std

Coefficient(s):
omega     alpha1     gamma1      beta1      delta      shape
0.017701   0.079995  -0.049126   0.903756   2.000000   3.912207

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega   0.017701    0.003124    5.667 1.45e-08 ***
alpha1  0.079995    0.011125    7.191 6.45e-13 ***
gamma1 -0.049126    0.042062   -1.168    0.243
beta1   0.903756    0.009674   93.421  < 2e-16 ***
delta   2.000000    0.222578    8.986  < 2e-16 ***
shape   3.912207    0.214974   18.198  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6552.334    normalized:  -1.154393

Description:
Thu Jul 05 09:39:37 2018 by user: fanj5

Title:
GARCH Modelling

Call:
garchFit(formula = ~aparch(2, 2), data = y, include.mean = FALSE,
trace = F)

Mean and Variance Equation:
data ~ aparch(2, 2)
<environment: 0x00000000281219a8>
[data = y]

Conditional Distribution:
norm

Coefficient(s):
omega     alpha1     alpha2     gamma1     gamma2      beta1      beta2
0.024472   0.052137   0.077632   0.135524  -0.086024   0.017708   0.825545
delta
2.000000

Std. Errors:
based on Hessian

Error Analysis:
Estimate  Std. Error  t value Pr(>|t|)
omega   0.024472    0.003233    7.570 3.73e-14 ***
alpha1  0.052137    0.008303    6.280 3.39e-10 ***
alpha2  0.077632    0.009063    8.566  < 2e-16 ***
gamma1  0.135524    0.063229    2.143   0.0321 *
gamma2 -0.086024    0.047027   -1.829   0.0674 .
beta1   0.017708    0.026265    0.674   0.5002
beta2   0.825545    0.025758   32.050  < 2e-16 ***
delta   2.000000    0.175473   11.398  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
-6961.148    normalized:  -1.226418

Description:
Thu Jul 05 09:39:39 2018 by user: fanj5



# Chapter 3: Multivariate Volatility Models¶

• 3.3: EWMA estimation
• 3.5: OGARCH estimation
• 3.7: DCC estimation
• 3.9: Comparison of EWMA, OGARCH, DCC
In [2]:
# Download stock prices in R
# Listing 3.1
# Last updated 2011
#
#

library(tseries)
library(zoo)

y = diff(log(p))*100       # 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 = length(y[,1])

In [3]:
# EWMA in R
# Listing 3.3
# Last updated August 2016
#
#

## 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:TT){                # loop though sample
S = lambda*S+(1-lambda)*t(y[i-1]) %*% y[i-1]
EWMA[i,] = c(S)[c(1,4,2)]     # convert matrix to vector
}
EWMArho = EWMA[,3]/sqrt(EWMA[,1]*EWMA[,2]) # calculate correlations

print(tail(EWMArho))

[1] 0.2296842 0.2236347 0.2218277 0.2270681 0.2054204 0.1992979
[1] 0.1931994 0.1914987 0.2177147 0.3436632 0.3122741 0.3159703

In [4]:
# OGARCH in R
# Listing 3.5
# Last updated 2011
#
#

library(gogarch)

res = gogarch(y,formula = ~garch(1,1),garchlist = c(include.mean=FALSE))
OOrho = ccor(res)

## OOrho is a vector of correlations

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(1, 1)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               1 1
Max GARCH Order:           1
Maximum Order:             1
Conditional Dist:          norm
h.start:                   2
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              1.547425

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -4.100332e-18 4.100332e-18    0.0    FALSE
omega   1.000000e-06 1.000000e+02    0.1     TRUE
alpha1  1.000000e-08 1.000000e+00    0.1     TRUE
gamma1 -1.000000e+00 1.000000e+00    0.1    FALSE
beta1   1.000000e-08 1.000000e+00    0.8     TRUE
delta   0.000000e+00 2.000000e+00    2.0    FALSE
skew    1.000000e-01 1.000000e+01    1.0    FALSE
shape   1.000000e+00 1.000000e+01    4.0    FALSE
Index List of Parameters to be Optimized:
omega alpha1  beta1
2      3      5
Persistence:                  0.9

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     7394.2603: 0.100000 0.100000 0.800000
1:     7365.8037: 0.0803059 0.100888 0.792824
2:     7347.8635: 0.0795917 0.117609 0.805475
3:     7321.6535: 0.0413903 0.134177 0.810643
4:     7298.5740: 0.0319866 0.144851 0.850116
5:     7285.1600: 0.0102903 0.122105 0.877909
6:     7268.7501: 0.0229233 0.0850644 0.893040
7:     7264.9496: 0.0167475 0.0802718 0.903717
8:     7264.7320: 0.0168769 0.0808174 0.904203
9:     7264.5623: 0.0161560 0.0809933 0.904186
10:     7264.2379: 0.0158076 0.0821189 0.905089
11:     7264.1123: 0.0145297 0.0826200 0.905342
12:     7263.9872: 0.0145185 0.0840036 0.905526
13:     7263.9680: 0.0140866 0.0838013 0.905872
14:     7263.9632: 0.0142765 0.0841415 0.905430
15:     7263.9601: 0.0143946 0.0839492 0.905329
16:     7263.9595: 0.0144066 0.0837914 0.905519
17:     7263.9591: 0.0143520 0.0836056 0.905672
18:     7263.9591: 0.0143828 0.0837819 0.905502
19:     7263.9591: 0.0143892 0.0837837 0.905506
20:     7263.9591: 0.0143906 0.0837584 0.905515
21:     7263.9590: 0.0143849 0.0837265 0.905558
22:     7263.9590: 0.0143739 0.0836913 0.905596
23:     7263.9590: 0.0143724 0.0836782 0.905612
24:     7263.9590: 0.0143724 0.0836788 0.905611

Final Estimate of the Negative LLH:
LLH:  9742.057    norm LLH:  1.71636
omega     alpha1      beta1
0.03441505 0.08367885 0.90561099

R-optimhess Difference Approximated Hessian Matrix:
omega    alpha1     beta1
omega  -173124.0 -166164.6 -232610.5
alpha1 -166164.6 -238704.5 -287921.6
beta1  -232610.5 -287921.6 -382203.1
attr(,"time")
Time difference of 0.03500414 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.2190139 secs

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(1, 1)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               1 1
Max GARCH Order:           1
Maximum Order:             1
Conditional Dist:          norm
h.start:                   2
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              0.6862081

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.610204e-16 1.610204e-16    0.0    FALSE
omega   1.000000e-06 1.000000e+02    0.1     TRUE
alpha1  1.000000e-08 1.000000e+00    0.1     TRUE
gamma1 -1.000000e+00 1.000000e+00    0.1    FALSE
beta1   1.000000e-08 1.000000e+00    0.8     TRUE
delta   0.000000e+00 2.000000e+00    2.0    FALSE
skew    1.000000e-01 1.000000e+01    1.0    FALSE
shape   1.000000e+00 1.000000e+01    4.0    FALSE
Index List of Parameters to be Optimized:
omega alpha1  beta1
2      3      5
Persistence:                  0.9

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     7725.2866: 0.100000 0.100000 0.800000
1:     7722.7430: 0.0906144 0.100515 0.798493
2:     7717.7086: 0.0890546 0.107219 0.805069
3:     7711.8394: 0.0718987 0.112195 0.811658
4:     7705.1823: 0.0656395 0.111399 0.829622
5:     7683.1700: 0.0375483 0.0835911 0.878395
6:     7682.8651: 0.0379565 0.0841896 0.880197
7:     7682.1325: 0.0364673 0.0830985 0.880799
8:     7681.1254: 0.0355957 0.0822047 0.884477
9:     7673.8776: 0.0191189 0.0496766 0.927414
10:     7671.5877: 0.0199738 0.0507118 0.928342
11:     7671.0284: 0.0185256 0.0514562 0.928234
12:     7670.2058: 0.0171226 0.0525377 0.930975
13:     7669.9692: 0.0143792 0.0524118 0.932739
14:     7669.8295: 0.0143091 0.0529748 0.933629
15:     7669.6873: 0.0147430 0.0529735 0.932667
16:     7669.6308: 0.0158092 0.0535238 0.930921
17:     7669.6250: 0.0156257 0.0535815 0.930832
18:     7669.6036: 0.0157060 0.0537751 0.930802
19:     7669.5844: 0.0156238 0.0540582 0.930498
20:     7669.5581: 0.0158414 0.0546699 0.929953
21:     7669.5426: 0.0159019 0.0551397 0.929251
22:     7669.5330: 0.0160279 0.0557897 0.928722
23:     7669.5310: 0.0158553 0.0557384 0.928911
24:     7669.5310: 0.0159022 0.0557496 0.928843
25:     7669.5310: 0.0159051 0.0557528 0.928847
26:     7669.5310: 0.0158849 0.0557380 0.928883
27:     7669.5310: 0.0158817 0.0557444 0.928882
28:     7669.5310: 0.0158816 0.0557403 0.928885

Final Estimate of the Negative LLH:
LLH:  5532.095    norm LLH:  0.9746467
omega      alpha1       beta1
0.007478336 0.055740288 0.928885329

R-optimhess Difference Approximated Hessian Matrix:
omega     alpha1      beta1
omega  -4690966 -1292219.9 -1616284.4
alpha1 -1292220  -480341.6  -527426.7
beta1  -1616284  -527426.7  -632146.9
attr(,"time")
Time difference of 0.111011 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.5050211 secs

In [5]:
# DCC in R
# Listing 3.7
# Last updated 2011
#
#

library(ccgarch)

## estimate univariate GARCH models to get starting values

f1 = garchFit(~ garch(1,1), data=y[,1],include.mean=FALSE)
f1 = f1@fit$coef f2 = garchFit(~ garch(1,1), data=y[,2],include.mean=FALSE) f2 = f2@fit$coef

## create vectors and matrices of starting values
a = c(f1[1], f2[1])
A = diag(c(f1[2],f2[2]))
B = diag(c(f1[3], f2[3]))
dccpara = c(0.2,0.6)

## estimate the model
dccresults=dcc.estimation(inia=a,iniA=A,iniB=B,ini.dcc=dccpara,dvar=y,model="diagonal")

## Parameter estimates and their robust standard errors in dcc.results$out DCCrho = dccresults$DCC[,2]

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(1, 1)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               1 1
Max GARCH Order:           1
Maximum Order:             1
Conditional Dist:          norm
h.start:                   2
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              0.6426198

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.669935e-16 1.669935e-16    0.0    FALSE
omega   1.000000e-06 1.000000e+02    0.1     TRUE
alpha1  1.000000e-08 1.000000e+00    0.1     TRUE
gamma1 -1.000000e+00 1.000000e+00    0.1    FALSE
beta1   1.000000e-08 1.000000e+00    0.8     TRUE
delta   0.000000e+00 2.000000e+00    2.0    FALSE
skew    1.000000e-01 1.000000e+01    1.0    FALSE
shape   1.000000e+00 1.000000e+01    4.0    FALSE
Index List of Parameters to be Optimized:
omega alpha1  beta1
2      3      5
Persistence:                  0.9

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     7785.0138: 0.100000 0.100000 0.800000
1:     7782.7078: 0.0911728 0.100995 0.801500
2:     7778.9615: 0.0886434 0.105868 0.808642
3:     7773.6590: 0.0729607 0.108289 0.817176
4:     7767.9361: 0.0661686 0.106552 0.833774
5:     7752.5835: 0.0449532 0.0815642 0.871867
6:     7749.6566: 0.0379184 0.0820200 0.884751
7:     7748.6911: 0.0267339 0.0792960 0.893871
8:     7743.7402: 0.0287565 0.0687111 0.903849
9:     7743.5044: 0.0274129 0.0646857 0.906726
10:     7742.7252: 0.0255079 0.0638911 0.911419
11:     7742.6496: 0.0225478 0.0622858 0.915285
12:     7742.5207: 0.0231552 0.0585416 0.918734
13:     7742.5193: 0.0227706 0.0583703 0.918441
14:     7742.4685: 0.0229251 0.0585204 0.918580
15:     7742.4524: 0.0226959 0.0589777 0.918545
16:     7742.4512: 0.0226569 0.0589680 0.918503
17:     7742.4501: 0.0226838 0.0590145 0.918481
18:     7742.4484: 0.0226899 0.0590679 0.918378
19:     7742.4451: 0.0227854 0.0592013 0.918214
20:     7742.4411: 0.0229070 0.0593988 0.917811
21:     7742.4381: 0.0230717 0.0596895 0.917488
22:     7742.4367: 0.0229720 0.0598279 0.917437
23:     7742.4364: 0.0230546 0.0598621 0.917283
24:     7742.4362: 0.0230995 0.0599810 0.917158
25:     7742.4362: 0.0230828 0.0599360 0.917211
26:     7742.4362: 0.0230828 0.0599363 0.917210

Final Estimate of the Negative LLH:
LLH:  5232.497    norm LLH:  0.9218635
omega     alpha1      beta1
0.00953228 0.05993628 0.91721034

R-optimhess Difference Approximated Hessian Matrix:
omega     alpha1      beta1
omega  -4089457 -1029107.8 -1299358.6
alpha1 -1029108  -357736.4  -387217.8
beta1  -1299359  -387217.8  -465578.8
attr(,"time")
Time difference of 0.0467999 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.206605 secs

Series Initialization:
ARMA Model:                arma
Formula Mean:              ~ arma(0, 0)
GARCH Model:               garch
Formula Variance:          ~ garch(1, 1)
ARMA Order:                0 0
Max ARMA Order:            0
GARCH Order:               1 1
Max GARCH Order:           1
Maximum Order:             1
Conditional Dist:          norm
h.start:                   2
llh.start:                 1
Length of Series:          5676
Recursion Init:            mci
Series Scale:              0.9482962

Parameter Initialization:
Initial Parameters:          $params Limits of Transformations:$U, $V Which Parameters are Fixed?$includes
Parameter Matrix:
U            V params includes
mu     -1.405184e-17 1.405184e-17    0.0    FALSE
omega   1.000000e-06 1.000000e+02    0.1     TRUE
alpha1  1.000000e-08 1.000000e+00    0.1     TRUE
gamma1 -1.000000e+00 1.000000e+00    0.1    FALSE
beta1   1.000000e-08 1.000000e+00    0.8     TRUE
delta   0.000000e+00 2.000000e+00    2.0    FALSE
skew    1.000000e-01 1.000000e+01    1.0    FALSE
shape   1.000000e+00 1.000000e+01    4.0    FALSE
Index List of Parameters to be Optimized:
omega alpha1  beta1
2      3      5
Persistence:                  0.9

--- START OF TRACE ---
Selected Algorithm: nlminb

R coded nlminb Solver:

0:     7308.8679: 0.100000 0.100000 0.800000
1:     7268.3526: 0.0784678 0.100994 0.791633
2:     7244.9928: 0.0766506 0.119819 0.804935
3:     7228.3507: 0.0342297 0.137620 0.809635
4:     7200.9599: 0.0287924 0.156200 0.851632
5:     7172.8137: 0.00853985 0.126221 0.880435
6:     7159.9675: 0.0214295 0.109813 0.877659
7:     7159.9243: 0.0176985 0.0921549 0.888492
8:     7157.1129: 0.0177346 0.0882042 0.898246
9:     7156.1513: 0.0112938 0.0869484 0.906475
10:     7155.8125: 0.0136246 0.0945070 0.899532
11:     7155.6586: 0.0136505 0.0961369 0.895966
12:     7155.2713: 0.0141273 0.0930951 0.898394
13:     7155.2709: 0.0139671 0.0919137 0.898938
14:     7155.2341: 0.0141919 0.0915743 0.899452
15:     7155.2259: 0.0142920 0.0909305 0.899522
16:     7155.2252: 0.0142645 0.0909255 0.899624
17:     7155.2250: 0.0142567 0.0909970 0.899547
18:     7155.2249: 0.0142851 0.0910750 0.899481
19:     7155.2247: 0.0142181 0.0911105 0.899555
20:     7155.2247: 0.0142468 0.0911595 0.899466
21:     7155.2247: 0.0142497 0.0911605 0.899468
22:     7155.2247: 0.0142487 0.0911571 0.899468
23:     7155.2247: 0.0142460 0.0911550 0.899475
24:     7155.2247: 0.0142483 0.0911483 0.899476
25:     7155.2247: 0.0142463 0.0911444 0.899481
26:     7155.2247: 0.0142465 0.0911447 0.899482
27:     7155.2247: 0.0142462 0.0911446 0.899482
28:     7155.2247: 0.0142458 0.0911447 0.899482
29:     7155.2247: 0.0142452 0.0911400 0.899487
30:     7155.2247: 0.0142444 0.0911405 0.899488

Final Estimate of the Negative LLH:
LLH:  6853.895    norm LLH:  1.207522
omega     alpha1      beta1
0.01280948 0.09114050 0.89948793

R-optimhess Difference Approximated Hessian Matrix:
omega    alpha1     beta1
omega  -1229502.2 -404258.1 -584098.8
alpha1  -404258.1 -203252.1 -251079.7
beta1   -584098.8 -251079.7 -343472.3
attr(,"time")
Time difference of 0.03119993 secs

--- END OF TRACE ---

Time to Estimate Parameters:
Time difference of 0.2292051 secs
****************************************************************
*  Estimation has been completed.                              *
*  The outputs are saved in a list with components:            *
*    out    : the estimates and their standard errors          *
*    loglik : the value of the log-likelihood at the estimates *
*    h      : a matrix of estimated conditional variances      *
*    DCC    : a matrix of DCC estimates                        *
*    std.resid : a matrix of the standardised residuals        *
*    first  : the results of the first stage estimation        *
*    second : the results of the second stage estimation       *
****************************************************************

In [6]:
# Sample statistics in R
# Listing 3.9
# Last updated 2011
#
#

matplot(cbind(EWMArho,DCCrho,OOrho),type='l',las=1,lty=1:3,col=1:3,ylab="")
mtext("Correlations",side=2,line=0.3,at=1,las=1,cex=0.8)
legend("bottomright",c("EWMA","DCC","OO"),lty=1:3,col=1: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)

[1] 0.7978846 1.7549833 2.0627128 2.3378028 2.6652142 3.3670901


# Chapter 5: Implementing Risk Forecasts¶

• 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 [26]:
# Download stock prices in R
# Listing 5.1
# Last updated August 2016
#
#

library(tseries)
library(zoo)

## convert prices of first two stocks to returns and adjust length
y1=tail(diff(log(coredata(prices[,1]))),4100)
y2=tail(diff(log(coredata(prices[,2]))),4100)
TT=length(y1)
y=cbind(y1,y2)
value = 1000 # portfolio value
p = 0.01     # probability

In [27]:
# Univariate HS in R
# Listing 5.3
# Last updated August 2016
#
#

ys = sort(y1) # sort returns
op = TT*p     # p percent smallest
VaR1 = -ys[op]*value

print(VaR1)

[1] 18.65861

In [28]:
# Multivariate HS in R
# Listing 5.5
# Last updated 2011
#
#

w = matrix(c(0.3,0.7)) # vector of portfolio weights
yp = y %*% w           # obtain portfolio returns
yps = sort(yp)
VaR2 = -yps[op]*value

print(VaR2)

[1] 19.81446

In [29]:
# Univariate ES in R
# Listing 5.7
# Last updated 2011
#
#

ES1 = -mean(ys[1:op])*value

print(ES1)

[1] 23.87292

In [30]:
# Normal VaR in R
# Listing 5.9
# Last updated 2011
#
#

sigma = sd(y1) # estimate volatility
VaR3 = -sigma * qnorm(p) * value

print(VaR3)

[1] 15.588

In [31]:
# Portfolio normal VaR in R
# Listing 5.11
# Last updated 2011
#
#

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)*value

print(VaR4)

[1] 18.13829

In [32]:
# Student-t VaR in R
# Listing 5.13
# Last updated August 2016
#
#

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) *  value

print(VaR5)

   sigma
17.88518

In [33]:
# Normal ES in R
# Listing 5.15
# Last updated June 2018
#
#

sigma = sd(y1)
ES2 = sigma*dnorm(qnorm(p))/p * value

print(ES2)

[1] 17.85862

In [34]:
# Direct integration ES in R
# Listing 5.17
# Last updated 2011
#
#

VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$value/p*value print(ES)  [1] 17.85862  In [35]: # MA normal VaR in R # Listing 5.19 # Last updated June 2018 # # WE=20 for (t in seq(TT-5,TT)){ t1=t-WE+1 window= y1[t1:t] # estimation window sigma=sd(window) VaR6 = -sigma * qnorm(p) * value print(VaR6) }  [1] 16.0505 [1] 16.1491 [1] 18.85435 [1] 18.88212 [1] 16.23053 [1] 16.16976  In [36]: # EWMA VaR in R # Listing 5.21 # Last updated August 2016 # # lambda = 0.94; s11 = var(y1[1:30]); # initial variance for (t in 2:TT){ s11 = lambda * s11 + (1-lambda) * y1[t-1]^2 } VaR7 = -qnorm(p) * sqrt(s11) * value print(VaR7)  [1] 16.75344  In [37]: # Two-asset EWMA VaR in R # Listing 5.23 # Last updated 2011 # # s = cov(y) # initial covariance for (t in 2:TT){ 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) * value print(VaR8)  [1] 20.50363  In [38]: # GARCH in R # Listing 5.25 # Last updated 2011 # # library(fGarch) g = garchFit(~garch(1,1),y1,include.mean=F,trace=F) omega = g@fit$matcoef[1,1]
alpha = g@fit$matcoef[2,1] beta = g@fit$matcoef[3,1]
sigma2 = omega + alpha*y[TT]^2 + beta*g@h.t[TT] # calc sigma2 for t+1
VaR9 = -sqrt(sigma2) * qnorm(p) * value

print(VaR9)

[1] 17.05734


# 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 [39]:
# 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)

$Call [1] 13.49852$Put
[1] 1.27641

$Delta.Call [1] 0.8395228$Delta.Put
[1] -0.1604772

$Gamma [1] 0.01723826  # 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)  1. 0.0693609158042818 2. 0.817775198724121 3. 0.942621732363477 4. 0.269381876336411 5. 0.169348123250529 6. 0.0338956224732101 7. 0.178785004187375 8. 0.641665365546942 9. 0.0228777434676886 10. 0.00832482660189271 1. -0.27229604424923 2. -0.315348711467784 3. -0.628255236517538 4. -0.106463884872094 5. 0.428014802202354 6. -0.777719582147445 7. -1.29388229761362 8. -0.779566508059429 9. 0.0119517585652984 10. -0.152416237822168 1. -0.54686521181825 2. 0.325766227984669 3. -0.310344756918247 4. 1.64011963529838 5. -0.60065788608405 6. -0.504397302337928 7. 0.146674867812759 8. 0.421606955584244 9. -1.10969020069546 10. 0.674102188748512 In [43]: # Price bond in R # Listing 7.5 # Last updated August 2016 # # yield=c(5.00, 5.69, 6.09, 6.38, 6.61, 6.79, 6.94, 7.07, 7.19, 7.30) # yield curve TT=length(yield) 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:TT))) # calculate price print(P)  [1] 9.913206  In [44]: # Simulate yields in R # Listing 7.7 # Last updated August 2016 # # 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=TT,ncol=S) for (i in 1:S) ysim[,i]=yield+r[i] matplot(ysim,type='l')  In [45]: # Simulate bond prices in R # Listing 7.9 # Last updated August 2016 # # SP = vector(length=S) for (i in 1:S){ # S simulations SP[i] = sum(cc/((1+ysim[,i]/100)^(TT))) } 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)))  In [46]: # Black-Scholes valuation in R # Listing 7.11 # Last updated August 2016 # # P0 = 50 # initial spot price sigma = 0.2 # annual volatility r = 0.05 # annual interest TT = 0.5 # time to expiration X = 40 # strike price f = bs(X,P0,r,sigma,TT) # analytical call price ## this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2) print(f)  $Call
[1] 11.08728

$Put [1] 0.09967718$Delta.Call
[1] 0.9660259

$Delta.Put [1] -0.03397407$Gamma
[1] 0.01066378


In [47]:
# Black-Scholes simulation in R
# Listing 7.13
# Last updated August 2016
#
#

set.seed(12)                      # set seed

S = 1e6                           # number of simulations
F = P0*exp(r*TT)                  # futures price
ysim = rnorm(S,-0.5*sigma^2*TT,sigma*sqrt(TT)) # sim returns, lognorm corrected
F=F*exp(ysim)                     # sim futures price
SP = F-X                          # payoff
SP[SP<0] = 0                      # set negative outcomes to zero
fsim = SP*exp(-r*TT)              # discount
call_sim = mean(fsim)             # simulated price

print(call_sim)

[1] 11.08709

In [48]:
# Option density plots in R
# Listing 7.15
# Last updated 2011
#
#

par(mfrow=c(1,2), pty="s")

hist(F,probability=TRUE,ylim=c(0,0.06))

x=seq(min(F),max(F),length=100)
lines(x, dnorm(x, mean = mean(F), sd = sd(SP)))
hist(fsim,nclass=100,probability=TRUE)

In [49]:
# Simulate VaR in R
# Listing 7.17
# Last updated 2011
#
#

set.seed(1)                 # set seed

S = 1e7                     # number of simulations
s2 = 0.01^2                 # daily variance
p = 0.01                    # probability
r = 0.05                    # annual riskfree rate
P = 100                     # price today
ysim = rnorm(S,r/365-0.5*s2,sqrt(s2)) # sim returns
Psim = P*exp(ysim) # sim future prices
q = sort(Psim-P) # simulated P/L
VaR1 = -q[p*S]

print(VaR1)

[1] 2.289809

In [50]:
# Simulate option VaR in R
# Listing 7.19
# Last updated August 2016
#
#

TT = 0.25;                             # time to expiration
X = 100;                               # strike price
sigma = sqrt(s2*250);                  # annual volatility
f = bs(X,P,r,sigma,TT)                 # analytical call price
fsim = bs(X,Psim,r,sigma,TT-(1/365))   # sim option prices
q = sort(fsim$Call-f$Call)             # simulated P/L
VaR2 = -q[p*S]

print(VaR2)

[1] 1.214623

In [51]:
# Example 7.3 in R
# Listing 7.21
# Last updated August 2016
#
#

X1 = 100
X2 = 110
f1 = bs(X1,P,r,sigma,TT)
f2 = bs(X2,P,r,sigma,TT)
f2sim = bs(X2,Psim,r,sigma,TT-(1/365))
f1sim = bs(X1,Psim,r,sigma,TT-(1/365))
q = sort(f1sim$Call+f2sim$Put+Psim-f1$Call-f2$Put-P);
VaR3 = -q[p*S]

print(VaR3)

[1] 1.494695

In [52]:
# Simulated two-asset returns in R
# Listing 7.23
# Last updated 2011
#
#

library (MASS)

set.seed(12)  # set seed

mu = c(r/365,r/365)                               # return mean
Sigma = matrix(c(0.01,0.0005,0.0005,0.02),ncol=2) # covariance matrix
y = mvrnorm(S,mu,Sigma)  # simulated returns

In [53]:
# Two-asset VaR in R
# Listing 7.25
# Last updated 2011
#
#

K=2
P = c(100,50)                                 # prices
x = c(1,1)                                    # number of assets
Port = P %*% x                                # portfolio at t
Psim = matrix(t(matrix(P,K,S)),ncol=K)*exp(y) # simulated prices
PortSim = Psim %*% x                          # simulated portfolio value
q = sort(PortSim-Port[1,1])                   # simulated P/L
VaR4 = -q[S*p]

print(VaR4)

[1] 25.96307

In [54]:
# A two-asset case in R with an option
# Listing 7.27
# Last updated August 2016
#
#

f = bs(P[2],P[2],r,sigma,TT)
fsim = bs(P[2],Psim[,2],r,sigma,TT-(1/365))
q = sort(fsim$Call+Psim[,1]-f$Call-P[1]);
VaR5 = -q[p*S]

print(VaR5)

[1] 20.79232


# Chapter 8: Backtesting and Stress Testing¶

• 8.3: Setting up backtest
• 8.5: Running backtest for EWMA/MA/HS/GARCH VaR
• 8.7: Backtesting analysis for EWMA/MA/HS/GARCH VaR
• 8.9: Bernoulli coverage test
• 8.11: Independence test
• 8.13: Running Bernoulli/Independence test on backtests
• 8.15: Running backtest for EWMA/HS ES
• 8.17: Backtesting analysis for EWMA/HS ES
In [55]:
# Load data in R
# Listing 8.1
# Last updated August 2016
#
#

library(tseries)
library(zoo)

y=diff(log(p)) # get returns
y=coredata(y)  # strip date information

In [56]:
# Set backtest up in R
# Listing 8.3
# Last updated August 2016
#
#

TT = length(y)               # number of obs for y
WE = 1000                    # estimation window length
p = 0.01                     # probability
l1 = WE*p                    # HS observation
value = 1;                   # portfolio value
VaR = matrix(nrow=TT,ncol=4) # matrix for forecasts
## EWMA setup
lambda = 0.94;
s11 = var(y[1:30]);
for(t in 2:WE) s11=lambda*s11+(1-lambda)*y[t-1]^2

library(fGarch)

In [57]:
# Running backtest in R
# Listing 8.5
# Last updated August 2016
#
#

for (t in (WE+1):TT){
t1 = t-WE;         # start of the data window
t2 = t-1;	         # end of the data window
window = y[t1:t2]  # data for estimation

s11=lambda*s11+(1-lambda)*y[t-1]^2
VaR[t,1] = -qnorm(p) * sqrt(s11) * value          # EWMA

VaR[t,2] = - sd(window) * qnorm(p)*value          # MA

ys = sort(window)
VaR[t,3] = -ys[l1]*value                          # HS

g=garchFit(formula=~garch(1,1),window,
trace=FALSE,include.mean=FALSE)
par=g@fit\$matcoef
s4=par[1]+par[2]*window[WE]^2+par[3]*g@h.t[WE]
VaR[t,4] = -qnorm(p) * sqrt(s4) * value           # GARCH(1,1)
}

In [58]:
# Backtesting analysis in R
# Listing 8.7
# Last updated June 2018
#
#

W1=WE+1
for (i in 1:4){
VR = sum(y[W1:TT]< -VaR[W1:TT,i])/(p*(TT-WE))
s = sd(VaR[W1:TT,i])
cat(i,"VR",VR,"VaR vol",s,"\n")
}
matplot(cbind(y[W1:TT],VaR[W1:TT,]),type='l',col=1:5,las=1,ylab="",lty=1:5)
legend("topleft",legend=c("Returns","EWMA","MA","HS","GARCH"),lty=1:5,col=1:5,bty="n")

1 VR 2.074423 VaR vol 0.01211406
2 VR 1.860565 VaR vol 0.006622265
3 VR 1.240376 VaR vol 0.01229172
4 VR 1.668092 VaR vol 0.01175914

In [59]:
# Bernoulli coverage test in R
# Listing 8.9
# Last updated August 2016
#
#

bern_test=function(p,v){
lv=length(v)
sv=sum(v)

al=log(p)*sv+log(1-p)*(lv-sv)
bl=log(sv/lv)*sv +log(1-sv/lv)*(lv-sv)
return(-2*(al-bl))
}

In [60]:
# Independence test in R
# Listing 8.11
# Last updated June 2018
#
#

ind_test=function(V){
J=matrix(ncol=4,nrow=length(V))
for (i in 2:length(V)){
J[i,1]=V[i-1]==0 & V[i]==0
J[i,2]=V[i-1]==0 & V[i]==1
J[i,3]=V[i-1]==1 & V[i]==0
J[i,4]=V[i-1]==1 & V[i]==1
}
V_00=sum(J[,1],na.rm=TRUE)
V_01=sum(J[,2],na.rm=TRUE)
V_10=sum(J[,3],na.rm=TRUE)
V_11=sum(J[,4],na.rm=TRUE)
p_00=V_00/(V_00+V_01)
p_01=V_01/(V_00+V_01)
p_10=V_10/(V_10+V_11)
p_11=V_11/(V_10+V_11)
hat_p=(V_01+V_11)/(V_00+V_01+V_10+V_11)
al = log(1-hat_p)*(V_00+V_10) + log(hat_p)*(V_01+V_11)
bl = log(p_00)*V_00 + log(p_01)*V_01 + log(p_10)*V_10 + log(p_11)*V_11
return(-2*(al-bl))
}

In [61]:
# Backtesting S&P 500 in R
# Listing 8.13
# Last updated August 2016
#
#

W1=WE+1
ya=y[W1:TT]
VaRa=VaR[W1:TT,]
m=c("EWMA","MA","HS","GARCH")
for (i in 1:4){
q= y[W1:TT]< -VaR[W1:TT,i]
v=VaRa*0
v[q,i]=1
ber=bern_test(p,v[,i])
ind=ind_test(v[,i])
cat(i,m[i],'Bernoulli',ber,1-pchisq(ber,1),"independence",ind,1-pchisq(ind,1),"\n")
}

1 EWMA Bernoulli 41.6257 1.105309e-10 independence 0.4411056 0.5065893
2 MA Bernoulli 27.90392 1.27491e-07 independence 17.44182 2.962367e-05
3 HS Bernoulli 2.535439 0.1113159 independence 11.48482 0.0007016709
4 GARCH Bernoulli 17.55348 2.793385e-05 independence 0.33451 0.5630154

In [62]:
# Backtest ES in R
# Listing 8.15
# Last updated August 2016
#
#

VaR = matrix(nrow=TT,ncol=2) # VaR forecasts for 2 models
ES = matrix(nrow=TT,ncol=2) # ES forecasts for 2 models
for (t in (WE+1):TT){
t1 = t-WE;
t2 = t-1;
window = y[t1:t2]

s11 = lambda * s11  + (1-lambda) * y[t-1]^2
VaR[t,1] = -qnorm(p) * sqrt(s11) * value  # EWMA
ES[t,1] = sqrt(s11) * dnorm(qnorm(p)) / p

ys = sort(window)
VaR[t,2] = -ys[l1]*value                  # HS
ES[t,2] = -mean(ys[1:l1]) * value
}

In [63]:
# Backtest ES in R
# Listing 8.17
# Last updated August 2016
#
#

ESa = ES[W1:TT,]
VaRa = VaR[W1:TT,]
for (i in 1:2){
q = ya <= -VaRa[,i]
nES = mean(ya[q] / -ESa[q,i])
cat(i,"nES",nES,"\n")
}

1 nES 1.223518
2 nES 1.053689


# Chapter 9: Extreme Value Theory¶

• 9.1: Calculation of tail index from returns
In [64]:
# Hill estimator in R
# Listing 9.1
# Last updated 2011
#
#

ysort = sort(y)                             # sort the returns
CT = 100                                    # set the threshold
iota = 1/mean(log(ysort[1:CT]/ysort[CT+1])) # get the tail index

print(iota)

# END

[1] 2.62971