# copulaExamplesCarmona.ssc
#
# author: Eric Zivot
# created: January 28, 2004
# updated: March 4, 2005
#
# comments
# 1. Copula examples from Carmona (2004) chapter 2

module(finmetrics)

#
# 2.1 Multivariate Data
#

#
# Brazil and Columbia coffee data
# note: dates don't match data in Carmona's book. Carmona's data starts
# in January 10, 1986 and ends in January 1, 1999
# FM data starts in March 10, 1993 and ends in January 1, 1999.
#
colIds(columbia.coffee)
colIds(brazil.coffee)
columbia.coffee[1:5,]
brazil.coffee[1:5,]
start(brazil.coffee)
end(brazil.coffee)

# plot prices
plot(columbia.coffee[,"value"],brazil.coffee[,"value"],
reference.grid=F,main="Coffee futures",plot.args=list(col=1:2))
legend(0,280,legend=c("Columbia","Brazil"),lty=c(1,1),col=1:2)

# create bivariate timeSeries of log returns
# this series matches output from page 50
coffee.ts = seriesMerge(brazil.coffee[,"log.return"],
columbia.coffee[,"log.return"],)
colIds(coffee.ts) = c("brazil","columbia")
coffee.ts[1:5,]

# plot log return data
seriesPlot(coffee.ts,one.plot=F,type="l",
main="Log returns on coffee futures",)

# extract data as a matrix and create objects used by Carmona
coffee.mat = as.matrix(seriesData(coffee.ts))
BCofLRet = coffee.mat[,"brazil"]
CCofLRet = coffee.mat[,"columbia"]

#
# 2.1.1 Density Estimation
#

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

# Note: carmona's data UTIL.index is not in FM. cannot reproduce figure 2.1
# on page 53

#
# 2.2 Multivariate Normal Distributiono
# 

# pg. 59
# generate bivariate normal data with mean zero, sd=1 and rho=0.18
set.seed(0)
TSAMPLE <- rmvnorm(n=128, mean=rep(0,2), sd=rep(1,2), rho=.18)

# plot bivariate density

TDENS <- kde2d(TSAMPLE[,1],TSAMPLE[,2])
persp(TDENS$x,TDENS$y,TDENS$z)
wireframe(z ~ x*y, con2tr(TDENS),
   aspect = c(1, 0.5), screen = list(z=20, x=-60), zoom = 1.2)

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

#
# 2.2.4  Let's Have Some Coffee
#

# pg. 60, figure 2.4: plot coffee prices
# note: Carmona's plot starts in 1986 whereas FM data starts in 
# March of 1993
plot(columbia.coffee[,"value"],brazil.coffee[,"value"],
reference.grid=F,main="Coffee futures",plot.args=list(col=1:2))
legend(0,280,legend=c("Columbia","Brazil"),lty=c(1,1),col=1:2)

# pg. 61, Fig 2.5: scatterplot of coffee returns
plot(BCofLRet,CCofLRet,xlab="Brazil returns",
ylab="Columbia returns",
main="Daily Returns on Coffee Futures: 1993 - 1999")

# consider scatterplot with zero values removed
# note: Carmona finds PNZ = 0.408 which is based on larger sample
NZ <- (BCofLRet != 0 & CCofLRet != 0)
PNZ <- mean(NZ)
PNZ

BLRet <- BCofLRet[NZ]
CLRet <- CCofLRet[NZ]
plot(BLRet,CLRet,
xlab="Brazil returns",
ylab="Columbia returns",
main="Returns on Coffee Futures with Zero Returns Removed")

# compare scatterplots with and without zeros removed

par(mfrow=c(1,2))
	plot(BCofLRet,CCofLRet,xlab="Brazil returns",
	ylab="Columbia returns",
	main="Zeros Retained")
	plot(BLRet,CLRet,
	xlab="Brazil returns",
	ylab="Columbia returns",
	main="Zeros Removed")
par(mfrow=c(1,1))

#
# pg. 62, Section 2.2.5: Is the joint distribution normal?
#

summaryStats(coffee.ts[NZ,])

XLIM <- c(-.4,.3)
YLIM <- c(-.2,.4)
Mu <- c(mean(BLRet), mean(CLRet))
Sigma <- var(cbind(BLRet,CLRet))
N <- length(BLRet)

# simulate bivariate normal data calibrated to coffee data
# and compare to actual data
set.seed(0)
CNsim <- rmvnorm(N,mean=Mu,cov=Sigma)
# pg. 63, Fig 2.6
par(mfrow=c(1,2))
	plot(BLRet,CLRet, xlim=XLIM, ylim=YLIM)
	title("Actual Returns")
	plot(CNsim,xlim=XLIM, ylim=YLIM)
	title("Simulated Normal Returns")
