i = 9000

source("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Final_function.R")

load("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500/outSP500.RDATA") # market proxy
#attributes(out.p)

#q=.50: method "stHS"
#q=.95,.99: method "stEVT"

## to extract vector of innovations (Z^I) based on
## AR(1)-GARCH(1,1) process with skew-t innovations

## Residuals from the AR(1)-GARCH(1,1) filter with skew-t innovations
zmatI <- out.p$zmat.st # n x 500 matrix, each row gives 500 residuals
# for estimation within one sliding window
## for CoVaR forecasts of L^I|L^j>VaRq based on CoVaR for Z's
muI <- out.p$mut.st # \mu_{t+1}^I: forecasts of conditional mean
sigI <- out.p$sigt.st # \sig_{t+1}^I: forecasts of volatility

n=dim(zmatI)[1] # number of estimation windows

## ===================================================
load("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500/outBA.RDATA") # ticker 'BA'

## Residuals from the AR(1)-GARCH(1,1) filter with skew-t innovations
zmatBA <- out.p$zmat.st # n x 500 matrix, each row gives 500 residuals
# for estimation within one sliding window

## VaR forecasts: using  AR(1)-GARCH(1,1) process with skew-t innovations and EVT
VaR.BA <- out.p$stEVT# n x 3 matrix; col1: q=.50, col2: q=.95, col3: q=.99
VaR.BA95 <- out.p$stEVT[,2] # .95-VaR forecasts for losses
VaR.BA99 <- out.p$stEVT[,3] # .99-VaR forecasts for losses

## \mu_{t+1}^j and \sig_{t+1}^j
mu.BA <- out.p$mut.st
sig.BA <- out.p$sigt.st

## VaRq forecasts (estimates) for Z^j
VarZ.BA95 <- (VaR.BA95-mu.BA)/sig.BA # q=.95
VarZ.BA99 <- (VaR.BA99-mu.BA)/sig.BA # q=.99```

i = 9000
zI=zmatI[i,]; zj=zmatBA[i,]; varZj=c(VarZ.BA95[i], VarZ.BA99[i]);q=c(.95,.99);
sim_y = as.matrix(cbind(zj, zI))
plot(sim_y)

fit = New_cov_update(sim_y, prob = 0.85) # our EVT method

tmp95 = NA; tmp99 = NA; tmpE95 = NA; tmpE99 = NA; tmpS95 = NA; tmpS99 =NA;

##### fitted EVT parameters #####
#alpha1, alpha2, rho, nu, sd1, sd2)
print(fit$EVT)

### CoVaR calculation
CoVaR = function(x, q, fit, Var_j){
  diag_inv_EVT = matrix(c(1/fit$EVT[5], 0,0,1/fit$EVT[6]),2,2)
  tmp_value_EVT = diag_inv_EVT %*% c(Var_j, x)
  z_EVT =  tmp_value_EVT[2] / tmp_value_EVT[1] - fit$EVT[3]
  theorm_EVT = 1-Thrm4_1(z_EVT, fit$EVT[1], fit$EVT[2], fit$EVT[3], fit$EVT[4])
  return(theorm_EVT - q)
}

try(tmp95 <- as.numeric(uniroot(CoVaR, q = 1-q[1], fit = fit, Var_j = varZj[1], lower = 0.1, upper = 100, tol = 8e-6)$root, print =1))
if (!is.na(tmp95)){
  try(tmp99 <- as.numeric(uniroot(CoVaR, q = 1-q[2], fit = fit, Var_j = varZj[2], lower = 0.1, upper = 100, tol = 8e-6)$root, print =1))
}

tmpE95 <- quantile(sim_y[sim_y[,1]>varZj[1],][,2],prob=q[1])
tmpE99 <- NaN
try(tmpE99 <- quantile(sim_y[sim_y[,1]>varZj[2],][,2],prob=q[2],na.rm = TRUE))

fit_skewR = mst.fit(y = -sim_y, plot.it=F) # Return fit using bivariate skew-t distribution 
#### calculate marginal density from bivarite-t fitting
xi = fit_skewR$dp$beta
Omega = fit_skewR$dp$Omega
alpha = fit_skewR$dp$alpha
nu = fit_skewR$dp$df
Omega_11 = Omega[1,1] - Omega[1,2] / Omega[2,2] * Omega[2,1]
alpha1_bar = (alpha[1] + Omega[1,2]/Omega[1,1] * alpha[2]) / (1 + alpha[2]^2*Omega_11)
##  \hat{q} the estiamted  marginal probability for 95%
sn::pst(x = -varZj[1], xi = xi[1], omega = Omega_11, alpha = alpha1_bar, nu = nu) 

