# copulaExamples.ssc				Examples for copula chapter
#
# author: Eric Zivot
# created: January 28, 2004
# updated: March 30, 2006
# updated: July 10, 2007
#	incorporated changes in copula plotting functions

module(finmetrics)

#
# 1. Motivating example
#

# BMW and Siemans data from EVIS library
start(bmw)
end(bmw)

# follow carmona and remove zeros to avoid jumps in empirical CDFs
# determine joint occurance of zeros
zero.idx = (bmw == 0 & siemens == 0)
sum(zero.idx)/seriesLength(bmw)
# 5% of data have joint zeros
nz.idx = (seriesData(bmw) != 0 & seriesData(siemens) != 0)

# create bivariate timeSeries and plot 
german.ts = seriesMerge(bmw,siemens)
colIds(german.ts)
seriesPlot(german.ts,one.plot=F,type="l", ref.line=F)

par(mfrow=c(2,1))
	plot(bmw, main="Daily returns on BMW", reference.grid=F)
	abline(h=0)
	plot(siemens, main="Daily returns on Siemens", reference.grid=F)
	abline(h=0)
par(mfrow=c(1,1))

# test for normality
qqPlot(german.ts[nz.idx,],strip=colIds(german.ts))
summaryStats(german.ts[nz.idx,])
normalTest(german.ts[nz.idx,],method="jb")

# scatterplot with zeros
plot(seriesData(bmw),seriesData(siemens),
xlab="BMW",ylab="Siemens")
abline(h=0,v=0)

# scatterplot without zeros
plot(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]),
xlab="BMW",ylab="Siemens")
abline(h=0,v=0)
par(mfrow=c(1,1))

# simulate multivariate normal data calibrated to data
cor(german.ts[nz.idx,])[1,2]

# without zeros
summaryStats(german.ts[nz.idx,])
cor(german.ts[nz.idx,])

XLIM <- c(min(bmw[nz.idx]),max(bmw[nz.idx]))
YLIM <- c(min(siemens[nz.idx]),max(siemens[nz.idx]))
mu.hat = colMeans(german.ts[nz.idx,])
Sigma.hat = var(german.ts[nz.idx,])
nobs = numRows(german.ts[nz.idx,])

# simulate bivariate normal data calibrated to german.ts
# and compare to actual data
set.seed(0)
german.sim <- rmvnorm(nobs,mean=mu.hat,cov=Sigma.hat)
colIds(german.sim) = colIds(german.ts)
plot(german.sim,xlim=XLIM, ylim=YLIM)
abline(h=0,v=0)

par(mfrow=c(1,2))
	plot(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]),
	xlab="BMW",ylab="Siemens",
	xlim=XLIM, ylim=YLIM)
	title("Actual Returns")
	abline(h=0,v=0)
	plot(german.sim,xlim=XLIM, ylim=YLIM)
	title("Simulated Normal Returns")
	abline(h=0,v=0)
par(mfrow=c(1,1))


# Note: Carmona's function kdest is not in FM
# use MASS function kde2d to compute bivariate density estimate
library(MASS)
args(kde2d)

# kernel density estimate for coffee data
DENS = kde2d(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]))
wireframe(z ~ x*y, con2tr(DENS),
   aspect = c(1, 1), screen = list(z=20, x=-60), zoom = 1.2)


DENS = kde2d(german.sim[,1],german.sim[,2])
wireframe(z ~ x*y, con2tr(DENS),
   aspect = c(1, 1), screen = list(z=20, x=-60), zoom = 1.2)

# Marginals and more measures of dependence

# From Carmona pg. 425
eda.shape <- function(x)
{
  par(mfrow = c(2,2))
  hist(x)
  boxplot(x)
  iqd <- summary(x)[5] -summary(x)[2]
  plot(density(x,width=2*iqd),xlab="x",ylab="",type="l")
  qqnorm(x)
  qqline(x)
  par(mfrow=c(1,1))
}

