Machine Learning Guru

import numpy as np
import pandas as pd
from datetime import datetime, timedelta

# Создание массива цен акций за последние 30 дней
initial_price = 100.0
daily_returns = np.random.normal(0.0005, 0.01, 30)  # Среднее значение 0.05%, стандартное отклонение 1%
prices = initial_price * np.cumprod(1 + daily_returns)

# Создание массива дат для временного ряда
end_date = datetime.now().date()
start_date = end_date - timedelta(days=29)
dates = np.array([start_date + timedelta(days=i) for i in range(30)])

# Создание структурированного массива для хранения данных OHLCV
dtype = [('date', 'datetime64[D]'), ('open', 'f8'), ('high', 'f8'), 
         ('low', 'f8'), ('close', 'f8'), ('volume', 'i4')]

# Генерация случайных данных OHLCV
n_days = 30
opens = prices
closes = opens * (1 + np.random.normal(0, 0.005, n_days))
highs = np.maximum(opens, closes) * (1 + np.abs(np.random.normal(0, 0.003, n_days)))
lows = np.minimum(opens, closes) * (1 - np.abs(np.random.normal(0, 0.003, n_days)))
volumes = np.random.randint(100000, 1000000, n_days)

# Создание структурированного массива с финансовыми данными
ohlcv_data = np.zeros(n_days, dtype=dtype)
ohlcv_data['date'] = dates
ohlcv_data['open'] = opens
ohlcv_data['high'] = highs
ohlcv_data['low'] = lows
ohlcv_data['close'] = closes
ohlcv_data['volume'] = volumes

# Пример загрузки данных из CSV файла
stock_data = np.genfromtxt('stock_history.csv', delimiter=',', names=True, 
                           dtype=[('date', 'U10'), ('open', 'f8'), ('high', 'f8'), 
                                  ('low', 'f8'), ('close', 'f8'), ('volume', 'i4')])

# Пример импорта данных из pandas DataFrame
import pandas as pd
import yfinance as yf  # Популярная библиотека для загрузки данных с Yahoo Finance

# Загрузка исторических данных для Apple за последний год
apple_df = yf.download('AAPL', period='1y')
apple_np = apple_df.to_numpy()  # Конвертация DataFrame в numpy массив

import numpy as np

# Предположим, у нас есть временной ряд цен закрытия
close_prices = np.array([100.0, 102.5, 101.8, 103.4, 105.7, 104.9, 107.2])

# Расчет абсолютного изменения цены
price_changes = np.diff(close_prices)

# Расчет простой доходности
simple_returns = np.diff(close_prices) / close_prices[:-1]

# Расчет логарифмической доходности
log_returns = np.log(close_prices[1:] / close_prices[:-1])

# Расчет кумулятивной доходности
cumulative_simple_return = np.prod(1 + simple_returns) - 1
cumulative_log_return = np.sum(log_returns)

# Расчет среднегодовой доходности (предполагая дневные данные и 252 торговых дня в году)
annualized_return = (1 + cumulative_simple_return) ** (252 / len(simple_returns)) - 1

print(f"Дневные изменения цены: {price_changes}")
print(f"Простые доходности: {simple_returns}")
print(f"Логарифмические доходности: {log_returns}")
print(f"Кумулятивная доходность: {cumulative_simple_return:.4f} (простая), {cumulative_log_return:.4f} (лог)")
print(f"Годовая доходность: {annualized_return:.4f}")

# Расчет волатильности на основе логарифмических доходностей
daily_volatility = np.std(log_returns)
annualized_volatility = daily_volatility * np.sqrt(252)  # Масштабирование для годовой волатильности

# Расчет коэффициента Шарпа (предполагая безрисковую ставку 0%)
mean_daily_return = np.mean(log_returns)
annualized_return = mean_daily_return * 252
sharpe_ratio = annualized_return / annualized_volatility

# Расчет скользящей волатильности с окном 5 дней
def rolling_volatility(returns, window=5):
    vol = np.zeros(len(returns) - window + 1)
    for i in range(len(vol)):
        vol[i] = np.std(returns[i:i+window])
    return vol

rolling_vol = rolling_volatility(log_returns)

# Расчет максимальной просадки
def max_drawdown(prices):
    # Создаем массив накопленных максимумов
    running_max = np.maximum.accumulate(prices)
    # Вычисляем просадку как (текущая цена / предыдущий максимум - 1)
    drawdown = prices / running_max - 1
    # Возвращаем максимальную просадку
    return np.min(drawdown)

