ARMA模型的MLE估计

暂时用来记录，日后完整。

MLE

$$ \begin{equation} \begin{aligned} & \underset{\boldsymbol{a},\boldsymbol{b},\theta}{\mathrm{minimize}} & & \sum_{t=1}^n \rho_\theta(\epsilon_t)\\ & \mathrm{subject\quad to}&& {\epsilon}_t = {y}_t - a_0 - \sum_{i=1}^p a_i{y}_{t-i} - \sum_{j=1}^q b_j {\epsilon}_{t-q},\text{ for }t=1,\ldots,n,\\ & && \text{moduli of the roots are greater than one}. \end{aligned} \end{equation} $$

ARMA的精确MLE一般可以通过两种方式得到：

状态空间演化，Kalman 滤波的思想，可见 Maximum Likelihood Fitting of ARMA Models to Time Series with Missing Observations
新息迭代法，本质是协方差矩阵的 Cholesky 分解，重要的是新息迭代法本身就提供了一种求解精确MLE的迭代形式

精确似然所带来的优化问题是非常难求解的：

新息迭代，那么算法的复杂度会非常高
直接使用线性规划方法，例如内点法，需要对新息的非线性约束进行处理，导致优化的问题严重非凸，只能找到局部最优解
ARMA 稳定性和可逆性的要求是非常难满足的，要求AR和MA部分的根在单位圆之外，对于比较高阶的模型，没有闭式的求根公式

对于稳定性和可逆性的约束 $S_p\times S_q$ ： $$ \begin{aligned} & S_p=\left\{\phi \in \mathbb{R}^p \mid 1-\phi_1 z-\ldots-\phi_p z^p \neq 0 \forall z \in \mathbb{C} \text { s.t. }|z| \leq 1\right\} \\ & S_q=\left\{\theta \in \mathbb{R}^q \mid 1+\theta_1 z+\ldots+\theta_p z^q \neq 0 \forall z \in \mathbb{C} \text { s.t. }|z| \leq 1\right\} . \end{aligned} $$ 一般用Jones reparametrization进行处理，可以通过引入变量的形式将该约束变为无约束。

除了精确MLE，还有条件MLE，即仅仅使用ARMA的模型，而忽略掉前 $m(m = \max(p,q))$ 个数据的分布,在数据量比较少时偏差很大，而且一般也需要数值求解。

Jones reparametrization

Jones reparametrization 旨在处理平稳型和可逆性的非线性约束，而偏自相关函数与模型的可逆性和稳定性有着天然的关系：

A note on enforcing stationarity in autoregressive-moving average models：

Theorem. We have that $y^{(k)} \in C_k$ if and only if $\left|y_k^{(k)}\right|<1$ and $y^{(k-1)} \in C_{k-1}$, where $$ y_i^{(k-1)}=\left(y_i^{(k)}-y_k^{(k)} y_{k-i}^{(k)}\right) /\left\{1-\left(y_k^{(k)}\right)^2\right\} \quad(i=1, \ldots, k-1) $$ For an autoregressive process of order $n, y^{(k)}$ solves the Yule-Walker equations of order $k$ for $k \leqslant n$. Thus for testing $y^{(n)} \in C_n$, the $n$ conditions $\left|y_k^{(k)}\right|<1$ are obvious, since $r_k \equiv y_k^{(k)}$ are partial autocorrelations.

凑巧的是 reverse Levinson-Durbin 提供了从 PACF 推出 AR 系数的方法（双射），本质还是使用Levinson-Durbin进行迭代，这样，我们就可以将对复杂的非线性约束转化到对PACF的约束上（$[-1,1]$）,从而迭代的优化PACF，最后才用PACF表达出AR系数。

Jones reparametrization进一步的对PACF进行了处理，将其变为无约束优化： $$ a_k=\frac{1-\exp \left(-u_k\right)}{1+\exp \left(-u_k\right)}, \quad k=1, \ldots, p $$

同样，对MA部分也可以使用类似的 Jones reparametrization。

可以发现，尽管使用Jones reparametrization可以解决非线性约束的问题，但是我们使用 $$u_i\rightarrow PACF\rightarrow a_i \rightarrow L(\theta)$$ 的方法依然会将所有的参数杂糅在一起，而且每一步优化的迭代非常复杂,使用Jones reparametrization并没有解决目标函数的非凸，也并没有降低计算复杂度，只是解决了一个非线性约束。

新息迭代方法

使用完整的历史信息构建似然和迭代

迭代：

