# Hierarchical Risk Parity

## Introduction
Investment managers face the daily challenge of constructing portfolios that reflect their views and predictions on risks and returns. Hierarchical Risk Parity (HRP) is an innovative portfolio allocation method designed to address this challenge.

**Key Idea:** HRP combines machine learning techniques with financial theory to build diversified portfolios that aim to perform well under various market conditions.

### Strengths of HRP
* **Robustness:** Doesn't require the covariance matrix to be invertible, unlike traditional mean-variance optimization.
* **Handles Ill-Conditioned Data:** Can work with singular or nearly singular covariance matrices.
* **Reduced Risk:** Empirical studies have shown HRP often produces portfolios with lower out-of-sample variance compared to other methods like the Critical Line Algorithm (CLA).
* **Strong Performance:** HRP portfolios have historically outperformed traditional approaches in various market scenarios.

## HRP Algorithm: A Three-Stage Process
1. **Tree Clustering:** Groups assets based on their correlation, creating a hierarchical structure resembling a tree.
2. **Quasi-Diagonalization:** Reorganizes the covariance matrix to place highly correlated assets together, aiding in risk management.
3. **Recursive Bisection:** Allocates weights to assets within each cluster, balancing risk and return based on an inverse-variance approach.
### Tree Clustering
This stage groups assets with similar risk profiles together. It uses a distance metric based on pairwise correlations:

#### Algorithm
:::success
<font color=blue>Algorithm: Tree Clustering
</font>
Input: $T \times N$ matrix $X$ of observations (e.g., returns series of $N$ variables over $T$ periods)
Output: Hierarchical structure of clusters

1. Compute $N \times N$ correlation matrix $\rho$:
$\rho = (\text{corr}(X_i, X_j))_{i,j = 1,\dots,N}$
2. Compute $N \times N$ distance matrix $D$:
$D = (d_{ij})_{i,j = 1,...,N}$ where $d_{ij} = \sqrt{\frac{1 - \rho_{ij}}{2}}$
3. Compute pairwise Euclidean distances $\tilde{d}$:
$\tilde{d}_{ij} = \sqrt{\sum_{k=1}^N (d_{ki} - d_{kj})^2}$
4. Initialize clusters: $C = \{\{1\}, \{2\}, \dots, \{N\}\}$
5. While $|C| > 1$:
a. Find $(i^*, j^*) = \arg\min_{i,j, i \neq j} {\tilde{d}_{ij}}$
b. Create new cluster $u = C[i^*] ∪ C[j^*]$
c. Update distances: For each cluster $v$ in $C$, $v \neq C[i^*]$, $v \neq C[j^*]$: $\tilde{d}[u,v] = \min(\tilde{d}[C[i^*],v], \tilde{d}[C[j^*],v])$
d. Update $C$: $C = C \backslash \{C[i^*], C[j^*]\} \cup \{u\}$
e. Update $\tilde{d}$ matrix:
- Remove rows and columns for $C[i^*]$ and $C[j^*]$
- Add row and column for u
6. Return final cluster hierarchy
:::

#### Code Implementation

def correlDist(corr):
	# A distance matrix based on correlation, where 0<=d[i,j]<=1
	# This is a proper distance metric
	dist=((1-corr)/2.)**.5 # distance matrix
	return dist

dist=correlDist(corr)
condensed_dist = squareform(dist)
link=sch.linkage(condensed_dist,'single')

**Linkage Matrix:** This matrix describes the hierarchical clustering process, capturing which assets are merged at each step and their corresponding distance.
**Dendrogram Interpretation:** The dendrogram is a tree-like structure where the leaves represent the individual data points, and the branches represent the merges. The height of the merge indicates the dissimilarity between the merged clusters. 

#### Linkage Methods
* **Complete Linkage:** Maximum distance between points in two clusters. $$\max_{x \in C_i,\,y \in C_j}d(x, y)$$
* **Single Linkage:** Maximum distance between points in two clusters. $$\min_{x \in C_i,\,y \in C_j}d(x, y)$$
* **Ward's Method:** Minimizes the increase in within-cluster variance when merging clusters. $$\sum_{x \in C_i \bigcup C_j}d(x, \mu_{C_i \bigcup C_j})$$

