[QF] Probing Volatility I: Covariance Matrix & Eigen Decomposition

The Covariance Matrix

Exploring the mathematical applications when tackling volatility & correlational properties of securities in the market, as well as the market itself, by dissecting covariance matrices

From basic statistics, for \(n=1,2,...,N\) time indexes of \(r_{i_n}\) and \(r_{j_n}\) as returns data of asset i and j, their covariance \(\sigma_{ij}\) is calculated as:

• \(\sigma_{ij} = \frac{\sum_{1}^{N}(r_{i_n} - \bar{r_i})(r_{j_n} - \bar{r_j})}{N-1}\) for sample size
• \(\sigma_{ij} = \frac{\sum_{1}^{N}(r_{i_n} - \bar{r_i})(r_{j_n} - \bar{r_j})}{N}\) for population size

When it comes to investment as a whole, I find most of my work tackling almost every problems related with risk management, or volatility. Recall our previous brief overview on the essential topics of this research series, given a constructed returns table \(RET\), we obtain the de-meaned version: \(\overline{RET} = \begin{vmatrix} (r_{1_1} - \bar{r_{1}})&(r_{2_1} - \bar{r_{2}}) &\cdots &(r_{M_1} - \bar{r_{M}}) \\ (r_{1_2} - \bar{r_{1}}) &(r_{2_2} - \bar{r_{2}}) &\cdots &(r_{M_2} - \bar{r_{M}}) \\ \vdots&\vdots &\ddots &\vdots \\ (r_{1_N} - \bar{r_{1}})&(r_{2_N} - \bar{r_{2}}) &\cdots &(r_{M_N} - \bar{r_{M}}) \end{vmatrix}\)

A N x M matrix representing N interval (hourly/daily/weekly/etc) returns of M securities, we can construct a sample size Covariance Matrix \(\boldsymbol{C}\)

\[\boldsymbol{C} = \frac{\overline{RET}^\top \cdot \overline{RET}}{N-1} = \begin{vmatrix} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{1M} \\ \sigma_{21} & \sigma_2^2 & \cdots & \sigma_{2M}\\ \vdots &\vdots &\ddots &\vdots \\ \sigma_{M1} & \sigma_{M2} &\cdots &\sigma_M^2 \end{vmatrix}\]

where \(\boldsymbol{C}\) is a M x M square matrix such that: -\(\sigma_i^2\) = variances of individual securities \(i = 1,2,...,M\)

• \(\sigma_{ij}\) where \(i\neq j\) = covariances between two different assets i & j

Important points to highlight:

• Although we can find the standard deviation of individual securities \(i = 1,2,...,M\) as \(\sigma_i = \sqrt{\sigma_i^2}\) , there is no “standard deviation” of two securities, where the covariance is written as \(\sigma_{ij}\).
• The covariance value is an unbounded measurement value that describes a different perspective, a perspective of looking at the “movement relationship” between two different entities.
• The covariance value \(\sigma_{ij}\) encompasses/derives a relevant metric correlation \(\rho_{ij} = \frac{\sigma_{ij}}{\sigma_i \sigma_j}\), where such value, unlike covariance, is bounded between (-1,1).
• In other words, to which we will subsequently further explore this application, correlation is basically a standardized version of covariance - \(\boldsymbol{N(0,1)}\), such that the “standardized variances”, being the diagonal values of \(\boldsymbol{C}\) turn to 1.

---

For demonstration purposes moving forward, I have preprocessed & cleaned a returns table \(RET\), logarithmically calculated from daily close prices from Apil 2019 - April 2020 (encompassing the latest COVID-19 market crash), of 950-1000 securities on the U.S. Equities Market, including individual stocks in all industries/sectors as well as ETFs of domestic (SPY, XLK, XLV, etc…) & foreign securities (FXI, EWJ, etc… ). Data obtained with IEX Cloud API

#---| Initialize all modules needed |---#
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

#----| Load preprocess RET Table
RET = pd.read_csv("ret.csv",index_col=0,header=0)
RET

	AMG	UNH	EVR	AMD	NKE	NRG	EV	VRSN	SNPS	PKI	...	Z	PLNT	PEN	MSGS	PSTG	HPE	MTCH	SQ	TEAM	UA
