diff --git a/experiments/backtest.py b/experiments/backtest.py index 8373de8..03765b7 100644 --- a/experiments/backtest.py +++ b/experiments/backtest.py @@ -24,9 +24,8 @@ def data_folder(): def load_data() -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, pd.DataFrame]: prices = pd.read_csv(data_folder() / "prices.csv", index_col=0, parse_dates=True) spread = pd.read_csv(data_folder() / "spreads.csv", index_col=0, parse_dates=True) - volume = pd.read_csv(data_folder() / "volumes.csv", index_col=0, parse_dates=True) rf = pd.read_csv(data_folder() / "rf.csv", index_col=0, parse_dates=True).iloc[:, 0] - return prices, spread, volume, rf + return prices, spread, rf @dataclass @@ -37,9 +36,9 @@ class OptimizationInput: prices: pd.DataFrame mean: pd.Series - covariance: pd.DataFrame + chol: np.array + volas: np.array spread: pd.DataFrame - volume: pd.DataFrame quantities: np.ndarray cash: float risk_target: float @@ -59,7 +58,7 @@ def run_backtest( weights and then execute the trades at time t. """ - prices, spread, volume, rf = load_data() + prices, spread, rf = load_data() n_assets = prices.shape[1] lookback = 500 @@ -81,12 +80,15 @@ def run_backtest( .dropna() ) # At time t includes data up to t+1 covariance_df = returns.ewm(halflife=125).cov() # At time t includes data up to t - days = returns.index + indices = range(lookback, len(prices) - forward_smoothing) + days = [prices.index[t] for t in indices] covariances = {} + cholesky_factorizations = {} for day in days: covariances[day] = covariance_df.loc[day] + cholesky_factorizations[day] = np.linalg.cholesky(covariances[day].values) - for t in range(lookback, len(prices) - forward_smoothing): + for t in indices: start_time = time.perf_counter() day = prices.index[t] @@ -96,17 +98,18 @@ def run_backtest( prices_t = prices.iloc[t - lookback : t + 1] # Up to t spread_t = spread.iloc[t - lookback : t + 1] - volume_t = volume.iloc[t - lookback : t + 1] mean_t = means.loc[day] # Forecast for return t to t+1 covariance_t = covariances[day] # Forecast for covariance t to t+1 + chol_t = cholesky_factorizations[day] + volas_t = np.sqrt(np.diag(covariance_t.values)) inputs_t = OptimizationInput( prices_t, mean_t, - covariance_t, + chol_t, + volas_t, spread_t, - volume_t, quantities, cash, risk_target, @@ -185,7 +188,10 @@ def interest_and_fees( cash_interest = cash * (1 + rf) ** days_t_to_t_minus_1 - cash short_valuations = np.clip(quantities, None, 0) * prices short_value = short_valuations.sum() - shorting_fee = short_value * (1 + rf) ** days_t_to_t_minus_1 - short_value + short_spread = 0.05 / 360 + shorting_fee = ( + short_value * (1 + rf + short_spread) ** days_t_to_t_minus_1 - short_value + ) return cash_interest + shorting_fee @@ -249,12 +255,16 @@ def asset_weights(self): return self.valuations.div(self.portfolio_value, axis=0) @property - def turnover(self) -> float: + def daily_turnover(self) -> pd.Series: trades = self.quantities.diff() prices = load_data()[0].loc[self.history] valuation_trades = trades * prices relative_trades = valuation_trades.div(self.portfolio_value, axis=0) - return relative_trades.abs().sum(axis=1).mean() * self.periods_per_year + return relative_trades.abs().sum(axis=1) + + @property + def turnover(self) -> float: + return self.daily_turnover.mean() * self.periods_per_year @property def mean_return(self) -> float: @@ -274,7 +284,7 @@ def max_leverage(self) -> float: @property def sharpe(self) -> float: - risk_free = load_data()[3].loc[self.history] + risk_free = load_data()[2].loc[self.history] excess_return = self.portfolio_returns - risk_free return ( excess_return.mean() / excess_return.std() * np.sqrt(self.periods_per_year) diff --git a/experiments/scaling_large.py b/experiments/scaling_large.py new file mode 100644 index 0000000..dd7a6d6 --- /dev/null +++ b/experiments/scaling_large.py @@ -0,0 +1,160 @@ +import cvxpy as cp +import numpy as np +import pandas as pd +from utils import generate_random_inputs + + +def main(): + fitting = True + scenarios = get_scenarios(fitting=fitting) + res = [] + for n_assets, n_factors in scenarios: + print(f"Running scenario with {n_assets} assets and {n_factors} factors") + solvers = [cp.CLARABEL] if fitting else [cp.MOSEK, cp.CLARABEL] + for solver in solvers: + for _ in range(1): + problem = run_scaling(n_assets, n_factors, solver) + assert problem.status in { + cp.OPTIMAL, + cp.OPTIMAL_INACCURATE, + }, problem.status + + res.append( + { + "n_assets": n_assets, + "n_factors": n_factors, + "solve_time": problem.solver_stats.solve_time, + "solver": solver, + } + ) + + df = pd.DataFrame(res) + + df = df.groupby(["n_assets", "n_factors", "solver"]).mean().reset_index() + + if fitting: + # Estimate the scaling exponents as solve time \approx a * n_assets^b * n_factors^c + n_assets = df["n_assets"].values + n_factors = df["n_factors"].values + log_solve_time = np.log(df["solve_time"].values) + + a = cp.Variable() + b = cp.Variable() + c = cp.Variable() + + objective = cp.Minimize( + cp.sum_squares( + a + b * np.log(n_assets) + c * np.log(n_factors) - log_solve_time + ) + ) + problem = cp.Problem(objective) + problem.solve() + assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}, problem.status + + print( + f"Estimated scaling exponents: a={np.exp(a.value):.2f}, \ + b={b.value:.2f}, c={c.value:.2f}" + ) + + else: + df.set_index(["n_assets", "n_factors"], inplace=True) + df = df.pivot(columns="solver", values="solve_time") + df = df.map(lambda x: f"{x:.2f}") + df = df.loc[:, [cp.MOSEK, cp.CLARABEL]] + + # Reset column and row indices + df.reset_index(inplace=True) + df.columns.name = None + df.index.name = None + + print(df.to_latex(index=False)) + + +def run_scaling( + n_assets: int, n_factors: int, solver: str +) -> tuple[np.ndarray, float, cp.Problem]: + mean, F, covariance = generate_random_inputs(n_assets, n_factors) + factor_chol = np.linalg.cholesky(covariance) + factor_volas = np.diag(factor_chol) + + equal_weights = np.ones(n_assets) / n_assets + np.sqrt(equal_weights @ F @ covariance @ F.T @ equal_weights) + sigma_target = 0 + + # The risk constraint is soft. + # For each percentage point of risk, we need to compensate with + # 5 percentage points of return. + + rho_mean = np.percentile(np.abs(mean), 20, axis=0) * np.ones(n_assets) + rho_covariance = 0.02 + L_max = 1.6 + T_max = 50 / 252 + + risk_free = 0.0001 + w_lower = np.ones(n_assets) * (-0.05) + w_upper = np.ones(n_assets) * 0.1 + c_lower = -0.05 + c_upper = 1.0 + gamma_risk = 5.0 + + w_prev = np.ones(n_assets) / n_assets + c_prev = 0.0 + + w, c = cp.Variable(n_assets), cp.Variable() + + z = w - w_prev + T = cp.norm1(z) + L = cp.norm1(w) + + # worst-case (robust) return + mean_return = w @ mean + risk_free * c + return_uncertainty = rho_mean @ cp.abs(w) + return_wc = mean_return - return_uncertainty + + # worst-case (robust) risk + risk = cp.norm2((F @ factor_chol).T @ w) + factor_volas = cp.norm2(F @ factor_chol, axis=1) + + risk_uncertainty = rho_covariance**0.5 * factor_volas @ cp.abs(w) + risk_wc = cp.norm2(cp.hstack([risk, risk_uncertainty])) + + objective = return_wc - gamma_risk * cp.pos(risk_wc - sigma_target) + + constraints = [ + cp.sum(w) + c == 1, + c == c_prev - cp.sum(z), + c_lower <= c, + c <= c_upper, + w_lower <= w, + w <= w_upper, + L <= L_max, + T <= T_max, + ] + + problem = cp.Problem(cp.Maximize(objective), constraints) + problem.solve(solver=solver, verbose=False) + + assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}, problem.status + return problem + + +def get_scenarios(fitting=False): + if not fitting: + return [ + (100, 10), + (500, 30), + (500, 50), + (2000, 50), + (2000, 100), + (10000, 50), + (10000, 100), + ] + else: + # fine grid for fitting + assets = np.logspace(2, 3, 10, dtype=int) + pairs = [(a, int(a * f)) for a in assets for f in np.logspace(-2, -1, 10)] + return pairs + + +if __name__ == "__main__": + main() diff --git a/experiments/scaling_small.py b/experiments/scaling_small.py new file mode 100644 index 0000000..1e910e1 --- /dev/null +++ b/experiments/scaling_small.py @@ -0,0 +1,269 @@ +from functools import lru_cache +import time +from matplotlib import pyplot as plt +import cvxpy as cp +import numpy as np +from backtest import BacktestResult, OptimizationInput, Timing, load_data, run_backtest +from experiments.utils import get_solver + + +def scaling_markowitz( + inputs: OptimizationInput, +) -> tuple[np.ndarray, float, cp.Problem]: + n_assets = inputs.n_assets + latest_prices = inputs.prices.iloc[-1] + portfolio_value = inputs.cash + inputs.quantities @ latest_prices + + # The risk constraint is soft. + # For each percentage point of risk, we need to compensate with + # 5 percentage points of return. + + rho_mean = np.percentile(np.abs(inputs.mean.values), 20, axis=0) * np.ones( + inputs.n_assets + ) + rho_covariance = 0.02 + L_max = 1.6 + T_max = 50 / 252 + + w_lower = np.ones(inputs.n_assets) * (-0.05) + w_upper = np.ones(inputs.n_assets) * 0.1 + c_lower = -0.05 + c_upper = 1.0 + gamma_risk = 5.0 + + w_prev = (inputs.quantities * latest_prices / portfolio_value).values + c_prev = inputs.cash / portfolio_value + + w, c = cp.Variable(n_assets), cp.Variable() + + z = w - w_prev + T = cp.norm1(z) + L = cp.norm1(w) + + # worst-case (robust) return + mean_return = w @ inputs.mean.values + inputs.risk_free * c + return_uncertainty = rho_mean @ cp.abs(w) + return_wc = mean_return - return_uncertainty + + # worst-case (robust) risk + risk = cp.norm2(inputs.chol.T @ w) + risk_uncertainty = rho_covariance**0.5 * inputs.volas @ cp.abs(w) + risk_wc = cp.norm2(cp.hstack([risk, risk_uncertainty])) + + objective = return_wc - gamma_risk * cp.pos(risk_wc - inputs.risk_target) + + constraints = [ + cp.sum(w) + c == 1, + c == c_prev - cp.sum(z), + c_lower <= c, + c <= c_upper, + w_lower <= w, + w <= w_upper, + L <= L_max, + T <= T_max, + ] + + problem = cp.Problem(cp.Maximize(objective), constraints) + problem.solve(solver=get_solver()) + assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}, problem.status + return w.value, c.value, problem + + +def parameter_scaling_markowitz( + inputs: OptimizationInput, +) -> tuple[np.ndarray, float, cp.Problem]: + problem, param_dict, w, c = get_parametrized_problem( + inputs.n_assets, inputs.risk_target + ) + latest_prices = inputs.prices.iloc[-1] + portfolio_value = inputs.cash + inputs.quantities @ latest_prices + + param_dict["chol"].value = inputs.chol + param_dict["volas"].value = inputs.volas + param_dict["rho_mean"].value = np.percentile( + np.abs(inputs.mean.values), 20, axis=0 + ) * np.ones(inputs.n_assets) + param_dict["w_prev"].value = ( + inputs.quantities * inputs.prices.iloc[-1] / portfolio_value + ).values + param_dict["c_prev"].value = inputs.cash / portfolio_value + param_dict["mean"].value = inputs.mean.values + param_dict["risk_free"].value = inputs.risk_free + + problem.solve(solver=get_solver()) + assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}, problem.status + return w.value, c.value, problem + + +@lru_cache +def get_parametrized_problem( + n_assets: int, risk_target: float +) -> tuple[cp.Problem, dict, cp.Variable, cp.Variable]: + rho_covariance = 0.02 + L_max = 1.6 + T_max = 50 / 252 + + w_lower = np.ones(n_assets) * (-0.05) + w_upper = np.ones(n_assets) * 0.1 + c_lower = -0.05 + c_upper = 1.0 + gamma_risk = 5.0 + + w_prev = cp.Parameter(n_assets) + c_prev = cp.Parameter() + mean = cp.Parameter(n_assets) + risk_free = cp.Parameter() + rho_mean = cp.Parameter(n_assets) + chol = cp.Parameter((n_assets, n_assets)) + volas = cp.Parameter(n_assets, nonneg=True) + + w, c = cp.Variable(n_assets), cp.Variable() + + z = w - w_prev + T = cp.norm1(z) + L = cp.norm1(w) + + # worst-case (robust) return + mean_return = w @ mean + risk_free * c + abs_weight_var = cp.Variable(n_assets, nonneg=True) + return_uncertainty = rho_mean @ abs_weight_var + return_wc = mean_return - return_uncertainty + + # worst-case (robust) risk + risk = cp.norm2(chol.T @ w) + risk_uncertainty = rho_covariance**0.5 * volas @ abs_weight_var + risk_wc = cp.norm2(cp.hstack([risk, risk_uncertainty])) + + objective = return_wc - gamma_risk * cp.pos(risk_wc - risk_target) + + constraints = [ + cp.sum(w) + c == 1, + c == c_prev - cp.sum(z), + c_lower <= c, + c <= c_upper, + w_lower <= w, + w <= w_upper, + L <= L_max, + T <= T_max, + cp.abs(w) <= abs_weight_var, + ] + + problem = cp.Problem(cp.Maximize(objective), constraints) + + param_dict = { + "w_prev": w_prev, + "mean": mean, + "risk_free": risk_free, + "rho_mean": rho_mean, + "chol": chol, + "volas": volas, + "c_prev": c_prev, + } + return problem, param_dict, w, c + + +def plot_timings(timings: list[Timing], specifier: str = "") -> None: + # Stacked area plot of cvxpy, solver, and other times + colors = ["#1f77b4", "#ff7f0e", "#2ca02c"] + # darker than darkgrey + line_color = "#888888" + plt.figure() + plt.stackplot( + range(len(timings)), + [timing.cvxpy for timing in timings], + [timing.solver for timing in timings], + [timing.other for timing in timings], + labels=["CVXPY", "Solver", "Other"], + colors=colors, + ) + + # add light horizontal line for average solver time + average_cvxpy_time = np.mean([timing.cvxpy for timing in timings]) + average_solver_time = np.mean([timing.solver for timing in timings]) + average_other_time = np.mean([timing.other for timing in timings]) + + plt.axhline( + average_cvxpy_time, + color=line_color, + linestyle="--", + ) + + plt.axhline( + average_solver_time + average_cvxpy_time, + color=line_color, + linestyle="--", + ) + + plt.axhline( + average_other_time + average_solver_time + average_cvxpy_time, + color=line_color, + linestyle="--", + ) + + plt.xlabel("Day of backtest") + plt.xlim(0, len(timings)) + + plt.ylabel("Time (s)") + plt.legend() + plt.savefig(f"figures/timing{specifier}.png") + plt.show() + + +def main(from_checkpoint: bool = True) -> None: + annualized_target = 0.10 + sigma_target = annualized_target / np.sqrt(252) + + if not from_checkpoint: + print("Running scaling") + # scaling_markowitz_result = run_backtest(scaling_markowitz, sigma_target, verbose=True) + # scaling_markowitz_result.save(f"checkpoints/scaling_{annualized_target}.pickle") + + print("Running parameter scaling") + + n_assets = load_data()[0].shape[1] + start = time.perf_counter() + problem, param_dict, _, _ = get_parametrized_problem(n_assets, sigma_target) + try: + for p in param_dict.values(): + p.value = np.zeros(p.shape) + problem.solve(solver=get_solver()) + except cp.SolverError: + pass + end = time.perf_counter() + print(f"First call to get_parametrized_problem took {end-start} seconds") + + scaling_parametrized_markowitz_result = run_backtest( + parameter_scaling_markowitz, sigma_target, verbose=True + ) + scaling_parametrized_markowitz_result.save( + f"checkpoints/scaling_parametrized_{annualized_target}.pickle" + ) + else: + # scaling_markowitz_result = BacktestResult.load( + # f"checkpoints/scaling_{annualized_target}.pickle" + # ) + scaling_parametrized_markowitz_result = BacktestResult.load( + f"checkpoints/scaling_parametrized_{annualized_target}.pickle" + ) + + # plot_timings(scaling_markowitz_result.timings) + total_time = sum(t.total for t in scaling_parametrized_markowitz_result.timings) + cvxpy_time = sum(t.cvxpy for t in scaling_parametrized_markowitz_result.timings) + solver_time = sum(t.solver for t in scaling_parametrized_markowitz_result.timings) + other_time = sum(t.other for t in scaling_parametrized_markowitz_result.timings) + print(f"Total time {total_time}") + print(len(scaling_parametrized_markowitz_result.timings)) + print( + f"Average time {total_time / len(scaling_parametrized_markowitz_result.timings)}" + ) + print( + f"Average solver time {solver_time / len(scaling_parametrized_markowitz_result.timings)}" + ) + print(f"CVXPY time {cvxpy_time/total_time}") + print(f"Solver time {solver_time/total_time}") + print(f"Other time {other_time/total_time}") + plot_timings(scaling_parametrized_markowitz_result.timings, "_parametrized") + + +if __name__ == "__main__": + main() diff --git a/experiments/taming.py b/experiments/taming.py index 963fe29..22b0266 100644 --- a/experiments/taming.py +++ b/experiments/taming.py @@ -2,22 +2,20 @@ import pandas as pd import cvxpy as cp from backtest import BacktestResult, OptimizationInput, run_backtest +from utils import get_solver from markowitz import Data, Parameters, markowitz -def basic_markowitz(inputs: OptimizationInput) -> np.ndarray: +def basic_markowitz(inputs: OptimizationInput) -> tuple[np.ndarray, float, cp.Problem]: """Compute the basic Markowitz portfolio weights.""" - mu, Sigma = inputs.mean.values, inputs.covariance.values - w = cp.Variable(inputs.n_assets) c = cp.Variable() - objective = mu @ w + objective = inputs.mean.values @ w - chol = np.linalg.cholesky(Sigma) constraints = [ cp.sum(w) + c == 1, - cp.norm2(chol.T @ w) <= inputs.risk_target, + cp.norm2(inputs.chol.T @ w) <= inputs.risk_target, ] problem = cp.Problem(cp.Maximize(objective), constraints) @@ -26,7 +24,7 @@ def basic_markowitz(inputs: OptimizationInput) -> np.ndarray: return w.value, c.value, problem -def equal_weights(inputs: OptimizationInput) -> np.ndarray: +def equal_weights(inputs: OptimizationInput) -> tuple[np.ndarray, float, cp.Problem]: """Compute the equal weights portfolio.""" n_assets = inputs.prices.shape[1] w = np.ones(n_assets) / (n_assets + 1) @@ -34,7 +32,9 @@ def equal_weights(inputs: OptimizationInput) -> np.ndarray: return w, c, None -def weight_limits_markowitz(inputs: OptimizationInput) -> np.ndarray: +def weight_limits_markowitz( + inputs: OptimizationInput, +) -> tuple[np.ndarray, float, cp.Problem]: lb = np.ones(inputs.n_assets) * (-0.05) ub = np.ones(inputs.n_assets) * 0.1 @@ -49,21 +49,25 @@ def weight_limits_markowitz(inputs: OptimizationInput) -> np.ndarray: return markowitz(data, param) -def leverage_limit_markowitz(inputs: OptimizationInput) -> np.ndarray: +def leverage_limit_markowitz( + inputs: OptimizationInput, +) -> tuple[np.ndarray, float, cp.Problem]: data, param = get_basic_data_and_parameters(inputs) param.L_max = 1.6 return markowitz(data, param) -def turnover_limit_markowitz(inputs: OptimizationInput) -> np.ndarray: +def turnover_limit_markowitz( + inputs: OptimizationInput, +) -> tuple[np.ndarray, float, cp.Problem]: data, param = get_basic_data_and_parameters(inputs) param.T_max = 50 / 252 # Maximum TO per year return markowitz(data, param) -def robust_markowitz(inputs: OptimizationInput) -> np.ndarray: +def robust_markowitz(inputs: OptimizationInput) -> tuple[np.ndarray, float, cp.Problem]: data, param = get_basic_data_and_parameters(inputs) param.rho_mean = np.percentile(np.abs(inputs.mean.values), 20, axis=0) * np.ones( inputs.n_assets @@ -90,7 +94,7 @@ def get_basic_data_and_parameters( idio_mean=np.zeros(n_assets), factor_mean=inputs.mean.values, risk_free=0, - factor_covariance_chol=np.linalg.cholesky(inputs.covariance.values), + factor_covariance_chol=inputs.chol, idio_volas=np.zeros(n_assets), F=np.eye(n_assets), kappa_short=np.zeros(n_assets), @@ -118,7 +122,7 @@ def get_basic_data_and_parameters( return data, param -def main(from_checkpoint: bool = False): +def main(from_checkpoint: bool = True) -> None: annualized_target = 0.10 if not from_checkpoint: @@ -244,9 +248,5 @@ def generate_table( ) -def get_solver(): - return cp.MOSEK if cp.MOSEK in cp.installed_solvers() else cp.CLARABEL - - if __name__ == "__main__": main() diff --git a/experiments/utils.py b/experiments/utils.py index 75d8251..e6cbda4 100644 --- a/experiments/utils.py +++ b/experiments/utils.py @@ -1,4 +1,5 @@ import numpy as np +import cvxpy as cp import pandas as pd @@ -28,6 +29,24 @@ def synthetic_returns( return synthetic_returns +def generate_random_inputs( + n_assets: int, n_factors: int +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + rng = np.random.default_rng(1) + mean = rng.normal(0, 0.05, size=n_assets) + + A = rng.uniform(-1, 1, size=(n_factors, n_factors)) + covariance = A @ A.T + + loadings = rng.normal(0, 1, size=(n_assets, n_factors)) + + return mean, loadings, covariance + + +def get_solver(): + return cp.MOSEK if cp.MOSEK in cp.installed_solvers() else cp.CLARABEL + + if __name__ == "__main__": prices = pd.read_csv("data/prices.csv", index_col=0, parse_dates=True) synthetic_returns = synthetic_returns(