max_dd = max_drawdown(close_prices)

print(f"Дневная волатильность: {daily_volatility:.4f}")
print(f"Годовая волатильность: {annualized_volatility:.4f}")
print(f"Коэффициент Шарпа: {sharpe_ratio:.4f}")
print(f"Максимальная просадка: {max_dd:.4f}")

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import pandas as pd
from datetime import datetime, timedelta

# Генерация тестовых данных
initial_price = 100.0
daily_returns = np.random.normal(0.0005, 0.01, 100)  # 100 дней
close_prices = initial_price * np.cumprod(1 + daily_returns)

# Даты
end_date = datetime.now().date()
dates = [end_date - timedelta(days=i) for i in range(len(close_prices))][::-1]

# Реализация индикаторов
def moving_average(data, window):
    weights = np.ones(window) / window
    ma_valid = np.convolve(data, weights, mode='valid')
    padded_ma = np.concatenate([np.full(window - 1, np.nan), ma_valid])
    return padded_ma

def exponential_moving_average(data, span):
    alpha = 2 / (span + 1)
    ema = np.zeros_like(data)
    ema[0] = data[0]
    for i in range(1, len(data)):
        ema[i] = alpha * data[i] + (1 - alpha) * ema[i - 1]
    return ema

def macd(data, fast=12, slow=26, signal=9):
    fast_ema = exponential_moving_average(data, fast)
    slow_ema = exponential_moving_average(data, slow)
    macd_line = fast_ema - slow_ema
    signal_line = exponential_moving_average(macd_line, signal)
    histogram = macd_line - signal_line
    return macd_line, signal_line, histogram

def relative_strength_index(returns, period=14):
    gains = np.where(returns > 0, returns, 0)
    losses = np.where(returns < 0, -returns, 0)

    avg_gain = np.zeros_like(returns)
    avg_loss = np.zeros_like(returns)

    avg_gain[:period] = np.mean(gains[:period])
    avg_loss[:period] = np.mean(losses[:period])

    for i in range(period, len(returns)):
        avg_gain[i] = (avg_gain[i - 1] * (period - 1) + gains[i]) / period
        avg_loss[i] = (avg_loss[i - 1] * (period - 1) + losses[i]) / period

    rs = avg_gain / np.where(avg_loss != 0, avg_loss, 1e-9)
    rsi = 100 - (100 / (1 + rs))
    return rsi

# Расчет индикаторов
ma_20 = moving_average(close_prices, 20)
ema_12 = exponential_moving_average(close_prices, 12)
macd_line, signal_line, histogram = macd(close_prices)
rsi = relative_strength_index(np.diff(close_prices) / close_prices[:-1])

# Добавляем NaN для соответствия длине графика
ma_20_padded = np.concatenate([np.full(len(close_prices) - len(ma_20), np.nan), ma_20])
ema_12_padded = ema_12  # Уже совпадает по длине

# === Визуализация ===
fig = plt.figure(figsize=(14, 10))
gs = GridSpec(3, 1, figure=fig, height_ratios=[3, 1, 1])

# 1. Основной график: цена, MA, EMA
ax1 = fig.add_subplot(gs[0])
ax1.plot(dates, close_prices, label='Цена закрытия', color='black')
ax1.plot(dates, ma_20_padded, label='MA(20)', linestyle='--', color='blue')
ax1.plot(dates, ema_12_padded, label='EMA(12)', linestyle='-.', color='orange')
ax1.set_title('Цена закрытия + MA(20) + EMA(12)')
ax1.grid(True)
ax1.legend()

# 2. MACD
ax2 = fig.add_subplot(gs[1], sharex=ax1)
ax2.bar(dates[-len(histogram):], histogram, label='MACD Гистограмма', color='gray', alpha=0.5)
ax2.plot(dates[-len(macd_line):], macd_line, label='MACD Line', color='blue')
ax2.plot(dates[-len(signal_line):], signal_line, label='Signal Line', color='red', linestyle='--')
ax2.axhline(0, linestyle=':', color='black', alpha=0.5)
ax2.set_title('MACD')
ax2.grid(True)
ax2.legend()

