Prediction methods

1 Differential equation model

1.1 Introduction

• 特点

当描述实际对象的某些特征随时间（或空间）而演变的过程、分析它的变化规律、预测它的未来形态、研究它的控制手段时，通常需要建立对象的 动态微分方程模型。

微分模型求解的结果就是问题的答案，该答案是 唯一 的。

• 典型的模型：
  • 传染病的预测模型
  • 经济增长预测模型
  • 兰彻斯特（Lanchester）战争预测模型
  • 药物在体内的分布于排除预测模型
  • 人口的预测模型
  • 烟雾的扩散与消失模型

模型的基本规律随着时间的增长趋势呈指数形式，根据变量的个数建立微分方程。

• 优点

短、中、长期的预测都能适用，既能反应 内部规律 以及 事物的内在关系，也能分析两个因素的 相关关系，精度相应的比较高，另外对模型的改进也比较容易理解和实现。

• 缺点

虽然反应的是内部规律，但由于方程的建立是以局部规律的独立性假定为基础，故中长期预测 偏差有点大，且 解较难得到。

1.2 Case

• Problem

  硫黄岛位于东京以南 660 英里的海面上，是日军的重要空军基地。美军在1945 年2 月开始进攻，激烈的战斗持续了一个月，双方伤亡惨重，日方守军21500 人全部阵亡或被俘，美方投入兵力73000 人，伤亡20265 人，战争进行到28 天时美军宣布占领该岛，实际战斗到36 天才停止。美军的战地记录有按天统计的战斗减员和增援情况。日军没有后援，战地记录则全部遗失。

• Solution Code

  dxy=@(t,x)[-0.0544*x(2)+54000*(t>=0 & t<1)+6000*(t>=2 & t<3)+13000*(t>=5 & t<6)
      -0.0106*x(1)];  %用匿名函数定义方程右端项，这里用逻辑语句定义分段函数
  [t,xy]=ode45(dxy,[0:36],[0,21500])
  subplot(211), plot(t,xy(:,1),'r*',t,xy(:,2),'gD')
  xlabel('时间t'),  ylabel('人数'), legend('美军','日军')
  subplot(212),  plot(xy(:,1),xy(:,2))  %画微分方程组的轨线
  xlabel('美军人数x'),  ylabel('日军人数y')  

  Result:

  

  Fig. 1-1 Result of case 1

2 Grey model

2.1 Introduction

• 主要特点

  模型使用的不是原始数据序列，而是生成的数据序列，其核心体系是 灰色模型（Grey model, GM）， 即对原始数据作 累加生成 得到近似的 **指数规律 **再进行建模的方法。

• 优点

  • 不需要很多的数据，一般只需要4个数据，就能解决历史数据少、序列的完整性及可靠性低的问题；
  • 能利用微分方程来充分挖掘系统的本质，精度高；
  • 能将无规律的原始数据进行生成得到规律性较强的生成序列，预算简便，易于检验，不考虑分布规律，不考虑变化趋势。

• 缺点

  只适用于 中短期的预测，只适合 指数增长的预测。

2.2 Categories

2.2.1 GM(1, 1) forecasting model

GM(1, 1) denotes the grey model is a first order difference equation and only include one variable.

Not finished, to be continued…

2.2.2 GM(2, 1), DGM and Verhulst model

GM(1, 1) 模型适用于具有较强 指数规律 的序列，只能描述单调的变化过程，对于非单调的摆动发展序列或 有饱和的 S 形序列，可以考虑建立 GM(2, 1), DGM and Verhulst model。

• 灰色 Verhulst 预测模型

  Verhulst 模型主要用来描述具有剥和状态的过程，即 S 形过程，常用于

  • 人口预测
  • 生物生长
  • 繁殖预测
  • 产品经济寿命预测

3 Difference equation

3.1 Introduction

在利用差分方程建模研究实际问题时，常需要根据统计数据用 最小二乘法 来 **拟合 **出差分方程的 系数。其系统稳定性讨论要用到代数方程的求根。

3.1.1 一阶自回归（AR(1)）

由于时间序列一般存在自相关，故最简单的预测方法为，使用过去值 $y_{t-1}$ 来预测当前值 $y_{t}$ ，即一阶自回归模型（AR(1)）

$y_t = \beta_0 + \beta_1 y_{t-1} + \varepsilon_t \ \ (t = 2, \dots, T)$

