Portfolio Optimization¶

import datetime
import pathlib
import pandas as pd
from pandas_datareader import data

ONE_DAY = datetime.timedelta(days=1)


class HistoricalDataLoader:
    def __init__(
        self,
        tickers: list[str],  # List of tickers for stocks
        operation_start_date: datetime.date,  # Start date of operation
        cache_dir: pathlib.Path,  # Directory path of stored data
        num_days_sampling=10,  # Number of days of sampling period (10 business days for default)
        num_days_operation=5,  # Number of days of operation period (5 business days for default)
        use_cache: bool = True,  # Set `False` to download all data from Stooq
        stdout: bool = False,  # Set `True` for standard output
    ):
        self._tickers = tickers
        self._num_days_sampling = num_days_sampling
        self._num_days_operation = num_days_operation
        self._operation_start_date = operation_start_date
        self._cache_dir = cache_dir
        self._stdout = stdout
        self._history = self._construct_history(use_cache)
        self._sample_end_date: None | datetime.date = None
        self._operation_end_date: None | datetime.date = None

    def _construct_history(self, use_cache: bool) -> pd.DataFrame:
        """Construct historical data of all stocks"""
        history = pd.DataFrame()
        for ticker in self._tickers:
            individual_history = self._fetch_save_individual_history(ticker, use_cache)
            history = pd.concat([history, individual_history], axis=1).reindex(
                index=individual_history.index
            )
        self._print_deletion_periods(history)
        # Delete data on non-business days
        return history.dropna(how="any")

    def _print_deletion_periods(self, history: pd.DataFrame) -> None:
        """Print non-business days"""
        mask = history.isnull().any(axis=1)
        i = 1
        while i < len(mask) - 1:
            if mask.iloc[i] and not mask.iloc[i - 1]:
                del_start_date = mask.index[i].to_pydatetime().date()  # type: ignore
                del_end_date = ""
                for j in range(i, len(mask) - 1):
                    if not mask.iloc[j + 1]:
                        del_end_date = mask.index[j].to_pydatetime().date()  # type: ignore
                        i = j
                        break
                if self._stdout:
                    print(f"deletion {del_start_date}-{del_end_date}, nan is found.")
            i += 1

    def _calc_sample_end_date(self) -> datetime.date:
        """Calculate the end date of the sampling of historical data"""
        check_date = self._operation_start_date
        # The end date of the sampling is the business day closest to the start of operation before the start of operation.
        # If you go back 14 days from the day before the start of the operation to the past, you should find at least one business day. (e.g., Chinese stocks around Chinese New Year)
        for _ in range(14):
            check_date -= ONE_DAY
            key = check_date.strftime("%Y-%m-%d")
            if key in self._history.index:
                return check_date
        raise RuntimeError(f"could not find a valid end date. {check_date=}")

    def _calc_operation_end_date(self) -> datetime.date:
        """Calculate the end date of the operation period if the historical data include the operation period"""
        hist = self.extract_operation_period()
        if len(hist) == 0:
            raise RuntimeError("historical data do not include operation period.")
        return hist.index[0].to_pydatetime().date()  # type: ignore

    def _exist_operation_start_date(self) -> bool:
        """Return whether the historical data contain data on the start date of the operation"""
        return self._operation_start_date.strftime("%Y-%m-%d") in self._history.index

    def exist_operation_duration(self) -> bool:
        """Return whether the historical data contain data for the entire operational period"""
        if not self._exist_operation_start_date():
            return False
        date_count = 0
        date = self._operation_start_date
        while date_count < self._num_days_operation:
            if date.strftime("%Y-%m-%d") in self._history.index:
                date_count += 1
            date += ONE_DAY
            if date > self._history.index[0].to_pydatetime().date():  # type: ignore
                return False
        return True

    @property
    def num_days_sampling(self) -> int:
        """Return number of days of sampling period"""
        return self._num_days_sampling

    @property
    def num_days_operation(self) -> int:
        """Return number of days of operation period"""
        return self._num_days_operation

    @property
    def tickers(self) -> list[str]:
        """Return list of tickers for stocks"""
        return self._tickers

    @property
    def history(self) -> pd.DataFrame:
        """Return historical data for the entire period"""
        return self._history

    @property
    def operation_end_date(self) -> datetime.date:
        """Return the end date of the operation period if the historical data include the operation period"""
        if self._operation_end_date is None:
            self._operation_end_date = self._calc_operation_end_date()
        return self._operation_end_date

    @property
    def sample_end_date(self) -> datetime.date:
        """Return the end date of the sampling period"""
        if self._sample_end_date is None:
            self._sample_end_date = self._calc_sample_end_date()
        return self._sample_end_date

    def _load_cache(self, filepath: pathlib.Path, num_rows: int | None = None):
        """Load historical data from saved file"""
        hist = pd.read_csv(filepath, index_col=0, nrows=num_rows)
        hist.index = pd.to_datetime(hist.index)  # type: ignore
        if self._stdout:
            print(f"Cache: {filepath.name} from {hist.index[-1]} to {hist.index[0]}")  # type: ignore
        return hist

    def _download_history(self, ticker: str, start_date: datetime.date | None = None):
        """Download historical data from Stooq since `start_date`"""
        if self._stdout:
            print(f"Download: {ticker}, from {start_date}")
        hist = data.DataReader(ticker, "stooq", start_date)
        if len(hist) == 0:
            return pd.DataFrame()
        else:
            return hist.drop(["Open", "High", "Low", "Volume"], axis=1)

    def _save_cache(self, df: pd.DataFrame, filepath: pathlib.Path):
        """Save historical data to file"""
        df.index.name = None
        df.to_csv(filepath)

    def _fetch_save_individual_history(
        self, ticker: str, use_cache: bool = True
    ) -> pd.DataFrame:
        """Obtain and store historical data for individual stock"""
        filepath = self._cache_dir.joinpath(ticker + ".csv")
        start_date: datetime.date | None = None
        hist = pd.DataFrame()
        if filepath.exists() and use_cache:
            hist = self._load_cache(filepath)
            # Start date of the required download period
            start_date = hist.index[0].date() + ONE_DAY  # type: ignore
        # +2 as a buffer to allow for a weekend
        delta_days = datetime.timedelta(days=self._num_days_operation + 2)
        if start_date is None or start_date <= self._operation_start_date + delta_days:
            hist = pd.concat([self._download_history(ticker, start_date), hist])
            self._save_cache(hist, filepath)
        return hist.rename(columns={"Close": ticker})  # Consider only the closing price

    def extract_sampling_period(self) -> pd.DataFrame:
        """Return historical data of the sampling period"""
        i_end = self._history.index.get_loc(self.sample_end_date.strftime("%Y-%m-%d"))
        ret = self._history.iloc[i_end : i_end + self._num_days_sampling]  # type: ignore
        return ret

    def extract_operation_period(self) -> pd.DataFrame:
        """Return historical data of the operational period (if included in historical data)"""
        key = self._operation_start_date.strftime("%Y-%m-%d")
        i_start = self._history.index.get_loc(key) if key in self._history.index else 0
        if not self.exist_operation_duration():
            print("The historical data do not include the operation period.")
            return pd.DataFrame()
        ret = self._history.iloc[i_start - self._num_days_operation + 1 : i_start + 1]  # type: ignore
        return ret

