本文作者：何百圣 哈尔滨工业大学（威海） 经济管理学院 数量金融方向
MLAT系列文章为校内课程作业，以blog的形式记录作业。笔者的课程任务是12gradient_boosting，本篇的Create datasets和Feature_engineering属于数据处理阶段。
文中对于流程做了梳理，同时对部分关键点做了介绍，使之更加易于理解，为后来者做参考。
Imports & Settings
import pandas as pd
import numpy as np
from pathlib import Path
import seaborn as sns
import pandas_datareader.data as web
from statsmodels.regression.rolling import RollingOLS
import statsmodels.api as sm
idx = pd.IndexSlice
DATA_STORE = Path('E:\machine learning for algorithmic trading\wiki.h5')1. H5文件是层次数据格式第5代的版本（Hierarchical Data Format，HDF5），用来储存大型文件，是美国超级计算与应用中心研发的文件格式。这一点csdn上其他博主介绍的比较详细。
我发现里面可以存入许多dataframe，最后数据处理完成后文件内是这样的，可以看到每一个dataframe都很大。
2. idx = pd.IndexSlice 是多重索引的用法， 这里的两重索引一重为date一重为ticker。具体再另一篇文档里有介绍。
Quandl Wiki Prices
这里使用的股票数据来源于Quandl community，不过数据已经在2018年停用了，所以这里实际上使用的是历史数据。
df = (pd.read_csv('E:\machine learning for algorithmic trading\wiki.csv',
#parse意为解析，parse_dates参数是将date列解析为日期格式
parse_dates=['date'],
#index_col参数将date列和ticker列设为index
index_col=['date', 'ticker'],
#该参数一般与parse_date参数连用，用于规范格式
infer_datetime_format=True)
.sort_index())
print(df.info(null_counts=True))
#存储数据
with pd.HDFStore(DATA_STORE) as store:
store.put('quandl/wiki/prices', df)
3. 目前df的行列信息如图，可以看到有千万行，memory_usage为1.4GB，这样的大小是没法用csv或xlsx来保存的，这也是使用h5文件的原因。adj_xxx意为经过调整的价格，考虑了股息红利等因素。我们使用的数据也是adj_close。
df.head()
Metadata on US-traded companies
接下来获得美国上市公司的metadata，但是由于NASDAQ已经不允许利用url自动爬取数据，所以只能去官网手动下载三个交易所AMEX, NASDAQ, NYSE的数据，然后concat一下。
nasdaq = pd.read_csv(r'E:/machine learning for algorithmic trading/nasdaq.csv')
amex = pd.read_csv(r'E:/machine learning for algorithmic trading/amex.csv')
nyse = pd.read_csv(r'E:/machine learning for algorithmic trading/nyse.csv')
df2 = pd.concat([nasdaq,amex,nyse]).dropna(how = 'all',axis = 1)
df2 = df2.rename(columns=str.lower).set_index('symbol')
#去掉重复的index行
df2 = df2[~df2.index.duplicated()]
df2 = df2.drop(['% change','country','volume','net change'],axis = 1)
print(df2.info())
可以看到关于上市公司的信息，例如上市时间，市值等。
#存储数据
with pd.HDFStore(DATA_STORE) as store:
store.put('us_equities/stocks', df2)
Get Data
START = 2000
END = 2018
with pd.HDFStore(DATA_STORE) as store:
#prices只选取了2000-2018的adi_close列
prices = (store['quandl/wiki/prices']
.loc[idx[str(START):str(END), :], 'adj_close']
.unstack('ticker'))
#stocks选取市值，上市年份，行业
stocks = store['us_equities/stocks'].loc[:, ['marketcap', 'ipoyear', 'sector']]
#这里stocks已经去除了重复的index4. unstack是展开成dataframe，stack是堆叠成series。这里unstack的参数是ticker，也就是根据ticker展开。
prices.head()
stocks.head()
#保证选取的公司既有财务数据也有价格数据，但仍有nan值存在
shared = prices.columns.intersection(stocks.index)
stocks = stocks.loc[shared,:]
prices = prices.loc[:,shared]Create monthly return series
为了减少训练时间，我们将日频数据改为月度数据。
#resample重采样
monthly_prices = prices.resample('M').last()
monthly_prices.info()
5. resample意为重采样，参数控制采样频率，last意为选最后一天，此处变化很多，可选sum等，很有意思。
monthly_prices.head()
outlier_cutoff = 0.01
data = pd.DataFrame()
lags = [1, 2, 3, 6, 9, 12]
for lag in lags:
data[f'return_{lag}m'] = (monthly_prices
.pct_change(lag)
.stack()
.pipe(lambda x: x.clip(lower=x.quantile(outlier_cutoff),
upper=x.quantile(1-outlier_cutoff)))
#normalize
.add(1)
.pow(1/lag)
.sub(1)
)
data = data.swaplevel().dropna()
data.info()6. outlier_cutoff是 winsorize level，用clip()函数实现
data.tail()
#经过swaplevel交换index后
data = data.swaplevel().dropna()
data.head()
这时候我们可以单独抽出一个公司来看其数据格式
data.loc[idx[ 'AA',:],:].head()
data.info()
Drop stocks with less than 10 yrs of returns
通过grouby\size可以看出有些公司不足十年，这些公司不是可分析数据
注意：此时的数据为月度
#由于部分公司数据不足十年，剔除掉
data.groupby(level = 'ticker').size()
min_obs = 120
nobs = data.groupby(level='ticker').size()
keep = nobs[nobs>min_obs].index
data = data.loc[idx[keep,:], :]
data.info()
data.describe()
通过热力图可以看出收益率之间的关系
sns.clustermap(data.corr('spearman'),annot = True,center = 0,cmap = 'Blues' )
data.index.get_level_values('ticker').nunique()
#此时还剩下1580个公司
15807. get_level_values应用于多重索引结构，在多重索引中不能直接取单一index
Rolling Factor Betas
应用Fama-French五因子模型，分别代表：market risk, size, value, operating profitability, and investment。
利用 pandas-datareader来获取数据。
factors = ['Mkt-RF', 'SMB', 'HML', 'RMW', 'CMA']
factor_data = web.DataReader('F-F_Research_Data_5_Factors_2x3', 'famafrench', start='2000')[0].drop('RF', axis=1)
factor_data.index = factor_data.index.to_timestamp()
factor_data = factor_data.resample('M').last().div(100)
factor_data.index.name = 'date'
factor_data.head()
factor_data.info()
与data.return_1m合并
factor_data = factor_data.join(data['return_1m']).sort_index()
factor_data.info()
factor_data.head()
接下来利用RollingOLS对factor_data进行拟合，自变量为五因子，因变量为return_1m
T = 24
betas = (factor_data.groupby(level='ticker',
group_keys=False)
.apply(lambda x: RollingOLS(endog=x.return_1m,
exog=sm.add_constant(x.drop('return_1m', axis=1)),
window=min(T, x.shape[0]-1))
#Use params_only to skip all calculations except parameter estimation
.fit(params_only=True)
.params
.drop('const', axis=1)))
#实际上betas每一列已经不是原来的值了，已经变成了估计出的系数对betas的参数作describe()处理
betas.describe().join(betas.sum(1).describe().to_frame('total'))
cmap = sns.diverging_palette(10, 220, as_cmap=True)
sns.clustermap(betas.corr(), annot=True, cmap=cmap, center=0)8. 这是一个很有意思的写法，变量cmap（color of map）获得了一个调色板，应用于clustermap中，as_cmap为True意为返回一个调色板，为False则是一个数值列表，每个数字代表不同的颜色。annot、center等为常规写法，不做介绍。
将betas中取得的参数结果与data合并
data = (data
.join(betas
.groupby(level='ticker')
.shift()))
data.info()9. 这里的shift()写法有些奇怪，按理说不应该shift，因为原过程不存在任何错位现象。经查阅发现，groupby形成的结果是DataFrameGroupby类型，shift()后转化为DataFrame类型，并且由于RollingOLS是用大量数据估计而成，所以shift一下影响不大，但是格式转化上的作用比较明显。
（此处存疑，并不十分确定）
由于shift合并，故而存在nan值，选择用平均值插入替换nan值。
data.loc[:, factors] = data.groupby('ticker')[factors].apply(lambda x: x.fillna(x.mean()))
data.head()
Momentum factors
We can use these results to compute momentum factors based on the difference between returns over longer periods and the most recent monthly return, as well as for the difference between 3 and 12 month returns as follows:
for lag in [2,3,6,9,12]:
data[f'momentum_{lag}'] = data[f'return_{lag}m'].sub(data.return_1m)
data[f'momentum_3_12'] = data[f'return_12m'].sub(data.return_3m)其实就是一个简单用sub()实现的相减过程。
Date Indicators
dates = data.index.get_level_values('date')
data['year'] = dates.year
data['month'] = dates.month
Lagged returns
滞后回报率
for t in range(1, 7):
data[f'return_1m_t-{t}'] = data.groupby(level='ticker').return_1m.shift(t)Target: Holding Period Returns
持有期回报率
for t in [1,2,3,6,12]:
data[f'target_{t}m'] = data.groupby(level='ticker')[f'return_{t}m'].shift(-t)
10. 各种回报率的理解
cols = ['target_1m',
'target_2m',
'target_3m',
'return_1m',
'return_2m',
'return_3m',
'return_1m_t-1',
'return_1m_t-2',
'return_1m_t-3']
data[cols].dropna().sort_index().head(10)
#dropna 删除了收益滞后项导致的nan值
不难看出，对于A公司来说，2001-05的收益为-0.140220，这一回报率也成为了前一期的持有一月回报率，所以实际上持有期回报率（HPR）是一种预测的概念，用未来数据算出。
这一数值也成为了return_1m_t-1在2001-06的数据，所以lagged return实际上就是历史收益率的概念，用历史数据得出。
Create age proxy
用IPO年份来作为age proxy，因为年份由小到大，所以age越小，代表公司越老
data = (data
.join(pd.qcut(stocks.ipoyear, q=5, labels=list(range(1, 6)))
.astype(float)
.fillna(0)
.astype(int)
.to_frame('age')))
#按年份从大到，所以ipo最早的公司为1
data.age = data.age.fillna(-1)
data.head()
此时的data已经有许多column了，故而显示不全
Create dynamic size proxy
使用 marketcap来create size proxy.
size_factor = (monthly_prices
.loc[data.index.get_level_values('date').unique(),
data.index.get_level_values('ticker').unique()]
.sort_index(ascending=False)
.pct_change()
.fillna(0)
.add(1)
.cumprod())
size_factor.head()11. 这里有一个细节，sort_index时使用了ascending = False，所以相当于使第一行为1，也就是今天的市值，根据pct_change向前推演出市值变化，因为此时每一个公司的上市时间不一定， 所以这样更具有可比性。数据结果如下：
msize = (size_factor
.mul(stocks
.loc[size_factor.columns, 'marketcap'])).dropna(axis=1, how='all')
data['msize'] = (msize
.apply(lambda x: pd.qcut(x, q=10, labels=list(range(1, 11)))
.astype(int), axis=1)
.stack()
.swaplevel())
data.msize = data.msize.fillna(-1)Combine data
到了最后一步了，合并数据，此时data的数据已经很丰富了
data = data.join(stocks[['sector']])
data.sector = data.sector.fillna('Unknown')
data.info()
存储数据，完事了
with pd.HDFStore(DATA_STORE) as store:
store.put('engineered_features', data.sort_index().loc[idx[:, :'2018-03-01'], :])
print(store.info())
1150
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