### Quasi-Diagonalization
This step reorders the rows and columns of the covariance matrix. The goal is to place the largest elements along the diagonal, which helps to visualize clusters and relationships between assets.

**Key Benefit:** Quasi-diagonalization places similar investments together and dissimilar investments further apart. This can be helpful in understanding risk concentrations within a portfolio.

#### Algorithm
:::success
<font color=blue>Algorithm: Quasi-Diagonalization</font>
Input:
- Linkage matrix `Y` from hierarchical clustering
- Original covariance matrix `Sigma`

Output:
- Reordered covariance matrix `Sigma'`

1. Initialize order list `L = []`
2. Function `QuasiDiagonalize(cluster)`:
* If `cluster` is a leaf node (original asset):
* Append cluster to `L`
* Else:
* `QuasiDiagonalize(cluster.left_child)`
* `QuasiDiagonalize(cluster.right_child)`
3. Start with the root cluster (last row of $Y$):
* `QuasiDiagonalize((Y[N-1, 1], Y[N-1, 2]))`
4. Reorder `Sigma` based on `L`:
`Sigma' = Sigma[L, :][:, L]`
5. Return `Sigma'`
:::

#### Code Implementation

def getQuasiDiag(link):
    # Sort clustered items by distance
    link=link.astype(int)
    sortIx=pd.Series([link[-1,0],link[-1,1]])
    numItems=link[-1,3] # number of original items
    while sortIx.max()>=numItems:
        sortIx.index=range(0,sortIx.shape[0]*2,2) # make space
        df0=sortIx[sortIx>=numItems] # find clusters
        i=df0.index;j=df0.values-numItems
        sortIx[i]=link[j,0] # item 1
        df0=pd.Series(link[j,1],index=i+1)
        sortIx=sortIx.append(df0) # item 2
        sortIx=sortIx.sort_index() # re-sort
        sortIx.index=range(sortIx.shape[0]) # re-index
    return sortIx.tolist()

## Recursive Bisection
This is where the actual asset weights are determined. The algorithm recursively divides the asset clusters into two equal subsets and rebalances the weights using an inverse-variance approach.

**Constraints:** The recursive bisection algorithm is flexible and can easily incorporate constraints on individual asset weights or groups of assets.

:::success
<font color=blue>Algorithm: Recursive Bisection</font>
Input:
- Covariance matrix $V$
- Quasi-diagonalized order of assets $L$

Output:
- Portfolio weights $w$

1. Initialize:
* $L = \{L_0\}$, where $L_0 = \{1, \dots, N\}$
* $w_n = 1$ for $n = 1, \dots, N$
2. While any $|L_i| > 1$ for $L_i \in L$:
* For each $L_i \in L$ with $|L_i| > 1$:
* Bisect $L_i$ into $L^{(1)}_i$ and $L^{(2)}_i$:
* $L^{(1)}_i$ is the first $\text{int}(\frac{|L_i|}{2})$ elements of $|L_i|$
* $L^{(2)}_i$ is the remaining elements of $L_i$
* Compute inverse-variance weights for each subset $j = 1, 2$:
* $V^{(j)}_i$ is the covariance matrix between the constituents of the $L^{(j)}_i$ bisection
* $\tilde{W}^{(j)}_i = \frac{\text{diag}[V^{(j)}_i ]^{−1}}{\text{tr}[\text{diag}[V^{(j)}_i]^{−1}]}$
* $\tilde{v}^{(j)}_i ≡ (\widetilde{W}^{(j)}_i)^T V^{(j)}_i \widetilde{W}^{(j)}_i$
* Compute split factor: $\alpha = 1 - \frac{\tilde{v}^{(1)}_i}{\tilde{v}^{(1)}_i + \tilde{v}^{(2)}_i}$ $\in [0,1]$
* Update weights:
* For $n \in L^{(1)}_i$: $w_n \leftarrow \alpha w_n$
* For $n \in L^{(2)}_i$: $w_n \leftarrow (1 - \alpha) w_n$
* Update $L$:
$L = L \backslash \{L_i\} \cup \{L^{(1)}_i, L^{(2)}_i\}$
3. Return $W$
:::

### Code Implementation

def getIVP(cov,**kargs):
	# Compute the inverse-variance portfolio
	ivp=1./np.diag(cov)
	ivp/=ivp.sum()
	return ivp