# show distributions
eda.shape(seriesData(bmw[nz.idx]))
eda.shape(seriesData(siemens[nz.idx]))

# mean excess plots for left and upper tails
par(mfrow=c(2,2))
	meplot(-bmw[nz.idx])
	title("BMW, lower tail")
	meplot(-siemens[nz.idx])
	title("Siemens, lower tail")
	meplot(bmw[nz.idx])
	title("BMW, upper tail")
	meplot(siemens[nz.idx])
	title("Siemens, upper tail")
par(mfrow=c(1,1))


# MLE of gpd shape parameters for lower and upper tails
par(mfrow=c(2,2))
	shape.plot(-bmw[nz.idx])	
	legend(0,0.75,"BMW, lower tail")
	shape.plot(-siemens[nz.idx])	
	legend(0,0.75,"Siemens, lower tail")
	shape.plot(bmw[nz.idx])	
	legend(0,0.6, "BMW, upper tail")
	shape.plot(siemens[nz.idx])
	legend(0,0.6,"Siemens, upper tail")
par(mfrow=c(1,1))

# fit 2-tailed gpd distribution by ml
# note: qq plots should be linear on the left
gpd.bmw.fit2 = gpd.tail(bmw[nz.idx], upper = 0.015, lower = -0.015)
gpd.bmw.fit2
# note - estimates of xi are similar for lower and upper tails

gpd.siemens.fit2 = gpd.tail(siemens[nz.idx], upper = 0.01, lower = -0.01,
upper.method="lmom")
gpd.siemens.fit2
# note - estimates of xi are similar for lower and upper tails
# upper.method="lmom" since mle did not converge

par(mfrow=c(2,2))
	tailplot(gpd.bmw.fit2,tail="lower")
	title("BMW Lower tail Fit")
	tailplot(gpd.bmw.fit2,tail="upper")
	title("BMW Upper tail Fit")
	tailplot(gpd.siemens.fit2,tail="lower")
	title("Siemens Lower tail Fit")
	tailplot(gpd.siemens.fit2,tail="upper")
	title("Siemens Upper tail Fit")
par(mfrow=c(1,1))

# MC simulations from fitted gpd distributions assuming no dependence structure
# use gpd.2q to compute simulation from fitted 2-tailed gpd model
args(gpd.2q)
nobs = numRows(bmw[nz.idx])
set.seed(123)
bmw.gpd.sim <- gpd.2q(runif(nobs),gpd.bmw.fit2)
siemens.gpd.sim <- gpd.2q(runif(nobs),gpd.siemens.fit2)

# do QQ-plot of actual data against simulated data

par(mfrow=c(1,2))
	qqplot(seriesData(bmw[nz.idx]),bmw.gpd.sim,
	xlab="Actual returns", ylab="Simulated returns")
	abline(0,1)
	title("BMW")
	qqplot(seriesData(siemens[nz.idx]),siemens.gpd.sim,
	xlab="Actual returns", ylab="Simulated returns")
	abline(0,1)
	title("Siemens")
par(mfrow=c(1,1))

# scatterplot of simulated returns with no dependence structure
plot(bmw.gpd.sim,siemens.gpd.sim)
abline(h=0,v=0)


#
# 2 estimate dependence functions
#

# compute Kendall's tau using splus function cor.test
# works on timeSeries!
?cor.test
german.tau = cor.test(bmw[nz.idx],siemens[nz.idx],
             method="k")$estimate
german.tau

# compute Spearman's rho using splus function cor.test
german.rho = cor.test(bmw[nz.idx],siemens[nz.idx],
             method="spearman")$estimate
german.rho

#
# 3 Copula constructor functions
#

# create normal copula object
# note: this is an SV4 object
ncop.7 = normal.copula(delta=0.7)
class(ncop.7)
slotNames(ncop.7)
inherits(ncop.7, what="copula")
slotNames(ncop.7)
ncop.7@parameters
ncop.7@param.lowbnd
ncop.7@param.upbnd
ncop.7@message
ncop.7

