全文链接：https://tecdat.cn/?p=33885
本文描述了帮助客户使用马尔可夫链蒙特卡洛（MCMC）方法通过贝叶斯方法估计基本的单变量随机波动模型，就像Kim等人（1998年）所做的那样（点击文末“阅读原文”获取完整代码数据）。
定义模型以及从条件后验中抽取样本的函数的代码也在Python脚本中提供。
%matplotlib inline from __future__ import division ...... from src import sv
来自Kim等人（1998年）的经典单变量随机波动性模型，在此之后简称KSC，如下所示:
这里，yt代表某个资产的修正后平均收益，ht为对数波动率
示例
我们将对1981年10月1日至1985年6月28日期间的英镑/美元汇率（查看文末了解数据免费获取方式）进行建模。
ex = pd.read_excel('es.xls') dta = np.l...... .iloc[1:] endg = (dta['...... ean()) * 100
准拟然估计
估计该模型参数的一种方法是Harvey等人（1994年）的“准拟然估计法”，其中将log(ε^2_t)用与均值和方差相同的高斯随机变量来近似替换。
mod_QSV = sv.QL...... ())
贝叶斯估计
KSC提供了一种使用贝叶斯技术估计该模型的替代方法；他们将log(ε^2_t)用高斯混合分布近似表示，使得：
其中 st 是一个指示随机变量，定义为 P(st=i)=qi, i=1,…,K (K 是混合成分数目)。定义了 (qi,mi,v2i) 表示组成高斯分布的值如下所示。
# 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η 的先验分布
我们考虑共轭先验分布：
其中我们将 σr=5 和 Sσ=0.01×σr=0.05。
? 的先验分布
定义 ??=(?+1)/2，我们对 ?? 指定一个先验分布：
正如在 KSC 中讨论的那样，该先验分布在 (?1,1) 上支持随机波动性过程的平稳性。
设置 ?(1)=20 和 ?(2)=1.5 意味着 E[?]=0.86。
最后：
μ 的先验分布
KSC 建议对 μ 设定一个模糊的先验分布（或者也可以稍微具有信息的先验分布，比如 μ～N(0,10)）。
从条件后验中采样
KSC 表明，在上述指定的先验条件下，我们可以按照以下方式从条件后验中采样：
采样 σ2η
条件后验分布为：
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
采样μ?
条件后验分布为：
def draw_pos...... * (1 - phi)**2 + ...... ) return norm.r...... 2_mu**0.5)
采样htTt=1?
在混合指示符（用于生成时变观测方程矩阵）和参数条件下，可以通过通常的模拟平滑器对状态进行采样。
采样stTt=1?
每个指示变量st只能取有限个离散值（因为它是一个指示变量，表示时间t时哪个混合分布处于活动状态）。KSC表明，可以从以下概率质量函数独立地采样混合指示符：
其中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))
由于参数ση?控制潜在随机波动率过程的方差，低估将抑制样本中波动率过程的变化。如下图所示
fig, ax = plt.subplots(f...... ax.legend();
点击标题查阅往期内容
【视频】随机波动率SV模型原理和Python对标普SP500股票指数预测|数据分享
01
02
03
04
最后，我们可能对参数的完全条件后验分布感兴趣。以下是这些分布，以及后验均值和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);