def getClusterVar(cov,cItems):
	# Compute variance per cluster
	cov_=cov.loc[cItems,cItems] # matrix slice
	w_=getIVP(cov_).reshape(-1,1)
	cVar=np.dot(np.dot(w_.T,cov_),w_)[0,0]
	return cVar

def getRecBipart(cov,sortIx):
    # Compute HRP alloc
    w=pd.Series(1,index=sortIx)
    cItems=[sortIx] # initialize all items in one cluster
    while len(cItems)>0:
        cItems=[i[j:k] for i in cItems for j,k in ((0,len(i)/2),\
            (len(i)/2,len(i))) if len(i)>1] # bi-section
        for i in xrange(0,len(cItems),2): # parse in pairs
            cItems0=cItems[i] # cluster 1
            cItems1=cItems[i+1] # cluster 2
            cVar0=getClusterVar(cov,cItems0)
            cVar1=getClusterVar(cov,cItems1)
            alpha=1-cVar0/(cVar0+cVar1)
            w[cItems0]*=alpha # weight 1
            w[cItems1]*=1-alpha # weight 2
    return w

### Covariance Predictors
We can see that only the method of predicting the covariance matrix directly affects the calculation of each asset weight.
Here we will introduce two difference covariance predictors: Rolling Window and Exponentially Weighted Moving Average (EWMA)
#### Rolling Window

def rolling_window(returns, memory, min_periods=20):
	rollw=pd.DataFrame()
	min_periods = max(min_periods, 1)
	times = returns.index
	assets = returns.columns
	returns = returns.values

	Sigmas = np.zeros((returns.shape[0], returns.shape[1], returns.shape[1]))
	Sigmas[0] = np.outer(returns[0], returns[0])

	for t in range(1, returns.shape[0]):
		alpha_old = 1 / min(t + 1, memory)
		alpha_new = 1 / min(t + 2, memory)

		if t >= memory:
			Sigmas[t] = alpha_new / alpha_old * Sigmas[t - 1] + alpha_new * (
             np.outer(returns[t], returns[t])
             - np.outer(returns[t - memory], returns[t - memory])
          )
		else:
			Sigmas[t] = alpha_new / alpha_old * Sigmas[t - 1] + alpha_new * (
                np.outer(returns[t], returns[t])
            )

	Sigmas = Sigmas[min_periods - 1 :]
	times = times[min_periods - 1 :]
	rollw=pd.concat([pd.DataFrame(Sigmas[t], index=assets, columns=assets) for t in range(len(times))], keys=[times[t].date() for t in range(len(times))])
	return rollw

#### Exponentially Weighted Moving Average (EWMA)

## Example: US Finance Stocks and Futures Portfolio

import matplotlib.pyplot as mpl
import scipy.cluster.hierarchy as sch
from scipy.spatial.distance import squareform
import numpy as np
import pandas as pd
import yfinance as yf  # Yahoo Finance API
from datetime import timedelta ,date
from dateutil.relativedelta import relativedelta

To showcase the effectiveness of HRP, we constructed a portfolio consisting of 35 diversified assets:

* **US Stocks:** `AAPL`, `ABT`, `ADBE`, `AMAT`, `AMZN`, `AVY`, `BALL`, `BAX`, `BDX`, `CMI`, `CPB`, `CSX`, `GILD`, `HAS`, `INSM`, `KO`, `MCD`, `MMM`, `MSFT`, `NVDA`, `PFE`, `TGT`, `TJX`, `TSM`, `WFC`, `XOM`, `YUM`
* **Futures Contracts:** `GC=F` (Gold), `HG=F` (Copper), `SI=F` (Silver), `^FVX` (CBOE Volatility Index), `^GSPC` (S&P 500), `^NDX` (Nasdaq 100), `^TNX` (10-Year Treasury Yield), `^TYX` (30-Year Treasury Yield)

The portfolio was rebalanced monthly using an 8-month lookback period, spanning from September 2000 to May 2024.

