# 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)