自回归模型¶

原理讲解¶

到目前为止，我们强调的是以回归模型来推断因果关系。如果从一个客户的角度仅关心某变量（比如股价）的未来值，则可以用该变量的过去值来预测其未来值(因为时间序列一般存在自相关)，这种模型被称为“单变量时间序列”。此时，可以不必理会因果关系，只考虑相关关系即可。比如，看到街上有人带伞，可以预测今天要下雨，但行人带伞并不导致下雨。对于样本数据 $y _ { 1 } , y _ { 2 } ,… y _T$，最简单的预测方法为一阶自回归 $AR(1)$:

其中，扰动项 $\varepsilon_t$ 为白噪声，满足零期望 $E \left( \varepsilon _ { t } \right) = 0$，同方差 $Var \left( \varepsilon _ { t } \right) = \sigma _ \varepsilon ^ { 2 }$，且无自相关 $Cov \left( \varepsilon _ { t } , \varepsilon _ { s } \right) =0$，$t \ne s$。假设$| \beta _ { 1 } | < 1$，则 {$y _ { t }$}为渐近独立的平稳过程。由于 $y _ { t - 1 }$ 依赖于 {$\varepsilon _ { t - 1 }$,…$\varepsilon _ { 1 }$}，而 $\varepsilon _ { t }$ 与 {$\varepsilon _ { t - 1 }$,…$\varepsilon _ { 1 }$} 不相关，故 $y _ { t - 1 }$ 与 $\varepsilon _ { t }$ 不相关，因此用OLS估计方程（1）是一致的。然而，使用 OLS 只能用观测值 $t = 2$,…,T 进行回归，将损失一个样本容量。

为了提高估计效率，考虑最大似然估计。为此，进一步假设 $\varepsilon_{t}$ 为独立同分布，且服从正态分布 $N\left(0, \sigma_{\varepsilon}^{2}\right)$。由于 ${y_t}$ 为平稳过程，其期望与方差均不随时间而变。对方程（1）两边同时取期望可得

经移项整理可知，{$y_t$} 的无条件期望为 $\frac{\beta_{0}}{1-\beta_{1}}$。对方程（1）两边同时取方差可得

故 $y_t$ 的无条件方差为 $\frac{\sigma_{\varepsilon}^{2}}{1-\beta_{1}^{2}}$。因此，$y_1$ 服从正态分布 $N\left(\frac{\beta_{0}}{1-\beta_{1}}, \frac{\sigma_{\varepsilon}^{2}}{1-\beta_{1}^{2}}\right)$，其（无条件）密度函数为

在给定 $y_1$ 的条件下，$y_2$ 的条件分布为 $y_{2} | y_{1} \sim N\left(\boldsymbol{\beta}_{0}+\boldsymbol{\beta}_{1} y_{1}, \sigma_{\varepsilon}^{2}\right)$，其条件密度为

则 $y_1$ 与 $y_2$ 的联合分布密度为

依此类推得 $y_{t} | y_{t-1} \sim N\left(\beta_{0}+\beta_{1} y_{t-1}, \sigma_{\varepsilon}^{2}\right), t=2, \cdots, T$，写出整个样本数据 $\left\{y_{1}, y_{2}, \cdots, y_{T}\right\}$ 的联合概率密度（即似然函数）为

取对数可得对数似然函数：

代入 $f_{y_{1}}\left(y_{1}\right)$ 与 $f_{y_{1} | y_{t-1}}\left(y_{t} | y_{t-1}\right)$ 的表达式可得

寻找最优参数 $\left(\beta_{0}, \beta_{1}, \sigma_{\varepsilon}^{2}\right)$，使得 $ln L$ 最大化，这个估计量被称为“精确最大似然估计量”。由于 $ln L$ 为 $\left(\beta_{0}, \beta_{1}, \sigma_{\varepsilon}^{2}\right)$）的非线性函数，故需要使用迭代法进行数值计算。

如果样本容量 $T$ 较大，则第一个观测值对似然函数的贡献较小，可以忽略。此时，相当于考虑在第一个观测值 $y_1$ 给定的情况下，{$y_2,\cdots,y_T$} 的条件分布，则对数似然函数简化为

最大化上式所得到的估计量称为“条件最大似然估计量”（ conditional MLE）。

更一般地，可以考虑 $p$ 阶自回归模型，记为 $AR(p)$：

对 $AR(p)$ 的估计与 $AR(1)$ 类似。在使用“条件 MLE”时，考虑在给定 {$y_1,\cdots,y_p$} 的情况下，{$y_p+1,\cdots,y_T$} 的条件分布。由于条件 MLE 等价于 OLS，而后者并不依赖于正态性假定，故条件 MLE 的一致性也不依赖于正态性假定。

