30  关联规则

30.1 安装有关软件包

原书附带的代码有点问题，因此我们使用一个比较知名的Python的包mlxtend。

请执行anaconda prompt（在Windows开始菜单中有），或者Terminal（对苹果用户）。

运行:

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple mlxtend

上述代码将会从清华大学的镜像中安装mlxtend。

30.2 基本概念

Support 支持度

所谓support，指的是“事件发生的概率”，就是概率论中的那个概率。因此，显然发生1个事件，或者同时发生多个事件的概率，这里都是support。

例如：发生A事件

\[ Support(A) = \frac{A发生的次数}{所有事件的数量} \]

A和B同时发生

\[ Support(A \rightarrow B) = \frac{A，B同时发生的次数}{所有事件的数量} \]

Confidence 支持度

所谓confidence，指的是“事件发生的条件概率”，同样就是概率论中的概念。

例如，发生A的情况下，发生B的概率。

\[ Confidence(A \rightarrow B) = P(B|A) = \frac{Support(A \cap B)}{Support(A)} \]

这里的概念和概率论中的完全一致，不需要展开。

Lift 提升度

衡量A的发生，对B发生的概率的提升程度。大于1为有提升（如2种商品总是一起买），等于1为无变化，小于1为有下降（如两种商品往往选了A就不会选B）。

30.3 购物清单示例

读取数据，手动进行一下计算。

数据：超市购买记录，数据位于case6/tr.xlsx

# 数据没有题头行，所以使用header=None
df = pd.read_excel('data/case6/tr.xlsx',header=None).set_index(0)
df
1 2 3 4
0
I1 西红柿 排骨 鸡蛋 NaN
I2 西红柿 茄子 NaN NaN
I3 鸡蛋 袜子 NaN NaN
I4 西红柿 排骨 茄子 NaN
I5 西红柿 排骨 袜子 酸奶
I6 鸡蛋 茄子 酸奶 NaN
I7 排骨 鸡蛋 茄子 NaN
I8 土豆 鸡蛋 袜子 NaN
I9 西红柿 排骨 鞋子 土豆

数据中，每一行是一个购物单（即一个事务），那么总事务的数量是9。

我们计算一下：

  1. Support(排骨, 西红柿)：同时购买排骨和西红柿的概率。
  2. Confidence(排骨, 西红柿)：购买了排骨的情况下，会购买西红柿的概率。

目测：

  1. 同时购买了排骨和西红柿的单子有4个，[1,4,5,9]，那么Support(排骨, 西红柿) = 4/9 = 0.4444
  2. 购买了排骨的单子有5个，[1,4,5,7,9]，那么Confidence(排骨, 西红柿) = 4/5 = 0.8

计算的结果和课本p119的表格相同。

30.4 用代码完成

这里使用来自mlxtend包的Apriori，同样是计算元素之间的Support和Confidence，但其计算速度比你自己写代码穷举要更快。

这也是我们一直的要求：尽量不要自己写算法，而是使用现成的代码。

mlxtend 安装见本章开头

30.4.1 数据预处理：

首先要把数据转为 one-hot 编码的形式（教材 p117 表格6-1）:

每一行是一个购物单，列是商品名称，数据是01变量，表示这个购物单中是否有这项商品，有则为1，反之为0。

# 后面的转换函数不接受“字符串和np.nan同时存在”的数据，因此需要把缺失值先进行填充
# 我们把缺失值填充为任何字符串，比如'x'
tmp_df = df.fillna('x')
tmp_df
1 2 3 4
0
I1 西红柿 排骨 鸡蛋 x
I2 西红柿 茄子 x x
I3 鸡蛋 袜子 x x
I4 西红柿 排骨 茄子 x
I5 西红柿 排骨 袜子 酸奶
I6 鸡蛋 茄子 酸奶 x
I7 排骨 鸡蛋 茄子 x
I8 土豆 鸡蛋 袜子 x
I9 西红柿 排骨 鞋子 土豆

使用TransactionEncoder，可以把上面的DataFrame转为one-hot格式的数据。

但是有关的函数不接收DataFrame格式的数据，我们可以用tmp_df.values，取得数据对应的ndarray的版本。

tmp_df.values
array([['西红柿', '排骨', '鸡蛋', 'x'],
       ['西红柿', '茄子', 'x', 'x'],
       ['鸡蛋', '袜子', 'x', 'x'],
       ['西红柿', '排骨', '茄子', 'x'],
       ['西红柿', '排骨', '袜子', '酸奶'],
       ['鸡蛋', '茄子', '酸奶', 'x'],
       ['排骨', '鸡蛋', '茄子', 'x'],
       ['土豆', '鸡蛋', '袜子', 'x'],
       ['西红柿', '排骨', '鞋子', '土豆']], dtype=object)

对数上述数据，先fit，再transform即可

from mlxtend.preprocessing import TransactionEncoder
te = TransactionEncoder()