Wold 表示 $$ \gamma_k=\sigma^2 \sum_{j=0}^{\infty} \psi_j \psi_{j+k}, \quad k=0,1,2, \ldots $$ $$ \psi_j= \begin{cases}1, & j=0, \\ b_j+\sum_{k=1}^p a_k \psi_{j-k}, & j=1,2, \ldots\end{cases} $$ Toeplitz 协方差矩阵建模 $$ Y_t= \begin{cases}X_t / \sigma, & t=1,2, \ldots, m \\ A(\mathscr{B}) X_t / \sigma, & t=m+1, \ldots\end{cases} $$ $$ E\left(Y_s Y_t\right)= \begin{cases}\sigma^{-2} \gamma_{t-s}, & 1 \leq s \leq t \leq m, \\ \sigma^{-2}\left[\gamma_{t-s}-\sum_{j=1}^p a_j \gamma_{t-s-j}\right], & 1 \leq s \leq m<t, \ \sum_{j=0}^q b_j b_{j+t-s}, & t \geq s>m,\end{cases} $$ Innovations Algorithm，Durbin-Levinson Algorithm，递归实现 Cholesky 分解的算法 $$ \left\{\begin{array}{l} \nu_0=E Y_1^2 \\ \theta_{n, n-k}=\left[E\left(Y_{n+1} Y_{k+1}\right)-\sum_{j=0}^{k-1} \theta_{k, k-j} \theta_{n, n-j} \nu_j\right] / \nu_k \\ \quad 0 \leq k \leq n-1 \\ \nu_n=E Y_{n+1}^2-\sum_{k=0}^{n-1} \theta_{n, n-k}^2 \nu_k \end{array}\right. $$ Cholesky 分解下的精确高斯似然：

利用Cholesky 分解思想引入残差构造似然 $$ L\left(\beta, \sigma^2\right)=\frac{\exp \left(-\frac{1}{2 \sigma^2} \sum_{k=1}^N Z_k^2 / \nu_{k-1}\right)}{(2 \pi)^{N / 2}\left(\sigma^{2 N} \nu_0 \nu_1 \cdots \nu_{N-1}\right)^{1 / 2}} $$ 正交投影残差 $$ S(\beta)=\sum_{k=1}^N Z_k^2 / \nu_{k-1} $$ $$ l(\beta) =\frac{1}{N} \ln \left(\nu_0 \nu_1 \cdots \nu_{N-1}\right)+\ln [S(\beta) / N] $$ $$ \hat{\sigma}^2=\frac{1}{N} S(\hat{\beta}) $$

解释新息迭代算法为Cholesky分解

用 Innovations Algorithm 实现 Cholesky 分解 的过程。我们以 一个长度为 3 的高斯向量 $Y = (Y_1, Y_2, Y_3)^T$ 为例，协方差矩阵为 $\Gamma_3$。

假设的协方差矩阵

我们假设 $Y$ 是平稳的高斯序列，有以下协方差结构： $$ \Gamma = \begin{bmatrix} \gamma_0 & \gamma_1 & \gamma_2 \\ \gamma_1 & \gamma_0 & \gamma_1 \\ \gamma_2 & \gamma_1 & \gamma_0 \end{bmatrix} $$ 例如取数值： $$ \gamma_0 = 1,\quad \gamma_1 = 0.5,\quad \gamma_2 = 0.25 $$ 于是： $$ \Gamma = \begin{bmatrix} 1 & 0.5 & 0.25 \\ 0.5 & 1 & 0.5 \\ 0.25 & 0.5 & 1 \end{bmatrix} $$ 我们要构造： $$ \Gamma = C D C^T $$

其中：

$C$：单位对角下三角矩阵，包含预测系数 $\theta_{i,j}$
$D = \mathrm{diag}(\nu_0, \nu_1, \nu_2)$：预测误差方差（新息方差）

然后我们就能得到 Cholesky 因子 $L = C D^{1/2}$，满足 $\Gamma = L L^T$

第一步：计算 $\nu_0 = \mathbb{E}[Y_1^2] = \gamma_0 = 1$

第二步：计算 $\theta_{1,1}$ 和 $\nu_1$ $$ \theta_{1,1} = \frac{\mathbb{E}[Y_2 Y_1]}{\nu_0} = \frac{\gamma_1}{\nu_0} = \frac{0.5}{1} = 0.5 $$ $$ \nu_1 = \mathbb{E}[Y_2^2] - \theta_{1,1}^2 \nu_0 = \gamma_0 - (0.5)^2 \cdot 1 = 1 - 0.25 = 0.75 $$ 第三步：计算 $\theta_{2,2}$、$\theta_{2,1}$，以及 $\nu_2$

按递推公式先计算： $$ \theta_{2,1} = \frac{\mathbb{E}[Y_3 Y_1]}{\nu_0} = \frac{\gamma_2}{\nu_0} = \frac{0.25}{1} = 0.25 $$ 再计算： $$ \theta_{2,2} = \frac{\gamma_1 - \theta_{1,1} \cdot \theta_{2,1} \cdot \nu_0}{\nu_1} = \frac{0.5 - (0.5)(0.25)(1)}{0.75} = \frac{0.5 - 0.125}{0.75} = \frac{0.375}{0.75} = 0.5 $$ 然后计算： $$ \nu_2 = \gamma_0 - \theta_{2,1}^2 \nu_0 - \theta_{2,2}^2 \nu_1 = 1 - 0.25^2 \cdot 1 - 0.5^2 \cdot 0.75 = 1 - 0.0625 - 0.1875 = 0.75 $$

得到 $C$ 和 $D$: $$ C = \begin{bmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.25 & 0.5 & 1 \end{bmatrix}, \quad D = \mathrm{diag}(1,\ 0.75,\ 0.75) $$

可见，$C$ 正是我们在高斯假设中得到的对新息线性变换的矩阵。

得到 Cholesky 因子 $L = C D^{1/2}$:

$$ D^{1/2} = \mathrm{diag}(1,\ \sqrt{0.75},\ \sqrt{0.75}) \approx \mathrm{diag}(1,\ 0.866,\ 0.866) $$

逐行乘：

$$ L=C D^{1 / 2}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.25 & 0.5 & 1 \end{array}\right]\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 0.866 & 0 \\ 0 & 0 & 0.866 \end{array}\right]=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0.5 & 0.866 & 0 \\ 0.25 & 0.433 & 0.866 \end{array}\right] $$

