Fit partition of unity copulas in R

Rnw

\documentclass{article}

\usepackage{etex} % Avoid an error due to a lack of registers
\usepackage[USenglish]{babel} % Defines the language of macros as well
\usepackage{mathptmx,amsmath}
%\usepackage[utf8]{inputenc}
\usepackage{bm}
\usepackage{caption}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{multicol}
\usepackage{dsfont}
%\usepackage[demo]{graphicx}
\usepackage[font=small,skip=2pt]{caption}
%\usepackage{subcaption}
\usepackage[round]{natbib}



\begin{document}


<<cache=FALSE,echo=FALSE,message=FALSE>>=
library("knitr")
library("tikzDevice")
opts_chunk$set(
   dev='cairo_pdf',
   dev.args=list(family='DejaVu Sans')
   #out.width='5in', 
   #fig.width=5,
   #fig.height=5/sqrt(2),
   #fig.path='figures-knitr/',
   #fig.align='center',
)
@


\section{Data}

This R-script follows the approach for fitting a partition-of-unity copula from \cite{partofun}.

Data set treated in \cite{cottin2014bernstein} is created and their ranks and relative ranks are calculated:

<<data, echo=TRUE>>=

# data
mdata <- cbind(c(0.468,9.951,0.866,6.731,1.421,2.040,2.967,1.200,0.426,1.946,
                 0.676,1.184,0.960,1.972,1.549,0.819,0.063,1.280,0.824,0.227)
             , c(0.966,2.679,0.897,2.249,0.956,1.141,1.707,1.008,1.065,1.162,
                 0.918,1.336,0.933,1.077,1.041,0.899,0.710,1.118,0.894,0.837)
             )
param.faktor <- 1
rdata <- apply(mdata,2,rank)/param.faktor
param.m <- dim(rdata)[1]
udata <- rdata/(1+param.m)

@
Have a look at the data:
<<datatable, echo=FALSE>>=

#cbind(i=1:20, x_i=mdata[,1], y_i=mdata[,2], r_1i=rdata[,1],r_2i=rdata[,2])

kable(cbind(i=1:20, x_i=mdata[,1], y_i=mdata[,2], r_1i=rdata[,1],r_2i=rdata[,2]), digits=2)

@
Also plot the original data and the relative rank data.
<<plotorig, echo=FALSE, fig.lp="fig:", fig.cap = 'graph of original data (left) and rank vectors (right)'>>=

op <- par(mfrow = c(1, 2), pty = "s")  
plot(mdata,xlim=c(0,12),ylim=c(0,3), col="red")

#library(ggplot2)
#ggplot(as.data.frame(mdata), aes(V1, V2)) +
#    geom_point(shape=1) +     # Use hollow circles
#    expand_limits(y=0) +
#    expand_limits(y=3) +
#    expand_limits(x=12)
plot(udata,xlim=c(0,1),ylim=c(0,1), col="red")
par(op)
rm(op)

@
Get empirical correlations for rank data and for relative ranks
<<empcorr>>=

cor(rdata)
cor(udata)

@

\section{Simulation}

Set number of simulations and copula parameters:
<<paramset, echo=1:4, warnings=FALSE>>=

sims.n <- 5000
par.a <- 17
par.b <- 22
bin.par.a <- 22
bin.par.b <- 27

library(copula)

@
The algorithm now is:
<<algo1, dev='cairo_pdf', label='algorithm1'>>=

#simulation
# 1) choose pair(s) of empirical ranks randomly
rsims.index <- ceiling(runif(sims.n)*dim(rdata)[1])
rsims <- rdata[rsims.index,]

# 2) Generate independent random number(s) z
usims <- matrix(runif(2*sims.n),nrow=sims.n,ncol=2) 
# (upper/lower Fréchet only use the left coloumn of random variables)

# 3a) Upper Fréchet shuffle
Z.upper <- (rsims-usims[,c(1,1)])/param.m

# 3b) Lower Fréchet shuffle
Z.lower <- cbind(rsims[,1]-usims[,1],rsims[,2]+usims[,1]-1)/param.m #WTF?? #only dimension 2