##  \hat{q} the estiamted  marginal probability for 99%
sn::pst(x = -varZj[2], xi = xi[1], omega = Omega_11, alpha = alpha1_bar, nu = nu)

### CoVaR calculation
CoVaR_Skew = function(x, q, fit_skewR, Var_j){
  xi = fit_skewR$dp$beta
  Omega = fit_skewR$dp$Omega
  alpha = fit_skewR$dp$alpha
  nu = fit_skewR$dp$df
  joint = sn::pmst(x = c(-Var_j, -x), xi = xi, Omega = Omega, alpha = alpha, nu = nu)[1]
  Omega_11 = Omega[1,1] - Omega[1,2] / Omega[2,2] * Omega[2,1]
  alpha1_bar = (alpha[1] + Omega[1,2]/Omega[1,1] *alpha[2]) / (1 + alpha[2]^2*Omega_11)
  marginal = sn::pst(x = -Var_j,xi = xi[1], omega = Omega_11, alpha = alpha1_bar, nu = nu)
  return(joint - q * marginal)
}

try(tmpS95 <- as.numeric(uniroot(CoVaR_Skew, q = 1-q[1], fit_skewR = fit_skewR,  Var_j = varZj[1], lower = 0.01, upper = 100, tol = 8e-6, check.conv = T, trace = 1)$root))
try(tmpS99 <- as.numeric(uniroot(CoVaR_Skew, q = 1-q[2], fit_skewR = fit_skewR,  Var_j = varZj[2], lower = 0.01, upper = 100, tol =  8e-6, check.conv = T, trace = 1)$root))

#### re-fit marginal density 
fit_st = st.mple(y = -zj)
##  \hat{q} the estiamted  marginal probability for 95%
sn::pst(x = -varZj[1], xi = fit_st$dp[1], omega = fit_st$dp[2], alpha = fit_st$dp[3], nu = fit_st$dp[4])

##  \hat{q} the estiamted  marginal probability for 99%
sn::pst(x = -varZj[2], xi = fit_st$dp[1], omega = fit_st$dp[2], alpha = fit_st$dp[3], nu = fit_st$dp[4])

CoVaR_Skew = function(x, q, fit_skewR, Var_j, fit_st){
  joint = sn::pmst(x = c(-Var_j, -x), xi = xi, Omega = Omega, alpha = alpha, nu = nu)[1]
  marginal = sn::pst(x = -Var_j, xi = fit_st$dp[1], omega = fit_st$dp[2], alpha = fit_st$dp[3], nu = fit_st$dp[4])
  return(joint - q * marginal)
}
try(tmpS95_re <- as.numeric(uniroot(CoVaR_Skew, q = 1-q[1], fit_skewR = fit_skewR,  Var_j = varZj[1], fit_st=fit_st, lower = 0.01, upper = 100, tol = 8e-6, check.conv = T, trace = 1)$root))
try(tmpS99_re <- as.numeric(uniroot(CoVaR_Skew, q = 1-q[2], fit_skewR = fit_skewR,  Var_j = varZj[2], fit_st = fit_st, lower = 0.01, upper = 100, tol =  8e-6, check.conv = T, trace = 1)$root))

### output
tmp = c(tmp95, tmp99, tmpE95, tmpE99, tmpS95, tmpS99, tmpS95_re, tmpS99_re)
### EVT95 EVT99 Emp95 Emp99 FP95 FP99 FP95_new FP99_new
print(tmp)

condition = zj > VarZ.BA95[i]
condition99 = zj > VarZ.BA99[i]

#### In-sample test for 95% 
c(sum(zI[condition] > tmp[1])/sum(condition)*100, sum(zI[condition] > tmp[3])/sum(condition)*100
, sum(zI[condition] > tmp[5])/sum(condition)*100, sum(zI[condition] > tmp[7])/sum(condition)*100)

#### In-sample test for 99% 
c(sum(zI[condition99] > tmp[2])/sum(condition99)*100,sum(zI[condition99] >tmp[4])/sum(condition99)*100,
  sum(zI[condition99] > tmp[6])/sum(condition99)*100, sum(zI[condition99] > tmp[8])/sum(condition99)*100)