<jupyter_start><jupyter_text>III - Неопределенность оптимального портфеля. Оптимизация CVaR. Целью работы является оценка неопределенности оптимального портфеля для нормального многомерного распределения и двух способов вычисления оптимального портфеля (оптимизация в модели Марковица и оптимизация CVaR)  Подготовка модели  Загружаем тикеры с фоднового рынка NASDAQ<jupyter_code>import pandas as pd
import numpy as np

DATA_PATH    = './downloader-data'
TICKERS_PATH = './tickers'

def read_tickers(stock_markets) -> pd.DataFrame:
    '''Returns pandas dataframe containing all tickers for the @stock_markets '''
    ticker_files = [f'{TICKERS_PATH}/{sm}.csv' for sm in stock_markets]
    tickers = pd.concat([pd.read_csv(tf) for tf in ticker_files], ignore_index=True)
    return tickers


stock_markets = ['NASDAQ']
tickers = read_tickers(stock_markets)
tickers<jupyter_output><empty_output><jupyter_text>Загружаем исторические данные за 2021 год для полученных тикеров<jupyter_code>def read_historical_data(tickers):
    '''Returns dict {ticker : historical data for the ticker}'''
    data_for_ticker = {}
    for index, (ticker, name) in tickers.iterrows():
        try:
            data = pd.read_csv(f'{DATA_PATH}/{ticker}.csv')
            if len(data) > 100:
                data_for_ticker[ticker] = data
        except:
            pass
    return data_for_ticker


data_for_ticker = read_historical_data(tickers)<jupyter_output><empty_output><jupyter_text>Пример данных для AAPL:<jupyter_code>data_for_ticker['AAPL']<jupyter_output><empty_output><jupyter_text>Выбираем 20 активов<jupyter_code>assets = ['MDLZ', 'MSFT', 'NXPI', 'PCAR', 'INTC',
          'NVDA', 'ILMN', 'DXCM', 'ROST', 'LULU',
          'PRFT', 'FTNT', 'WIRE', 'NFLX', 'GAIN',
          'SYNA', 'NYMTP', 'ANAT', 'TRNS', 'SFBS']

assets_data = {}
for asset in assets:
    assets_data[asset] = data_for_ticker[asset]<jupyter_output><empty_output><jupyter_text>Оцениваем математические ожидания доходностей<jupyter_code>import pandas as pd
import numpy as np

def add_logret(ticker_data: pd.DataFrame, by_column='Close') -> pd.DataFrame:
    '''Returns @ticker_data with calculated "logret" and "-logret" columns'''
    ticker_data = ticker_data.assign(logret=np.log(ticker_data[by_column]).diff())
    ticker_data['-logret'] = ticker_data['logret'].mul(-1)
    return ticker_data


def get_logret_mean_std(tickers, data_map, by_column='Close') -> pd.DataFrame:
    '''Returns @result pd.DataFrame such that @result.loc[ticker] == [logret_mean, logret_std]'''
    result = pd.DataFrame(data=[], columns=['ticker', 'logret_mean', 'logret_std'])
    result.set_index('ticker', inplace=True)

    for ticker in tickers:
        ticker_data = data_map[ticker]
        ticker_logret = np.log(ticker_data[by_column]).diff()
        result.loc[ticker] = [ticker_logret.mean(), ticker_logret.std()]
    
    return result
  

for ticker in assets_data.keys():
    assets_data[ticker] = add_logret(assets_data[ticker])

estims = get_logret_mean_std(assets, assets_data)
estims<jupyter_output><empty_output><jupyter_text>Находим матрицу выборочных ковариаций доходностей<jupyter_code>def get_covariation_matrix(tickers,                          
                           data_map=data_for_ticker,
                           by_column='logret'):
    columns = pd.DataFrame()
    for ticker in tickers:
        columns[ticker] = data_map[ticker][by_column][1:]
        
    matrix = columns.corr()
    return matrix

cov_matrix = get_covariation_matrix(assets, data_map=assets_data)
cov_matrix
import plotly.express as px

fig = px.imshow(cov_matrix, title='Correlation Matrix',
                color_continuous_scale=px.colors.diverging.RdYlGn,
                zmin=-1, zmax=1)
fig.show()<jupyter_output><empty_output><jupyter_text>Проверяем вырожденность матрицы и число обусловленности<jupyter_code>det = np.linalg.det(cov_matrix)
print(f'det = {det}')

