在这一部分，我们尝试使用一些机器学习模型，来进行股票收益率的预测。机器学习研究如何从数据中产生“模型”的算法，即“学习算法”。有了学习算法，我们把经验数据提供给它，它就能基于这些数据产生模型，在面对新的情况时，模型会给我们提供相应的判断。在我们要讨论的股票预测的例子中，我们根据历史数据研究一个“模型”，这个模型在得到了今天的开盘价、收盘价等信息之后，能够预测明天股票收盘价是涨还是跌。

我们使用的是有监督学习，即根据有标签的历史数据来训练模型，这里的标签是根据要预测的股票价格的变动给出的。当我们要预测明天的对数收益率时，预测的是连续数值，可以归为一个回归问题，此时可以用我们之前介绍过的均方根误差作为预测是否准确的评价标准。当我们要预测明天收盘价的涨跌时，预测的是离散值，可以归为一个分类问题，此时可以用准确率或精度（Accuracy），即分类正确的样本数占样本总数的比例，作为评价指标。其他的评价指标还包括查准率，查全率，F1分数等。

在这一章节，我们会尝试scikit-learn中包含的线性回归、k近邻、随机森林模型和支持向量机，以及XGBoost模型，尝试预测对数收益率和涨跌。

1. 准备工作

在这一部分，除了使用之前已经安装和使用过的scikit-learn之外，我们还需要安装XGBoost。XGBoost最初是一个研究项目，由当时在华盛顿大学的陈天奇负责。近几年，由于受到在许多机器学习竞赛中多次获奖，XGBoost逐渐成为机器学习领域具有较强实用性的算法工具选择之一。XGBoost在Windows上的安装比较麻烦，一种比较简单的方法是先下载XGBoost对应自己电脑的Python和Windows版本的.whl文件，然后在命令行中利用pip命令安装wheel文件。首先到这个网站下载wheel文件，然后找到下载的文件所在的路径，运行pip命令：pip install xgboost-0.90-***.whl。XGBoost在Mac上的安装可以参考这篇帖子。

安装完成之后，我们导入这一章节需要用到的packages：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
%config InlineBackend.figure_format = 'svg'
import warnings
warnings.filterwarnings('ignore')  # 忽略一些warnings
import sklearn
print(sklearn.__version__)
import xgboost
print(xgboost.__version__)

这里我使用的scikit-learn版本为0.20.3，XGBoost版本为0.90。与上一章节类似，我们载入数据，并把数据分成训练、验证和测试集：

def split_data():
    df = pd.read_csv('1.csv', dtype={'trade_date': str})
    df['trade_date'] = pd.to_datetime(df['trade_date'], format='%Y%m%d')
    df = df.set_index('trade_date', drop=False)  # 指明drop=False保留trade_dt这一列
    df = df.sort_index()  # 按照时间顺序排序

    # 计算log return
    df['close1'] = df['close'].shift(-1)  # 1天之后的收盘价
    df['log_ret'] = np.log(df['close1'] / df['close'])  # 1天之后的log return

    df = df.dropna()

    train = df['2010-01-01':'2017-12-31']  # 我们使用2010-2017年的数据作为训练集
    cv = df['2018-01-01': '2018-12-31']  # 我们使用2018年的数据作为验证集
    train_cv = df['2010-01-01':'2018-12-31']  # 训练集和验证集
    test = df['2019-01-01': '2019-12-31']  # 我们使用2019年的数据作为测试集

    return df, train, cv, train_cv, test

df, train, cv, train_cv, test = split_data()
print('Train dataset shape: %s' % str(train.shape))
print('Validation dataset shape: %s' % str(cv.shape))
print('Test dataset shape: %s' % str(test.shape))

输出为：

Train dataset shape: (1874, 13)
Validation dataset shape: (243, 13)
Test dataset shape: (243, 13)

定义均方根误差的评价指标：

def get_rmse(y_true, y_pred):
    '''
    Compute root mean squared error (RMSE)
    '''
    # 把输入的变量转成numpy的array格式
    y_true, y_pred = np.asarray(y_true), np.asarray(y_pred)
    return np.sqrt(((y_true - y_pred) ** 2).mean())