最终结果：验证 $$ \Gamma = L L^T \Rightarrow \text{构造出的 } L \text{ 正是 } \Gamma \text{ 的 Cholesky 分解} $$

现有的arma参数求解器

Matlab

EstMdl = estimate(Mdl,y)

使用Box书中的方法根据请款产生初始参数估计，作为迭代起点
确保初始迭代解满足约束要求
设置负对数似然目标函数
根据稳定性、可逆性、平稳性等设置约束条件arimaNonLineqrConstraints
调用fmincon进行优化
检查优化情况，改进数值问题匹配模型
计算后续结果，输出模型

构造了专用的负对数似然计算函数，内部 internal.econ.arimaxMex 使用C语言快速迭代的计算残差，上述优化问题中关于残差的第一条约束已经隐藏进似然函数的计算中。

非线性约束arimaNonLineqrConstraints就是直白的特征根在单位圆的外部，然后使用内点法障碍函数处理该非线性约束，所以最后该约束未必是满足的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


初始化：
  Estimate parameters of an ARIMAX model by one of the following approaches:

  (1) If the model has a regression component OR no MA component, then 
      estimate the constant, AR, and predictor coefficients by OLS regression. 

     After the ARX component is estimated, a subsequent estimation of any 
     MA coefficients is made by applying the auto-covariance approach outlined 
      in Box, Jenkins, and Reinsel (Appendix A6.2, pages 220-21), to the 
     filtered residuals obtained in the first step, modified to account 
      for any user-specified coefficient values.

  (2) If the model has no regression component AND an MA component, then 
      estimate the AR and MA coefficients by the auto-covariance method 
     outlined in Box, Jenkins, and Reinsel (Appendix A6.2, pages 220-21), 
    modified to account for any user-specified coefficient values.

目标函数：
F  = @(X) nLogLike(X, YData, XData, E, V, Mdl, AR.Lags, MA.Lags, maxPQ, T, isDistributionT, options, ...
                   userSpecifiedY0, userSpecifiedE0, userSpecifiedV0, trapValue);

优化：
 [coefficients, nLogL, exitflag] = fmincon(F, X0, conStruct.A, conStruct.b, ...
                 conStruct.Aeq, conStruct.beq, conStruct.lb, conStruct.ub, conStruct.nonLinear, options);

fmincon 默认调用内点法

Mathematica

R

Other

重要定理

AR模型协方差递推公式：

