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) 有时被称为 随机差分方程。)
Example 32.1
例如,
某个家庭的劳动收入对数,或
某个经济体的货币需求对数。
在两种情况任意一种下,(32.1) 显示当前值是通过上一个值的线性函数以及一个独立同分布的冲击
(我们使用
一旦我们指定一个初始条件
为了使分析变得更简单,我们将假设
过程
是 独立同分布 且符合标准正态分布,初始条件
从正态分布 中抽取,初始条件
独立于 。
32.2.1. 移动平均表示#
从时间
如果我们一直追溯到零时,我们得到
方程 (32.2) 显示
参数,
初始条件
和从时间
到现在的冲击 。
在整个过程中,符号
32.2.2. 分布动态#
这个模型的一个好处是很容易追踪一系列分布
让我们看看我们如何做到这一点。
首先我们指出,对于每个时间
这是从 (32.2) 显见的,因为独立正态随机变量的线性组合是正态分布的。
鉴于
设
我们可以从 (32.2) 确定这些值,或者我们可以使用以下递归表达式:
通过分别取等式两侧的期望和方差,从 (32.1) 得到这些表达式。
在计算第二个表达式时,我们利用了
(这是根据我们的假设和 (32.2) 得出的。)
给定 (32.2) 中的动态和初始条件
下面的代码利用以上所得来追踪边缘分布序列
参数是
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()
32.3. 平稳性和渐近稳定性#
当我们使用模型来研究现实世界时,通常希望我们的模型具有清晰准确的预测。
对于动态问题,准确的预测与稳定性有关。
例如,如果一个动态模型预测通货膨胀总是收敛到某种稳态,那么这个模型提供了一个明确的预测。
(预测可能是错误的,但即便如此,它也是有帮助的,因为我们可以判断模型的质量。)
注意,在上图中,序列
如果我们进一步研究未来的情况,这一点就更加明显:
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()
此外,极限不依赖于初始条件。
例如,另一个密度序列也会收敛到相同的极限。
事实上,可以很容易地证明,只要
为了证明这一点,我们只需查看前两个矩的动态, 如 (32.3) 中所给出的。
当
(可以通过阅读我们的 一维动力学讲座来了解确定性收敛的背景。)
因此
我们可以使用以下代码确认这对于上面的序列是有效的。
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()
请注意,根据上述参数,我们看到序列
我们看到,至少对于这些参数,AR(1) 模型具有很强的稳定性特性。
32.3.1. 平稳分布#
让我们试着更深入地理解极限分布
平稳分布是 AR(1) 过程更新规则的一个“不动点”。
换句话说,如果
另一种针对当前情况的说法是:如果一个在
(当然,我们假设
事实上,可以证明
因此,当
32.4. 遍历性#
不同的作者使用遍历性这一概念有不同的方式。
在当前情况中可以理解为:即使
特别是,时间序列的平均值收敛于平稳分布下的期望值。
实际上,可以证明,只要
只要右侧的积分是有限且良定义的。
注意:
在 (32.6) 中,收敛性以概率1成立。
由 [Meyn and Tweedie, 2009] 编写的教科书是关于遍历性的经典参考。
Example 32.2
如果我们考虑恒等函数
即,时间序列样本均值收敛于平稳分布的均值。
出于多种原因,遍历性非常重要。
例如,(32.6) 可用来测试理论。
在这个方程中,我们可以使用观测数据来评估 (32.6) 的左侧。
我们可以使用理论的 AR(1) 模型来计算右侧。
如果
32.5. 练习#
Exercise 32.1
假设
随机变量的
当这个随机变量是
这里
根据 (32.6), 对于任何
当
通过模拟验证一系列
Solution to Exercise 32.1
这是一个解法:
from numba import njit
from scipy.special import factorial2
@njit
def sample_moments_ar1(k, m=100_000, mu_0=0.0, sigma_0=1.0, seed=1234):
np.random.seed(seed)
sample_sum = 0.0
x = mu_0 + sigma_0 * np.random.randn()
for t in range(m):
sample_sum += (x - mu_star)**k
x = a * x + b + c * np.random.randn()
return sample_sum / m
def true_moments_ar1(k):
if k % 2 == 0:
return std_star**k * factorial2(k - 1)
else:
return 0
k_vals = np.arange(6) + 1
sample_moments = np.empty_like(k_vals)
true_moments = np.empty_like(k_vals)
for k_idx, k in enumerate(k_vals):
sample_moments[k_idx] = sample_moments_ar1(k)
true_moments[k_idx] = true_moments_ar1(k)
fig, ax = plt.subplots()
ax.plot(k_vals, true_moments, label="真实矩")
ax.plot(k_vals, sample_moments, label="样本矩")
ax.legend()
plt.show()
Exercise 32.2
编写一个你自己的一维核密度估计,用于从样本中估计概率密度。
将其写为一个类,该类在初始化时接受数据
对于
编写该类,使得带宽默认遵循 Silverman 的规则(参见此页面中的“经验规则”讨论)。通过以下步骤测试你编写的类:
从分布
中模拟数据在适当的范围内绘制核密度估计
在同一图形上绘制
的密度
给定如下分布
使用
对结果进行总结。(你认为这是这些分布的良好估计吗?)
Solution to Exercise 32.2
这是一个解答:
K = norm.pdf
class KDE:
def __init__(self, x_data, h=None):
if h is None:
c = x_data.std()
n = len(x_data)
h = 1.06 * c * n**(-1/5)
self.h = h
self.x_data = x_data
def f(self, x):
if np.isscalar(x):
return K((x - self.x_data) / self.h).mean() * (1/self.h)
else:
y = np.empty_like(x)
for i, x_val in enumerate(x):
y[i] = K((x_val - self.x_data) / self.h).mean() * (1/self.h)
return y
def plot_kde(ϕ, x_min=-0.2, x_max=1.2):
x_data = ϕ.rvs(n)
kde = KDE(x_data)
x_grid = np.linspace(-0.2, 1.2, 100)
fig, ax = plt.subplots()
ax.plot(x_grid, kde.f(x_grid), label="估计值")
ax.plot(x_grid, ϕ.pdf(x_grid), label="真实密度")
ax.legend()
plt.show()
from scipy.stats import beta
n = 500
parameter_pairs= (2, 2), (2, 5), (0.5, 0.5)
for α, β in parameter_pairs:
plot_kde(beta(α, β))
我们可以看到,当基础分布是平滑时,核密度估计是有效的,但在其他情况下则效果不佳。
Exercise 32.3
在本讲中我们讨论了以下事实:对于
其中
通过仿真确认这一点,至少是近似的。设定
首先,使用上述真实分布绘制
其次,在同一个图形上(用不同的颜色),绘制
从分布
中生成 次 抽样使用方程
更新全部的数据使用
的结果样本通过核密度估计产生的估计密度。
尝试用
Solution to Exercise 32.3
这是我们的解答
a = 0.9
b = 0.0
c = 0.1
μ = -3
s = 0.2
n = 2000 # 样本数
# 理论上的分布
μ_next = a * μ + b
s_next = np.sqrt(a**2 * s**2 + c**2)
# 理论密度
ψ = norm(μ, s)
ψ_next = norm(μ_next, s_next)
# 使用仿真值来估计ψ_{t+1}
x_draws = ψ.rvs(n) # 从ψ生成样本
x_draws_next = a * x_draws + b + c * np.random.randn(n) # 更新规则
# 对仿真结果进行核密度估计
kde = KDE(x_draws_next)
# 绘制结果
x_grid = np.linspace(μ - 1, μ + 1, 100)
fig, ax = plt.subplots()
ax.plot(x_grid, ψ.pdf(x_grid), label="$\psi_t$")
ax.plot(x_grid, ψ_next.pdf(x_grid), label="$\psi_{t+1}$")
ax.plot(x_grid, kde.f(x_grid), label="estimate of $\psi_{t+1}$")
ax.legend()
plt.show()
正如预测的那样,模拟分布与理论分布大致吻合。