Introduction

This notebook is an implementation of Jón Daníelsson's Financial Risk Forecasting (Wiley, 2011) in MATLAB 2018a, with annotations and introductory examples. The introductory examples (Appendix) are similar to Appendix C in the original book.

The following code has been adapted from the original, implemented in MATLAB 2016a. Occasional lines of code which require different syntax to run in MATLAB 2016a are also noted in trailing comments.

Bullet point numbers correspond to the 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 MATLAB

Created in MATLAB 2018a (June 2018)

  • M.1: Entering and Printing Data
  • M.2: Vectors, Matrices and Sequences
  • M.3: Importing Data (to be updated)
  • M.4: Basic Summary Statistics
  • M.5: Calculating Moments
  • M.6: Basic Matrix Operations
  • M.7: Statistical Distributions
  • M.8: Statistical Tests
  • M.9: Time Series
  • M.10: Loops and Functions
  • M.11: Basic Graphs
  • M.12: Miscellaneous Useful Functions
In [1]:
% Entering and Printing Data
% Listing M.1
% Last updated June 2018
%
%

x = 10;             % assign x the value 10, silencing output print with ;
disp(x)             % display x
    10


In [2]:
% Vectors, Matrices and Sequences
% Listing M.2
% Last updated June 2018
%
%

y = [1,3,5,7,9]            % lists are denoted by square brackets

y(3)                       % calling 3rd element (MATLAB indices start at 1)

size(y)                    % shows that y is 1 x 5 (a row vector, by default)

length(y)                  % as expected, y has length 5

v = nan(2,3)               % fill a 2 x 3 matrix with NaN values

size(v)                    % as expected, v is size (2,3)

