Chapter 5. Implementing Risk Forecasts

p = read.csv('stocks.csv')
y = apply(log(p),2,diff)
portfolio_value = 1000 
p = 0.01

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

y1 = y[,1]                     # select one asset
ys = sort(y1)                  # sort returns
op = ceiling(length(y1)*p)     # p percent smallest, rounded up
VaR1 = -ys[op]*portfolio_value
print(VaR1)

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

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

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

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

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

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

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

sigma = sqrt(t(w) %*% cov(y) %*% w)[1] # portfolio volatility
VaR4 = -sigma * qnorm(p)*portfolio_value
print(VaR4)

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

library(QRM)
scy1 = (y1)*100                     # scale the returns 
res = fit.st(scy1) 
sigma1 = unname(res$par.ests['sigma']/100)  # rescale the volatility
nu = unname(res$par.ests['nu'])
VaR5 = - sigma1 * qt(df=nu,p=p) *  portfolio_value
print(VaR5)

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 = sd(y1)
ES2 = sigma*dnorm(qnorm(p))/p * portfolio_value
print(ES2)

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

VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$value/p*portfolio_value
print(ES)

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

WE=20
for (t in seq(length(y1)-5,length(y1))){
  t1=t-WE+1
  window= y1[t1:t] # estimation window
  sigma=sd(window)
  VaR6 = -sigma * qnorm(p) * portfolio_value
  print(VaR6)
}

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

lambda = 0.94
s11 = var(y1)             # initial variance, using unconditional
for (t in 2:length(y1)){ 
  s11 = lambda * s11  + (1-lambda) * y1[t-1]^2
}
VaR7 = -qnorm(p) * sqrt(s11) * portfolio_value
print(VaR7)

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 

s = cov(y)                        # initial covariance
for (t in 2:dim(y)[1]){
  s = lambda*s + (1-lambda)*y[t-1,] %*% t(y[t-1,])
}
sigma = sqrt(t(w) %*% s %*% w)[1] # portfolio vol
VaR8 = -sigma * qnorm(p) * portfolio_value
print(VaR8)

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 

library(rugarch)
spec = ugarchspec(variance.model = list( garchOrder = c(1, 1)),
                  mean.model = list( armaOrder = c(0,0),include.mean = FALSE))
res = ugarchfit(spec = spec, data = y1, solver = "hybrid")
omega = res@fit$coef['omega']
alpha = res@fit$coef['alpha1']
beta = res@fit$coef['beta1']
sigma2 = omega + alpha * tail(y1,1)^2 + beta * tail(res@fit$var,1)  
VaR9 = -sqrt(sigma2) * qnorm(p) * portfolio_value 
names(VaR9)="VaR"
print(VaR9)

[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