其中，扰动项 $\varepsilon_t$ 为白噪声，故无自相关，即 $\rm{Cov}(\varepsilon_t,\ \varepsilon_s) = 0, \forall t \neq s$。假设自回归系数 $|\beta_1| < 1 $，则 $ {y_t}$ 为渐进独立的平稳过程。

3.1.2 高阶自回归（AR§）

在 AR(1) 模型中，假设扰动项为无自相关，故可用 $OLS$ 进行一致的估计。然而， 如果模型为 AR(2)，但却被误设为 AR(1)，则意味着二阶滞后项 $\beta_2 y_{t-2}$ 被纳入扰动项：

$y_t = \beta_0 + \beta_1 y_{t-1} + (\beta_2 y_{t-2} + \varepsilon_t)$

其中，由于扰动项为 $(\beta_2 y_{t-2} + \varepsilon_t) $，故扰动项与解释变量 $y_{t-1}$ 相关，因为 $\rm{Cov}(y_{t-1},\ y_{t-2}) \neq 0$。此时，扰动项存在自相关，$OLS$ 不再一致，需引入 $\beta_2 y_{t-2}$ 才能得到一致估计。

另外，从预测的角度来看，更高阶的滞后项也可能包含有用的信息。为此，更一般地，考虑 $ p$ 阶自回归模型，记为 $AR(p)$：

$y_t = \beta_0 + \beta_1 y_{t-1} + \dots + \beta_p y_{t-p} + \varepsilon_t$

其中，扰动项 $\varepsilon_t$ 为白噪声，无自相关，故 $OLS$ 为一致估计。然而，通常我们并不知道滞后期 $p$。

估计 $\hat{p}$ 的方法：

• 设一个最大滞后期 $p_{\max}$，然后令 $\hat{p} = p_{\max} $ 进行估计，并对最后一个滞后期系数的显著性进行 $t$ 检验。

  如果接受该系数为 0，则令 $\hat{p} = p_{\max} -1 $，重新进行估计，再对最后一个滞后期的系数进行 $t$ 检验。

• 使用信息准则，选择 $\hat{p}$ 使得 $AIC$ 或 $BIC$ 最小化，分别记为 $\hat{p}_{AIC}$ 与 $\hat{p}_{BIC}$。

3.1.3 自回归分布滞后模型（ADL）

在自回归 $AR(p)$ 中，为提高预测力或解释力，可引入其他解释变量，构成 自回归分布滞后模型（Autoregressive Distributed Lag Model, ADL(p, 1) or ARDL(p, q)）。

$y_t = \beta_o + \beta_1 y_{t-1} + \dots + \beta_p y_{t-p} + \gamma_1 x_{t-1} + \dots + \gamma_q x_{t-q} + \varepsilon_t$

其中，$p$ 为解释变量 $y$ 的自回归阶数，$q$ 为解释变量 $x$ 的滞后阶数。假定扰动项 $\varepsilon_t$ 为白噪声，则 $OLS$ 为一致估计。

3.1.4 误差修正模型（ECM）

ADL 是一种动态模型。从经济理论而言，相关变量之间可能存在长期的均衡关系，而变量的短期变动则是向着这个长期均衡关系的调整。ECM 正是这一思想在计量经济学中的体现。

考虑 ADL(1, 1) 模型：

$y_t = \beta_o + \beta_1 y_{t-1} + \gamma_1 x_{t-1} + \varepsilon_t$

其中，$|\beta_1| < 1$，故为平稳过程。假设经济理论 认为 $(y,\ x)$ 之间存在长期均衡关系：

$y = \phi + \theta x$

其中，$\phi$ 和 $\theta$ 为待定参数。对方程两边求期望，并令 $y^\star = \rm{E}(y_t) = \rm{E}(y_{t-1})$， $x^\star = \rm{E}(x_t) = \rm{E}(x_{t-1})$，可得：

$y^\star = \frac{\beta_0}{1-\beta_1} + \frac{\gamma_1}{1-\beta_1} x^\star$

所以，$\phi = \frac{\beta_0}{1-\beta_1}$， $\theta = \frac{\gamma_0}{1-\beta_1}$。其中 $ \theta$ 即为长期程数，衡量当 $x$ 永久性变化 1 单位时，将导致 $y$ 的永久变化幅度。在方程（5）两边同时减去 $y_{t-1}$ 得：

$\Delta y_t = (\beta_1 -1)(y_{t-1} - \phi - \theta x_{t-1}) + \varepsilon_t$

