Simulating Correlated Stock Returns in Python (SciPy)

Question

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?

Answer 1

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.