par(mfrow=c(1,1))

#
# pg. 63, section 2.3: 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))
}

# pg. 64, Fig. 2.7
eda.shape(BLRet)
eda.shape(CLRet)


# look at tail plots - note tail option is not supported in FinMetrics
# want about 8% of data above upper threshold and below lower threshold
# compare with EVIS function shape

# pg. 65, Fig. 2.8. 
# MLE of gpd shape parameters for left tail
# note: results don't match Carmona due to different sample sizes
shape.plot(-BLRet)	
# MLE of gpd shape parameters for right tail
# note: results don't match Carmona due to different sample sizes
shape.plot(BLRet)	

# compare with EVIS function shape

# pg. 66, Fig. 2.9 fit gpd distribution by ml
# note: qq plots should be linear on the left
# results don't match Carmona
B.est <- gpd.tail(BLRet, upper = 0.04, lower = -0.045)
class(B.est)
names(B.est)
# note: no summary method; however a plot method
B.est
# no standard errors are computed - this is a bug
plot(B.est)
# note - plots are only given for upper tail estimates

# compare with EVIS function gpd
# both gpd.tail and gpd return objects of class "gpd" so they should
# give the same results - they do!


# pg. 67, Fig. 2.10: examine tail plots on log scale
par(mfrow=c(1,2))
	tailplot(B.est,tail="lower")
	title("Lower tail Fit")
	tailplot(B.est,tail="upper")
	title("Upper tail Fit")
par(mfrow=c(1,1))

# do the same for the columbia data
# results don't match Carmona
par(mfrow=c(1,2))
	shape.plot(-CLRet)	# lower tail
	shape.plot(CLRet)	# upper tail
par(mfrow=c(1,1))

# note: qq plots should be linear on the left
C.est <- gpd.tail(CLRet, upper = 0.0375, lower = -0.03)
C.est

# examine tail plots on log scale
par(mfrow=c(1,2))
	tailplot(C.est,tail="lower")
	title("Lower tail Fit")
	tailplot(C.est,tail="upper")
	title("Upper tail Fit")
par(mfrow=c(1,1))

# pg. 67, MC simulations
# use gpd.2q to compute simulation from fitted 2-tailed gpd model
# note: gdp.1q is used for 1-tail gpd model
args(gpd.2q)
set.seed(123)
BLRet.sim <- gpd.2q(runif(length(BLRet)),B.est)
CLRet.sim <- gpd.2q(runif(length(CLRet)),C.est)

# do QQ-plot of actual data against simulated data
# pg. 68, Fig. 2.11
qqplot(BLRet,BLRet.sim)
abline(0,1)
qqplot(CLRet,CLRet.sim)
abline(0,1)

# scatterplot of simulated returns and actual returns
# pg. 68, Fig. 2.12
plot(BLRet.sim,CLRet.sim,
main="Simulated returns")
plot(BLRet,CLRet,
main="Actual returns")


#
# 2.3.2 More measures of dependence
#

# compute Kendall's tau using splus function cor.test
# strangely, result matches Carmona, pg. 69 even though sample size
# is different
?cor.test
coffee.tau = cor.test(BLRet,CLRet,method="k")$estimate
coffee.tau

# compute Spearman's rho using splus function cor.test
# strangely, result matches Carmona, pg. 69 even though sample size
# is different

coffee.rho = cor.test(BLRet,CLRet,method="spearman")$estimate
coffee.rho

#
# 2.4 Copulas and Random Simulations
#

# pg. 70, Fig 2.13
# compute cdf for coffee data based on fitted 2-tail gdp models
# these cdf estimates are uniformly distributed but dependent
U <- gpd.2p(BLRet, B.est)
V <- gpd.2p(CLRet, C.est)
plot(U,V,main="Dependence between coffee returns")
# compute empirical copula
EMPCOP <- empirical.copula(U,V)
class(EMPCOP)
slotNames(EMPCOP)
# more on empirical.copula below in section on fitting copulas

# do the same for the simulated data
U.sim <- gpd.2p(BLRet.sim, B.est)
V.sim <- gpd.2p(CLRet.sim, C.est)
plot(U.sim,V.sim, main="Dependence between simulated returns")

# remark: could one also use the empirical CDF and do the same thing?

# 
# 2.4.2 First Examples of Copula Families
#

# create normal copula object
# note: this is an SV4 object
ncop.7 = normal.copula(0.7)
class(ncop.7)
slotNames(ncop.7)
# 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.