class Stocks():
  def __init__(self, arg_stocks_list=[], arg_begin='2015-01-01', arg_end='2023-10-01'):
    """
    Initialize the Stocks class.

    Parameters:
    - arg_stocks_list: List of stock symbols.
      For example:
        ['2330.TW','2337.TW','2357.TW','2454.TW','3231.TW','3443.TW'] or
        ['TSM','AAPL','GOOG','MSFT','AMZN','GOOGL']
    - arg_begin: Start date for the stock data.
    - arg_end: End date for the stock data.

    Returns:
    - None: This function doesn't return anything but initializes the class instance.
    """
    self.stocks_list = arg_stocks_list  # List of stock symbols
    self.stocks_list.sort()
    self.start_date = arg_begin  # Start date
    self.end_date = arg_end  # End date
    self.N = len(self.stocks_list)  # Number of stocks
    self.data = None
    self.set_data()
    self.M = len(self.data)  # The number of rows in the data
    self.rtns = pd.DataFrame()
    self.log_rtns = pd.DataFrame()
    self.set_rtns()
  def set_data(self):
    self.data = pd.DataFrame()  # Initialize an empty DataFrame to store data
    # Iterate over each stock in the list of stocks
    for stock in self.stocks_list:
      # Download the stock data using the yfinance library
      df = yf.download(stock, start=self.start_date, end=self.end_date)
      # Select only the adjusted closing price from the downloaded data
      df = df[['Adj Close']]
      # Rename the column to the stock symbol for clarity
      df = df.rename(columns = {'Adj Close': stock})
      # Check if the main DataFrame is empty
      if self.data.empty:
        # If it is, assign the downloaded data to the main DataFrame
        self.data = df
      else:
        # If it's not, concatenate the downloaded data to the main DataFrame
        self.data = pd.concat([self.data, df], axis=1)
    return
  def set_rtns(self):
    self.rtns = self.data.pct_change() + 1
    self.log_rtns = np.log(self.data.pct_change() + 1)
    self.rtns.fillna(1, inplace=True)
    self.log_rtns.fillna(1, inplace=True)
    return
  def plotComparePrice(self,use_log_rtns=True):
    mpl.figure(figsize=(20,12))
    mpl.title('Normalized Prices', fontsize=20)
    mpl.xticks(fontsize=16)
    mpl.yticks(fontsize=16)
    if use_log_rtns:
      mpl.ylabel('Logarithmic Price', fontsize=16)
      mpl.plot(1+ np.log(self.data/self.data.loc[self.data.index[0]]),label=self.data.columns)
    else:
      mpl.ylabel('Price', fontsize=16)
      mpl.plot((self.data/self.data.loc[self.data.index[0]]),label=self.data.columns)
    mpl.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=12)
    mpl.savefig('ComparePrice.png', format='png', bbox_inches='tight')
    mpl.clf();mpl.close()
    return

stocks = Stocks(stocks_list,start_date ,end_date)
stocks.plotComparePrice()


---

