[Python] Inverse Covariance를 통한 수익률 유사도 측정

import FinanceDataReader as fdr
fdr.__version__

import numpy as np
import pandas as pd
from scipy import linalg

from datetime import datetime
import pytz

from sklearn.datasets import make_sparse_spd_matrix
from sklearn.covariance import GraphLassoCV, ledoit_wolf
from sklearn.preprocessing import StandardScaler
from sklearn import cluster, manifold

import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

np.random.seed(0)

# 한국거래소 상장종목 전체
# 용도 : 코드와 종목명 가져오기
df_krx = fdr.StockListing('KRX')

# 코스피 종목 추출
# 용도 : DataReader 사용시 종목별 기간 조회 건수가 상이하다.
#        코스피 데이터 날짜를 기준으로 가져오기 위해 사용
strt_dt = '2018-07-01' # 시작일 지정
end_dt = '2018-12-31'  # 종료일 지정
kospi_df = fdr.DataReader('KS11', strt_dt,end_dt)
kospi_df.reset_index(inplace = True)

def stock_reader(kospi_df, code_list, df='p', n=0):

    if n == 0:
        n = len(code_list)

    stock_df = pd.DataFrame()
    stock_df['Date'] = kospi_df['Date']
    print("동기간 KOSPI 생성일수 : ", len(kospi_df['Date']))
    normal_cnt = 0
    err_cnt = 0
    code_nm_list = []
    symbol_list = []

    for code in code_list:
        stock = df_krx[df_krx.Symbol == code]
        code_nm = list(stock.Name)[0]
        symbol = list(stock.Symbol)[0]


        try:
            temp = fdr.DataReader(code, strt_dt, end_dt)

            # 데이터일수가 시장보다 작으면 skip(최근 상장 데이터로 판단)
            if len(temp) < len(kospi_df['Date']):
                err_cnt += 1
                print("skip : (",err_cnt,")", code, code_nm, strt_dt, end_dt, ", 건수 : ", len(temp))
                continue

            temp.reset_index(inplace = True)
            temp_df = pd.merge(temp[['Date','Close', 'Open']], kospi_df[['Date']], on='Date', how='right')

            if df == 'v':
                stock_df[code_nm] = temp_df.Close - temp_df.Open #variation
            elif df == 'p':
                stock_df[code_nm] = temp_df.Close # price

            normal_cnt += 1
            code_nm_list.append(code_nm)
            symbol_list.append(symbol)
            print("정상 : (",normal_cnt,")", code, code_nm, strt_dt, end_dt, ", 건수 : ", len(temp), "->", len(stock_df))
        except:
            err_cnt += 1
            print("skip : (",err_cnt,")", code, code_nm, strt_dt, end_dt, ", 건수 : ", len(temp))

        if normal_cnt == n:
            print('총', n,'개 생성 설정 / ', normal_cnt, '개 생성 완료')
            break # n개 종목 생성시 종료

        # 데이터 정렬
        stock_df.sort_values('Date', ascending=True, inplace=True) # ascending=True 오름차순, False 내림차순
        # 결측치 채우기
        stock_df.fillna(method='ffill', inplace=True)

    return stock_df, code_nm_list, symbol_list

n = 100 # 생성할 종목수 지정
df_krx_list = df_krx['Symbol'].head(n*2) # 임시로 2배까지 루프

stock_df, code_nm_list, symbol_list = stock_reader(kospi_df, df_krx_list, 'v', n) # v :variation, p : price(close)

stock_df.dropna(inplace=True)
X = stock_df.set_index('Date')
X /= X.std(axis=0) # 표준화

#estimate inverse covariance
graph = GraphLassoCV(verbose=10) # 그리드서치를 통해 유의한 alphas가 있는 경우 지정
graph.fit(X)

gl_cov = graph.covariance_
gl_prec = graph.precision_
gl_alphas =graph.cv_alphas_
gl_scores = np.mean(graph.grid_scores, axis=1)

plt.figure()
sns.heatmap(gl_prec) # vmin=-2, vmax=2

plt.figure()
plt.plot(gl_alphas, gl_scores, marker='o', color='b', lw=2.0, label='GraphLassoCV')
plt.title("Graph Lasso Alpha Selection")
plt.xlabel("alpha")
plt.ylabel("score")
plt.legend()

names = np.array(code_nm_list)
symbols = np.array(symbol_list)

#cluster using affinity propagation
_, labels = cluster.affinity_propagation(gl_cov) # 공분산 데이터로 클러스터링
num_labels = np.max(labels)

for i in range(num_labels+1):
    print("Cluster %i: %s" %((i+1), ', '.join(names[labels==i])))

