###########################################
# Topic 7                                 # 
# Copulas                                 # 
###########################################
rm(list = ls())
require(zoo)
require(xts)
require(knitr)
require(ggplot2)
require(copula)
require(psych)
require(MASS)
source("Block2Functions.R")


# 1. Loading data for WIG20 stocks
####################################


# setwd(.)

filename <- "wig20.Rdata"
load(filename)

# filename  <- "http://web.sgh.waw.pl/~mrubas/EFII/Dane/wig.csv"
# data      <- read.csv(filename, row.names = "X")
# data      <- zoo(data,as.Date(rownames(data)))

dates     <- index(data) 
names(data)

# selecting the sample
startDate <- as.Date("2005-01-01")
endDate   <- as.Date("2050-01-01")
y         <- window(data, start=startDate, end=endDate)

# selecting a pair of stocks
y          <- na.omit(y[,c("pko","kgh")])

# log returns of stocks
dy  <- 100*diff(log(y))

# porfolio weights
w   <- c(0.5,0.5)

# 2. Multivariate normal distribution   
######################################
p     <- 0.05
H     <- 10

mu    <- colMeans(dy); m   <- mu%*%w;
Sig   <- cov(dy);      s   <- sqrt(t(w)%*%Sig%*%w)

VaRmnorm     <- qnorm(p)*s + m             
ESmnorm      <- m - s*dnorm(qnorm(p))/p    

VaRHmnorm    <- sqrt(1:H)*qnorm(p)*s + (1:H)*m    
ESHmnorm     <- (1:H)*m - sqrt(1:H)*s*dnorm(qnorm(p))/p    

load("VaRESH.RData")

kable(data.frame(h=1:H, unorm = -VaRH[,"VaRHnorm"],mnorm = -VaRHmnorm),digits=3)
kable(data.frame(h=1:H, unorm = -ESH[,"ESHnorm"],mnorm = -ESHmnorm),digits=3)

# 3. Is correlation a good description of dependence
#############################################################################

# Simulated data from normal distribution 
dyN    <- mvrnorm(n = dim(dy)[1], mu=mu, Sigma = Sig)

cov(dy)
cov(dyN)

# Comparison: data vs simulated data from normal distribution
par(mfrow=c(1,2), cex = 0.75, bty="l")
plot(coredata(dy), xlab="asset 1", ylab="asset 2", main = "Data")
plot(dyN, xlab="asset 1", ylab="asset 2", main = "Simulated data from Normal dist.")

# 4. Marginal and joint distribtions             
####################################
require(psych)

rho0 <- 0.5
Sig0 <- matrix(c(1, rho0, rho0, 1),2,2)    
mu0  <- c(0,0)

z0  <- mvrnorm(n = 1000, mu=mu0, Sigma = Sig0)
pairs.panels(z0)

# marginal distributions (should be from Uniform[0,1])       
u0 <- pnorm(z0)                                                
pairs.panels(u0)

# arbitrary changing marginals, i.e gamma / t-Student marginal distributions + normal Copula
x0a <- qchisq(u0[,1],df=5)      
x0b <- qt(u0[,2],df=3)
x0  <- cbind(x0a,x0b)
pairs.panels(x0)


# 5. Introduction: create your own series with copulas and marginals                    
######################################################################

require(copula)

# STEP 1. defining copula
Cop1 <- ellipCopula(family = "normal", dim = 2, param = 0.5)
Cop2 <- ellipCopula(family = "t", dim = 2, param = 0.5, df = 3)
Cop3 <- archmCopula(family = "clayton", param = 3)
Cop4 <- archmCopula(family = "frank", param = 5)
Cop5 <- archmCopula(family = "gumbel", param = 2)
Cop   <- Cop1

# random draws from copula, i.e. for (U,V)
tempUV <- rCopula(1000,Cop)
colnames(tempUV) <- c("U", "V")

# windows()
par(mfrow=c(2,2), cex = 0.7, bty="l", mar=c(3,3,3,3))
persp(Cop, pCopula, main = "copula: C(U,V)")
persp(Cop, dCopula, main = "copula density: c(U,V)")
contour(Cop, dCopula, main = "contour for c(U,V)")
plot(tempUV,pch=19, main = "Scatter plot: U vs. V")


