Using copula to account for endogeneity in regression

Hello, all,
I have been working for last couple of weeks on using copula model to account for dependency between an independent variable and error. The model I used for a start was from Park and Gupta (2012).
I simulated a very simple endogenous regression model in R, but I am not able to recover the true parameter (coefficient of x1).
Can anybody help? @James_Savage @economicurtis are you able to help?

set.seed(123)

Generate independent variables x1
x1 ← rnorm(250)

Generate error term
error ← rnorm(250)

Generate endogenous variable z (related to x1 and error)
z ← 1.29 * x1 + error

Generate dependent variable y
beta0 ← 2
beta1 ← 0.5

y ← beta0 + beta1 * x1 + z

model = lm(y ~x1)
summary(model)

data ← list(N = n,x=x1, y = y,m = 2)

functions {
  vector normal_copula(vector u, vector v, real rho) {
    int N = size(u);
    vector[N] copula_values;

    for (i in 1:N) {
      copula_values[i] = (0.5 * rho * (-2. * u[i] * v[i] + u[i] * u[i] * rho + v[i] * v[i] * rho)) /
          (-1. + rho * rho) - 0.5 * log(1. - rho * rho);
    }

    return copula_values;
  }
}

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

parameters {
    real beta0;
  real beta;             // Coefficient for x
  real<lower=0> sigma_e; // Standard deviation of the error term
  real<lower=-1, upper=1> rho;   // Correlation parameter between x and the correlated error component
  vector[N] u;           // Uncorrelated error component
}

transformed parameters {
  vector[N] error;       // Error term
  vector[N] mu;          // Mean predicted value

  mu = beta0 + beta * x;         // Mean predicted value
  error = sigma_e * (rho * normal_copula(x, u, rho) + sqrt(1 - rho^2) * u);  // Error term using normal copula
}

model {
  // Priors
    beta0 ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma_e ~ cauchy(0, 5);
  rho ~ uniform(-1, 1);

  // Likelihood
  y ~ normal(mu + error, sigma_e);  // Likelihood of y

  // Priors for the uncorrelated error component
  u ~ normal(0, 1);  // Assuming standard normal distribution for u
}

****
@economicurtis