node_model = manifold.LocallyLinearEmbedding(n_components=2, n_neighbors=6, eigen_solver='dense')
embedding = node_model.fit_transform(X.T).T  # 일일 가격차 표준화값( (close-open)/std)을 기준으로 embedding

select_label = 14 # 라벨 사용자 지정

plot_temp = pd.DataFrame()
plot_temp['symbols'] = list(symbols)
plot_temp['names'] = list(names)
plot_temp['gl_prec'] = list(gl_prec)
plot_temp['em0'] = list(embedding[0])
plot_temp['em1'] = list(embedding[1])
plot_temp['labels'] = list(labels)
result_df = plot_temp
plot_temp = plot_temp[plot_temp['labels'] == select_label]
# 임베딩 변경
em0_array = plot_temp['em0']
em1_array = plot_temp['em1']
embedding_temp = np.vstack([em0_array, em1_array])

group_df, group_names, group_symbols = stock_reader(kospi_df, list(plot_temp['symbols']) , 'p', n) # v :variation, p : price(close)
group_df.set_index('Date', inplace=True)

group_df.plot(figsize=(15,7))
plt.ylabel('Price');
plt.xlabel('Time');
plt.show()

prec_temp = pd.DataFrame()
final_df = pd.DataFrame()
total_prec = ()

for x in range(0,len(plot_temp)):

    prec_temp['gl_prec'] =list(plot_temp['gl_prec'].iloc[x])
    prec_temp['labels'] = list(labels)
    final_df = prec_temp[prec_temp['labels'] == select_label]
    if len(total_prec) == 0:
        total_prec = np.array(final_df['gl_prec'])
    else:
        total_prec = np.vstack([total_prec, np.array(final_df['gl_prec'])])

# 일일 가격변화 표준값 X = (Close - Open) / std
# 클러스터링 그룹(labels) : X의 공분산값을 기준으로 음의 제곱 유클리드 거리 계산, (GraphLassoCV, affinity_propagation)
# 종목의 위치(embedding) : 일별 가격변동의 표준화값을 기준으로 2차원 embedding (LocallyLinearEmbedding)
# 연결강도(partial_corr) : X의 정밀도값을 기준으로 결정(그룹내에서 실제 유사도가 얼마나 큰지) (GraphLassoCV)

plt.figure()
plt.clf()
ax = plt.axes([0.,0.,1.,1.])
plt.axis('off')

partial_corr = total_prec
d = 1 / np.sqrt(np.diag(partial_corr))    
non_zero = (np.abs(np.triu(partial_corr, k=1)) > 0.02)  #connectivity matrix

#plot the nodes
plt.scatter(plot_temp['em0'],plot_temp['em1'], s = 100*d**2, c = plot_temp['labels'], cmap = plt.cm.get_cmap("Spectral"))

#plot the edges
start_idx, end_idx = np.where(non_zero)
# segments = [[np.array(plot_temp['em0'].iloc[:start])[0], np.array(plot_temp['em1'].iloc[:stop])[0]] for start, stop in zip(start_idx, end_idx)]
segments = [[embedding_temp[:,start], embedding_temp[:,stop]] for start, stop in zip(start_idx, end_idx)]
values = np.abs(partial_corr[non_zero])
lc = LineCollection(segments, zorder=0, cmap=plt.cm.hot_r, norm=plt.Normalize(0,0.7*values.max()))
lc.set_array(values)
lc.set_linewidths(5*values)
ax.add_collection(lc)

#plot the labels
for index, (x_names, x_labels, (x,y)) in enumerate(zip(np.array(plot_temp['names']),np.array(plot_temp['labels']), np.array(plot_temp[['em0','em1']]))):
    plt.text(x,y,x_names,size=12)

#generate plots
plt.figure()
plt.clf()
ax = plt.axes([0.,0.,1.,1.])
plt.axis('off')

partial_corr = gl_prec
d = 1 / np.sqrt(np.diag(partial_corr))    
non_zero = (np.abs(np.triu(partial_corr, k=1)) > 0.02)  #connectivity matrix

#plot the nodes
plt.scatter(embedding[0], embedding[1], s = 100*d**2, c = labels, cmap = plt.cm.get_cmap("Spectral"))

#plot the edges
start_idx, end_idx = np.where(non_zero)
segments = [[embedding[:,start], embedding[:,stop]] for start, stop in zip(start_idx, end_idx)]
values = np.abs(partial_corr[non_zero])
lc = LineCollection(segments, zorder=0, cmap=plt.cm.hot_r, norm=plt.Normalize(0,0.7*values.max()))
lc.set_array(values)
lc.set_linewidths(5*values)
ax.add_collection(lc)

#plot the labels
# for index, (name, label, (x,y)) in enumerate(zip(names, labels, embedding.T)):
#     plt.text(x,y,name,size=12)