ポートフォリオ最適化¶

import datetime
import pandas as pd

def load_stock_prices() -> pd.DataFrame:
    return pd.read_csv(
        "../../../storage/portfolio/dummy_stock_price.csv",
        index_col="Date",
        parse_dates=True,
    )

import amplify
import numpy as np

def calculate_return_rates(prices: np.ndarray, num_days_operation: int) -> np.ndarray:
    return (prices[num_days_operation:] - prices[:-num_days_operation]) / prices[
        :-num_days_operation
    ]

def optimize_portfolio(
    historical_data: pd.DataFrame,
    num_days_operation: int,
    gamma: float = 20,
    max_w: int = 20,
    constraint_weight: float = 1.0,
    timeout: datetime.timedelta = datetime.timedelta(seconds=5),
):
    # 銘柄名のリスト
    stock_names = list(historical_data.columns)
    
    # 投資比率(%) を表す変数 `w_i` を作成 (0 以上 `max_w` 以下の値を取る整数変数)
    gen = amplify.VariableGenerator()
    w = gen.array("Integer", len(stock_names), bounds=(0, max_w))
    
    # 制約条件を作成 (w の総和は 100)
    constraint = amplify.equal_to(w.sum(), 100)
    
    # 目的関数を作成 (収益率の平均を最大化)
    w_ratio = w / 100  # w (単位は %) を実数に変換したもの
    
    # num_days_operation 営業日運用した場合の銘柄ごとの収益率を計算する
    return_rates = calculate_return_rates(
        historical_data.to_numpy(), num_days_operation
    )
    
    # ポートフォリオの収益率を定式化
    # それぞれの銘柄の (収益率の平均) * (投資割合) を足し合わせる
    portfolio_return_rate = (w_ratio * np.mean(return_rates, axis=0)).sum()
    
    # ポートフォリオの共分散 (二次元配列) を計算
    # 配列の i 行 j 列は銘柄 i と銘柄 j の収益率の共分散を表す
    covariance_matrix = np.cov(return_rates, rowvar=False)
    
    # ポートフォリオのリスクを定式化
    # 全体の収益率の分散を表す w についての二次多項式
    portfolio_variance = w_ratio @ covariance_matrix @ w_ratio  # type: ignore
    
    # 目的関数の定式化 (gamma はリスク回避度を表すパラメータ)
    objective = -portfolio_return_rate + 0.5 * gamma * portfolio_variance
    
    # 最適化モデルを作成
    model = amplify.Model(objective, constraint_weight * constraint)
    
    # ソルバークライアントの作成とソルバーの設定
    client = amplify.FixstarsClient()
    client.parameters.timeout = timeout
    # ローカル環境等で使用する場合は Fixstars Amplify AE のアクセストークンを入力してください
    # client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
    
    # 最適化を実行
    result = amplify.solve(model, client)
    
    # 実行結果を解析
    if len(result) == 0:
        raise RuntimeError("No feasible solution found")
    
    # 全ての銘柄それぞれにいくつ (%) 投資するかを得る
    w_values = w.evaluate(result.best.values)
    
    # 投資比率が 0 より大きい銘柄のみを抽出してポートフォリオを作成
    portfolio = {
        stock_name: int(w_value)
        for stock_name, w_value in zip(stock_names, w_values)
        if w_value > 0
    }
    
    # 得られたポートフォリオの収益率とリスクを計算
    return_rate = portfolio_return_rate.evaluate(result.best.values)
    variance = portfolio_variance.evaluate(result.best.values)
    return portfolio, return_rate, variance

stock_prices = load_stock_prices()
stock_prices

import matplotlib.pyplot as plt
    
plt.plot(stock_prices["salmon"], color="salmon")
plt.plot(stock_prices["darkslategray"], color="darkslategray")
plt.plot(stock_prices["hotpink"], color="hotpink")
plt.show()

stock_prices_history = stock_prices.loc["2023":"2023"]

portfolio, return_rate, variance = optimize_portfolio(
    stock_prices_history, num_days_operation=20
)

import matplotlib
    
# カラーマップ
colors = tuple(matplotlib.colormaps.get_cmap("Set3")(range(12)))
    
# パイチャートの描画
patches, texts, autotexts = plt.pie(  # type: ignore
    list(portfolio.values()),
    labels=list(portfolio.keys()),
    radius=1.5,
    autopct="%.f%%",
    colors=colors,
    labeldistance=0.8,
    wedgeprops={"linewidth": 1.0, "edgecolor": "white"},
    pctdistance=0.5,
)
for text in texts:
    text.set_horizontalalignment("center")
plt.show()