w = repmat([1,2,3]', 2, 3) % repeats matrix twice by rows, thrice by columns

s = 1:10                   % s is a list of integers from 1 to 10 inclusive
y =

     1     3     5     7     9


ans =

     5


ans =

     1     5


ans =

     5


v =

   NaN   NaN   NaN
   NaN   NaN   NaN


ans =

     2     3


w =

     1     1     1
     2     2     2
     3     3     3
     1     1     1
     2     2     2
     3     3     3


s =

     1     2     3     4     5     6     7     8     9    10


In [3]:
% Importing Data
% Listing M.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 csvread to read it
%%
%% Example:
%% data = csvread('data.csv', 1, 0);
%% the two numbers behind are the row offset and column offset
%% so here we ignore the first row (ie. the header)

In [4]:
% Basic Summary Statistics
% Listing M.4
% Last updated June 2018
%
%

y = [3.14,15,9.26,5];

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
corrcoef(y)    % corr matrix = [1] for single vector
sort(y)        % sorting in ascending order
log(y)         % natural log
ans =

   32.4000


ans =

   2.1807e+03


ans =

    15


ans =

    3.1400


ans =

   11.8600


ans =

    8.1000


ans =

    7.1300


ans =

   27.7224


ans =

   27.7224


ans =

     1


ans =

    3.1400    5.0000    9.2600   15.0000


ans =

    1.1442    2.7081    2.2257    1.6094


In [5]:
% Calculating Moments
% Listing M.5
% Last updated June 2018
%
%

mean(y)      % mean
var(y)       % variance
std(y)       % unbiased standard deviation, by default
skewness(y)  % skewness
kurtosis(y)  % kurtosis
ans =

    8.1000


ans =

   27.7224


ans =

    5.2652


ans =

    0.4700


ans =

    1.7153


In [6]:
% Basic Matrix Operations
% Listing M.6
% Last updated June 2018
%
%

z = [1, 2; 3, 4]  % z is a 2 x 2 matrix (Note the use of ; as row separator)
x = [1, 2]        % x is a 1 x 2 matrix

%% Note: z * x is undefined since the two matrices are not conformable

z * x'            % this evaluates to a 2 x 1 matrix

vertcat(z,x)      % "stacking" z and x vertically
horzcat(z,x')     % "stacking z and x' horizontally

%% Note: dimensions must match along the combining axis)
z =

     1     2
     3     4


x =

     1     2


ans =

     5
    11


ans =

     1     2
     3     4
     1     2


ans =

     1     2     1
     3     4     2


In [7]:
% Statistical Distributions
% Listing M.7
% Last updated June 2018
%
%


q = -3:1:3                       % specify a set of values

p = 0.1:0.1:0.9                  % specify a set of probabilities

norminv(p, 0, 1)                 % element-wise inverse Normal quantile

tcdf(q, 4)                       % element-wise cdf under Student-t(4)

chi2pdf(q, 2)                    % element-wise pdf under Chisq(2)

%% One can also obtain pseudorandom samples from distributions

x = trnd(5, 100, 1);             % Sampling 100 times from t dist with 5 df

y = normrnd(0, 1, 100, 1);       % Sampling 50 times from a standard normal 

%% Given sample data, we can also obtain MLE estimates of distribution parameters:

res = fitdist(x, "Normal")       % Fitting x to normal dist
q =

    -3    -2    -1     0     1     2     3


p =

    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000


ans =

   -1.2816   -0.8416   -0.5244   -0.2533         0    0.2533    0.5244    0.8416    1.2816


ans =

    0.0200    0.0581    0.1870    0.5000    0.8130    0.9419    0.9800


ans =

         0         0         0    0.5000    0.3033    0.1839    0.1116


res = 

  NormalDistribution

  Normal distribution
       mu = 0.146605   [-0.155625, 0.448836]
    sigma =  1.52317   [1.33736, 1.76943]


In [8]:
% Statistical Tests
% Listing M.8
% Last updated June 2018
%
%

x = trnd(5, 500, 1);                    % Create hypothetical dataset x

[h1, p1, jbstat] = jbtest(x)            % Jarque-Bera test for normality

[h2, p2, lbstat] = lbqtest(x,'lags',20) % Ljung-Box test for serial correlation
Warning: P is less than the smallest tabulated value, returning 0.001.
> In jbtest (line 136)

h1 =

     1


p1 =

   1.0000e-03


jbstat =

  398.3104


h2 =

  logical

   0


p2 =

    0.1588


lbstat =

   26.2154


In [9]:
% Time Series
% Listing M.9
% Last updated June 2018
%
%

x = trnd(5, 60, 1); % Create hypothetical dataset x

subplot(1,2,1)
autocorr(x, 20)     % autocorrelation for lags 1:20
subplot(1,2,2)
parcorr(x,20)       % partial autocorrelation for lags 1:20

In [10]:
% Loops and Functions
% Listing M.10
% Last updated June 2018
%
%

%% For loops

for i = 3:7        % iterates through [3,4,5,6,7]
    i^2  
end

%% If-else loops

X = 10;

if (rem(X,3) == 0)
    disp("X is a multiple of 3")
else 
    disp("X is not a multiple of 3")
end
    
%% Functions (example: a simple excess kurtosis function)
%% NOTE: in MATLAB, functions can be defined in 2 locations:
%% 1) in a separate file (e.g. excess_kurtosis.m in this case) in the workspace
%% 2) in the same file as the rest of the code, BUT at the end of the file

%% function k = excess_kurtosis(x, excess) 
%%     if nargin == 1                      % if there is only 1 argument
%%         excess = 3;                     % set excess = 3
%%     end                                 % this is how optional param excess is set
%%     m4 = mean((x-mean(x)).^4);
%%     k = m4/(std(x)^4) - excess;
%% end
ans =

     9


ans =

    16


ans =

    25


ans =

    36


ans =

    49

X is not a multiple of 3

In [11]:
% Basic Graphs
% Listing M.11
% Last updated June 2018
%
%

y = normrnd(0, 1, 50, 1);
z = trnd(4, 50, 1);

subplot(2,2,1)
bar(y)            % bar plot
subplot(2,2,2)
plot(y)           % line plot
subplot(2,2,3)
histogram(y)      % histogram
subplot(2,2,4)
scatter(y,z)      % scatter plot

In [12]:
% Miscellaneous Useful Functions
% Listing M.12
% Last updated June 2018
%
%

%% Convert objects from one type to another with int8() etc
%% To check type, use isfloat(object), isinteger(object) and so on
x = 8.0;
isfloat(x)
x = int8(x);
isinteger(x)
ans =

  logical

   1


ans =

  logical

   1


Chapter 1: Financial Markets, Prices and Risk

  • 1.1/1.2: Loading hypothetical stock prices, converting to returns, plotting returns
  • 1.3/1.4: Summary statistics for returns timeseries
  • 1.5/1.6: Autocorrelation function (ACF) plots, Ljung-Box test
  • 1.7/1.8: Quantile-Quantile (QQ) plots
  • 1.9/1.10: Correlation matrix between different stocks
