import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import pandas_datareader.data as web
import datetime
import backtrader as bt
%matplotlib inline
import pyfolio as pf
import quantstats
import sys
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.stats import rankdata
from scipy.stats import stats
from scipy.optimize import minimize

start = '2018-06-18'
end = '2021-05-21'

vti = web.DataReader("VTI", 'yahoo', start, end)['Adj Close'].to_frame("vti")
vea = web.DataReader("VEA", 'yahoo', start, end)['Adj Close'].to_frame("vea")
vwo = web.DataReader("VWO", 'yahoo', start, end)['Adj Close'].to_frame("vwo")
tlt = web.DataReader("TLT", 'yahoo', start, end)['Adj Close'].to_frame("tlt")
ltpz = web.DataReader("LTPZ", 'yahoo', start, end)['Adj Close'].to_frame("ltpz")
vclt = web.DataReader("VCLT", 'yahoo', start, end)['Adj Close'].to_frame("vclt")
iau = web.DataReader("IAU", 'yahoo', start, end)['Adj Close'].to_frame("iau")
dbc = web.DataReader("DBC", 'yahoo', start, end)['Adj Close'].to_frame("dbc")
xlb = web.DataReader("XLB", 'yahoo', start, end)['Adj Close'].to_frame("xlb")
xlc = web.DataReader("XLC", 'yahoo', start, end)['Adj Close'].to_frame("xlc")
xle = web.DataReader("XLE", 'yahoo', start, end)['Adj Close'].to_frame("xle")
xlf = web.DataReader("XLF", 'yahoo', start, end)['Adj Close'].to_frame("xlf")
xli = web.DataReader("XLI", 'yahoo', start, end)['Adj Close'].to_frame("xli")
xlk = web.DataReader("XLK", 'yahoo', start, end)['Adj Close'].to_frame("xlk")
xlp = web.DataReader("XLP", 'yahoo', start, end)['Adj Close'].to_frame("xlp")
xlu = web.DataReader("XLU", 'yahoo', start, end)['Adj Close'].to_frame("xlu")
xlv = web.DataReader("XLV", 'yahoo', start, end)['Adj Close'].to_frame("xlv")
xly = web.DataReader("XLY", 'yahoo', start, end)['Adj Close'].to_frame("xly")
xlre = web.DataReader("XLRE", 'yahoo', start, end)['Adj Close'].to_frame("xlre")

price_df = pd.concat([vti, vea, vwo, tlt, ltpz, vclt, iau, dbc, xlb, xlc, xle, xlf, xli, xlk, xlp, xlu, xlv, xly, xlre], axis=1)
return_df = price_df.pct_change().dropna(axis=0)

return_df

# correlation matrix
corr = return_df.corr()

# distance matrix
d_corr = np.sqrt(0.5*(1-corr))

link = linkage(d_corr, 'single')
Z = pd.DataFrame(link)

fig = plt.figure(figsize=(25, 10))
dn = dendrogram(Z)
# plt.show()

Z

def get_quasi_diag(link):
    
    # sort clustered items by distance
    
    link = link.astype(int)
    
    # get the first and the second item of the last tuple
    sort_ix = pd.Series([link[-1,0], link[-1,1]]) 
    
    # the total num of items is the third item of the last list
    num_items = link[-1,3]
    
    # if the max of sort_ix is bigger than or equal to the max_items
    while sort_ix.max() >= num_items:

        sort_ix.index = range(0, sort_ix.shape[0]*2, 2) # make space
        
        df0 = sort_ix[sort_ix >= num_items] # find clusters
        
        # df0 contains even index and cluster index
        i = df0.index
        j = df0.values - num_items 
        
        sort_ix[i] = link[j,0] # item 1
        
        df0  = pd.Series(link[j, 1], index=i+1)
        
        sort_ix = sort_ix.append(df0) # item 2
        sort_ix = sort_ix.sort_index()
        
        sort_ix.index = range(sort_ix.shape[0])
    
    return sort_ix.tolist()

sort_ix = get_quasi_diag(link)
sort_ix

def get_cluster_var(cov, c_items):
    cov_ = cov.iloc[c_items, c_items] # matrix slice
    # calculate the inverse-variance portfolio
    ivp = 1./np.diag(cov_)
    ivp/=ivp.sum()
    w_ = ivp.reshape(-1,1)
    c_var = np.dot(np.dot(w_.T, cov_), w_)[0,0]
    return c_var

def get_rec_bipart(cov, sort_ix):
    # compute HRP allocation
    # intialize weights of 1
    w = pd.Series(1, index=sort_ix)
    
    # intialize all items in one cluster
    c_items = [sort_ix]
    while len(c_items) > 0:
        # bisection
        c_items = [i[int(j):int(k)] for i in c_items for j,k in 
                   ((0,len(i)/2),(len(i)/2,len(i))) if len(i)>1]
        
        # parse in pairs
        for i in range(0, len(c_items), 2):
            
            c_items0 = c_items[i] # cluster 1
            c_items1 = c_items[i+1] # cluter 2
            
            c_var0 = get_cluster_var(cov, c_items0)
            c_var1 = get_cluster_var(cov, c_items1)
            
            alpha = 1 - c_var0/(c_var0+c_var1)
            
            w[c_items0] *= alpha
            w[c_items1] *=1-alpha
            
    return w

cov = return_df.cov()
weights = get_rec_bipart(cov, sort_ix)

new_index = [return_df.columns[i] for i in weights.index]
weights.index = new_index

weights

weights[0:19]

plt.figure(figsize = (12, 8))
plt.bar(list(weights.index), weights)