import amplify
import numpy as np
from typing import Any
import matplotlib.pyplot as plt


def calc_return_cov_rate(
    h: HistoricalDataLoader,  # Historical data
    purchase_prices: np.ndarray,  # Purchase prices of each stock
    is_sampling_period: bool,  # Boolean whether to use historical data of the sampling period
) -> tuple[np.ndarray | None, np.ndarray | None]:
    """Calculate the average return and covariance matrix for each stock over the operation period"""
    if is_sampling_period:
        prices_end = h.extract_sampling_period().to_numpy()
        prices_start = np.roll(prices_end, -h.num_days_operation + 1, axis=0)
        stop = h.num_days_sampling - h.num_days_operation + 1
    else:
        prices_end = h.extract_operation_period().to_numpy()
        # If the data of the operation period is not included in the stored data
        if len(prices_end) == 0:
            return None, None
        # The price of each stock at the start of the operation is the purchase price
        prices_start = purchase_prices[np.newaxis, :]
        # One return data set (i.e., no covariance matrix)
        stop = 1
    return_rates = ((prices_end - prices_start) / prices_start)[0:stop]
    if len(return_rates) == 1 or len(return_rates[0]) == 1:
        return return_rates.mean(axis=0), None
    return return_rates.mean(axis=0), np.cov(return_rates, rowvar=False)


