事前準備として、ノイズ除去で使用する関数をモジュール化しておく。
　ここでは、MarPat.pyとしてある。

　yahoofinanceから、100のリターン行列を作り、ノイズ除去無し、ノイズ除去あり、ノイズ除去とデトーニングありの各場合の効率的フロンティアを計算する。

Yahoo Financeからticker symbolを使い、終値を読み込み、リターン行列を作成する

import numpy as np
import pandas as pd
from datetime import datetime
from pandas_datareader import data as pdr
import yfinance as yfin

yfin.pdr_override()

ticker_symbols=['AAPL','ABNB','ABSI','ACN','ADBE','ADI','ADP','ADSK','AEP','AFG',
                'ALGN','AMAT','AMD','AMGN','AMZN','ANSS','AON','ASML','AVGO','AXP','AZN',
                'BA','BABA','BAC','BKNG','BKR','BP','BUD','BX','C','CAT','CDNS','CEG','CHTR',
                'CMCSA','COST','CPRT','CRM','CSCO','CSGP','CSX','CTAS','CTSH','CVX','DLTR',
                'DXCM','EBAY','EXC','FANG','FAST','FTNT','GEHC','GOOG','GOOGL','IBM',
                'IDXX','ILMN','INTC','INTU','ISRG','JD','KDP','KHC','KLAC','LCID',
                'LULU','MAR','MCHP','MDLZ','MELI','META','MRNA','MRVL','MSFT','NFLX',
                'NXPI','NVDA','ODFL','ORLY','PANW','PAYX','PCAR','PYPL','QCOM','REGN','POST',
                'SGEN','SIRI','SNPS','TMUS','TTD','TXN','VRSK','VRTX','WBA','WBD','WDAY',
                'WEL','ZM','ZS']
    
startdate= datetime(2023,1,1)
enddate = datetime(2023,8,1)

s_data= pdr.get_data_yahoo(ticker_symbols, startdate, enddate,interval='1h')['Adj Close']

return_month = s_data.pct_change()
return_mat=return_month.dropna()

　intervalは、一日が1d、1時間が1h、30分が30m、15分が15mで変えられる。

ノイズ除去のために、Marchenko-Patstur分布の当てはめを行う。

　リターン行列の平均ベクトり、分散ベクトル、標準偏差、共分散行列、相関行列をとる。

return_mean=return_mat.mean()
return_var=return_mat.var()
return_std=return_mat.std()
return_cov=return_mat.cov()
return_corr=return_mat.corr()

import MarPat as MP
import matplotlib.pyplot as plt
from scipy.optimize import minimize

q=return_month.shape[0]/return_month.shape[1]
eVal0,eVec0=MP.getPCA(return_corr)
bw=MP.findOptimalBWidth(np.diag(eVal0))
bW=bw['bandwidth']

eMax0,var0=MP.findMaxEval(np.diag(eVal0),q,bWidth=bW)
nFact0=eVal0.shape[0]-np.diag(eVal0)[::-1].searchsorted(eMax0)

効率的フロンティアを計算する

空売り禁止として、ポートフォリオの重みを$${(0,1)}$$とする。

 #ポートフォリオの期待収益率の分散 
def min_func_var(weights):
    return np.dot(weights.T, np.dot(cov, weights))

from scipy.optimize import minimize

bWidth=bW

mean_high = return_mean.max().round(3)
mean_low = return_mean.min().round(3)
trets = np.linspace(mean_low, mean_high, 50)

n_assets =return_month.shape[1]
x0 = [1. / n_assets] * n_assets
bounds = tuple((0,1) for i in range(n_assets))

tvols0 = []
cov=return_cov
for tret in trets:
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1},
               {'type': 'eq', 'fun': lambda x: np.sum(return_mean*x) - tret})
    res = minimize(fun=min_func_var, x0=x0, method='SLSQP', bounds=bounds,constraints=constraints)
    tvols0.append(np.sqrt(res['fun']))
tvols0 = np.array(tvols0)


tvols =[]
cov=MP.deNoiseCov(return_cov, q,bWidth)
for tret in trets:
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1},
               {'type': 'eq', 'fun': lambda x: np.sum(return_mean*x) - tret})
    res = minimize(fun=min_func_var, x0=x0, method='SLSQP', bounds=bounds,constraints=constraints)
    tvols.append(np.sqrt(res['fun']))
    
    
tvols = np.array(tvols)

　次にノイズ除去した共分散行列からデトーニングで市場成分を抜いた行列は得意行列となっているので、ポートフォリオ計算では次元を元に戻す必要があることに注意する。

dn_corr=MP.cov2corr(cov)
eVal1,eVec1=MP.getPCA(dn_corr)
dt_corr=MP.detonedCorr(dn_corr,eVal1, eVec1)
dt_cov=MP.corr2cov(dt_corr,np.diag(cov)**.5)
cov=dt_cov
weights = []


for tret in trets:
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1},
               {'type': 'eq', 'fun': lambda x: np.sum(return_mean*x) - tret})
    res = minimize(fun=min_func_var, x0=x0, method='SLSQP', bounds=bounds,constraints=constraints)
    weights.append(res['x'])

weights=np.array(weights)
weights.shape[0]
eVal1,eVec1=MP.getPCA(dt_corr)
tvols1=[]

for i in range(weights.shape[0]):
    weight=np.dot(eVec1,weights[i])
    tvols1.append(np.sqrt(min_func_var(weight)))
    

tvols1 = np.array(tvols1)

グラフ化

import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.scatter(tvols0,trets,c=trets/tvols0, marker='o',label='cov')
plt.scatter(tvols, trets, c=trets/tvols, marker='x',label='denoised cov')
plt.scatter(tvols1,trets,c=trets/tvols1, marker='*',label='detoned denoised cov')
plt.grid(True)
plt.ylim([0.0001,None])
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.legend()