In [1]:
import numpy as np
The model we are looking is is the standard one-(systemic)-factor model where a Gaussian vector $X_i$ is driven by a Standard Gaussian vector as follows $$ X_i = \rho * Z_0 + \sqrt{1-\rho^2} * Z_i $$
By definition the $Z_i$ satisfy $E[Z_iZ_j]=\delta_{ij}$ where $\delta_{ij}$ is the well known Kronecker Delta. Simple algebra then yields that the covariance of the $X_i$ is as follows $$ Cov[X_i X_j] = E[X_i X_j] = \delta_{ij} + (1-\delta_{ij}) \cdot 2 \rho^2 $$ which is a complicated way of writing that the $X_i$ are of unity variance, and that the covariance between different $X_i$ is $2\rho^2$.
Here we fill in the model parameters: the dimension of the space $d$, the correlation parameter $\rho$, and the number of draws $N$
In [2]:
d = 4
rho = 0.5
N = 1000
In order to implement this efficiently all the above equations must be encoded into a vector / matrix notation so that we can avoid loops in the intepreted Python code.
We start with the vector $z=(z_0, \ldots, z_d)$ where $d$ is the number of unique random variables $X_{i=1\ldots d}$. As we generate the whole series of samples in one go $z$ is actually a matrix $z_{in}$, with the index $i$ denoting the random variable, and $n=0...N-1$ denoting the number of the draw.
We will now generate this matrix where the rows z[n][.] correspond to one draw of the whole set, and where the columns z[.][i] correspond to all draws of the variable $Z_i$. In particular, z[.][0] are all draws of the systemic factor $Z_0$.
In [3]:
z = np.random.standard_normal((N, d+1))
#z
We now need a matrix $M$ that translates single draw $z$ into a single draw $x$. Note that this very same matrix will, when applied to the matrix containg all draws of $Z$ generate a matrix containing all draws of $X$.
As a reminder, we have $$ z = (z_0, z_1, \ldots, z_d) $$ and we want $$ x = (\rho z_0 + \bar{\rho} z_1, \ldots, \rho z_0 + \bar{\rho} z_d) $$ where $\bar{\rho}=\sqrt{1-\rho^2}$ is a simple shorthand.
The matrix that transforms $z$ into $x$ has $d+1$ rows and $d$ cols, with the first row being equal to $\rho$ and the first sub-diagonal being $\bar{\rho}$ (note that we have row-vectors, so they will be multiplied onto the matrix from the left). The matrix looks as follows $$ m = \begin{pmatrix} \rho & \rho & \ldots & \rho \\ \bar{\rho} \\ & \bar{\rho} \\ & & \ddots \\ & & & \bar{\rho} \end{pmatrix} $$ and we generate it by stacking a vector containing $\rho$ everywhere onto a matrix containing $\bar{\rho}$ on the diagonal
In [4]:
ma = rho * np.ones(d)
mb = sqrt(1-rho*rho) * np.identity(d)
m = np.vstack((ma,mb))
The $x$ are now generated by simple left multiplication $$ x = z \cdot m $$
$x$ is like $z$ a matrix where the rows x[n][.] correspond to one draw of the whole set, and where the columns x[.][i-1] correspond to all draws of the variable $X_i$. Note the shift in index because $i=1\ldots d$!
In [5]:
x = np.dot(z,m)
We are now testing whether we produce the correct covariance matrix. As a reminder: we want $1$ on the diagonal, the $2\rho^2$ on all off-diagonal elements. The term rowvar=0 indicates our dataformat, ie all observations of one draw in one row, and the bias=1 term means we are using the biased estimator that is normalised by $N$ instead of $N-1$
In [6]:
print "2 rho^2 = %f" % (rho*rho*2)
cov = np.cov(x, rowvar=0, bias=1)
cov
When $\rho$ is not constant our matrix looks as follows $$ m = \begin{pmatrix} \rho_1 & \rho_2 & \ldots & \rho_d \\ \bar{\rho_1} \\ & \bar{\rho_2} \\ & & \ddots \\ & & & \bar{\rho_d} \end{pmatrix} $$ which means we can no longer create it from building blocks but we have to run a short loop. Everything else remains the same.
In [7]:
rhovec = np.random.uniform(0,1,d) # random factor loadings, just to have some numbers to play
m2 = np.zeros((d+1,d))
for i in range(0, d):
rhoi = rhovec[i]
m2[0][i] = rhoi
m2[i+1][i] = sqrt(1-rhoi*rhoi)
rhovec, m2
Out[7]:
In [8]:
x2 = np.dot(z,m2)
In [9]:
cov = np.cov(x2, rowvar=0, bias=1)
cov
Out[9]:
(c) Stefan Loesch / oditorium 2014; all rights reserved (license)
In [10]:
import sys
print(sys.version)