# Markowitz Portfolio Numerical Experiments

## Data and Backtests
The experiment uses a dataset of daily adjusted closing prices of stocks in 18 ETFs and the 3-month U.S. Treasury bill rate. This dataset is divided into two periods:

* **In-sample data:** From January 1st, 2018, to June 30th, 2023, used for model training.
* **Out-of-sample data:** From July 1st, 2023, to June 31st, 2024, used for evaluating the model's performance on unseen data.

The 18 ETFs chosen cover a broad range of asset classes, including:

| Ticker | Name | Description | Asset Class |
|--------|------|-------------|-------------|
| AGG | iShares Core U.S. Aggregate Bond ETF | Invests in U.S. investment-grade bonds | Fixed Income |
| DBC | Invesco DB Commodity Index Tracking Fund | Tracks various commodities | Commodities |
| GLD | SPDR Gold Shares | Invests in gold bullion | Commodities |
| IBB | iShares Biotechnology ETF | Tracks biotechnology stocks | Equity |
| ITA | iShares U.S. Aerospace & Defense ETF | Tracks U.S. aerospace and defense companies | Equity |
| PBJ | Invesco Dynamic Food & Beverage ETF | Tracks U.S. food and beverage companies | Equity |
| IEF | iShares 7-10 Year Treasury Bond ETF | Invests in U.S. Treasury bonds (7-10 years) | Fixed Income |
| VNQ | Vanguard Real Estate ETF | Invests in REITs | Real Estate |
| VTI | Vanguard Total Stock Market ETF | Tracks broad U.S. stock market | Equity |
| XLB | Materials Select Sector SPDR Fund | Tracks materials sector stocks | Equity |
| XLE | Energy Select Sector SPDR Fund | Tracks energy sector stocks | Equity |
| XLF | Financial Select Sector SPDR Fund | Tracks financial sector stocks | Equity |
| XLI | Industrial Select Sector SPDR Fund | Tracks industrial sector stocks | Equity |
| XLK | Technology Select Sector SPDR Fund | Tracks technology sector stocks | Equity |
| XLP | Consumer Staples Select Sector SPDR Fund | Tracks consumer staples sector stocks | Equity |
| XLU | Utilities Select Sector SPDR Fund | Tracks utilities sector stocks | Equity |
| XLV | Health Care Select Sector SPDR Fund | Tracks healthcare sector stocks | Equity |
| XLY | Consumer Discretionary Select Sector SPDR Fund | Tracks consumer discretionary sector stocks | Equity |

### Import Packages and Functions

import matplotlib.pyplot as plt
import pandas as pd
import pandas_datareader as pdr
import yfinance
import datetime
import cvxportfolio as cvx
import numpy as np
import cvxpy as cp
import seaborn as sn
import quandl

### Load the Data

# Data Source: Yahoo Finance and FRED
def load_data():
    tickers = ["AGG", "DBC", "GLD", "IBB", "ITA", "PBJ", "IEF", "VNQ", "VTI", "XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY"]
    start_date = '2007-01-01'
    end_date = pd.Timestamp.today()

    returns = pd.DataFrame(dict([(ticker, yfinance.download(ticker, start=start_date, end=end_date)['Adj Close'].pct_change()) for ticker in tickers]))
    returns["USDOLLAR"] = quandl.get('FRED/DTB3', start_date=start_date, end_date=end_date)/(250*100)
    returns = returns.fillna(method='ffill').iloc[1:]

    
    return start_date, end_date, returns

start_date, end_date, returns = load_data()
market_data = cvx.UserProvidedMarketData(returns=returns, cash_key='USDOLLAR')


## Estimating Return and Covariance

### Return Prediction
To simulate a realistic yet proprietary return prediction method, we generate synthetic returns using the following model:

$$
\hat{r}_t = \alpha(r_t+ \epsilon_t)
$$

