32. AR(1) 过程#

32.1. 概览#

在这个讲座中,我们将研究一类非常简单的随机模型,称为 AR(1) 过程。

这些简单的模型在经济研究中一次又一次地被用来表示诸如

  • 劳动收入

  • 股息

  • 生产力等序列的动态。

AR(1) 过程不仅有广泛的应用,还可以帮助我们理解很多非常重要的概念。

让我们从一些函数库导入开始:

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5)  # 设置默认图像大小

import matplotlib as mpl
FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf"
mpl.font_manager.fontManager.addfont(FONTPATH)
plt.rcParams['font.family'] = ['Source Han Serif SC']

32.2. AR(1) 模型#

AR(1) 模型(一阶自回归模型)形式是:

(32.1)#\[X_{t+1} = a X_t + b + c W_{t+1}\]

其中 \(a, b, c\) 是标量参数

(方程 (32.1) 有时被称为 随机差分方程。)

Example 32.1

例如,\(X_t\) 可能是

  • 某个家庭的劳动收入对数,或

  • 某个经济体的货币需求对数。

在两种情况任意一种下,(32.1) 显示当前值是通过上一个值的线性函数以及一个独立同分布的冲击 \(W_{t+1}\) 来演变的。

(我们使用 \(t+1\) 作为 \(W_{t+1}\) 的下标,因为在时间 \(t\) 时这个随机变量还未被观察到。)

一旦我们指定一个初始条件 \(X_0\),我们就可以用(32.1) 生成一个时间序列 \(\{ X_t\}\)

为了使分析变得更简单,我们将假设

  • 过程 \(\{ W_t \}\)独立同分布 且符合标准正态分布,

  • 初始条件 \(X_0\) 从正态分布 \(N(\mu_0, v_0)\) 中抽取,

  • 初始条件 \(X_0\) 独立于 \(\{ W_t \}\)

32.2.1. 移动平均表示#

从时间 \(t\) 向后迭代,我们得到

\[ X_t = a X_{t-1} + b + c W_t = a^2 X_{t-2} + a b + a c W_{t-1} + b + c W_t = a^3 X_{t-3} + a^2 b + a^2 c W_{t-2} + b + c W_t = \cdots \]

如果我们一直追溯到零时,我们得到

(32.2)#\[X_t = a^t X_0 + b \sum_{j=0}^{t-1} a^j + c \sum_{j=0}^{t-1} a^j W_{t-j}\]

方程 (32.2) 显示 \(X_t\) 是一个明确定义的随机变量,其值取决于

  • 参数,

  • 初始条件 \(X_0\)

  • 从时间 \(t=1\) 到现在的冲击 \(W_1, \ldots W_t\)

在整个过程中,符号 \(\psi_t\) 将用来指代这个随机变量 \(X_t\) 的概率密度。

32.2.2. 分布动态#

这个模型的一个好处是很容易追踪一系列分布 \(\{ \psi_t \}\),这些分布对应于时间 序列 \(\{ X_t\}\) 。具体来说,我们可以在每个日期 \(t\) 上追踪观测到的 \(X_t\) 的边缘分布。

让我们看看我们如何做到这一点。

首先我们指出,对于每个时间 \(t\)\(X_t\) 是正态分布的。

这是从 (32.2) 显见的,因为独立正态随机变量的线性组合是正态分布的。

鉴于 \(X_t\) 是正态分布的,如果我们能确定它的一阶矩和二阶矩,就可以知道完整的分布 \(\psi_t\)

\(\mu_t\)\(v_t\) 分别表示 \(X_t\) 的均值和方差。

我们可以从 (32.2) 确定这些值,或者我们可以使用以下递归表达式:

(32.3)#\[\mu_{t+1} = a \mu_t + b \quad ext{和} \quad v_{t+1} = a^2 v_t + c^2\]

通过分别取等式两侧的期望和方差,从 (32.1) 得到这些表达式。

在计算第二个表达式时,我们利用了 \(X_t\)\(W_{t+1}\) 的独立性。

(这是根据我们的假设和 (32.2) 得出的。)

给定 (32.2) 中的动态和初始条件 \(\mu_0, v_0\),我们得到 \(\mu_t, v_t\),因此

\[ \psi_t = N(\mu_t, v_t) \]

下面的代码利用以上所得来追踪边缘分布序列 \(\{ \psi_t \}\)

参数是

a, b, c = 0.9, 0.1, 0.5

mu, v = -3.0, 0.6  # 初始条件 mu_0, v_0