con = np.linalg.cond(cov_matrix, p='fro')
print(f'con = {con}')<jupyter_output>det = 0.002219356005352476
con = 55.688725974909254<jupyter_text>1. Истинный оптимальный портфель в модели Марковица с заданным отношением к риску.  С заданным отношением к риску  подобираем константу b таким образом, что истинный оптимальный CVaR портфель совпадает с истинным оптимальным портфелем п.1. Значение константы смотри в упражнениях к теме.$$ b = \frac{1}{\sqrt{2 \pi }} \frac{1}{(1 - \beta)} exp(-(\Phi^{-1}(\beta))^2 / 2)$$<jupyter_code>import math
from scipy.stats import norm

beta = 0.95
b = (1 / math.sqrt(2 * math.pi)) * (1 / (1 - beta)) * np.exp(-((norm.ppf(beta))**(-1))**2 / 2) 
print(f'b = {b}')<jupyter_output>b = 6.632540979594675<jupyter_text>Решаем задачу оптимизации $$ -E(x)+ с \cdot σ(x) \rightarrow min $$при условиях:$$ x_1 + x_2 + ... +x_N = 1 $$$$ x_i >= 0 $$где$$ E(x)= E_1x_1 + E_2 x_2 + ... + E_Nx_N $$$$ σ^2(x)=\Sigma\Sigma σ_{i,j} x_i \cdot x_j $$<jupyter_code>import math 
import numpy as np
from scipy.optimize import minimize

def get_E(x, vector_E):
    E = sum([(E_i * x_i) for E_i, x_i in zip(vector_E, x)]) 
    return E


def get_sigma(x, matrix_cov):
    sigma_squared = 0
    
    for i in range(len(x)):
        for j in range(len(x)):
            simga_i_j = matrix_cov.iloc[i].iloc[j]
            sigma_squared += simga_i_j * x[i] * x[j]
    
    sigma = math.sqrt(sigma_squared)
    return sigma
            
            
def target_function(x, vector_E, matrix_cov, minusE):
    E = get_E(x, vector_E)
    sigma = get_sigma(x, matrix_cov)
    if minusE:
        result = -E + b * sigma
    else:
        result = E + b * sigma
    return result


def find_optimal(vector_E=estims['logret_mean'], matrix_cov=cov_matrix, minusE=True):
    x0 = np.array([1/len(assets)] * len(assets))
    solution = minimize(target_function, x0, args=(vector_E, matrix_cov, minusE), 
                        method='SLSQP', 
                        constraints=[{'type': 'eq',  'fun': lambda x: sum(x) - 1}],
                        bounds=[(0, 1)] * len(assets))
    if not solution.success:
        raise Exception(opt.message)
    return solution
    
true_opt_solution = find_optimal()<jupyter_output><empty_output><jupyter_text>Веса портфеля:<jupyter_code>import plotly.express as px

print(f'Сумма весов: {sum(true_opt_solution.x)}')
fig = px.bar(x=assets, y=true_opt_solution.x)
fig.show()<jupyter_output>Сумма весов: 1.0000000000000002<jupyter_text>Значение целевой функции:<jupyter_code>print(true_opt_solution.fun)
def get_portfolio_coordinates(weights,
                              _tickers,
                              _logret_mean,
                              _logret_std):
    portfolio_logret_mean = 0
    portfolio_logret_std = 0
     
    for i in range(len(_tickers)):
        portfolio_logret_mean += weights[i] * _logret_mean[i]
        for j in range(len(_tickers)):
            portfolio_logret_std += weights[i] * weights[j] * _logret_std.loc[_tickers[i], _tickers[j]]
            
    portfolio_logret_std = math.sqrt(portfolio_logret_std)
    return portfolio_logret_std, portfolio_logret_mean 
import plotly.express as px
from sklearn.preprocessing import normalize

ticker_colour = estims['logret_mean'] / estims['logret_std']
ticker_size = (ticker_colour - ticker_colour.min()) / ticker_colour.max()
fig = px.scatter(estims, x='logret_std', y='logret_mean',
                 hover_name=estims.index, 
                 color_continuous_scale=px.colors.diverging.RdYlGn,
                 color=ticker_colour,
                 size=ticker_size)