# 3. RSI
ax3 = fig.add_subplot(gs[2], sharex=ax1)
ax3.plot(dates[-len(rsi):], rsi, label='RSI', color='purple')
ax3.axhline(30, linestyle='--', color='red', alpha=0.3)
ax3.axhline(70, linestyle='--', color='green', alpha=0.3)
ax3.set_title('RSI (14)')
ax3.set_ylim(0, 100)
ax3.grid(True)
ax3.legend()

plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Создаем временной ряд с пропущенными значениями (np.nan)
prices_with_gaps = np.array([100.0, 101.2, np.nan, 103.5, np.nan, np.nan, 105.2, 106.1])

# Определение маски для обнаружения пропущенных значений
missing_mask = np.isnan(prices_with_gaps)
valid_mask = ~missing_mask

# 1. Заполнение предыдущими значениями (forward fill)
def forward_fill(arr):
    # Создаем маску пропущенных значений
    mask = np.isnan(arr)
    
    # Индексы непропущенных значений
    idx = np.where(~mask)[0]
    
    # Для каждого пропущенного значения находим ближайший предыдущий непропущенный индекс
    # np.searchsorted с side='right' даст нам индекс первого элемента, который больше искомого
    # Затем вычитаем 1, чтобы получить последний элемент, который меньше или равен искомому
    fill_idx = np.searchsorted(idx, np.arange(len(arr)), side='right') - 1
    
    # Учитываем случай, когда первые элементы могут быть пропущены
    valid_fill_idx = fill_idx >= 0
    out = np.copy(arr)
    
    # Заполняем только пропущенные значения
    out[mask & valid_fill_idx] = arr[idx[fill_idx[mask & valid_fill_idx]]]
    
    return out

# 2. Линейная интерполяция
def linear_interpolate(arr):
    # Создаем копию исходного массива
    result = np.copy(arr)
    
    # Находим индексы непропущенных значений
    valid_indices = np.where(~np.isnan(arr))[0]
    
    # Для каждого участка между непропущенными значениями выполняем интерполяцию
    for i in range(len(valid_indices) - 1):
        start_idx = valid_indices[i]
        end_idx = valid_indices[i + 1]
        
        if end_idx - start_idx > 1:  # Есть пропуски между двумя валидными значениями
            # Рассчитываем шаг изменения
            step = (arr[end_idx] - arr[start_idx]) / (end_idx - start_idx)
            
            # Заполняем пропуски
            for j in range(start_idx + 1, end_idx):
                result[j] = arr[start_idx] + step * (j - start_idx)
    
    return result

# Применяем различные методы заполнения пропусков
filled_forward = forward_fill(prices_with_gaps)
filled_interpolate = linear_interpolate(prices_with_gaps)

# Для сравнения также используем встроенные методы numpy
# Заполнение константой (например, средним значением)
mean_value = np.nanmean(prices_with_gaps)
filled_mean = np.where(missing_mask, mean_value, prices_with_gaps)

print("Исходный ряд с пропусками:", prices_with_gaps)
print("После forward fill:", filled_forward)
print("После линейной интерполяции:", filled_interpolate)
print("После заполнения средним:", filled_mean)

# Визуализация результатов
plt.figure(figsize=(12, 6))
time_points = np.arange(len(prices_with_gaps))

plt.plot(time_points[valid_mask], prices_with_gaps[valid_mask], 'o-', label='Исходные данные')
plt.plot(time_points, filled_forward, '--', label='Forward Fill')
plt.plot(time_points, filled_interpolate, '-.', label='Линейная интерполяция')
plt.plot(time_points, filled_mean, ':', label='Заполнение средним')

plt.xlabel('Время')
plt.ylabel('Цена')
plt.title('Сравнение методов заполнения пропусков в финансовом временном ряду')
plt.legend()
plt.grid(True)

# Генерируем временной ряд с выбросами
np.random.seed(42)
normal_returns = np.random.normal(0.0005, 0.01, 100)  # Нормальные доходности
# Добавляем несколько выбросов
normal_returns[25] = 0.15  # Сильный положительный выброс
normal_returns[50] = -0.12  # Сильный отрицательный выброс
normal_returns[75] = 0.08  # Умеренный положительный выброс

# Обнаружение выбросов с использованием z-score
def detect_outliers_zscore(data, threshold=3):
    mean = np.mean(data)
    std = np.std(data)
    z_scores = np.abs((data - mean) / std)
    return z_scores > threshold

