Frank Copula

[ants007]

Hello all , I am trying to work on Frank Copula. Here is my artifical data set that generates correlation between x and y in R using frank copula

library(VineCopula)

Set the copula parameter
theta ← -10

Set the marginals (marginal distributions for x and y)
marginals ← c(“norm”, “norm”)

Create the Frank copula object
copula ← frankCopula(param = theta)

Generate random variables from the copula
n ← 500 # Number of data points
u ← rCopula(n, copula)

Apply inverse transformation to get simulated values for x and y
x ← qnorm(u[, 1])
y ← qnorm(u[, 2])

Plot the simulated data
plot(x, y, pch = 16, xlab = “x”, ylab = “y”, main = “Simulated Data”)

Compute and print the sample correlation
correlation ← cor(x, y)
print(paste(“Sample Correlation:”, correlation))

I am using the following Stan Code to recover the value of theta which is -10, but I am having no success, can some one please help
@spinkney I know you have done some work on this can you help please?
@andre.pfeuffer can you look into this too, please?

---

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

functions {
real frank_copula(real u, real v, real theta) {
  return log(theta) + log1m(exp(-theta)) - theta * (u + v)
         - 2 * log(1 - exp(-theta) - (1 - exp(-theta * u)) * (1 - exp(-theta * v)));
}

}

data {
  int<lower=0> N;  // Number of observations
  vector[N] x;     // First variable
  vector[N] y;     // Second variable
}

transformed data {
  vector[N] x_scaled = (x - min(x)) / (max(x) - min(x));
  vector[N] y_scaled = (y - min(y)) / (max(y) - min(y));
}

parameters {
  real theta;  // Copula parameter
}

model {
  // Priors
  // No constraint on theta

  // Generate the copula
  vector[N] z;
  for (i in 1:N) {
    z[i] = frank_copula(x_scaled[i], y_scaled[i], theta);
  }

  // Likelihood
  for (i in 1:N) {
    target += z[i];
  }
}

generated quantities {

}
****

[spinkney]

I believe it needs to handle the case when \theta < 0. You can easily vectorize this ldpf too (see below). I don’t think you want to do the scaling either. Just input u directly.

mod_out <- mod$sample(
  data = list(N = n,
              x = u[, 1],
              y = u[, 2]),
  parallel_chains = 4
)

functions {
real frank_copula(real u, real v, real theta) {
  return log(theta * (1 - exp(-theta))) - theta * (u + v) -
  log( (1 - exp(-theta) - (1 - exp(-theta * u)) * (1 - exp(-theta * v)))^2);
}

real frank_copula_lpdf(vector u, vector v, real theta) {
  int N = num_elements(u);
  return N * log(theta * (1 - exp(-theta))) - theta * sum(u + v) -
  sum(log( (1 - exp(-theta) - (1 - exp(-theta * u)) .* (1 - exp(-theta * v)))^2));
}

}
data {
  int<lower=0> N;  // Number of observations
  vector[N] x;     // First variable
  vector[N] y;     // Second variable
}
parameters {
  real theta;  // Copula parameter
}
model {
  // Priors
  // No constraint on theta
   target += frank_copula_lpdf(x, y, theta);
}