全文链接:https://tecdat.cn/?p=33885

原文出处:拓端数据部落公众号

本文描述了帮助客户使用马尔可夫链蒙特卡洛(MCMC)方法通过贝叶斯方法估计基本的单变量随机波动模型,就像Kim等人(1998年)所做的那样。

定义模型以及从条件后验中抽取样本的函数的代码也在Python脚本中提供。

%matplotlib inline
from __future__ import division
......

from src import sv

来自Kim等人(1998年)的经典单变量随机波动性模型,在此之后简称KSC,如下所示:

image.png

这里,yt代表某个资产的修正后平均收益,ht为对数波动率

示例

我们将对1981年10月1日至1985年6月28日期间的英镑/美元汇率进行建模。

image.png

ex = pd.read_excel('es.xls')
dta = np.l......
.iloc[1:]

endg = (dta['......
ean()) * 100

准拟然估计

估计该模型参数的一种方法是Harvey等人(1994年)的“准拟然估计法”,其中将log(ε^2_t)用与均值和方差相同的高斯随机变量来近似替换。

mod_QSV = sv.QL......
())

image.png

贝叶斯估计

KSC提供了一种使用贝叶斯技术估计该模型的替代方法;他们将log(ε^2_t)用高斯混合分布近似表示,使得:

image.png
其中 st 是一个指示随机变量,定义为 P(st=i)=qi, i=1,…,K (K 是混合成分数目)。定义了 (qi,mi,v2i) 表示组成高斯分布的值如下所示。

image.png