回顾上一章节中的对数收益率的预测结果：

oneday_rmse = 0.028567
mov_avg_rmse = 0.020553
arima_rmse = 0.020049
print('one-day prediction:\t %f' % oneday_rmse)
print('mov_avg prediction:\t %f' % mov_avg_rmse)
print('arima prediction:\t %f' % arima_rmse)

输出为：

one-day prediction:     0.028567
mov_avg prediction:     0.020553
arima prediction:     0.020049

2. 预测对数收益率

首先我们继续尝试预测对数收益率。

2.1 线性回归

我们首先使用scikit-learn中的Linear Regression，熟悉一下scikit-learn的使用。利用一个长度为N的滑动时间窗内的历史对数收益率，对下一个对数收益率进行线性回归。给定数据集，其中，线性回归试图学得一个线性模型以尽可能准确地预测目标值。在这里，我们的线性模型的输入变量为前天的对数收益率，输出目标为明天的对数收益率，我们要做的事情就是拟合出一条误差最小的直线，然后预测这个直线对应的下一个点，如下图所示：

我们继续使用上一章节定义的通用函数：

def get_preds_rolling(df, target_col, N, pred_func, offset):
    # 我们在一个DataFrame df上应用rolling函数，例如train_cv
    # target_col是我们要预测的column，例如这里是close
    # N是时间窗的大小，也就是说我们会使用t-1,t-2, t-N天的数据
    # pred_func是我们真正使用的预测函数，基于滑动窗内的数据预测下一个数据点，之后会给出不同的定义
    # offset是返回结果时，我们提取出真正要预测的部分，例如输入train_cv时，我们要提取cv的部分，即offset=train.shape[0]

    # 直接调用rolling得到的结果是t, t-1, t-(N-1)的结果
    pred_list = df[target_col].rolling(window = N).apply(pred_func)

    # 因此我们将预测结果整体往后位移一位，最前面补上一个np.nan
    pred_list = np.concatenate((np.array([np.nan]), np.array(pred_list[:-1])))

    # 我们返回的结果是对应offset之后到末尾的部分
    return pred_list[offset:]

线性回归在scikit-learn中对应LinearRegression：

# 导入我们要使用的函数
from sklearn.linear_model import LinearRegression

# 定义我们的特定的预测函数pred_func
def pred_func(x):
    X_train = np.array(range(len(x)))  # X_train对应横坐标上：[0,1,2,3……，len(x)-1] 这里我们假设数据是按交易日等间隔分布的
    y_train = np.array(x)  # y_train对应纵坐标上具体要回归的数值
    X_train = X_train.reshape(-1, 1)  # 需要转成[[0], [1], [2], ……]的格式，下同
    y_train = y_train.reshape(-1, 1)
    regr = LinearRegression(fit_intercept=True)  # 定义模型，设置拟合截距项
    regr.fit(X_train, y_train)  # 拟合模型
    pred = regr.predict(np.array([len(x)]).reshape(-1, 1))  # 预测下一个数值，即横坐标对应len(x)-1
    return pred[0][0]

# 进行一个简单的试验，让大家了解一下pred_func的工作工程
x = [2, 3, 4, 5]
print(pred_func(x))

输出结果为：

6.000000000000001

然后我们在验证集上对参数N进行搜索：

cv_rmse_list = []
target_col = 'log_ret'
offset = train.shape[0]
for N in np.arange(2, 31):
    y_pred = get_preds_rolling(train_cv, target_col, N, pred_func, offset)
    y_true = cv[target_col].values
    # 评估RMSE
    cv_rmse_list.append(get_rmse(y_true, y_pred))

画出RMSE随N变化的曲线：

plt.figure()
plt.plot(np.arange(2, 31), cv_rmse_list, 'x-')
plt.grid()
plt.xlabel('N')
plt.ylabel('RMSE')
plt.xlim([2, 30])
plt.xticks([2, 5, 10, 15, 20]);

结果如图所示：

根据验证集上的表现，我们选取N=30：

