如何模拟含指定shape、rate与phi的近似伽马分布AR(1)时间序列？

问题背景 #

需要生成满足以下条件的时间序列：

• 服从AR(1)过程，自相关系数为指定的phi
• 边际分布近似匹配指定shape和rate参数的伽马分布
• 接受通过调整参数来实现近似，无需严格的伽马AR(1)过程

核心思路：调整创新项的伽马参数 #

直接用arima.sim传入目标伽马参数作为创新项，得到的序列边际分布会偏离目标（因为AR(1)的边际矩由创新项矩和自相关系数共同决定）。我们可以通过匹配边际分布的均值和方差，反推出创新项的伽马参数，从而让生成的序列近似符合目标伽马分布。

推导过程

目标伽马分布的均值和方差：

• 均值：$\mu = \frac{\text{shape}}{\text{rate}}$
• 方差：$\sigma^2 = \frac{\text{shape}}{\text{rate}^2}$

对于平稳AR(1)过程 $X_t = \phi X_{t-1} + \varepsilon_t$（$\varepsilon_t$ 为伽马分布的创新项），其边际矩满足：

• $E[X_t] = \frac{E[\varepsilon_t]}{1-\phi}$
• $\text{Var}[X_t] = \frac{\text{Var}[\varepsilon_t]}{1-\phi^2}$

联立方程求解创新项的伽马参数（$\text{shape}\varepsilon$、$\text{rate}\varepsilon$）：

1. $\text{rate}_\varepsilon = \frac{\text{rate}}{1+\phi}$
2. $\text{shape}_\varepsilon = \text{shape} \times \frac{1-\phi}{1+\phi}$

实现代码 #

library(ggplot2)

# 目标参数
n <- 500
target_shape <- 1.5
target_rate <- 5
phi <- 0.3

# 计算调整后的创新项参数
rate_epsilon <- target_rate / (1 + phi)
shape_epsilon <- target_shape * (1 - phi) / (1 + phi)

# 生成目标伽马分布样本
target_gamma <- rgamma(n = n, shape = target_shape, rate = target_rate)

# 生成近似伽马分布的AR(1)序列
ar1_gamma <- arima.sim(
  list(order = c(1,0,0), ar = phi),
  n = n,
  rand.gen = rgamma,
  shape = shape_epsilon,
  rate = rate_epsilon
)

# 密度对比
ggplot() +
  geom_density(aes(x = target_gamma), fill = "#FF8C00", alpha = 0.5) +
  geom_density(aes(x = ar1_gamma), fill = "#A034F0", alpha = 0.5) +
  labs(title = "目标伽马分布与近似AR(1)序列密度对比", x = "数值", y = "密度") +
  theme_minimal()

说明 #

• 该方法通过匹配一二阶矩，让序列的边际分布近似目标伽马分布，自相关结构严格符合指定的AR(1)过程
• 若需要更精准的密度匹配，可以尝试通过模拟迭代微调参数，或使用Copula方法构建依赖结构，但上述方法已能满足大多数近似需求

内容的提问来源于stack exchange，提问作者user111024