ARIMA与SARIMA状态空间形式转换编程求助：差分扩展问题

嗨，我正好在时间序列建模中经常处理这类状态空间转换的问题，咱们一步步拆解来解决你的困惑！首先你提到的Koopman著作第54页的ARMA(p,q)状态空间形式，咱们先快速回顾，再扩展到差分、季节差分的场景。

1. 先锚定ARMA(p,q)的基础状态空间（对应Koopman第54页） #

Koopman书中的ARMA(p,q)通常用创新形式（Innovation Form），把所有滞后的观测值和噪声项都纳入状态向量，方便后续扩展。假设我们的ARMA(p,q)模型是：
$$
y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + ... + \phi_p y_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q \epsilon_{t-q}
$$
其中$\epsilon_t \sim N(0, \sigma^2)$是白噪声。

对应的状态向量一般是$s_t = [y_t, y_{t-1}, ..., y_{t-p+1}, \epsilon_t, \epsilon_{t-1}, ..., \epsilon_{t-q+1}]^T$，维度为$p+q$。此时状态空间的两个核心方程为：

• 状态方程：$s_t = F s_{t-1} + G \epsilon_t$
  转移矩阵$F$的结构很直观：
  • 第一行：前p个元素是AR系数$\phi_1$到$\phi_p$，后q个元素是MA系数$\theta_1$到$\theta_q$
  • 第2到p行：是“下移矩阵”，比如第2行只有第2列（对应$y_{t-1}$）为1，其余为0，用来把上一时刻的$y_{t-1}$移到当前状态的$y_{t-2}$位置
  • 第p+1到p+q行：同样是下移矩阵，第p+2行只有第p+1列（对应$\epsilon_{t-1}$）为1，用来把上一时刻的$\epsilon_{t-1}$移到当前状态的$\epsilon_{t-2}$位置
• 观测方程：$y_t = H s_t$，其中$H$是行向量，第一个元素为1，其余为0（直接提取状态向量里的$y_t$）

2. 扩展到ARIMA(p,d,q)：嵌入差分的逆过程 #

ARIMA的核心是对原始序列$y_t$做d阶差分得到平稳序列$z_t = (1-L)^d y_t$，而$z_t$服从ARMA(p,q)。我们不需要单独建模$z_t$再转换，而是把**差分的逆过程（积分）**直接嵌入状态向量。

举个具体例子：ARIMA(p,1,q)（d=1） #

一阶差分$z_t = y_t - y_{t-1}$，逆过程就是$y_t = y_{t-1} + z_t$。$z_t$的ARMA(p,q)模型为：
$$
z_t = \phi_1 z_{t-1} + ... + \phi_p z_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q \epsilon_{t-q}
$$

此时构造状态向量需要包含：

• 原始观测的滞后项：$y_t$
• 差分后序列的滞后项：$z_t, z_{t-1}, ..., z_{t-p+1}$
• 噪声的滞后项：$\epsilon_t, \epsilon_{t-1}, ..., \epsilon_{t-q+1}$

状态向量$s_t = [y_t, z_t, z_{t-1}, ..., z_{t-p+1}, \epsilon_t, ..., \epsilon_{t-q+1}]^T$，对应的状态方程：

• 第一行（$y_{t+1}$）：$y_{t+1} = y_t + z_{t+1}$，而$z_{t+1}$用ARMA方程替换，得到$y_{t+1} = y_t + \phi_1 z_t + ... + \phi_p z_{t-p+1} + \epsilon_{t+1} + \theta_1 \epsilon_t + ... + \theta_q \epsilon_{t-q+1}$
• 第二行（$z_{t+1}$）：就是上面的$z_{t+1}$的ARMA表达式
• 第3到p+1行：下移$z$的滞后项（比如第3行只有第2列为1，把$z_t$移到$z_{t-1}$的位置）
• 第p+2到p+q+1行：下移$\epsilon$的滞后项

观测方程依然是$y_t = H s_t$，$H$的第一个元素为1。

对于d>1的情况 #

