n <- 100000
x <- matrix(NA, nrow <- 3, ncol=n)
x[,1] <- c(1, 1, 1)

for(i in 2:n) {
x[1,i]<- rexp(1, rate = 1+x[2,i-1] * x[3,i-1])
x[2,i]<- rexp(1, rate = 2+x[1,i] * x[3,i-1])
x[3,i]<- rexp(1, rate = 3+x[1,i] * x[2,i])
}