date
4/1/2019	0.039181	-0.006981	0.026032	0.032385	0.012040	-0.016616	0.023294	0.023945	0.020798	0.012684	...	0.016274	0.014590	-0.023260	-0.000205	0.031175	0.030005	-0.006202	0.018514	0.013082	-0.001591
4/2/2019	-0.000180	-0.004613	-0.001178	0.014687	-0.010142	0.013787	0.006523	0.004935	-0.000510	-0.001846	...	0.024340	0.012116	0.005624	0.014667	0.013696	-0.005044	0.007438	-0.009479	0.015079	0.010034
4/3/2019	0.018240	0.005834	-0.002468	0.081451	0.001185	-0.000472	0.006959	0.006454	-0.000766	0.012041	...	0.027533	0.009448	-0.012190	0.004362	0.025559	0.001895	0.011228	0.018998	-0.002078	0.005764
4/4/2019	0.006065	0.006285	0.007814	0.002409	0.009544	-0.002365	0.010466	-0.010475	-0.018567	-0.007841	...	-0.001615	-0.007042	-0.016749	-0.011583	-0.026876	0.013785	-0.037143	-0.033387	-0.045579	0.032385
4/5/2019	0.012367	0.005603	0.001917	-0.003789	0.001406	-0.004033	0.006134	0.016571	0.011473	0.011991	...	0.003763	0.010125	0.004551	0.006717	0.004384	0.004966	-0.005629	0.006820	0.010020	-0.006597
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3/26/2020	0.101575	0.085379	0.053125	0.062323	0.064807	0.071551	0.073436	0.052101	0.053364	0.052086	...	0.022066	0.039924	0.079247	0.010383	0.109337	0.114583	-0.031277	0.067172	0.062938	0.061196
3/27/2020	-0.008481	-0.051996	-0.001504	-0.019558	-0.012774	0.045344	-0.022240	-0.011378	-0.028454	-0.054757	...	-0.072186	-0.038137	-0.023598	0.009727	-0.098243	-0.069807	-0.017452	-0.049201	-0.040977	-0.051534
3/30/2020	0.069726	0.035772	0.003863	0.027109	0.025504	-0.037936	0.039860	0.080825	0.038871	0.022604	...	-0.023468	-0.036764	0.049167	-0.048365	-0.023118	0.023152	0.064909	0.030647	-0.007692	0.010759
3/31/2020	-0.021578	-0.007590	-0.013585	-0.051007	-0.031409	-0.042728	-0.046351	-0.045906	-0.004029	-0.016207	...	-0.028464	0.002055	-0.011524	-0.064122	-0.008097	-0.034416	-0.009494	-0.048808	-0.028017	-0.042508
4/1/2020	-0.051881	-0.049568	-0.009599	-0.040840	-0.043348	-0.015527	-0.062040	-0.026046	-0.022219	-0.060932	...	-0.135993	-0.139947	-0.049035	-0.014964	-0.089231	-0.013479	-0.069436	-0.112428	-0.023664	-0.078700

We will use Pandas to perform calculation for our covariance matrix, as it is automatically normalized by (N-1) for sample population. We confirm this by performing our own calculation to check for sanity

def sample_covariance(ret):
    demeaned_RET = pd.DataFrame({
        s:(ret[s]-ret[s].mean()) for s in ret.columns})
    return ((demeaned_RET.T.dot(demeaned_RET))/(len(ret)-1))

raw_cov = sample_covariance(RET)
pd_cov = RET.cov()

False not in [np.allclose(raw_cov[i],pd_cov[i]) for i in raw_cov.columns]

True

They indeed do match. We proceed onto the next step.

Eigen Decomposition

This topic is a part of Linear Algebra, an extremely powerful spectrum within mathematics, that elevated my perspective from reality to abstraction, specifically in modeling nature using infinite-dimensional abstract vector space. For further academic readings for ground-up studies, I recommend watching this MIT lecture series.

From a linear algebraic perspective, specifically in our approach, a matrix acts as a vector space, defined by a system of equations, containing subspaces such as row space, column space, etc… as well as an operator that perform transformation, or mapping, of vectors from one subspace to another.

Example: Casting a shadow of an object = mapping a 3D entity to its 2D form

In this work, we are exploring a mapping from N-dimensional space to another N-dimensional space, specifically working with a square matrix. The main point of Eigen Decomposition is simply to find invariant subspaces under transformation by a square matrix. The analysis of such extracted eigen pairs is otherwise known as Principal Component Analysis (PCA), commonly used in other non-finance topics for dimensionality reduction, clustering data features, etc…

NOTE: Unlike the common usages such as compression of images’ pixels, or classification of different plant types, we are utilizing PCA to extract eigen components from a covariance matrix calculated from time series datasets of financial instruments. Always keep this context in mind when approaching mathematically, abstract or applied.

Given our M x M covariance matrix \(\boldsymbol{C}\) of M securities, as square matrix calculated from returns of N time indexes, we seek to find \(e_i\) and \(\lambda_i\) , respectively as the eigen vectors and eigen values, where for \(i = 1,2,...,M\), comes with such eigen pair that satisfy the condition: \(\boldsymbol{C} \cdot e_i = \lambda_i e_i\), such that:

