Simulate an AR(1) process that approximates a gamma distribution

I’d like to simulate an AR(1) time series that approximates a gamma distribution rather than a normal distribution. I’d like the result to have a specified rate and shape along with a AR(1) process with a given phi. I have seen this post and understand that a strict gamma simulation is problematic, but I'd be satisfied with a result that approximates the shape and rate. Some crude adjustment of shape and rate as a function of phi is likely enough for my needs. Any idea?

E.g.,

n <- 500
foo <- rgamma(n = n, shape = 1.5 , rate = 5)

phi <- 0.3
bar <- arima.sim(list(order = c(1,0,0), ar = phi), n = n, 
                 rand.gen = rgamma, shape = 1.5,
                 rate = 5)

ggplot() + 
  geom_density(aes(x=foo),fill="#FF8C00",alpha=0.5) + 
  geom_density(aes(x=bar),fill="#A034F0",alpha=0.5)

Borrowing the function from this CV answer, we can first sample from the desired distribution, then reorder the samples to approximate the expected ACF from the AR process.

set.seed(721893540)

n <- 500
foo <- rgamma(n = n, shape = 1.5 , rate = 5)
phi <- 0.3

# match the autocorrelation for the first 4 lags, when the
# autocorrelation (in the limit as n -> Inf) first drops below 0.01
alpha <- phi^(0:ceiling(log(0.01)/log(0.3)))

# reorder foo to get the desired ACF (ARIMA(0,0,1), phi = 0.3)
bar <- acf.reorder(foo, alpha)

The ACF of bar will match alpha more closely than would be expected from a randomly generated AR(1):

# autocorrelation function of a random Gaussian AR(1):
ACF <- acf(arima.sim(list(order = c(1,0,0), ar = phi), n), plot = 0)$acf

data.frame(
  normal = ACF[1:length(alpha)],
  gamma = acf(bar, length(alpha) - 1, plot = 0)$acf
)
#>         normal      gamma
#> 1  1.000000000 1.00000000
#> 2  0.334034734 0.29929104
#> 3  0.060370711 0.09063563
#> 4  0.001948240 0.03818762
#> 5 -0.004354808 0.01196619

If this behavior is undesirable, one option is to set alpha to the ACF from a randomly generated AR(1) (until it stops decreasing monotonically).

alpha <- ACF[1:which.max(diff(ACF) > 0)]
bar2 <- acf.reorder(foo, alpha)

data.frame(
  normal = ACF[1:length(alpha)],
  gamma = acf(bar, length(alpha) - 1, plot = 0)$acf,
  gamma2 = acf(bar2, length(alpha) - 1, plot = 0)$acf
)
#>         normal      gamma       gamma2
#> 1  1.000000000 1.00000000  1.000000000
#> 2  0.334034734 0.29929104  0.333494498
#> 3  0.060370711 0.09063563  0.060757048
#> 4  0.001948240 0.03818762 -0.017663549
#> 5 -0.004354808 0.01196619  0.008826011

Including acf.reorder for convenience.

acf.reorder <- function(x, alpha) {
  tol <- 1e-5
  maxIter <- 10L
  n <- length(x)
  xx <- sort(x)
  y <- rnorm(n)
  w0 <- w <- alpha1 <- alpha
  m <- length(alpha)
  i1 <- sequence((m - 1):1)
  i2 <- sequence((m - 1):1, 2:m)
  i3 <- cumsum((m - 1):1)
  tol10 <- tol/10
  iter <- 0L
  x <- xx[rank(filter(y, w, circular = TRUE))]
  SSE0 <- Inf
  f <- function(ww) {
    sum((c(1, diff(c(0, cumsum(ww[i1]*(ww[i2]))[i3]))/sum(ww^2)) - alpha1)^2)
  }
  ACF <- function(x) acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m]
  
  while ((SSE <- sum((ACF(x) - alpha)^2)) > tol) {
    if (SSE < SSE0) {
      SSE0 <- SSE
      w <- w0
    }
    if ((iter <- iter + 1L) == maxIter) break
    w1 <- w0
    a <- 0
    sse0 <- Inf
    
    while (max(abs(alpha1 - a)) > tol10) {
      a <- c(1, diff(c(0, cumsum(w1[i1]*(w1[i2]))[i3]))/sum(w1^2))
      
      if ((sse <- sum((a - alpha1)^2)) < sse0) {
        sse0 <- sse
        w0 <- w1
      } else {
        # w0 failed to converge; try optim
        w1 <- optim(w0, f, method = "L-BFGS-B")$par
        a <- c(1, diff(c(0, cumsum(w1[i1]*(w1[i2]))[i3]))/sum(w1^2))
        if (sum((a - alpha1)^2) < sse0) w0 <- w1
        break
      }
      
      w1 <- (w1*alpha1/a + w1)/2
    }
    
    x <- xx[rank(filter(y, w0, circular = TRUE))]
    alpha1 <- (alpha1*alpha/ACF(x) + alpha1)/2
  }
  
  xx[rank(filter(y, w, circular = TRUE))]
}