# Обнаружение выбросов с использованием метода межквартильного размаха (IQR)
def detect_outliers_iqr(data, factor=1.5):
    q25 = np.percentile(data, 25)
    q75 = np.percentile(data, 75)
    iqr = q75 - q25
    lower_bound = q25 - factor * iqr
    upper_bound = q75 + factor * iqr
    return (data < lower_bound) | (data > upper_bound)

# Применяем методы обнаружения выбросов
outliers_zscore = detect_outliers_zscore(normal_returns)
outliers_iqr = detect_outliers_iqr(normal_returns)

# Обработка выбросов путем винсоризации (обрезания до границ)
def winsorize(data, limits):
    if isinstance(limits, tuple):
        lower_limit, upper_limit = limits
    else:
        lower_limit = upper_limit = limits
    
    lower_percentile = np.percentile(data, lower_limit * 100)
    upper_percentile = np.percentile(data, 100 - upper_limit * 100)
    
    winsorized_data = np.copy(data)
    winsorized_data[data < lower_percentile] = lower_percentile winsorized_data[data > upper_percentile] = upper_percentile
    
    return winsorized_data

# Применяем винсоризацию с порогом 5%
winsorized_returns = winsorize(normal_returns, 0.05)

# Визуализация результатов
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(normal_returns, label='Исходные доходности')
plt.plot(np.where(outliers_zscore, normal_returns, np.nan), 'ro', label='Выбросы (Z-score)')
plt.title('Обнаружение выбросов с использованием Z-score')
plt.legend()
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(normal_returns, label='Исходные доходности')
plt.plot(np.where(outliers_iqr, normal_returns, np.nan), 'go', label='Выбросы (IQR)')
plt.title('Обнаружение выбросов с использованием IQR')
plt.legend()
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(normal_returns, label='Исходные доходности')
plt.plot(winsorized_returns, label='Винсоризованные доходности')
plt.title('Результат винсоризации')
plt.legend()
plt.grid(True)

plt.tight_layout()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# Генерация синтетического временного ряда
n_points = 365 * 2  # 2 года дневных данных
time = np.arange(n_points)

# Компоненты ряда
trend = 0.05 * time  # Линейный тренд
season_annual = 10 * np.sin(2 * np.pi * time / 365)  # Годовая сезонность
season_quarterly = 5 * np.sin(2 * np.pi * time / (365 / 4))  # Квартальная сезонность
noise = np.random.normal(0, 3, n_points)  # Случайный шум

# Финальный временной ряд
ts = trend + season_annual + season_quarterly + noise

# Функция для расчёта скользящего среднего
def moving_average(data, window):
    weights = np.ones(window) / window
    return np.convolve(data, weights, mode='same') 

# Выделение тренда методом скользящего среднего
window_size = 365  # Окно в один год
trend_ma = moving_average(ts, window_size)

# Детрендирование — убираем тренд из исходного ряда
detrended = ts - trend_ma

# Периодограмма для поиска доминирующих частот
frequencies, power = signal.periodogram(detrended)

# Поиск доминирующих частот
dominant_freq_idx = np.argsort(power)[-5:]  # 5 самых "мощных" частот
dominant_periods = 1 / frequencies[dominant_freq_idx]
dominant_periods = dominant_periods[dominant_periods < n_points / 2]  # Убираем слишком длинные

print("Доминирующие периоды в данных:", dominant_periods)

# Визуализация
plt.figure(figsize=(15, 10))

# Исходный ряд
plt.subplot(4, 1, 1)
plt.plot(time, ts, label='Исходный ряд')
plt.title('Исходный временной ряд')
plt.grid(True)
plt.legend()

# Тренд
plt.subplot(4, 1, 2)
plt.plot(time, trend_ma, 'r-', label='Тренд (Скользящее среднее)')
plt.plot(time, ts, 'b-', alpha=0.3, label='Оригинал')
plt.title('Выделенный тренд')
plt.grid(True)
plt.legend()

# Детрендированный ряд
plt.subplot(4, 1, 3)
plt.plot(time, detrended, label='Детрендированный ряд')
plt.title('Ряд после удаления тренда (сезонность + шум)')
plt.grid(True)
plt.legend()

