I am trying to simulate correlation in order to price a correlation swap (via Monte-Carlo). For simplicity, let's assume we have 2 assets, and everything is correlated with $\rho$, and there is no drift or anything. Say, the swap has $M$ observation periods, and we simulate $N$ paths per asset.
So we start by using a multivariate normal distribution with $\mu=\begin{bmatrix}0\\0\end{bmatrix}$ and $\Sigma=\begin{bmatrix}1&\rho\\\rho&1\end{bmatrix}$. Now we generate $M$ correlated pairs of normal variables $R = \begin{bmatrix}r_{11}&\ldots & r_{1M}\\r_{21}&\ldots&r_{2M}\end{bmatrix} \in \mathbb{R}^{2\times M}$. We do this $N$ times, such that we obtain $R_1,\ldots,R_N$. For each sample $R_i$ we now compute the realised sample correlation $\bar{\rho}_i$ between its 2 assets.
Finally, we can compute the average of all realised sample correlations: $$\bar{\rho} = \frac1N \cdot\sum_{i=1}^N \bar{\rho}_i$$
What I have noticed is, the average realised sample correlation is disproportionally often lower than the target correlation, i.e. $$\bar{\rho} < \rho$$ I would have expected $\bar{\rho} \approx \rho$.
Below are a few sample plots, which highlight this effect (with 3 out of 4 paths lower than 50%, and 1 path around 50%):
Source Code to generate graph: zerobin
Question: Does anybody know why this is the case? Or am I generating random variables incorrectly?
