R语言中的马尔可夫区制转移(Markov regime switching)模型

金融分析师通常关心市场何时“发生变化”:几个月或者几年内市场的典型行为可以立即转变为非常不同的行为。

由Kaizong Ye,Sherry Deng撰写

金融分析师通常关心检测市场何时“发生变化”:几个月或至几年内市场的典型行为可以立即转变为非常不同的行为。

投资者希望及时发现这些变化,以便可以相应地调整其策略,但是这样做可能很困难。

×

1.1 马尔可夫过程

       马尔可夫过程(Markov process)是一类随机过程。由俄国数学家A.A.马尔可夫于1907年提出。该过程具有如下特性:在已知目前状态(现在)的条件下,它未来的演变(将来)不依赖于它以往的演变 (过去 )。例如森林中动物头数的变化构成——马尔可夫过程。在现实世界中,有很多过程都是马尔可夫过程,如液体中微粒所作的布朗运动、传染病受感染的人数、车站的候车人数等,都可视为马尔可夫过程。(这里虽然我也不清楚这些现象到底是不是,姑且就认为是吧!

马尔科夫过程中最核心的几个概念:过去,现在,将来。其中最核心的在于“现在”如何理解。

在马尔可夫性的定义中,"现在"是指固定的时刻,但实际问题中常需把马尔可夫性中的“现在”这个时刻概念推广为停时(见随机过程)。例如考察从圆心出发的平面上的布朗运动,如果要研究首次到达圆周的时刻 τ以前的事件和以后的事件的条件独立性,这里τ为停时,并且认为τ是“现在”。如果把“现在”推广为停时情形的“现在”,在已知“现在”的条件下,“将来”与“过去”无关,这种特性就叫强马尔可夫性。具有这种性质的马尔可夫过程叫强马尔可夫过程。在相当一段时间内,不少人认为马尔可夫过程必然是强马尔可夫过程。首次提出对强马尔可夫性需要严格证明的是J.L.杜布。直到1956年,才有人找到马尔可夫过程不是强马尔可夫过程的例子。马尔可夫过程理论的进一步发展表明,强马尔可夫过程才是马尔可夫过程真正研究的对象。

这段话实在是太过于抽象了,不好理解,心里有数就行,因为这里的过去、现在、将来和我们生活中是有所差别的,不太好理解!

所以:一个马尔科夫过程就是指过程中的每个状态的转移只依赖于之前的 n个状态,这个过程被称为 n阶马尔科夫模型,其中 n是影响转移状态的数目。最简单的马尔科夫过程就是一阶过程,每一个状态的转移只依赖于其之前的那一个状态,这也是后面很多模型的讨论基础,很多时候马尔科夫链、隐马尔可夫模型都是只讨论一阶模型,甚至很多文章就将一阶模型称之为马尔科夫模型,现在我们知道一阶只是一种特例而已了

对于一阶马尔科夫模型,则有:

如果第 i 时刻上的取值依赖于且仅依赖于第 i−1 时刻的取值,即


从这个式子可以看出,xi 仅仅与 xi-1有关,二跟他前面的都没有关系了,这就是一阶过程。
 

总结:马尔科夫过程指的是一个状态不断演变的过程,对其进行建模后称之为马尔科夫模型,在一定程度上,马尔科夫过程和马尔科夫链可以打等号的。

1.2 马尔科夫性(无后效型)

在马尔科夫过程中,在给定当前知识或信息的情况下,过去(即当前以前的历史状态)对于预测将来(即当前以后的未来状态)是无关的。这种性质叫做无后效性。简单地说就是将来与过去无关,值与现在有关,不断向前形成这样一个过程。

1.3 马尔可夫链

时间和状态都是离散的马尔可夫过程称为马尔可夫链,简记为Xn=X(n),n=0,1,2…马尔可夫链是随机变量X1,X2,X3…的一个数列。

这种离散的情况其实草是我们所讨论的重点,很多时候我们就直接说这样的离散情况就是一个马尔科夫模型。

(1)关键概念——状态空间

马尔可夫链是随机变量X1,X2,X3…Xn所组成的一个数列,每一个变量Xi 都有几种不同的可能取值,即他们所有可能取值的集合,被称为“状态空间”,而Xn的值则是在时间n的状态。

(2)关键概念——转移概率(Transition Probability)

马尔可夫链可以用条件概率模型来描述。我们把在前一时刻某取值下当前时刻取值的条件概率称作转移概率。

上面是一个条件概率,表示在前一个状态为s的条件下,当前状态为t的概率是多少。

(3)关键概念——转移概率矩阵

很明显,由于在每一个不同的时刻状态不止一种,所以由前一个时刻的状态转移到当前的某一个状态有几种情况,那么所有的条件概率会组成一个矩阵,这个矩阵就称之为“转移概率矩阵”。

一、什么是隐马尔可夫模型HMM(Hidden Markov Model, HMM)

1.1 从一个简单的例子说起

既然说它是在马尔科夫模型的基础之上发展来的,而这必然有关系了,先看一个简单的例子来描述这个情景。

在某些情况下马尔科夫过程不足以描述我们希望发现的模式。回到之前那个天气的例子,一个隐居的人可能不能直观的观察到天气的情况,但是有一些海藻。民间的传说告诉我们海藻的状态在某种概率上是和天气的情况相关的。在这种情况下我们有两个状态集合,一个可以观察到的状态集合(海藻的状态)和一个隐藏的状态(天气的状况)。我们希望能找到一个算法可以根据海藻的状况和马尔科夫假设来预测天气的状况。

其中,隐藏状态的数目和可以观察到的状态的数目可能是不一样的。在语音识别中,一个简单的发言也许只需要80个语素来描述,但是一个内部的发音机制可以产生不到80或者超过80种不同的声音。同理,在一个有三种状态的天气系统(sunny、cloudy、rainy)中,也许可以观察到四种潮湿程度的海藻(dry、dryish、damp、soggy)。在此情况下,可以观察到的状态序列和隐藏的状态序列是概率相关的。于是我们可以将这种类型的过程建模为一个隐藏的马尔科夫过程一个和这个马尔科夫过程概率相关的并且可以观察到的状态集合

看了这个例子就明白了,实际上就说如果某一个状态我们是没办法直接观察到的,但是我们可以通过观察与这个状态紧密相关的另一个我们可以观察的状态来推导那个观察不到的状态,就是这样一个过程。

2.2 HMM较为正式的定义

隐马尔可夫模型(hidden Markov model, HMM)是另一种常用的概率模型。在该模型中,观测值(可观测状态)是依据一定的概率由某些不可见的内部状态(隐状态)决定,而这些状态之间服从某种马尔科夫模型。

下面给出几个关键概念:

(1)隐状态,观测状态

(2)发射概率(emission probability)

上面的可观测状态 yi 是观测到的数据,它的取值根据条件概率


得到,这个表示在隐状态Zi的条件之下得到yi的概率.这一概率称作发射概率(emission probability)

同时,由于隐状态 zi 满足马尔科夫转移概率,即

1.3 HMM 定义的数学表达

前面都是一些文字性的描述,那么使用数学语言怎么描述HMM呢?

对于HMM模型,首先我们假设,Q是所有可能的隐藏状态的集合,,V是所有可能的观测状态的集合,即:

其中,N是可能的隐藏状态数,M是所有的可能的观察状态数。

对于一个长度为T的序列,I 对应的状态序列, O是对应的观察序列,即:

 

其中,任意一个隐藏状态 it∈Q,任意一个观察状态ot∈V.

 

二、隐马尔可夫模型的核心知识点

2.1 HMM的两个假设

(1) 齐次马尔科夫链假设——本质上就是隐藏状态是一个马尔科夫链

即任意时刻的隐藏状态只依赖于它前一个隐藏状态。当然这样假设有点极端,因为很多时候我们的某一个隐藏状态不仅仅只依赖于前一个隐藏状态,可能是前两个或者是前三个。但是这样假设的好处就是模型简单,便于求解。如果在时刻 t 的隐藏状态是,在时刻t+1的隐藏状态是, 则从时刻 t到时刻 t+1的HMM状态转移概率 aij 可以表示为:

这样 aij 可以组成马尔科夫链的状态转移矩阵A:

注意:这里的“状态转移矩阵A”实际上就是马尔科夫链的“概率转移矩阵”,它的大小是NxN的方阵,N表示的是隐藏状态的可能取值数目。

(2)观测独立性假设。

即任意时刻的观察状态只仅仅依赖于当前时刻的隐藏状态,这也是一个为了简化模型的假设。如果在时刻 t 的隐藏状态是 it=qj, 而对应的观察状态为 ot=vk, 则该时刻观察状态 vk在隐藏状态 qj下生成的概率为 bj(k),满足:

这样 bj(k) 可以组成观测状态生成的概率矩阵B:

注意:这里的意思是,在某一个时刻,观测状态仅仅与这一时刻的隐藏状态有关,而与之前的观测状态是无关的,实际上也就是说观测序列并没有组成一个马尔科夫链哦!

所以这里的矩阵B称之为观测值转移矩阵。它是一个N行M列的矩阵,N表示隐藏状态的取值种类数,M表示观测值的状态取值种类数。

除此之外,我们需要一组在时刻t=1的隐藏状态概率分布 π:

注意:这里的初始概率分布是一个N维向量,N表示隐藏状态的取值种类数。

2.2 HMM模型的三要素

一个HMM模型,可以由隐藏状态初始概率分布π, 状态转移概率矩阵A观测状态概率矩阵B决定,这称之为隐马尔科夫模型的三要素。π,A决定隐藏状态序列,B决定观测序列。因此,HMM模型可以由一个三元组λ表示如下:



RHmmCRAN不再可用,因此我想使用其他软件包复制功能实现马尔可夫区制转换(Markov regime switching)模型从而对典型的市场行为进行预测,并且增加模型中对参数的线性约束功能。

library(SIT)
load.packages('quantmod')
y=returns
ResFit = HMM(y, nStates=2)
DimObs = 1


课程

R语言数据分析挖掘必知必会

从数据获取和清理开始,有目的的进行探索性分析与可视化。让数据从生涩的资料,摇身成为有温度的故事。

立即参加

matplot(fb$Gamma, type='l', main='Smoothed Probabilities', ylab='Probability')
		legend(x='topright', c('State1','State2'),  fill=1:2, bty='n')
 fm2 = fit(mod, verbose = FALSE)

使用logLik在迭代69处收敛:125.6168

	#*****************************************************************
	# #添加一些数据,看看模型是否能够识别状态
	#****************************************************************** 
	bear2  = rnorm( 100, -0.01, 0.20 )
	bull3 = rnorm( 100, 0.10, 0.10 )
	bear3  = rnorm( 100, -0.01, 0.25 )
	true.states = c(true.states, rep(2,100),rep(1,100),rep(2,100))
	y = c( bull1, bear,  bull2, bear2, bull3, bear3 )
 

DimObs = 1

 
	plot(data, type='h', x.highlight=T)
		plota.legend('Returns + Detected Regimes')
#*****************************************************************
# 加载历史价格
#****************************************************************** 
data = env()
getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
 
price = Cl(data$SPY)
	open = Op(data$SPY)
ret = diff(log(price))
	ret = log(price) - log(open)
 
atr = ATR(HLC(data$SPY))[,'atr']
 
fm2 = fit(mod, verbose = FALSE)

使用logLik在迭代30处收敛:18358.98

print(summary(fm2))
Initial state probabilties model pr1 pr2 pr3 pr4 0 0 1 0
 
Transition matrix toS1 toS2 toS3 toS4 fromS1 9.821940e-01 1.629595e-02 1.510069e-03 8.514403e-45 fromS2 1.167011e-02 9.790209e-01 8.775478e-68 9.308946e-03 fromS3 3.266616e-03 8.586650e-47 9.967334e-01 1.350529e-69 fromS4 3.608394e-65 1.047516e-02 1.922545e-130 9.895248e-01
 
Response parameters Resp 1 : gaussian Resp 2 : gaussian Re1.(Intercept) Re1.sd Re2.(Intercept) Re2.sd St1 2.897594e-04 0.006285514 1.1647547 0.1181514 St2 -6.980187e-05 0.008186433 1.6554049 0.1871963 St3 2.134584e-04 0.005694483 0.4537498 0.1564576 St4 -4.459161e-04 0.015419207 2.7558362 0.7297283
 
 	Re1.(Intercept)	Re1.sd	Re2.(Intercept)	Re2.sd
St1	0.000289759401378951	0.00628551404616354	1.16475474419891	0.118151350440916
St2	-6.98018749098021e-05	0.00818643307634358	1.65540488736983	0.187196307284941
St3	0.000213458358141314	0.00569448330115608	0.453749781945066	0.156457606460757
St4	-0.00044591612667264	0.0154192070819596	2.75583620018895	0.72972830143278
probs = posterior(fm2)
 
print(head(probs))
rownames(x)	state	S1	S2	S3	S4
1	3	0	0	1	0
2	3	0	0	1	0
3	3	0	0	1	0
4	3	0	0	1	0
5	3	0	0	1	0
6	3	0	0	1	0
layout(1:3)
plota(temp, type='l', col='darkred')
	plota.legend('Market Regimes', 'darkred')
layout(1:4)


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds