NumPy 2: arrays, broadcasting, vectorización, álgebra lineal y generadores aleatorios.

Del bucle al array

Hasta ahora, una serie de precios era una lista de Python y una rentabilidad, un bucle for (lección 0.3.4). Funciona, pero no escala: cinco años de cierres diarios son 1.300 números; cinco años de barras de un minuto de 500 acciones son casi 250 millones. NumPy resuelve el problema con una estructura, el ndarray, y una forma de pensar: operar con columnas enteras en lugar de elemento a elemento.

Un ndarray se diferencia de una lista en tres cosas:

Los cierres de SINT, el ETF ficticio del curso, se leen con pandas (lección 0.5.2) y se pasan a NumPy con to_numpy():

import numpy as np
import pandas as pd

precios = pd.read_csv("sint_diario.csv", parse_dates=["fecha"],
                      index_col="fecha")
cierres = precios["cierre"].to_numpy()
print(type(cierres).__name__, cierres.dtype, cierres.shape)
print(cierres[:5])
print(cierres[-1])
cierres[-1]

ndarray float64 (1303,)
[99.14 99.25 97.38 98.56 99.05]
142.22
np.float64(142.22)

La última línea muestra un detalle de NumPy 2: el repr de un escalar incluye su tipo (np.float64(142.22)), mientras que print muestra solo el número. Es un cambio de NumPy 2.0 (NEP 51) que sorprende al leer notebooks antiguos: el valor es el mismo; solo cambia cómo se representa.

Vectorización: rentabilidades sin bucles

La rentabilidad simple de un día es el cierre de hoy entre el de ayer, menos uno. Con un array, «hoy» son todos los cierres salvo el primero (cierres[1:]) y «ayer», todos salvo el último (cierres[:-1]). Dividir los dos arrays divide cada par de elementos alineados:

r_bucle = []
for i in range(1, len(cierres)):
    r_bucle.append(cierres[i] / cierres[i - 1] - 1)
r_bucle = np.array(r_bucle)

r = cierres[1:] / cierres[:-1] - 1
print(r.shape)
print(np.array_equal(r, r_bucle))

(1302,)
True

El resultado es idéntico bit a bit: las mismas divisiones en el mismo orden. La diferencia está en quién ejecuta el bucle. En la versión vectorizada, el bucle existe, pero lo recorre código compilado de NumPy sobre un bloque de memoria contiguo, sin crear un objeto de Python por cada número. Medido con timeit sobre un millón de precios, en el equipo en que se preparó esta lección la versión vectorizada fue unas 50 veces más rápida que el bucle; la cifra exacta depende de la máquina, el orden de magnitud no. El código para medirlo en tu equipo:

import timeit

largo = np.tile(cierres, 800)  # algo más de un millón de precios
timeit.timeit(lambda: largo[1:] / largo[:-1] - 1, number=10)

Las rentabilidades logarítmicas, $\ln(P_t / P_{t-1})$, se suman en el tiempo en lugar de multiplicarse. NumPy tiene funciones universales (ufuncs) como np.log, np.diff o np.cumprod que actúan sobre todo el array:

r_log = np.diff(np.log(cierres))
total_simple = np.prod(1 + r) - 1
total_log = np.expm1(r_log.sum())
print(f"Rentabilidad total: {total_simple:.4%}")
print(f"Con logaritmos:     {total_log:.4%}")
print(f"Último / primero:   {cierres[-1] / cierres[0] - 1:.4%}")

Rentabilidad total: 43.4537%
Con logaritmos:     43.4537%
Último / primero:   43.4537%

Las tres cuentas llegan al mismo número por caminos distintos. Comprobar un resultado por dos vías es un hábito que ahorra muchos errores en un backtest.

Las comparaciones también se vectorizan: r < -0.02 devuelve un array de bool del mismo tamaño, que sirve de máscara para seleccionar, contar o construir una señal con np.where:

caidas = r < -0.02
print(f"Sesiones con caída > 2 %: {caidas.sum()}")
print(f"Peor sesión: {r.min():.2%} (posición {r.argmin()})")
print(f"Media tras las caídas: {r[1:][caidas[:-1]].mean():.3%}")

# 1 si el cierre está por encima de su media de 200 sesiones
media_200 = np.convolve(cierres, np.ones(200) / 200, mode="valid")
senal = np.where(cierres[199:] > media_200, 1, 0)
print(senal.shape, senal.mean().round(3))

Sesiones con caída > 2 %: 45
Peor sesión: -4.81% (posición 338)
Media tras las caídas: -0.077%
(1104,) 0.745

True cuenta como 1 al sumar, así que caidas.sum() es el número de días que cumplen la condición. La tercera línea es un pequeño estudio de evento: la rentabilidad media del día siguiente a una caída fuerte; el desfase r[1:] frente a caidas[:-1] hace que solo se use la información disponible al cierre del día de la caída. Con 45 casos, un −0,077 % no permite concluir nada; el método sí sirve. La media móvil con np.convolve funciona, pero es la última vez que se calcula así en el curso: pandas lo hace en una línea y con las fechas alineadas (lección 0.5.3).

Broadcasting: operar con formas distintas

cartera_sint.csv tiene los cierres de cuatro activos ficticios: RVAR (renta variable), RFIJ (renta fija), ORO y TECH. En NumPy son una matriz con una fila por sesión y una columna por activo, y las rentabilidades se calculan igual que con un solo activo, porque el corte [1:] actúa sobre las filas:

cartera = pd.read_csv("cartera_sint.csv", parse_dates=["fecha"],
                      index_col="fecha")
activos = cartera.columns.to_list()
P = cartera.to_numpy()
R = P[1:] / P[:-1] - 1
medias = R.mean(axis=0)
print(activos, P.shape, R.shape, medias.shape)
print((medias * 252).round(4))

['RVAR', 'RFIJ', 'ORO', 'TECH'] (1303, 4) (1302, 4) (4,)
[0.1089 0.0053 0.0019 0.1232]

axis=0 recorre las filas y deja una media por columna. Para restar a cada activo su media, la matriz (1302, 4) y el vector (4,) tienen que operar juntos. NumPy lo permite con broadcasting: compara las formas empezando por la derecha y, en cada dimensión, acepta que sean iguales o que una de las dos valga 1 (o no exista, que cuenta como 1). La dimensión de tamaño 1 se «estira» sin copiar datos.

X = R - medias
print(X.shape)
print(np.abs(X.mean(axis=0)).max() < 1e-15)

Z = X / R.std(axis=0, ddof=1)   # rentabilidades estandarizadas
print(Z.std(axis=0, ddof=1).round(6))

(1302, 4)
True
[1. 1. 1. 1.]

Tras restar la media, cada columna tiene media cero salvo errores de redondeo minúsculos (por eso se compara con una tolerancia y no con == 0; la coma flotante se vio en la lección 0.3.1). Dividir por la desviación típica de cada columna estandariza las cuatro series: ya se pueden comparar en la misma escala.

La tercera fila de la tabla es el error típico. La rentabilidad de RVAR es un vector (1302,), y restarla a la matriz para obtener el rendimiento relativo de cada activo frente a la renta variable falla:

R - R[:, 0]

ValueError: operands could not be broadcast together with shapes (1302,4) (1302,)

La solución es dar a ese vector la forma de una columna, (1302, 1), con None (o np.newaxis) en el índice:

relativa = R - R[:, 0][:, None]
print(relativa.shape)
print((relativa.mean(axis=0) * 252).round(4))

(1302, 4)
[ 0.     -0.1036 -0.107   0.0143]

La primera columna es cero (RVAR menos sí misma). Cuando un broadcasting falla, imprime las formas (.shape) de los dos operandos y compáralas por la derecha: casi siempre falta o sobra una dimensión.

Álgebra lineal: covarianzas y volatilidad de una cartera

La matriz de covarianzas resume cómo se mueven los activos juntos. Con las rentabilidades ya centradas (X), es un producto de matrices: $\Sigma = X^\top X / (n - 1)$. El operador @ es el producto matricial de NumPy, y np.cov hace la misma cuenta:

n = R.shape[0]
S_manual = X.T @ X / (n - 1)
S = np.cov(R, rowvar=False)
print(np.allclose(S_manual, S))

S_anual = S * 252
volatilidades = np.sqrt(np.diag(S_anual))
correlaciones = np.corrcoef(R, rowvar=False)
for nombre, vol in zip(activos, volatilidades):
    print(f"{nombre:<5}{vol:7.2%}")
print(correlaciones.round(2))

True
RVAR  16.83%
RFIJ   4.96%
ORO   15.11%
TECH  28.27%
[[ 1.   -0.19  0.12  0.85]
 [-0.19  1.    0.19 -0.14]
 [ 0.12  0.19  1.    0.02]
 [ 0.85 -0.14  0.02  1.  ]]

rowvar=False le dice a np.cov que las variables (los activos) están en las columnas. La diagonal de la matriz son las varianzas de cada activo; anualizadas (×252 sesiones) y con raíz cuadrada, las volatilidades. La correlación de 0,85 entre TECH y RVAR y la negativa de RFIJ con los dos son las que hacen interesante combinar los activos.

La varianza de una cartera con pesos $w$ es la forma cuadrática $w^\top \Sigma w$. Con pesos iguales:

w = np.full(4, 0.25)
vol_cartera = np.sqrt(w @ S_anual @ w)
media_vols = w @ volatilidades
print(f"Media ponderada de volatilidades: {media_vols:.2%}")
print(f"Volatilidad de la cartera:        {vol_cartera:.2%}")

Media ponderada de volatilidades: 16.29%
Volatilidad de la cartera:        11.68%

La cartera es bastante menos volátil que la media de sus componentes (la figura de la cabecera compara las seis cifras): las correlaciones bajas o negativas compensan parte de los movimientos. Esa diferencia es la diversificación, y solo aparece si se usa la matriz completa; la media ponderada de volatilidades la ignora.

np.linalg resuelve sistemas de ecuaciones lineales, invierte matrices y calcula autovalores. Un ejemplo con consecuencia práctica: la cartera de mínima varianza (totalmente invertida) tiene pesos proporcionales a $\Sigma^{-1}\mathbf{1}$. En lugar de invertir la matriz, se resuelve el sistema $\Sigma x = \mathbf{1}$ con np.linalg.solve, que es más estable numéricamente:

x = np.linalg.solve(S_anual, np.ones(4))
w_mv = x / x.sum()
for nombre, peso in zip(activos, w_mv):
    print(f"{nombre:<5}{peso:7.1%}")
print(f"Volatilidad: {np.sqrt(w_mv @ S_anual @ w_mv):.2%}")

RVAR   17.7%
RFIJ   85.7%
ORO     1.0%
TECH   -4.4%
Volatilidad: 4.40%

El resultado concentra la cartera en renta fija y pone un peso negativo (una posición corta) en TECH: es lo que pide una optimización que solo mira la varianza. Por qué estas carteras son tan sensibles a la estimación de $\Sigma$ y cómo se acotan se ve en la lección 0.7.3; aquí basta con saber que @ y np.linalg hacen el trabajo.

Números aleatorios reproducibles

Simulaciones de Monte Carlo, bootstrap de rentabilidades, particiones de validación: el trading cuantitativo usa aleatoriedad en todas partes y necesita que sea reproducible. NumPy ofrece dos interfaces:

• La heredada: np.random.seed(42) seguido de np.random.normal(...). Fija la semilla de un generador global que comparte todo el programa, incluidas las librerías que importes.
• La actual: rng = np.random.default_rng(42) crea un Generator propio (algoritmo PCG64) que se pasa como argumento a quien lo necesite. La política de NumPy (NEP 19) recomienda esta para todo código nuevo; la heredada se mantiene solo por compatibilidad.

El problema del estado global se ve con dos funciones que comparten ese generador:

def simular_cartera():
    return np.random.normal(0.0004, 0.01, size=3).round(4)

def otra_libreria():
    np.random.normal(size=5)  # consume números del mismo generador

np.random.seed(42)
print(simular_cartera())
np.random.seed(42)
otra_libreria()
print(simular_cartera())

[ 0.0054 -0.001   0.0069]
[-0.0019  0.0162  0.0081]

Misma semilla, resultados distintos: entre la semilla y la simulación, otra parte del programa gastó números. En un proyecto real, ese «otra_libreria» puede ser una función de validación cruzada o un paquete que ni sabes que usa np.random. Con un Generator explícito, cada componente consume solo el suyo:

def simular_rentabilidades(rng, dias, mu=0.0004, sigma=0.01):
    return rng.normal(mu, sigma, size=dias)

a = simular_rentabilidades(np.random.default_rng(2026), 252)
b = simular_rentabilidades(np.random.default_rng(2026), 252)
print(np.array_equal(a, b))

# corrientes independientes para procesos en paralelo
hijos = np.random.default_rng(2026).spawn(3)
print([round(float(g.normal()), 4) for g in hijos])

