hrp_hbdc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6eca0855",
   "metadata": {},
   "source": [
    "# Hierarchical Risk Parity - A HBDC Research Note"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a024862",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Risk-Parity (RP) allocates proportional to inverse variance: $w_k = \\frac{\\nu_k^{-1}}{\\sum_k\\nu_k^{-1}}$. This completely ignores correlation information.  \n",
    "\n",
    "Mean-Variance-Optimization takes into account full covariance matrix. $w_k = \\Sigma^{-1}\\mu$. However, MVO weights are unstable (big changes for small changes in covariances) and frequently extreme (corner solutions).\n",
    "\n",
    "Hierarchical Risk-Parity sits in between these two extremes. It takes into account dependency information in somewhat \"coarse\" way by grouping assets that are similar in the sense of a distance metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "071e7d98",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TaskLocalRNG()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Distributions, Random, Statistics, LinearAlgebra\n",
    "Random.seed!(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab9a68ee",
   "metadata": {},
   "source": [
    "## An Example\n",
    "\n",
    "HRP is best understood through a numerical example. Consider three assets with the following covariance matrix of returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "id": "8e9cf276",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 0.0225      0.00900343  0.00946224\n",
       " 0.00900343  0.04        0.0137452\n",
       " 0.00946224  0.0137452   0.0225"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Σ = [0.0225 0.00900343 0.00946224; 0.00900343 0.04 0.0137452; 0.00946224 0.0137452 0.0225]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bdb73f0",
   "metadata": {},
   "source": [
    "We need a distance metric since neither covariance nor correlation is a proper distance metric. It can be shown that $d(i,j) = \\sqrt{(1-\\rho_{ij})/2}$ is a proper distance metric. We convert the covariance matrix into a distance matrix by first extracting the correlation matrix and then converting each correlation entry into a distance measure with the formula above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "d4591f83",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cov2corr (generic function with 1 method)"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function cov2corr(C)\n",
    "    σinv = diagm(1.0 ./ sqrt.(diag(C)))\n",
    "    return σinv * C * σinv\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "id": "bf830aa2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 1.0       0.300114  0.420544\n",
       " 0.300114  1.0       0.458173\n",
       " 0.420544  0.458173  1.0"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = cov2corr(Σ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "d017bc1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " 0.0       0.59156   0.538264\n",
       " 0.59156   0.0       0.520493\n",
       " 0.538264  0.520493  0.0"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = sqrt.(0.5*(1.0.-R))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f5ad6ed",
   "metadata": {},
   "source": [
    "Under risk-parity the allocations to three assets would only depend on the diagonal entries in $\\Sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "id": "92d01acf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ivp (generic function with 3 methods)"
      ]
     },
     "execution_count": 188,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ivp(C::Matrix{Float64})\n",
    "    iv = 1.0 ./ diag(C)\n",
    "    return iv ./ sum(iv)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "id": "8f110c03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 0.3902439024390244\n",
       " 0.21951219512195122\n",
       " 0.3902439024390244"
      ]
     },
     "execution_count": 189,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wrp = ivp(Σ)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3109b47",
   "metadata": {},
   "source": [
    "For this portfolio, under hierarchical risk-parity, we first find the pair that is closest in distance and treat it as one asset. From the distance matrix D, we see that assets 2 and 3 are the closest. We group them together into a portfolio (say, asset 4). The mix of 2 and 3 in this portfolio is done as per risk-parity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "id": "3c48453b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.36\n",
       " 0.64"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w2, w3 = ivp(Σ[2:3, 2:3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed40804e",
   "metadata": {},
   "source": [
    "Variance of this asset 4 can be found from the covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "id": "c3c7880d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "clustervar (generic function with 1 method)"
      ]
     },
     "execution_count": 191,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function clustervar(C, cluster)\n",
    "    V = C[cluster, cluster]\n",
    "    u = ivp(V)\n",
    "    return dot(u' * V, u)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "1cf962a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02073378816"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v4 = clustervar(Σ, [2,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8185571",
   "metadata": {},
   "source": [
    "Now we are left with two assets 1 and 4. We conclude the allocation process by simply allocating based on risk parity between assets 1 and 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "17df1c8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.47957370941607536\n",
       " 0.5204262905839246"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w1, w4 = ivp([Σ[1,1], v4])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9098c66",
   "metadata": {},
   "source": [
    "But, we know that asset 4 is just a combination of 36% of asset 2 and 64% of asset 3. So we cascade the allocation to 4 down to 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "id": "97d10662",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.18735346461021288\n",
       " 0.3330728259737118"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w2, w3 = w4 * [w2, w3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bf2c910",
   "metadata": {},
   "source": [
    "This gives the final HRP allocation of"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "id": "92d199ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Float64}:\n",
       " 0.47957370941607536\n",
       " 0.18735346461021288\n",
       " 0.3330728259737118"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "whrp = [w1, w2, w3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12aeda18",
   "metadata": {},
   "source": [
    "Compare this with RP. Asset 1 is dissimilar to both assets 2 and 3, so it receives a relatively higher allocation in HRP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "id": "5db3441d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Matrix{Float64}:\n",
       " 0.390244  0.479574\n",
       " 0.219512  0.187353\n",
       " 0.390244  0.333073"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[wrp whrp]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b9e0317",
   "metadata": {},
   "source": [
    "The recursive allocation can be visualised on a tree. Varaince of the portfolio represented by each node is in brackets. The final allocations are 48%, 52%x36%, 52%x64% for assets 1, 2 and 3 respectively.\n",
    "![hrp](hrp_tree.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb893b74",
   "metadata": {},
   "source": [
    "## Summary of HRP methodology\n",
    "\n",
    "1. Compute agglomorative clustering of assets in the given portfolio. That is, at every stage cluster or group pair of assets or clusters from the previous stage that are closest in the distance metric. This gives a binary tree. For example, in a portfolio of 4 assets numbered 1 to 4, we could have ((3,2),1,4) in the first stage, then (((3,2),1),4) in the next stage.\n",
    "\n",
    "2. Allocate progressively using RP from the top of the tree. In the above example, allocate between ((3,2),1) and 4 using RP. Then, allocate between (3,2) and 1 using RP. And finally, allocate between 3 and 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d897479",
   "metadata": {},
   "source": [
    "## Alternative Methodology\n",
    "\n",
    "It is possible to do clustering and recursive allocation in one pass. The idea is to start at the root of the tree, with all assets in one portfolio. At every stage we split the portfolio into two by removing the asset that is farthest from all other assets in the distance metric. We do this by computing for each asset its minimum distance to other assets. Then we find the asset that has the maximum of the minimum distances and remove it from the portfolio. We allocate using RP between this asset and the remaining portfolio. This process is repeated till there is only one asset remaining in the portfolio. This essentailly recreates agglomerative clustering in the reverse direction.\n",
    "\n",
    "Code below follows this method."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "445bdd86",
   "metadata": {},
   "source": [
    "## Helper Funcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "a0d1bd17",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ivp (generic function with 3 methods)"
      ]
     },
     "execution_count": 199,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ivp(v::Vector{Float64})\n",
    "    iv = 1.0 ./ v\n",
    "    return iv ./ sum(iv)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bcf4f36",
   "metadata": {},
   "source": [
    "Split a given vector of indices (id) into an index (k) that is farthest from the rest of the indices (id)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "0cba518c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "split (generic function with 1 method)"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function split(D::Matrix{Float64}, id::Union{UnitRange{Int64}, Vector{Int64}})\n",
    "    if length(id) == 1\n",
    "        return nothing, id\n",
    "    end\n",
    "    tmp = []\n",
    "    for i in id\n",
    "        push!(tmp, minimum([D[i,k] for k in setdiff(id, i)]))\n",
    "    end\n",
    "    k = id[argmax(tmp)]\n",
    "    id = setdiff(id, k)\n",
    "    return id, k\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ed47057",
   "metadata": {},
   "source": [
    "Compute HRP weights by recursive RP allocation between splits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "id": "9aa6bdef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hrp (generic function with 2 methods)"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function hrp(Σ, D)\n",
    "    id = 1:size(D,2)\n",
    "    w = fill(1.0, size(D,2))\n",
    "    while length(id) > 1\n",
    "        id, k = split(D, id)\n",
    "        v1 = clustervar(Σ, id)\n",
    "        v2 = Σ[k,k]\n",
    "        α = ivp([v1,v2])[1]\n",
    "        w[id] .= w[id] * α\n",
    "        w[k] = w[k] * (1. - α) \n",
    "    end\n",
    "    return w\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da090dd5",
   "metadata": {},
   "source": [
    "## A Bigger Example\n",
    "Generate fake data for a portfolio of 10 assets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "id": "1cf62105",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10×10 Matrix{Float64}:\n",
       " 0.0       0.55444   0.596279  0.871755  …  0.383599  0.851367  0.484048\n",
       " 0.55444   0.0       0.526299  0.850921     0.583789  0.629766  0.673496\n",
       " 0.596279  0.526299  0.0       0.739134     0.551224  0.662068  0.636274\n",
       " 0.871755  0.850921  0.739134  0.0          0.89054   0.537583  0.768226\n",
       " 0.764139  0.81803   0.651709  0.610314     0.75588   0.686476  0.759294\n",
       " 0.682054  0.806884  0.673949  0.819175  …  0.521068  0.766023  0.628537\n",
       " 0.896789  0.820001  0.71908   0.562226     0.783188  0.442969  0.838506\n",
       " 0.383599  0.583789  0.551224  0.89054      0.0       0.824285  0.598282\n",
       " 0.851367  0.629766  0.662068  0.537583     0.824285  0.0       0.73743\n",
       " 0.484048  0.673496  0.636274  0.768226     0.598282  0.73743   0.0"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = rand(LKJ(10,1), 1)[1]\n",
    "σ = rand(0.1:0.01:0.2, 10)\n",
    "Σ = diagm(σ)*R*diagm(σ) \n",
    "D = sqrt.((1.0 .- R)*0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "id": "885d88af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10×2 Matrix{Float64}:\n",
       " 0.0740929  0.0508197\n",
       " 0.0664989  0.0909951\n",
       " 0.12248    0.13186\n",
       " 0.0830661  0.128253\n",
       " 0.142048   0.183589\n",
       " 0.12248    0.103041\n",
       " 0.0740929  0.0607706\n",
       " 0.106694   0.0731804\n",
       " 0.0664989  0.0772713\n",
       " 0.142048   0.10022"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wrp = ivp(Σ)\n",
    "whrp = hrp(Σ, D)\n",
    "[wrp whrp]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.1",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}