where:
* $r_t$ is the actual (realized) return of the asset at time $t$
* $\hat{r}_t$ is the predicted return of the asset at time $t$.
* $\epsilon_t$ is zero-mean Gaussian noise with variance $\sigma_\epsilon^2I$, representing the uncertainty in the prediction.
* We set the noise standard deviation $\sigma_\epsilon=0.02$ (14%), reflecting a moderately noisy forecast.
* $\alpha$ is a scaling factor derived from the square of the Information Coefficient (IC), a measure of the correlation between predicted and actual returns.
* The IC ranges from -1 to 1, with higher values indicating better predictive skill.
* A higher IC indicates a stronger relationship and a more skilled forecast.
* We set the annualized information ratio (IR), which is the IC divided by its tracking error, to 0.15. This implies a typical forecast error of approximately ±0.3% per day.

**Note:** The scaling factor $\alpha$ is chosen to minimize the mean squared error (MSE) between predicted and actual returns:

$$
\alpha = \frac{\sigma_r^2}{\sigma_r^2 + \sigma^2_\epsilon}
$$

where $\sigma_r^2=0.0005$ is the variance of realized returns, corresponding to a daily volatility of roughly 2%.

### Covariance Prediction
* We estimate the covariance matrix of asset returns using an Exponentially Weighted Moving Average (EWMA) model. This model assigns greater weight to recent observations, allowing it to adapt to changing market conditions.
* We use a half-life of 125 trading days, meaning that the influence of an observation halves after 125 days. The corresponding decay factor $\beta$ is calculated as: $$\beta = 1 − \frac{\ln 2}{\text{halflife}} \approx 0.994$$

# Return Estimate
def synthetic_returns(returns, information_ratio, forward_smoothing):
    rng = np.random.default_rng(1)
    returns = returns.rolling(forward_smoothing).mean().shift(-(forward_smoothing - 1))
    var_r = returns.var()

    alpha = information_ratio**2
    var_eps = var_r * (1 - alpha) / alpha
    noise = rng.normal(0, np.sqrt(var_eps), size=returns.shape)
    synthetic_returns = alpha * (returns + noise)
    return synthetic_returns

r_hat = synthetic_returns(returns, information_ratio=0.15, forward_smoothing=5).shift(1).dropna().iloc[:, :-1]

# Covariance estimate
covariance = returns.shift(1).ewm(halflife=125).cov().dropna()

## Single-Period Markowitz Optimization
We explore various Markowitz portfolio optimization strategies, targeting an annualized volatility of 10% ($\sqrt{250}\sigma = 0.1$), corresponding to a daily variance of 0.00004.

### Initial Setup
We begin with an initial capital of $1,000,000 held entirely in cash. The portfolio is rebalanced daily to target weights, considering previous returns and executing trades at closing prices. To accurately model borrowing costs, we use the daily equivalent of the 3-Month Treasury Bill Secondary Market Rate as the borrow fee.

# Risk constraint
risk_target = (cvx.FullCovariance(covariance) <= 0.00004)

# Transaction and holding cost models
HALF_SPREAD = 10E-4  

# Annualized borrow fee
r_hat_with_cash = market_data.returns.rolling(window=250).mean().shift(1).dropna()
BORROW_FEE = r_hat_with_cash.iloc[:, -1]

tcost_model = cvx.TcostModel(a=HALF_SPREAD, b=None)
hcost_model = cvx.HcostModel(short_fees=BORROW_FEE)

# Market simulator (daily trading)
market_sim = cvx.MarketSimulator(
    market_data=market_data,
    costs=[tcost_model, hcost_model]
)

# Initial portfolio (uniform on non-cash assets, with $1M total capital)
init_portfolio = pd.Series(index=market_data.returns.columns, data=1000000/18)
init_portfolio.USDOLLAR = 0

### Equally Weighted Portfolio (Benchmark)
As a benchmark, we use an equally weighted portfolio, allocating equal capital to each of the 18 ETFs.

# Define equally weight setup
w_b = pd.Series(index=market_data.returns.columns, data=1)
w_b.USDOLLAR = 0.
w_b /= sum(w_b)
target_weights = w_b
rebalancing_times = pd.date_range('2018-01-01', '2024-06-30', freq='M')

spo_policy_equal_weight = cvx.PeriodicRebalance(target=target_weights, rebalancing_times=rebalancing_times)

### Basic Markowitz
The classic Markowitz optimization seeks to maximize expected returns while staying within the risk constraint:

\begin{split}
\text{minimize} &\quad \mu^T w \\
\text{subject to} &\quad w^T\Sigma w \le (\sigma^{tar})^2, \ 1^Tw = 1
\end{split}