In [13]:
% Download S&P 500 data in MATLAB
% Listing 1.2
% Last updated August 2016
%
%

price = csvread('index.csv', 1, 0);
y=diff(log(price)); % calculate returns

plot(y)             % plot returns


% END

In [14]:
% Sample statistics in MATLAB
% Listing 1.4
% Last updated June 2018
%
%

%% the function sacf uses Kevin Sheppard's MFE toolbox
%% download at https://www.kevinsheppard.com/MFE_Toolbox

mean(y)
std(y)
min(y)
max(y)
skewness(y)
kurtosis(y)
sacf(y,1,[],0)
sacf(y.^2,1,[],0)
[h,pValue,stat]=jbtest(y);
[h,pValue,stat]=lbqtest(y,'lags',20);          
[h,pValue,stat]=lbqtest(y.^2,'lags',20);

%% NOTE: in MATLAB 2018a, some functions require name-value pairs
%% e.g. MATLAB 2016a: [h,pValue,stat] = lbqtest(y,20)
ans =

   2.5816e-04


ans =

    0.0100


ans =

   -0.1020


ans =

    0.1067


ans =

    0.1526


ans =

   16.9812


ans =

    0.0148


ans =

    0.1992

Warning: P is less than the smallest tabulated value, returning 0.001.
> In jbtest (line 136)

In [15]:
% ACF plots and the Ljung-Box test in MATLAB
% Listing 1.6
% Last updated August 2016
%
%

%% subplots here are just for ease of visualization

subplot(1,2,1)
autocorr(y, 20)
subplot(1,2,2)
autocorr(y.^2, 20)

In [16]:
% QQ plots in MATLAB
% Listing 1.8
% Last updated 2011
%
%

%% subplots here are just for ease of visualization

subplot(1,2,1)
qqplot(y)
subplot(1,2,2)
qqplot(y, fitdist(y,'tLocationScale'))

In [17]:
% Download stock prices in MATLAB
% Listing 1.10
% Last updated 2011
%
%

price = csvread('stocks.csv', 1, 0);
y=diff(log(price));
corr(y) % correlation matrix
ans =

    1.0000    0.2297    0.2126
    0.2297    1.0000    0.1451
    0.2126    0.1451    1.0000


Chapter 2: Univariate Volatility Modelling

  • 2.1/2.2: GARCH and t-GARCH estimation
  • 2.3/2.4: APARCH estimation
In [18]:
% ARCH and GARCH estimation in MATLAB
% Listing 2.2
% Last updated June 2018
%
%

p = csvread('index.csv', 1, 0);

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

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

tarch(y,1,0,0);              % ARCH(1)
tarch(y,4,0,0);              % ARCH(4)
tarch(y,4,0,1);              % GARCH(4,1)
tarch(y,1,0,1);              % GARCH(1,1)
tarch(y,1,0,1,'STUDENTST');  % t-GARCH(1,1)

In [19]:
% Advanced ARCH and GARCH estimation in MATLAB
% Listing 2.4
% Last updated August 2016
%
%

aparch(y,1,1,1);              % APARCH(1,1)
aparch(y,2,2,1);              % APARCH(2,1)
aparch(y,1,1,1,'STUDENTST');  % t-APARCH(1,1)

Chapter 3: Multivariate Volatility Models

  • 3.1/3.2: Loading hypothetical stock prices
  • 3.3/3.4: EWMA estimation
  • 3.5/3.6: OGARCH estimation (unavailable as of June 2018)
  • 3.7/3.8: DCC estimation (unavailable as of June 2018)
  • 3.9/3.10: Comparison of EWMA, OGARCH, DCC (unavailable as of June 2018)
In [1]:
% Download stock prices in MATLAB
% Listing 3.2
% Last updated August 2016
%
%

p = csvread('stocks.csv',1,0);
p = p(:,[1,2]); % consider first two stocks

y = diff(log(p))*100; % convert prices to returns
y(:,1)=y(:,1)-mean(y(:,1)); % subtract mean
y(:,2)=y(:,2)-mean(y(:,2));
T = length(y);

In [2]:
% EWMA in MATLAB
% Listing 3.4
% Last updated June 2018
%
%