te_data = te.fit(tmp_df.values).transform(tmp_df.values) 
te_data
array([[ True, False,  True, False, False,  True, False, False,  True],
       [ True, False, False,  True, False,  True, False, False, False],
       [ True, False, False, False,  True, False, False, False,  True],
       [ True, False,  True,  True, False,  True, False, False, False],
       [False, False,  True, False,  True,  True,  True, False, False],
       [ True, False, False,  True, False, False,  True, False,  True],
       [ True, False,  True,  True, False, False, False, False,  True],
       [ True,  True, False, False,  True, False, False, False,  True],
       [False,  True,  True, False, False,  True, False,  True, False]])

te.columns_：可以获得对应的列名

te.columns_
['x', '土豆', '排骨', '茄子', '袜子', '西红柿', '酸奶', '鞋子', '鸡蛋']

把数据转为DataFrame，就能看得更清楚：

data = pd.DataFrame(te_data,columns=te.columns_).drop('x',axis=1) 
data * 1
土豆 排骨 茄子 袜子 西红柿 酸奶 鞋子 鸡蛋
0 0 1 0 0 1 0 0 1
1 0 0 1 0 1 0 0 0
2 0 0 0 1 0 0 0 1
3 0 1 1 0 1 0 0 0
4 0 1 0 1 1 1 0 0
5 0 0 1 0 0 1 0 1
6 0 1 1 0 0 0 0 1
7 1 0 0 1 0 0 0 1
8 1 1 0 0 1 0 1 0

这个表格就是教材p117的表6-1，只是我们采用的是现成的代码去处理。

对比一下原始数据的第一行，和one-hot数据的第一行，都表达了同一个含义：第一张购物单是[西红柿, 排骨, 鸡蛋]

df.head(1)
1 2 3 4
0
I1 西红柿 排骨 鸡蛋 NaN

30.4.2 计算关联

按教材的做法，找出最小支持度为0.2的关系。

首先计算所有商品以及商品集的support。即某个商品集出现在全部清单中的频率。

from mlxtend.frequent_patterns import apriori

# min_support： 支持度不低于0.2
# use_colnames：把列名（商品名）作为itemsets
freq_data = apriori(data, min_support=0.2,use_colnames=True).sort_values('support',ascending=False)
freq_data
support itemsets
1 0.555556 (排骨)
4 0.555556 (西红柿)
6 0.555556 (鸡蛋)
2 0.444444 (茄子)
8 0.444444 (西红柿, 排骨)
3 0.333333 (袜子)
0 0.222222 (土豆)
5 0.222222 (酸奶)
7 0.222222 (茄子, 排骨)
9 0.222222 (鸡蛋, 排骨)
10 0.222222 (西红柿, 茄子)
11 0.222222 (茄子, 鸡蛋)
12 0.222222 (袜子, 鸡蛋)

显而易见，itemsets就是商品集，support就是支持度（出现在购物单中的次数），如support(西红柿，排骨) = 0.4444，和我们前面算的一样。

做个简单的筛选，选择商品数量大于等于2的行

freq_data[freq_data.itemsets.apply(len) >=2]
support itemsets
8 0.444444 (西红柿, 排骨)
7 0.222222 (茄子, 排骨)
9 0.222222 (鸡蛋, 排骨)
10 0.222222 (西红柿, 茄子)
11 0.222222 (茄子, 鸡蛋)
12 0.222222 (袜子, 鸡蛋)

计算出频率之后，可以进一步计算关联规则。

association_rules函数：

  1. 第一个参数就是前面计算频率数据
  2. metric：你希望用什么指标来做筛选，比如这里用confidence来筛选
  3. min_threshold：对于metric中指定的那个指标，你筛选出值不低于多少

对于课本的例子，置信度不低于0.4，那么后2个参数就是metric='confidence',min_threshold=0.5

from mlxtend.frequent_patterns import association_rules

rules = association_rules(freq_data, metric='confidence',min_threshold=0.5)
rules
antecedents consequents antecedent support consequent support support confidence lift leverage conviction
0 (西红柿) (排骨) 0.555556 0.555556 0.444444 0.800000 1.44 0.135802 2.222222
1 (排骨) (西红柿) 0.555556 0.555556 0.444444 0.800000 1.44 0.135802 2.222222
2 (茄子) (排骨) 0.444444 0.555556 0.222222 0.500000 0.90 -0.024691 0.888889
3 (茄子) (西红柿) 0.444444 0.555556 0.222222 0.500000 0.90 -0.024691 0.888889
4 (茄子) (鸡蛋) 0.444444 0.555556 0.222222 0.500000 0.90 -0.024691 0.888889
5 (袜子) (鸡蛋) 0.333333 0.555556 0.222222 0.666667 1.20 0.037037 1.333333

选择对应的列，再排序一下，就是课本p124的结果。

rules.iloc[:,[0,1,4,5]].sort_values(['support','confidence'],ascending=False)
antecedents consequents support confidence
0 (西红柿) (排骨) 0.444444 0.800000
1 (排骨) (西红柿) 0.444444 0.800000
5 (袜子) (鸡蛋) 0.222222 0.666667
2 (茄子) (排骨) 0.222222 0.500000
3 (茄子) (西红柿) 0.222222 0.500000
4 (茄子) (鸡蛋) 0.222222 0.500000

