Chapter 5. Implementing Risk Forecasts

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));
y=[y1 y2];
T=length(y1);
value = 1000;                 % portfolio value
p = 0.01;                     % probability

using CSV, DataFrames;
p = CSV.read("stocks.csv", DataFrame);
y1 = diff(log.(p[:,1]));
y2 = diff(log.(p[:,2]));
y = hcat(y1,y2);
T = size(y,1); 
value = 1000; # portfolio value 
p = 0.01;     # probability

ys = sort(y1);   % sort returns
op = ceil(T*p);  % p percent smallest, rounded up to meet VaR probability requirement
VaR1 = -ys(op)*value

ys = sort(y1)            # sort returns
op = ceil(Int, T*p)      # p percent smallest, rounding up
VaR1 = -ys[op] * value
println("Univariate HS VaR ", Int(p*100), "%: ", round(VaR1, digits = 3), " USD")

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

w = [0.3; 0.7]          # vector of portfolio weights
yp = y * w              # portfolio returns
yps = sort(yp)
VaR2 = -yps[op] * value
println("Multivariate HS VaR ", Int(p*100), "%: ", round(VaR2, digits = 3), " USD")

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

using Statistics;
ES1 = -mean(ys[1:op]) * value
println("ES: ", round(ES1, digits = 3), " USD")

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

using Distributions;
sigma = std(y1);      # estimate volatility
VaR3 = -sigma * quantile(Normal(0,1),p) * value
println("Normal VaR", Int(p*100), "%: ", round(VaR3, digits = 3), " USD")

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

sigma = sqrt(w'*cov(y)*w)   # portfolio volatility
VaR4 = -sigma * quantile(Normal(0,1), p) * value
println("Portfolio normal VaR", Int(p*100), "%: ", round(VaR4, digits = 3), " USD")

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

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

sigma = std(y1)
ES2 = sigma * pdf(Normal(0,1), (quantile(Normal(0,1), p))) / p * value
println("Normal ES: ", round(ES2, digits = 3), " USD")

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

using QuadGK;
VaR = -quantile(Normal(0,1), p)
integrand(x) = x * pdf(Normal(0,1), x)
ES3 = -sigma * quadgk(integrand, -Inf, -VaR)[1] / p * value
println("Normal integrated ES: ", round(ES3, digits = 3), " USD")

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

WE = 20
for t in T-5:T
    t1 = t-WE
    window = y1[t1+1:t] # estimation window
    sigma = std(window)
    VaR6 = -sigma*quantile(Normal(0,1),p)*value
    println("MA Normal VaR", Int(p*100), "% using observations ", t1, " to ", t, ": ",
        round(VaR6, digits = 3), " USD")
end

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 

lambda = 0.94
s11 = var(y1) # initial variance
for t in 2:T
    s11 = lambda * s11 + (1-lambda) * y1[t-1]^2 
end
VaR7 = -sqrt(s11) * quantile(Normal(0,1), p) * value
println("EWMA VaR ", Int(p*100), "%: ", round(VaR7, digits = 3), " USD")

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 

s = cov(y) # initial covariance
for t in 2:T
    s = lambda * s + (1-lambda) * y[t-1,:] * (y[t-1,:])'
end
sigma = sqrt(w'*s*w) # portfolio vol
VaR8 = -sigma * quantile(Normal(0,1), p) * value
println("Two-asset EWMA VaR ", Int(p*100), "%: ", round(VaR8, digits = 3), " USD")

[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 

using ARCHModels;
garch1_1 = fit(GARCH{1,1}, y1; meanspec = NoIntercept);
garch_VaR_in = VaRs(garch1_1, :0.01)
cond_vol = predict(garch1_1, :volatility)     # 1-day-ahead conditional volatility
garch_VaR_out = -cond_vol * quantile(garch1_1.dist, p) * value
println("GARCH VaR ", Int(p*100), "%: ", round(garch_VaR_out, digits = 3), " USD")