• \(e_{i} = \begin{vmatrix} e_{i_1}\\ e_{i_2}\\ \vdots\\ e_{i_M} \end{vmatrix}\) , where we can construct an M x M matrix of the linearly independent eigen vectors as \(\boldsymbol{E} = \begin{vmatrix} e_1 &e_2 &\cdots &e_M \end{vmatrix} =\begin{vmatrix} e_{1_1} & e_{2_1} &\cdots &e_{M_1} \\ e_{1_2} & e_{2_2} &\cdots &e_{M_2} \\ \vdots &\vdots &\ddots &\vdots \\ e_{1_M} &e_{2_M} &\cdots &e_{M_M} \end{vmatrix}\)

• \(\lambda_i\) is a real value, each associated with \(e_i\), to which we can diagonalize it into another M x M square matrix: \(\boldsymbol{V} = \begin{vmatrix} \lambda_1 &0 &\cdots &0 \\ 0 &\lambda_2 &\cdots &0 \\ \vdots &\vdots &\ddots &\vdots \\ 0 &0 &\cdots &\lambda_M \end{vmatrix}\)

Due to the nature of \(\boldsymbol{C}\) being not just as a square matrix, but also a Hermitian Matrix, our values will be real & the matrix is diagonalizable such that: \(\boldsymbol{C} = \boldsymbol{E^\top} \cdot \boldsymbol{V} \cdot \boldsymbol{E} = \begin{vmatrix} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{1M} \\ \sigma_{21} & \sigma_2^2 & \cdots & \sigma_{2M}\\ \vdots &\vdots &\ddots &\vdots \\ \sigma_{M1} & \sigma_{M2} &\cdots &\sigma_M^2 \end{vmatrix}\)

To observe the significance of to why this is powerful when it comes to applying it on the covariance matrix, we will first highlight the significance of the eigen pairs (\(e_i\) & \(\lambda_i\))

• Eigen vectors \(e_i\) are normalized independent vectors from each other, meaning that they are orthogonal to each other with \(\left \| e_i \right \| = 1\)
• Since there are M eigen vectors \(e_i\) extracted from an M x M matrix, their orthogonality means that they represent the basis of M-dimensional vector space.
• Therefore, in other words, we perceive the eigen vectors as the directional vectors basis, the axis that builds the entire subspace that described by such matrix, or in this case the covariance matrix \(\boldsymbol{C}\)
• Each eigen value \(\lambda_i\) simply represents the directional variance of its associated eigen vector \(e_i\) , independent from one another, such that together they explain the total volatility of M securities, without any specified allocations.

---

We will first visualize the meaning of such extract eigen pairs by giving an example of two securities - SPY(S&P500 ETF) and IEF(7-10 Year Treasury ETF). Plotting their returns against each other in 2-Dimensional Space and we have:

ax = RET.plot.scatter("SPY","IEF",alpha=0.2)
ax.set_aspect("equal")

Observe the relationship between two assets, exhibiting negative correlation, reflecting behaviorially through the financial market (i.e.: As the market goes down, demand falls for equities and rises in treasury, especially during COVID-19 market crash ).

RET[["SPY","IEF"]].corr()

	SPY	IEF
SPY	1.000000	-0.516561
IEF	-0.516561	1.000000

If we plot the cummulative returns between the two assets, we can distinctly observe their correlational properties as exhibited through the scatter plot:

(1+RET[["SPY","IEF"]]).cumprod().plot()

Using sklearn PCA module to extract our eigen pairs. Below is our pre-coded function to perform PCA, working on dataframe being inputted & return the eigen vectors & eigen values:

#----| Perform Eigen Decomposition on Covariance Matrix of any given returns table, standardized/normalized or not
#----- dynamic for any NxM size of returns matrix
def EIGEN(ret,_simple=False):
    _pca = PCA().fit(ret)
    eVecs = _pca.components_
    eVals = _pca.explained_variance_
    _eigenNames = ["eigen_"+str(e+1) for e in range(len(eVals))]
    _eigenValues = pd.Series(eVals,index=_eigenNames,name="eigenValues")
    _eigenVectors = pd.DataFrame(eVecs.T,index=ret.columns,columns=_eigenNames)
    return _eigenVectors,_eigenValues

With each eigen vector magnified by 3 times the square root of its associated eigen value (since each eigen value represents the directional variance of such eigen vector, as stated above), graphing such vectors on top of the scatter plot, each in both directions, using matplotlib we obtain:

ax = plt.figure().gca()
ret = RET[["SPY","IEF"]]
x,y = list(ret.columns)
vectors,lambdas = EIGEN(ret)
ax.scatter(ret[x], ret[y], alpha=0.1)
ax.set_xlabel(x)
ax.set_ylabel(y)
_start = ret.mean(axis=0)
_colors = pd.Series(["red","green"],index=lambdas.index)
for e in lambdas.index:
    v = vectors[e] * (3*np.sqrt(lambdas[e]))
    ax.plot([_start[x]-v[x],_start[x]+v[x]],
              [_start[y]-v[y],_start[y]+v[y]],color=_colors[e],linewidth=2,label=e)