然而，通常我们并不知道滞后期 $p$。因此，如何估计 $\hat{p}$ 在实践中有着重要的意义。

估计滞后期 $p$ 有两种方法：

1. 由大到小的序贯 $t$ 规则
2. 使用信息准则：AIC、BIC 准则等

在实践中，可以结合以上两种方法来确定 $\hat{p}$。如果二者不一致，为了保守起见（即尽量避免遗漏变量偏差），可取二者滞后阶数的大者。

最大似然估计原理¶

如果回归模型存在非线性，可能不方便使用 OLS，这时常常采用最大似然估计法（MLE）或非线性最小二乘法（NLS）。

假设随机向量 $y$ 的概率密度函数为 $f(y; \theta)$，其中 $\theta$ 为 $K$ 维未知参数向量，通过抽取随机样本 ${y_1,\cdots,y_n}$ 来估计 $\theta$。假设${y_1,\cdots,y_n}$ 为独立同分布（iid），则样本数据的联合密度函数为 $f(y_1; \theta)f(y_2; \theta) \cdots f(y_n; \theta)$。

在抽样之前，${y_1,\cdots,y_n}$ 被视为随机向量。抽样之后，${y_1,\cdots,y_n}$ 就有了特定的样本值。因此，可以将样本的联合密度函数视为在 ${y_1,\cdots,y_n}$ 给定的情况下，未知参数 $\theta$ 的函数。定义“似然函数”（likelihood function）为

由此可知，似然函数与联合密度函数完全相等，只是 $\theta$ 与 ${y_1,\cdots,y_n}$ 的角色互换，即把 $\theta$ 作为自变量，而视 ${y_1,\cdots,y_n}$ 为已给定的。为了运算方便，常把似然函数取对数，将乘积的形式转化为求和的形式

“最大似然估计法"（Maximum Likelihood Estimation，简记 MLE 或 ML）来源于一个简单而深刻的思想：给定样本取值后，该样本最有可能来自参数 $\theta$ 为何值的总体。换言之，寻找 $\hat{\boldsymbol{\theta}}_{\mathrm{ML}}$，使得观测到样本数据的可能性最大，即最大化对数似然函数（loglikelihoodfunction）：

在数学上，常把最大似然估计量 $\hat{\boldsymbol{\theta}}_{\mathrm{ML}}$ 写为

其中，“argmax"（即argument of the maximum）表示能使 $ln L(y; \theta)$ 最大化的 $\theta$ 取值。假设存在唯一内点解，则该无约束极值问题的一阶条件为

原理实现¶

statsmodels 库实现¶

详情请参考：AutoReg 模型

matlab实现¶

使用 OLS 估计代码¶

ARxx.m

fminsearch 函数（最大似然估计）¶

使用¶

fminsearch 函数用来求解多维无约束的线性优化问题，用 derivative-free 的方法找到多变量无约束函数的最小值

语法：

x = fminsearch(fun,x0)   
     x = fminsearch(fun,x0,options)    
     [x,fval] = fminsearch(...)
     [x,fval,exitflag] = fminsearch(...)
     [x,fval,exitflag,output] = fminsearch(...)

解释：

fminsearch 能够从一个初始值开始，找到一个标量函数的最小值。通常被称为无约束非线性优化

• x = fminsearch(fun,x0) 从 x0 开始，找到函数 fun 中的局部最小值 x，x0 可以是标量，向量，矩阵。fun 是一个函数句柄

• x = fminsearch(fun,x0,options) 以优化参数指定的结构最小化函数，可以用 optimset 函数定义这些参数。（见 matlab help）

• [x,fval] = fminsearch(...)返回在结果 x 出的目标函数的函数值

• [x,fval,exitflag] = fminsearch(...) 返回 exitflag 值来表示 fminsearch 退出的条件：

  • 1--函数找到结果 x
  • 0--函数最大功能评价次数达到，或者是迭代次数达到
  • -1--算法由外部函数结束

• [x,fval,exitflag,output] = fminsearch(...) 返回一个结构输出 output，包含最优化函数的信息：

  • output.algorithm 使用的优化算法
  • output.funcCount 函式计算次数
  • output.iterations 迭代次数
  • output.message 退出信息

优化参数选项options¶

可以通过 optimset 函数设置或改变这些参数。其中有的参数适用于所有的优化算法，有的则只适用于大型优化问题，另外一些则只适用于中型问题。

首先描述适用于大型问题的选项。这仅仅是一个参考，因为使用大型问题算法有一些条件。对于 fminunc（fminsearch）函数来说，必须提供梯度信息。