spo_policy_base = cvx.SinglePeriodOpt(
    objective=cvx.ReturnsForecast(r_hat),
    constraints=[risk_target],
    include_cash_return=False
)

# Run backtest for Basic Markowitz and Equally weight
results_base = market_sim.run_multiple_backtest(
    h=[init_portfolio]*2,
    start_time='2018-01-03',  end_time='2023-06-30',
    policies=[spo_policy_base, spo_policy_equal_weight]
)

# Plot comparison
fig = plt.figure(figsize=(10,8))
results_base[0].v.plot(label="Basic Markowitz")
results_base[1].v.plot(label="Equally Weighted Portfolio")
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.title('Basic Markowitz vs Equally Weighted Portfolio')
plt.legend()


### Taming Markowitz
The classic Markowitz model, while theoretically elegant, often faces challenges in real-world applications due to its sensitivity to estimation errors and its tendency to produce overly concentrated portfolios. To address these limitations, we explore several techniques to "tame" Markowitz by incorporating additional constraints and adjustments to the objective function:

#### Markowitz with Regularization
We start by incorporating regularization terms into the objective function. These terms penalize both trading costs and holding costs, where the penalties are scaled by factors $\gamma^{hold}$ and $\gamma^{trade}$, respectively. This encourages the optimizer to seek solutions that not only maximize returns but also minimize unnecessary trading and the costs associated with maintaining portfolio positions. In this initial exploration, we set both penalty factors to 1.

gamma_trade, gamma_hold = 1., 1.
spo_policy_reg = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) 
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[risk_target],
    include_cash_return=False
)

#### Additional Constraints
Next, we introduce several types of constraints to further refine the portfolio optimization process:

##### Weight Limits
We impose constraints on the maximum and minimum allocation to each asset, setting a long position limit of 40% and a short position limit of -25%. This promotes diversification and helps mitigate the risk of overly concentrated positions.

max_weight = cvx.MaxWeights(0.4)
min_weight = cvx.MinWeights(-0.25)
spo_policy_w = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) 
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[max_weight, min_weight, risk_target],
    include_cash_return=False
)

##### Leverage Limit
We set a leverage limit to restrict the total exposure of the portfolio. A maximum leverage of 2 is chosen, consistent with the original paper.

l_limit = cvx.LeverageLimit(2)
spo_policy_l = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) 
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[l_limit, risk_target],
    include_cash_return=False
)

##### Turnover Limit
To manage transaction costs, we restrict the portfolio turnover to 25%. This limits the amount of trading activity and helps preserve returns that might otherwise be eroded by excessive trading fees.

turnover_limit = cvx.TurnoverLimit(0.25)
spo_policy_t = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) 
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[turnover_limit, risk_target],
    include_cash_return=False
)

#### Robustness to Forecasts
We also acknowledge the inherent uncertainty in predicting returns and covariance. To address this, we introduce robustness measures:

* **Robust Return Forecast:** We adjust the expected return forecast by adding a penalty term that accounts for the uncertainty in the forecast. This penalty is proportional to the absolute value of the portfolio weights, scaled by a factor `rho_mean`, which represents the 20th percentile of the absolute value of the return forecast error.
* **Robust Covariance Forecast:** We augment the covariance matrix with an additional term that accounts for the uncertainty in the covariance estimate. This term is proportional to the squared weighted sum of individual asset volatilities, scaled by a factor `rho_covariance`.

rho_mean = np.percentile(np.abs(r_hat), 20, axis=0) * np.ones(r_hat.shape[1])
rho_covariance = 0.1
rho_mean = pd.Series(rho_mean, index=r_hat.columns)
variance = returns.ewm(halflife=125).var().shift(1).dropna().iloc[:, :-1]
risk_target_robust = risk_target_robost = ((cvx.FullCovariance(covariance) + rho_covariance * cvx.RiskForecastError(variance)) <= 0.00004)

spo_policy_robust = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) - cvx.ReturnsForecastError(deltas = rho_mean)
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[risk_target_robust],
    include_cash_return=False
)

##### Markowitz with All Constraints and Robustness
Finally, we combine all the aforementioned constraints and robustness adjustments into a single optimization problem. The resulting strategy, referred to as "Markowitz with All Constraints," aims to find an optimal portfolio that balances return maximization with risk management, transaction cost control, and robustness to forecast errors.