# q_i, m_i, v_i^2
ksc_aras = np.array([......
)

在给定 stTt=1 的条件下,每个时间段的观测方程是由一个高斯噪声项定义的。

通过设置 K=7 是对 logε2t 进行很好近似的方法,Omori et al. (2007) 将其扩展到 K=10。

class TLDT(sm.t......
Model):
    """
    时变局部线性确定性趋势
  ......

        # 转换为对数平方,带有偏移量
        endog = n.logenog**2+ offset

        # 初始化基本模型
        super(TVLLDT, self)._......
tationary')

        # 设置观测方程的时变数组
        self['o......
.nobs))

        # 设置状态空间矩阵的固定分量
        self['d......
0] = 1

    def update......
7036, v_i^2)
        self['o......
rams[1]
        self['state_cov', 0, 0] = params[2]

先验分布

为了计算模型,我们需要为参数 θ 的先验分布进行特定的指定。下面的先验规范取自于 KSC。

σ2η 的先验分布

我们考虑共轭先验分布:

image.png
其中我们将 σr=5 和 Sσ=0.01×σr=0.05。

ϕ 的先验分布

定义 ϕ∗=(ϕ+1)/2,我们对 ϕ∗ 指定一个先验分布:

image.png
正如在 KSC 中讨论的那样,该先验分布在 (−1,1) 上支持随机波动性过程的平稳性。

设置 ϕ(1)=20 和 ϕ(2)=1.5 意味着 E[ϕ]=0.86。

最后:

image.png

μ 的先验分布

KSC 建议对 μ 设定一个模糊的先验分布(或者也可以稍微具有信息的先验分布,比如 μ∼N(0,10))。

从条件后验中采样

KSC 表明,在上述指定的先验条件下,我们可以按照以下方式从条件后验中采样:

采样 σ2η

条件后验分布为:

image.png

def draw_po......
or_params=(5, 0.05)):
    sigma_r, S_sigma = prior_params

    v1 = sig......
i * (states[0, :-1] - mu))**2)
    delta1 = S_sigma + tmp1 + tmp

    return ingammars(v1,scal=deta1)
采样 ϕ

我们可以应用 Metropolis 步骤:从 N(ϕ^,Vθ) 中生成一个提议值 ϕ∗

def g(phi, ......


    # 先验分布对非平稳过程给予零权重
    if np.abs(phi) >= 1:
        ret......
2) / 2 * sigma2
    tmp2 = 0.5 * np.log(1 - phi**2)

    return n......

    V_phi = sigma2 / tmp2

    proposal ......
om.uniform() else phi
采样μ̂

条件后验分布为:

image.png

def draw_pos......
 * (1 - phi)**2 + ......
)

    return norm.r......
2_mu**0.5)
采样htTt=1̂

在混合指示符(用于生成时变观测方程矩阵)和参数条件下,可以通过通常的模拟平滑器对状态进行采样。

采样stTt=1̂

每个指示变量st只能取有限个离散值(因为它是一个指示变量,表示时间t时哪个混合分布处于活动状态)。KSC表明,可以从以下概率质量函数独立地采样混合指示符:

image.png

其中fN(y∗t∣a,b)表示均值为a,方差为b的高斯随机变量在y∗t处的概率密度。

def (mod states):
    resid = od.nog[:, 0] - states[0]

    # 构建均值 (nobs x 7), 方差 (7,), 先验概率 (7,)
    means = ks_aram......
0]

    # 调整维度以便广播计算
    resid = np.repe......
[None, :], mod.nobs,
                                    axis=0)

    # 计算对数似然 (nobs x 7)
    loglikelihoods = -0.5 * ((resi......
* variances))

    # 得到(与后验(对数))成比例的值 (nobs x 7)
    posterior_kernel = log......
ilities)

    # 归一化得到实际后验概率
    tmp = logsumxp(psterir_kernl,axis=1)
    posterior_probabilitie......
d, states)

    # 从后验中抽取样本
    varaes = np.radom.niorm(ize=od.obs)
......
    sample = np.argmax(tmp, axis=1)

    return sample

MCMC

下面我们进行10,000次迭代以从后验中进行抽样。在下面展示结果时,我们将舍弃前5,000次迭代作为燃烧期,并且在剩下的5,000次迭代中,我们只保存每十次迭代的结果。然后从剩下的500次迭代中计算结果。


# 设置模型和模拟平滑器
md = TVLLT(eog)
mo.(0, sothr_stateTrue)
sim = md.siutin_sother()

# 模拟参数
nitertons = 10000
brn = 5000
tin = 10

# 存储轨迹
trae_sooted = np.eros((_iteations+ 1 mod.nobs))......

trce_sim2 = np.ers((n_iteations + 1, 1))

# 初始值 (p. 367)
trce_miing[0] = 0
[0] = 0.95
trace_sigma2[0] = 0.5
# 迭代
for s in range(1, n_teations + 1):
    # 更新模型参数
    mod.updat_ming(tace_mixing[s-1])......
    # 模拟平滑
    sim.smuate()......


    # 抽取混合指标
    trac_miing[s] = drawmixngmod states)
    
    # 抽取参数
    tra_phi[s] = (mod, sates, trace_phi[s-1], trace_mu[s-1], trace_sigma2[s-1])......

结果

下面我们给出参数的后验均值。我们还展示了相应的QMLE估计值。这些估计值与 ϕ 和 β 的后验均值相似,但是对于 ση² 的QMLE估计值约为贝叶斯方法的一半,可能表明准拟然方法的一个缺点。

# 参数的后验均值
menphi = n.men(trae_hi[burn:thin])......

print('  beta          = %.5f' % npexp(rs_LSVparams[2] / 2))

image.png

由于参数ση²控制潜在随机波动率过程的方差,低估将抑制样本中波动率过程的变化。如下图所示

fig, ax = plt.subplots(f......

ax.legend();

output_26_0.png

最后,我们可能对参数的完全条件后验分布感兴趣。以下是这些分布,以及后验均值和QMLE估计值。

fig, axes = plt.subplots(1, 3, ......

axes[0].set(title=r'$\phi$', ylim=ylim)
axes[0].legend(loc='upper left')
......
axes[2].set(title=r'$\beta$', ylim=ylim);

output_28_0.png


rate reset presents banks with irrbb challenge_1365585930_standard_1536x1536.webp

最受欢迎的见解

1.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%bb%98%e5%88%b6%e7%94...)HAR-RV-J与递归神经网络(RNN)混合模型预测和交易大型股票指数的高频波动率

2.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%94%9f%e5%ad%98%e5%88...)WinBUGS对多元随机波动率模型:贝叶斯估计与模型比较

3.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e5%a6%82%e4%bd%95%e5%9c...)波动率的实现:ARCH模型与HAR-RV模型

4.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e4%b8%ad%e4%bd%bf%e7%94...)R语言ARMA-EGARCH模型、集成预测算法对SPX实际波动率进行预测

5.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%94%9f%e5%ad%98%e5%88...)使用R语言随机波动模型SV处理时间序列中的随机波动率

6.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80ggplot2%e8%af%af%e5%b7%a...)R语言多元COPULA GARCH 模型时间序列预测

7.[](https://tecdat.cn/r-%e8%af%ad%e8%a8%80%e7%bb%98%e5%88%b6%e5%8...)R语言基于ARMA-GARCH过程的VAR拟合和预测

8.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e5%a6%82%e4%bd%95%e6%89...)R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型

9.R语言对S&P500股票指数进行ARIMA + GARCH交易策略


拓端tecdat
195 声望46 粉丝