GSoC 2017 项目：哈密顿间接推断

这是我 2017 年 Google Summer of Code 项目的总结。相关代码库包含了各种模型估计的示例。除了这个代码库之外，我还作为 GSOC 2017 项目的一部分，参与了 HamiltonianABC 及其分支的开发。另请参阅 Bayesian_Examples.jl。

GSOC 2017 项目：哈密顿蒙特卡罗和伪贝叶斯间接似然
完整数据的参数贝叶斯间接似然
项目的第一阶段
项目的第二阶段
随机波动率模型
我在 GSOC 期间遇到的问题
1. 1 困难的辅助模型
2. 2 类型稳定性问题
未来工作
参考文献

GSOC 2017 项目：哈密顿蒙特卡罗和伪贝叶斯间接似然

今年夏天，我有机会参加了 Google Summer of Code 项目。我的项目使用 Julia 语言，主要目标是实现间接推断 (A. A. Smith 1993; A. Smith 2008) 来克服在使用基于似然的估计方法时通常出现的难题（例如难以处理或计算成本很高的似然函数）。预期哈密顿蒙特卡罗将带来更高效的采样过程。

在 Tamás K. Papp 的指导下，我完成了使用间接推断 (II) 和哈密顿蒙特卡罗进行贝叶斯估计方法的重大修订。我还熟悉了使用 git、创建问题、创建代码库等操作。

在这里，我将介绍这些方法并提供一些背景信息，并更详细地讨论一个示例。

通常，当我们面临难以处理的似然函数或计算成本极高的似然函数时，可以选择使用替代的辅助模型来提取和估计感兴趣的参数。这些替代模型应该更容易处理。Drovandi 等人回顾了一系列参数贝叶斯间接推断 (pBII) 方法，我重点关注了 Gallant 和 McCulloch (2009) 提出的完整数据的参数贝叶斯间接似然 (pdBIL) 方法。pdBIL 使用辅助模型的似然函数来代替难以处理的似然函数。pdBIL 不比较汇总统计量，而是按以下方式工作

首先生成数据，一旦我们有了数据，就可以估计辅助模型的参数。然后，将估计的参数与观察/生成的数据一起代入辅助似然函数。之后，我们可以将此似然函数用于我们选择的贝叶斯方法，例如 MCMC。

为了总结该方法，首先我们有参数向量 $\theta$ 和观察到的数据 y。我们希望计算 $\ell(\theta|y)$ 的似然函数，但它难以处理或计算成本很高。在这种情况下，使用 pdBIL，我们必须找到一个辅助模型 (A)，以以下方式近似真实似然函数

首先，我们必须从之前提出的参数 $\theta$ 的数据生成过程中生成点，用 x* 表示。
然后，我们计算 x 下辅助似然的 MLE 以获得用 $\phi$ 表示的参数。

\phi(x^{\star}) = \argmax_{\phi} (x^{\star}|\phi)

在这些参数 $\phi$ 下，我们现在可以计算 $\ell_{A}(y|\phi)$ 的似然函数。期望辅助似然函数尽可能接近真实似然函数，在捕捉模型和生成数据的相关方面的意义上。

项目的第一阶段

在我的项目的第一阶段，我使用 pdBIL 编写了 Drovandi 等人提出的两个模型的代码。在计算辅助模型的似然函数后，我使用了随机游走 Metropolis-Hastings MCMC 从目标分布中采样，得到玩具模型。在这个阶段，我使用的方法是众所周知的。复制 Drovandi 等人的玩具模型的目的是找出我们以后可能面临的哪些问题，并提出一个可用的接口。此阶段的结果是 HamiltonianABC（与 Tamás K. Papp 合作）。

项目的第二阶段

在第一阶段之后，我阅读了 Betancourt (2017) 并对 Tamás K. Papp 的 DynamicHMC.jl 进行了代码修订，包括检查代码并将其与论文进行比较。除了使用哈密顿蒙特卡罗方法之外，ForwardDiff 包的前向模式自动微分也是此阶段的另一个主要因素。该项目的创新之处在于找到一种方法将所有组件组合在一起，以从中获得有效的估计。最大的问题是定义类型稳定的函数以加速采样过程。

随机波动率模型

在第二阶段之后，我为 DynamicHMC.jl 编写了经济模型的代码。随机波动率模型就是其中之一。在下一节中，我将介绍设置过程。

Ornstein-Ulenbeck 随机波动率模型的连续时间版本描述了时间 t 时的回报如何具有零均值，并且其波动率受其方差的连续时间 Ornstein-Ulenbeck 过程控制。金融产品的价值的巨大波动意味着波动率过程是变化的。这就是为什么我们需要模型中的随机元素。由于我们只能访问离散时间的数据，因此对模型进行离散化是自然的。