spo_policy_all = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) -  cvx.ReturnsForecastError(deltas = rho_mean)
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model,
    constraints=[risk_target_robust, turnover_limit, l_limit, max_weight, min_weight],
    include_cash_return=False
)

#### Comparison
We now compare the performance of the various Markowitz strategies we've defined, along with the benchmark equally weighted portfolio.

# Run backtest for all strategies
results = market_sim.run_multiple_backtest(
    h=[init_portfolio] * 5,
    start_time='2018-01-03',  
    end_time='2023-06-30',
    policies=[spo_policy_l, spo_policy_w, spo_policy_t, spo_policy_robost, spo_policy_all]
)

# Print summary statistics
print(results)

# Plot performance
plt.figure(figsize=(12, 8))
results.plot()
plt.title('Performance Comparison of Markowitz Variants')
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.legend(strategies.keys())
plt.tight_layout()
plt.show()

# Calculate and display additional metrics
additional_metrics = results.estimate_additional_metrics()
print(additional_metrics)

| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Mininal Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight | 8.2% |13.4% | 0.48 | 0.1 | 1 |-15.6% |
| Basic Markowitz | -19.6% | 11.0% | -1.93 | 73.5 | 0.87 |-68.2% |
| Markowitz with regularization | -8.5% | 11.6% | -0.86| 77.0 |16.23|-43.2% |
| Leverage-limited |7.7% |11.0%| 0.56 | 55.8 | 1.98 | -14.6% |
| Weight-limited |11.1%| 10.8% |0.88 |50.8 | 3.56 | -9.5% |
| Turnover-limited |12.0%| 11.0% |0.95 | 22.9 |4.83 |-13.3% |
| Robust | 10.1% |7.8%| 1.41 | 34.3 | 1.16 | -7.5% |
| Include all constrains | 17.9% |8% | 2.04 | 19.7 | 1.143 | -5.2% |


#### Key Findings
* **Leverage's Impact:** Strategic use of leverage can enhance risk-adjusted returns and reduce portfolio volatility, but excessive leverage poses risks.
* **Weight Constraints' Role:** Imposing limits on individual asset weights can further improve returns, but may lead to unintended high leverage.
* **Turnover Considerations:** Limiting portfolio turnover can boost returns, but doesn't necessarily address leverage concerns.
* **Robust Optimization's Promise:** The most effective approach involves incorporating robustness to forecast uncertainty, leading to superior performance, higher Sharpe ratios, and better risk management compared to other strategies.

### Markowitz++: A Comprehensive Strategy
The Markowitz++ strategy enhances the basic Markowitz model by integrating all the previously discussed techniques to address its practical limitations. By incorporating constraints on leverage, turnover, and individual asset weights, along with robust adjustments to forecasts, Markowitz++ seeks to create a more practical and resilient portfolio.

#### Initial Parameters and Soft Constraints
The initial parameters for Markowitz++ are:

* Holding cost penalty $\gamma^{hold}$: 1
* Trading cost penalty $\gamma^{trade}$: 1
* Target volatility $\sigma^{tar}$: 10% annualized
* Minimum weight $w^{min}$: -25%
* Maximum weight $w^{max}$: 40%
* Target leverage $L^{tar}$: 2
Target turnover $T^{tar}$: 25%

Soft constraints are used to gently enforce limits on leverage, turnover, and weights, with penalty factors of 0.005 each.

$$
\gamma^{w^{max}}= 0.005, \ \gamma^{w^{min}}=0.005, \ \gamma^{lev} = 0.005, \ \gamma^{turn} = 0.005
$$

gamma_trade, gamma_hold = 1., 1.
spo_policy_markowitz_plus = cvx.SinglePeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) - cvx.ReturnsForecastError(deltas = rho_mean)
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model
        - 0.005 * cvx.SoftConstraint(cvx.TurnoverLimit(0.25))
        - 0.005 * cvx.SoftConstraint(cvx.LeverageLimit(2))
        - 0.005 * cvx.SoftConstraint(cvx.MaxWeights(0.40))
        - 0.005 * cvx.SoftConstraint(cvx.MinWeights(-0.25)),
    constraints=[risk_target_robust],
    include_cash_return=False
)

