如何自主实现AR(p)时间序列过程的模拟算法？

首先得给你点个赞——你已经摸到了AR(p)过程的两种核心视角：一种是基于Wold分解的无限白噪声线性组合，另一种是基于递推关系的有限初始值模拟。咱们来逐个拆解，帮你理清哪种思路更适合实际编程：

关于你的第一种思路：无限白噪声表示 #

你写的推导其实是理论上完全正确的——平稳AR(p)过程确实可以表示为无限阶MA过程（Wold定理），也就是所有过去白噪声的加权和：
$$ X_t = \sum_{k=0}^\infty \psi_k Z_{t-k} $$
其中$\psi_k$是脉冲响应系数，可以通过递推$\psi_0=1$，$\psi_k = \phi_1\psi_{k-1}+\dots+\phi_p\psi_{k-p}$（$k\geq1$，$\psi_k=0$当$k<0$）得到。

但实际模拟时这个思路不可行，因为你没法生成无限个过去的白噪声$Z_{-1}, Z_{-2}, ...$。除非你截断这个无限和（比如只取前M项，M足够大），但这本质上还是近似，不如直接用递推法高效准确。

关于你的第二种思路：有限初始值递推 #

这才是实际模拟AR(p)过程的标准路径！核心是先确定前p个初始值$X_1, X_2, ..., X_p$，然后用AR(p)的递推公式生成后续值。这里的关键是初始值的选择，不同选择会影响模拟序列的平稳性：

可选的初始值方案 #

• 方案1：使用平稳分布生成初始值
  平稳AR(p)过程的联合分布是多元正态（如果$Z_t$是正态白噪声），协方差矩阵$\Gamma$满足Yule-Walker方程：
  $$ \Gamma(k) = \phi_1\Gamma(k-1)+\dots+\phi_p\Gamma(k-p) \quad (k\geq1) $$
  其中$\Gamma(0)=\text{Var}(X_t)$，$\Gamma(k)=\text{Cov}(X_t,X_{t-k})$。你可以通过Cholesky分解$\Gamma$，从标准正态分布生成样本，再转换为符合平稳分布的初始值$X_1,...,X_p$。这种方法最准确，但需要计算协方差矩阵，适合对精度要求高的场景。

• 方案2：设置预热期（Burn-in）
  这是最常用的工程化方法：

  1. 随便设一组初始值（比如全0，或者直接用前p个白噪声$Z_1,...,Z_p$）；
  2. 先模拟远多于你需要的序列长度（比如需要N个样本，先模拟N+1000个）；
  3. 丢弃前1000个“预热”样本，只保留后面的N个。
     这么做的原因是，随着递推次数增加，初始值的影响会被逐渐稀释，后面的样本会趋近于平稳分布的样本。这种方法简单易实现，效果也足够好。

• 方案3：简单初始值（仅用于快速测试）
  如果只是做快速验证，可以直接设$X_1=X_2=...=X_p=0$，或者用$X_k=Z_k$（$k=1..p$）。但要注意，前p个样本可能不符合平稳分布，序列开头会有“瞬态”，适合对精度要求不高的场景。

具体算法步骤（以预热期方案为例） #

假设你要模拟长度为N的AR(p)序列，步骤如下：

1. 验证平稳性：先确认AR(p)的特征方程$1-\phi_1z-\dots-\phi_pz^p=0$的所有根都在单位圆外（这是平稳的充要条件）；
2. 生成白噪声：生成$N+M$个独立同分布的白噪声样本$Z_1,Z_2,...,Z_{N+M}$（M是预热长度，比如取1000）；
3. 初始化序列：设$X_1=X_2=...=X_p=0$（或其他初始值）；
4. 递推生成序列：对$t=p+1$到$t=N+M$，计算：
   $$ X_t = \phi_1X_{t-1} + \phi_2X_{t-2} + \dots + \phi_pX_{t-p} + Z_t $$
5. 截取有效样本：取$X_{M+1}, X_{M+2}, ..., X_{M+N}$作为最终的模拟序列。

代码示例（Python） #

比如模拟一个AR(2)过程，参数$\phi_1=0.5, \phi_2=-0.2$，白噪声是标准正态分布：

import numpy as np

def simulate_ar(p, phi, n_samples, burn_in=1000, noise_mean=0, noise_var=1):
    # 验证平稳性：检查特征根是否都在单位圆外
    roots = np.roots(np.concatenate([[1], -np.array(phi)]))
    if any(np.abs(roots) >= 1):
        raise ValueError("AR(p) process is not stationary!")
    
    # 生成白噪声序列
    total_length = n_samples + burn_in
    z = np.random.normal(noise_mean, np.sqrt(noise_var), total_length)
    
    # 初始化AR序列
    x = np.zeros(total_length)
    # 递推生成后续值
    for t in range(p, total_length):
        x[t] = np.dot(phi, x[t-p:t]) + z[t]
    
    # 丢弃预热期，返回有效样本
    return x[burn_in:]

# 测试：模拟AR(2)序列
ar_series = simulate_ar(p=2, phi=[0.5, -0.2], n_samples=1000)

总结 #

• 你的第一种思路是理论上的Wold分解，适合理解AR(p)的本质，但不适合实际模拟；
• 第二种思路的有限初始值递推是实际模拟的标准方法，搭配预热期方案既简单又能保证序列的平稳性。

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