LargeScale – 当设为'on'时使用大型算法，若设为'off'则使用中型问题的算法。

适用于大型和中型算法的参数：

• Diagnostics – 打印最小化函数的诊断信息。

• Display – 显示水平。选择'off'，不显示输出；选择'iter'，显示每一步迭代过程的输出；选择'final'，显示最终结果。打印最小化函数的诊断信息。

• GradObj – 用户定义的目标函数的梯度。对于大型问题此参数是必选的，对于中型问题则是可选项。

• MaxFunEvals – 函数评价的最大次数。

• MaxIter – 最大允许迭代次数。

• TolFun – 函数值的终止容限。

• TolX – x 处的终止容限。

只用于大型算法的参数：

• Hessian – 用户定义的目标函数的 Hessian 矩阵。

• HessPattern –用于有限差分的 Hessian 矩阵的稀疏形式。若不方便求 fun 函数的稀疏 Hessian 矩阵 H，可以通过用梯度的有限差分获得的 H 的稀疏结构（如非零值的位置等）来得到近似的 Hessian 矩阵 H。若连矩阵的稀疏结构都不知道，则可以将 HessPattern 设为密集矩阵，在每一次迭代过程中，都将进行密集矩阵的有限差分近似（这是缺省设置）。这将非常麻烦，所以花一些力气得到 Hessian 矩阵的稀疏结构还是值得的。

• MaxPCGIter – PCG 迭代的最大次数。

• PrecondBandWidth – PCG 前处理的上带宽，缺省时为零。对于有些问题，增加带宽可以减少迭代次数。

• TolPCG – PCG 迭代的终止容限。

• TypicalX – 典型 x 值。

只用于中型算法的参数：

• DerivativeCheck – 对用户提供的导数和有限差分求出的导数进行对比。

• DiffMaxChange – 变量有限差分梯度的最大变化。

• DiffMinChange - 变量有限差分梯度的最小变化。

• LineSearchType – 一维搜索算法的选择。

• exitflag:描述退出条件

  • exitflag>0 表示目标函数收敛于解 x 处。
  • exitflag=0 表示已经达到函数评价或迭代的最大次数。
  • exitflag<0 表示目标函数不收敛。

• output:该参数包含下列优化信息：

  • output.iterations – 迭代次数。
  • output.algorithm – 所采用的算法。
  • output.funcCount – 函数评价次数。
  • output.cgiterations – PCG 迭代次数（只适用于大型规划问题）。
  • output.stepsize – 最终步长的大小（只用于中型问题）。
  • output.firstorderopt – 一阶优化的度量：解 x 处梯度的范数。

• fminunc（fminsearch）为中型优化算法的搜索方向提供了 4 中算法，由 options 中的参数 HessUpdate 控制

  • HessUpdate=‘bfgs’（默认值），为拟牛顿的 BFGS 法
  • HessUpdate='dfp'为拟牛顿 DFP 法
  • HessUpdate=‘steepdesc’最速下降法

• fminunc（fminsearch）中为中型优化算法的步长一维搜索提供了两种算法，由 options 中参数 LineSearchType 控制

  • LineSearchType='quadcubic'混合的二次和三次多项式插值
  • LineSearchType=‘cubicpoly’三次多项式插值

使用 fminunc 和 fminsearch 都可能会得到局部最优解。

fminsearch函数的缺陷¶

1. 依赖于初值
2. 不同目标函数需要选择不同的算法
3. 仅是局部最优算法，并不是全局最优
4. 迭代法与解析解的速度不是一个量级

fminsearch函数实例¶

1. 求解 $z=(x+1)^2+(y-1)^2$ 的最小值

三维图像显示：

可能不太好从图像里看出极值，这里再举一个例子：

$z=-\sqrt{(x+1)^2+(y-1)^2}$ 三维图

2. 求解 $\ln L=-\frac{1}{2} \ln 2 \pi-\frac{1}{2} \ln \left[\sigma_{\varepsilon}^{2} /\left(1-\beta_{1}^{2}\right)\right]-\frac{\left[y_{1}-\left(\beta_{0} /\left(1-\beta_{1}\right)\right)\right]^{2}}{2 \sigma_{\varepsilon}^{2} /\left(1-\beta_{1}^{2}\right)}-\frac{T-1}{2} \ln 2 \pi-\frac{T-1}{2} \ln \sigma_{\varepsilon}^{2}-\sum_{t=2}^{T} \frac{\left(y_{t}-\beta_{0}-\beta_{1} y_{t-1}\right)^{2}}{2 \sigma_{\varepsilon}^{2}}$ 的最大值

最终结果为：