定理（Yule－Walker方程）$A R(p)$ 序列的自协方差函数满足 $$ \boldsymbol{\gamma}_n=\Gamma_n \boldsymbol{a}_n, \quad \gamma_0=\boldsymbol{\gamma}_n^T \boldsymbol{a}_n+\sigma^2, \quad n \geq p, $$ 即 $$ \begin{aligned} \gamma_k & =a_1 \gamma_{k-1}+a_2 \gamma_{k-2}+\cdots+a_p \gamma_{k-p}, \quad k \geq 1 \\ A(\mathscr{B}) \gamma_k & =0, \quad k \geq 1 \\ A(\mathscr{B}) \gamma_0 & =\gamma_0-a_1 \gamma_1-a_2 \gamma_2-\cdots-a_p \gamma_p=\sigma^2 \end{aligned} $$ 特别地，当 $n=p$ 时 $$ \Gamma_p\left(\begin{array}{c} a_1 \\ a_2 \\ \ldots \\ a_p \end{array}\right)=\left(\begin{array}{c} \gamma_1 \\ \gamma_2 \\ \ldots \\ \gamma_p \end{array}\right) $$ 记 $\phi_0=1, \phi_1=-a_1, \ldots, \phi_p=-a_p$ ，则 $A(z)=\sum_{j=0}^p \phi_j z^j, ~ \mathrm{AR}$ 模型可写成 $\sum_{j=0}^p \phi_j X_{t-j}=\varepsilon_t$ 。 Yule－Walker方程可写成 $$ \left(\begin{array}{cccc} \gamma_0 & \gamma_1 & \cdots & \gamma_p \\ \gamma_1 & \gamma_0 & \cdots & \gamma_{p-1} \\ \vdots & \vdots & & \vdots \\ \gamma_p & \gamma_{p-1} & \cdots & \gamma_0 \end{array}\right)\left(\begin{array}{c} \phi_0 \\ \phi_1 \\ \vdots \\ \phi_p \end{array}\right)=\left(\begin{array}{c} \sigma^2 \\ 0 \\ \vdots \\ 0 \end{array}\right) $$

使用Yule-Walker系数系数不止可以对AR模型进行处理，对于平稳序列的线性预测都可以使用Yule-Walker系数。

定理（Y－W系数的最小相位性） 如果实数列 $\gamma_k, k=0,1, \ldots, n$ 使得 $$ \Gamma_{n+1} \triangleq\left(\begin{array}{cccc} \gamma_0 & \gamma_1 & \cdots & \gamma_n \\ \gamma_1 & \gamma_0 & \cdots & \gamma_{n-1} \\ \vdots & \vdots & & \vdots \\ \gamma_n & \gamma_{n-1} & \cdots & \gamma_0 \end{array}\right)>0 $$ 则解出的 $n$ 阶Yuler－Walker系数 $\boldsymbol{a}_n$ 满足如下最小相位条件： $$ 1-\sum_{j=1}^n a_{n j} z^j \neq 0, \quad|z| \leq 1 . $$ 最小相位性就是以 $\boldsymbol{a}_n$ 为系数的 $\mathrm{AR}(p)$ 模型能表示成因果性线性平稳列的充分必要条件。一般的线性平稳列的自协方差列正定，所以其任意 $n$ 阶Yuler－Walker系数都满足最小相位条件。

下面的Levinson-Durbin 递推定理提供了一种递归计算Y-W系数的方法。