# 3c) rook copula
Z.rook <- (rsims-usims)/param.m

# 3d) Gauss copula
Z.Gauss <- (rsims-1+rCopula(sims.n,normalCopula(0.8, dim=2)) )/param.m
  
# 4) for all three of them: Generate a pair (x,y) according to Lemma 1
getij <- function(Z, a, b, model=c("binom","nbinom","poisson"),mkplot=TRUE) {
  if (model == "binom") 
    d<-ceiling(sweep(Z,2,c(a,b),"*")) 
  else if (model == "nbinom") 
    d<-floor(sweep(Z/(1-Z),2,c(a,b),"*")) 
  else if (model == "poisson") 
    d<-floor(sweep(-log(1-Z),2,c((log(a+1)-log(a)),(log(b+1)-log(b))),"/")) 
  else
    d <- NULL
  if (mkplot)
    plot(d,xlim=c(0,40),ylim=c(0,40))
  return(d)
}

@

Some plots of the patchwork

<<algoresult, dev='cairo_pdf', echo=FALSE>>=

op <- par(mfrow = c(2, 2), pty = "s")  
plot(Z.upper,xlim=c(0,1),ylim=c(0,1),pch=21,cex=.1)
plot(Z.lower,xlim=c(0,1),ylim=c(0,1),pch=21,cex=.1)
plot(Z.rook,xlim=c(0,1),ylim=c(0,1),pch=21,cex=.1)
plot(Z.Gauss,xlim=c(0,1),ylim=c(0,1),pch=21,cex=.1)#,col="red",bg="red", type="p",lwd=1,pch=21)
par(op)
rm(op)

@

Some more plots from step 4. of algorithm 1 \eqref{algorithm1}: pairs $(i,j)$ for which the parameters $p_{ij}$ are positive

<<moreplots, dev='cairo_pdf', echo=FALSE, fig.lp="fig:", fig.cap = c('negative binomial copula, $a=18$, $b=22$','binomial copula, $a=22$, $b=27$','Poisson copula, $a=17$, $b=22$')>>=

#### for negative Binomial:
op <- par(mfrow = c(2, 2), pty = "s")  
nBin.ij.upper <- getij(Z.upper, par.a, par.b,model="nbinom") 
nBin.ij.lower <- getij(Z.lower, par.a, par.b,model="nbinom") 
nBin.ij.rook  <- getij(Z.rook, par.a, par.b,model="nbinom") 
nBin.ij.Gauss <- getij(Z.Gauss, par.a, par.b,model="nbinom") 
par(op)
rm(op)

#### for Binomial:
op <- par(mfrow = c(2, 2), pty = "s")  
Bin.ij.upper <- getij(Z.upper, bin.par.a, bin.par.b,model="binom") 
Bin.ij.lower <- getij(Z.lower, bin.par.a, bin.par.b,model="binom") 
Bin.ij.rook  <- getij(Z.rook,  bin.par.a, bin.par.b,model="binom") 
Bin.ij.Gauss <- getij(Z.Gauss, bin.par.a, bin.par.b,model="binom") 
par(op)
rm(op)
#### for Poisson:
op <- par(mfrow = c(2, 2), pty = "s")  
Pois.ij.upper <- getij(Z.upper, par.a, par.b,model="poisson") 
Pois.ij.lower <- getij(Z.lower, par.a, par.b,model="poisson") 
Pois.ij.rook  <- getij(Z.rook,  par.a, par.b,model="poisson") 
Pois.ij.Gauss <- getij(Z.Gauss, par.a, par.b,model="poisson") 
par(op)
rm(op)
@

Furthermore we get:

<<step5, dev='cairo_pdf'>>=

POIcopula <- function(n, ij, a, b, model=c("binom","nbinom","poisson")) {
  if (model == "binom") 
    return(cbind(qbeta(runif(n),ij[,1],rep(1,n)*a+1-ij[,1])
             , qbeta(runif(n),ij[,2],rep(1,n)*b+1-ij[,2])))
  else if (model == "nbinom") 
    return(cbind(qbeta(runif(n),ij[,1]+1,rep(1,n)*a+1)
             , qbeta(runif(n),ij[,2]+1,rep(1,n)*b+1)))
  else if (model == "poisson") 
    return(cbind(1-exp(-qgamma(runif(n),shape=ij[,1]+1,scale=1/(1+rep(1,n)*a)))
             , 1-exp(-qgamma(runif(n),shape=ij[,2]+1,scale=1/(1+rep(1,n)*b)))))
  
}