ax.set_aspect("equal")
ax.legend()

REMARK: The eigen components are orthogonal to each other, encapsulating the magnitudes of standard deviation in both dimensions. We interpret M eigen values as independently explained variances that, together, describe the total variance of our M selected securities.Such mathematical approach can take us further to any M-dimensional space such that we can describe the entire variance of M securities, in such abstract vector space that we can’t visualize, down to the independent components.

Plotly 3D Interactive Visualization

For the sake of illustration & fancy observation, below are functions written to actually plot such results for 3 Dimensional Space (M=3) using Plotly, a powerful interactive data visualization module:

import plotly.graph_objects as go
#----| Plot (with labels) & return result of 3 selected assets with each other with given table
def EIGEN_3D(ret):
    assert(ret.shape[1] == 3),"3 assets only"
    vectors,lambdas = EIGEN(ret)
    _start = ret.mean(axis=0)
    x,y,z = list(ret.columns)
    explained_variance = lambdas/lambdas.sum()
    _names = pd.Series([(i+" - "+str(
                    (explained_variance[i]*100).round(2))+"%") for i in explained_variance.index],
                           index=explained_variance.index)
    result = {"vectors":vectors,"lambdas":lambdas,
              "explained_variance":explained_variance}

    _DATA_ = []
    #----| SCATTER PLOT OF RETURNS DATA |----#
    _DATA_.append(go.Scatter3d(x=RET[x],y=RET[y],z=RET[z],mode="markers",marker=dict(
                                size=5,opacity=0.2),name=" - ".join([ret.index[0],ret.index[-1]])))

    #----| LINE PLOT OF EIGEN COMPONENTS |----#
    for e in lambdas.index:
        v = vectors[e] * (3*np.sqrt(lambdas[e]))
        _DATA_.append(go.Scatter3d(x=[_start[x]-v[x],_start[x]+v[x]],
                                   y=[_start[y]-v[y],_start[y]+v[y]],
                                   z=[_start[z]-v[z],_start[z]+v[z]],
                                   mode="lines",line=dict(width=6),name=_names[e]))

    fig = go.Figure(_DATA_)
    fig.update_layout(
        scene=dict(
        xaxis_title=x,
        yaxis_title=y,
        zaxis_title=z,aspectmode="auto"))
    result["figure"] = fig
    return result

Observe the results for SPY(S&P500 ETF), IEF(7-10 Year Treasury ETF) & GLD(SPDR Gold ETF), representing fundamentally different asset classes (Stocks, T-Bills & Gold respectively)

result = EIGEN_3D(RET[["SPY","IEF","GLD"]])
result["figure"].show()

Eigen Portfolios Introduction

Before we conclude this first section, I would like to shed light on some concepts that come up once observing this result:

• Given we just explore PCA application on finding the independent/orthogonal eigen vectors that when each magnified by its respective eigen value, describes the shape of the universe’s returns, with each eigen vector \(e_i\) being normalized into a unit vector, such that \(\left \| e_i \right \| = 1\)
• For any allocations weight of M amount of securities, whereas \(w = \begin{vmatrix} w_1\\ w_2\\ \vdots\\ w_M \end{vmatrix}\) such that \(-1 \leq w_{i}\leq 1\) and \(\sum_{1}^{M}w_i = 1\) we are simply finding a allocation vector in M-dimensional space!

We then pose the question:

What if we align our allocations vector in the direction of an eigen vector?

• By simply normalize any eigen vector \(e_{i} = \begin{vmatrix} e_{i_1}\\ e_{i_2}\\ \vdots\\ e_{i_M} \end{vmatrix}\) by its components summation, \(w_{e_i} = \frac{e_i}{e_i \cdot \vec{\boldsymbol{1}}} = \frac{e_i}{\sum_{n=1}^{M}e_{i_n}}\), we obtain an eigen portfolio

This is actually a field of study to which there exist many different research papers, feel free to look them up. Although, commonly agreed among the community, empirically proven, the largest eigen portfolio (eigen_1) associating with the largest eigen value is referred to as the “market portfolio” when M securities are picked as securities that comprises the market, which is explained in the next chapter. Though we have understood, as well as visualize, the meanings & functionalities of eigen decomposition on a symmetric M x M Covariance Matrix, derived from any given N x M returns table of M securities, we have only learned the ABC of such mathematical language. There arises many questions regarding the invariant nature of such eigen components, implementations & utilities on how we could use such knowledge to benefit our investment decision making. Moving forward to the subsequent chapters, we will start diving deeper into the nature of such eigen components, thus covariance & correlations, in a much heavier quantitative way, specifically with more details on portfolio & risk management and touching base on optimization as a whole.