Ornstein-Ulenbeck 随机波动率模型的离散时间版本

y_{t} = x_{t} + \epsilon_{t} \quad \text{其中}\quad \epsilon_{t} ∼ \Chi^{2}(1)

x_{t} = \rho * x_{t-1} + \sigma * \nu_{t} \quad \text{其中}\quad \nu_{t} ∼ \mathcal N(0, 1)

离散时间版本被用作数据生成过程。其中 $y_t$ 表示收益的对数， $x_{t}$ 是方差的对数，而 $\epsilon_{t}$ 和 $\nu_{t}$ 是未观测的噪声项。

对于辅助模型，我们使用了两个回归。第一个回归是对一阶差分进行的AR(2)过程，第二个也是对原始变量进行的AR(2)过程，以捕捉水平。

"""
    lag_matrix(xs, ns, K = maximum(ns))

Matrix with differently lagged xs.
"""
function lag_matrix(xs, ns, K = maximum(ns))
    M = Matrix{eltype(xs)}(length(xs)-K, maximum(ns))
    for i ∈ ns
        M[:, i] = lag(xs, i, K)
    end
    M
end

"first auxiliary regression y, X, meant to capture first differences"
function yX1(zs, K)
    Δs = diff(zs)
    lag(Δs, 0, K), hcat(lag_matrix(Δs, 1:K, K), ones(eltype(zs), length(Δs)-K), lag(zs, 1, K+1))
end

"second auxiliary regression y, X, meant to capture levels"
function yX2(zs, K)
    lag(zs, 0, K), hcat(ones(eltype(zs), length(zs)-K), lag_matrix(zs, 1:K, K))
end

一阶差分的AR(2)过程可以概括为

给定一个序列Y，它是一阶差分的一阶差分。也就是Y在时间t的“变化的变化”。离散函数的二阶差分可以解释为连续函数的二阶导数，即函数在时间t某一点的“加速度”。在这个模型中，我们想要捕捉收益的对数的“加速度”。

原始变量的AR(2)过程需要捕捉 $\rho$ 的影响。事实证明，ρ在对一阶差分进行的AR(2)过程中影响较弱。这就是我们需要第二个辅助模型的原因。

我现在将描述使用动态哈密顿蒙特卡罗方法估计随机波动率模型中感兴趣参数所需的步骤。首先，我们需要一个可调用的Julia对象，它以DiffResult类型返回对数密度和梯度。之后，我们编写一个计算密度的函数，然后使用ForwardDiff包在包装函数中计算其梯度。

随机波动率模型所需的包

using ArgCheck
using Distributions
using Parameters
using DynamicHMC
using StatsBase
using Base.Test
using ContinuousTransformations
using DiffWrappers
import Distributions: Uniform, InverseGamma

首先，我们定义一个结构体。此结构体应包含观测数据、先验、冲击和对参数执行的转换，但组件可能会因估计模型而异。

struct StochasticVolatility{T, Prior_ρ, Prior_σ, Ttrans}
    "observed data"
    ys::Vector{T}
    "prior for ρ (persistence)"
    prior_ρ::Prior_ρ
    "prior for σ_v (volatility of volatility)"
    prior_σ::Prior_σ
    "χ^2 draws for simulation"
    ϵ::Vector{T}
    "Normal(0,1) draws for simulation"
    ν::Vector{T}
    "Transformations cached"
    transformation::Ttrans
end

在指定数据生成函数以及特定模型的一些辅助函数和附加函数后（整个模块可以在src文件夹中找到），我们可以使模型结构体可调用，返回对数密度。由于我们对参数进行了转换，因此需要logjac。

function (pp::StochasticVolatility)(θ)
    @unpack ys, prior_ρ, prior_σ, ν, ϵ, transformation = pp
    ρ, σ = transformation(θ)
    logprior = logpdf(prior_ρ, ρ) + logpdf(prior_σ, σ)
    N = length(ϵ)

    # Generating xs, which is the latent volatility process

    xs = simulate_stochastic(ρ, σ, ϵ, ν)
    Y_1, X_1 = yX1(xs, 2)
    β₁ = qrfact(X_1, Val{true}) \ Y_1
    v₁ = mean(abs2,  Y_1 - X_1*β₁)
    Y_2, X_2 = yX2(xs, 2)
    β₂ = qrfact(X_2, Val{true}) \ Y_2
    v₂ = mean(abs2,  Y_2 - X_2*β₂)
    # We work with first differences
    y₁, X₁ = yX1(ys, 2)
    log_likelihood1 = sum(logpdf.(Normal(0, √v₁), y₁ - X₁ * β₁))
    y₂, X₂ = yX2(ys, 2)
    log_likelihood2 = sum(logpdf.(Normal(0, √v₂), y₂ - X₂ * β₂))
    logprior + log_likelihood1 + log_likelihood2 + logjac(transformation, θ)