比如d=2，二阶差分$z_t = y_t - 2y_{t-1} + y_{t-2}$，逆过程是$y_t = 2y_{t-1} - y_{t-2} + z_t$。此时状态向量需要加入$y_{t-1}$，状态方程的第一行就变成$y_{t+1} = 2y_t - y_{t-1} + z_{t+1}$，以此类推，把d阶积分的递归关系嵌入状态方程即可。

3. 扩展到SARIMA(p,d,q)(P,D,Q)_s：加入季节项 #

SARIMA的复杂度在于同时有普通差分、季节差分，以及普通ARMA和季节ARMA项。我们可以把它拆成两步处理：

1. 对$y_t$做d阶普通差分+D阶季节差分，得到平稳序列$z_t = (1-L)^d (1-L^s)D y_t$
2. $z_t$服从乘积型ARMA模型：$(1-\phi_1 L - ... - \phi_p L^p)(1-\Phi_1 L^s - ... - \Phi_P L^{Ps}) z_t = (1+\theta_1 L + ... + \theta_q L^q)(1+\Theta_1 L^s + ... + \Theta_Q L^{Qs}) \epsilon_t$

此时状态向量需要额外加入：

• 季节差分的逆过程项：比如$y_{t-s}, y_{t-2s}, ...$（对应D阶季节积分）
• 季节AR(P)的滞后项：$z_{t-s}, z_{t-2s}, ..., z_{t-Ps}$
• 季节MA(Q)的滞后噪声项：$\epsilon_{t-s}, ..., \epsilon_{t-Qs}$

示例：SARIMA(1,1,1)(1,1,1)_12（月度数据） #

• 普通差分d=1：$w_t = y_t - y_{t-1}$
• 季节差分D=1：$z_t = w_t - w_{t-12}$
• z_t的模型：$(1-\phi L)(1-\Phi L^{12}) z_t = (1+\theta L)(1+\Theta L^{12}) \epsilon_t$

状态向量可以设为$s_t = [y_t, y_{t-1}, y_{t-12}, z_t, z_{t-1}, z_{t-12}, \epsilon_t, \epsilon_{t-1}, \epsilon_{t-12}]^T$，然后根据递归关系写出状态方程：

• $y_{t+1} = y_t + w_{t+1} = y_t + z_{t+1} + w_{t-11}$（而$w_{t-11} = y_{t-11} - y_{t-12}$，如果需要可以把$y_{t-11}$也加入状态向量）
• $z_{t+1} = \phi z_t + \Phi z_{t-11} + \epsilon_{t+1} + \theta \epsilon_t + \Theta \epsilon_{t-11}$

观测方程还是直接提取$y_t$即可。

4. 实用技巧：简化推导的小方法 #

如果你觉得手动推导矩阵太繁琐，可以用以下技巧：

• 当q≤p时，MA项可以被吸收到状态方程的噪声项中，用可控形式的状态空间，只需要包含观测值的滞后项，不需要单独的噪声滞后项
• 用“扩展状态向量”的思路：不管是差分还是季节项，只要把需要的滞后变量都放进状态向量，再根据递归关系写出转移矩阵F，就不会出错

最后给个极简示例：ARIMA(1,1,1)的状态空间 #

让你更直观地看到转换后的结果：
ARIMA(1,1,1)的模型展开为：
$$y_t = y_{t-1} + \phi(y_{t-1} - y_{t-2}) + \epsilon_t + \theta \epsilon_{t-1}$$

状态向量设为$s_t = [y_t, y_{t-1}, \epsilon_t]^T$，对应的状态方程：
$$
\begin{bmatrix} y_{t+1} \ y_t \ \epsilon_{t+1} \end{bmatrix} =
\begin{bmatrix} 1+\phi & -\phi & \theta \ 1 & 0 & 0 \ 0 & 0 & 0 \end{bmatrix}
\begin{bmatrix} y_t \ y_{t-1} \ \epsilon_t \end{bmatrix} +
\begin{bmatrix} 1 \ 0 \ 1 \end{bmatrix} \epsilon_{t+1}
$$

观测方程：
$$
y_t = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} y_t \ y_{t-1} \ \epsilon_t \end{bmatrix}
$$

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