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.
##
## 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
16.25
115.785
9
3.1
  1. 3.1
  2. 9
5.41666666666667
4.15
9.905833
9.905833
1
  1. 3.1
  2. 4.15
  3. 9
1.131402
1.423108
2.197225
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.905833
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
13
24
12
131
242
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.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 [13]:
# Download S&P500 data in R
# Listing 1.1
# Last updated June 2018
#
#

library(tseries)
library(zoo)


 
price = zoo(read.csv('index.csv', header=TRUE, sep=','))
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
#
#

p = ts(read.csv('stocks.csv',header=TRUE,sep=','))

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)

p = zoo(read.csv('index.csv', header=TRUE, sep=','))
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.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 [2]:
# Download stock prices in R
# Listing 3.1
# Last updated 2011
#
#

library(tseries)
library(zoo)

p = zoo(read.csv('stocks.csv',header=TRUE,sep=','))
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(head(EWMArho))
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.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 [26]:
# Download stock prices in R
# Listing 5.1
# Last updated August 2016
#
#

library(tseries)
library(zoo)

prices = zoo(read.csv('stocks.csv',header=TRUE,sep=','))
## 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.1: Loading hypothetical stock prices, converting to returns
  • 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)

p = zoo(read.csv('index.csv',header=TRUE,sep=','))
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

All rights reserved, Jon Danielsson, 2011-2018.