I am trying to find conditional probability of the form P(X<x|Y=y) for two jointly distributed random variables based on the copula estimate from training data. I use R package copula but can not figure out the best way to do it.

What I do now - estimate empirical copula parameters on training data, generate 100000 outcomes from this distribution, construct rank-transformed data for testing data based on training data, find number of cases X<x within Y=y+/-eps for each outcome in testing data based on generated distribution. The code for doing thing is below.

Could you please advice whether there is better way of finding P(X<x|Y=y) for testing data based on the training data?


    t.cop0 <- tCopula(0.5,dim=2,dispstr='un',df=1.7)
    gendata <- rCopula(300,t.cop0)
    train <- gendata[1:199,]
    test <- gendata[200:300,]

    ptrain <- pobs(train)
    tau <- cor(train,method='kendall')[2]
    t.cop <- tCopula(tau,dim=2,dispstr='un',df=3)
    fit.mpl <- fitCopula(t.cop,ptrain,method='mpl',estimate.variance=FALSE)
    empiricalCopula <- tCopula(fit.mpl@estimate[1],dim=2,dispstr='un',df=fit.mpl@estimate[2])

    p1 <- sapply(as.numeric(test[,1]),function(q)rank(c(q,train[,1]))[1]/nrow(train+2))
    p2 <- sapply(as.numeric(test[,2]),function(q)rank(c(q,train[,2]))[1]/nrow(train+2))
    ptest <- cbind(p1,p2)

    e <- rCopula(100000,empiricalCopula)
    eps <- .1
    cp <- sapply(1:nrow(ptest),function(i)
            sum(e[,2]<=ptest[i,2] & e[,1]>=(ptest[i,1]-eps) & e[,1]<=(ptest[i,1]+eps))/
            sum(e[,1]>=(ptest[i,1]-eps) & e[,1]<=(ptest[i,1])+eps))
  • $\begingroup$ And how to find, for example $P(X>a|Y>b)$ in terms of spcopula package, without calculating double integrals of copula density? $\endgroup$ – user64186 Dec 20 '14 at 0:24
  • $\begingroup$ the df=1.7 was used in the second line. The df is the degree of freedom It must be integer. $\endgroup$ – Nick Apr 8 '17 at 22:58

Given a bivariate copula $C(u,v)$, a bivariate CDF $H(x,y)$ and marginal CDFs $F(x)$ and $G(y)$ with $H(x,y)=C(F(x),G(y))$. Then, what you are targeting at is (with some abuse of notation):

$$P(X < x | Y=y) = \frac{P(X < x, Y=y)}{P(Y=y)} = \frac{\frac{\partial}{\partial v} C(F(x), G(y)) \cdot g(y)}{g(y)} = \frac{\partial}{\partial v} C(F(x), G(y))$$

The partial derivatives of copulas are available from the spcopula package on r-forge. Your R call for $P(X < 0.2 | Y=0.3)$ assuming for instance a Gumbel copula with parameter 5 and standard normal margins would then look like:

ddvCopula(c(pnorm(0.2), pnorm(0.3)), gumbelCopula(5))
[1] 0.3834244
  • $\begingroup$ Thank you so much Ben! This is exactly what I was looking for. $\endgroup$ – Alex Mar 21 '14 at 15:30
  • $\begingroup$ Do you mean $P(X \leq x | Y=y)$? $\endgroup$ – Kiran K. Oct 17 '15 at 12:58
  • 1
    $\begingroup$ @Ben, how to install the your spcopula package? $\endgroup$ – Nick Apr 9 '17 at 14:28

in R package VineCopula:

# for Gumbel copula with parameter =5
A = BiCop(family=4, par=5, tau = NULL, check.pars = TRUE)
u1 = 0.03
u2 = 0.3
B = BiCopHfunc1(u1, u2, A)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.