Datawhale 数学建模导论第八章学习笔记

时间序列预测与投资模型

“杂七杂八”的知识：

• 时间序列含义：按照时间排下来的序列。例如股票。
• 平稳型时间序列：数值上：均值，方差和协方差不会随着事件变化的序列；数据呈现上：1.序列的趋势线是一条水平线；2.序列不会某段分布密集，某段分布稀疏；3.序列不会某段波动很大或很小。

• 时间序列的分解模型的概念：一种用来将时间序列数据分解成不同组成部分的方法。

• 时间系列的组成：
• 长期趋势（T）：现象在较长时期内持续发展变化的一种趋向或状态，（通常）表现为一条光滑曲线趋势线。
• 季节波动（S）：由于季节的变化引起的现象发展水平的规则变动，（通常）表现为周期相对短一些的周期曲线。
• 循环波动（I）：在某段时间内，不具严格规则的周期性连续变动，（通常）表现为周期更长的周期曲线。
• 不规则波动（噪声）（C）：由于众多偶然因素对时间序列造成的影响。
• 适用场景：长达几千/上万条的时间序列。（短期的序列eg:十几条的，就可以用之前的回归模型）

• 分解模型的目的：是将时间序列拆解成上述几个成分，这样可以更好地理解时间序列数据中不同影响因素的贡献，帮助分析趋势、预测未来发展趋势或者提取季节性调整数据。

• 常见的分解方法：

1. 加法模型：加法指的是时间序分的组成是相互独立的，四个成分都有相同的量纲。
   加法模型假设时间序列是趋势、噪声、季节性和随机性的加和：[ Y_t = T_t + S_t + C_t + I_t ]

2. 乘法模型： 乘法模型输出部分和趋势项有相同的量纲，季节项和循环项是比例数，不规则变动项为独立随机变量序列，服从正态分布。
   乘法模型假设时间序列是趋势、噪声、季节性和随机性的乘积：[ Y_t = T_t \times S_t \times C_t \times I_t ]

在实际应用中，可以使用多种统计方法或者时间序列分析工具（如STL分解、Holt-Winters季节性分解等）来执行分解操作，以便于后续的分析和预测工作。

• 时间序列的建模：

1. 参数学习：通过模型参数分析序列的特征，并基于“领域知识”分析序列特点挖掘一些；
2. 预测：

• 预测周期和预测精度是一对冲突概念，如果追求高预测精度，预测周期就不能太长，因为预测时间长，难以考虑环境变化与突发事件对序列产生的影响，所做的预测也就没有意义了。
• 长期预测：只能做趋势预测；
• 短期预测

移动平均法与指数平滑法

• 适用场景：１.大趋势的建模；２.短期的趋势外推。

移动平均法

• 移动平均法的概念：用一组最近的实际数据值来预测未来时间序列的一种常用方法，是一种简单平滑预测技术。
• 适用场景：１.短期预测；２.序列需求既不快速增长／下降，且不存在季节性因素时，能有效消除预测中的随机波动。
• 移动平均法的分类（根据预测时使用的各元素的权重不同）：１.简单移动平均；２.加权移动平均。
• 基本思想：根据时间序列资料，逐项推移，依次计算包含一定项数的序列平均值，以反映长期趋势的方法。

• 二次移动平均法的一般预测公式推导

设时间集为{t|1≤t≤N, t为正整数}, 需要预测的随机序列为{Yt, 1≤t≤N}, N1(2≤N1≤N-2) 为二次移动平均法中的第一次移动平均步长, N2(2≤N2≤N-N1+1) 为二次移动平均法中的第二次移动平均步长 (2) 。则经过第一次和第二次移动平均后的序列可分别表示为:

通过式 (1) 和式 (2) 可以计算得到二次移动平均的两次平均值序列, 见表1。在表1中, 对随机序列{Yt, 1≤t≤N}而言, Yt为其在t时刻的随机变量, {Yt, 1≤t≤N}的期望值序列呈线性递增或递减态势, 则从时间节点t来看, t+T时刻的随机变量Yt+T的期望可表示为:0

其中, at为线性模型的截距项, bt为线性模型的斜率, E (Yt+T|t) 为基于t时刻的信息集下Yt+T的期望值。则t+T时刻的随机变量Yt+T可表示为:

 表1 二次移动平均法中两次平均值的计算：

注:NA表示没有相关数据。0其中, εt+T满足E (εt+T) =0, Var (εt+T) =σε2且cov (εm, εs) =0 (m≠s时) 。为对at和bt做出估计, 首先需计算Mt(1)和Mt(2), 再计算Mt(1)和Mt(2)的期望, 最后用Mt(1)和Mt(2)的估计值^Mt(1)和^Mt(2)代替各自的期望求解at和bt的估计值^at和^bt。分别计算如下:

对上式中Mt(2)两边求期望得:

 假定从随机序列{Yt, 1≤t≤N}抽取了一个样本量为n的连续等间隔观测值序列yt(1≤t≤n≤N∪N1+N2-1≤n) , 那么对E (Mt(1)) 和E (Mt(2)) 的估计值可分别表示为:

其中,

^at表示线性模型截距的样本估计值, ^bt表示线性模型斜率的样本估计值。联立式 (7) 和式 (8) , 可得: 

式 (9) 和式 (10) 是更一般的二次移动平均法的截距和斜率计算式。当N1=N2时, 可简化为其特殊形式:

进一步, 在式 (4) 的基础上, 结合式 (9) 和式 (10) , Yt+T均值的预测值可表示为:

当N1=N2时, 结合式 (11) 和式 (12) , Yt+T均值的预测值可表示为:

为了及时捕捉当期信息下截距和斜率对未来预测值的影响, 一般取T=1。此情境下, 二次移动平均法预测的均方误差的计算公式为:

通过最小化式 (15) 的取值, 即可获得二次移动平均法下的一组 (N1, N2) , 该组两次移动平均步长即为所求二次移动平均的最优移动步长。

(参考文献：二次移动平均法中两次移动平均步长关系研究  席建国）

指数平滑法

• 基本思想：将权重按照指数级进行衰减
• 指数平滑法公式：

设初始序列为yt,yt-1,yt-2,…,则t时刻的一次指数平滑值为：

其中，St(1)为一次指数平滑值，α为权重系数（也称为平滑系数），且0<α<1。由式(1)，得原始值yt,yt-1,yt-2的权重分别为α，α（1-α）,α（1-α）2。变换式（1)

则预测公式为：

二次指数平滑值为：

预测模型为：

其中，St(2)为二次指数平滑值，Ft+T为T期的预测值，T为预测超前期数，T=1,2,3……,n,at,bt为第T期的预测系数。

进一步采用三次指数平滑法进行预测，三次指数平滑值为：

预测模型为：

预测系数为：

其中，St(3)为三次指数平滑值。

（参考文献：基于三次指数平滑法对防疫医疗物资备库的研究  唐杰 ）

注意：时间序列中如果应用移动平均，预测序列的数据量会少一个窗口长度；而应用指数平滑法的时候，趋势线的长度和原始序列的长度是对齐的。

ARIMA系列模型

杂七杂八的知识点：

• ARIMA模型实际上是由多个模型组合而来，最起初的模型也都是针对平稳时间序列而言的。
• 如果原始序列非平稳但差分以后平稳也是可以变换使用这些模型的。

AR（自回归）模型

• 基本思想：用同一变数例如x的之前各期，亦即𝑥1x1​至𝑥𝑡−1xt−1​来预测本期𝑥𝑡xt​的表现，并假设它们为一线性关系；描述了当前值与历史值之间的相关关系，用变量自身的历史数据对当前数据进行预测。
• 适用场景：必须满足弱平稳性的要求，且必须具有自相关性，自相关系数小于0.5则不适用。
• 模型公式：对于p阶自回归模型（AR），其递推公式形如：

  

 这个方程的形式本质上是自己与自己的历史做回归。比如，如果这个自回归模型为三项（也可以记为AR(3)的形式），那么就是以当天序列值为因变量，前面三天的序列值为自变量构建回归模型。在构建回归方程的过程中还会引入一个白噪声项，也就是取值服从标准正态分布的一个随机序列。

注意：时间序列中自回归的思想在后面也很有用。本质上自己和自己的历史去做回归也不一定局限在线性的模型形式，也可以用多项式去做一个广义的回归，还可以用支持向量机等构建一个机器学习模型。

MA（移动平均）模型

• 模型主要特点：当前的输出（或时间序列值）被视为过去白噪声误差的加权和。
• 模型重要特征：它的自协方差函数自协方差函数和自相关系数自相关系数表现出特定的模式。
• 自协方差函数在某个滞后阶数后趋于零，表现出q阶截尾的特性，而自相关系数则表现出q阶截尾的特性。这与自回归（AR）模型形成对比，后者的自相关系数表现出拖尾的特性。
• 模型适用场景：时间序列的数据是平稳的，即它们的均值和方差保持不变，并且每个时间点的数据都是独立的。
• 在实际应用中，MA模型常与其他模型结合使用，如自回归滑动平均（ARMA）模型和自回归移动平均（ARIMA）模型，以适应更复杂的时间序列分析需求。
• 模型公式：对于q阶移动平均模型（MA），其递推公式形如：

  

 ARMA模型和ARIMA模型

 ARMA模型：