def testPerformance():
	# data
	start_date = date(2000,1,1)
	end_date = date(2024,6,1)
	td = (end_date - start_date)/timedelta(days=1)
	stocks_list=['AAPL', 'ABT', 'ADBE', 'AMAT', 'AMZN', 'AVY', 'BALL', 'BAX', 'BDX', 'CMI', 'CPB', 'CSX','GC=F', 'GILD', 'HAS','HG=F', 'INSM', 'KO', 'MCD', 'MMM', 'MSFT', 'NVDA', 'PFE','SI=F', 'TGT', 'TJX', 'TSM', 'WFC', 'XOM', 'YUM', '^FVX', '^GSPC', '^NDX', '^TNX', '^TYX']
	stocks = Stocks(stocks_list,start_date ,end_date)
	weight_history=pd.DataFrame()
	tm=8
	assetHRP=1
	assetEW=1
	assetRandom_1=1
	assetRandom_2=1
	assetRandom_3=1
	xl=[start_date + relativedelta(months=+(tm))]
	y_HRP=[1]
	y_EW=[1]
	y_Random_1=[1]
	y_Random_2=[1]
	y_Random_3=[1]
	np.random.seed(0)
	rollw=rolling_window(stocks.rtns, memory=168, min_periods=168)
	for i in range(int((td/30)-tm-4)):
		x=stocks.rtns[start_date + relativedelta(months=+i) : start_date + relativedelta(months=+(tm+i))]
		
		dt=start_date + relativedelta(months=+(tm+i))
		temp=rollw[dt:dt]
		while temp.empty:
			dt=dt + relativedelta(days=+1)
			temp=rollw[dt:dt]
		cov=pd.DataFrame(temp.values, index=temp.columns, columns=temp.columns)

		corr=x.corr()
		# cluster
		dist=correlDist(corr)
		condensed_dist = squareform(dist)
		link=sch.linkage(condensed_dist,'single')
		sortIx=getQuasiDiag(link)
		sortIx=corr.index[sortIx].tolist() # recover labels
		# Capital allocation
		hrp=getRecBipart(cov,sortIx)
		weight_history=pd.concat([weight_history,pd.DataFrame([hrp],index=[start_date + relativedelta(months=+(tm+i))])])
		tempHRP=hrp.copy()
		tempEW=0
		tempRandom_1=hrp.copy()
		rd_w_1=random_weight(len(stocks.rtns.columns))
		tempRandom_2=hrp.copy()
		rd_w_2=random_weight(len(stocks.rtns.columns))
		tempRandom_3=hrp.copy()
		rd_w_3=random_weight(len(stocks.rtns.columns))

		k=len(stocks.rtns.columns)-1
		for key in stocks_list:
			rtn=float(stocks.rtns[start_date + relativedelta(months=+(tm+i)): start_date + relativedelta(months=+(tm+1+i))].cumprod(axis = 0).iloc[-1][key])
			tempHRP[key]*=rtn
			tempRandom_1[key]=rd_w_1[k]*rtn
			tempRandom_2[key]=rd_w_2[k]*rtn
			tempRandom_3[key]=rd_w_3[k]*rtn
			k-=1
			tempEW+=rtn
		tempEW/=len(stocks_list)
		assetHRP*=tempHRP.sum()
		assetEW*=tempEW
		assetRandom_1*=tempRandom_1.sum()
		assetRandom_2*=tempRandom_2.sum()
		assetRandom_3*=tempRandom_3.sum()
		xl.append(start_date + relativedelta(months=+(tm+i)))
		y_HRP.append(assetHRP)
		y_EW.append(assetEW)
		y_Random_1.append(assetRandom_1)
		y_Random_2.append(assetRandom_2)
		y_Random_3.append(assetRandom_3)
	plotPerformance(xl,y_HRP.copy(),y_EW.copy(),y_Random_1.copy(),y_Random_2.copy(),y_Random_3.copy())
	plotHistoryWeight(weight_history)
	summary_stat(y_HRP.copy())
	summary_stat(y_EW.copy())
	summary_stat(y_Random_1.copy())
	summary_stat(y_Random_2.copy())
	summary_stat(y_Random_3.copy())
	plotDrawdown(xl,y_HRP.copy(),y_EW.copy(),y_Random_1.copy(),y_Random_2.copy(),y_Random_3.copy())
	return

def random_weight(arg_len):
	rd = np.random.random(size=arg_len-1)
	rd.sort()
	rd_w=rd.copy()
	rd_w=np.append(rd_w,1)
	for k in range(len(rd)):
		rd_w[k+1]=rd_w[k+1]-rd[k]
	return rd_w

### Tree Clustering
The dendrogram below illustrates the hierarchical clustering of assets, revealing natural groupings based on their correlation structure.

def plotDendrogram(path,linkage,labels=None):
	if labels is None:labels=[]
	p = len(labels)
	mpl.figure(figsize=(12,6))
	mpl.title('Hierarchical Clustering', fontsize=20)
	mpl.ylabel('stock symbol', fontsize=16)
	mpl.xlabel('distance', fontsize=16)

	# call dendrogram to get the returned dictionary 
	# (plotting parameters can be ignored at this point)
	R = sch.dendrogram(
		linkage,
		truncate_mode='lastp',  # show only the last p merged clusters
		p=p,  # show only the last p merged clusters
		no_plot=True, orientation='right'
	)
	# create a label dictionary
	temp = {R["leaves"][ii]: labels[ii] for ii in range(len(R["leaves"]))}
	def llf(xx):
		return "{}".format(temp[xx])

	sch.dendrogram(
		linkage,
		truncate_mode='lastp',  # show only the last p merged clusters
		p=p,  # show only the last p merged clusters
		color_threshold=0.35, # If color_threshold is None or ‘default’,the threshold is set to 0.7*max(Z[:,2])
		leaf_label_func=llf,
		# leaf_rotation=60.,
		leaf_font_size=12.,
		show_contracted=True,  # to get a distribution impression in truncated branches
		orientation='right', 
	)
	mpl.savefig(path, format='png', bbox_inches='tight')
	mpl.clf();mpl.close() # reset pylab
	return

