一.背景
1.1 数据源
现有2000年到2019年的数据，其中y列表示收入，x1~x13表示与收入相关联的13个特征值。
1.2 处理目的
a.分析、识别影响收入y的关键属性；
b.预测2020年、2021年的收入；
二. 分析过程
2.1 读取数据
import pandas as pd
import numpy as np
import sys
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn')#设置样式
plt.rcParams['font.sans-serif']=['SimHei']# 支持中文
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
import warnings
warnings.filterwarnings('ignore')
data = pd.read_csv('data/data.csv')
2.2 相关性分析
使用heatmap来画图呈现数据相关性。使用mask为得出的相关性矩阵corr画遮罩。
corr = data.corr(method='pearson')
corr
plt.figure(figsize=(12,10))
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True #cmap是设置热图的颜色
cmap = sns.diverging_palette(220, 10, as_cmap=True)
#绘制热图
g = sns.heatmap(corr, mask=mask, cmap=cmap, square=True, annot=True, fmt='0.2f')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
#plt.savefig("corr.png")
plt.show()
最终得出的相关性矩阵和相关关系图如下。
从上面的相关性呈现来看，可以得出以下结论：
a. 特征值x11与收入y的关系不显著，呈现负相关。
b. 其余属性与收入呈现高度正相关。
c. x2,x3,x7与多条属性存在共线性。
本次分析目的在于找出与收入y显著相关的特征值，然而上述结果显示，数十个指标与收入y都呈现显著性相关，这些属性之间存在着信息的重复，因此接下来需要运用Lasso回归来对上述特征值做进一步筛选。
2.3 Lasso回归-选取相关变量
使用Lasso回归，处理共线性，有效降维。
from sklearn.linear_model import Lasso
lasso = Lasso(1000) # 调用Lasso()函数，设置λ的值为1000
lasso.fit(data.iloc[:,0:13],data['y'])
tar = lasso.coef_ != 0 # 返回一个相关系数是否为零的布尔数组
data_a = data.iloc[:,0:13]
new_reg_data = data_a.iloc[:, tar ] # 返回相关系数非零的数据
new_reg_data 即为Lasso回归选中的新的显著性特征变量X。下图为Lasso回归选出的特征值。
对于x1，x3，x4，x5，x6，x7，x8，x13这些新选取的相关性特征值，y收入值与他们分别的相关性系数如下。
2.4 构建灰色预测模型
白色系统：指一个系统的的内部特征是完全可以认知的，即系统的信息完全公开。
黑色系统：指一个系统的内部信息对外界来说是一无所知的，只能通过它与外界的联系来加以观测研究。
灰色系统：一部分的信息是已知的，另一部分是未知的，系统内各因素间有不确定的关系，GM（1，1）就是灰色系统，也是最常用的模型系统。
灰色预测模型（Gray Forecast Model）是通过少量的、不完全的信息，建立数学模型并做出预测的一种预测方法。对单调性曲线有较好的预测。
定义一个灰色预测函数
def GM11(x0):
import numpy as np
x1 = x0.cumsum() #1-AGO序列
z1 = (x1[:len(x1)-1] + x1[1:])/2.0 #紧邻均值（MEAN）生成序列
z1 = z1.reshape((len(z1),1))
B = np.append(-z1, np.ones_like(z1), axis = 1)
Yn = x0[1:].reshape((len(x0)-1, 1))
[[a],[b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Yn) #计算参数
f = lambda k: (x0[0]-b/a)*np.exp(-a*(k-1))-(x0[0]-b/a)*np.exp(-a*(k-2)) #还原值
delta = np.abs(x0 - np.array([f(i) for i in range(1,len(x0)+1)]))
C = delta.std()/x0.std()
P = 1.0*(np.abs(delta - delta.mean()) < 0.6745*x0.std()).sum()/len(x0)
return f, a, b, x0[0], C, P #返回灰色预测函数、a、b、首项、方差比、小残差概率
new_reg_data 为上一个步骤选取的特征值。给接下来需要预测的两个年份填上空值。
new_reg_data.index = range(2000,2020)
#给接下来两个年份填上空值
new_reg_data.loc[2020] = None
new_reg_data.loc[2021] = None
l = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']
for i in l:
f = GM11(new_reg_data.loc[range(2000,2020),i].values)[0]
new_reg_data.loc[2020,i] = f(len(new_reg_data)-1) # 2020年预测结果
new_reg_data.loc[2021,i] = f(len(new_reg_data)) # 2021年预测结果
new_reg_data[i] = new_reg_data[i].round(2) # 保留两位小数
y = list(data['y'].values) # 提取财政收入列，合并至新数据框中
y.extend([np.nan,np.nan])
new_reg_data['y'] = y
new_reg_data_GM11 = new_reg_data
new_reg_data_GM11 为初步预测的y值。结果如下图。
2.5 构建支持向量机回归模型预测
调取模型，选取特征值与Y值。
from sklearn.svm import LinearSVR
data_train = new_reg_data_GM11.loc[range(2000,2020)].copy()
#数据标准化
data_mean = data_train.mean()
data_std = data_train.std()
data_train = (data_train-data_mean)/data_std
#选取测试集数据集
feature =['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']
x_train = data_train[feature].values
y_train = data_train['y'].values
linearsvr.predict(x) ，这是y是预测后归一化的值；
这里x,y都做过数据的归一化处理，根据公式 y归一值=(y归一前的值-y.mean())/y.std() ,因此有 y归一前的值=y归一值*y.std()+y.mean()
#调用LinearSVR()
linearsvr = LinearSVR()
linearsvr.fit(x_train,y_train)
# 给所有数据做归一化处理，刚才训练集选取的是2000-2019年数据，现在是全部数据集2000-2021年
x = ((new_reg_data_GM11[feature] - data_mean[feature])/data_std[feature]).values
new_reg_data_GM11['y_pred'] = linearsvr.predict(x) * data_std['y'] + data_mean['y']
最终预测的结果表格为new_reg_data_GM11
2.6 比较y-实际值与y-预测值的差异
画出预测结果图
fig = new_reg_data_GM11[['y','y_pred']].plot(subplots = True, style=['b-o','r-*']) # 画出预测结果图
plt.show()
结果如下