# Run backtest
results_markowitz_plus = market_sim.run_backtest(
    h=[init_portfolio],
    start_time='2018-01-03',  end_time='2023-06-30',
    policies=[spo_policy_markowitz_plus]
)

print(results_markowitz_plus)

| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Markowitz++ | 22.8% | 10% | 2.05 | 21.9| 1.65 | -7.9% |

### Parameter Tuning
To further refine the strategy, we perform a grid search to optimize the penalty factors for holding costs, trading costs, leverage, turnover, and weight constraints. The goal is to maximize the in-sample Sharpe ratio while adhering to realistic constraints on turnover, leverage, and volatility.

#### Criterias of Adjusting a Parameter
* Increase in in-sample Sharpe Ratio
* In-sample annualized turnover $\leq 50 \%$
* In-sample maximum leverage $\leq 2$
* In-sample annualized volatility $\leq 15\%$

#### Tuning Process

import itertools

def evaluate_strategy(params):
    gamma_hold, gamma_trade, gamma_lev, gamma_turn, gamma_w_min, gamma_w_max = params
    
    strategy = cvx.SinglePeriodOpt(
        objective = cvx.ReturnsForecast(r_hat) - rho_covariance * cvx.ReturnsForecastError(deltas = rho_mean)
            - gamma_trade * tcost_model
            - gamma_hold * hcost_model
            - gamma_turn * cvx.SoftConstraint(cvx.TurnoverLimit(0.25))
            - gamma_lev * cvx.SoftConstraint(cvx.LeverageLimit(1.6))
            - gamma_w_max * cvx.SoftConstraint(cvx.MaxWeights(0.50))
            - gamma_w_min * cvx.SoftConstraint(cvx.MinWeights(-0.05)),
        constraints=[risk_target_robust],
        include_cash_return=False
    )
    
    result = market_sim.run_multiple_backtest(
        h=[init_portfolio],
        start_time='2018-01-03',  end_time='2023-06-30',
        policies=[strategy]
    )
    
    metrics = result[0]
    return metrics.sharpe_ratio, metrics.turnover.mean(), metrics.leverage.mean(), metrics.volatility


param_grid = {
    'gamma_hold': [0.5, 1, 2],
    'gamma_trade': [1, 3, 5],
    'gamma_lev': [0.0000001, 0.00005,0.001,0.01,0.25,1],
    'gamma_turn': [0.0000001, 0.00005,0.001,0.01,0.25,1],
    'gamma_w_min': [0.0000001, 0.00005,0.001,0.01,0.25,1],
    'gamma_w_max': [0.0000001, 0.00005,0.001,0.01,0.25,1]
}

best_params = None
best_sharpe = -float('inf')

for params in itertools.product(*param_grid.values()):
    sharpe, turnover, leverage, volatility = evaluate_strategy(params)
    
    if (sharpe > best_sharpe and 
        turnover <= 0.5 and 
        leverage <= 2 and 
        volatility <= 0.15):
        best_sharpe = sharpe
        best_params = params

print("Best parameters:", best_params)
print("Best Sharpe ratio:", best_sharpe)

Best parameters: (0.5, 3, 0.0005, 0.25, 0.00005, 0.0000001)
Best Sharpe ratio: 2.05

#### Backtest Comparison
Let's compare the performance of all strategies, including the tuned Markowitz++:

| -- | Return | Volatility| Sharpe | Turnover | Leverage | Mininal Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight | 8.2% |13.4% | 0.48 | 0.1 | 1 |-15.6% |
| Basic Markowitz | -19.6% | 11.0% | -1.93 | 73.5 | 0.87 |-68.2% |
| Markowitz with regularization | -8.5% | 11.6% | -0.86| 77.0 |16.23|-43.2% |
| Leverage-limited |7.7% |11.0%| 0.56 | 55.8 | 1.98 | -14.6% |
| Weight-limited |11.1%| 10.8% |0.88 |50.8 | 3.56 | -9.5% |
| Turnover-limited |12.0%| 11.0% |0.95 | 22.9 |4.83 |-13.3% |
| Robust | 10.1% |7.8%| 1.41 | 34.3 | 1.16 | -7.5% |
| Include all constrains | 17.9% |8% | 2.04 | 19.7 | 1.143 | -5.2% |
| Markowitz++ | 22.8% | 10% | 2.05 | 21.9| 1.65 | -7.9% |
| Tuned Markowitz++| 23.6% | 10.1% | 2.18 | 22.0 | 1.67 | -7.6% |


