iPython Cookbook - Monte Carlo I

Generating a Monte Carlo simulation based on a one (systemic) factor model

In [1]:
import numpy as np

The Model

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


Generating the $Z_i$

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))

Generating the matrix $m$

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))

Generating the $X_i$

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)

  File "<ipython-input-6-f7eb3b16b3f0>", line 1
    print "2 rho^2 = %f" % (rho*rho*2)
SyntaxError: invalid syntax


Generating the matrix $m$ for non-constant $\rho$

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

(array([ 0.06650426,  0.32545442,  0.17522132,  0.18735068]),
 array([[ 0.06650426,  0.32545442,  0.17522132,  0.18735068],
        [ 0.99778614,  0.        ,  0.        ,  0.        ],
        [ 0.        ,  0.94555773,  0.        ,  0.        ],
        [ 0.        ,  0.        ,  0.98452907,  0.        ],
        [ 0.        ,  0.        ,  0.        ,  0.98229309]]))

In [8]:
x2 = np.dot(z,m2)

In [9]:
cov = np.cov(x2, rowvar=0, bias=1)

array([[  1.03617347e+00,   8.42507497e-04,  -4.61986301e-02,
       [  8.42507497e-04,   9.73906504e-01,   1.71766657e-02,
       [ -4.61986301e-02,   1.71766657e-02,   9.78382916e-01,
       [ -3.99929059e-02,   3.35121828e-02,   5.88462278e-02,

Licence and version

(c) Stefan Loesch / oditorium 2014; all rights reserved (license)

In [10]:
import sys

3.4.0 (default, Apr 11 2014, 13:05:11) 
[GCC 4.8.2]