以下是分布序列:

from scipy.stats import norm

sim_length = 10
grid = np.linspace(-5, 7, 120)

fig, ax = plt.subplots()

for t in range(sim_length):
    mu = a * mu + b
    v = a**2 * v + c**2
    ax.plot(grid, norm.pdf(grid, loc=mu, scale=np.sqrt(v)),
            label=fr"$\psi_{t}$",
            alpha=0.7)

ax.legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=1)

plt.show()
_images/ec1937a7e974c38bb96b0cacf27b2e2f49d36496015d2616ab738d581254da88.png

32.3. 平稳性和渐近稳定性#

当我们使用模型来研究现实世界时,通常希望我们的模型具有清晰准确的预测。

对于动态问题,准确的预测与稳定性有关。

例如,如果一个动态模型预测通货膨胀总是收敛到某种稳态,那么这个模型提供了一个明确的预测。

(预测可能是错误的,但即便如此,它也是有帮助的,因为我们可以判断模型的质量。)

注意,在上图中,序列 \(\{ \psi_t \}\) 似乎正在收敛到一个极限分布,这表明存在某种稳定性。

如果我们进一步研究未来的情况,这一点就更加明显:

def plot_density_seq(ax, mu_0=-3.0, v_0=0.6, sim_length=40):
    mu, v = mu_0, v_0
    for t in range(sim_length):
        mu = a * mu + b
        v = a**2 * v + c**2
        ax.plot(grid,
                norm.pdf(grid, loc=mu, scale=np.sqrt(v)),
                alpha=0.5)

fig, ax = plt.subplots()
plot_density_seq(ax)
plt.show()
_images/53713afb5d43eb37fa8a8508e6fb1dad7a5abb060428622e9e8b996138aeb5cd.png

此外,极限不依赖于初始条件。

例如,另一个密度序列也会收敛到相同的极限。

fig, ax = plt.subplots()
plot_density_seq(ax, mu_0=4.0)
plt.show()
_images/11646ff605e31d15004f525f702e66f2a15619a0b53ac65338c19b8c098bb439.png

事实上,可以很容易地证明,只要 \(|a| < 1\),不管初始条件如何,都会发生这种收敛。

为了证明这一点,我们只需查看前两个矩的动态, 如 (32.3) 中所给出的。

\(|a| < 1\) 时,这些序列会收敛到各自的极限

(32.4)#\[\mu^* := \frac{b}{1-a} \quad 和 \quad v^* = \frac{c^2}{1 - a^2}\]

(可以通过阅读我们的 一维动力学讲座来了解确定性收敛的背景。)

因此

(32.5)#\[\psi_t \to \psi^* = N(\mu^*, v^*) \quad 当 t \to \infty\]

我们可以使用以下代码确认这对于上面的序列是有效的。

fig, ax = plt.subplots()
plot_density_seq(ax, mu_0=4.0)

mu_star = b / (1 - a)
std_star = np.sqrt(c**2 / (1 - a**2))  # v_star的平方根
psi_star = norm.pdf(grid, loc=mu_star, scale=std_star)
ax.plot(grid, psi_star, 'k-', lw=2, label=r"$\psi^*$")
ax.legend()

plt.show()
_images/76aa123a72f2b3168d14618c5e98fa814f90331d309f69f425aff25ccd62a992.png

请注意,根据上述参数,我们看到序列 \(\{ \psi_t \}\) 收敛到 \(\psi^*\)

我们看到,至少对于这些参数,AR(1) 模型具有很强的稳定性特性。

32.3.1. 平稳分布#

让我们试着更深入地理解极限分布 \(\psi^*\)

平稳分布是 AR(1) 过程更新规则的一个“不动点”。

换句话说,如果 \(\psi_t\) 是平稳的,那么对所有 \(j\)\(\psi_{t+j} = \psi_t\)\(\mathbb N\) 时成立。

另一种针对当前情况的说法是:如果一个在 \(\mathbb R\) 上的概率密度 \(\psi\) 对 AR(1) 过程是平稳的,则有

\[ X_t \sim \psi \quad \implies \quad a X_t + b + c W_{t+1} \sim \psi \]

\(\psi^*\)(32.5) 中具有这一性质——验证这一点是留给读者的练习。

(当然,我们假设 \(|a| < 1\) 从而 \(\psi^*\) 是 良定义的。)

事实上,可以证明 \(\mathbb R\) 上没有其他分布具有这一性质。