# Периодограмма
plt.subplot(4, 1, 4)
plt.semilogy(frequencies, power, label='Периодограмма')
plt.title('Спектр мощности (периодограмма, логарифмическая шкала)')
plt.axvline(1 / 365, color='r', linestyle='--', label='Годовой цикл')
plt.axvline(4 / 365, color='g', linestyle='--', label='Квартальный цикл')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# Генерируем синтетические доходности для 5 активов
np.random.seed(42)
n_assets = 5
n_days = 252  # Один торговый год

# Создаем случайную матрицу корреляций
random_corr = np.random.uniform(-0.5, 0.9, size=(n_assets, n_assets))
# Обеспечиваем симметричность
random_corr = (random_corr + random_corr.T) / 2
# Устанавливаем единицы по диагонали
random_corr[np.diag_indices_from(random_corr)] = 1

# Генерируем доходности с заданной корреляционной структурой
# используя разложение Холецкого
std_devs = np.random.uniform(0.01, 0.05, n_assets)  # Годовые волатильности
means = np.random.uniform(0.03, 0.15, n_assets)  # Годовые доходности

# Преобразуем матрицу корреляций в матрицу ковариаций
cov_matrix = np.diag(std_devs) @ random_corr @ np.diag(std_devs)

# Генерируем коррелированные доходности
daily_means = means / 252
daily_covs = cov_matrix / 252
cholesky = np.linalg.cholesky(daily_covs)
uncorrelated_returns = np.random.normal(0, 1, size=(n_days, n_assets))
correlated_returns = uncorrelated_returns @ cholesky.T + daily_means

# Рассчитываем эмпирическую матрицу корреляций
empirical_corr = np.corrcoef(correlated_returns.T)

# Рассчитываем эмпирическую матрицу ковариаций
empirical_cov = np.cov(correlated_returns.T)

print("Теоретическая матрица корреляций:")
print(random_corr)
print("\nЭмпирическая матрица корреляций:")
print(empirical_corr)

# Визуализация матрицы корреляций
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.imshow(random_corr, cmap='coolwarm', vmin=-1, vmax=1)
plt.colorbar(label='Корреляция')
plt.title('Теоретическая матрица корреляций')
plt.xticks(np.arange(n_assets))
plt.yticks(np.arange(n_assets))

plt.subplot(1, 2, 2)
plt.imshow(empirical_corr, cmap='coolwarm', vmin=-1, vmax=1)
plt.colorbar(label='Корреляция')
plt.title('Эмпирическая матрица корреляций')
plt.xticks(np.arange(n_assets))
plt.yticks(np.arange(n_assets))

plt.tight_layout()

import numpy as np
import matplotlib.pyplot as plt

# Генерируем коррелированные лог-доходности для 5 активов
np.random.seed(42)
n_days = 1000
n_assets = 5

# Создаём случайную ковариационную матрицу
random_weights = np.random.random(n_assets)
random_weights /= random_weights.sum()
cov_matrix = np.cov(np.random.randn(n_assets, n_assets))
cov_matrix = cov_matrix @ cov_matrix.T  # делаем положительно определенной

# Генерируем коррелированные доходности
means = np.array([0.001, 0.0005, 0.0004, 0.0008, 0.0007])  # средние дневные доходности
correlated_returns = np.random.multivariate_normal(means, cov_matrix, size=n_days)

# Расчёт годовых доходностей и волатильностей
annual_returns = np.mean(correlated_returns, axis=0) * 252
annual_volatilities = np.std(correlated_returns, axis=0) * np.sqrt(252)
annual_cov_matrix = np.cov(correlated_returns, rowvar=False) * 252  # годовая ковариационная матрица

print("Годовые доходности активов:", annual_returns)
print("Годовые волатильности активов:", annual_volatilities)

# Определяем функции для расчётов по портфелю
def portfolio_return(weights, returns):
    return np.sum(weights * returns)

def portfolio_volatility(weights, cov_matrix):
    return np.sqrt(weights.T @ cov_matrix @ weights)

def sharpe_ratio(weights, returns, cov_matrix):
    ret = portfolio_return(weights, returns)
    vol = portfolio_volatility(weights, cov_matrix)
    return ret / vol if vol != 0 else 0

# Генерируем случайные портфели
n_portfolios = 10000
all_weights = np.random.random((n_portfolios, n_assets))
all_weights = all_weights / np.sum(all_weights, axis=1)[:, np.newaxis]

