# Simulating Correlated Stock Returns in Python (SciPy)

I'm looking to generate stock returns with inter-stock correlation in Python. However, the output is not behaving properly and may have accidental temporal correlation causing issues.

This code is designed to generate num_paths of correlated stock returns given a Panda's Series of annualized returns, a DataFrame of constant covariances and a DateIndex of dates (date_index).

from pandas import DataFrame, concat
from scipy.stats import multivariate_normal

def correlated_returns(num_paths, returns, covariances, date_index, periods_per_year=1):
period_returns = (1 + returns) ** (1 / periods_per_year) - 1 if periods_per_year != 1 else returns
mn = multivariate_normal(period_returns, covariances / periods_per_year, allow_singular=True)

digits = len(str(num_paths))
paths = [DataFrame(mn.rvs(size=len(date_index)), index=date_index, columns=returns.index) for _ in range(num_paths)]
keys = [f'Run {str(run_num).zfill(digits)}' for run_num in range(num_paths)]
return concat(paths, axis='columns', keys=keys, names=['Run', 'Returns'])


My test code based on the following covariance matrix shows a few oddities in the results which might be related.

correlation = 0.2  # inter-stock correlation 0.18
annualized_return = 7 / 100  # Simulated return for each stock
stocks = [f'Stock {i}' for i in range(simulated_stocks)]
constituent_weights = DataFrame(1 / simulated_stocks, date_index, stocks)

returns = Series(annualized_return, stocks)
volatilities = Series(volatility, stocks)
correlations = DataFrame(correlation, stocks, stocks)
fill_diagonal(correlations.values, 1)
covariances = correlations.mul(volatilities, axis='index').mul(volatilities, axis='columns')


Over a large number of simulations using equal constituent_weights (above) the index_returns appear to be negatively auto-correlated (using autocorr in Pandas). The annualized_returns of the index come in less than 7% expected when the correlation is greater than 0, but equal to 7% when the correlation is 0. Very wierd.

for simulation in range(simulations):
runs = correlated_returns(1, returns, covariances, date_index, frequency_scale)
return_history = runs['Run 0']

index_returns = return_history.mul(constituent_weights).sum(axis='columns') \
.div(constituent_weights.sum(axis='columns'))

annualized_return = (1 + index_returns[1:]).prod() ** (1 / simulation_years) - 1


Am I misusing multivariate_normal somehow?

• Quick thanks to the a couple people on this site that helped me realize that I must have a code error rather than discovering new finance :) Nov 14, 2021 at 15:12

I have not run your code, but this is an example where Jensen's inequality holds. Suppose you look at returns of a single asset over two periods, $$R_1$$ and $$R_2$$. Returns here mean gross returns, such as 1.07 in your case. If those returns are independent, we have $$\mathrm{E}(R_1R_2)$$ equals $$\mathrm{E}(R_1)\mathrm{E}(R_2)$$. However, this product is not what you look at. You annualize return by taking the $$n$$-th root, which is a concave function of total return. Define $$f$$ to be your annualization function. Then, by Jensen's inequality, $$\mathrm{E}(f(R_1R_2)) \le f(\mathrm{E}(R_1R_2))$$. If you assume a particular distribution for the returns, you may quantify this gap (e.g. for lognormally distributed prices, the geometric mean will be the arithmetic mean minus half the variance).

A numerical example (with R code):

ans <- NULL
for (n in c(1e3, 1e4, 1e5, 1e6, 1e7)) {

R1 <- pmax(0, rnorm(n,mean = 1.07, sd = 0.2))
R2 <- pmax(0, rnorm(n,mean = 1.07, sd = 0.2))

ans <- rbind(ans,
c(n,
mean(R1), mean(R2), mean(R1*R2),
mean((R1*R2)^(1/2))))

}
colnames(ans) <- c("trials", "E[R1]", "[ER2]", "E[R1R2]", "E[f(R1R2)]")
ans
##        trials    E[R1]    [ER2]  E[R1R2] E[f(R1R2)]
## [1,]     1000 1.065672 1.076422 1.148000   1.061579
## [2,]    10000 1.066393 1.074114 1.145864   1.060501
## [3,]   100000 1.069893 1.069176 1.143758   1.059819
## [4,]  1000000 1.070105 1.069802 1.144723   1.060283
## [5,] 10000000 1.069981 1.070046 1.144939   1.060348