### Out-of-Sample Backtest

results_out = market_sim.run_multiple_backtest(
    h=[init_portfolio]*9,
    start_time='2023-07-01',
    end_time='2024-06-30',
    policies=[spo_policy_base,cvx.PeriodicRebalance(target=target_weights, rebalancing_times=rebalancing_times),\
    spo_policy_l, spo_policy_w, spo_policy_t,spo_policy_all,spo_policy_soft,spo_policy_robost,spo_policy_soft_tun])

print(results_out)

plt.figure(figsize=(12, 8))
results_out[0].v.plot(label = "Basic")
results_out[1].v.plot(label = "Equal weight")
results_out[2].v.plot(label = "Leverage-limited")
results_out[3].v.plot(label = "Weight-limit")
results_out[4].v.plot(label = "Turnover-limit")
results_out[5].v.plot(label = "All constrains")
results_out[6].v.plot(label = "Markowitz++")
results_out[7].v.plot(label = "Robust")
results_out[8].v.plot(label = "Tune Markowitz++")
plt.title('Out-of-Sample Performance Comparison (without Robust)')
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.legend()

| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight |11.7% | 8.8% | 0.73 | 0.1 | 1 | -8.7% |
| Basic Markowitz | -22.0% | 8.2% | -3.33 | 68.2 | 0.814 | -20.4% |
| Weight-limited | -0.0% | 9.0%| -0.59 | 46.7| 3.546 |-7.4% |
| Leverage-limited |-7.3% | 9.3% | -1.36| 55.0| 1.98 |-13.2%|
| Turnover-limited | 10.1% |8.4% | 0.57 | 22.5 | 3.712 | -4.6% |
| Robust |2.8% | 6.6% | 0.39 | 34.2 |1.27| -4.8%|
| All | 7.9% |5.9%|0.44 | 20.2| 1.34 | -3.0% |
| Markowitz++ | 8.8% |5.9%| 0.59 | 20.1 | 1.32 | -2.8% |
| Tuned Markowitz++| 13.5% |6.0% | 1.36 | 20.7 | 1.32 |-3.7% |


## Multi-Period Markowitz Optimization
While single-period Markowitz optimization focuses on the immediate next period, multi-period optimization takes a longer view. It aims to find an optimal sequence of portfolio adjustments over multiple periods, considering future trading costs and potential market changes.

### Strategy and Evaluation
We explore multi-period optimization by varying the planning horizon (the number of periods considered in each optimization). The optimization process incorporates the same objective and constraints as the tuned Markowitz++ strategy, but now with the added dimension of time.

lookforward = [2, 3, 5, 10, 20]
policies ={}
for d in lookforward:
    policies[d] = \
    cvx.MultiPeriodOpt(
        objective = cvx.ReturnsForecast(r_hat) - rho_covariance * cvx.ReturnsForecastError(deltas = rho_mean)
            - gamma_trade * tcost_model
            - gamma_hold * hcost_model
            - 0.25 * cvx.SoftConstraint(turnover_limit)
            - 0.05 * cvx.SoftConstraint(l_limit)
            - 0.00001 * cvx.SoftConstraint(max_weight)
            - 0.00005 * cvx.SoftConstraint(min_weight),
        constraints=[risk_target_robost],
        include_cash_return=False,planning_horizon=d)
        
results_multi={}

import warnings
warnings.filterwarnings('ignore')
results_multi.update(dict(zip(policies.keys(), market_sim.run_multiple_backtest(h=[init_portfolio] * len(policies), start_time='2018-01-03',  end_time='2023-06-30',
                                              policies= policies.values(), parallel=True))))