# STEP 2. adding marginal distributions
MVD1  <- mvdc(copula=Cop, margins=c("unif","unif"), paramMargins=list(list(min=0, max=1), list(min=0, max=1)) ) 
MVD2  <- mvdc(copula=Cop, margins=c("norm","norm"), paramMargins=list(list(mean=0, sd=1), list(mean=0, sd=1)) )  
MVD3  <- mvdc(copula=Cop, margins=c("t","t"), paramMargins=list(list(df=3), list(df=3)))  
#MVD4  <- mvdc(copula=Cop, margins=c("gamma","gamma"), paramMargins=list(list(shape=2, scale=1), list(shape=2, scale=1)) )  
#MVD5  <- mvdc(copula=Cop, margins=c("gamma","t"), paramMargins=list(list(shape=2, scale=1), list(df=3)) )  
MVD   <- MVD2

# random draws from multivariate distribution, i.e. for (X,Y)
tempXY  <- rMvdc(1000,MVD)                                                                                   # generating draws from copula
colnames(tempXY) <- c("X", "Y")

xy_lim = c(-3,3)
#windows()
par(mfrow=c(2,2), cex = 0.7, bty="l",mar=c(3,3,3,3))
persp(MVD, pMvdc, xlim=xy_lim, ylim=xy_lim, main = "cumulative density: F(X,Y)")
persp(MVD, dMvdc, xlim=xy_lim, ylim=xy_lim, main = "density: f(X,Y)")
contour(MVD, dMvdc, xlim = xy_lim, ylim=xy_lim, main = "Contour for f(X,Y)")
plot(tempXY,pch=19, main = "Scatter plot: X vs Y",xlim = xy_lim, ylim=xy_lim)


# 6. Estimate parameters of normal copula with normal marginals 
#################################################################

# Step 1: Parameters of normal distribution for marginals
m1  <- mean(dy[,1]); s1 <- sd(dy[,1]); U  <- pnorm((dy[,1]-m1)/s1); 
m2  <- mean(dy[,2]); s2 <- sd(dy[,2]); V  <- pnorm((dy[,2]-m2)/s2)
UV  <- as.matrix(cbind(U,V))

par(mfrow=c(1,2), cex = 0.7, bty="l")
plot(coredata(dy),main = "Scatter plot for data")
plot(UV, main = "Scatter plot for marginal distr.")

# Step 2: estimating the parameter of copula 
Cop <- normalCopula()

sVal   <- fitCopula(Cop, UV, method = "itau")@estimate   # method of moments for the starting value
CopEst <- fitCopula(Cop, UV, method = "ml", start=sVal)  # ML estimation

CopEst
cor(dy)

# Step 3: Simulating the data for (X,Y) to calculate VaR
H    = 10
Nsim = 10000
p    = 0.05     # tolerance level

# a. Copula
UVsim <- rCopula(Nsim*H,CopEst@copula)
XYsim <- cbind(qnorm(UVsim[,1])*s1+m1, qnorm(UVsim[,2])*s2+m2)

par(mfrow=c(1,2), cex = 0.7, bty="l")
plot(coredata(dy),  main='Returns')
plot(XYsim[1:dim(dy)[1],],col='red', main='Simulated')

# portfolio return
Rdraws     <- matrix(XYsim%*%w,H,Nsim)
temp       <- RdrawsToVaRES(RDraws,p)
VaRHcopN   <- temp$VaR
ESHcopN    <- temp$ES

# Analytical from multivariate normal vs simulated from norm marginals +norm copula
kable(data.frame(h=1:H, copNorm = -VaRHcopN,mNorm = -VaRHmnorm),digits=3)
kable(data.frame(h=1:H, copNorm = -ESHcopN, mNorm = -ESHmnorm),digits=3)


# 6a. Normal Copula: 1 step procedure 
#######################################
rho  <- cor(dy)[2,1]
Thet1<- c(m1,s1,m2,s2,rho) # starting values