Use `color_threshold=0.35` let the height less than 0.35 be same color, which also meaning they are highly correlated assets


### Quasi-Diagonalization
By reordering the correlation matrix based on the clustering, we obtain a quasi-diagonalized matrix. This visually highlights blocks of highly correlated assets, aiding in risk assessment.

def plotCorrMatrix(path,corr,labels=None):
	# Heatmap of the correlation matrix
	if labels is None:labels=[]
	mpl.figure(figsize=(8,6))
	mpl.pcolor(corr)
	mpl.colorbar()
	mpl.yticks(np.arange(.5,corr.shape[0]+.5),labels)
	mpl.xticks(np.arange(.5,corr.shape[0]+.5),labels,rotation = 60.)
	mpl.savefig(path)
	mpl.clf();mpl.close() # reset pylab
	return
plotCorrMatrix('HRP_corr.png',corr,labels=corr.columns)

**Correlation matrix**

sortIx=getQuasiDiag(link)
sortIx=corr.index[sortIx].tolist() # recover labels
df0=corr.loc[sortIx,sortIx] # reorder
plotCorrMatrix('HRP_corr_Q_diag.png',df0,labels=df0.columns)

**Quasi-diagonalisation of the correlation matrix**

### Recursive Bisection and Performance:
The final HRP portfolio weights were determined through recursive bisection.
Here we use Covariance Predictors-Rolling Window to Predict the covariance matrix of rebalance month.

**History Weight**

def plotHistoryWeight(w_history):
	tot = np.zeros(len(w_history))
	mpl.figure(figsize=(30,13))
	mpl.title('History of weight values', fontsize=30)
	mpl.xticks(fontsize=16)
	mpl.yticks(fontsize=16)
	mpl.ylabel('Weights', fontsize=24)
	for col_idx in w_history:
		mpl.fill_between(x=w_history.index.tolist(), y1=w_history[col_idx].values+tot,y2=tot,step='post',label=col_idx)
		tot += w_history[col_idx].values
	mpl.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=12)
	mpl.savefig('HistoryWeight.png', format='png', bbox_inches='tight')
	mpl.clf();mpl.close()
	return

**With Rolling Window**

**With EWMA**
![]()

**Compare Performance**
The resulting portfolio's performance is compared against equal-weighted (EW), Buy and hold portfolio (B&H) and randomly weighted portfolios in the figures below.

def plotPerformance(d,target1,target2,target3,target4,target5):
	mpl.figure(figsize=(12,6))
	mpl.plot(d, target1, 'r',label='HRP')
	mpl.plot(d, target2, 'b',label='EW')
	mpl.plot(d, target3, 'g',label='Random_1')
	mpl.plot(d, target4, 'm',label='Random_2')
	mpl.plot(d, target5, 'c',label='Random_3')
	mpl.title('Performance', fontsize=20)
	mpl.xticks(rotation = 60.)
	mpl.ylabel('cumprod', fontsize=16)
	mpl.legend()
	mpl.savefig('Performance.png', format='png', bbox_inches='tight')
	mpl.clf();mpl.close()
	return


**History Drawdown**

def historyDrawdown(target):
	history_drawdown=[]
	max_price=target[0]
	for price in target:
		if price > max_price:
			max_price=price
		history_drawdown.append((price - max_price) *100 / max_price)
	return history_drawdown

def plotDrawdown(d,target1,target2,target3,target4,target5):
	drawdown_HRP=historyDrawdown(target1)
	drawdown_EW=historyDrawdown(target2)
	drawdown_Random_1=historyDrawdown(target3)
	drawdown_Random_2=historyDrawdown(target4)
	drawdown_Random_3=historyDrawdown(target5)

	mpl.figure(figsize=(12,6))
	mpl.title('History Drawdown(%)', fontsize=16)
	mpl.plot(d, drawdown_HRP, 'r',label='HRP Drawdown(%)')
	mpl.plot(d, drawdown_EW, 'b',label='EW Drawdown(%)')
	mpl.plot(d, drawdown_Random_1, 'g',label='Random_1 Drawdown(%)')
	mpl.plot(d, drawdown_Random_2, 'm',label='Random_2 Drawdown(%)')
	mpl.plot(d, drawdown_Random_3, 'c',label='Random_3 Drawdown(%)')
	mpl.ylabel('Drawdown(%)', fontsize=16)
	mpl.legend()
	mpl.savefig('HistoryDrawdown.png', format='png', bbox_inches='tight')
	mpl.clf();mpl.close()
	return 