其中，$(y_{t-1} - \phi - \theta x_{t-1})$ 衡量上一期对均衡条件 $y = \phi + \theta x$ 的偏离（误差），而 $(\beta_1 -1)(y_{t-1} - \phi - \theta x_{t-1})$ 为根据上期的误差做的反向修正，成为 误差修正项。

• ECM 优点

  一般 ADL 模型都可转化为 ECM 模型。ECM 模型经济含义十分明确，而且可以分别考察长期效益（长期均衡关系）与短期效应（误差修正效应）。

4 Markov model

4.1 Defination of Markov Chain

现实生活中有很多这样的现象，某一系统在已知现在情况的条件下，系统未来时刻的情况只与现在有关，而与过去的历史无直接关系。

如，研究一个商店的累计销售额，如果现在时刻的累计销售额已知，未来某一时刻的累计销售额与现在时刻以前的任一时刻累计销售额无关。

描述这类 随机现象 的数学模型称为 马尔科夫模型，简称马氏模型

• 定义 1

设 $\{\xi_n,\ n = 1,\ 2,\dots \}$ 是一个随机序列，状态空间 $E$ 为有限或可列集，对于任意的正整数 $m,\ n $，若 $i,\ j,\ ,\ i_k \in E\ (k=1, \dots, n-1)$，有：

$P\{\xi_{n+m} = j | \xi_n = i,\ \xi_{n-1} = i_{n-1}, \dots, \xi_1 = i_1\} = P\{\xi_{n+m} = j | \xi_n = i\}$

则称 $\{\xi_{n} , n = 1,\ 2, \dots\}$ 是一个马尔可夫链。式 (9) 也被称为马氏性。

式 (9) 中，可以证明，对于 $m = 1$ 时成立，则它对 任意 的正整数 $m$ 都成立。因此，只要当 $m=1$ 式，式 (9) 成立，就可以称随机序列 $\{\xi_n,\ n=1,\ 2, \dots\}$ 具有 马氏性，即 $\{\xi_{n} , n = 1,\ 2, \dots\}$ 是一个马尔可夫链。

• 定义 2

设 $\{\xi_n,\ n = 1,\ 2,\dots \}$ 是一个马尔科夫链，如果式 (9) 右边的条件概率与 $n$ 无关，即：

$P\{\xi_{n+m} = j | \xi_n = i \} = p_{ij}(m)$

则称 $\{\xi_n,\ n = 1,\ 2,\dots \}$ 为时齐的马尔可夫链。称 $p_{ij}(m)$ 为系统由状态 $i$ 经过 $m$ 个时间间隔（或 $m$ 步）转移到状态 $j$ 的转移概率。式 (10) 称为时齐性，它的含义式系统由状态 $i$ 到 状态 $j$ 的转移概率 只依赖与时间间隔的长短，与起始时刻无关。本章介绍的马尔可夫链都是时齐的，因此省略 时齐 二字。

4.2 转移概率矩阵及柯尔莫哥洛夫定理

4.2.1 转移概率矩阵

对于马尔科夫链 $\{\xi_n,\ n = 1,\ 2,\dots \}$ ，称以 $m$ 步转移概率 $p_{ij}(m)$ 为元素的矩阵 $P (m) = p_{ij}(m)$ 为马尔科夫链的 $ m$ 步转移概率矩阵。当 $m = 1$时，记 $P(1) = P$ 称为马尔可夫链的 一步转移矩阵，或简称 转移矩阵。

转移矩阵具有如下性质：

