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



Implementation

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



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)



Testing

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




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



Appendix

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




Out[7]:

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




Out[9]:

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



Licence and version



In [10]:

import sys
print(sys.version)




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