target_col = 'log_ret'
# 对测试集使用我们选出的最好的N
N = 30
offset = train_cv.shape[0]
# 这次应用在整个df上进行滑动预测
y_pred = get_preds_rolling(df, target_col, N, pred_func, offset)
y_true = test[target_col].values
# 评估RMSE
lin_reg_rmse = get_rmse(y_true, y_pred)
print('Results of linear regression prediction')
print('RMSE: %f' % lin_reg_rmse)

输出为：

Results of linear regression prediction
RMSE: 0.021567

跟之前的最好的ARIMA结果0.020049相比，线性回归并没有在这个特定数据集上取得提高，我们接下来继续尝试使用更多的输入特征和更复杂的方法。

2.2 构造输入特征

2.2.1 构造基本的输入特征

我们用过去5天的open,high,low,close,vol和amount作为基本的输入特征。大家也可以尝试对vol和amount进行log变换后作为输入特征，这也是一种常见的数据预处理步骤。

basic_features = ['open', 'high', 'low', 'close', 'vol', 'amount']  # 原始数据中的特征
input_features = basic_features[:]  # 我们要输入给模型的特征
lagged_days = 5
lagged_df = df.copy()  # 复制一份df
for feature in basic_features:
    for day in range(1, lagged_days):
        col = feature+'_lag_'+str(day)  # 我们把前几天的数值作为新的特征
        lagged_df[col] = lagged_df[feature].shift(day)
        input_features.append(col)
print(input_features)

输出结果为我们可以使用的输出特征：

['open', 'high', 'low', 'close', 'vol', 'amount', 'open_lag_1', 'open_lag_2', 'open_lag_3', 'open_lag_4', 'high_lag_1', 'high_lag_2', 'high_lag_3', 'high_lag_4', 'low_lag_1', 'low_lag_2', 'low_lag_3', 'low_lag_4', 'close_lag_1', 'close_lag_2', 'close_lag_3', 'close_lag_4', 'vol_lag_1', 'vol_lag_2', 'vol_lag_3', 'vol_lag_4', 'amount_lag_1', 'amount_lag_2', 'amount_lag_3', 'amount_lag_4']

然后我们用TA-Lib手动创建一些特征，看看之后能否对我们的预测效果带来提升。

import talib
extended_input_features = input_features[:]  # 在input_features上添加
# 加上过去5天的指数滑动平均和布林带指标，大家也可以尝试加入更多的输入特征
timeperiod = 5
lagged_df['ema5'] = talib.EMA(lagged_df['close'], timeperiod=timeperiod)
extended_input_features.append('ema5')

nbdevup = 2  # 上下都是2倍标准差
nbdevdn = 2
lagged_df['upperband5'], lagged_df['middleband5'], lagged_df['lowerband5'] = talib.BBANDS(lagged_df['close'], timeperiod=timeperiod, nbdevup=nbdevup, nbdevdn=nbdevdn)
extended_input_features.append('upperband5')
extended_input_features.append('middleband5')
extended_input_features.append('lowerband5')
lagged_df = lagged_df.dropna()  # 删除最开始有缺失值的几行
print(extended_input_features)

输出为：

['open', 'high', 'low', 'close', 'vol', 'amount', 'open_lag_1', 'open_lag_2', 'open_lag_3', 'open_lag_4', 'high_lag_1', 'high_lag_2', 'high_lag_3', 'high_lag_4', 'low_lag_1', 'low_lag_2', 'low_lag_3', 'low_lag_4', 'close_lag_1', 'close_lag_2', 'close_lag_3', 'close_lag_4', 'vol_lag_1', 'vol_lag_2', 'vol_lag_3', 'vol_lag_4', 'amount_lag_1', 'amount_lag_2', 'amount_lag_3', 'amount_lag_4', 'ema5', 'upperband5', 'middleband5', 'lowerband5']

在新的DataFrame上重新划分训练集、验证集和测试集。

train = lagged_df['2010-01-01':'2017-12-31']  # 我们使用2010-2017年的数据作为训练集
cv = lagged_df['2018-01-01': '2018-12-31']  # 我们使用2018年的数据