You see that the total return (E[R1R2]) matches the expectation (1.07 times 1.07 is 1.1449), but the average annualized return is lower than 0.07.

Update: Jensen's inequality isn't strict (and even if it were, the gap could still be trivial). The gap will depend on the distribution of returns. Your example should be close to the lognormal case (and it would get closer with more periods). So the gap is single-period mean minus half the variance. With 500 assets and no correlation, variance will drop to almost zero, so arithmetic and geometric mean should be almost the same. As correlation increases, so will variance, and hence the estimate of the geometric mean (i.e. the annualized return) will decrease.

Here is a code example with 12 periods and increasing levels of correlation [I'd rather keep the first example, because I really simple examples :-) ]. The function randomReturns is taken from the NMOF package, which I maintain. It returns a matrix of size 12 times 500 of normally-distributed variates.

library("NMOF")
rhos <- c(0, 0.1, 0.3, 0.9)  ## levels of correlation to test
na <- 500                    ## number of assets
n <- 1e4                     ## number of trials
w <- rep(1/na, na)           ## equal weights

ans <- NULL
for (rho in rhos) {
message(rho)

results.geom <- numeric(n)
results.TR <- numeric(n)
for (i in seq_len(n)) {
R <- randomReturns(na = na,
ns = 12,
mean = 0.07,
sd = 0.2,
rho = rho)
results.geom[i] <- prod(R %*% w + 1)^(1/12) - 1
results.TR[i]   <- prod(R %*% w + 1) - 1
}

## compute E: expected annualized
##            return under lognormality
C <- array(rho, dim = c(na, na))
diag(C) <- 1
E <- 0.07 - 0.5 * sum(0.2^2 * C)/na^2

ans <- rbind(ans,
c(rho, mean(results.TR), mean(results.geom), E))

}
colnames(ans) <- c("rho", "E[ΣR]", "E[f(ΣR)]", "ln-expected")
ans
##      rho  E[R1R2] E[f(R1R2)] ln-expected
## [1,] 0.0 1.252762 0.06998791    0.069960
## [2,] 0.1 1.258257 0.06843583    0.067964
## [3,] 0.3 1.253640 0.06471616    0.063972
## [4,] 0.9 1.244869 0.05353727    0.051996


After 10000 samples, the total return is roughly as expected (1.07^12 - 1 = 1.252), for all levels of correlation. The average annualized return is reasonably close to what would be expected under lognormality.

• The numbers here are the right order of magnitude to match my issue for sure, but I'm not sure quite how to understand the finance of what you are saying. Am I calculating the annualized return incorrectly at the end (I should be taking the mean of the annual results rather than the mean of the annualized return)? I'm still confused then why would the inter-stock correlation matter to the final result? Nov 16, 2021 at 21:40
• Specifically, if I understand the above shouldn't I be getting a mean annualized return of less than 7% even in the case where the inter-stock correlation is zero? Nov 16, 2021 at 22:02
• It has nothing to do with finance: it's all statistics :-) People are sometimes surprised that there is no unbiased estimator for standard deviation, unless you make assumptions about and adjustments for the distribution of the data. You could take the average total return (which will be unbiased) and then annualize it. That is, compute f(E(R)), not E(f(R)). I have updated the answer. Nov 17, 2021 at 8:44
• "With 500 assets and no correlation, variance will drop to almost zero, so arithmetic and geometric mean should be almost the same." Many thanks. This is what I was missing. The correlation changes the index volatility which changes the annual returns. Nov 17, 2021 at 11:17