首先，我们先导入需要的库和数据，并进行数据清洗和平滑处理。此处我们采用了简单指数平滑法进行平滑处理。
```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# 导入数据
data = pd.read_csv('data.csv', index_col=0)
# 数据清洗
data.dropna(inplace=True)
# 简单指数平滑法平滑处理
alpha = 0.8
data['close'] = data['close'].ewm(alpha=alpha).mean()
# 数据标准化
data = (data - data.mean()) / data.std()
# 绘制数据图
plt.plot(data)
plt.show()
```
接下来，我们对数据进行灰色马尔可夫链建模，并计算模型参数。
```python
# 灰色马尔可夫链建模
def GM11(x0):
x1 = np.cumsum(x0)
z1 = (x1[:-1] + x1[1:]) / 2.0
B = np.append(-z1.reshape(-1, 1), np.ones_like(z1).reshape(-1, 1), axis=1)
Y = x0[1:].reshape(-1, 1)
[[a], [b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)
X = np.zeros_like(x0)
X[0] = x0[0]
for i in range(1, len(x0)):
X[i] = (x0[0] - b/a) * np.exp(-a*(i-1)) - (x0[0] - b/a) * np.exp(-a*i)
return X
# 计算模型参数
X0 = data['close'].values
X1 = np.array([GM11(X0[i:i+5]) for i in range(len(X0)-4)])
P = np.zeros((len(X1), len(X1)))
for i in range(len(X1)):
for j in range(len(X1)):
if i >= j:
P[i][j] = np.sum(X1[i] == X1[j]) / len(X1[i])
```
接下来，我们对模型预测的结果进行检验，包括残差检查、关联度检验和后验差检验。
```python
# 残差检查
e = np.abs(X0[4:] - X1[:, 4])
plt.plot(e)
plt.show()
# 关联度检验
r = np.corrcoef(data['close'].values[4:], X1[:, 4])[0][1]
print('关联度：', r)
# 后验差检验
delta = np.abs(X0[4:] - np.dot(P, X0[:-4]))
C = delta.std() / X0.std()
P_value = 1.0 - 2.0 / (len(X0) - 1)
if C < 0.35 and P_value > 0.05:
print('后验差比值：{:.2f}，P值：{:.2f}，模型精度等级：好'.format(C, P_value))
elif C < 0.5 and P_value > 0.05:
print('后验差比值：{:.2f}，P值：{:.2f}，模型精度等级：合格'.format(C, P_value))
else:
print('后验差比值：{:.2f}，P值：{:.2f}，模型精度等级：不合格'.format(C, P_value))
```
接下来，我们根据模型的预测结果划分系统状态，并检验所得序列是否具有马氏性。
```python
# 划分系统状态
s = np.zeros(len(X0))
s[0] = 1
for i in range(1, len(X0)):
if X0[i] > X1[:, i-1].max():
s[i] = np.argmin(X1[:, i-1]) + 2
else:
s[i] = np.argmin(X1[:, i-1]) + 1
# 检验序列是否具有马氏性
N = 3
T = len(X0)
F = np.zeros((N, N))
for i in range(1, T):
F[int(s[i-1]-1)][int(s[i]-1)] += 1
for i in range(N):
if sum(F[i]) != 0:
F[i] /= sum(F[i])
print('状态转移概率矩阵：\n', F)
```
接下来，我们对灰色马尔科夫链模型进行预测，得到未来的状态概率分布和预测值。
```python
# 预测未来的状态概率分布
T0 = len(X0)
future_n = 5
future_s = np.zeros((future_n, T0+future_n))
future_s[:, 0] = [i+1 for i in range(future_n)]
for i in range(future_n):
for j in range(T0, T0+i):
future_s[i][j+1-T0] = np.argmin(X1[:, j-1]) + 1
for i in range(T0+1, T0+future_n):
F[int(future_s[:, i-T0-1]-1)][int(future_s[:, i-T0]-1)] += 1
for i in range(N):
if sum(F[i]) != 0:
F[i] /= sum(F[i])
print('未来5天的状态概率分布：\n', F)
# 预测未来的值
future_X = np.zeros(future_n)
for i in range(future_n):
if future_s[i][-1] == 1:
future_X[i] = X1[:, -1].min()
else:
future_X[i] = X1[:, -1][future_s[i][-1]-2]
future_X = future_X * data['close'].std() + data['close'].mean()
print('未来5天的预测值：\n', future_X)
```
接下来，我们用加权灰色马尔可夫链模型进行建模，并计算加权灰色马尔可夫链的状态转移概率矩阵，对加权灰色马尔科夫链模型进行预测，得到未来的预测值。
```python
# 加权灰色马尔可夫链建模
def WGM11(x0, weight):
x1 = np.cumsum(x0)
z1 = (x1[:-1] + x1[1:]) / 2.0
B = np.append(-(z1*weight).reshape(-1, 1), weight.reshape(-1, 1), axis=1)
Y = x0[1:].reshape(-1, 1)
[[a], [b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)
X = np.zeros_like(x0)
X[0] = x0[0]
for i in range(1, len(x0)):
X[i] = (x0[0] - b/a) * np.exp(-a*(i-1)) - (x0[0] - b/a) * np.exp(-a*i)
return X
# 计算加权灰色马尔可夫链的状态转移概率矩阵
X0 = data['close'].values
T = len(X0)
N = 3
W = np.zeros((N, T))
W[:, 0] = [1, 0, 0]
for i in range(1, T):
if X0[i] > X1[:, i-1].max():
s = np.argmin(X1[:, i-1])
for j in range(N):
if j == s:
W[j][i] = 0.5
else:
W[j][i] = 0.25
else:
s = np.argmin(X1[:, i-1])
for j in range(N):
if j == s:
W[j][i] = 0.75
else:
W[j][i] = 0.125
P = np.zeros((N, N))
for i in range(N):
for j in range(N):
P[i][j] = np.sum(W[i][1:] * (s[1:] == j+1)) / np.sum(W[i][1:])
print('加权灰色马尔可夫链的状态转移概率矩阵：\n', P)
# 预测未来的状态概率分布
future_n = 5
future_W = np.zeros((N, T+future_n))
future_W[:, 0] = [1, 0, 0]
for i in range(1, T+future_n):
if i <= T:
if X0[i] > X1[:, i-1].max():
s = np.argmin(X1[:, i-1])
for j in range(N):
if j == s:
future_W[j][i] = 0.5
else:
future_W[j][i] = 0.25
else:
s = np.argmin(X1[:, i-1])
for j in range(N):
if j == s:
future_W[j][i] = 0.75
else:
future_W[j][i] = 0.125
else:
for j in range(N):
future_W[j][i] = np.sum(future_W[:, i-T-1] * P[:, j])
print('未来5天的状态概率分布：\n', future_W[:, -future_n:])
# 预测未来的值
future_X = np.zeros(future_n)
for i in range(future_n):
if future_W[:, -1][0] > future_W[:, -1][1]:
future_X[i] = X1[0, -1]
else:
future_X[i] = X1[1, -1]
future_X = future_X * data['close'].std() + data['close'].mean()
print('未来5天的预测值：\n', future_X)
```
最后，我们可视化以上所有的预测结果。
```python
# 可视化预测结果
plt.plot(np.arange(len(X0)), X0, label='real')
plt.plot(np.arange(len(X0)-4)+4, X1[-1], label='predict')
plt.plot(np.arange(len(X0), len(X0)+future_n), future_X, label='future')
plt.legend()
plt.show()
```
扫一扫
6446
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