可以很简单地解释一下数据，如第一行：对于每张购物单，西红柿和排骨同时出现的频率是0.4444，西红柿出现的情况下，排骨出现的概率是0.8。

30.4.3 测试

from mlxtend.frequent_patterns import fpgrowth

fpgrowth(data, min_support=0.2,use_colnames=True)
support itemsets
0 0.555556 (鸡蛋)
1 0.555556 (西红柿)
2 0.555556 (排骨)
3 0.444444 (茄子)
4 0.333333 (袜子)
5 0.222222 (酸奶)
6 0.222222 (土豆)
7 0.444444 (西红柿, 排骨)
8 0.222222 (鸡蛋, 排骨)
9 0.222222 (西红柿, 茄子)
10 0.222222 (茄子, 排骨)
11 0.222222 (茄子, 鸡蛋)
12 0.222222 (袜子, 鸡蛋)

30.5 国际股票指数关联示例

问：国际股票指数（包括我国的上证指数）之间，大跌有没有关联性？如美股跳水，A股会跳水吗？

用关联规则的话说，就是把每一国的指数大跌作为一个事件，那么这些事件的support和confidence（概率和条件概率）是什么？

df = pd.read_excel('data/case6/国际股票价格指数日交易数据表.xlsx',
                    converters={'Trddt':pd.to_datetime})
df.head()
Indexcd Trddt Opnidx Highidx Lowidx Clsidx
0 DJI 2010-01-04 10430.69 10604.97 10430.69 10583.96
1 DJI 2010-01-05 10584.56 10584.56 10522.52 10572.02
2 DJI 2010-01-06 10564.72 10594.99 10546.55 10573.68
3 DJI 2010-01-07 10571.11 10612.37 10505.21 10606.86
4 DJI 2010-01-08 10606.40 10619.40 10554.33 10618.19 
index_info_df = pd.read_excel('data/case6/IDX_Gidxinfo.xlsx')
index_info_df.head()
Indexcd Idxnme
0 DJI 美国道琼斯工业指数
1 FCHI 法国CAC40指数
2 FTSE 英国富时100指数
3 GDAXI 德国DAX指数
4 HSI 香港恒生指数 
data = pd.merge(index_info_df,df,on='Indexcd')
data.head()
Indexcd Idxnme Trddt Opnidx Highidx Lowidx Clsidx
0 DJI 美国道琼斯工业指数 2010-01-04 10430.69 10604.97 10430.69 10583.96
1 DJI 美国道琼斯工业指数 2010-01-05 10584.56 10584.56 10522.52 10572.02
2 DJI 美国道琼斯工业指数 2010-01-06 10564.72 10594.99 10546.55 10573.68
3 DJI 美国道琼斯工业指数 2010-01-07 10571.11 10612.37 10505.21 10606.86
4 DJI 美国道琼斯工业指数 2010-01-08 10606.40 10619.40 10554.33 10618.19

计算每只股票的跌幅。

注意我们不知道原始数据有没有排序，所以我们这里先排序一下，再进行计算。

data['ret'] = data.sort_values(['Indexcd','Trddt']).groupby('Indexcd')['Clsidx'].pct_change()
data['jump'] = (data['ret'] < -0.005) * 1 # 跌幅是否大于0.5%？
data.head()
Indexcd Idxnme Trddt Opnidx Highidx Lowidx Clsidx ret jump
0 DJI 美国道琼斯工业指数 2010-01-04 10430.69 10604.97 10430.69 10583.96 NaN 0
1 DJI 美国道琼斯工业指数 2010-01-05 10584.56 10584.56 10522.52 10572.02 -0.001128 0
2 DJI 美国道琼斯工业指数 2010-01-06 10564.72 10594.99 10546.55 10573.68 0.000157 0
3 DJI 美国道琼斯工业指数 2010-01-07 10571.11 10612.37 10505.21 10606.86 0.003138 0
4 DJI 美国道琼斯工业指数 2010-01-08 10606.40 10619.40 10554.33 10618.19 0.001068 0

取得我们需要的列：指数名称，日期，和是否跳水，再长数据转宽数据，然后去掉含有na的行（即只保留所以指数共同的交易日）

这就是我们要用来计算关联规则的材料

oh_data = data[['Idxnme', 'Trddt', 'jump']].pivot(
    index='Trddt', columns='Idxnme', values='jump').dropna().astype(int)

oh_data.iloc[:5,:5]
Idxnme 中国沪深300指数 俄罗斯Micex指数 印度孟买30指数 台湾加权指数 富时新加坡海峡时报指数
Trddt
2010-01-12 0 1 1 0 1
2010-01-13 1 0 0 1 1
2010-01-14 0 0 0 0 0
2010-01-15 0 0 0 0 0
2010-01-19 0 0 1 1 0

按课本的例子，选择支持度不低于0.08，支持度不低于0.9的指数集合。

首先计算支持度（频率数据），最小支持度是0.08。