# illustrate functions
ncop.7
Kendalls.tau(ncop.7)
Spearmans.rho(ncop.7)

contour.plot(ncop.7)
persp.dcopula(ncop.7) # fig 2.14
contour.dcopula(ncop.7)
persp.pcopula(ncop.7) # fig 2.14
contour.pcopula(ncop.7) # same as contour.plot?

# pg. 73, Fig 2.14
# surface plot of Gaussian copula with rho=0.7
# first plot CDF, then density
persp.pcopula(ncop.7)
persp.dcopula(ncop.7)

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

# illustrate functions
gcop.1.5
Kendalls.tau(gcop.1.5)
Spearmans.rho(gcop.1.5)

contour.plot(gcop.1.5)
persp.dcopula(gcop.1.5) # fig 2.15
contour.dcopula(gcop.1.5)
persp.pcopula(gcop.1.5) # fig 2.15
contour.pcopula(gcop.1.5) # same as contour.plot?

# pg. 73, Fig 2.14
# surface plot of Gaussian copula with rho=0.7
# first plot CDF, then density
# note: error in CDF plot: z-axis should be "CDF"
persp.pcopula(gcop.1.5)
persp.dcopula(gcop.1.5)


#
# 2.4.3 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.copula)
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=1.4, fx = N(3,4^2),
# and fy = Student-t with df=3
BIV1 <- bivd(cop=gumbel.copula(1.4), marginX="norm", marginY="t", 
param.marginX=c(3,4), param.marginY=3)
class(BIV1)
slotNames(BIV1)
# 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.

# pg. 76, Fig. 2.16 show bivariate density
persp.dbivd(BIV1)
contour.dbivd(BIV1)
# show bivariate CDF
persp.pbivd(BIV1)
contour.pbivd(BIV1)

# pg 77, Fig 2.17
# note: can't do figure below because FM uses trellis plotting function
# contourplot

par(mfrow=c(1,2))
	contour.dbivd(BIV1)
	title("Density of the Bivariate Distribution BIV1")
	contour.pbivd(BIV1)
	title("CDF of the Bivariate Distribution BIV1")
par(mfrow=c(1,1))

#
# 2.2.2 Fitting copulas
#

# create empirical copula object for coffee data
?empirical.copula
args(empirical.copula)
EMPCOP <- empirical.copula(U,V)
class(EMPCOP)
slotNames(EMPCOP)
# 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.

# contour plot of CDF
contour.plot(EMPCOP)
tail.index(EMPCOP)
Kendalls.tau(EMPCOP)
Spearmans.rho(EMPCOP)
persp.dcopula(EMPCOP)
contour.dcopula(EMPCOP)
persp.pcopula(EMPCOP)
contour.pcopula(EMPCOP)


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

# fit Gumbel copula to coffee data
ESTC <- fit.copula(EMPCOP,family="gumbel")
# note hessian failed to invert
class(ESTC)
names(ESTC)
methods(class="copula.fit")
ESTC$copula
IC(ESTC)

# pg. 79
# compare Kendall's tau and Spearman's rho for empirical copula
# and fitted copula
# note: results are not the same as Carmona
Kendalls.tau(EMPCOP)
# note: cannot pass ESTC directly to Kendalls.tau
Kendalls.tau(ESTC$copula)
Spearmans.rho(EMPCOP)
# note: cannot pass ESTC directly to Spearmans.rho
Spearmans.rho(ESTC$copula)



#
# 2.4.5 Monte Carlo Simulations with Copulas
#


# simulate 100 observations from bivariate normal copula
ncop.7.sim = rcopula(ncop.7,1000)
class(ncop.7.sim)
names(ncop.7.sim)
plot(ncop.7.sim$x,ncop.7.sim$y)

# pg. 79
# simulate bivariate data from Gumbel copula fitted to 
# coffee data with gpd marginals using the function rcopula

N <- length(BLRet)
# note: cannot pass ESTC directly as done by Carmona; must pass ESTC$copula
SD <- rcopula(ESTC$copula,N)
Xsim <- gpd.2q(SD$x, B.est)
Ysim <- gpd.2q(SD$y, C.est)

# Fig. 2.19
# compare simulated data with actual data
par(mfrow=c(1,2))
	plot(BLRet,CLRet,main="Actual Returns")
	plot(Xsim,Ysim,main="Simulated Returns")
par(mfrow=c(1,1))

# do the same plot but overlay acutal and simulated data
plot(BLRet,CLRet)
points(Xsim,Ysim,col=5)