df_future = stock_prices.loc["2024":"2024"]
num_days_operation = 20
    
historical_return_rates = calculate_return_rates(
    stock_prices_history.to_numpy(), num_days_operation
)
max_profit_stock: str = stock_prices_history.columns[
    historical_return_rates.mean(axis=0).argmax()
]
max_profit_portfolio = {max_profit_stock: 100}
uniform_ratio_portfolio = {stock_name: 1 for stock_name in df_future.columns}
    
    
def calculate_portfolio_return_rates(portfolio: dict[str, int]):
    return_rates = calculate_return_rates(df_future.to_numpy(), num_days_operation)
    ratio_array = (
        np.array([portfolio.get(stock_name, 0) for stock_name in df_future.columns])
        / 100
    )  # 各銘柄への投資割合の配列
    return (ratio_array * return_rates).sum(
        axis=1
    )  # 各運用開始日ごとに、ポートフォリオ全体の収益率を計算する
    
    
optimized_return_rates = calculate_portfolio_return_rates(portfolio)
    
max_profit_return_rates = calculate_portfolio_return_rates(max_profit_portfolio)
    
uniform_return_rates = calculate_portfolio_return_rates(uniform_ratio_portfolio)
    
print(
    f"optimized:  max return rate = {np.max(optimized_return_rates) * 100:.2f}%, "
    f"mean return rate = {np.mean(optimized_return_rates) * 100:.2f}%, "
    f"min return rate = {np.min(optimized_return_rates) * 100:.2f}%, "
    f"variance = {np.var(optimized_return_rates):.5f}"
)
print(
    f"max profit: max return rate = {np.max(max_profit_return_rates) * 100:.2f}%, "
    f"mean return rate = {np.mean(max_profit_return_rates) * 100:.2f}%, "
    f"min return rate = {np.min(max_profit_return_rates) * 100:.2f}%, "
    f"variance = {np.var(max_profit_return_rates):.5f}"
)
print(
    f"uniform:    max return rate = {np.max(uniform_return_rates) * 100:.2f}%, "
    f"mean return rate = {np.mean(uniform_return_rates) * 100:.2f}%, "
    f"min return rate = {np.min(uniform_return_rates) * 100:.2f}%, "
    f"variance = {np.var(uniform_return_rates):.5f}"
)
    
bins = np.linspace(-40, 40, 50)
plt.hist(
    optimized_return_rates * 100,
    label="optimized",
    bins=bins,  # type: ignore
    color="royalblue",
    alpha=0.8,
    zorder=3,
)
plt.hist(
    max_profit_return_rates * 100,
    label="max profit",
    bins=bins,  # type: ignore
    color="coral",
    alpha=0.8,
    zorder=1,
)
plt.hist(
    uniform_return_rates * 100,
    label="uniform ratio",
    bins=bins,  # type: ignore
    color="gold",
    alpha=0.8,
    zorder=2,
)
plt.legend()
plt.xlabel("return rate (%)")
plt.show()

TAX_RATE = 0.2
rng = np.random.default_rng()
    
    
def get_portfolio(
    prices: pd.DataFrame,
    start_date: datetime.date,
    num_days_backward: int,
    num_days_operation: int,
) -> dict[str, int]:
    """過去データを用いて最適化したポートフォリオを取得する
    
    Args:
        prices (pd.DataFrame): 株価の時系列データ
        start_date (datetime.date): 運用開始日
        num_days_backward (int): 最適化に使用する過去データの日数
    
    Returns:
        dict[str, int]: 銘柄名をキー、投資比率を値とする辞書
    """
    
    # 運用開始日の前の営業日を取得
    previous_date = start_date - datetime.timedelta(days=1)
    
    # 運用開始日からさかのぼって `num_days_backward` 日分の株価データを取得
    stock_price_history = prices.loc[: str(previous_date)].iloc[-num_days_backward:]
    
    # 運用開始前のデータを用いてポートフォリオを作成
    portfolio, _, _ = optimize_portfolio(stock_price_history, num_days_operation)
    
    return portfolio
    
    