### Key Observations
* **Risk Management:** HRP demonstrates lower volatility and smaller maximum drawdowns, particularly during the 2008 financial crisis and the 2020 pandemic, suggesting superior risk-adjusted performance.
* **Competitive Returns:**
- HRP achieved a CAGR of 16.87%, slightly outperforming the equally weighted (EW) portfolio with a CAGR of 16.26%
- HRP’s volatility was 15.78%, slightly lower than the EW portfolio’s 15.88% , lower than all random portfolios.
- HRP’s Sharpe ratio was higher than EW and all random portfolios.
- HRP experienced a maximum drawdown of 36.14%, which was slightly better than the EW portfolio’s 37.18%, better than Random Portfolio 2 and 3.
* **Summary Statistics:** The table below summarizes key performance metrics:

def summary_stat(target):
	max_drawdown = 0
	max_price=target[0]
	for price in target:
		if price > max_price:
			max_price=price
		drawdown = (price - max_price) / max_price
		if drawdown < max_drawdown:
			max_drawdown = drawdown
	pct_rtn=(np.array(target[1:])/np.array(target[:-1]))-1
	print('Compound annual growth rate',(target[-1]/target[0])**(12/len(target))-1)
	print('Vol.(annual)', (pct_rtn.std())*(12**(1/2)))
	print('Sharpe ratio(annual)',(pct_rtn.mean())*(12**(1/2))/(pct_rtn.std()))
	print('max_drawdown',max_drawdown)
	print('Sortino ratio(annual)',(pct_rtn.mean())*(12**(1/2))/(pct_rtn[pct_rtn<0].std()))
	print('Calmar ratio(annual)',(pct_rtn.mean())*(12**(1/2))/max_drawdown)
	return

| | HRP with Rolling window | EW | Random 1 | Random 2 | Random 3 | 1/n B&H | HRP with EWMA |
|:----------------------:|:-----------------------:|:-------:|:--------:|:--------:|:--------:|:-------:|:-------------:|
| CAGR | 16.87% | 16.26% | 18.62% | 15.36% | 16.35% | | |
| Volatility(annualized) | 15.78% | 15.88% | 16.46% | 17.10% | 17.18% | | |
| Sharpe ratio | 1.077 | 1.038 | 1.131 | 0.929 | 0.976 | | |
| Maximum drawdown | -36.14% | -37.18% | -31.14% | -37.82% | -37.11% | | |

---

## **Reference**
* [Marcos López de Prado](https://www.quantresearch.org/index.html)
* [Building Diversified Portfolios that Outperform Out-of-Sample](https://doi.org/10.3905/jpm.2016.42.4.059)
* [HRP.py.txt](https://www.quantresearch.org/HRP.py.txt)
* [HRP_MC.py.txt](https://www.quantresearch.org/HRP_MC.py.txt)
* [Hierarchical Construction of Investment Portfolios Using Clustered Machine Learning](https://patents.google.com/patent/US20180089762A1)
* [Asset Allocation - Hierarchical Risk Parity](https://www.mathworks.com/videos/asset-allocation-hierarchical-risk-parity-1577957794663.html)
* [Covariance/Correlation Matrix HRP-Clustering](https://youtu.be/wtd3K-Ubr1g?si=rQomxWT2Xinw8UXo)
* [Towards Robust Portfolios](https://www.quantresearch.org/FIVE_VROBUST_Index__Research-EN.pdf)
* [Hierarchical Risk Parity on RAPIDS: An ML Approach to Portfolio Allocation](https://developer.nvidia.com/blog/hierarchical-risk-parity-on-rapids-an-ml-approach-to-portfolio-allocation/)
* [Asset Allocation - Hierarchical Risk Parity](https://youtu.be/e21MfMe5vtU?si=hS23MDNGDQlSmD-A)
* [Covariance Predictor](https://hackmd.io/@Howard531/B1DphEqXR)