聚类、因子分析-python

<随手练习>

人工智能股票实证分析

一：研究背景和意义：

近几年人工智能的兴起，导致大量人工智能企业的诞生或崛起，故本文从4个维度选取15个主要财务指标对159家人工智能股票进行聚类分析和因子分析，主要将公司分成三大类和采用因子分析对目前公司进行综合排名，并对排名靠前的三只股票进行ARMA模型预测其未来一段时间的股票收盘价。

?

二：指标选取

???评估公司基本状况最为核心的财务指标就是公司的偿债能力、营运能力、盈利能力和发展能力。由此本文选取如下15个指标作为最能反映这些能力的指标进行分析：

表1 指标系数表

指标类型	指标名称	指标描述	编号
偿债能力指标	流动比率	比率越高，企业资产的变现能力越强，短期偿债能力亦越强	X1
	速度比率	衡量企业流动资产中可以立即变现用于偿还流动负债的能力	X2
	现金比率	现金比率越高，说明变现能力越强	X3
	资产负债率	反映债权人所提供的资本占全部资本的比例	X4
营运能力指标	应收账款周转率	衡量企业应收账款周转速度及管理效率的指标	X5
	存货周转率	衡量企业销售能力及存货管理水平的综合性指标	X6
	流动资产周转率	评价企业资产利用率的一个重要指标	X7
	总资产周转率	衡量资产投资规模与销售水平之间配比情况的指标	X8
盈利能力指标	总资产净利润率	反映的是公司运用全部资产所获得利润的水平	X9
	净资产收益率	衡量公司运用自有资本的效率	X10
	营业利润率	反映企业营业收入创造净利润的能力	X11
发展能力指标	主营业务收入增长率	评价企业成长状况和发展能力的重要指标	X12
	总资产增长率	反映企业本期资产规模的增长情况	X13
	基本每股收益	反映企业的经营成果，衡量普通股的获利水平及投资风险	X14
	净利润增长率	衡量一个企业经营效益的重要指标	X15

?

三：实证分析(python)

数据提取/清洗：

本文涉及的人工智能股票的财务数据均爬取自公布在网易财经上披露的近期年度财务指标数据。其主要数据显示如下：

from urllib2 import Request
from bs4 import BeautifulSoup

# 忽略未经核实的ssl证书认证
ssl._create_default_https_context=ssl._create_unverified_context()

req = Request.request.Request(xiaoqu_links,headers=hds[random.randint(0,len(hds)-1)])
source_code = Request.request.urlopen(req,timeout=10).read()
plain_text = source_code.decode(‘UTF-8‘) ?#,errors=‘ignore‘)
soup = BeautifulSoup(plain_text,"html.parser")

融资数据：

销售数据：

?

本文将对已有人工智能股票数据进行清洗和异常值删除后，共得到159只有效的股票数据，并提取其主要财务数据(如表1所示)，最终数据结果如下：

?

?

聚类分析：

?

标准化数据：由于不同指标数据间的量纲不一致，为了消除其影响，故将数据集进行Z标准化:

?

from sklearn import preprocessing as pc

?

X = pc.scale(X_basic[[‘x1‘,‘x2‘,‘x3‘,‘x4‘,‘x5‘,‘x6‘,

?

‘x7‘,‘x8‘,‘x9‘,‘x10‘,‘x11‘,‘x12‘,‘x13‘,‘x14‘,‘x15‘]])

?

为了更好地了解不同公司的特征财务状况，对标准化处理后的样本数据用上述15个财务指标为变量，采用KMeans聚类方法将公司聚成3类(自己指定),最终聚类结果如下：

?

from sklearn.cluster import KMeans

?

# 进行聚类，聚成 3 类
n = 3
estimator = KMeans(n_clusters=n)
estimator.fit(X)
label_pred = estimator.labels_
length = len(label_pred)

?

?

从上图知共分成三类，其中第三类含有博信股份、鲁亿通等5家企业，第一类含有奥飞娱乐，宝德股份等11家企业，最后均为第二类，也是企业最多的一个类。对各类指标的均值在EXCEL进行雷达图，其结果如下图所示：

从上述雷达图可以看出：从各类均值指标雷达直观看出：

这3类人工智能公司在指标x12至x15上均没有显著性差异，即企业间的发展能力指标基本相差不远，这与2017年人工智能在国内被重视和高速发展实际相符。第一类企业的营运能力是最好的，但是盈利能力介于其他两列中间。第二类是发展最稳定的企业群体，它的各项指标间波动均不大，是发展最稳定和盈利最好的企业，应该是属于稳定型企业。第三类企业的债偿能力指标是最好的，但是同时盈利能力指标是最差的。