# create Gumbel copula object
# note: this is an SV4 object
gcop.2 = gumbel.copula(2)
class(gcop.2)
slotNames(gcop.2)

gcop.2 = ev.copula(family="gumbel", param=2)
class(gcop.2)
gcop.2


# available methods are:
# print, pcopula, dcopula, rcopula, 
# contour.plot, Kendalls.tau, Spearmans.rho.
# available function operating on normal copula objects
# persp.dcopula, persp.pcopula, contour.dcopula, contour.pcopula.

#
# Visualization functions
#

# normal copula with delta=0.7
contour.plot(ncop.7)
persp.dcopula(ncop.7) 
contour.dcopula(ncop.7)
persp.pcopula(ncop.7) 
contour.pcopula(ncop.7) 

# create 4 panel plot
# note: the functions junk, junk2, junk3, and junk4 are modified versions
# of the functions persp.dcopula, contour.dcopula, persp.pcopula, and
# contour.pcopula that allow for easier printing

print.trellis(junk2(ncop.7), split = c(1,1,2,2), more = T)  # bot left
print.trellis(junk(ncop.7), split = c(1,2,2,2), more = T)   # top left
print.trellis(junk4(ncop.7), split = c(2,1,2,2), more = T)  # bot right
print.trellis(junk3(ncop.7) ,split = c(2,2,2,2))            # top right

# EZ: 7/10/057. plotting functions have been updated so that they
# optionally return trellis objects for easy printing

