灰色预测模型GM python 灰色预测模型gm(1,1) 过程

灰色预测模型

  灰色预测是对既含有已知信息又含有不确定信息的系统进行预测，就是对在一定范围内变化的、与时间有关的灰色过程进行预测。
  灰色预测对原始数据进行生成处理来寻找系统变动的规律，并生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

---

目录

• 灰色预测模型
• 一、GM(1,1)模型简介
• 二、GM(1,1)原理
• 三、准指数规律的检验
• 四、GM(1,1)模型的评价
• 五、模型扩展（★）

---

一、GM(1,1)模型简介

  GM(1,1)是最简单的灰色预测模型，它是使用原始的离散非负数据列，通过一次累加生成削弱随机性的较有规律的新的离散数据列，然后通过建立微分方程模型，得到在离散点处的解经过累减生成的原始数据的近似估计值，从而预测原始数据的后续发展。（本文中只探究 GM(1,1) 模型，第一个 1 表示微分方程是一阶的，后面的 1 表示只有一个变量）

---

二、GM(1,1)原理

$x^{(0)}=(x^{(0)}(1)，x^{(0)}(2),\cdots,x^{(0)}(n))$ 是最初非负数据列，对其一次累加得到新的生成数据 $x^{(1)}$

$x^{(1)}=(x^{(1)}(1)，x^{(1)}(2),\cdots,x^{(1)}(n))$，其中：$x^{(1)}(m)=\sum_{i=1}^{m}{x^{(0)}(i)},m=1,2,\cdots,n$。

$z^{(1)}$ 为数列 $x^{(1)}$ 的紧邻生成数列，即 $z^{(1)}=(z^{(1)}(1)，z^{(1)}(2),\cdots,z^{(1)}(n))$，其中: $z^{(1)}(m)=\delta x^{(1)}(m)+(1-\delta)x^{(1)}(m-1),m=2,3,\cdots$,$n$ 且 $\delta =0.5$。

  

     

---

  我们称方程

$x^{(0)}(k)+az^{(1)}(k)=b$ 为 GM(1,1) 模型的基本形式（k=2,3,…,n），其中，

$b$ 表示灰作用量。

$-a$ 表示发展系数。下面引入矩阵形式：

$u=(a,b)^T,\quad Y=\left[ \begin{array}{c} x^{(0)}(2)\\ x^{(0)}(3)\\ \vdots \\ x^{(0)}(n) \end{array}\right],\quad B=\left[ \begin{array}{c} -z^{(1)}(2)\quad 1\\ -z^{(1)}(3)\quad 1\\ \quad\vdots \quad\qquad \vdots \\ -z^{(1)}(n)\quad 1 \end{array}\right]$

$x^{(0)}(k)+az^{(1)}(k)=b$ 可表示为：$\bf Y=Bu$
  利用最小二乘法得到参数 a,b 的估计值为：
$\color{fuchsia}{ \hat u=\left( \begin{array}{c} \hat a\\ \hat b\\ \end{array} \right)=\bf (B^TB)^{-1}B^TY}$

  实际上就是将 $x^{(0)}$ 序列视为因变量 $y$, $z^{(1)}$ 序列是为自变量 $x$，进行回归。
$\color{red}{\bf x^{(0)}(k)=-az^{(1)}(k)+b \quad\Rightarrow\quad y=kx+b}$

---

引入最小二乘法(OLS)

  最小二乘法定义：
$\hat{y_i} =kx_i+b,\quad \hat{k},\hat{b}=\mathop{arg\min}\limits_{k,b}(\sum_{i=1}^{n} (y_i-\hat{y_i})^2)= \mathop{arg\min}\limits_{k,b} \sum_{i=1}^{n}{(y_i-kx_i-b)^2}$

  我们令：
$L=\sum_{i=1}^{n}(y_i-kx_i-b)^2=[y_1-kx_1-b,y_2-kx_2-b,\cdots,y_n-kx_n-b] \left[ \begin{array}{c} y_1-kx_1-b\\ y_2-kx_2-b\\ \vdots\\ y_n-kx_n-b\\ \end{array} \right]$

  并且令：
$矩阵 Y = \left[ \begin{array}{c} y_1\\ y_2\\ \vdots\\ y_n\\ \end{array} \right],\quad X=\left[ \begin{array}{c} 1 \quad x_1\\ 1 \quad x_2\\ \vdots \text{ }\quad\text{ }\vdots\\ 1 \quad x_n\\ \end{array} \right],\quad \beta= \left[ \begin{array}{c} b\\ k\\ \end{array} \right]$

  那么：
$X\beta=\left[ \begin{array}{c} b+kx_1\\ b+kx_2\\ \vdots\\ b+kx_n\\ \end{array} \right],\quad Y-X\beta=\left[ \begin{array}{c} y_1-kx_1-b\\ y_2-kx_2-b\\ \vdots\\ y_n-kx_n-b\\ \end{array} \right]$

  所以有：