def calc_profit_risk(
    h: HistoricalDataLoader,  # Historical data
    w_i_percent: amplify.PolyArray | np.ndarray,  # Investment ratio (%)
    purchase_prices: np.ndarray,  # Purchase prices of stocks
    is_sampling_period=True,  # Boolean whether to use historical data of the sampling period
) -> tuple[Any, Any]:
    """Calculate and return the return and variance for the entire given portfolio"""
    w_i = w_i_percent / 100
    return_rates, cov_rates = calc_return_cov_rate(
        h, purchase_prices, is_sampling_period
    )
    if return_rates is None:
        return None, None
    profit = amplify.einsum("i,i->", w_i, return_rates)
    if cov_rates is None:
        return profit, None
    risk = amplify.einsum("i,j,ij->", w_i, w_i, cov_rates)
    return profit, risk


def generate_objective(
    h: HistoricalDataLoader,  # Historical data
    w_i_percent: amplify.PolyArray,  # Investment ratio (%)
    purchase_prices: np.ndarray,  # Purchase prices of stocks
    gamma: float,  # Gamma in mean and variance model
) -> amplify.Poly:
    """Construct and return an objective function based on the mean/variance model and an additive term"""
    profit, risk = calc_profit_risk(h, w_i_percent, purchase_prices)
    return -profit + gamma * 0.5 * risk


def generate_constraint(w_i_percent: amplify.PolyArray) -> amplify.Constraint:
    """Construct and return the constraint that results in a total investment ratio of 100%"""
    return amplify.equal_to(w_i_percent, 100)


def solve(
    h: HistoricalDataLoader,  # Historical data
    purchase_prices: np.ndarray,  # Purchase prices of stocks
    gamma=1,  # Gamma in mean and variance model
    constraint_weight: float = 1,  # Weight for constraint condition
    timeout_ms: int = 5000,  # Solving timeout (ms)
    num_solves: int = 3,  # Number of serial runs for solve (https://amplify.fixstars.com/ja/docs/amplify/v1/serial.html)
    plot_optim_history=False,
) -> np.ndarray:
    """Optimize portfolio based on historical data, estimated purchase price of each stock, and investment budget"""
    gen = amplify.VariableGenerator()

    # q: Investment ratio (%)
    # For diversification, limit ratio of each stock to 50%
    q = gen.array("Integer", len(h.tickers), bounds=(0, 50))

    objective = generate_objective(h, q, purchase_prices, gamma)
    constraint = generate_constraint(q)
    model = amplify.Model(objective, constraint_weight * constraint)

    client = amplify.FixstarsClient()
    client.parameters.timeout = datetime.timedelta(milliseconds=timeout_ms)
    client.parameters.outputs.num_outputs = 0  # Return all solutions found
    # client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # If you use Amplify in a local environment or Google Colaboratory, enter your Amplify API token.

    # Solve the problem
    result = amplify.solve(model, client, num_solves=num_solves)

    if len(result) == 0:
        print("no solution is found")

    if plot_optim_history:
        optim_history(result)

    return q.evaluate(result.best.values)


def optim_history(result: amplify.Result) -> None:
    """Display optimization history"""
    ax = plt.figure().add_subplot()
    # Obtain the time of each solution and the value of the objective function
    for i in range(result.num_solves):
        times = [solution.time.total_seconds() for solution in result.split[i]]
        objectives = [solution.objective for solution in result.split[i]]
        ax.plot(times, objectives, marker="*")
    ax.set_xlabel("elapsed time in seconds")
    ax.set_ylabel("objective value")
    ax.grid(True)
    plt.show()

# Directory path for storing historical data
history_data_dir = pathlib.Path("../../../storage/portfolio")

# Start date of operation
date_of_operation = datetime.date.fromisoformat("2024-04-01")

# NASDAQ100 component tickers (May 18, 2024)
nasdaq100_stocks = ["ADBE", "ADP", "ABNB", "GOOGL", "GOOG", "AMZN", "AMD", "AEP", "AMGN", "ADI",
                     "ANSS", "AAPL", "AMAT", "ASML", "AZN", "TEAM", "ADSK", "BKR", "BIIB", "BKNG",
                     "AVGO", "CDNS", "CDW", "CHTR", "CTAS", "CSCO", "CCEP", "CTSH", "CMCSA", "CEG",
                     "CPRT", "CSGP", "COST", "CRWD", "CSX", "DDOG", "DXCM", "FANG", "DLTR", "DASH",
                     "EA", "EXC", "FAST", "FTNT", "GEHC", "GILD", "GFS", "HON", "IDXX", "ILMN",
                     "INTC", "INTU", "ISRG", "KDP", "KLAC", "KHC", "LRCX", "LIN", "LULU", "MAR",
                     "MRVL", "MELI", "META", "MCHP", "MU", "MSFT", "MRNA", "MDLZ", "MDB", "MNST",
                     "NFLX", "NVDA", "NXPI", "ORLY", "ODFL", "ON", "PCAR", "PANW", "PAYX", "PYPL",
                     "PDD", "PEP", "QCOM", "REGN", "ROP", "ROST", "SIRI", "SBUX", "SNPS", "TTWO",
                     "TMUS", "TSLA", "TXN", "TTD", "VRSK", "VRTX", "WBA", "WBD", "WDAY", "XEL", "ZS"]  # fmt: skip