因此,当 \(|a| < 1\) 时,AR(1) 模型恰好有一个平稳概率密度,那就是 \(\psi^*\)

32.4. 遍历性#

不同的作者使用遍历性这一概念有不同的方式。

在当前情况中可以理解为:即使 \(\{X_t\}\) 不是独立同分布的,大数定律也是有效的。

特别是,时间序列的平均值收敛于平稳分布下的期望值。

实际上,可以证明,只要 \(|a| < 1\),我们就有

(32.6)#\[\frac{1}{m} \sum_{t = 1}^m h(X_t) \to \int h(x) \psi^*(x) dx \quad 当 m \to \infty\]

只要右侧的积分是有限且良定义的。

注意:

Example 32.2

如果我们考虑恒等函数 \(h(x) = x\),我们得到

\[ \frac{1}{m} \sum_{t = 1}^m X_t \to \int x \psi^*(x) dx \quad \text{当 } m \to \infty \]

即,时间序列样本均值收敛于平稳分布的均值。

出于多种原因,遍历性非常重要。

例如,(32.6) 可用来测试理论。

在这个方程中,我们可以使用观测数据来评估 (32.6) 的左侧。

我们可以使用理论的 AR(1) 模型来计算右侧。

如果 \(\frac{1}{m} \sum_{t = 1}^m X_t\) 即使在大量观测下也不接近 \(\psi^*(x)\),那么我们的理论便有可能是错误的,我们将需要修订它。

32.5. 练习#

Exercise 32.1

假设 \(k\) 是一个自然数。

随机变量的 \(k\) 阶中心矩定义为

\[ M_k := \mathbb E [ (X - \mathbb E X )^k ] \]

当这个随机变量是 \(N(\mu, \sigma^2)\) 时,已知

\[\begin{split} M_k = \begin{cases} 0 & 如果 k 是奇数 \\ \sigma^k (k-1)!! & 如果 k 是偶数 \end{cases} \end{split}\]

这里 \(n!!\)双阶乘

根据 (32.6), 对于任何 \(k \in \mathbb N\)

\[ \frac{1}{m} \sum_{t = 1}^m (X_t - \mu^* )^k \approx M_k \]

\(m\) 是较大时。

通过模拟验证一系列 \(k\),使用本讲中的默认参数。

Exercise 32.2

编写一个你自己的一维核密度估计,用于从样本中估计概率密度。

将其写为一个类,该类在初始化时接受数据 \(X\) 和带宽 \(h\),并提供一个方法 \(f\),使得

\[ f(x) = \frac{1}{hn} \sum_{i=1}^n K \left( \frac{x-X_i}{h} \right) \]

对于 \(K\),使用高斯核(\(K\) 是标准正态密度函数)。

编写该类,使得带宽默认遵循 Silverman 的规则(参见此页面中的“经验规则”讨论)。通过以下步骤测试你编写的类:

  1. 从分布 \(\phi\) 中模拟数据 \(X_1, \ldots, X_n\)

  2. 在适当的范围内绘制核密度估计

  3. 在同一图形上绘制 \(\phi\) 的密度

给定如下分布 \(\phi\) 类型:

使用 \(n=500\)

对结果进行总结。(你认为这是这些分布的良好估计吗?)

Exercise 32.3

在本讲中我们讨论了以下事实:对于 \(AR(1)\) 过程

\[ X_{t+1} = a X_t + b + c W_{t+1} \]

其中 \(\{ W_t \}\) 独立同分布且标准正态,

\[ \psi_t = N(\mu, s^2) \implies \psi_{t+1} = N(a \mu + b, a^2 s^2 + c^2) \]

通过仿真确认这一点,至少是近似的。设定

  • \(a = 0.9\)

  • \(b = 0.0\)

  • \(c = 0.1\)

  • \(\mu = -3\)

  • \(s = 0.2\)

首先,使用上述真实分布绘制 \(\psi_t\)\(\psi_{t+1}\)

其次,在同一个图形上(用不同的颜色),绘制 \(\psi_{t+1}\) 如下:

  1. 从分布 \(N(\mu, s^2)\) 中生成 \(n\)\(X_t\) 抽样

  2. 使用方程 \(X_{t+1} = a X_t + b + c W_{t+1}\) 更新全部的数据

  3. 使用 \(X_{t+1}\) 的结果样本通过核密度估计产生的估计密度。

尝试用 \(n=2000\) 并确认基于抽样实验的 \(\psi_{t+1}\) 估计会收敛到理论分布。