?

因子分析：

1：进行因子分析，为验证样本的可行性，首先对样本数据进行KMO和巴特利特检验，结果如下表所示：

from factor_analyzer import factor_analyzer,Rotator

# KMO值

print round(factor_analyzer.calculate_kmo(X_basic)[1],5)

# 巴特利特球形度值

print round(factor_analyzer.calculate_bartlett_sphericity(X_basic)[1],5)

#0.63664
#0.0

从上图可知：KMO检验中，其KMO的概率值为0.63371,变量间的偏相关性较一般强，可以进行因子分析。巴特利特球形度检验检验中，

其伴随概率为0.00小于显著水平5%，说明这些变量可以进行因子分析。

?

2：初始化主成分分析参数：

# 主成分提取参数设置:
# ?n_factors = 4:4个主成分
# ?method = ‘principal‘：提取方法为主成分提取方法
# ?rotation = ‘varimax‘：旋转方法为最大方差法
fa = factor_analyzer.FactorAnalyzer(n_factors=4,method=‘principal‘,rotation=‘varimax‘)
fa.fit(X_basic)

?

3：因子分析中各主成分贡献率，以及累计贡献率结果如下表：

# 公因子方差
print fa.get_communalities()
print ‘n‘
# 贡献率
print fa.get_factor_variance()

###其中方差贡献率结果中：特征值，方差贡献率，累计方差贡献率

从上述的总方差解释中，可以看出总体主要成分的累积贡献度为70.92%，整体贡献度

算高。其次各主成分的贡献度占比依次是22.19%，20.83%，14.369%，13.54%。同时这4个特征根反映出来的信息量占全部信息量的70.92%高于70%，

满足要求，因此选取这4个特征根进行分析。

?

4：因子分析中的旋转后的因子载荷矩阵如下：

# 旋转后的因子载荷矩阵
rotator = Rotator()
print rotator.fit_transform(fa.loadings_)

从上表因子载荷矩阵中，可以知道：

（1）主成分1与x9,x10,x11,x14,x15间相关系数较大，主要体现的是企业的盈利能力指标，故可以称作为“盈利因子”

（2）主成分2与x1,x2,x3,4间相关系数较大，主要体现的是企业的债偿能力指标，故可以称作为“债偿因子”

（3）主成分3与x12,x13间相关系数较大，主要体现的是企业的发展能力指标，故可以称作为“发展因子”

（4）主成分4与x5,x6,x7,x8间相关系数较大，主要体现的是企业的营运能力指标，故可以称作为“营运因子”

?

?

5：成分得分系数矩阵：

?

# 成分得分系数矩阵
corr = np.matrix(corr)
load_matrix = np.matrix(load_matrix)
print np.dot(np.linalg.inv(corr),load_matrix)

?

?

6：各个主成分得分和综合得分：

?

F=(成分1*贡献率1+成分2*贡献率2+….+成分4*贡献率4)/累计贡献率：

?

?means=datas.mean()