• 对一切 $i,\ j \in E, 0 \leq p_{ij}(m) \leq 1 $，
• 对一切 $i \in E, \sum_\limits{j \in E} p_{ij}(m) = 1 $，-
• 对一切 i,\ j \in E, p_{ij}(0) = \delta_{ij} = \left\{\begin{array}{k}1,\ (i = j) \\ 0,\ (i \neq j) \end{array}\right.

当实际问题可以用马尔可夫链来描述时，首先要确定它的 状态空间及参数集合，然后确定它的 一步转移概率。该概率的确定，可以由问题的 内在规律 得到，也可以由过去的经验给出，还可以根据 预测数据 来估计。案例可见 4.4 Cases 的 Case 1。

4.2.2 柯尔莫哥洛夫定理

• 柯尔莫哥洛夫- 开普曼定理

设 $\{\xi_n,\ n = 1,\ 2,\dots \}$ 是一个马尔可夫链，其状态空间 $E = \{1,\ 2,\dots\}$，则对任意正整数 $m,\ n$，有：

$p_{ij}(n+m) = \sum_{k \in E} p_{ik}(n)p_{kj}(m)$

其中：$i,\ j \in E$。

• 定理

设 $P$ 是一步马尔科夫链转移矩阵（ $P$ 的行向量是概率向量），$P^{(0)}$ 是初始分布行向量，则第 $n$ 步的 概率分布 为：

$P^{(n)} = P^{(0)} P^n$

4.3 转移概率的渐近性质—极限概率分布

4.3.1 Introduction

• 随着 $n$ 的增大，$P^n$ 是否会趋于某以固定矩阵？

先考虑一个简单的例子：

转移矩阵 $ P = \begin{bmatrix}
0.5 & 0.5 \
0.7 & 0.3 \
\end{bmatrix}$，当 $n \rightarrow + \infty$，有：

$P^n \rightarrow \begin{bmatrix} \frac{7}{12} & \frac{5}{12}\\ \frac{7}{12} & \frac{5}{12} \\ \end{bmatrix}$

若取

\begin{align*} u = \begin{bmatrix} \frac{7}{12} & \frac{5}{12}\\ \frac{7}{12} & \frac{5}{12} \\ \end{bmatrix} \end{align*}

则 $uP = u$，$u^T$ 为矩阵 $P^T$ 的对应于特征值为 $\lambda = 1$ 的特征（概率）向量，$u$ 也称为 $P$ 的不动点向量。

• 哪些转移矩阵具有不动点向量？

  **定义：**马尔可夫链的转移矩阵 $P$ 是正则 $\Leftrightarrow \ \exist\ k \in N^+$，使 $P^k$ 任一元素都为正

  **定理：**若 $P$ 是一个马尔科夫链的正则阵，则：

  • $P$ 有 唯一 的不动点向量 $W,\ W$ 的每个分量都为正；
  • $P$ 的 $n$ 次幂 $P^n$（ $n$ 为正整数）随 $n$ 的增加趋于矩阵 $\overline{W},\ \overline{W}$ 的每一行向量均等于不动点向量 $W$。

设时齐马尔科夫链的状态空间为 $E$，如果对于所有 $i,\ j \in E$，转移概率 $p_{ij}(n)$ 存在极限：

$\lim_{n\to \infty}p_{ij}(n) = \pi_j\quad (\text{不依赖于 }i)$

或：

$P(n) = P^n \underset{(n\to \infty)}{\longrightarrow} \begin{bmatrix} \pi_1 & \pi_2 & \dots & \pi_j & \dots \\ \pi_1 & \pi_2 & \dots & \pi_j & \dots \\ \vdots & \vdots & \ddots & \vdots & \ddots \\ \pi_1 & \pi_2 & \ddots & \pi_j & \dots \\ \vdots & \vdots & \ddots & \vdots & \ddots \\ \end{bmatrix}$

则称此链具有 遍历性。若 \sum_\limits j \pi_j = 1，则同时称 $\vec {\pi} = [\pi_1,\ \pi_2, \dots]$ 为链的 极限分布。

• 就有限链的遍历性给出一个充分条件：

  **定理：**设时齐马尔可夫链 $\{\xi_n,\ n = 1,\ 2,\dots \}$ 的状态空间为 $E = \{a_1,\dots,\ a_N\},\ P = (p_{ij})$ 是它的一步转移矩阵，如果存在正整数 $m$，使对任意的 $a_i,\ a_j \in E$，都有：

  $p_{ij}(m) > 0,\ i, j = 1, 2, \dots, N$

  则此链具有 遍历性；且有极限分布 $\vec \pi = [\pi_i, \dots, \pi_N]$，他是方程组：

  $\pi = \pi P\quad \mathrm{or\quad} \pi_j = \sum_{i=1}^N \pi_i p_{ij}, \ j = 1, \dots,\ N$

  的满足条件 \pi_j > 0,\ \sum_\limits{j=1}^N \pi_j = 1 的唯一解。

4.4 Cases

• Case 1:

某计算机机房的一台计算机经常出故障，研究者每隔15 分钟观察一次计算
机的运行状态，收集了24 小时的数据（共作97 次观察）。用1 表示正常状态，用0 表示不正常状态，所得的数据序列如下：

1110010011111110011110111111001111111110001101101
111011011010111101110111101111110011011111100111

**解：**设 $X_n (n = 1,\dots,97)$ 为第 $n$ 个时段的计算机状态，可以认为它是一个时齐马氏链，状态空间 $E = \{0,\ 1\}$，编写如下 $Matlab$ 程序：

a1='1110010011111110011110111111001111111110001101101';
a2='111011011010111101110111101111110011011111100111';
a=[a1 a2];
f00=length(findstr('00',a))
f01=length(findstr('01',a))
f10=length(findstr('10',a))
f11=length(findstr('11',a))

Result:

求得 96 次的转移情况是：

$\left\{ \begin{array}{c} 0 \rightarrow 0,\ \ 8 \text{次}; \qquad\ 0 \rightarrow 1,\ 18 \text{次} \\ 1 \rightarrow 0,\ 18 \text{次}; \qquad 1 \rightarrow 1,\ 52 \text{次} \end{array} \right.$

因此，一步转移概率可用频率近似地表示为：

$\left\{ \begin{array}{c} P_{00} = P\{X_{n+1} = 0 | X_n = 0 \} \approx \frac{8}{8 + 18} = \frac{4}{13} \\ P_{01} = P\{X_{n+1} = 1 | X_n = 0 \} \approx \frac{18}{8 + 18} = \frac{9}{13} \\ P_{10} = P\{X_{n+1} = 0 | X_n = 1 \} \approx \frac{18}{18 + 52} = \frac{9}{35} \\ P_{11} = P\{X_{n+1} = 1 | X_n = 1 \} \approx \frac{52}{18 + 52} = \frac{26}{35} \end{array} \right.$

• Case 2:

若顾客的购买是无记忆的，即已知现在顾客购买情况，未来顾客的购买情况不受过去购买历史的影响，而只与现在购买情况有关。现在市场上供应 $A,\ B,\ C$ 三个不同厂家生产的 50 克袋装味精，用 $\xi_n = 1,\ \xi_n, \ \xi_n = 3$ 分别表示”顾客第 $n$ 次购买 $A,\ B,\ C$ 厂的味精”。显然，$\{\xi_n, n = 1,\ 2,\dots\}$ 是一个马氏链。若已知第一次顾客购买三个厂味精的概率依次为 $0.2,\ 0.4,\ 0.4$。又知道一般顾客购买的倾向由表2给出。求

1）顾客第四次购买各家味精的概率；

2）预测经过长期的多次购买滞后，顾客的购买倾向如何？

	A	B	C
A	0.8	0.1	0.1
B	0.5	0.1	0.4
C	0.5	0.3	0.2

解：

1）第一次购买的概率分布为：$P^{(1)} = [0.2,\ 0.4,\ 0.4]$，一步状态转移矩阵为：

$P = \begin{bmatrix}0.8 & 0.1 & 0.1 \\ 0.5 & 0.1 & 0.4 \\ 0.5 & 0.3 & 0.2\\\end{bmatrix}$，则顾客第四次购买各家味精的概率为：

$P^{(4)} = P^{(1)} P^3 = [0.7004, 0.136, 0.1636]$

2）这个马尔可夫链的转移矩阵 $P$ 满足 $p_{ij}(m) > 0,\ i,\ j = 1, 2,\dots, N$，可以求出其极限概率分布。为此，解下列方程组：

$\left\{ \begin{array}{l} p_1 = 0.2p_1 + 0.8p_2 + 0.1p_3,\\ p_2 = 0.8p_1 + 0.3p_3, \\ p_3 = 0.2p_2 + 0.6p_3, \\ p_1 + p_2 + p_3 = 1, \end{array} \right.$

求得 $p_1 = \frac{5}{7},\ p_2 = \frac{11}{84},\ p_3 = \frac{13}{84}$。这说明，无论第一次顾客购买的情况如何，经过长期多次购买后，$A$ 厂产的味精占有市场的 $\frac{5}{7}$，$B,\ C$ 两厂的产品分别占有市场的 $\frac{11}{84}$ 和 $\frac{13}{84}$。

Relative codes:

format rat
p=[0.8 0.1 0.1;0.5 0.1 0.4;0.5 0.3 0.2];
a=[p'-eye(3);ones(1,3)];
b=[zeros(3,1);1];
p_limit=a\b

或者利用求转移矩阵 $P$ 的转置矩阵 $P^T$ 的特征值 1 对应的特征(概率)向量，求得极
限概率。编写程序如下：

clc,clear
p=[0.8 0.1 0.1;0.5 0.1 0.4;0.5 0.3 0.2];
p=sym(p');
[x,y]=eig(p)
y=diag(y);y=double(y);
ind=find(y==max(y));
p=x(:,ind)/sum(x(:,ind))