print.trellis(contour.dcopula(ncop.7, value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(persp.dcopula(ncop.7, value="trellis"), split = c(1,2,2,2), more = T)   # top left
print.trellis(contour.pcopula(ncop.7, value="trellis"), split = c(2,1,2,2), more = T)  # bot right
print.trellis(persp.pcopula(ncop.7, value="trellis") ,split = c(2,2,2,2))            # top right


# normal copula with delta=0.0.0001
normal.cop.0 = normal.copula(delta=0.0001)
print.trellis(contour.dcopula(normal.cop.0, value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(persp.dcopula(normal.cop.0, value="trellis"), split = c(1,2,2,2), more = T)   # top left
print.trellis(contour.pcopula(normal.cop.0, value="trellis"), split = c(2,1,2,2), more = T)  # bot right
print.trellis(persp.pcopula(normal.cop.0, value="trellis") ,split = c(2,2,2,2))            # top right

# Gumbel copula with delta=2
contour.plot(gcop.2)
persp.dcopula(gcop.2) 
contour.dcopula(gcop.2)
persp.pcopula(gcop.2) 
contour.pcopula(gcop.2) 

print.trellis(contour.dcopula(gcop.2,value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(persp.dcopula(gcop.2, value="trellis"), split = c(1,2,2,2), more = T)   # top left
print.trellis(contour.pcopula(gcop.2, value="trellis"), split = c(2,1,2,2), more = T)  # bot right
print.trellis(persp.pcopula(gcop.2, value="trellis") ,split = c(2,2,2,2))            # top right

# Clayton copula with delta = 2

ks.cop.2 = kimeldorf.sampson.copula(delta=2)
print.trellis(contour.dcopula(ks.cop.2, value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(persp.dcopula(ks.cop.2, value="trellis"), split = c(1,2,2,2), more = T)   # top left
print.trellis(contour.pcopula(ks.cop.2, value="trellis"), split = c(2,1,2,2), more = T)  # bot right
print.trellis(persp.pcopula(ks.cop.2, value="trellis") ,split = c(2,2,2,2))            # top right


# plot contours of indep, min and max copulas
max.cop = normal.copula(0.999)
indep.cop = normal.cop.0
min.cop = kimeldorf.sampson.copula(delta=0.001)

print.trellis(contour.dcopula(max.cop, value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(contour.dcopula(indep.cop,value="trellis"), split = c(1,2,2,2), more = T)  # top left
print.trellis(contour.dcopula(min.cop,value="trellis"), split = c(2,2,2,2))            # top right




#
# 4 Measures of dependence
#

# normal copula with delta = 0.7
Kendalls.tau(ncop.7)
Spearmans.rho(ncop.7)
tail.index(ncop.7)

# gumbel copula with delta = 2
Kendalls.tau(gcop.2)
Spearmans.rho(gcop.2)
tail.index(gcop.2)

# gumbel copula with delta = 10
gcop.10 = gumbel.copula(delta=10)
Kendalls.tau(gcop.10)
Spearmans.rho(gcop.10)
tail.index(gcop.10)

#
# 5 Simulation of uniform variables with copula joint CDF 
#

set.seed(123)
normal.7.sim = rcopula(ncop.7,2000)
normal.0.sim = rcopula(normal.cop.0,2000)
gumbel.2.sim = rcopula(gumbel.cop.2,2000)
ks.2.sim = rcopula(ks.cop.2,2000)

par(mfrow=c(2,2))
	plot(normal.0.sim, main="Independent Copula",
	ylab="u",xlab="v")
	plot(normal.7.sim, main="Normal Copula, delta=0.7",
	ylab="u",xlab="v")
	plot(gumbel.2.sim, main="Gumbel Copula, delta=2",
	ylab="u",xlab="v")
	plot(ks.2.sim, main="Clayton Copula, delta=2",
	ylab="u",xlab="v")
par(mfrow=c(1,1))



#
# 6 Copulas and General Bivariate Distributions
#

# use bivd to create bivariate distribution from copula
# and two marginals
# note: In FM help the function bivd is called bivd.copula, however
# bivd.copula does not exist.
args(bivd)

# pg. 75-76
# the marginal densities must computed with available splus functions
# starting with "d". For example, if marginX="norm" then there must be a 
# function available called "dnorm" that computes density values for the
# normal distribution. The parameters in param.marginX must 
# match the optional arguments to the "d" dist function. For example, 
# dnorm has optional argumens mean and sd. Therefore, param.marginX is
# a vector with values for mean and sd.

# create bivariate density using gumbel copula with beta=2, fx = N(3,4^2),
# and fy = Student-t with df=3
gumbel.biv <- bivd(cop=gumbel.copula(2), Xmarg="norm", Ymarg="t", 
                   param.Xmarg=c(3,4), param.Ymarg=3)
class(gumbel.biv)
slotNames(gumbel.biv)
gumbel.biv@copula
# note: no show method

# The following functions are implemented for all child classes that 
# inherit from the class "bivd": 
# pbivd, dbivd, rbivd, persp.dbivd, persp.pbivd, contour.dbivd, 
# contour.pbivd.

# note: make junk functions for multiple plotting
# continue making changes to plotting functions
persp.dbivd(gumbel.biv)
contour.dbivd(gumbel.biv)
persp.pbivd(gumbel.biv)
contour.pbivd(gumbel.biv)

# use junk functions for 4 panel plot

print.trellis(persp.dbivd(gumbel.biv, value="trellis"), split = c(1,1,2,2), more = T)  # bot left
print.trellis(contour.dbivd(gumbel.biv, value="trellis"), split = c(1,2,2,2), more = T)   # top left
print.trellis(persp.pbivd(gumbel.biv, value="trellis"), split = c(2,1,2,2), more = T)  # bot right
print.trellis(contour.pbivd(gumbel.biv, value="trellis") ,split = c(2,2,2,2))            # top right

# compute joint orthant probabilities

pbivd(gumbel.biv,
      x = c(-2, 0),
      y = c(-2, 0))


# create random simulation from bivariate distribution

set.seed(123)
gumbel.biv.sim = rbivd(gumbel.biv,n=2000)
class(gumbel.biv.sim)
names(gumbel.biv.sim)
U.x = pnorm(gumbel.biv.sim$x,mean=3,sd=4)
V.y = pt(gumbel.biv.sim$y,df=3)

par(mfrow=c(2,2))
	hist(gumbel.biv.sim$x, xlab="x")
	hist(gumbel.biv.sim$y, xlab="y")
	plot(as.matrix(gumbel.biv.sim))
	plot(U.x,V.y,xlab="u",ylab="v")
par(mfrow=c(1,1))

#
# 7 simulations from arbitrary distributions
#

# simulate bivariate data from Gumbel copula fitted to 
# return data with gpd marginals using the function rcopula

n.obs <- seriesLength(bmw[nz.idx])
set.seed(123)
u.data <- rcopula(ncop.7,n.obs)
bmw.cop.sim <- gpd.2q(u.data$x, gpd.bmw.fit2)
siemens.cop.sim <- gpd.2q(u.data$y, gpd.siemens.fit2)

# Fig. 2.19
# compare simulated data with actual data
par(mfrow=c(1,2))
	plot(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]),
	     main="Actual Returns")
	plot(bmw.cop.sim,siemens.cop.sim,main="Simulated Returns")
par(mfrow=c(1,1))

par(mfrow=c(1,2))
	plot(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]),
	xlab="BMW",ylab="Siemens",
	xlim=XLIM, ylim=YLIM)
	title("Actual Returns")
	abline(h=0,v=0)
	plot(bmw.cop.sim,siemens.cop.sim, 
	xlab="BMW",ylab="Siemens",
	xlim=XLIM, ylim=YLIM)
	title("Simulated Returns from Copula")
	abline(h=0,v=0)
par(mfrow=c(1,1))


# do the same plot but overlay acutal and simulated data
plot(seriesData(bmw[nz.idx]),seriesData(siemens[nz.idx]),
	 xlab="BMW",ylab="Siemens")
points(bmw.cop.sim,siemens.cop.sim, col=5)


#
# 8 Fitting copulas
#

# compute empirical copula for return data

# compute cdf for returns data based on fitted 2-tail gdp models
# these cdf estimates are uniformly distributed but dependent
U.bmw.gpd <- gpd.2p(bmw[nz.idx], gpd.bmw.fit2)
V.siemens.gpd <- gpd.2p(siemens[nz.idx], gpd.siemens.fit2)
plot(U.bmw.gpd,V.siemens.gpd)


# do the same for the simulated data
U.sim <- gpd.2p(bmw.gpd.sim, gpd.bmw.fit2)
V.sim <- gpd.2p(siemens.gpd.sim, gpd.siemens.fit2)
plot(U.sim,V.sim)


empcop.bs <- empirical.copula(U.bmw.gpd,V.siemens.gpd)
class(empcop.bs)
slotNames(empcop.bs)

# The "empirical.copula" class of objects currently 
# has methods for the following generic functions: 
# Afunc, LAMBDA, pcopula, contour.plot, tail.index, 
# Kendalls.tau, Spearmans.rho.
# Furthermore, the following functions are implemented for this class: 
# persp.dcopula, persp.pcopula, contour.dcopula, contour.pcopula.


plot(empcop.bs@x,empcop.bs@y) 
contour.pcopula(empcop.bs) 
persp.dcopula(empcop.bs) 
contour.dcopula(empcop.bs)

# put 4 plots on same graph
print.trellis(junk2(empcop.bs), split = c(1,1,2,2), more = T)  # bot left
print.trellis(xyplot(U.bmw.gpd~V.siemens.gpd), split = c(1,2,2,2), more = T)   # top left
print.trellis(junk4(empcop.bs), split = c(2,1,2,2), more = T)  # bot right
print.trellis(junk3(empcop.bs) ,split = c(2,2,2,2))            # top right

Kendalls.tau(empcop.bs)
Spearmans.rho(empcop.bs)
tail.index(empcop.bs)

# fit.copula is used for parametric fitting of copulas
?fit.copula
args(fit.copula)

# fit normal copula to german returns
cop.normal.fit = fit.copula(empcop.bs, family="normal",
                            plot=T)
class(cop.normal.fit)
methods(class="copulaFit")
cop.normal.fit
# compute Kendall's tau and Spearmans rho
Kendalls.tau(cop.normal.fit$copula)
Spearmans.rho(cop.normal.fit$copula)

# simulate return data from fitted copula
# do the method in two ways

normal.gpd.biv <- bivd(cop=cop.normal.fit$copula, 
                       Xmarg="gpd", Ymarg="t", 
                   param.Xmarg=c(3,4), param.Ymarg=3)



# fit Gumbel copula to return data
cop.gumbel.fit <- fit.copula(empcop.bs,family="gumbel",
                             plot=T)
cop.gumbel.fit

# compare fits
compare.copulaFit(cop.normal.fit, cop.gumbel.fit)

#
# 9 Risk Management with Copulas
#

# simulate portfolio returns from bivariate model with
# normal margins and Clayton copula:
# X ~ N(0.00034,0.01476^2)
# Y ~ N(0.00021, 0.01140^2)
# C ~ Clayton(delta=2) 

# create bivd object
clayton.biv <- bivd(cop=kimeldorf.sampson.copula(delta=2), 
                    Xmarg = "norm", Ymarg = "norm", 
                    param.Xmarg=c(0.00034, 0.0147), 
                    param.Ymarg=c(0.00021, 0.01140))
# simulate from bivariate distn and compute log-return
set.seed(123)
xy.sim = rbivd(clayton.biv,10000)
R.sim = log(0.7*exp(xy.sim$x) + 0.3*exp(xy.sim$y))
VaR.95.est = quantile(-R.sim,probs=0.95) 
ES.95.est = -mean(R.sim[-R.sim > VaR.95.est])
VaR.95.est
ES.95.est

# repeat analysis with normal distn
rho.hat = cor(as.matrix(xy.sim))[1,2]
normal.biv <- bivd(cop=normal.copula(delta=rho.hat), 
                    Xmarg = "norm", Ymarg = "norm", 
                    param.Xmarg=c(0.00034, 0.0147), 
                    param.Ymarg=c(0.00021, 0.01140))
# simulate from bivariate distn and compute log-return
set.seed(123)
xy.sim.norm = rbivd(normal.biv,10000)
R.sim.norm = log(0.7*exp(xy.sim.norm$x) + 0.3*exp(xy.sim.norm$y))
VaR.95.est.norm = quantile(R.sim.norm,probs=0.05) 
-VaR.95.est.norm
ES.95.est.norm = mean(R.sim.norm[R.sim.norm < VaR.95.est.norm])
-ES.95.est.norm

# compute VaR on 2 asset portfolio
args(VaR.exp.sim)
VaR.out = VaR.exp.sim(n=10000,Q=c(0.01,0.05),
                      copula=cop.normal.fit$copula,
                      x.est=gpd.bmw.fit2,
                      y.est=gpd.siemens.fit2,
                      lambda1=0.7,
                      lambda2=0.3)
VaR.out


# compute cumulative probabilities from fitted bivariate distributions
# compare joint tail probabilities of copula fit with bivariate normal

# compute 95% and 99% quantiles
bmw.highq = gpd.2q(p=c(0.95,0.99), 
                   obj=gpd.bmw.fit2, linear=T)
bmw.highq
bmw.lowq = gpd.2q(p=c(0.05,0.01), 
                   obj=gpd.bmw.fit2, linear=T)
bmw.lowq
siemens.highq = gpd.2q(p=c(0.95,0.99), 
                   obj=gpd.siemens.fit2, linear=T)
siemens.highq
siemens.lowq = gpd.2q(p=c(0.05,0.01), 
                   obj=gpd.siemens.fit2, linear=T)
siemens.lowq

gpdjoint.2p(x=bmw.lowq, y=siemens.lowq, 
            copula=cop.normal.fit$copula, 
            x.est=gpd.bmw.fit2,
            y.est=gpd.siemens.fit2)