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)#Xt+1=aXt+b+cWt+1

其中 a,b,c 是标量参数

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

Example 32.1

例如,Xt 可能是

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

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

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

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

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

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

  • 过程 {Wt}独立同分布 且符合标准正态分布,

  • 初始条件 X0 从正态分布 N(μ0,v0) 中抽取,

  • 初始条件 X0 独立于 {Wt}

32.2.1. 移动平均表示#

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

Xt=aXt1+b+cWt=a2Xt2+ab+acWt1+b+cWt=a3Xt3+a2b+a2cWt2+b+cWt=

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

(32.2)#Xt=atX0+bj=0t1aj+cj=0t1ajWtj

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

  • 参数,

  • 初始条件 X0

  • 从时间 t=1 到现在的冲击 W1,Wt

在整个过程中,符号 ψt 将用来指代这个随机变量 Xt 的概率密度。

32.2.2. 分布动态#

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

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

首先我们指出,对于每个时间 tXt 是正态分布的。

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

鉴于 Xt 是正态分布的,如果我们能确定它的一阶矩和二阶矩,就可以知道完整的分布 ψt

μtvt 分别表示 Xt 的均值和方差。

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

(32.3)#μt+1=aμt+bextvt+1=a2vt+c2

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

在计算第二个表达式时,我们利用了 XtWt+1 的独立性。

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

给定 (32.2) 中的动态和初始条件 μ0,v0,我们得到 μt,vt,因此

ψt=N(μt,vt)

下面的代码利用以上所得来追踪边缘分布序列 {ψ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. 平稳性和渐近稳定性#

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

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

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

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

注意,在上图中,序列 {ψ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)#μ:=b1av=c21a2

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

因此

(32.5)#ψtψ=N(μ,v)t

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

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

请注意,根据上述参数,我们看到序列 {ψt} 收敛到 ψ

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

32.3.1. 平稳分布#

让我们试着更深入地理解极限分布 ψ

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

换句话说,如果 ψt 是平稳的,那么对所有 jψt+j=ψtN 时成立。

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

XtψaXt+b+cWt+1ψ

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

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

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

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

32.4. 遍历性#

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

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

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

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

(32.6)#1mt=1mh(Xt)h(x)ψ(x)dxm

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

注意:

Example 32.2

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

1mt=1mXtxψ(x)dx当 m

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

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

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

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

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

如果 1mt=1mXt 即使在大量观测下也不接近 ψ(x),那么我们的理论便有可能是错误的,我们将需要修订它。

32.5. 练习#

Exercise 32.1

假设 k 是一个自然数。

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

Mk:=E[(XEX)k]

当这个随机变量是 N(μ,σ2) 时,已知

Mk={0kσk(k1)!!k

这里 n!!双阶乘

根据 (32.6), 对于任何 kN

1mt=1m(Xtμ)kMk

m 是较大时。

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

Exercise 32.2

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

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

f(x)=1hni=1nK(xXih)

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

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

  1. 从分布 ϕ 中模拟数据 X1,,Xn

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

  3. 在同一图形上绘制 ϕ 的密度

给定如下分布 ϕ 类型:

使用 n=500

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

Exercise 32.3

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

Xt+1=aXt+b+cWt+1

其中 {Wt} 独立同分布且标准正态,

ψt=N(μ,s2)ψt+1=N(aμ+b,a2s2+c2)

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

  • a=0.9

  • b=0.0

  • c=0.1

  • μ=3

  • s=0.2

首先,使用上述真实分布绘制 ψtψt+1

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

  1. 从分布 N(μ,s2) 中生成 nXt 抽样

  2. 使用方程 Xt+1=aXt+b+cWt+1 更新全部的数据

  3. 使用 Xt+1 的结果样本通过核密度估计产生的估计密度。

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