Notebook Computing Optimal Submatrix of an Invalid Correlation Matrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing Optimal Submatrix of an Invalid Correlation Matrix\n",
    "\n",
    "\n",
    "### Goal: Given invalid correlation matrix C,find the optimal submatrix that has the largest similarity to the corresponding principal submatrix.\n",
    "\n",
    "### Key Mathematical Operations:\n",
    "\n",
    "* $C \\in R^{n \\times n}$ : fake correlation matrix\n",
    "* [$J(C)$](https://www.researchgate.net/publication/280050954_The_most_simple_methodology_to_create_a_valid_correlation_matrix_for_risk_management_and_option_pricing_purposes): function computing the optimal valid correlation matrix representing $C$ \n",
    "* $\\hat{C} = J(C),\\in R^{n×n}$\n",
    "* $similarity ∶= 1 - \\frac{2}{n(n-1)} \\times \\sum_{i<j}|\\frac{C_{ij}-\\hat{C}_{ij}}{C_{ij}}|$\n",
    "* $change(i) := \\sum_{j}|C_{ij}-\\hat{C}_{ij}|$\n",
    "* $i^* = max_{_{i=1,2,3,...,n}}change(i)$\n",
    "* $j^* = i^*$\n",
    "* $C_{sub\\setminus i^*,j^*} := \\{ C_{ij} \\}_{i\\neq i^*,j\\neq j^*},\\in R^{ n-1 \\times n-1}$\n",
    "* $\\hat{C}_{sub\\setminus i^*,j^*} := J(C_{sub\\setminus i^*,j^*}), \\in R^{ n-1 \\times n-1}$\n",
    "\n",
    "### Algorithm\n",
    "$C_{optimal}=C$ <br />\n",
    "$\\hat{C}_{optimal}=\\hat{C}$ <br />\n",
    "$similarity_{optimal}=similarity$ <br />\n",
    "$for \\quad m = n-1,n-2,...,5,4,3:$ <br />\n",
    "$\\quad \\quad C=C_{sub\\setminus i^*,j^*}$ <br />\n",
    "$\\quad \\quad \\hat{C}=\\hat{C}_{sub\\setminus i^*,j^*}$ <br />\n",
    "$\\quad \\quad if \\quad similarity > similarity_{optimal}:$ <br />\n",
    "$\\quad \\quad \\quad \\quad C_{optimal}=C$ <br />\n",
    "$\\quad \\quad \\quad \\quad \\hat{C}_{optimal}=\\hat{C}$ <br />\n",
    "$\\quad \\quad \\quad \\quad similarity_{optimal}=similarity$ <br />\n",
    "$\\quad \\quad \\quad \\quad if \\quad similarity_{optimal} = 1:$ <br />\n",
    "$\\quad \\quad \\quad \\quad \\quad \\quad stop$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Scripting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as LA\n",
    "from statsmodels.stats.moment_helpers import cov2corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $J(C): find\\_valid\\_corr$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_valid_corr(C,search_size=1000,phi_size=100000,max_attempt=500,step_unit=1e-17,disp=False):\n",
    "    '''\n",
    "    Compute the optimal valid correlation matrix representing matrix C\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    C: numpy array, symmetric square matrix whose diagonal values equal to 1 with any other element from -1 to 1\n",
    "    search_size: int, number of valid correlation matrices for C in the search space\n",
    "    phi_size: int, number of standard normal draws\n",
    "    max_attempt: int, maximum attempts to construct a positive definite seed matrix by slightly changing eigenvalues \n",
    "    step_unit: float, at each attempt, increase non-positive eigenvalues by one step_unit \n",
    "    disp: boolean, if true print information on seed matrix, eigenvalues, and upper matrix after Cholesky decomposition\n",
    "    \n",
    "    Return:\n",
    "    -------\n",
    "    optimal_C_hat: numpy array, optimal valid correlation matrix representing C\n",
    "    '''\n",
    "    if C.shape[0] != C.shape[1]:\n",
    "        raise ValueError('Matrix is not square matrix')\n",
    "    if np.sqrt(np.sum((C.T - C)**2)) > (C.shape[0]**2)*1e-5:\n",
    "        raise ValueError('Matrix is not symmetric')\n",
    "    \n",
    "    val_,vec = LA.eigh(C)\n",
    "    if (val_>=-1e-08).all(): # positive is defined as >= -1e-8 due to numerical precision\n",
    "        print('C is already a PSD matrix')\n",
    "        return C\n",
    "    \n",
    "    chol_done = False\n",
    "    attempt = 0\n",
    "    while (not chol_done) and attempt <= max_attempt:\n",
    "        try:\n",
    "            #print('{}th attempt'.format(attempt))\n",
    "            val = val_.copy()\n",
    "            val[val<0] = attempt*step_unit # set negative eigenvalues into 0\n",
    "                               # 0 here represented by 1e-16, this is to make psd matrix a little more positive\n",
    "                               # to allow Cholesky Decomposition to work\n",
    "            B_star = vec*(np.sqrt(val)) # take square root for every eigen values\n",
    "            B = B_star/np.linalg.norm(B_star,axis=1,keepdims=True) # normalize row vectors\n",
    "            seed_C = np.matmul(B,B.T)\n",
    "            \n",
    "            U = LA.cholesky(seed_C) # return upper matrix\n",
    "            chol_done = True\n",
    "        except:\n",
    "            if attempt == max_attempt:\n",
    "                raise RuntimeError('Failed to find Seed Matrix')\n",
    "            attempt += 1\n",
    "    if disp:\n",
    "        print('seed correlation matrix:')\n",
    "        print(seed_C)\n",
    "        print('eigen values of seed matrix:')\n",
    "        print(LA.eigh(seed_C)[0])\n",
    "        print('chol_done:',chol_done,'at {}th attempt'.format(attempt))\n",
    "        print('Upper Matrix:')\n",
    "        print(U)\n",
    "    \n",
    "    optimal_C_hat = None\n",
    "    optimal_error = np.inf \n",
    "    for _ in range(search_size):\n",
    "        phi = np.random.randn(phi_size,C.shape[0])\n",
    "        R = phi.dot(U)\n",
    "        C_hat = cov2corr(np.cov(R,rowvar=False))\n",
    "        error = np.sum((C_hat - C)**2)\n",
    "        if error < optimal_error and (LA.eigh(C_hat)[0]>=-1e-8).all(): # positive is defined as >= -1e-5 due to computing reason\n",
    "            optimal_C_hat = C_hat\n",
    "            optimal_error = error\n",
    "        else:\n",
    "            continue\n",
    "    print('element-wise error:',optimal_error)\n",
    "    print('distance error:',error_measure(C,optimal_C_hat,'d'))\n",
    "    print('similarity:',error_measure(C,optimal_C_hat,'s'))\n",
    "    return optimal_C_hat\n",
    "    #return optimal_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $ Similarity: error\\_measure$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_measure(C,C_hat,method='d'):\n",
    "    '''\n",
    "    Measure the difference/similarity between Matrix C and its simulated counterparty C_hat\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    C: numpy array, symmetric square matrix whose diagonal values equal to 1 with any other element from -1 to 1\n",
    "    C_hat: numpy array, symmetric square matrix whose diagonal values equal to 1 with any other element from -1 to 1\n",
    "    method: str, distance/d to compute the square root of sum of square of element-wise difference\n",
    "                 similarity/s to compute the 1 - sum of element-wise change of the upper matrix devided by number of elements \n",
    "    Return:\n",
    "    -------\n",
    "    float, distance or similarity\n",
    "    '''\n",
    "    if C.shape != C_hat.shape:\n",
    "        raise ValueError('Simulated Matrix does not have the same dimension the original Matrix')\n",
    "    if C.shape[0] != C.shape[1]:\n",
    "        raise ValueError('Matrices are not square matrix')\n",
    "    if np.sqrt(np.sum((C.T - C)**2)) > (C.shape[0]**2)*1e-6:\n",
    "        raise ValueError('Matrix C is not symmetric')\n",
    "    if np.sqrt(np.sum((C_hat.T - C_hat)**2)) > (C_hat.shape[0]**2)*1e-6:\n",
    "        raise ValueError('Matrix C_hat is not symmetric')\n",
    "    \n",
    "    if method in ('distance','d'):\n",
    "        return np.sqrt(np.sum((C - C_hat)**2))\n",
    "    elif method in ('similarity','s'):\n",
    "        n = C.shape[0]\n",
    "        return 1-np.sum(np.abs(np.triu((C - C_hat),k=1)/C))/(n*(n-1)/2)\n",
    "    else:\n",
    "        raise NameError('method can only be distance(d)/similarity(s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $ Algorithm: find\\_optimal\\_subcorr$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_optimal_subcorr(C,stopping_dim=3):\n",
    "    '''\n",
    "    Compute the optimal valid submatrix that has the largest similarity to the corresponding principal submatrix of a given matrix\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    C: numpy array, symmetric square matrix whose diagonal values equal to 1 with any other element from -1 to 1\n",
    "    stopping_dim: int, the lowest dimension of interest to search on\n",
    "    \n",
    "    Return:\n",
    "    -------\n",
    "    opt_C: numpy array, the principal submatrix \n",
    "    opt_C_hat: numpy array, optimal valid submatrix that has the largest similarity to the corresponding principal submatrix\n",
    "    '''\n",
    "    cur_C = C.copy()\n",
    "    print('Dimension: '+str(C.shape[0]))\n",
    "    cur_C_hat = find_valid_corr(cur_C)\n",
    "    opt_C = cur_C\n",
    "    opt_C_hat = cur_C_hat\n",
    "    opt_sim = error_measure(cur_C,cur_C_hat,'s')\n",
    "    for m in range(C.shape[0]-1,stopping_dim-1,-1):\n",
    "        print('Dimension: '+str(m))\n",
    "        change_vec = np.sum(np.abs(cur_C - cur_C_hat),axis=1)\n",
    "        ind_max = change_vec.argmax()\n",
    "        remain_ind = [i for i in range(cur_C.shape[0]) if i != ind_max]\n",
    "        cur_C = cur_C[np.ix_(remain_ind,remain_ind)]\n",
    "        cur_C_hat = find_valid_corr(cur_C)\n",
    "        sim = error_measure(cur_C,cur_C_hat,'s')\n",
    "        if sim > opt_sim:\n",
    "            opt_C = cur_C\n",
    "            opt_C_hat = cur_C_hat\n",
    "            opt_sim = sim\n",
    "            if opt_sim > 1-1e-10:\n",
    "                break\n",
    "    print('optimal similarity {} at dimension {}'.format(opt_sim,opt_C.shape[0]))\n",
    "    return opt_C, opt_C_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}