• 只需要考虑内生变量而无需考虑其他外生变量。
• 适用场景：序列是平稳序列或者差分后是平稳序列。
•  ARMA模型组成：自回归（AR）、差分（I）和移动平均（MA）三个部分。

1. 自回归模型：

   对于p阶自回归模型（AR），其递推公式形如：

   

2.  移动平均模型：

   对于q阶移动平均模型（MA），其递推公式形如：

   

3.  差分模型：

   

   注意：时间序列可能自身不平稳，但是差分一次以后可能就平稳了。差分一次可能还不平稳，但差分两次就平稳了。所以定义了差分模型。但一般来讲ARMA模型的差分次数不应该超过两次，超过两次的话应该考虑建模方法的问题了。

 ARMIA模型

• 模型基本形式：由自回归模型阶数p、差分阶数d和移动平均阶数q就可以确定ARIMA模型的基本形式，我们将其简记为ARIMA(p,d,q)

  

 对于模型最优参数选择，我们使用AIC准则法（也叫赤池信息准则）分析。

AIC准则法（也叫赤池信息准则）分析：

• AIC统计量的定义：

  

  L表示模型的最大似然函数
• 当两个模型之间存在较大差异时，差异主要体现在似然函数项，当似然函数差异不显著时，模型复杂度则起主要作用，从而参数个数少的模型是较好的选择。一般而言，当模型复杂度提高时，似然函数L也会增大，从而使AIC变小，但是参数过多时，根据奥卡姆剃刀原则，模型过于复杂从而AIC增大容易造成过拟合现象。AIC不仅要提高似然项，而且引入了惩罚项，使模型参数尽可能少，有助于降低过拟合的可能性。
• 另一个常用的准则叫贝叶斯信息准则。贝叶斯决策理论是主观贝叶斯派归纳理论的重要组成部分。是在不完全情报下，对部分未知的状态用主观概率估计，然后用贝叶斯公式对发生概率进行修正，最后再利用期望值和修正概率做出最优决策。它的公式为：

  

 实际上，相对于AIC准则，BIC准则可能更加常用一些。

灰色系统模型

• 适用场景：小中期的序列预测，并且适合有指数上升趋势的序列波动。
• 模型原理：我们先从最基本的GM(1,1)模型说起。若已知数据列𝑥(0)x(0)，进行一次累加生成新的数列。我们需要针对已知数据𝑥(0)=(𝑥(0)𝑡(1),𝑥(0)𝑡(2),...,𝑥(0)𝑡(𝑛))x(0)=(x(0)t​(1),x(0)t​(2),...,x(0)t​(n))进行累加，通过累加生成一阶累加序列：

  

 对序列均值化可以得到：

 虽然是个离散的差分模型，我们当它连续，建立灰微分方程：

 以及对应的白化微分方程：

 实际上这两个方程是等价的，通过求解白化微分方程并使用最小二乘法去拟合参数，可以得到方程的解为：

 利用最小二乘法的矩阵模式，求解这个模型，可以得到方程

然后向后差分就可以得到原数据的预测值。

另一个更一般的模型是𝐺𝑀（2,1）GM（2,1）模型。对于原始序列，得到一次累加序列𝑥1x1和一次差分序列𝑥0′x0′，然后我们可以得到：

 这就是GM（2，1）的灰色方程。白化微分方程为：

 同样去解这个方程即可。但求解过程较为复杂，所以我们以GM（1,1）为例实现灰色预测模型的建模：

import numpy as np
import math
import matplotlib.pyplot as plt
history_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]
def GM11(history_data,forcast_steps):
  n = len(history_data) # 确定历史数据体量
  X0 = np.array(history_data) # 向量化
  # 级比检验的部分可以自行补充
  lambda0=np.zeros(n-1)
  for i in range(n-1):
​    if history_data[i]:
​      lambda0[i]=history_data[i+1]/history_data[i]
​    if lambda0[i]<np.exp(-2/(n+1)) or lambda0[i]>np.exp(2/n+2):
​      print("GM11模型失效")
​      return -1
  #累加生成
  history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
  X1 = np.array(history_data_agg)
  #计算数据矩阵B和数据向量Y
  B = np.zeros([n-1,2])
  Y = np.zeros([n-1,1])
  for i in range(0,n-1):