# Instantiate historical data class
nasdaq100_stocks_history = HistoricalDataLoader(
    nasdaq100_stocks,
    date_of_operation,
    history_data_dir,
    stdout=True,
)

nasdaq100_stocks_history.extract_sampling_period()["ADBE"]

rng = np.random.default_rng(0)


def calc_purchase_prices(h: HistoricalDataLoader, low=1.0, high=1.01) -> np.ndarray:
    """Calculate the purchase prices of the subject stocks. Add random prices in the range of [low, high] to each closing price on the day before the operation start date."""
    # Closing prices on the day before the operation start date (i.e., the last day of the sampling period)
    close = h.extract_sampling_period()[h.sample_end_date : h.sample_end_date].iloc[0]
    return close.to_numpy() * rng.uniform(low, high, size=len(close))


# The purchase prices at start of operation
purchase_prices_stocks = calc_purchase_prices(nasdaq100_stocks_history)

solution = solve(nasdaq100_stocks_history, purchase_prices_stocks)

import math

TAX_RATE = 0.2


def evaluate(
    sol: np.ndarray,
    purchase_prices: np.ndarray,
    h: HistoricalDataLoader,
    show_historical=False,
    tax_rate=TAX_RATE,
) -> float | None:
    """Evaluate the investment results based on the portfolio and return the overall portfolio rate of return based on historical data and (if data exist) the rate of return over the investment period."""

    portfolio = {
        ticker: f"{int(sol[i])}%" for i, ticker in enumerate(h.tickers) if sol[i] > 0
    }
    print(f"portfolio: {portfolio}")

    # Rate of return/risk for the sampling period
    h_profit, h_risk = calc_profit_risk(h, sol, purchase_prices)
    if show_historical:
        print(f"historical profit: {float(h_profit)*100:.1f} %")
        if h_risk is not None:
            # Display standard deviation value as risk
            print(f"historical risk: {math.sqrt(float(h_risk))*100:.1f} %")

    op_profit, _ = calc_profit_risk(h, sol, purchase_prices, is_sampling_period=False)

    if op_profit is None:
        return None

    # Rate of return after operation (Only if data during the investment period is given, `tax_rate` on earnings is also considered)
    print(f"operational profit: {float(op_profit)*100:.1f} %")
    return (
        float(op_profit) * (1 - tax_rate) if float(op_profit) > 0 else float(op_profit)
    )


# Evaluate the optimal portfolio found in the upper cell
return_rate = evaluate(
    solution,
    purchase_prices_stocks,
    nasdaq100_stocks_history,
    show_historical=True,
)

nasdaq100_index_history = HistoricalDataLoader(
    ["^NDX"],  # Ticker of NASDAQ 100 Index
    date_of_operation,
    history_data_dir,
)

# The purchase price at start of operation
purchase_prices_index = calc_purchase_prices(nasdaq100_index_history)

# The investment ratio is set to 100% since only one stock (index) is used
return_rate = evaluate(
    np.array([100]),
    purchase_prices_index,
    nasdaq100_index_history,
    show_historical=True,
)