定理（Levinson-Durbin 递推） 如果 $\Gamma_{n+1}$ 正定，则 $\gamma_k, k=0,1, \ldots, n$ 的 $1,2, \ldots, n, n+1$ 阶Yuler－Walker系数 $\left\{a_{i j}, i=1, \ldots, n+1, j=1, \ldots, i\right\}$ 和均方误差 $\sigma_k^2$ 可以如下递推计算： $$ \begin{aligned} \sigma_0^2 & =\gamma_0 \\ a_{1,1} & =\gamma_1 / \gamma_0 \\ \sigma_k^2 & =\sigma_{k-1}^2\left(1-a_{k, k}^2\right) \\ a_{k+1, k+1} & =\frac{\gamma_{k+1}-a_{k, 1} \gamma_k-a_{k, 2} \gamma_{k-1}-\cdots-a_{k, k} \gamma_1}{\gamma_0-a_{k, 1} \gamma_1-a_{k, 2} \gamma_2-\cdots-a_{k, k} \gamma_k} \\ a_{k+1, j} & =a_{k, j}-a_{k+1, k+1} a_{k, k+1-j}, \quad 1 \leq j \leq k \end{aligned} $$ 其中 $$ \sigma_k^2=E\left(X_{k+1}-\left(a_{k, 1} X_k+\cdots+a_{k, k} X_1\right)\right)^2 $$ 是用 $X_k, X_{k-1}, \ldots, X_1$ 预测 $X_{k+1}$ 的均方误差。
注. 关于 $\sigma_k^2$ ，记 $\boldsymbol{X}_k=\left(X_k, X_{k-1}, \ldots, X_1\right)^T$ ，有： $$ \begin{aligned} \sigma_k^2 & =E\left[X_{k+1}-\boldsymbol{a}_k^T \boldsymbol{X}_k\right]^2 \\ & =E\left[\left(X_{k+1}-\boldsymbol{a}_k^T \boldsymbol{X}_k\right) X_{k+1}\right]-E\left[\left(X_{k+1}-\boldsymbol{a}_k^T \boldsymbol{X}_k\right) \boldsymbol{a}_k^T \boldsymbol{X}_k\right] \\ & =\gamma_0-\boldsymbol{a}_k^T\left(\gamma_1 \cdots \gamma_k\right)^T-0 \quad \text { (根据Y-W方程) } \\ & =\gamma_0-a_{k, 1} \gamma_1-\cdots-a_{k, k} \gamma_k \end{aligned} $$ 这是 $a_{k+1, k+1}$ 的递推公式的分母，所以 $a_{k+1, k+1}$ 的递推公式也可以写成 $$ a_{k+1, k+1}=\frac{\gamma_{k+1}-a_{k, 1} \gamma_k-a_{k, 2} \gamma_{k-1}-\cdots-a_{k, k} \gamma_1}{\sigma_k^2} $$ 注意 $\sigma_k^2$ 是用 $k$ 个历史值预报第 $k+1$ 个的最小均方误差线性预测的均方误差。

定义（偏相关系数） 如果 $\Gamma_n$ 正定，称 $a_{n, n}$ 为 $\left\{X_t\right\}$ 或 $\left\{\gamma_k\right\}$ 的 $n$ 阶偏（自）相关系数。序列 $\left\{a_{n, n}, n=1,2, \ldots\right\}$ 称为 $\left\{X_t\right\}$ 或 $\left\{\gamma_k\right\}$ 的 $n$ 阶偏（自）相关函数。

偏自相关是 $X_1$ 和 $X_{n+1}$ 之间的如下意义下的偏相关系数： $$ \begin{aligned} a_{n, n}=\operatorname{Corr}[ & X_1-L\left(X_1 \mid X_2, \ldots, X_n\right) \\ & \left.X_{n+1}-L\left(X_{n+1} \mid X_2, \ldots, X_n\right)\right] \end{aligned} $$ 即 $a_{n, n}$ 为 $X_1$ 和 $X_{n+1}$ 扣除 $X_2, \ldots, X_n$ 的线性影响后的相关系数。
设 $\left\{X_t\right\}$ 是AR $(p)$ 序列。其自协方差函数正定。由Yule－Walker方程知其 $n$ 阶 $(n \geq p) \mathrm{Y}-\mathrm{W}$ 系数为 $$ \begin{aligned} \boldsymbol{a}_n & =\left(a_1, \ldots, a_p, 0, \ldots, 0\right)^T \ & =\left(a_{n, 1}, a_{n, 2}, \ldots, a_{n, n}\right)^T, \quad n \geq p \end{aligned} $$ 其偏相关系数满足 $$ a_{n, n}= \begin{cases}a_p \neq 0, & n=p \\ 0, & n>p\end{cases} $$ 称此性质为 $\operatorname{AR}(p)$ 序列的相关系数 $p$ 后截尾。反之，如果一个零均值平稳列偏相关系数 $p$ 后截尾，则它必是AR $(p)$ 序列（见下面的定理）。

定理（AR序列的偏相关函数条件） 设零均值平稳列 $\left\{X_t\right\}$ 的自协方差函数 $\left\{\gamma_k\right\}$ 是正定序列，则 $\left\{X_t\right\}$ 是 $A R(p)$ 序列的充分必要条件是，它的偏相关系数 $\left\{a_{n, n}\right\} p$ 后截尾。

实际操作中，我们使用经验估计量估计协方差，再通过Levinson-Durbin 递推获得偏自相关系数，偏自相关系数在MLE的重参数化中具有重大作用。

论文中的一些结论

Improved Maximum Likelihood Estimation of ARMA Models:

A closed form expression of the ARMA exact likelihood function was firstly given in The exact likelihood function for a mixed autoregressive-moving average process. Afterwards, the focus shifted to finding expressions of the exact likelihood being more suitable for its computation An Algorithm for the Exact Likelihood of a Mixed Autoregressive-Moving Average Process,Computation of the exact likelihood function of an arima process