Chapter 5. Implementing Risk Forecasts

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

import numpy as np
from scipy import stats
p = np.loadtxt('stocks.csv',delimiter=',',skiprows=1)
p = p[:,[0,1]]      # consider two stocks
y1 = np.diff(np.log(p[:,0]), n=1, axis=0)
y2 = np.diff(np.log(p[:,1]), n=1, axis=0)
y = np.stack([y1,y2], axis = 1)
T = len(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)

from math import ceil
ys = np.sort(y1) # sort returns
op = ceil(T*p)    # p percent smallest
VaR1 = -ys[op - 1] * value
print(VaR1)

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 = np.squeeze(np.matmul(y, w)) # portfolio returns
yps = np.sort(yp)
VaR2= -yps[op - 1] * value
print(VaR2)

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

ES1 = -np.mean(ys[:op]) * value
print(ES1)

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

sigma = np.std(y1, ddof=1) # estimate volatility
VaR3 = -sigma * stats.norm.ppf(p) * value 
print(VaR3)

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

sigma = np.sqrt(np.mat(w)*np.mat(np.cov(y,rowvar=False))*np.transpose(np.mat(w)))[0,0]
VaR4 = -sigma * stats.norm.ppf(p) * value 
print(VaR4)

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 = stats.t.fit(scy1)
sigma = res[2]/100      # rescale volatility
nu = res[0]
VaR5 = -sigma*stats.t.ppf(p,nu)*value
print(VaR5)

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

sigma = np.std(y1, ddof=1)
ES2 = sigma * stats.norm.pdf(stats.norm.ppf(p)) / p * value
print(ES2)

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

from scipy.integrate import quad
VaR = -stats.norm.ppf(p)
integrand = lambda q: q * stats.norm.pdf(q)
ES = -sigma * quad(integrand, -np.inf, -VaR)[0] / p * value
print(ES)

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 in range(T-5,T+1):
    t1 = t-WE
    window = y1[t1:t]      # estimation window
    sigma = np.std(window, ddof=1)
    VaR6 = -sigma*stats.norm.ppf(p)*value
    print (VaR6)

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)

lmbda = 0.94
s11 = np.var(y1[0:30], ddof = 1)     # initial variance
for t in range(1, T):
    s11 = lmbda*s11 + (1-lmbda)*y1[t-1]**2
VaR7 = -np.sqrt(s11)*stats.norm.ppf(p)*value
print(VaR7)

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 = np.cov(y, rowvar = False)
for t in range(1,T):
    s = lmbda*s+(1-lmbda)*np.transpose(np.asmatrix(y[t-1,:]))*np.asmatrix(y[t-1,:])
sigma = np.sqrt(np.mat(w)*s*np.transpose(np.mat(w)))[0,0]
VaR8 = -sigma * stats.norm.ppf(p) * value
print(VaR8)

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)

from arch import arch_model 
am = arch_model(y1, mean = 'Zero', vol='Garch', p=1, o=0, q=1, dist='Normal', rescale = False)
res = am.fit(update_freq=5, disp = "off")
omega = res.params.loc['omega']
alpha = res.params.loc['alpha[1]']
beta = res.params.loc['beta[1]']
sigma2 = omega + alpha*y1[T-1]**2 + beta * res.conditional_volatility[-1]**2
VaR9 = -np.sqrt(sigma2) * stats.norm.ppf(p) * value
print(VaR9)