all_returns = np.array([portfolio_return(w, annual_returns) for w in all_weights])
all_volatilities = np.array([portfolio_volatility(w, annual_cov_matrix) for w in all_weights])
all_sharpe_ratios = all_returns / all_volatilities

# Находим лучшие портфели
max_sharpe_idx = np.argmax(all_sharpe_ratios)
max_sharpe_weights = all_weights[max_sharpe_idx]
max_sharpe_return = all_returns[max_sharpe_idx]
max_sharpe_volatility = all_volatilities[max_sharpe_idx]

min_vol_idx = np.argmin(all_volatilities)
min_vol_weights = all_weights[min_vol_idx]
min_vol_return = all_returns[min_vol_idx]
min_vol_volatility = all_volatilities[min_vol_idx]

# Вывод результатов
print("\nПортфель с максимальным коэффициентом Шарпа:")
print(f"Доходность: {max_sharpe_return:.4f}")
print(f"Волатильность: {max_sharpe_volatility:.4f}")
print(f"Коэффициент Шарпа: {all_sharpe_ratios[max_sharpe_idx]:.4f}")
print("Веса активов:", max_sharpe_weights)

print("\nПортфель с минимальной волатильностью:")
print(f"Доходность: {min_vol_return:.4f}")
print(f"Волатильность: {min_vol_volatility:.4f}")
print(f"Коэффициент Шарпа: {all_sharpe_ratios[min_vol_idx]:.4f}")
print("Веса активов:", min_vol_weights)

# Визуализация эффективной границы
plt.figure(figsize=(12, 8))
plt.scatter(all_volatilities, all_returns, c=all_sharpe_ratios, cmap='viridis', alpha=0.5)
plt.colorbar(label='Коэффициент Шарпа')

# Отдельные активы
plt.scatter(annual_volatilities, annual_returns, c='red', marker='*', s=200, label='Отдельные активы')

# Портфель с максимальным Шарпом
plt.scatter(max_sharpe_volatility, max_sharpe_return, c='gold', marker='X', s=200, label='Максимальный Шарп')

# Портфель с минимумом волатильности
plt.scatter(min_vol_volatility, min_vol_return, c='green', marker='P', s=200, label='Минимум волатильности')

plt.title('Эффективная граница портфелей')
plt.xlabel('Волатильность (годовая)')
plt.ylabel('Доходность (годовая)')
plt.legend()
plt.grid(True)
plt.show()

# Пример использования метода Монте-Карло для моделирования цен акций
# и оценки стоимости опциона колл

# Параметры модели
S0 = 100  # Начальная цена акции
K = 105   # Цена исполнения опциона
T = 1.0   # Срок до экспирации (в годах)
r = 0.05  # Безрисковая ставка
sigma = 0.2  # Волатильность
n_simulations = 10000  # Количество симуляций
n_steps = 252  # Количество шагов (торговых дней) в году

# Геометрическое броуновское движение для моделирования цен акций
def simulate_gbm(S0, mu, sigma, T, n_steps, n_simulations):
    dt = T / n_steps
    nudt = (mu - 0.5 * sigma**2) * dt
    sigdt = sigma * np.sqrt(dt)
    
    # Генерируем случайные приращения для всех симуляций и шагов сразу
    Z = np.random.normal(0, 1, size=(n_simulations, n_steps))
    
    # Вычисляем логарифмические доходности
    log_returns = nudt + sigdt * Z
    
    # Начинаем с логарифма начальной цены
    log_price_paths = np.zeros((n_simulations, n_steps + 1))
    log_price_paths[:, 0] = np.log(S0)
    
    # Накапливаем доходности для получения логарифмических цен
    log_price_paths[:, 1:] = log_price_paths[:, 0, np.newaxis] + np.cumsum(log_returns, axis=1)
    
    # Преобразуем обратно в цены
    price_paths = np.exp(log_price_paths)
    
    return price_paths

# Моделируем пути цен акций
price_paths = simulate_gbm(S0, r, sigma, T, n_steps, n_simulations)

# Рассчитываем стоимость опциона колл по формуле: max(S_T - K, 0) * exp(-r*T)
option_payoffs = np.maximum(price_paths[:, -1] - K, 0)
option_price_mc = np.mean(option_payoffs) * np.exp(-r * T)

# Рассчитываем теоретическую цену по формуле Блэка-Шоулза для сравнения
from scipy.stats import norm