$\begin{array}{l} L=(Y-X\beta)^T(Y-X\beta) \\ \quad=(Y^T-\beta^TX^T)(Y-X\beta)\\ \quad=Y^TY-Y^TX\beta-\beta^TX^TY+\beta^TX^TX\beta \\ \end{array}\\ 则 \text{ }\hat \beta=\left[ \begin{array}{c} \hat b\\ \hat k\\ \end{array} \right]=\mathop{arg\min}\limits_{\beta}{(L)}=\mathop{argmin}\limits_{\beta}=(Y^TY-Y^TX\beta-\beta^TX^TY+\beta^TX^TX\beta)\\$

  对矩阵求导：
$\frac {dL}{dt}=-X^TY-X^TY+2X^TX\beta=0\Rightarrow X^TX\beta=X^TY\\$

  所以：
$\color{fuchsia}{\bf\hat\beta=(X^TX)^{-1}X^TY}$

---

利用 OLS 估计的回归结果可以得出 $\hat a$ 和 $\hat b$，即 $x^{(0)}{(k)}=-\hat az^{(1)}(k)+\hat b \quad (k=2,3,\cdots,n)$

$x^{(0)}{(k)}=-\hat az^{(1)}(k)+\hat b \Rightarrow \color {blue} { x^{(1)}(k)-x^{(1)}(k-1) \color{black}{\text{ } = \text{ }-\hat a } \color{fuchsia}{ z^{(1)}(k)}} \color{black}{\text{ }+\text{ }} \color{red}{\hat b}$

  对于以上式子:
$等式左边：\color {blue} { x^{(1)}(k)-x^{(1)}(k-1)}\color{#000}{ = \int_{k-1}^k {\frac{dx^{(1)}(t)}{dt}} \,dt}\quad(牛顿莱布尼茨公式)\\ 等式右边：\color {fuchsia} {z^{(1)}(k)} \color {black} {\text{ } =\text{ } \frac{x^{(1)}(k)+x^{(1)}(k-1)}{2}}\approx \int_{k-1}^k x^{(1)}(t)dt\quad(定积分几何意义)$

$\color{blue}{\int_{k-1}^k {\frac{dx^{(1)}(t)}{dt}} \,dt}\color{black}{\text{ }\approx\text{ }-\hat a}\color{fuchsia}{\int_{k-1}^k x^{(1)}(t)dt}\color{black}{\text{ }+\text{ }}\color{red}{\int_{k-1}^k \hat b dt}\color{black}{=\int_{k-1}^k\left[-\hat ax^{(1)}(t)+\hat b\right]}dt$

$微分方程：\bf\frac{dx^{(1)}(t)}{dt}=-\hat ax^{(1)}(t)+\hat b\quad被称为\text{ }GM(1,1)\text{ }模型的白化方程$

$\bf x^{(0)}(k)+az^{(1)}(k)=b被称为\text{ }GM(1,1)\text{ }的灰色微分方程$

---

白化方程求解：

$\frac{dx^{(1)}(t)}{dt}=-\hat ax^{(1)}(t)+\hat b$，如果取初始值 $\hat x^{(1)}(t)\mid_{t=1}=x^{(0)}(1)$，可求出对应的解为：
$x^{(1)}(t)=\left[ x^{(0)}(1)-\frac{\hat b}{\hat a}\right]e^{-\hat a(t-1)}+\frac{\hat b}{\hat a}\\ 所以\quad x^{(1)}(m+1)=\left[ x^{(0)}(1)-\frac{\hat b}{\hat a}\right]e^{-\hat a(m)}+\frac{\hat b}{\hat a},\quad m=1,2,\cdots,n-1$

$x^{(1)}(m)=\displaystyle\sum_{i=1}^{m}x^{(0)}(i)，m=1,2,\cdots,n$，所以可以得到：
$\bf x^{(0)}(m+1)=x^{(1)}(m+1)-x^{(1)}(m)=(1-e^{\hat a})\left[ x^{(0)}(1)-\frac{\hat b}{\hat a}\right]e^{-\hat a(m)}+\frac{\hat b}{\hat a},\quad m=1,2,\cdots,n-1$

$m≥n$

---

$f(x)=C_1e^{C_2(x-1)}$，其中这里的指数规律阵对 $x^{(1)}(k)$

---

三、准指数规律的检验

1. 数据具有准指数规律是使用灰色系统建模的理论基础。
2. 累加 $r$ 次的序列为 $x^{(r)}=(x^{(r)}(1)，x^{(r)}(2),\cdots,x^{(r)}(n))$，定义级比 $\sigma(k)=\frac{x^{(r)}(k)}{x^{(r)}(k-1)},k=2,3,\cdots,n$。
3. 如果 $\forall k,\sigma(k) \in[a,b]$，且区间长度 $\delta=b-a＜0.5$，则称累加 $r$
4. 具体到 GM(1,1) 模型中，我们只需判断累加一次后的序列 $x^{(1)}=(x^{(1)}(1),x^{(1)}(2),\cdots,x^{(1)}(n))$
5. 根据上述公式：
   序列 $x^{(1)}$ 的级比 $\sigma(k)=\frac{x^{(1)}(k)}{x^{(1)}(k-1)}=\frac{x^{(0)}(k)+x^{(1)}(k-1)}{x^{(1)}(k-1)}=\frac{x^{(0)}(k)}{x^{(1)}(k-1)}+1$
   定义 $\rho(k)=\frac{x^{(0)}(k)}{x^{(1)}(k-1)}$ 为原始序列 $x^{(0)}$ 的光滑比，注意到： $\rho(k)=\frac{x^{(0)}(k)}{x^{(0)}(1)+x^{(0)}(2)+\cdots+x^{(0)}(k-1)}$
   假设 $x^{(0)}$ 为非负序列，那么随着 $k$ 增加，最终 $\rho(k)$ 会逐渐接近 0，因此要是的具有 $x^{(1)}$ 具有准指数规律，即 $\forall k$，区间长度 $\delta$ ＜0.5，只需要保证 $\rho(k)∈(0,0.5)$，此时序列 $x^{(1)}$ 的级比 $\sigma(k)∈(1,1.5)$

---

注意：一般前两期：$\rho(2)$ 和 $\rho(3)$

---

四、GM(1,1)模型的评价

  使用 GM(1,1) 模型对未来的数据进行预测时，首先需要检验 GM(1,1) 模型对原数据的拟合程度（对原始数据的还原效果）。一般有两种方法：残差检验和级比偏差检验。

---

残差检验：

$\varepsilon(k)=x^{(0)}(k)-\hat x^{(0)}(k)，k=2,3,\cdots,n$
  相对残差： $\varepsilon_r(k)=\frac{|x^{(0)}-\hat x^{(0)}(k)|}{ x^{(0)}(k)}\times100\%，k=2,3,\cdots,n$

  平均相对残差：
$\overline \varepsilon_r = \frac{1}{n-1} \sum_{k=2}^{n}|\varepsilon_r(k)|$

• 如果 $\overline \varepsilon_r＜ 20\%$，则认为 GM(1,1) 对原始数据的拟合达到了一般要求。
• 如果 $\overline \varepsilon_r＜ 10\%$，则认为 GM(1,1) 对原始数据的拟合达效果非常不错。（10%不绝对）

---

级比偏差检验：

$x^{(0)}(k-1)$ 和$x^{(0)}(k)$ 计算出原始数据的级比 $\sigma(k)$：
$\sigma(k)=\frac{x^{(0)}(k)}{x^{(0)}(k-1)},\quad(k=2,3,\cdots,n)$

$(-\hat a)$ 计算出相应的级比偏差和平均级比偏差:
$\eta(k)=\left|1-\frac{1-0.5\hat a}{1+0.5\hat a}\frac{1}{\sigma(k)}\right|,\quad \bar \eta=\sum_{k=2}^{n}\frac{\eta(k)}{n-1}$

• 如果 $\bar \eta＜ 0.2$，则认为 GM(1,1) 对原始数据的拟合达到了一般要求。
• 如果 $\bar \eta＜ 0.1$，则认为 GM(1,1) 对原始数据的拟合达效果非常不错。

---

$\eta(k)$ 越小，说明 $x^{(0)}(k)$ 和 $\hat x^{(0)}(k)$ 越接近。特别地，当 $\eta(k)=0$ 时，可以得到：
$\sigma(k)=\frac{1-0.5\hat a}{1+0.5\hat a}=\frac{\hat x^{(0)}(k)}{x^{(0)}(k-1)}$

---

五、模型扩展（★）

1. 数据是以年份度量的非负数据（如果是月份或者季度数据就要用 时间序列模型);
2. 数据能经过准指数规律的检验（除了前两期外，后面至少90%的期数的光滑比要低于0.5，规定的 90% 不绝对) ;
3. 数据的期数较短且和其他数据之间的关联性不强（小于等于10，也不能太短了，比如只有 3 期数据），要是数据期数较长，一般用传统的时间序列模型比较合适。
4. 在传统的 GM(1,1) 模型的基础上，每预测一次，将预测的数据作为已知数据进行下一次预测，那么这种模型为 新信息 GM(1,1) 模型。在新信息 GM(1,1) 模型的基础上，去掉最老信息 $x^{(0)}(1)$ ，那么这种模型为 新陈代谢 GM(1,1) 模型。应对比传统GM(1,1)、新信息 GM(1,1) 、新陈代谢 GM(1,1) 三种模型的预测效果，抉择使用。

---