Chapter 3

Example

library(copula)
sigma <-0.3
kopula.norm <- normalCopula(param = sigma, dim = 2)
par(mfrow=c(2,2))
contour(kopula.norm, dCopula,  nlevels=50)
contour(kopula.norm, pCopula,  nlevels=20)
persp(kopula.norm,dCopula, xlab="u1", ylab="u2", zlab="c(u1, u2)")
persp(kopula.norm,pCopula, xlab="u1", ylab="u2", zlab="C(u1, u2)")

Example

library(statmod)
mvd<- mvdc(kopula.norm, c("invgauss", "norm"),
list(list(mean=350,dispersion=4.5), list(mean = 0, sd =0.2)))

par(mfrow=c(2,2))
persp(mvd, dMvdc, col=4, xlim = c(0, 2), ylim=c(-0.5, 0.5),
xlab="x1", ylab="x2", zlab="f(x1, x2)")
persp(mvd, pMvdc, col=4, xlim = c(0, 2), ylim=c(-0.5, 0.5),
xlab="x1", ylab="x2", zlab="f(x1, x2)")
contour(mvd, dMvdc, col=4, xlim = c(0, 2), ylim=c(-0.5, 0.5), 
xlab="x1", ylab="x2",nlevels=30)
contour(mvd, pMvdc, col=4, xlim = c(0, 2), ylim=c(-0.5, 0.5), 
xlab="x1", ylab="x2",nlevels=30)

Example

library(fitdistrplus)

sim <- rMvdc(1000,mvd)

fit.X1<-fitdist(sim[,1], dinvgauss, start= list(mean=350,dispersion=4.5), method="mle")
(mean.invgauss<-coef(fit.X1)[1]) 

##     mean 
## 167.3121

(dispersion.invgauss<-coef(fit.X1)[2])

## dispersion 
##   4.527281

#plot(fit.X1)

fit.X2<-fitdist(sim[,2], "norm", method="mle")
(mean.norm<-coef(fit.X2)[1]) 

##       mean 
## 0.01122929

(sd.norm<-coef(fit.X2)[2])

##        sd 
## 0.2029674

u.1<-pinvgauss(sim[,1],mean=mean.invgauss,dispersion=dispersion.invgauss)
u.2<-pnorm(sim[,2], mean=mean.norm, sd=sd.norm)
U<-cbind(u.1, u.2)

#Figure
par(mfrow=c(1,2))
plot(sim, xlim = c(0, 2), ylim=c(-0.5, 0.5),xlab=expression(x[1]), ylab=expression(x[2]))
plot(U, xlim = c(0, 1), ylim=c(0, 1) ,xlab=expression(u[1]), ylab=expression(u[2]))

fit.normal<- fitCopula(normalCopula(dim = 2),U, method="ml")
summary(fit.normal)

## Call: fitCopula(copula, data = data, method = "ml")
## Fit based on "maximum likelihood" and 1000 2-dimensional observations.
## Normal copula, dim. d = 2 
##       Estimate Std. Error
## rho.1   0.2895      0.028
## The maximized loglikelihood is 44.93 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        6        6

(sigma.fit<-coef(fit.normal))

##     rho.1 
## 0.2895408

gofCopula(normalCopula(sigma.fit), sim)

## 
##  Parametric bootstrap-based goodness-of-fit test of Normal copula,
##  dim. d = 2, with 'method'="Sn", 'estim.method'="mpl":
## 
## data:  x
## statistic = 0.017616, parameter.rho.1 = 0.29468, p-value = 0.517