%% create a matrix to hold covariance matrix for each t
EWMA = nan(T,3); 
lambda = 0.94;
S = cov(y);                 % initial (t=1) covar matrix
EWMA(1,:) = S([1,4,2]);     % extract var and covar
for i = 2:T                 % loop though the sample
    S = lambda*S+(1-lambda)* y(i-1,:)'*y(i-1,:);
    EWMA(i,:) = S([1,4,2]); % convert matrix to vector
end
EWMArho = EWMA(:,3)./sqrt(EWMA(:,1).*EWMA(:,2)); % calculate correlations

In [3]:
% OGARCH in MATLAB
% Listing 3.6
% Last updated August 2016
%
%

[par, Ht] = o_mvgarch(y,2, 1,1,1);
Ht = reshape(Ht,4,T)';
%% Ht comes from o_mvgarch as a 3D matrix, this transforms it into a 2D matrix
OOrho = Ht(:,3) ./ sqrt(Ht(:,1) .* Ht(:,4));

%% OOrho is a vector of correlations

In [4]:
% DCC in MATLAB
% Listing 3.8
% Last updated August 2016
%
%

[p, lik, Ht] = dcc(y,1,1,1,1);
Ht = reshape(Ht,4,T)';
DCCrho = Ht(:,3) ./ sqrt(Ht(:,1) .* Ht(:,4));

%% DCCrho is a vector of correlations

In [5]:
% Correlation comparison in MATLAB
% Listing 3.10
% Last updated June 2018
%
%

plot([EWMArho,OOrho,DCCrho])
legend('EWMA','DCC','OGARCH','Location','SouthWest')

Chapter 4: Risk Measures

  • 4.1/4.2: Expected Shortfall (ES) estimation under normality assumption
In [6]:
% ES in MATLAB
% Listing 4.2
% Last updated August 2016
%
%

p = [0.5,0.1,0.05,0.025,0.01,0.001];
VaR = -norminv(p)
ES = normpdf(norminv(p))./p
VaR =

         0    1.2816    1.6449    1.9600    2.3263    3.0902


ES =

    0.7979    1.7550    2.0627    2.3378    2.6652    3.3671


Chapter 5: Implementing Risk Forecasts

  • 5.1/5.2: Loading hypothetical stock prices, converting to returns
  • 5.3/5.4: Univariate HS Value at Risk (VaR)
  • 5.5/5.6: Multivariate HS VaR
  • 5.7/5.8: Univariate ES VaR
  • 5.9/5.10: Normal VaR
  • 5.11/5.12: Portfolio Normal VaR
  • 5.13/5.14: Student-t VaR (unavailable as of June 2018)
  • 5.15/5.16: Normal ES VaR
  • 5.17/5.18: Direct Integration Normal ES VaR
  • 5.19/5.20: MA Normal VaR
  • 5.21/5.22: EWMA VaR
  • 5.23/5.24: Two-asset EWMA VaR
  • 5.25/5.26: GARCH(1,1) VaR
In [7]:
% Download stock prices in MATLAB
% Listing 5.2
% Last updated August 2016
%
%

stocks = csvread('stocks.csv',1,0);
p1 = stocks(:,1);             % consider first two stocks
p2 = stocks(:,2); 

y1=diff(log(p1));             % convert prices to returns
y2=diff(log(p2));
y1=y1(length(y1)-4100+1:end); % length adjustment
y2=y2(length(y2)-4100+1:end); % length adjustment
y=[y1 y2];
T=length(y1);
value = 1000;                 % portfolio value
p = 0.01;                     % probability

In [8]:
% Univariate HS VaR in MATLAB
% Listing 5.4
% Last updated 2011
%
%

ys = sort(y1);   % sort returns
op = T*p;        % p percent smallest
VaR1 = -ys(op)*value
VaR1 =

   18.6586


In [9]:
% Multivariate HS VaR in MATLAB
% Listing 5.6
% Last updated 2011
%
%

w = [0.3; 0.7];    % vector of portfolio weights
yp = y*w;          % portfolio returns
yps = sort(yp);
VaR2 = -yps(op)*value
VaR2 =

   19.8145


In [10]:
% Univariate ES in MATLAB
% Listing 5.8
% Last updated 2011
%
%

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

   23.8729


In [11]:
% Normal VaR in MATLAB
% Listing 5.10
% Last updated 2011
%
%

sigma = std(y1); % estimate volatility
VaR3 = -sigma * norminv(p) * value
VaR3 =

   15.5880