def simulate_operation(
    tickers: list[str],  # Tickers of stocks
    operation_date: datetime.date | list[datetime.date],  # Start date of operation
    num_operation_rounds: int,  # Number of operation rounds
) -> tuple[list[datetime.date], list[float]]:
    """Operation simulation of buying and selling over `num_operation_rounds` time's investment rounds."""

    # Scan historical data and store the start and end dates of each operation round in `operation_dates`
    operation_dates: list[datetime.date] = []
    if isinstance(operation_date, list):
        operation_dates = operation_date
    else:
        i_operation_round = 0
        while i_operation_round < num_operation_rounds:
            while True:
                h = HistoricalDataLoader(tickers, operation_date, history_data_dir)
                #  If `operation_date` is not a business day, revalidate one day later as a business day
                if not h.exist_operation_duration():
                    operation_date += ONE_DAY
                    break
                i_operation_round += 1
                operation_dates.append(operation_date)
                operation_dates.append(h.operation_end_date)
                operation_date += datetime.timedelta(days=h.num_days_sampling)
                break

    # List to include the amount of assets at the beginning and end of each investment round
    fund_evolution: list[float] = []
    funds = 1.0  # Initial asset value
    print(f"Initial funds: {funds:.2f}")

    for i, date in enumerate(operation_dates):
        # If `i` is odd, `date` is the end date of the operation
        if i % 2 == 1:
            continue

        h = HistoricalDataLoader(tickers, date, history_data_dir)

        print("============================================")
        print(f"Operation round: {int(i / 2) + 1} of {num_operation_rounds}, {date}")

        # Add the amount of assets at the start of the current round to `fund_evolution`
        fund_evolution.append(funds)

        # Purchase prices of all stocks at the start of current round
        prices = calc_purchase_prices(h)

        # If there is more than one candidate of stocks, find the optimal portfolio; if one candidate, buy 100% of that stock
        if len(prices) > 1:
            solution = solve(h, prices)
        else:
            solution = np.array([100])

        # Evaluate the investment results in current round and calculate the amount of assets after increase/decrease
        operational_return = evaluate(solution, prices, h)
        if operational_return is not None:
            funds *= 1.0 + operational_return
            print(f"Current funds: {funds:.2f}")

        # Add the amount of assets at the end of the current round to `fund_evolution`
        fund_evolution.append(funds)

    return operation_dates, fund_evolution

start_date = datetime.date.fromisoformat(
    "2024-01-01"
)  # Start date of operation of the first round
num_operation_rounds = 5  # Number of operation rounds

operation_dates, optim_portfo_evolution = simulate_operation(
    nasdaq100_stocks, start_date, num_operation_rounds
)

_, index_op_evolution = simulate_operation(["^NDX"], start_date, num_operation_rounds)

round_start_date = operation_dates[0]
round_end_date = operation_dates[-1]

h = HistoricalDataLoader(["^NDX"], round_start_date, history_data_dir)

# Purchase price on the start date of the first investment round
prices = calc_purchase_prices(h)

# Index on the last day of the last operational round
i_end = h.history.index.get_loc(round_end_date.strftime("%Y-%m-%d"))

# Taxed funds after sale on the last day of the last investment round
return_rate = h.history.iloc[i_end].iloc[0] / prices[0]
if return_rate > 1:
    return_rate = (return_rate - 1) * (1 - TAX_RATE) + 1
index_no_op_evolution = [1, return_rate]

import matplotlib.pyplot as plt
from matplotlib import dates as mdates

ax = plt.figure().add_subplot()
color_stocks = [0.2, 0.2, 1]
color_index = [1, 0.2, 0.2]

# Plot "optim portfo", "index-op" for all operational rounds
operation_dates = np.array(operation_dates)
ax.plot(operation_dates, optim_portfo_evolution, linestyle=":", color=color_stocks)
ax.plot(operation_dates, index_op_evolution, linestyle=":", color=color_index)
for i, date in enumerate(operation_dates):
    text = "E"
    if i % 2 == 0:
        text = "S"
        (optimal,) = ax.plot(
            operation_dates[i : i + 2],
            optim_portfo_evolution[i : i + 2],
            color=color_stocks,
            marker="*",
        )
        (index,) = ax.plot(
            operation_dates[i : i + 2],
            index_op_evolution[i : i + 2],
            color=color_index,
            marker="*",
        )
    ax.annotate(text, (operation_dates[i], optim_portfo_evolution[i]), fontsize=10)
    ax.annotate(text, (operation_dates[i], index_op_evolution[i]), fontsize=10)

# Plot "index-no-op" (Only on the first day of the first operation round and the last day of the last operation round)
(no_transact,) = ax.plot(
    [operation_dates[0], operation_dates[-1]],
    index_no_op_evolution,
    marker="o",
    color=color_index,
)

ax.legend(
    [optimal, index, no_transact],
    ["optim portfo", "index-op", "index-no-op"],
    loc="lower right",
)
ax.annotate(
    " S: Start date of round\n E: End date",
    (operation_dates[0], max(optim_portfo_evolution) * 0.99),
    fontsize=10,
)
ax.set_xlabel("Date", fontsize=10)
ax.set_ylabel("Total asset", fontsize=10)
ax.tick_params(labelsize=10)
ax.set_xlim((operation_dates[0] - ONE_DAY, operation_dates[-1] + ONE_DAY))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%m/%d"))

plt.show()