# From the above, we knew that when look ahead period = 5, it has maximum the sharp ratio, lowest volitility and the highest return.
mpo_policy = cvx.MultiPeriodOpt(
    objective = cvx.ReturnsForecast(r_hat) - rho_covariance * cvx.ReturnsForecastError(deltas = rho_mean)#+ cvx.CashReturn(cash_returns = r_hat_cash)
        - gamma_trade * tcost_model
        - gamma_hold * hcost_model
        - 0.25 * cvx.SoftConstraint(turnover_limit)
        - 0.05 * cvx.SoftConstraint(l_limit)
        - 0.00001 * cvx.SoftConstraint(max_weight)
        - 0.00005 * cvx.SoftConstraint(min_weight),
        constraints=[risk_target_robost],
    include_cash_return=False,planning_horizon=5)

# Run backtests for the Multi-period Markowitz++ and Equal Weight policies
results_multi = market_sim.run_multiple_backtest(
    h=[init_portfolio]*3,
    start_time='2018-01-03',  end_time='2023-06-30',
    policies=[spo_policy_soft_tun, mpo_policy, cvx.PeriodicRebalance(target=target_weights, rebalancing_times=rebalancing_times)])
)

# Print statistics
print(results_multi)

### In-Sample and Out-of-Sample Performance
We compare the performance of:

* **Tuned Markowitz++ (Single-Period):** The best-performing single-period strategy.
* **Multi-Period Markowitz++:** The strategy with the optimal planning horizon.
* **Equally Weighted:** The benchmark.

#### In-Sample Results

plt.figure(figsize=(12, 8))
results_multi[0].v.plot(label='Tune Markowitz++')
results_multi[1].v.plot(label='MPO Markowitz')
results_multi[2].v.plot(label='Equal Weight')
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.title('In-Sample Performance: Multi-period Markowitz++ vs Equal Weight')
plt.legend()


| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight | 8.2% |13.4% | 0.48 | 0.1 | 1 |-15.6% |
| Tuned Markowitz++| 23.6% | 10.1% | 2.18 | 22.0 | 1.67 | -7.6% |
| Multi-period Markowitz++ |24.6% | 9.2%| 2.27 | 25| 1.70| -5.2% |

### Out-of-Sample Results

# Run out-of-sample backtests
results_out_of_sample = market_sim.run_multiple_backtest(
    h=[init_portfolio]*3,
    start_time='2023-07-01',
    end_time='2024-06-30',
    policies=[spo_policy_soft_tun, mpo_policy, cvx.PeriodicRebalance(target=target_weights,
    rebalancing_times=rebalancing_times)]
)

# Print statistics for out-of-sample backtests
print(results_out_of_sample)

# Plot out-of-sample performance
plt.figure(figsize=(12, 8))
results_out_of_sample[0].v.plot(label='Tune Markowitz++')
results_out_of_sample[1].v.plot(label='MPO Markowitz')
results_out_of_sample[2].v.plot(label='Equal Weight')
plt.xlabel('Date')
plt.ylabel('Portfolio Value')
plt.title('Out-of-Sample Performance Comparison')
plt.legend()


| Strategy | Return | Volatility| Sharpe | Turnover | Leverage | Drawdown |
| -------- | -------- | -------- | -------- | -------- | -------- | -------- |
| Equal weight |11.7% | 8.8% | 0.73 | 0.1 | 1 | -8.7% |
| Tuned Markowitz++ | 13.5% | 6.0% | 1.36 |20.7| 1.32 |-3.7% |
| Multi-period Markowitz++ |14.3% | 7.7%| 1.56 | 24.9 | 1.78| -3.5% |

### Observations
* **Multi-Period Advantage:** The multi-period optimization strategy consistently outperforms both the single-period and equally weighted strategies, indicating the value of incorporating a longer planning horizon into portfolio decisions.
* **Superior Risk-Adjusted Returns:** The multi-period approach delivers higher Sharpe ratios, suggesting better risk-adjusted performance compared to the other strategies.
* **Trade-Offs:** While the multi-period strategy exhibits slightly higher volatility and turnover, it compensates with higher returns, showcasing a favorable risk-return trade-off.

## Reference
* S. Boyd, K. Johansson, R. Kahn, P. Schiele, and T. Schmelzer Markowitz Portfolio Construction at Seventy, Journal of Portfolio Management, 50(8), 117–160, July 2024.
* [markowitz-reference](https://github.com/cvxgrp/markowitz-reference)
* [cvxportfolio](https://github.com/cvxgrp/cvxportfolio)