Is it possible to generate two correlated exponential variables with constraint on each other preserving their means? I wanted to generate two negatively correlated exponentially distributed variables (remission time $t_r$ and survival time $t$, for example) with the constraint that $t>t_r$. I wanted that the mean of $t_r$ be around $1$ and the mean of $t$ be around $2$. I tried the following codes:
rho   <- -0.5
mu    <- rep(0,2)
Sigma <- matrix(rho, nrow=2, ncol=2) + diag(2)*(1 - rho)

library(MASS)

compute.tr.t <- function(req.n, paccept) {
  req.n      <- round(req.n / paccept)
  rawvars    <- mvrnorm(req.n, mu=mu, Sigma=Sigma)
  pvars      <- pnorm(rawvars)
  tr         <- qexp(pvars[,1], 1/1)
  t          <- qexp(pvars[,2], 1/2)
  keep       <- which(t > tr)
  return(data.frame(t=t[keep],tr=tr[keep]))
}

req.n   <- n
paccept <- 1
res     <- data.frame()

while (req.n > 0) {
  new.res <- compute.tr.t(req.n, paccept)
  res     <- rbind(res, new.res)
  req.n   <- n - nrow(res)
  paccept <- nrow(new.res) / n# updated paccept according to last step
}

But here the means of both $t_r$ and $t$ changes and I also do not get the desired correlation. 
mean(res$tr)
[1] 0.4660927
mean(res$t)
[1] 2.859441
print(cor(res$tr,res$t))
[1] -0.237159

Even if I accept the reduced correlation, is it possible to generate two exponentially distributed random variables preserving their means and the constraint that $t>t_r$? 
I asked this in stackoverflow and I was suggested to ask this here. Is it that the distributions of $t_r$ and $t$ do not remain exponential anymore?   
 A: I offer a candidate solution which satisfies your requirements on the means. You will have to determine whether this satisfies your other desires.
Let $mean_{tr}$ = desired mean of tr, and $mean_t$ = desired mean of t.
Solution: 


*

*Draw all exponentials below independently

*Generate tr per an exponential random variable having mean of $mean_{tr}$

*For each sample value of tr, generate corresponding sample t = tr + [exponential random variable having mean of $mean_t - mean_{tr}]$


It might be easier to view the solution as follows: If tr is an n by 1 vector of exponential random numbers having mean of $mean_{tr}$, and $X$ is a vector of exponential random numbers having mean of $mean_t - mean_{tr}$, then $t = tr + X$ will be a vector of sample values of t.
tr is exponentially distributed. t is exponentially distributed conditional on tr, or equivalently is a shifted exponential, in which the shift is by the amount tr, which of course differs in numerical value for each random variate of tr.
Applying this solution to your stated values of $mean_{tr} = 1, mean_t = 2$, results in $corr(tr,t) = \sqrt{2}/2 = 0.707...$
Edit: In an attempt to satisfy all your requirements, including negative correlation between tr and t, I generated tr and X according to a Gaussian copula (rather than independently), but otherwise as per the solution above.  For tr and t having respective means of 1 and 2, the smallest correlation I was able to achieve was 0.42, using a correlation of essentially -1 in the Gaussian for the copula.
Alternatively, if you were willing to forgo the requirement that t > tr, then your desired correlation between tr and t of -0.5 could be achieved by using a Gaussian copula to generate tr and t (without the "t= tr + X" method in my solution above) with correlation of approximately -0.724 in the Gaussian for the copula.
I don't think it's possible to simultaneously satisfy all your requirements with the numerical values you have provided for mean(t), mean(tr), and corr(tr,t). I'd be very interested if someone can show that it is possible. Or perhaps someone can prove that it is not possible.
However, it is possible to meet all your requirements for other values.  For instance, if mean(tr) = 1 and mean(tr) > 5.9, the corr(tr,t) = -0.5 can be achieved using my "t = tr + X" method combined with the Gaussian copula, per the method described above, and that satisfies all your requirements (except that mean(t) must be > 5.9).  Similarly, if you only require corr(tr,t) be negative, but not that it equal -0.5, then if mean(tr) = 1, this approach will work if mean(t) > 2.6.  Perhaps use of this approach with some copula other than Gaussian would further "improve" things. I leave that as an exercise to you.