fig.show()
# _x, _y = get_portfolio_coordinates(true_opt_solution.x, assets, estims['logret_mean'], estims['logret_std'])<jupyter_output><empty_output><jupyter_text>2. Оценка неопределенности оптимального портфеля в модели Марковица с заданным отн. к риску.   2.1 Задаём число наблюдений T=30. С помощью генератора многомерного нормального распределения создаём выборку размера Т из нормального распределения с вектором математических ожиданий  E=(E1, E2, …, EN) и матрицей ковариаций (σi,j).<jupyter_code>T = 30
sample_raw = np.random.multivariate_normal(estims['logret_mean'], cov_matrix, T)
sample = pd.DataFrame(columns=assets, data=sample_raw)
sample<jupyter_output><empty_output><jupyter_text>2.2 По построенной выборке делаем оценку Eest вектора математических ожиданий ...<jupyter_code>estE = sample.mean()
estE<jupyter_output><empty_output><jupyter_text>... и оценку (σesti,j) матрицы ковариаций.<jupyter_code>estCov = sample.corr()
estCov
import plotly.express as px

fig = px.imshow(estCov, title='Correlation Matrix',
                color_continuous_scale=px.colors.diverging.RdYlGn, 
                zmin=-1, zmax=1)
fig.show()<jupyter_output><empty_output><jupyter_text>2.3 Используя эти оценки решаем задачу оптимизации<jupyter_code>opt_solution = find_optimal(estE, estCov)<jupyter_output><empty_output><jupyter_text>Находим и фиксируем веса портфеля ...<jupyter_code>print(f'Сумма весов: {sum(opt_solution.x)}')
fig = px.bar(x=assets, y=opt_solution.x)
fig.show()<jupyter_output>Сумма весов: 1.0000000000000029<jupyter_text>и значение целевой функции<jupyter_code>print(opt_solution.fun)<jupyter_output>2.918238358413867<jupyter_text>2.4 Сравниаем два портфеля: истинный (п.1) и выборочный (п.2.3). Оцениваем относительную ошибку в определении весов портфеля в норме Manhattan (L1 норма Минковского). Делайте выводы и сравнение в системе координат (σ, E).<jupyter_code>from scipy.spatial.distance import cityblock
cityblock(true_opt_solution.x, opt_solution.x)<jupyter_output><empty_output><jupyter_text>2.5 Повторите эксперимент S=40 раз и оцените среднюю относительную ошибку по S повторениям эксперимента. Сделайте выводы.  Сделайте сравнение в системе координат (σ, E).<jupyter_code>import copy

S = 40
errors = []
x_vectors = pd.DataFrame(columns = [f'x_{i}' for i in range(len(assets))], data=[])


for iteration in range(0, S):
    T = 30
    sample_raw = np.random.multivariate_normal(estims['logret_mean'], cov_matrix, T)
    sample = pd.DataFrame(columns=assets, data=sample_raw)
    
    estE = sample.mean()
    estCov = sample.cov()
    
    opt_solution = find_optimal(estE, estCov)
        
    errors.append(cityblock(true_opt_solution.x, opt_solution.x))
    x_vectors.loc[iteration] = (copy.deepcopy(opt_solution.x))

    
print(f'Mean Erorr = {np.mean(errors)}')
print(f'Mean X:')
x_vectors.mean()
import plotly.express as px

fig = px.bar(x=assets, y=x_vectors.mean(), title='Mean X')
fig.show()<jupyter_output><empty_output><jupyter_text>2.6  Предположите, что нам известны точные значения математических ожиданий E=(E1, E2, …, EN). Повторяем пп. 2.2-2.5. используя оценку только матрицы ковариаций<jupyter_code>import copy

S = 40
errors = []
x_vectors = pd.DataFrame(columns = [f'x_{i}' for i in range(len(assets))], data=[])


for iteration in range(0, S):
    T = 30
    sample_raw = np.random.multivariate_normal(estims['logret_mean'], cov_matrix, T)
    sample = pd.DataFrame(columns=assets, data=sample_raw)
    
    estE = sample.mean()
    estCov = sample.cov()
    
    opt_solution = find_optimal(estims['logret_mean'], estCov)
        
    errors.append(cityblock(true_opt_solution.x, opt_solution.x))
    x_vectors.loc[iteration] = (copy.deepcopy(opt_solution.x))

    
print(f'Mean Erorr = {np.mean(errors)}')
print(f'Mean X:')
x_vectors.mean()
import plotly.express as px