def black_scholes_call(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

option_price_bs = black_scholes_call(S0, K, T, r, sigma)

print(f"Цена опциона колл (Монте-Карло): {option_price_mc:.4f}")
print(f"Цена опциона колл (Блэк-Шоулз): {option_price_bs:.4f}")

# Визуализация нескольких путей и распределения конечных цен
plt.figure(figsize=(15, 10))

plt.subplot(2, 1, 1)
# Отображаем только первые 100 путей для наглядности
plt.plot(np.arange(n_steps + 1) / n_steps, price_paths[:100].T, alpha=0.1, color='blue')
plt.axhline(K, color='r', linestyle='--', label=f'Цена исполнения (K={K})')
plt.title(f'Моделирование путей цены акции (показаны первые 100 из {n_simulations})')
plt.xlabel('Время (лет)')
plt.ylabel('Цена акции')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.hist(price_paths[:, -1], bins=50, density=True, alpha=0.7)
plt.axvline(S0 * np.exp(r * T), color='g', linestyle='-', 
            label=f'Ожидаемая цена ({S0 * np.exp(r * T):.2f})')
plt.axvline(K, color='r', linestyle='--', label=f'Цена исполнения (K={K})')
plt.title('Распределение цен акции в момент экспирации')
plt.xlabel('Цена акции')
plt.ylabel('Плотность вероятности')
plt.legend()
plt.grid(True)

plt.tight_layout()

import numpy as np
import matplotlib.pyplot as plt

# Генерация синтетических данных для демонстрации
np.random.seed(42)
n_samples = 100
x = np.sort(np.random.uniform(0, 1, n_samples))
y_true = 3 * x + 2
y = y_true + np.random.normal(0, 0.1, n_samples)  # Добавляем шум

# Формирование признаков для линейной регрессии
X = np.vstack([np.ones(n_samples), x]).T  # Добавляем столбец единиц для свободного члена

# Аналитическое решение методом наименьших квадратов
# w = (X^T X)^(-1) X^T y
w = np.linalg.inv(X.T @ X) @ X.T @ y

# Полученные коэффициенты регрессии
intercept, slope = w
print(f"Истинные коэффициенты: [2, 3]")
print(f"Оцененные коэффициенты: [{intercept:.4f}, {slope:.4f}]")

# Предсказание и расчет коэффициента детерминации R^2
y_pred = X @ w
r2 = 1 - np.sum((y - y_pred)**2) / np.sum((y - y.mean())**2)
print(f"Коэффициент детерминации R^2: {r2:.4f}")

# Визуализация результатов
plt.figure(figsize=(10, 6))
plt.scatter(x, y, alpha=0.7, label='Наблюдения')
plt.plot(x, y_true, 'r--', label='Истинная зависимость')
plt.plot(x, y_pred, 'g-', label='Предсказание модели')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Линейная регрессия с использованием NumPy')
plt.legend()
plt.grid(True)
plt.show()

def create_time_windows(data, window_size, horizon=1):
    """
    Создает матрицу входных данных X и целевых значений y для прогнозирования временных рядов.
    
    Args:
        data: одномерный массив временного ряда
        window_size: размер окна (количество предыдущих наблюдений для прогноза)
        horizon: горизонт прогнозирования (количество шагов вперёд)
    
    Returns:
        X: матрица признаков, где каждая строка содержит window_size последовательных наблюдений
        y: массив целевых значений для прогнозирования
    """
    n_samples = len(data) - window_size - horizon + 1
    X = np.zeros((n_samples, window_size))
    y = np.zeros(n_samples)
    
    for i in range(n_samples):
        X[i, :] = data[i:i+window_size]
        y[i] = data[i+window_size+horizon-1]
    
    return X, y

# Пример использования
# Генерируем синтетический временной ряд (простая синусоида с шумом)
n_points = 200
time = np.arange(n_points)
signal = np.sin(0.05 * time) + 0.2 * np.random.normal(0, 1, n_points)

# Создаем окна для прогнозирования на 1 шаг вперед с использованием 10 предыдущих значений
window_size = 10
horizon = 1
X, y = create_time_windows(signal, window_size, horizon)

print(f"Размер матрицы признаков: {X.shape}")
print(f"Размер вектора целей: {y.shape}")

plt.figure(figsize=(12, 6))
plt.plot(time, signal, label='Исходный сигнал')
plt.plot(time[window_size:window_size+len(y)], y, 'ro', markersize=3, alpha=0.6, label='Целевые значения')
plt.title('Временной ряд и целевые значения для прогнозирования')
plt.legend()
plt.grid(True)
plt.show()

class SimpleNeuralNetwork:
    def __init__(self, input_size, hidden_size, output_size, learning_rate=0.01):
        # Инициализация весов
        self.W1 = np.random.randn(input_size, hidden_size) * 0.01
        self.b1 = np.zeros((1, hidden_size))
        self.W2 = np.random.randn(hidden_size, output_size) * 0.01
        self.b2 = np.zeros((1, output_size))
        self.learning_rate = learning_rate
    
    def sigmoid(self, x):
        return 1 / (1 + np.exp(-x))
    
    def sigmoid_derivative(self, x):
        return x * (1 - x)
    
    def forward(self, X):
        # Прямой проход
        self.z1 = np.dot(X, self.W1) + self.b1
        self.a1 = self.sigmoid(self.z1)
        self.z2 = np.dot(self.a1, self.W2) + self.b2
        self.a2 = self.z2  # Линейная активация для регрессии
        return self.a2
    
    def backward(self, X, y, y_pred):
        # Обратное распространение ошибки
        m = X.shape[0]
        
        # Вычисление градиентов
        dz2 = (y_pred - y) / m
        dW2 = np.dot(self.a1.T, dz2)
        db2 = np.sum(dz2, axis=0, keepdims=True)
        
        dz1 = np.dot(dz2, self.W2.T) * self.sigmoid_derivative(self.a1)
        dW1 = np.dot(X.T, dz1)
        db1 = np.sum(dz1, axis=0, keepdims=True)
        
        # Обновление весов
        self.W1 -= self.learning_rate * dW1
        self.b1 -= self.learning_rate * db1
        self.W2 -= self.learning_rate * dW2
        self.b2 -= self.learning_rate * db2
    
    def train(self, X, y, epochs):
        losses = []
        for epoch in range(epochs):
            # Прямой проход
            y_pred = self.forward(X)
            
            # Вычисление функции потерь (среднеквадратичная ошибка)
            loss = np.mean((y_pred - y) ** 2)
            losses.append(loss)
            
            # Обратное распространение ошибки
            self.backward(X, y, y_pred)
            
            # Выводим прогресс каждые 100 эпох
            if epoch % 100 == 0:
                print(f"Эпоха {epoch}: потери = {loss:.6f}")
        
        return losses

# Пример использования для прогнозирования финансового временного ряда
# Генерируем синтетические данные
np.random.seed(42)
n_points = 200
X_raw = np.linspace(0, 4 * np.pi, n_points)
y_raw = np.sin(X_raw) + 0.1 * np.random.randn(n_points)

# Нормализация данных
X = (X_raw - np.mean(X_raw)) / np.std(X_raw)
y = (y_raw - np.mean(y_raw)) / np.std(y_raw)

# Создаем окна для прогнозирования
window_size = 5
X_windows, y_targets = create_time_windows(y, window_size)

# Разделение на обучающую и тестовую выборки
train_size = int(0.8 * len(X_windows))
X_train, X_test = X_windows[:train_size], X_windows[train_size:]
y_train, y_test = y_targets[:train_size], y_targets[train_size:]

# Создаем и обучаем нейронную сеть
nn = SimpleNeuralNetwork(input_size=window_size, hidden_size=10, output_size=1, learning_rate=0.03)
losses = nn.train(X_train, y_train.reshape(-1, 1), epochs=500)

# Оценка на тестовых данных
y_pred = nn.forward(X_test)
mse = np.mean((y_pred - y_test.reshape(-1, 1)) ** 2)
print(f"Среднеквадратичная ошибка на тестовых данных: {mse:.6f}")

# Визуализация результатов
plt.figure(figsize=(15, 10))

# График обучения
plt.subplot(2, 1, 1)
plt.plot(losses)
plt.title('Динамика функции потерь')
plt.xlabel('Эпоха')
plt.ylabel('MSE')
plt.grid(True)

# Сравнение предсказаний с реальными данными
plt.subplot(2, 1, 2)
plt.plot(y_test, label='Реальные значения')
plt.plot(y_pred, label='Предсказания')
plt.title('Сравнение предсказаний с реальными данными')
plt.xlabel('Индекс выборки')
plt.ylabel('Нормализованное значение')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()