Cop1    <- normalCopula(rho)
MVD1    <- mvdc(Cop1, c("norm","norm"),list(list(mean = 1, sd =1), list(mean = m2, sd =s2))) 
MVDEst1 <- fitMvdc(coredata(dy), MVD1, start = Thet1)
XYSim   <- rMvdc(Nsim*H,MVDEst1@mvdc) 

# 7. Selecting the best copula
#####################################

# Step 1: selecting marginals
require(rugarch)
require(moments)

mType <- 2

# Empirical marginal cdf 
if(mType==1){
  T  <- dim(dy)[1]
  UV <- matrix(0,T,2)
  UV[order(coredata(dy)[,1]),1] <- seq(1,T)/T;  
  UV[order(coredata(dy)[,2]),2] <- seq(1,T)/T;
}
# Normal distribution
if(mType==2){
  m1  <- mean(dy[,1]); s1 <- sd(dy[,1]); U  <- pdist("norm",(dy[,1]-m1)/s1); 
  m2  <- mean(dy[,2]); s2 <- sd(dy[,2]); V  <- pdist("norm",(dy[,2]-m2)/s2);
  UV  <- as.matrix(cbind(U,V))
}
# t-Student
if(mType==3){
  m1  <- mean(dy[,1]); s1 <- sd(dy[,1]); v1 <- 4 + 6/(kurtosis(dy[,1])-3); U  <- pdist("std",(dy[,1]-m1)/s1,shape = v1); 
  m2  <- mean(dy[,2]); s2 <- sd(dy[,2]); v2 <- 4 + 6/(kurtosis(dy[,2])-3); V  <- pdist("std",(dy[,2]-m2)/s2,shape = v2 );  
  UV  <- as.matrix(cbind(U,V))
}

par(mfrow=c(1,2), cex = 0.7, bty="l")
plot(coredata(dy),main = "Scatter plot for XY")
plot(UV, main = "Scatter plot for marginal distr.")


# Step 2. Fitting copula to UV
Cop1 <- normalCopula()
Cop2 <- tCopula(dispstr="un")
Cop3 <- gumbelCopula()
Cop4 <- claytonCopula()
Cop5 <- frankCopula()
Cop  <- Cop1

# ml estimation often breaks down...
# CopEst <- fitCopula(Cop, u, method = "itau")
sVal <- fitCopula(Cop, UV, method = "itau")@estimate
CopEst <- fitCopula(Cop, UV, method = "ml", start=sVal)
# For tCopula use the function below
# CopEst <- fitCopula(Cop2, UV, method = "itau.mpl", start=c(cor(dy)[1,2],5))
# CopEst <- fitCopula(Cop2, UV, method = "ml", start=c(cor(dy)[1,2],5))

CopEst@estimate
CopEst@loglik

# 8. VaR/ES
#####################################
H    = 10
Nsim = 10000
p    = 0.05     # tolerance level

# a. Copula
UVsim <- rCopula(Nsim*H,CopEst@copula)

# b. returning to obsevables
# Empirical marginal cdf 
if(mType==1){XYsim <- cbind(quantile(dy[,1],UVsim[,1]), quantile(dy[,2],UVsim[,2]))}
# Normal distribution
if(mType==2){XYsim <- cbind(qdist("norm",UVsim[,1])*s1+m1, qdist("norm",UVsim[,2])*s2+m2)}
# t-Student
if(mType==3){XYsim <- cbind(qdist("std",UVsim[,1],shape = v1)*s1+m1, qdist("std",UVsim[,2],shape = v2)*s2+m2)}

par(mfrow=c(1,2), cex = 0.7, bty="l")
plot(coredata(dy),  main='Returns')
plot(XYsim[1:dim(dy)[1],], col='red', main='Simulated')

# portfolio return
Rdraws     <- matrix(XYsim%*%w,H,Nsim)
temp       <- RdrawsToVaRES(RDraws,p)
VaRHcop    <- temp$VaR
ESHcop     <- temp$ES

VaRH <- cbind(VaRH,VaRHmnorm, VaRHcopN, VaRHcop)
ESH  <- cbind(ESH,ESHmnorm, ESHcopN, ESHcop)

save("VaRH","ESH", file = "VaRESH.RData")

kable(cbind(1:H,-VaRH), digits=2)
kable(cbind(1:H,-ESH), digits=2)