Hi

This is a question about copulas in Stan. I am trying to estimate a copula model which introduces some correlation between a Gaussian and a Poisson regression. This is my data generating process:

```
library(MASS)
s_size <- 600
x1 <- rnorm(s_size, mean = 1, sd = 3)
x2 <- sample(c(0, 1), s_size, replace = TRUE)
# correlation parameter
rho_pop <- .4
cov_matrix <- matrix(c(1, rho_pop, rho_pop, 1), ncol = 2)
cr_error <- mvrnorm(s_size, mu = c(0, 0), Sigma = cov_matrix)
# transform gaussian error to uniform
x_transf <- pnorm(cr_error)
# transform outcomes back to gaussian/Poisson
y_gauss <- qnorm(x_transf[,2], mean = (1 + .5 * x1 +.05 * x2), sd = .8)
y_pois <- qpois(x_transf[,1], lambda = exp(-1 + 0.3 * x1 - 0.5 * x2))
```

The Stan code I use is provided below. The copula code works fine when using simple outcomes (without a regression model). However, if I try to use it for the regression setup, I consistently estimate the wrong copula correlation parameter. The issue seems to arise when I introduce the transformed parameter step to transform the Gaussian and Poisson outcomes to uniform.

```
functions {
real normal_copula_llpdf(real u, real v, real rho) {
real rho_sq = square(rho);
return (0.5 * rho * (-2. * u * v + square(u) * rho + square(v) * rho)) / (-1. + rho_sq)
- 0.5 * log1m(rho_sq);
}
}
data {
int<lower=1> s_size; // number of responses
int<lower=0> K1; // Number of independent variables
int<lower=0> K2; // Number of independent variables
matrix[s_size, K1] mX1; // independent variables
matrix[s_size, K2] mX2; // independent variables
vector[s_size] y_gauss;
int y_pois[s_size];
}
parameters {
real mu_gauss;
vector[K1] betas_gauss;
real<lower=0> sigma_gauss;
real mu_pois;
vector[K2] betas_pois;
real<lower=-1,upper= 1> rho;
}
transformed parameters {
// transform gaussian outcomes to uniform by taking normal cdf
real<lower=0,upper=1> y_gauss_cop[s_size];
real<lower=0,upper=1> y_pois_cop[s_size];
for (i in 1:s_size){
y_gauss_cop[i] = normal_cdf(y_gauss[i], mu_gauss+mX1[i,]*betas_gauss,sigma_gauss);
y_pois_cop[i] = poisson_cdf(y_pois[i], exp(mu_pois + mX2[i,]*betas_pois));
}}
model {
rho ~ normal(0,1);
y_gauss ~ normal(mu_gauss + mX1*betas_gauss, sigma_gauss);
y_pois ~ poisson(exp(mu_pois + mX2*betas_pois));
// copula pdf
for (i in 1:s_size){
target += normal_copula_llpdf(y_gauss_cop[i], y_pois_cop[i], rho);
}
}
```

I should also mention Iâ€™m not 100 % sure if the mistake concerns the stan code or if it is a conceptual mistake (so perhaps Iâ€™m just formulating the wrong model). Anybody would have any idea about how to tackle this?