In [12]:
% Portfolio normal VaR in MATLAB
% Listing 5.12
% Last updated 2011
%
%

sigma = sqrt(w' * cov(y) * w); % portfolio volatility
VaR4 = - sigma * norminv(p) *  value
VaR4 =

   18.1383


In [13]:
% Student-t VaR in MATLAB
% Listing 5.14
% Last updated 2011
%
%

scy1=y1*100;          % scale the returns
res=mle(scy1,'distribution','tlocationscale');
sigma1 = res(2)/100;  % rescale the volatility
nu = res(3);
VaR5 = - sigma1 * tinv(p,nu) * value
VaR5 =

   17.8882


In [14]:
% Normal ES in MATLAB
% Listing 5.16
% Last updated June 2018
%
%

sigma = std(y1);
ES2=sigma*normpdf(norminv(p))/p * value
ES2 =

   17.8586


In [15]:
% Direct integration ES in MATLAB
% Listing 5.18
% Last updated 2011
%
%

VaR = -norminv(p);
ES = -sigma*quad(@(q) q.*normpdf(q),-6,-VaR)/p*value
ES =

   17.8580


In [16]:
% MA normal VaR in MATLAB
% Listing 5.20
% Last updated June 2018
%
%

WE=20;
for t=T-5:T
    t1=t-WE+1;
    window=y1(t1:t);  % estimation window
    sigma=std(window);
    VaR6 = -sigma * norminv(p) * value
end
VaR6 =

   16.0505


VaR6 =

   16.1491


VaR6 =

   18.8543


VaR6 =

   18.8821


VaR6 =

   16.2305


VaR6 =

   16.1698


In [17]:
% EWMA VaR in MATLAB
% Listing 5.22
% Last updated 2011
%
%

lambda = 0.94;	
s11 = var(y1(1:30)); % initial variance
for t = 2:T	
    s11 = lambda * s11  + (1-lambda) * y1(t-1)^2;
end
VaR7 = -norminv(p) * sqrt(s11) * value
VaR7 =

   16.7534


In [18]:
% Two-asset EWMA VaR in MATLAB
% Listing 5.24
% Last updated 2011
%
%

s = cov(y);               % initial covariance
for t = 2:T
    s = lambda * s +  (1-lambda) * y(t-1,:)' * y(t-1,:);
end
sigma = sqrt(w' * s * w); % portfolio vol
VaR8 = - sigma * norminv(p) * value
VaR8 =

   20.5036


In [19]:
% GARCH in MATLAB
% Listing 5.26
% Last updated August 2016
%
%

[parameters,ll,ht]=tarch(y1,1,0,1);
omega = parameters(1)
alpha = parameters(2)
beta = parameters(3)
sigma2 = omega + alpha*y1(end)^2 + beta*ht(end) % calc sigma2 for t+1
VaR9 = -sqrt(sigma2) * norminv(p) * value
____________________________________________________________
   Diagnostic Information

Number of variables: 3

Functions 
Objective:                            tarch_likelihood
Gradient:                             finite-differencing
Hessian:                              finite-differencing (or Quasi-Newton)
 

Algorithm selected
   quasi-newton


____________________________________________________________
   End diagnostic information
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           4         -14928.8                          19.7
     1          12         -14930.8     0.00888901           13.6  
     2          16         -14931.9              1           6.78  
     3          20         -14932.6              1           5.29  
     4          24         -14932.9              1           5.26  
     5          28         -14935.1              1           13.6  
     6          32         -14935.6              1            5.4  
     7          36         -14935.9              1            2.1  
     8          40         -14935.9              1          0.568  
     9          44         -14935.9              1          0.114  
    10          48         -14935.9              1         0.0223  
    11          52         -14935.9              1        0.00062  
    12          56         -14935.9              1       0.000892  

Local minimum possible.

fminunc stopped because the size of the current step is less than
the selected value of the step size tolerance.




omega =

   9.9400e-07


alpha =

    0.0607


beta =

    0.9173


sigma2 =

   5.3679e-05


VaR9 =

   17.0442


Chapter 6: Analytical Value-at-Risk for Options and Bonds

  • 6.1/6.2: Black-Scholes function definition
  • 6.3/6.4: Black-Scholes option price calculation example
In [ ]:
% Black-Scholes function in MATLAB
% Listing 6.2
% Last updated 2011
%
%

%% To run this code block in Jupyter notebook:
%% delete all lines above the line with file bs.m, then run

%%file bs.m
function  res = bs(K,P,r,sigma,T)
	d1 = (log(P./K)+(r+(sigma^2)/2)*T)./(sigma*sqrt(T));
	d2 = d1 - sigma*sqrt(T);
	res.Call = P.*normcdf(d1,0,1)-K.*exp(-r*T).*normcdf(d2,0,1);
	res.Put = K.*exp(-r*T).*normcdf(-d2,0,1)-P.*normcdf(-d1,0,1);
	res.Delta.Call = normcdf(d1,0,1);
	res.Delta.Put = res.Delta.Call -1;
	res.Gamma = normpdf(d1,0,1)./(P*sigma*sqrt(T));
end
In [20]:
% Black-Scholes in MATLAB
% Listing 6.4
% Last updated 2011
%
%

f=bs(90,100,0.05,0.2,0.5)
f = 

  struct with fields:

     Call: 13.4985
      Put: 1.2764
    Delta: [1x1 struct]
    Gamma: 0.0172


Chapter 7: Simulation Methods for VaR for Options and Bonds

  • 7.1/7.2: Plotting normal distribution transformation
  • 7.3/7.4: Random number generation from Uniform(0,1), Normal(0,1)
  • 7.5/7.6: Bond pricing using yield curve
  • 7.7/7.8: Yield curve simulations
  • 7.9/7.10: Bond price simulations
  • 7.11/7.12: Black-Scholes analytical pricing of call
  • 7.13/7.14: Black-Scholes Monte Carlo simulation pricing of call
  • 7.15/7.16: Option density plots
  • 7.17/7.18: VaR simulation of portfolio with only underlying
  • 7.19/7.20: VaR simulation of portfolio with only call
  • 7.21/7.22: VaR simulation of portfolio with call, put and underlying
  • 7.23/7.24: Simulated two-asset returns
  • 7.25/7.26: Two-asset portfolio VaR
  • 7.27/7.28: Two-asset portfolio VaR with a call
In [21]:
% Transformation in MATLAB
% Listing 7.2
% Last updated 2011
%
%

x=-3:0.1:3;
plot(x,normcdf(x))

In [22]:
% Various RNs in MATLAB
% Listing 7.4
% Last updated August 2016
%
%

rng default; % set seed

S=10;

rand(S,1)
randn(S,1)
trnd(4,S,1)
ans =

    0.8147
    0.9058
    0.1270
    0.9134
    0.6324
    0.0975
    0.2785
    0.5469
    0.9575
    0.9649


ans =

   -1.3499
    3.0349
    0.7254
   -0.0631
    0.7147
   -0.2050
   -0.1241
    1.4897
    1.4090
    1.4172


ans =

    1.2458
   -1.8803
    0.4892
    2.4722
    0.5898
    1.0088
    0.8056
   -0.2653
    0.2206
   -0.8372


In [23]:
% Price bond in MATLAB
% Listing 7.6
% Last updated 2011
%
%

yield = [5.00 5.69 6.09 6.38 6.61...
         6.79 6.94 7.07 7.19 7.30]; % yield curve
T = length(yield); 
r=0.07;                             % initial yield rate
Par=10;                             % par value
coupon=r*Par;                       % coupon payments
cc=zeros(1,T)+coupon;               % vector of cash flows 
cc(10)=cc(10)+Par;                  % add par to cash flows
P=sum(cc./((1+yield./100).^(1:T)))  % calculate price
P =

    9.9132


In [24]:
% Simulate yields in MATLAB
% Listing 7.8
% Last updated 2011
%
%

randn('state',123);         % set the seed
sigma = 1.5;                % daily yield volatility
S = 8;                      % number of simulations
r = randn(1,S)*sigma;       % generate random numbers

ysim=nan(T,S);
for i=1:S
    ysim(:,i)=yield+r(i);
end
ysim=repmat(yield',1,S)+repmat(r,T,1);
plot(ysim)

In [25]:
% Simulate bond prices in MATLAB
% Listing 7.10
% Last updated 2011
%
%

SP = nan(S,1);
for s = 1:S                                        % S simulations
    SP(s) = sum(cc./(1+ysim(:,s)'./100).^((1:T)));
end
SP = SP-(mean(SP) - P);                            % correct for mean
bar(SP)

In [26]:
% Black-Scholes valuation in MATLAB
% Listing 7.12
% Last updated 2011
%
%

P0 = 50;                      % initial spot price
sigma = 0.2;                  % annual volatility
r = 0.05;                     % annual interest
T = 0.5;                      % time to expiration
X = 40;                       % strike price
f = bs(X,P0,r,sigma,T)        % analytical call price

%% this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
f = 

  struct with fields:

     Call: 11.0873
      Put: 0.0997
    Delta: [1x1 struct]
    Gamma: 0.0107


In [27]:
% Black-Scholes simulation in MATLAB
% Listing 7.14
% Last updated 2011
%
%

randn('state',0);         % set seed
S = 1e6;                  % number of simulations
F = P0*exp(r*T);          % futures price
ysim=randn(S,1)*sigma*sqrt(T)-0.5*T*sigma^2; % sim returns, lognorm corrected
F=F*exp(ysim);            % sim futures price
SP = F-X;                 % payoff
SP(find(SP < 0)) = 0;     % set negative outcomes to zero
fsim =SP*exp(-r*T) ;      % discount
mean(fsim)                % simulated price
ans =

   11.0996


In [28]:
% Option density plots in MATLAB
% Listing 7.16
% Last updated 2011
%
%

subplot(1,2,1)
histfit(F);

subplot(1,2,2)
hist(fsim,100);

In [29]:
% Simulate VaR in MATLAB
% Listing 7.18
% Last updated 2011
%
%

randn('state',0);   % 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 = randn(S,1)*sqrt(s2)+r/365-0.5*s2; % sim returns
Psim = P*exp(ysim);  % sim future prices 
q = sort(Psim-P);  % simulated P/L
VaR1 = -q(S*p)
VaR1 =

    2.2906


In [30]:
% Simulate option VaR in MATLAB
% Listing 7.20
% Last updated 2011
%
%

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

    1.2150


In [31]:
% Example 7.3 in MATLAB
% Listing 7.22
% Last updated 2011
%
%

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

    1.4951


In [32]:
% Simulated two-asset returns in MATLAB
% Listing 7.24
% Last updated 2011
%
%

randn('state',12)                 % set seed
mu = [r/365 r/365]';              % return mean
Sigma=[0.01 0.0005; 0.0005 0.02]; % covariance matrix
y = mvnrnd(mu,Sigma,S);           % simulated returns

In [33]:
% Two-asset VaR in MATLAB
% Listing 7.26
% Last updated 2011
%
%

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

   25.9700


In [34]:
% A two-asset case in MATLAB with an option
% Listing 7.28
% Last updated 2011
%
%

f = bs(P(2),P(2),r,sigma,T);
fsim=bs(P(2),Psim(:,2),r,sigma,T-(1/365));
q = sort(fsim.Call+Psim(:,1)-f.Call-P(1)); 
VaR5 = -q(p*S)
VaR5 =

   20.8125


Chapter 8: Backtesting and Stress Testing

  • 8.1/8.2: Loading hypothetical stock prices, converting to returns
  • 8.3/8.4: Setting up backtest
  • 8.5/8.6: Running backtest for EWMA/MA/HS/GARCH VaR
  • 8.7/8.8: Backtesting analysis for EWMA/MA/HS/GARCH VaR
  • 8.9/8.10: Bernoulli coverage test
  • 8.11/8.12: Independence test
  • 8.13/8.14: Running Bernoulli/Independence test on backtests
  • 8.15/8.16: Running backtest for EWMA/HS ES
  • 8.17/8.18: Backtesting analysis for EWMA/HS ES
In [35]:
% Load data in MATLAB
% Listing 8.2
% Last updated August 2016
%
%

price = csvread('index.csv', 1, 0);

y=diff(log(price)); % get returns

In [36]:
% Set backtest up in MATLAB
% Listing 8.4
% Last updated August 2016
%
%

T = length(y);  % number of obs for return y
WE = 1000;      % estimation window length 
p = 0.01;       % probability
l1 = WE*p ;     % HS observation
value = 1;      % portfolio value
VaR = NaN(T,4); % matrix for forecasts
%% EWMA setup
lambda = 0.94;
s11 = var(y(1:30));
for t = 2:WE
    s11=lambda*s11+(1-lambda)*y(t-1)^2;
end

In [39]:
% Running backtest in MATLAB
% Listing 8.6
% Last updated June 2018
%
%

for t = WE+1:T
    t1 = t-WE;          % start of the data window
    t2 = t-1;           % end of data window
    window = y(t1:t2) ; % data for estimation
    s11 = lambda*s11  + (1-lambda)*y(t-1)^2;
    VaR(t,1) = -norminv(p) * sqrt(s11)  *value; % EWMA
    VaR(t,2) = -std(window)*norminv(p)*value; % MA
    ys = sort(window);
    VaR(t,3) = -ys(l1)*value; % HS
    [par,ll,ht]=tarch(window,1,0,1);
    h=par(1)+par(2)*window(end)^2+par(3)*ht(end);
    VaR(t,4) = -norminv(p)*sqrt(h)*value; % GARCH(1,1)
end

In [40]:
% Backtesting analysis in MATLAB
% Listing 8.8
% Last updated 2011
%
%

for i=1:4
    VR = length(find(y(WE+1:T)<-VaR(WE+1:T,i)))/(p*(T-WE)); 
    s = std(VaR(WE+1:T,i));         
    disp([i VR s])          
end
plot([y(WE+1:T) VaR(WE+1:T,:)])
    1.0000    2.0744    0.0121

    2.0000    1.8606    0.0066

    3.0000    1.2404    0.0123

    4.0000    1.6895    0.0115


In [ ]:
% Bernoulli coverage test in MATLAB
% Listing 8.10
% Last updated 2011
%
%

function res=bern_test(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)
	res=-2*(al-bl);
end
In [ ]:
% Independence test in MATLAB
% Listing 8.12
% Last updated 2011
%
%

function res=ind_test(V)
	T=length(V);
	J=zeros(T,4);
	for i = 2:T
		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;
	end	
	V_00=sum(J(:,1));
	V_01=sum(J(:,2));
	V_10=sum(J(:,3));
	V_11=sum(J(:,4));
	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;
	res= -2*(al-bl);
end
In [41]:
% Backtesting S&P 500 in MATLAB
% Listing 8.14
% Last updated June 2018
%
%

ya=y(WE+1:T);
VaRa=VaR(WE+1:T,:);
for i=1:4
	q=find(y(WE+1:T)<-VaR(WE+1:T,i));
	v=VaRa*0;
	v(q,i)=1;
	ber=bern_test(p,v(:,i));
	in=ind_test(v(:,i));
	disp([i,ber,1-chi2cdf(ber,1),in,1-chi2cdf(in,1)])
end
bl =

 -471.9091

    1.0000   41.6257    0.0000    0.4411    0.5066


bl =

 -432.8188

    2.0000   27.9039    0.0000   17.4418    0.0000


bl =

 -312.2446

    3.0000    2.5354    0.1113   11.4848    0.0007


bl =

 -400.7082

    4.0000   18.6034    0.0000    0.2984    0.5849


In [42]:
% Backtest ES in MATLAB
% Listing 8.16
% Last updated 2011
%
%

VaR = NaN(T,2);  % VaR forecasts for 2 models 
ES = NaN(T,2);   % ES forecasts for 2 models 
for t = WE+1:T
    t1 = t-WE; 
    t2 = t-1; 
    window = y(t1:t2) ;
    
    s11 = lambda * s11  + (1-lambda) * y(t-1)^2; 
    VaR(t,1) = -norminv(p) * sqrt(s11)  *value; % EWMA
    ES(t,1) = sqrt(s11) * normpdf(norminv(p)) / p;

    ys = sort(window);
    VaR(t,2) = -ys(l1) * value;          % HS
    ES(t,2) = -mean(ys(1:l1)) * value;  
end

In [43]:
% ES in MATLAB
% Listing 8.18
% Last updated 2011
%
%

VaRa = VaR(WE+1:T,:);
ESa = ES(WE+1:T,:);
for i = 1:2
	q = find(ya <= -VaRa(:,i));
	nES = mean(ya(q) ./ -ESa(q,i));
	disp([i,nES])
end
    1.0000    1.2235

    2.0000    1.0537


Chapter 9: Extreme Value Theory

  • 9.1/9.2: Calculation of tail index from returns
In [44]:
% Hill estimator in MATLAB
% Listing 9.2
% 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

% END
iota =

    2.6297



All rights reserved, Jon Danielsson, 2011-2018.