​    B[i][0] = -0.5*(X1[i] + X1[i+1])
​    B[i][1] = 1
​    Y[i][0] = X0[i+1]
  #计算GM(1,1)微分方程的参数a和b
  A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
  a = A[0][0]
  b = A[1][0]
  #建立灰色预测模型
  XX0 = np.zeros(n)
  XX0[0] = X0[0]
  for i in range(1,n):
​    XX0[i] = (X0[0] - b/a)*(1-math.exp(a))*math.exp(-a*(i))
  #模型精度的后验差检验
  e=sum(X0-XX0)/n
  #求历史数据平均值
  aver=sum(X0)/n
  #求历史数据方差
  s12=sum((X0-aver)**2)/n
  #求残差方差
  s22=sum(((X0-XX0)-e)**2)/n
  #求后验差比值
  C = s22 / s12  
  #求小误差概率
  cobt = 0
  for i in range(0,n):
​    if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12):
​      cobt = cobt+1
​    else:
​      cobt = cobt
  P = cobt / n
  f = np.zeros(forcast_steps)
  if (C < 0.35 and P > 0.95):
​    #预测精度为一级
​    print('往后各年预测值为：')
​    for i in range(0,forcast_steps):
​      f[i] = (X0[0] - b/a)*(1-math.exp(a))*math.exp(-a*(i+n))   
​    print(f)
  else:
​    print('灰色预测法不适用')
  return f
f=GM11(history_data,20)
plt.plot(range(11,31),f)
plt.plot(range(1,11),history_data)
plt.show()

 

从图可以看到，预测值相对于真实值的波动更加光滑。但从短期序列的预测情况而言，灰色预测方法的拟合程度还是相对比较高的。

灰色关联模型

• 灰色关联分析方法概念：根据因素之间发展趋势的相似或相异程度，亦即“灰色关联度”，作为衡量因素间关联程度的一种方法。
• 基本思想：确定参考列和比较列以后需要对数列进行无量纲化处理，然后计算灰色关联系数。
• 用均值处理法，处理每个属性的数据：

  

  灰色关联系数的定义如下：

 其中ρ不超过0.5643时分辨力最好，这里为了简洁，可以取之为0.5。灰色关联度为关联系数在样本上的平均值，计算出每个属性的灰色关联度以后就可以进行分析。

模型案例：

数据如下：

 可以写出以下代码：

#导入相关库
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# 解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
#导入数据
data=pd.read_excel('huiseguanlian.xlsx')
# print(data)
#提取变量名 x1 -- x7
label_need=data.keys()[1:]
# print(label_need)
#提取上面变量名下的数据
data1=data[label_need].values
print(data1)
[m,n]=data1.shape #得到行数和列数
data2=data1.astype('float')
data3=data2
ymin=0
ymax=1
for j in range(0,n):
  d_max=max(data2[:,j])
  d_min=min(data2[:,j])
  data3[:,j]=(ymax-ymin)*(data2[:,j]-d_min)/(d_max-d_min)+ymin
print(data3)
# 绘制 x1,x4,x5,x6,x7 的折线图
t=range(2007,2014)
plt.plot(t,data3[:,0],'*-',c='red')
for i in range(4):
  plt.plot(t,data3[:,2+i],'.-')
plt.xlabel('year')
plt.legend(['x1','x4','x5','x6','x7'])
plt.title('灰色关联分析')
plt.show()
# 得到其他列和参考列相等的绝对值
for i in range(3,7):
  data3[:,i]=np.abs(data3[:,i]-data3[:,0])
#得到绝对值矩阵的全局最大值和最小值
data4=data3[:,3:7]
d_max=np.max(data2)
d_min=np.min(data2)
a=0.5 #定义分辨系数
# 计算灰色关联矩阵
data4=(d_min+a*d_max)/(data4+a*d_max)
xishu=np.mean(data4, axis=0)
print(' x4,x5,x6,x7 与 x1之间的灰色关联度分别为：')
print(xishu)

 

上图为不同的属性与x1之间的灰色关联分析图。属性x4-x7与x1之间的灰色关联系数分别为0.95294652 0.92674346 0.9004367 0.80079348，相比于后面会讲到的皮尔逊相关系数，灰色关联分析更适合时间序列的关联分析探究。

股票序列与投资组合

马科维兹均值方差模型

