Simulación de rendimientos bursátiles correlacionados en Python (SciPy)

Question

Estoy buscando generar rendimientos de acciones con correlación entre acciones en Python. Sin embargo, la salida no se está comportando correctamente y puede tener correlación temporal accidental causando problemas.

Este código está diseñado para generar número_de_caminos de los rendimientos correlacionados de las acciones dada una Serie de Panda de devuelve , un DataFrame de constante covarianzas y un DateIndex de fechas ( índice_de_fecha ).

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

Mi código de prueba basado en la siguiente matriz de covarianza muestra algunas rarezas en los resultados que podrían estar relacionadas.

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

A lo largo de un gran número de simulaciones con igual pesos_constituyentes (arriba) el index_returns parecen estar negativamente autocorrelacionados (usando autocorr en Pandas). El rendimientos_anuales del índice son inferiores al 7% esperado cuando la correlación es mayor que 0, pero iguales al 7% cuando la correlación es 0. Muy extraño.

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

¿Estoy usando mal multivariante_normal ¿de alguna manera?

Answer 1

No he ejecutado su código, pero este es un ejemplo donde la desigualdad de Jensen se mantiene. Supongamos que usted mira los rendimientos de un solo activo durante dos períodos, $R_1$ y $R_2$ . Los rendimientos aquí significan rendimientos brutos, como 1,07 en su caso. Si esos rendimientos son independientes, tenemos $\mathrm{E}(R_1R_2)$ es igual a $\mathrm{E}(R_1)\mathrm{E}(R_2)$ . Sin embargo, este producto no es lo que se mira. La rentabilidad se anualiza tomando el $n$ -enésima root, que es una función cóncava de la rentabilidad total. Definir $f$ para ser su función de anualización. Entonces, por la desigualdad de Jensen, $\mathrm{E}(f(R_1R_2)) \le f(\mathrm{E}(R_1R_2))$ . Si se asume una distribución particular para los rendimientos, se puede cuantificar esta brecha (por ejemplo, para los precios distribuidos lognormalmente, la media geométrica será la media aritmética menos la mitad de la varianza).

Un ejemplo numérico (con código R):

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

Se ve que el rendimiento total ( E[R1R2] ) coincide con la expectativa (1,07 por 1,07 es 1,1449), pero la rentabilidad media anualizada es inferior a 0,07.

---

Actualización: La desigualdad de Jensen no es estricta (e incluso si lo fuera, la brecha podría seguir siendo trivial). La brecha dependerá de la distribución de los rendimientos. Su ejemplo debería estar cerca del caso lognormal (y se acercaría más con más períodos). Así que la brecha es la media de un solo período menos la mitad de la varianza. Con 500 activos y sin correlación, la varianza caerá casi a cero, por lo que la media aritmética y la geométrica deberían ser casi iguales. A medida que aumente la correlación, también lo hará la varianza y, por tanto, la estimación de la media geométrica (es decir, la rentabilidad anualizada) disminuirá.

Aquí hay un ejemplo de código con 12 periodos y niveles crecientes de correlación [prefiero quedarme con el primer ejemplo, porque realmente son ejemplos sencillos :-) ]. La función randomReturns se toma del NMOF paquete, que yo mantengo. Devuelve una matriz de tamaño 12 por 500 de variantes normalmente distribuidas.

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

Después de 10000 muestras, el rendimiento total es más o menos el esperado ( 1.07^12 - 1 = 1.252 ), para todos los niveles de correlación. El rendimiento medio anualizado se aproxima razonablemente a lo que cabría esperar según la lognormalidad.