stds=datas.std()
adict={}
for i in range(15):
????a=‘x‘+str(i+1)
????adict[i]=[means[a],stds[a]]
# 企业各主成分得分
def get_Fac_score(x,no,adict):
????arr = pr_matrix_score_t[no].tolist()[0]
????sum=0
????for i in range(15):
????????a=‘x‘+str(i+1)
????????sum=sum+(x[a]-adict[i][0])/(adict[i][1])*arr[i]
????return sum
datas[‘FAC1‘]= datas.apply(get_Fac_score,args=(0,adict),axis=1)
datas[‘FAC2‘]= datas.apply(get_Fac_score,args=(1,axis=1)
datas[‘FAC3‘]= datas.apply(get_Fac_score,args=(2,axis=1)
datas[‘FAC4‘]= datas.apply(get_Fac_score,args=(3,axis=1)
datas[‘FAC‘] = datas.apply(lambda x:(x[‘FAC1‘]*var[1][0]+x[‘FAC2‘]*var[1][1]+
???????????x[‘FAC3‘] * var[1][2]+x[‘FAC4‘]*var[1][3])/var[2][3],axis=1)
datas.to_excel(‘result_fac.xls‘,index=False)

结论：

其中类一的综合排名均靠前，说明企业的营运能力指标在企业排名中最为重要，其次类二的排名居中间位置，排名整体靠后的是类三，其中它的企业盈利能力指标是最差的。从本次因子分析中可以知道：企业的营运能力和盈利能力对企业排名有着明显的影响。

?

?

?

?

python 代码：

# kmodatas = pd.read_excel(‘datas.xls‘)datas_1 = pd.read_excel(‘stone_class.xls‘)datas=pd.merge(datas,datas_1)datas_1=datas.copy()# 读取指标X_basic = datas[[‘x1‘,‘x2‘,‘x3‘,‘x4‘,‘x5‘,‘x6‘,‘x7‘,‘x8‘,‘x9‘,‘x10‘,‘x11‘,‘x12‘,‘x13‘,‘x14‘,‘x15‘]]# 相关系数corr = X_basic.corr().valuesX_basic = pc.scale(    X_basic[[‘x1‘,‘x15‘]])def kmo():    # 相关系数    corr = X_basic.corr().values    # print corr    # 逆矩阵    corr_inv = np.linalg.inv(corr)    # 对角线取倒数，其他为0    S2 = np.diag(corr_inv)    S2 = np.linalg.inv(np.diag(S2))    # 反映像协方差矩阵    AIS = np.dot(S2,corr_inv)    AIS = np.dot(AIS,S2)    # 是映像协方差矩阵    IS = corr + AIS - 2 * S2    # 将矩阵AIS对角线上的元素开平方，并且将其余元素都变成0    Dai = np.diag(np.sqrt(np.diag(AIS)))    # IR是映像相关矩阵    IR = np.dot(np.linalg.inv(Dai),IS)    IR = np.dot(IR,np.linalg.inv(Dai))    # AIR是反映像相关矩阵    AIR = np.dot(np.linalg.inv(Dai),AIS)    AIR = np.dot(AIR,np.linalg.inv(Dai))    a = np.power((AIR - np.diag(np.diag(AIR))),2)    a = np.sum(a,axis=1)    AA = np.sum(a)    b = corr - np.identity(15)    b = np.power(b,2)    b = np.sum(b,axis=1)    BB = np.sum(b)    MSA = b / (b + a)    AIR = AIR - np.identity(15) + np.diag(MSA)    kmo = BB / (AA + BB)    return kmodef bart():    corr = X_basic.corr().values    # 计算结果有问题    # bart=st.bartlett(*corr)    detCorr = np.linalg.det(corr)    n = len(X_basic)    p = len(X_basic.columns)    statistic = -math.log(detCorr) * (n - 1 - (2 * p + 5) / 6)    df = p * (p - 1) / 2    # 双侧概率    pval = (1.0 - st.chi2.cdf(statistic,df)) * 2    return statistic,pval# kmo_result =  kmo()# barr_result = bart()# print kmo_result,barr_result#0.636636281781# (2129.4863886695084,0.0)# KMO值print round(factor_analyzer.calculate_kmo(X_basic)[1],5)# 巴特利特球形度值print round(factor_analyzer.calculate_bartlett_sphericity(X_basic)[1],5)# 主成分提取参数设置:#  n_factors = 4:4个主成分#  method = ‘principal‘：提取方法为主成分提取方法#  rotation = ‘varimax‘：旋转方法为最大方差法fa = factor_analyzer.FactorAnalyzer(    n_factors=4,method=‘principal‘,rotation=‘varimax‘)fa.fit(X_basic)# print ‘n‘# # 公因子方差# print fa.get_communalities()# print ‘n‘# # 贡献率var = fa.get_factor_variance()print var[0]print var[1]print var[2]# print ‘n‘# # 旋转后的因子载荷矩阵rotator = Rotator()load_matrix = rotator.fit_transform(fa.loadings_)# 成分得分系数矩阵corr = np.matrix(corr)load_matrix = np.matrix(load_matrix)pr_matrix_score = np.dot(np.linalg.inv(corr),load_matrix)pr_matrix_score_t = pr_matrix_score.Tmeans=datas.mean()stds=datas.std()adict={}for i in range(15):    a=‘x‘+str(i+1)    adict[i]=[means[a],stds[a]]# 企业各主成分得分def get_Fac_score(x,adict):    arr = pr_matrix_score_t[no].tolist()[0]    sum=0    for i in range(15):        a=‘x‘+str(i+1)        sum=sum+(x[a]-adict[i][0])/(adict[i][1])*arr[i]    return sumdatas[‘FAC1‘]= datas.apply(get_Fac_score,axis=1)datas[‘FAC2‘]= datas.apply(get_Fac_score,axis=1)datas[‘FAC3‘]= datas.apply(get_Fac_score,axis=1)datas[‘FAC4‘]= datas.apply(get_Fac_score,axis=1)datas[‘FAC‘] = datas.apply(lambda x:(x[‘FAC1‘]*var[1][0]+x[‘FAC2‘]*var[1][1]+                                    x[‘FAC3‘] * var[1][2]+x[‘FAC4‘]*var[1][3])/var[2][3],axis=1)datas.to_excel(‘result_fac.xls‘,index=False)