• 适用场景：解决最优投资组合选择问题。
• 基本思想：通过研究各资产的预期收益、方差和协方差来确定最优投资组合。
• 模型原理：
• 对于风险函数： {D \left( w\mathop{{}}\nolimits_{{1}}r\mathop{{}}\nolimits_{{b}}+w\mathop{{}}\nolimits_{{2}}r\mathop{{}}\nolimits_{{g}} \left) =w\mathop{{}}\nolimits^{{T}}{ \sum {w= \left[ \begin{array}{*{20}{l}}{ \sigma \mathop{{}}\nolimits_{{2}}^{{1}}}&{ \sigma \mathop{{}}\nolimits_{{12}}}\\{ \sigma \mathop{{}}\nolimits_{{12}}}&{ \sigma \mathop{{}}\nolimits_{{2}}^{{2}}}\end{array} \right] }}\right. \right. } 其中，𝜎21 𝜎22σ21​ σ22​分别代表波动率，𝜎212 σ212​ 代表协方差，
•  对于收益函数，若考虑在第二天就将比特币和黄金全部套现，那么套现折算减去当日购买价格和交易额，所得即为预期收益：

 那么模型形式为：

{s.t.{ \left\{ {\begin{array}{*{20}{l}}{w\mathop{{}}\nolimits_{{1}}+w\mathop{{}}\nolimits_{{w}} \le 1}\\{-B \le w\mathop{{}}\nolimits_{{1}} \le 1}\\{-G \le w\mathop{{}}\nolimits_{{2}} \le 1}\end{array}}\right. }}

这是一个多目标规划问题。为使问题简化，我们引入乘子：

{s.t.{ \left\{ {\begin{array}{*{20}{l}}{w\mathop{{}}\nolimits_{{1}}+w\mathop{{}}\nolimits_{{w}} \le 1}\\{-B \le w\mathop{{}}\nolimits_{{1}} \le 1}\\{-G \le w\mathop{{}}\nolimits_{{2}} \le 1}\end{array}}\right. }}

对这一问题进行求解即可得到对应的策略w。

最大夏普比率模型

1990年诺贝尔经济学奖得主威廉·夏普，认为当投资者建立一个风险投资组合，他们至少应该需要投资回报达到团队220718116页25无风险投资，或更多，在此基础上，他提出了夏普比率。夏普比率是可以同时考虑回报和风险的三个经典指标之一。夏普比率的目的是计算一个投资组合每单位总风险将产生多少超额回报。如果夏普比率为正，则表示基金的回报率高于波动风险；如果为负，则表示基金的操作风险大于回报率。夏普比率越高，投资组合就越好。

夏普比率是一种综合考虑风险和收益的指标，它被定义为：

 其中，𝑟𝑓rf​为无风险利率，通常取0.04作为市场估计值。此时优化问题变为：

{s.t.{ \left\{ {\begin{array}{*{20}{l}}{w\mathop{{}}\nolimits_{{1}}+w\mathop{{}}\nolimits_{{w}} \le 1}\\{-B \le w\mathop{{}}\nolimits_{{1}} \le 1}\\{-G \le w\mathop{{}}\nolimits_{{2}} \le 1}\end{array}}\right. }}

可以看到本质上这个问题还是把风险和收益二者进行了一个综合。此时求解问题是一个非有理函数，可通过数值方法得到其最优解。

风险平价模型

风险平价是由爱德华·钱博士在2005年提出的。风险平价是一种资产配置哲学，它为投资组合中的不同资产分配相等的风险权重。风险平价的本质实际上是假设各种资产的夏普比率在长期内往往是一致的，以找到投资组合的长期夏普比率的最大化。

对于比特币和黄金的风险贡献率，可以分别计算为： { \left\{ {\begin{array}{*{20}{l}}{P\mathop{{}}\nolimits_{{1}}=1-\frac{{w\mathop{{}}\nolimits_{{2}}^{{2}} \sigma \mathop{{}}\nolimits_{{2}}^{{2}}}}{{D \left( w,r \right) }}}\\{P\mathop{{}}\nolimits_{{2}}=1-\frac{{w\mathop{{}}\nolimits_{{2}}^{{1}} \sigma \mathop{{}}\nolimits_{{2}}^{{1}}}}{{D \left( w,r \right) }}}\end{array}}\right. } 为使两个序列风险尽可能一致，我们构造这样一个规划模型：

{s.t.{ \left\{ {\begin{array}{*{20}{l}}{w\mathop{{}}\nolimits_{{1}}+w\mathop{{}}\nolimits_{{w}} \le 1}\\{-B \le w\mathop{{}}\nolimits_{{1}} \le 1}\\{-G \le w\mathop{{}}\nolimits_{{2}} \le 1}\end{array}}\right. }}

最优的组合投资策略也就只需要对这个规划进行求解就可以了。

马尔可夫模型

• 适用场景：1.离散时间序列；2.连续数据的预测。
• 马尔可夫模型的概念：
• 一个马尔科夫链是离散时间的随机过程，系统的下一个状态仅仅依赖当前的所处状态，与在它之前发生的事情无关。写成表达式就是：

 马氏定理是指对于一个非周期马尔科夫链有状态转移矩阵P，有： {\mathop{{lim}}\limits_{{n \to \infty }}P\mathop{{}}\nolimits^{{n}}={ \left[ {\begin{array}{*{20}{l}}{ \pi \left( 1 \right) }&{ \pi \left( 2 \right) }&{L}&{ \pi \left( j \right) }\\{ \pi \left( 1 \right) }&{ \pi \left( 2 \right) }&{L}&{ \pi \left( j \right) }\\{M}&{M}&{L}&{M}\\{ \pi \left( 1 \right) }&{ \pi \left( 2 \right) }&{L}&{ \pi \left( j \right) }\end{array}\begin{array}{*{20}{l}}{L}\\{L}\\{L}\\{L}\end{array}} \right] }}

 而细致平稳定理是说：若对于非周期马尔可夫链，

我们就说𝜋𝑖(𝑥)πi​(x)是马尔可夫链的平稳分布。上式被称为细致平稳条件。这两个定理是马尔可夫模型中最根本的两个定理。

马尔可夫随机场就是概率无向图模型，它是一个可以用无向图表示联合概率分布。假设有一个联合概率分布P(Y)，其中Y代表一组随机变量，该联合概率分布可以由无向图来表示，图中的每一个节点表示的是Y中的一个随机变量，图中的每条边表示的两个随机变量之间的概率依赖关系，那么这个联合概率分布P(Y)怎么样才能构成一个马尔可夫随机场呢？答案是：联合概率分布P(Y)满足成对马尔可夫性、局部马尔可夫性和全局马尔可夫性这三个中的任意一个。

隐马尔可夫模型通常用一个五元组𝜆=(𝑁,𝑀,𝜋,𝐴,𝐵)λ=(N,M,π,A,B)定义：

• 模型的状态数𝑞𝑡∈{𝑆1,𝑆2,...,𝑆𝑁}qt​∈{S1​,S2​,...,SN​}。

• 模型观测值数𝑂𝑡∈{𝑉1,𝑉2,...,𝑉𝑀}Ot​∈{V1​,V2​,...,VM​}。

• 状态转移概率𝑎𝑖𝑗=𝑃(𝑞𝑡+1=𝑆𝑗∣𝑞𝑡=𝑆𝑖)aij​=P(qt+1​=Sj​∣qt​=Si​)。

• 观察概率𝑏𝑗𝑘=𝑃(𝑂𝑘=𝑉𝑘∣𝑞𝑘=𝑆𝑗)bjk​=P(Ok​=Vk​∣qk​=Sj​)。

• 初始状态概率𝜋=(𝜋1,𝜋2,...,𝜋𝑁),𝜋𝑖=𝑃(𝑞𝑖=𝑆𝑖)π=(π1​,π2​,...,πN​),πi​=P(qi​=Si​)。

在隐马尔可夫模型中存在这样几个经典问题：

• 已知模型参数计算某个序列的概率：使用前向后向算法。

• 已知参数寻找最可能产生某一序列的隐含状态：使用维特比算法。

• 已知输出寻找最可能的状态转移与概率：使用鲍姆-韦尔奇算法

• 马尔可夫模型的实现

隐马尔可夫模型（Hidden Markov Model, HMM）是一个强大的工具，用于模拟具有隐藏状态的时间序列数据。HMM广泛应用于多个领域，如语音识别、自然语言处理和生物信息学等。在处理HMM时，主要集中于三个经典问题：评估问题、解码问题和学习问题。三个问题构成了使用隐马尔可夫模型时的基础框架，使得HMM不仅能够用于模拟复杂的时间序列数据，还能够从数据中学习和预测。

1.评估问题

在隐马尔可夫模型（Hidden Markov Model, HMM）的应用中，评估问题是指确定一个给定的观测序列在特定HMM参数下的概率。简而言之，就是评估一个模型生成某个观测序列的可能性有多大。模型评估问题通常使用前向算法解决。前向算法是一个动态规划算法，它通过累积“前向概率”来计算给定观测序列的概率。前向概率定义为在时间点t观察到序列的前t个观测，并且系统处于状态i的概率。算法的核心是递推公式，它利用前一时刻的前向概率来计算当前时刻的前向概率。

import numpy as np

# 定义模型参数
states = {'Rainy': 0, 'Sunny': 1}
observations = ['walk', 'shop', 'clean']
start_probability = np.array([0.6, 0.4])
transition_probability = np.array([[0.7, 0.3], [0.4, 0.6]])
emission_probability = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])

# 观测序列，用索引表示
obs_seq = [0, 1, 2]  # 对应于 'walk', 'shop', 'clean'

# 初始化前向概率矩阵
alpha = np.zeros((len(obs_seq), len(states)))

# 初始化
alpha[0, :] = start_probability * emission_probability[:, obs_seq[0]]

# 递推计算
for t in range(1, len(obs_seq)):
    for j in range(len(states)):
        alpha[t, j] = np.dot(alpha[t-1, :], transition_probability[:, j]) * emission_probability[j, obs_seq[t]]

# 序列的总概率为最后一步的概率之和
total_prob = np.sum(alpha[-1, :])

print("Forward Probability Matrix:")
print(alpha)
print("\nTotal Probability of Observations:", total_prob)

2.解码问题

 在隐马尔可夫模型（Hidden Markov Model, HMM）中，解码问题是指给定一个观测序列和模型参数，找出最有可能产生这些观测的隐状态序列。这个问题的核心是如何从已知的观测数据中推断出隐含的状态序列，这在许多应用场景中非常有用，如语音识别、自然语言处理、生物信息学等。解决这一问题最常用的算法是维特比算法，一种动态规划方法，它通过计算并记录达到每个状态的最大概率路径，从而找到最可能的状态序列。

import numpy as np

def viterbi(obs, states, start_p, trans_p, emit_p):
    """
    Viterbi Algorithm for solving the decoding problem of HMM
    obs: 观测序列
    states: 隐状态集合
    start_p: 初始状态概率
    trans_p: 状态转移概率矩阵
    emit_p: 观测概率矩阵
    """
    V = [{}]
    path = {}

    # 初始化
    for y in states:
        V[0][y] = start_p[y] * emit_p[y][obs[0]]
        path[y] = [y]

    # 对序列从第二个观测开始进行运算
    for t in range(1, len(obs)):
        V.append({})
        newpath = {}

        for cur_state in states:
            # 选择最可能的前置状态
            (prob, state) = max((V[t-1][y0] * trans_p[y0][cur_state] * emit_p[cur_state][obs[t]], y0) for y0 in states)
            V[t][cur_state] = prob
            newpath[cur_state] = path[state] + [cur_state]

        # 不更新path
        path = newpath

    # 返回最终路径和概率
    (prob, state) = max((V[len(obs) - 1][y], y) for y in states)
    return (prob, path[state])

# 定义状态、观测序列及模型参数
states = ('Rainy', 'Sunny')
observations = ('walk', 'shop', 'clean')
start_probability = {'Rainy': 0.6, 'Sunny': 0.4}
transition_probability = {
   'Rainy' : {'Rainy': 0.7, 'Sunny': 0.3},
   'Sunny' : {'Rainy': 0.4, 'Sunny': 0.6},
}
emission_probability = {
   'Rainy' : {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
   'Sunny' : {'walk': 0.6, 'shop': 0.3, 'clean': 0.1},
}

# 应用维特比算法
result = viterbi(observations,
                 states,
                 start_probability,
                 transition_probability,
                 emission_probability)
print(result)

3.学习问题

理解隐马尔可夫模型（HMM）的模型学习问题关键在于确定模型参数，以最大化给定观测序列的出现概率。解决这一学习问题的常用方法是鲍姆-韦尔奇算法，这是一种迭代算法，通过交替执行期望步骤（E步骤）和最大化步骤（M步骤）来找到最大化观测序列概率的参数。E步骤计算隐状态的期望值，而M步骤则更新模型参数以最大化观测序列的概率。这一过程会持续重复，直至满足一定的收敛条件，如参数变化量低于特定阈值或达到预设的迭代次数。通过这种方式解决学习问题，我们可以获得一组能够很好解释给定观测数据的模型参数，这表明模型能够捕捉到观测数据中的统计规律，用于生成观测序列、预测未来观测值或识别新观测序列中的模式。

import numpy as np
from hmmlearn import hmm

# 假设我们有一组观测数据，这里我们随机生成一些数据作为示例
# 实际应用中，你应该使用真实的观测数据
n_samples = 1000
n_components = 3  # 假设我们有3个隐状态
obs_dim = 2  # 观测数据的维度，例如二维的观测空间

# 随机生成观测数据
np.random.seed(42)
obs_data = np.random.rand(n_samples, obs_dim)

# 初始化GaussianHMM模型
# 这里我们指定了n_components隐状态数量和covariance_type协方差类型
model = hmm.GaussianHMM(n_components=n_components, covariance_type='full', n_iter=100)

# 使用观测数据训练模型
# 注意：实际应用中的数据可能需要更复杂的预处理步骤
model.fit(obs_data)

# 打印学习到的模型参数
print("学习到的转移概率矩阵:")
print(model.transmat_)
print("\n学习到的均值:")
print(model.means_)
print("\n学习到的协方差:")
print(model.covars_)

深入理解隐马尔可夫模型（HMM）处理的三种经典问题——评估问题、解码问题和学习问题，可以将通过一个完整的示例来展示这些问题的应用和解决方案。如有一个简单的天气模型，其中的状态（隐藏状态）包括晴天（Sunny）和雨天（Rainy），观测（可见状态）包括人们的三种活动：散步（Walk）、购物（Shop）和清洁（Clean）。可以使用HMM来处理评估问题、解码问题和学习问题。

from hmmlearn import hmm
import numpy as np

# 定义模型参数
states = ["Rainy", "Sunny"]
n_states = len(states)

observations = ["walk", "shop", "clean"]
n_observations = len(observations)

start_probability = np.array([0.6, 0.4])

transition_probability = np.array([
    [0.7, 0.3],
    [0.4, 0.6],
])

emission_probability = np.array([
    [0.1, 0.4, 0.5],
    [0.6, 0.3, 0.1],
])

# 创建模型
model = hmm.MultinomialHMM(n_components=n_states)
model.startprob_ = start_probability
model.transmat_ = transition_probability
model.emissionprob_ = emission_probability
model.n_trials = 4

# 观测序列
obs_seq = np.array([[0], [1], [2]]).T  # 对应于观测序列 ['walk', 'shop', 'clean']

# 计算观测序列的概率
logprob = model.score(obs_seq)
print(f"Observation sequence probability: {np.exp(logprob)}")

# 继续使用上面的模型参数和观测序列

# 使用Viterbi算法找出最可能的状态序列
logprob, seq = model.decode(obs_seq, algorithm="viterbi")
print(f"Sequence of states: {', '.join(map(lambda x: states[x], seq))}")

# 假设我们只有观测序列，不知道模型参数
obs_seq = np.array([[0], [1], [2], [0], [1], [2]]).T  # 扩展的观测序列

# 初始化模型
model = hmm.MultinomialHMM(n_components=n_states, n_iter=100)
model.fit(obs_seq)

# 打印学习到的模型参数
print("Start probabilities:", model.startprob_)
print("Transition probabilities:", model.transmat_)
print("Emission probabilities:", model.emissionprob_)

• 条件随机场 

条件随机场(Conditional Random Field)是 马尔可夫随机场 + 隐状态的特例。

区别于生成式的隐马尔可夫模型，CRF是判别式的。CRF 试图对多个随机变量（代表状态序列）在给定观测序列的值之后的条件概率进行建模：

给定观测序列𝑋={𝑋1,𝑋2,...,𝑋𝑛}X={X1​,X2​,...,Xn​}，以及隐状态序列𝑌={𝑦1,𝑦2,...,𝑦𝑛}Y={y1​,y2​,...,yn​}的情况下，构建条件概率模型𝑃(𝑌∣𝑋)P(Y∣X)。若随机变量Y构成的是一个马尔科夫随机场，则 𝑃(𝑌∣𝑋)P(Y∣X)为CRF。

借助 sklearn_crfsuite 库实现

import sklearn_crfsuite

X_train = ...
y_train = ...

crf = sklearn_crfsuite.CRF(
    algorithm='lbfgs',
    c1=0.1,
    c2=0.1,
    max_iterations=100,
    all_possible_transitions=True
)
crf.fit(X_train, y_train)
# CPU times: user 32 s, sys: 108 ms, total: 32.1 s
# Wall time: 32.3 s
y_pred = crf.predict(X_test)
metrics.flat_f1_score(y_test, y_pred,
                      average='weighted', labels=labels)

本章小结

• 本章主要介绍了时间序列预测模型与投资选股模型。时间序列的核心就是如何预测，对不同体量的序列数据有不同的预测方法:

1. 小体量数据我们会用回归建模；
2. 中体量数据我们开始使用灰色系统建模；
3. 中大型体量数据我们便开始使用ARIMA系列模型建模；
4. 等到数据体量为大体量数据时，我们便会使用机器学习尤其是神经网络建模.

• 当我们能够相对准确预测未来时，我们认为还不够，因为我们最终的目的是根据我们预测的结果进行投资选股的决策，这就又是一个优化问题。在这个优化问题中，我们需要综合考虑风险和收益，想办法让收益最大风险最小，于是便有了马科维兹均值方差理论、最大夏普模型和风险平价模型。二者构成了在量化投资中的理论基石，是一个不可分的整体。