end

我们需要进行变换，因为参数位于 $\Re^{n}$ 的真子集内，但我们希望使用 $\Re^{n}$ 。ContinuousTransformation 包用于进行这些变换。我们保存这些变换，以便可调用对象保持类型稳定，从而使过程更快。

$\nu$ 和 $\epsilon$ 是随机变量，我们在变换后使用它们来模拟观测点。这样，模拟变量在参数上是连续的，后验是可微的。

给定定义的函数，我们现在可以开始估计和采样过程。

RNG = Base.Random.GLOBAL_RNG
# true parameters and observed data
ρ = 0.8
σ = 0.6
y = simulate_stochastic(ρ, σ, 10000)
# setting up the model
model = StochasticVolatility(y, Uniform(-1, 1), InverseGamma(1, 1), 10000)
# we start the estimation process from the true values
θ₀ = inverse(model.transformation, (ρ, σ))
# wrap for gradient calculations
fgw = ForwardGradientWrapper(model, θ₀)
# sampling
sample, tuned_sampler = NUTS_tune_and_mcmc(RNG, fgw, 5000; q = θ₀)

以下图表显示了参数的结果。

rho_plot

sigma_plot

通过分析上述图表，我们可以说后验值与真实值相当接近。同样值得一提的是，先验不影响后验值。

我在 GSOC 期间遇到的问题

1 困难的辅助模型

真实模型是由 Rayner 和 MacGillivray (2002) 描述的 g-and-k 分位数函数。
辅助模型是三成分正态混合模型。

我们在这个模型中遇到了严重的问题。

首先，我编写了有限成分正态混合模型的最大似然估计 (MLE) 代码，该代码根据观测数据和所需的混合数计算正态分布的均值、方差和权重。使用 g-and-k 分位数函数时，我遇到了所谓的“隔离”现象，这意味着一个观测点是异常值，权重为 1，其他观测点的权重为 $\theta$ ，导致方差等于 $\theta$ 。有一些方法可以解决隔离问题，但感兴趣的参数仍然没有收敛到真实值。这个模型还需要进一步研究。

2 类型稳定性问题

为了有效地使用自动微分方法，我必须将函数编写为类型稳定的，否则采样函数将需要花费数小时才能运行。请参阅以下示例。

这不是类型稳定的。

function simulate_stochastic(ρ, σ, ϵs, νs)
    N = length(ϵs)
    @argcheck N == length(νs)
    xs = Vector(N)
    for i in 1:N
        xs[i] = (i == 1) ? νs[1]*σ*(1 - ρ^2)^(-0.5) : (ρ*xs[i-1] + σ*νs[i])
    end
    xs + log.(ϵs) + 1.27
end

这是类型稳定的。

function simulate_stochastic(ρ, σ, ϵs, νs)
    N = length(ϵs)
    @argcheck N == length(νs)
    x₀ = νs[1]*σ*(1 - ρ^2)^(-0.5)
    xs = Vector{typeof(x₀)}(N)
    for i in 1:N
        xs[i] = (i == 1) ? x₀ : (ρ*xs[i-1] + σ*νs[i])
    end
    xs + log.(ϵs) + 1.27
end

未来工作

更复杂的模型。
解决三成分正态混合模型中的隔离问题。
每次迭代更新冲击。
优化。

参考文献

Betancourt, M. (2017)。Hamiltonian Monte Carlo 的概念介绍。
Drovandi, C. C.，Pettitt, A. N. 和 Lee, A. (2015)。使用参数辅助模型的贝叶斯间接推断。
Gallant, A. R. 和 McCulloch, R. E. (2009)。关于一般科学模型的确定及其在资产定价中的应用。
Martin, G. M.，McCabe, B. P. M.，Frazier, D. T.，Maneesoonthorn, W. 和 Robert, C. P. (2016)。状态空间模型中基于辅助似然的近似贝叶斯计算。
Rayner, G. D. 和 MacGillivray, H. L. (2002)。g-and-k 和广义 g-and-h 分布的数值最大似然估计。在：统计计算 12 57-75。
Smith, A. A. (2008)。“间接推断”。在：新帕尔格雷夫经济学词典，第二版（即将出版）。
Smith, A. A. (1993)。“使用模拟向量自回归估计非线性时间序列模型”。在应用计量经济学杂志 8.S1 中。

GSoC 2017 项目：哈密顿间接推断

2017 年 9 月 19 日 | Dorisz Albrecht