freq_data = apriori(oh_data, min_support=0.08,use_colnames=True).sort_values('support',ascending=False)
freq_data
/home/lee/anaconda3/lib/python3.9/site-packages/mlxtend/frequent_patterns/fpcommon.py:111: DeprecationWarning: DataFrames with non-bool types result in worse computationalperformance and their support might be discontinued in the future.Please use a DataFrame with bool type
  warnings.warn(
support itemsets
1 0.303378 (俄罗斯Micex指数)
0 0.298917 (中国沪深300指数)
7 0.288719 (法国CAC40指数)
6 0.288082 (日本日经225)
11 0.275972 (香港恒生指数)
... ... ...
75 0.080306 (中国沪深300指数, 韩国KOSPI指数, 台湾加权指数)
97 0.080306 (印度孟买30指数, 台湾加权指数, 韩国KOSPI指数)
114 0.080306 (香港恒生指数, 德国DAX指数, 富时新加坡海峡时报指数)
123 0.080306 (香港恒生指数, 德国DAX指数, 日本日经225)
86 0.080306 (美国道琼斯工业指数, 德国DAX指数, 俄罗斯Micex指数)

140 rows × 2 columns

用频率数据计算置信度，最小置信度是0.9。

rules = association_rules(freq_data, metric='confidence', min_threshold=0.9)

rules.iloc[:, [0, 1, 4, 5]].sort_values('confidence', ascending=False)
antecedents consequents support confidence
6 (美国道琼斯工业指数, 德国DAX指数, 英国富时100指数) (法国CAC40指数) 0.097514 0.968354
10 (印度孟买30指数, 德国DAX指数, 英国富时100指数) (法国CAC40指数) 0.087317 0.958042
4 (美国道琼斯工业指数, 英国富时100指数) (法国CAC40指数) 0.104525 0.937143
1 (美国道琼斯工业指数, 德国DAX指数) (法国CAC40指数) 0.113448 0.936842
7 (美国道琼斯工业指数, 英国富时100指数, 法国CAC40指数) (德国DAX指数) 0.097514 0.932927
11 (德国DAX指数, 日本日经225, 英国富时100指数) (法国CAC40指数) 0.085405 0.930556
2 (德国DAX指数, 俄罗斯Micex指数, 英国富时100指数) (法国CAC40指数) 0.109624 0.929730
8 (香港恒生指数, 德国DAX指数, 英国富时100指数) (法国CAC40指数) 0.095602 0.914634
9 (香港恒生指数, 英国富时100指数, 法国CAC40指数) (德国DAX指数) 0.095602 0.914634
0 (德国DAX指数, 英国富时100指数) (法国CAC40指数) 0.168260 0.913495
12 (法国CAC40指数, 日本日经225, 英国富时100指数) (德国DAX指数) 0.085405 0.911565
3 (印度孟买30指数, 德国DAX指数) (法国CAC40指数) 0.105163 0.906593
5 (美国道琼斯工业指数, 英国富时100指数) (德国DAX指数) 0.100701 0.902857

这就是课本p129的表6-7：

  1. support：antecedents和consequents同时下跌超过0.5%，占所有交易日的比例。
  2. confidence：如果antecedents都下跌超过0.5%，那么consequents也会跳水超过0.5%的概率。

那么第1行的解释是：

  1. (德国DAX指数, 美国道琼斯工业指数, 英国富时100指数) (法国CAC40指数)，这4个指数同时跌超过0.5%，占所有共同交易日9.75%。所有共同的交易日中，一起跳水的日子大概10%。
  2. 如果(德国DAX指数, 美国道琼斯工业指数, 英国富时100指数) 同时下跌超过0.5%，那么(法国CAC40指数)下跌超过0.5%的概率96.84%。前3个指数如果跳水，法国指数也几乎必然跳水。

如果要计算涨幅的关联，显然也完全类似，只要把data['ret'] < -0.005改一下即可

30.5.1 一对一规则

从样本中提取1对1的规则，只要把antecedents和consequents里面的元素数量同时为1的行选出来即可。

这里借用教材p128的数值，最小支持度0.1，最小置信度0.6。

freq_data = apriori(oh_data, min_support=0.1,use_colnames=True).sort_values('support',ascending=False)
rules = association_rules(freq_data, metric='confidence', min_threshold=0.6)
rules.iloc[:, [0, 1, 4, 5]].head()
/home/lee/anaconda3/lib/python3.9/site-packages/mlxtend/frequent_patterns/fpcommon.py:111: DeprecationWarning: DataFrames with non-bool types result in worse computationalperformance and their support might be discontinued in the future.Please use a DataFrame with bool type
  warnings.warn(
antecedents consequents support confidence
0 (德国DAX指数) (法国CAC40指数) 0.224984 0.836493
1 (法国CAC40指数) (德国DAX指数) 0.224984 0.779249
2 (法国CAC40指数) (英国富时100指数) 0.191842 0.664459
3 (英国富时100指数) (法国CAC40指数) 0.191842 0.775773
4 (德国DAX指数) (英国富时100指数) 0.184194 0.684834 
# antecedents和consequents中的元素数量同时为1的行
one_rule = np.logical_and(rules.antecedents.apply(len) == 1,rules.consequents.apply(len) == 1)
rules[one_rule].iloc[:,[0, 1,4,5]].reset_index(drop=True)
antecedents consequents support confidence
0 (德国DAX指数) (法国CAC40指数) 0.224984 0.836493
1 (法国CAC40指数) (德国DAX指数) 0.224984 0.779249
2 (法国CAC40指数) (英国富时100指数) 0.191842 0.664459
3 (英国富时100指数) (法国CAC40指数) 0.191842 0.775773
4 (德国DAX指数) (英国富时100指数) 0.184194 0.684834
5 (英国富时100指数) (德国DAX指数) 0.184194 0.744845
6 (富时新加坡海峡时报指数) (香港恒生指数) 0.144041 0.704050
7 (韩国KOSPI指数) (香港恒生指数) 0.139579 0.636628
8 (韩国KOSPI指数) (日本日经225) 0.138942 0.633721
9 (美国道琼斯工业指数) (法国CAC40指数) 0.130019 0.680000
10 (美国道琼斯工业指数) (德国DAX指数) 0.121096 0.633333

这个表即课本p128-p129的表6-6。

30.6 行业联动和轮动

对应课本第11章，p223

info_df = pd.read_excel('data/case11/指数基本信息表.xlsx').iloc[:,[0,1]]
info_df.columns = ['指数代码','指数名称']
info_df.head()
指数代码 指数名称
0 801010 农林牧渔
1 801020 采掘
2 801030 化工
3 801040 黑色金属
4 801050 有色金属 
price_df = pd.read_excel('data/case11/指数交易数据表20100104-20170307.xlsx',converters={'交易日期':pd.to_datetime})
price_df
指数代码 交易日期 收盘价
0 801010 2010-01-04 2058.27
1 801010 2010-01-05 2081.28
2 801010 2010-01-06 2048.70
3 801010 2010-01-07 1999.31
4 801010 2010-01-08 2018.81
... ... ... ...
54731 801890 2017-03-01 1673.23
54732 801890 2017-03-02 1666.43
54733 801890 2017-03-03 1667.51
54734 801890 2017-03-06 1692.56
54735 801890 2017-03-07 1694.80

54736 rows × 3 columns

和前面一样，我们把指数名称merge进价格序列。

price_df = pd.merge(price_df,info_df)
price_df
指数代码 交易日期 收盘价 指数名称
0 801010 2010-01-04 2058.27 农林牧渔
1 801010 2010-01-05 2081.28 农林牧渔
2 801010 2010-01-06 2048.70 农林牧渔
3 801010 2010-01-07 1999.31 农林牧渔
4 801010 2010-01-08 2018.81 农林牧渔
... ... ... ... ...
54731 801890 2017-03-01 1673.23 机械设备
54732 801890 2017-03-02 1666.43 机械设备
54733 801890 2017-03-03 1667.51 机械设备
54734 801890 2017-03-06 1692.56 机械设备
54735 801890 2017-03-07 1694.80 机械设备

54736 rows × 4 columns

查看每个行业指数的长度

length_df = price_df.groupby('指数代码').agg(len).reset_index()
length_df
指数代码 交易日期 收盘价 指数名称
0 801010 1741 1741 1741
1 801020 1741 1741 1741
2 801030 1741 1741 1741
3 801040 1741 1741 1741
4 801050 1741 1741 1741
5 801060 998 998 998
6 801070 998 998 998
7 801080 1741 1741 1741
8 801090 998 998 998
9 801100 998 998 998
10 801110 1741 1741 1741
11 801120 1741 1741 1741
12 801130 1741 1741 1741
13 801140 1741 1741 1741
14 801150 1741 1741 1741
15 801160 1741 1741 1741
16 801170 1741 1741 1741
17 801180 1741 1741 1741
18 801190 998 998 998
19 801200 1741 1741 1741
20 801210 1741 1741 1741
21 801220 998 998 998
22 801230 1741 1741 1741
23 801710 1741 1741 1741
24 801720 1741 1741 1741
25 801730 1741 1741 1741
26 801740 1741 1741 1741
27 801750 1741 1741 1741
28 801760 1741 1741 1741
29 801770 1741 1741 1741
30 801780 1741 1741 1741
31 801790 1741 1741 1741
32 801880 1741 1741 1741
33 801890 1741 1741 1741

某些行业指数短一点，可能是比较晚才推出。提取交易日长度为1741的行业。

full_len_codes = length_df.query("交易日期 == 1741")['指数代码']
full_len_codes
0     801010
1     801020
2     801030
3     801040
4     801050
7     801080
10    801110
11    801120
12    801130
13    801140
14    801150
15    801160
16    801170
17    801180
19    801200
20    801210
22    801230
23    801710
24    801720
25    801730
26    801740
27    801750
28    801760
29    801770
30    801780
31    801790
32    801880
33    801890
Name: 指数代码, dtype: int64

在交易信息price_df中，提取指数代码属于full_len_codes，用.isin()方法。

mask = price_df['指数代码'].isin(full_len_codes)
mask
0        True
1        True
2        True
3        True
4        True
         ... 
54731    True
54732    True
54733    True
54734    True
54735    True
Name: 指数代码, Length: 54736, dtype: bool
price_df = price_df[mask]
price_df
指数代码 交易日期 收盘价 指数名称
0 801010 2010-01-04 2058.27 农林牧渔
1 801010 2010-01-05 2081.28 农林牧渔
2 801010 2010-01-06 2048.70 农林牧渔
3 801010 2010-01-07 1999.31 农林牧渔
4 801010 2010-01-08 2018.81 农林牧渔
... ... ... ... ...
54731 801890 2017-03-01 1673.23 机械设备
54732 801890 2017-03-02 1666.43 机械设备
54733 801890 2017-03-03 1667.51 机械设备
54734 801890 2017-03-06 1692.56 机械设备
54735 801890 2017-03-07 1694.80 机械设备

48748 rows × 4 columns

简单检查一下，按代码分组后，是否全部分组的长度都是1741

all(price_df.groupby('指数代码').agg(len)['收盘价'] == 1741)
True

那么现在price_df中的全部指数，交易日的长度都一致了。

30.6.1 月频涨跌

首先计算月频下的涨跌幅和是否上涨。

注意因为我们不确定价格数据是不是完全按交易日期排序，保险起见我们手动排序一次。

  1. 在计算收益率之前，我们要对“每个指数代码内部按日期排序”，但因为不能分组后再组内排序，因为我们用多重排序来达到同样的效果。
  2. 排序后，按指数代码分组，取出收盘价计算日收益率，以及是否上涨的标志。
  3. 去掉含有na的行。
price_w_df = price_df.sort_values(['指数代码','交易日期']).set_index('交易日期').groupby('指数代码').resample('M').last()
price_w_df
指数代码 收盘价 指数名称
指数代码 交易日期
801010 2010-01-31 801010 2021.48 农林牧渔
2010-02-28 801010 2132.22 农林牧渔
2010-03-31 801010 2105.68 农林牧渔
2010-04-30 801010 2034.38 农林牧渔
2010-05-31 801010 1923.09 农林牧渔
... ... ... ... ...
801890 2016-11-30 801890 1731.66 机械设备
2016-12-31 801890 1623.81 机械设备
2017-01-31 801890 1594.32 机械设备
2017-02-28 801890 1667.32 机械设备
2017-03-31 801890 1694.80 机械设备

2436 rows × 3 columns

price_w_df = price_w_df.drop('指数代码',axis = 1).reset_index()
price_w_df
指数代码 交易日期 收盘价 指数名称
0 801010 2010-01-31 2021.48 农林牧渔
1 801010 2010-02-28 2132.22 农林牧渔
2 801010 2010-03-31 2105.68 农林牧渔
3 801010 2010-04-30 2034.38 农林牧渔
4 801010 2010-05-31 1923.09 农林牧渔
... ... ... ... ...
2431 801890 2016-11-30 1731.66 机械设备
2432 801890 2016-12-31 1623.81 机械设备
2433 801890 2017-01-31 1594.32 机械设备
2434 801890 2017-02-28 1667.32 机械设备
2435 801890 2017-03-31 1694.80 机械设备

2436 rows × 4 columns

price_w_df['月收益率'] = price_w_df['收盘价'].pct_change()
price_w_df['月up']  =  (price_w_df['月收益率'] > 0) * 1
price_w_df
指数代码 交易日期 收盘价 指数名称 月收益率 月up
0 801010 2010-01-31 2021.48 农林牧渔 NaN 0
1 801010 2010-02-28 2132.22 农林牧渔 0.054782 1
2 801010 2010-03-31 2105.68 农林牧渔 -0.012447 0
3 801010 2010-04-30 2034.38 农林牧渔 -0.033861 0
4 801010 2010-05-31 1923.09 农林牧渔 -0.054705 0
... ... ... ... ... ... ...
2431 801890 2016-11-30 1731.66 机械设备 0.059553 1
2432 801890 2016-12-31 1623.81 机械设备 -0.062281 0
2433 801890 2017-01-31 1594.32 机械设备 -0.018161 0
2434 801890 2017-02-28 1667.32 机械设备 0.045788 1
2435 801890 2017-03-31 1694.80 机械设备 0.016482 1

2436 rows × 6 columns

data = price_w_df.pivot(index = '交易日期',columns='指数名称',values = '月up').dropna().astype(int)
data.iloc[:,:7]
指数名称 交通运输 传媒 公用事业 农林牧渔 化工 医药生物 商业贸易
交易日期
2010-01-31 0 0 0 0 0 1 0
2010-02-28 1 1 1 1 1 1 1
2010-03-31 1 1 1 0 1 0 1
2010-04-30 0 0 0 0 0 1 0
2010-05-31 0 0 0 0 0 1 0
... ... ... ... ... ... ... ...
2016-11-30 1 1 1 1 1 1 1
2016-12-31 0 0 0 0 0 0 0
2017-01-31 1 0 0 0 0 0 0
2017-02-28 1 1 1 1 1 1 1
2017-03-31 0 1 1 1 1 1 0

87 rows × 7 columns

在数据中，2017年的3月只有几个交易日，因此可以去掉（教材p235）。

data = data.loc[:'2017-02']
data.iloc[:,:7]
指数名称 交通运输 传媒 公用事业 农林牧渔 化工 医药生物 商业贸易
交易日期
2010-01-31 0 0 0 0 0 1 0
2010-02-28 1 1 1 1 1 1 1
2010-03-31 1 1 1 0 1 0 1
2010-04-30 0 0 0 0 0 1 0
2010-05-31 0 0 0 0 0 1 0
... ... ... ... ... ... ... ...
2016-10-31 1 1 1 1 1 1 1
2016-11-30 1 1 1 1 1 1 1
2016-12-31 0 0 0 0 0 0 0
2017-01-31 1 0 0 0 0 0 0
2017-02-28 1 1 1 1 1 1 1

86 rows × 7 columns

上述表格即教材p238的图11-8。注意教材用的算法是最后一天的价格减去第一天价格，这个算法不正确，但在月度上和pct_change()只有一点差异。

30.6.2 月频关联度

接下来就是计算关联度。按课本的参数，最小支持度为0.47，最小置信度为0.9。

注意：

  1. 最小支持度（事件或者多个事件发生的概率）如果设置得太低，那么就会产生极大量的组合，运算时间可能极长，并且你的内存可能不足以完成运算。因此这个值一般先设置的得大一点。
  2. fpgrowth算法和apriori算法的目的相同，但速度更快，所以下面的代码也可以直接把apriori()函数改为fpgrowth()函数，参数一样。
from mlxtend.frequent_patterns import apriori, fpgrowth, association_rules

freq_data = apriori(data, min_support=0.47,use_colnames=True).sort_values('support',ascending=False).reset_index(drop=True)
freq_data
/home/lee/anaconda3/lib/python3.9/site-packages/mlxtend/frequent_patterns/fpcommon.py:111: DeprecationWarning: DataFrames with non-bool types result in worse computationalperformance and their support might be discontinued in the future.Please use a DataFrame with bool type
  warnings.warn(
support itemsets
0 0.616279 (医药生物)
1 0.604651 (家用电器)
2 0.593023 (有色金属)
3 0.581395 (建筑材料)
4 0.581395 (机械设备)
... ... ...
158 0.476744 (电气设备, 轻工制造)
159 0.476744 (通信, 电子元器件)
160 0.476744 (电子元器件, 综合)
161 0.476744 (纺织服装, 电子元器件)
162 0.476744 (电子元器件, 电气设备, 计算机, 轻工制造, 机械设备)

163 rows × 2 columns

在月度的频率上，单个板块上涨，或者多个板块同时上涨，其概率高于0.47的，在上表中。

获得全部的规则，最低置信度为0.9

rules = association_rules(freq_data, metric='confidence', min_threshold=0.9)
rules.iloc[:, [0, 1, 4, 5]]
antecedents consequents support confidence
0 (计算机) (电子元器件) 0.546512 0.940000
1 (电子元器件) (计算机) 0.546512 0.979167
2 (化工) (轻工制造) 0.534884 0.958333
3 (轻工制造) (化工) 0.534884 0.978723
4 (采掘) (有色金属) 0.523256 0.937500
... ... ... ... ...
396 (电气设备, 机械设备) (计算机, 电子元器件, 轻工制造) 0.476744 0.953488
397 (计算机, 轻工制造) (机械设备, 电气设备, 电子元器件) 0.476744 0.911111
398 (计算机, 机械设备) (电气设备, 电子元器件, 轻工制造) 0.476744 0.931818
399 (机械设备, 轻工制造) (计算机, 电气设备, 电子元器件) 0.476744 0.953488
400 (电气设备) (机械设备, 计算机, 电子元器件, 轻工制造) 0.476744 0.931818

401 rows × 4 columns

找到1对1的规则，并且查看置信度最高的前几个。

mask = rules.apply(lambda x: len(x.antecedents) == 1 and len(x.consequents) == 1,axis=1)
rules.iloc[:, [0, 1, 4, 5]][mask].sort_values('confidence',ascending=False).head(10)
antecedents consequents support confidence
1 (电子元器件) (计算机) 0.546512 0.979167
3 (轻工制造) (化工) 0.534884 0.978723
36 (电气设备) (机械设备) 0.500000 0.977273
31 (电气设备) (计算机) 0.500000 0.977273
2 (化工) (轻工制造) 0.534884 0.958333
7 (轻工制造) (计算机) 0.523256 0.957447
42 (综合) (化工) 0.500000 0.955556
25 (综合) (轻工制造) 0.500000 0.955556
64 (电气设备) (电子元器件) 0.488372 0.954545
0 (计算机) (电子元器件) 0.546512 0.940000

30.6.3 测试

本月涨跌和下月的涨跌，有没有关联？

data_s1 = data.shift(-1)
data_s1.columns = [x + 's1' for x in data_s1.columns]
data_s1.iloc[:7,:7]
交通运输s1 传媒s1 公用事业s1 农林牧渔s1 化工s1 医药生物s1 商业贸易s1
交易日期
2010-01-31 1.0 1.0 1.0 1.0 1.0 1.0 1.0
2010-02-28 1.0 1.0 1.0 0.0 1.0 0.0 1.0
2010-03-31 0.0 0.0 0.0 0.0 0.0 1.0 0.0
2010-04-30 0.0 0.0 0.0 0.0 0.0 1.0 0.0
2010-05-31 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2010-06-30 1.0 1.0 1.0 1.0 1.0 1.0 1.0
2010-07-31 1.0 1.0 1.0 1.0 1.0 1.0 1.0 
data2 = pd.concat([data,data_s1],axis=1).dropna().astype(int)
data2
交通运输 传媒 公用事业 农林牧渔 化工 医药生物 商业贸易 国防军工 家用电器 建筑材料 ... 综合s1 计算机s1 轻工制造s1 通信s1 采掘s1 银行s1 非银金融s1 食品饮料s1 餐饮旅游s1 黑色金属s1
交易日期
2010-01-31 0 0 0 0 0 1 0 0 0 1 ... 1 1 1 1 1 1 0 1 1 1
2010-02-28 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 0 1 1 1 0 1 0
2010-03-31 1 1 1 0 1 0 1 0 1 1 ... 0 0 0 0 0 0 0 0 0 0
2010-04-30 0 0 0 0 0 1 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
2010-05-31 0 0 0 0 0 1 0 0 0 0 ... 0 0 0 0 0 0 1 0 0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2016-09-30 0 0 0 0 0 0 0 0 1 1 ... 1 1 1 1 1 1 1 1 1 1
2016-10-31 1 1 1 1 1 1 1 1 1 1 ... 0 1 1 1 1 1 1 1 0 1
2016-11-30 1 1 1 1 1 1 1 1 1 1 ... 0 0 0 0 0 0 0 0 0 0
2016-12-31 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 1 1 1 1 0 1
2017-01-31 1 0 0 0 0 0 0 1 1 1 ... 1 1 1 1 1 1 1 1 1 1

85 rows × 56 columns

from mlxtend.frequent_patterns import apriori, fpgrowth, association_rules

freq_data = fpgrowth(data2, min_support=0.4,use_colnames=True).sort_values('support',ascending=False).reset_index(drop=True)
freq_data
/home/lee/anaconda3/lib/python3.9/site-packages/mlxtend/frequent_patterns/fpcommon.py:111: DeprecationWarning: DataFrames with non-bool types result in worse computationalperformance and their support might be discontinued in the future.Please use a DataFrame with bool type
  warnings.warn(
support itemsets
0 0.611765 (家用电器s1)
1 0.611765 (医药生物)
2 0.611765 (医药生物s1)
3 0.600000 (家用电器)
4 0.588235 (机械设备s1)
... ... ...
24699 0.400000 (餐饮旅游s1, 化工s1, 通信s1, 电气设备s1, 医药生物s1, 农林牧渔s1)
24700 0.400000 (餐饮旅游s1, 化工s1, 通信s1, 电气设备s1, 医药生物s1, 机械设备s1)
24701 0.400000 (餐饮旅游s1, 化工s1, 通信s1, 电气设备s1, 医药生物s1, 机械设备s1, 农...
24702 0.400000 (餐饮旅游s1, 通信s1, 轻工制造s1, 电气设备s1, 医药生物s1, 机械设备s1)
24703 0.400000 (非银金融, 机械设备, 化工)

24704 rows × 2 columns

freq_data2 = freq_data[freq_data.itemsets.agg(len) <= 2]
rules = association_rules(freq_data2, metric='confidence', min_threshold=0.9)
rules.iloc[:, [0, 1, 4, 5]]
antecedents consequents support confidence
0 (电子元器件s1) (计算机s1) 0.552941 0.979167
1 (计算机s1) (电子元器件s1) 0.552941 0.959184
2 (计算机) (电子元器件) 0.541176 0.938776
3 (电子元器件) (计算机) 0.541176 0.978723
4 (轻工制造s1) (化工s1) 0.541176 0.978723
... ... ... ... ...
121 (传媒) (通信) 0.435294 0.902439
122 (传媒) (电气设备) 0.435294 0.902439
123 (非银金融) (家用电器) 0.435294 0.902439
124 (建筑装饰) (家用电器) 0.435294 0.902439
125 (建筑装饰) (建筑材料) 0.435294 0.902439

126 rows × 4 columns

mask1 = rules.antecedents.apply(lambda x: not any(['s1' in y for y in x ]))
mask2 = rules.consequents.apply(lambda x: all(['s1' in y for y in x ]))

mask = np.logical_and(mask1,mask2)

rules.iloc[:, [0, 1, 4, 5]][mask]
antecedents consequents support confidence

答案是似乎没有。