fig = px.bar(x=assets, y=x_vectors.mean(), title='Mean X')
fig.show()<jupyter_output><empty_output><jupyter_text>3. Оценка неопределенности оптимального CVaR портфеля С заданным отношением к риску  подобираем константу b таким образом, что истинный оптимальный CVaR портфель совпадает с истинным оптимальным портфелем п.1. Значение константы смотри в упражнениях к теме.$$ b = \frac{1}{\sqrt{2 \pi }} \frac{1}{(1 - \beta)} exp(-(\Phi^{-1}(\beta))^2 / 2)$$<jupyter_code>beta = 0.95
b = (1 / math.sqrt(2 * math.pi)) * (1 / (1 - beta)) * np.exp(-((norm.ppf(beta))**(-1))**2 / 2) 
print(f'b = {b}')
def get_E(x, vector_E):
    E = sum([(E_i * x_i) for E_i, x_i in zip(vector_E, x)]) 
    return E


def get_sigma(x, matrix_cov):
    sigma_squared = 0
    
    for i in range(len(x)):
        for j in range(len(x)):
            simga_i_j = matrix_cov.iloc[i].iloc[j]
            sigma_squared += simga_i_j * x[i] * x[j]
    
    sigma = math.sqrt(sigma_squared)
    return sigma
            
            
def target_function_CVaR(x, vector_E, matrix_cov):
    E = get_E(x, vector_E)
    sigma = get_sigma(x, matrix_cov)
    result = E + b * sigma
    return result


def find_optimal_CVaR(vector_E, matrix_cov):
    x0 = np.array([1/len(assets)] * len(assets))
    solution = minimize(target_function_CVaR, x0, args=(vector_E, matrix_cov), 
                        method='SLSQP', 
                        constraints=[{'type': 'eq',  'fun': lambda x: sum(x) - 1}],
                        bounds=[(0, 1)] * len(assets))
    if not solution.success:
        raise Exception(opt.message)
    return solution<jupyter_output><empty_output><jupyter_text>Используя оценки из 2.2 решаем задачу оптимизации<jupyter_code>opt_solution_CVaR = find_optimal_CVaR(estE, estCov)<jupyter_output><empty_output><jupyter_text>Находим и фиксируем веса портфеля ...<jupyter_code>print(f'Сумма весов: {sum(opt_solution_CVaR.x)}')
fig = px.bar(x=assets, y=opt_solution_CVaR.x)
fig.show()<jupyter_output><empty_output><jupyter_text>и значение целевой функции<jupyter_code>print(opt_solution_CVaR.fun)<jupyter_output><empty_output><jupyter_text>3.2 Сравниаем два портфеля: истинный (п.1) и найденный в п.3.1. Оцениваем относительную ошибку в определении весов портфеля в норме Manhattan (L1 норма Минковского). Сравниваем с ошибкой портфеля из п. 2.3<jupyter_code>er1 = cityblock(true_opt_solution.x, opt_solution.x)
er2 = cityblock(true_opt_solution.x, opt_solution_CVaR.x)
print(f'Ошибка оптимального VaR портфеля(п.2.3) = {er1}')
print(f'Ошибка оптимального CVaR портфеля = {er2}')<jupyter_output><empty_output><jupyter_text>3.3 Повторите эксперимент S=40 раз и оцените среднюю относительную ошибку по S повторениям эксперимента. Сделайте выводы.  Сравните с ошибкой из п. 2.5.<jupyter_code>import copy

S = 40
errors = []
x_vectors = pd.DataFrame(columns = [f'x_{i}' for i in range(len(assets))], data=[])


for iteration in range(0, S):
    T = 30
    sample_raw = np.random.multivariate_normal(estims['logret_mean'], cov_matrix, T)
    sample = pd.DataFrame(columns=assets, data=sample_raw)
    
    estE = sample.mean()
    estCov = sample.cov()
    
    opt_solution_CVaR = find_optimal_CVaR(estE, estCov)
        
    errors.append(cityblock(true_opt_solution.x, opt_solution_CVaR.x))
    x_vectors.loc[iteration] = (copy.deepcopy(opt_solution_CVaR.x))

    
print(f'Mean Erorr = {np.mean(errors)}')
print(f'Mean X:')
x_vectors.mean()
import plotly.express as px

fig = px.bar(x=assets, y=x_vectors.mean(), title='Mean X')
fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12450446818517771, 'fraction_numerical': 0.01666330712893906, 'mean_word_length': 4.150173010380623, 'pattern_counts': {'":': 0, '<': 100, '<?xml version=': 0, '>': 105, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '273071', 'license': 'MIT', 'n_tokens_mistral': 6044, 'n_tokens_neox': 6071, 'n_words': 1108, 'path': 'Task4.ipynb', 'repo_name': 'aailyinh/fin-ops', 'size': '14883'}