def simulate_stock_trading(
    prices: pd.DataFrame,
    funds: float,
    start_date: datetime.date,
    num_days_operation: int,
    portfolio: dict[str, int],
    tax_rate=TAX_RATE,
) -> float:
    """与えられた運用日数とポートフォリオに基づいて株式売買をシミュレーションする
    
    Args:
        prices (pd.DataFrame): 株価の時系列データ
        funds (float): 運用資金
        start_date (datetime.date): 運用開始日
        num_days_operation (int): 運用日数
        portfolio (dict[str, int]): ポートフォリオ
        tax_rate (_type_, optional): 譲渡益税率
    
    Returns:
        float: 運用結果の資金
    """
    
    # 銘柄ごとの投資比率の配列に変換
    weights = np.array(
        [portfolio.get(stock_name, 0) / 100 for stock_name in prices.columns]
    )
    
    # 各銘柄の 1 株あたり購入額を前営業日の株価から決定
    previous_date = start_date - datetime.timedelta(days=1)
    start_prices = prices.loc[: str(previous_date)].iloc[-1].to_numpy()
    # 1% 程度の購入価格の増分を考慮
    start_prices = start_prices * rng.uniform(1.0, 1.01, size=len(prices.columns))
    
    # 各銘柄の 1 株あたり売却額を計算
    end_prices = prices.loc[str(start_date) :].iloc[num_days_operation - 1].to_numpy()
    # 1% 程度の売却金額の差分を考慮
    end_prices = end_prices * rng.uniform(0.99, 1.0, size=len(prices.columns))
    
    # 利益率を計算
    return_rate: float = (weights * (end_prices / start_prices)).sum()
    
    # 利益が出た場合は課税分を差し引く
    if return_rate > 1:
        return_rate = 1 + (1 - tax_rate) * (return_rate - 1)
    
    # 売却額 (= 購入額 x 利益率) を返却
    return funds * return_rate

def simulate_stick_operation(
    prices: pd.DataFrame,
    num_rounds: int,
    simulation_start_date: datetime.date,
    num_days_sampling: int,
    num_days_operation: int,
) -> list[tuple[datetime.date, float, datetime.date, float]]:
    """与えられた運用日数とサイクル数に基づいて株式売買をシミュレーションする
    
    Args:
        prices (pd.DataFrame): 株価の時系列データ
        num_rounds (int): ラウンド数
        simulation_start_date (datetime.date): シミュレーションの開始日
        num_days_sampling (int): 最適化に用いる過去データの日数
        num_days_operation (int): 運用日数
    
    Returns:
        list[tuple[datetime.date, float, datetime.date, float]]: _description_
    """
    
    # 開始資金
    current_funds = 1.0
    
    # (購入日, 購入額, 売却日, 売却額) を格納するためのリスト
    operation_history: list[tuple[datetime.date, float, datetime.date, float]] = []
    
    # シミュレーション開始日以降の株価の時系列データ
    prices_start = prices.loc[str(simulation_start_date) :]
    
    for i in range(num_rounds):
        # ラウンドの開始日と終了日
        start_date = prices_start.iloc[num_days_operation * i].name.date()  # type: ignore
        end_date = prices_start.iloc[
            num_days_operation * i + num_days_operation - 1
        ].name.date()  # type: ignore
    
        print(f"Round: {i+1}/{num_rounds}, {start_date} - {end_date}")
    
        # ポートフォリオの最適化
        portfolio = get_portfolio(
            prices, start_date, num_days_sampling, num_days_operation
        )
    
        # 株式売買のシミュレーション
        next_funds = simulate_stock_trading(
            prices, current_funds, start_date, num_days_operation, portfolio
        )
    
        # 運用履歴に追加
        operation_history.append((start_date, current_funds, end_date, next_funds))
    
        print(
            f"Profit: {next_funds / current_funds:.3f}, Funds: {current_funds:.3f} -> {next_funds:.3f}"
        )
    
        current_funds = next_funds
    
    return operation_history

operation_history = simulate_stick_operation(
    stock_prices, 10, datetime.date(2024, 1, 1), 100, 20
)
# 最終ラウンドの運用結果
operation_history[-1][3]

import itertools
import matplotlib.pyplot as plt
from matplotlib import dates as mdates
    
ax = plt.figure().add_subplot()
    
    
def plot(operation_history, color, label):
    for start_date, start_funds, end_date, end_funds in operation_history:
        (line,) = ax.plot(
            [start_date, end_date], [start_funds, end_funds], color=color, marker="o"
        )
    line.set_label(label)  # type: ignore
    
    for history1, history2 in itertools.pairwise(operation_history):
        _, _, end_date1, end_funds1 = history1
        start_date2, start_funds2, _, _ = history2
        ax.plot(
            [end_date1, start_date2],
            [end_funds1, start_funds2],
            color=color,
            linestyle=":",
        )
    
    
plot(operation_history, "C0", "optimized")
    
ax.legend(loc="lower right")
ax.set_xlabel("Date", fontsize=10)
ax.set_ylabel("Total asset", fontsize=10)
ax.tick_params(labelsize=10)
ax.xaxis.set_major_formatter(mdates.DateFormatter("%m/%d"))
    
plt.show()