# negative binomial
nBin.xy.upper <- POIcopula(sims.n, nBin.ij.upper, par.a, par.b, model="nbinom")
nBin.xy.lower <- POIcopula(sims.n, nBin.ij.lower, par.a, par.b, model="nbinom")
nBin.xy.rook  <- POIcopula(sims.n, nBin.ij.rook , par.a, par.b, model="nbinom")
nBin.xy.Gauss <- POIcopula(sims.n, nBin.ij.Gauss , par.a, par.b, model="nbinom")
# binomial
Bin.xy.upper <- POIcopula(sims.n, Bin.ij.upper, bin.par.a, bin.par.b, model="binom")
Bin.xy.lower <- POIcopula(sims.n, Bin.ij.lower, bin.par.a, bin.par.b, model="binom")
Bin.xy.rook  <- POIcopula(sims.n, Bin.ij.rook , bin.par.a, bin.par.b, model="binom")
Bin.xy.Gauss <- POIcopula(sims.n, Bin.ij.Gauss , bin.par.a, bin.par.b, model="binom")
# Poisson
Pois.xy.upper <- POIcopula(sims.n, Pois.ij.upper, par.a, par.b, model="poisson")
Pois.xy.lower <- POIcopula(sims.n, Pois.ij.lower, par.a, par.b, model="poisson")
Pois.xy.rook  <- POIcopula(sims.n, Pois.ij.rook , par.a, par.b, model="poisson")
Pois.xy.Gauss <- POIcopula(sims.n, Pois.ij.Gauss , par.a, par.b, model="poisson")

@

The corresponding plots are:

<<step4plots, dev='cairo_pdf', echo=FALSE>>=
copulaAndPoints <- function(data, redpoints) {
  plot(data,xlim=c(0,1),ylim=c(0,1))
  points(redpoints,col="red",bg="red", type="p",lwd=2,pch=21)  
}
op <- par(mfrow = c(2, 2), pty = "s")  
  copulaAndPoints(nBin.xy.upper, udata)
  copulaAndPoints(nBin.xy.lower, udata)
  copulaAndPoints(nBin.xy.rook, udata)
  copulaAndPoints(nBin.xy.Gauss, udata)
par(op)
rm(op)
op <- par(mfrow = c(2, 2), pty = "s")  
  copulaAndPoints(Bin.xy.upper, udata)
  copulaAndPoints(Bin.xy.lower, udata)
  copulaAndPoints(Bin.xy.rook, udata)
  copulaAndPoints(Bin.xy.Gauss, udata)
par(op)
rm(op)
op <- par(mfrow = c(2, 2), pty = "s")  
  copulaAndPoints(Pois.xy.upper, udata)
  copulaAndPoints(Pois.xy.lower, udata)
  copulaAndPoints(Pois.xy.rook, udata)
  copulaAndPoints(Pois.xy.Gauss, udata)
par(op)
rm(op)

@

Finally output the empirical correslations

<<empcorrelations, dev='cairo_pdf'>>=

cor(nBin.xy.upper)
cor(nBin.xy.lower)
cor(nBin.xy.rook)
cor(nBin.xy.Gauss)

cor(Bin.xy.upper)
cor(Bin.xy.lower)
cor(Bin.xy.rook)
cor(Bin.xy.Gauss)

cor(Pois.xy.upper)
cor(Pois.xy.lower)
cor(Pois.xy.rook)
cor(Pois.xy.Gauss)

cor(Z.upper)
cor(Z.lower)
cor(Z.rook)
cor(Z.Gauss)

@


\bibliographystyle{plainnat} 
%Plain oder alpha
\bibliography{C://Users/maendle/Dropbox/PHD/umlauttest} %..\PHD\umlauttest}

\end{document}