True
[0.1107, 0.544, 1.403]

spawn deriva generadores hijos estadísticamente independientes a partir de uno padre: es la forma correcta de repartir aleatoriedad entre procesos (lección 0.4.6) sin que dos de ellos repitan la misma secuencia. float(...) pasa el escalar de NumPy a un float de Python para que la lista se imprima sin np.float64.

Un uso típico en trading es el bootstrap: remuestrear con reposición las rentabilidades observadas para ver cuánto podría variar un estadístico. Por ejemplo, la rentabilidad anualizada de SINT:

rng = np.random.default_rng(7)
idx = rng.integers(0, len(r), size=(5_000, len(r)))
medias_boot = r[idx].mean(axis=1) * 252
p5, p95 = np.percentile(medias_boot, [5, 95])
print(f"Observada: {r.mean() * 252:.1%}")
print(f"Intervalo del 90 %: de {p5:.1%} a {p95:.1%}")

Observada: 8.6%
Intervalo del 90 %: de -4.1% a 21.7%

idx es una matriz de 5.000 × 1.302 índices aleatorios; r[idx] usa esa matriz como índice (fancy indexing) y devuelve 5.000 muestras de golpe, sin un solo bucle. El intervalo es muy ancho: cinco años de datos dicen poco de la rentabilidad esperada, como anticipaba la lección 0.1.2.

Valores ausentes: NaN

Los datos reales tienen huecos: una sesión sin cotización, un fallo del proveedor, un activo que todavía no existía. En NumPy, un hueco en un array de float es np.nan (not a number), y tiene dos propiedades peligrosas:

con_hueco = cierres[:6].copy()
con_hueco[3] = np.nan          # sesión sin cotización
r_hueco = con_hueco[1:] / con_hueco[:-1] - 1
print(r_hueco.round(4))
print(r_hueco.mean(), np.nanmean(r_hueco).round(5))
print(np.isnan(r_hueco).sum(), "rentabilidades perdidas")
print(np.nan == np.nan, r_hueco > 0)

[ 0.0011 -0.0188     nan     nan -0.0039]
nan -0.00722
2 rentabilidades perdidas
False [ True False False False False]

Primero, se propaga: un solo cierre ausente deja dos rentabilidades en nan y la media de todo el array en nan. np.nanmean, np.nansum o np.nanstd ignoran los ausentes, pero hay que pedirlo. Segundo, toda comparación con nan es falsa, incluso nan == nan; por eso se detecta con np.isnan, no con ==. Para una señal, eso significa que r > 0 devuelve False en un día sin dato, y la regla lo trata como «no comprar» sin avisar.

Ejercicios

1. Calcula la volatilidad anualizada de cada uno de los cuatro activos sin usar np.cov: solo con R, broadcasting y np.sqrt. Comprueba que coincide con volatilidades.

varianzas = ((R - R.mean(axis=0)) ** 2).sum(axis=0) / (len(R) - 1)
vol_manual = np.sqrt(varianzas * 252)
print(vol_manual.round(4))
print(np.allclose(vol_manual, volatilidades))

[0.1683 0.0496 0.1511 0.2827]
True

2. Con pesos del 60 % en RVAR y del 40 % en RFIJ (y cero en los demás), calcula la volatilidad de la cartera y compárala con la media ponderada de las dos volatilidades. ¿Por qué la diferencia es mayor que en la cartera equiponderada?

w_6040 = np.array([0.6, 0.4, 0.0, 0.0])
print(f"Cartera:        {np.sqrt(w_6040 @ S_anual @ w_6040):.2%}")
print(f"Media ponderada: {w_6040 @ volatilidades:.2%}")

Cartera:        9.92%
Media ponderada: 12.08%

3. Simula con default_rng(11) 1.000 trayectorias de un año (252 sesiones) de rentabilidades normales con la media y la desviación típica diarias de SINT. ¿En qué porcentaje de trayectorias el año cierra en pérdidas?

rng = np.random.default_rng(11)
anos = rng.normal(r.mean(), r.std(ddof=1), size=(1_000, 252))
final = np.prod(1 + anos, axis=1) - 1
print(f"Años en pérdidas: {(final < 0).mean():.1%}")
print(f"Mediana: {np.median(final):.1%}")

Años en pérdidas: 36.7%
Mediana: 6.7%