NumPy 2: arrays, broadcasting, vectorization, linear algebra and random generators.

From the loop to the array

Until now, a price series was a Python list and a return was a for loop (lesson 0.3.4). It works, but it does not scale: five years of daily closes are 1,300 numbers; five years of one-minute bars for 500 stocks are almost 250 million. NumPy solves the problem with a structure, the ndarray, and a way of thinking: operating on whole columns instead of element by element.

An ndarray differs from a list in three ways:

The closes of SINT, the course's fictitious ETF, are read with pandas (lesson 0.5.2) and converted to NumPy with 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)

The last line shows a NumPy 2 detail: the repr of a scalar includes its type (np.float64(142.22)), while print shows only the number. It is a NumPy 2.0 change (NEP 51) that surprises people reading old notebooks: the value is the same; only its representation changes.

Vectorization: returns without loops

The simple return of a day is today's close divided by yesterday's, minus one. With an array, "today" is all the closes except the first (cierres[1:]) and "yesterday", all except the last (cierres[:-1]). Dividing the two arrays divides each pair of aligned elements:

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

The result is identical bit for bit: the same divisions in the same order. The difference lies in who runs the loop. In the vectorized version, the loop exists, but it is executed by NumPy's compiled code over a contiguous block of memory, without creating a Python object for each number. Measured with timeit over a million prices, on the machine where this lesson was prepared the vectorized version was about 50 times faster than the loop; the exact figure depends on the machine, the order of magnitude does not. The code to measure it on your machine:

import timeit

largo = np.tile(cierres, 800)  # just over a million prices
timeit.timeit(lambda: largo[1:] / largo[:-1] - 1, number=10)

Log returns, $\ln(P_t / P_{t-1})$, add up over time instead of multiplying. NumPy has universal functions (ufuncs) such as np.log, np.diff or np.cumprod that act on the whole array:

r_log = np.diff(np.log(cierres))
total_simple = np.prod(1 + r) - 1
total_log = np.expm1(r_log.sum())
print(f"Total return:   {total_simple:.4%}")
print(f"With logs:      {total_log:.4%}")
print(f"Last / first:   {cierres[-1] / cierres[0] - 1:.4%}")

Total return:   43.4537%
With logs:      43.4537%
Last / first:   43.4537%

The three calculations reach the same number by different paths. Checking a result two ways is a habit that saves many errors in a backtest.

Comparisons are vectorized too: r < -0.02 returns a bool array of the same size, which serves as a mask to select, count or build a signal with np.where:

caidas = r < -0.02
print(f"Sessions with a drop > 2%: {caidas.sum()}")
print(f"Worst session: {r.min():.2%} (position {r.argmin()})")
print(f"Mean after the drops: {r[1:][caidas[:-1]].mean():.3%}")

# 1 if the close is above its 200-session average
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))

Sessions with a drop > 2%: 45
Worst session: -4.81% (position 338)
Mean after the drops: -0.077%
(1104,) 0.745

True counts as 1 when summing, so caidas.sum() is the number of days that meet the condition. The third line is a small event study: the mean return on the day after a sharp drop; the offset of r[1:] against caidas[:-1] ensures that only the information available at the close of the drop day is used. With 45 cases, a −0.077% allows no conclusion; the method is still useful. The moving average with np.convolve works, but it is the last time the course computes it this way: pandas does it in one line and with the dates aligned (lesson 0.5.3).

Broadcasting: operating on different shapes

cartera_sint.csv has the closes of four fictitious assets: RVAR (equities), RFIJ (fixed income), ORO (gold) and TECH. In NumPy they are a matrix with one row per session and one column per asset, and the returns are computed the same way as with a single asset, because the [1:] slice acts on the rows:

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 runs along the rows and leaves one mean per column. To subtract each asset's mean, the (1302, 4) matrix and the (4,) vector have to operate together. NumPy allows it with broadcasting: it compares the shapes starting from the right and, in each dimension, accepts that they are equal or that one of the two is 1 (or does not exist, which counts as 1). The dimension of size 1 is "stretched" without copying data.

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

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

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

After subtracting the mean, each column has zero mean except for tiny rounding errors (which is why the comparison uses a tolerance and not == 0; floating point was covered in lesson 0.3.1). Dividing by each column's standard deviation standardizes the four series: they can now be compared on the same scale.

The third row of the table is the common mistake. The return of RVAR is a (1302,) vector, and subtracting it from the matrix to get each asset's performance relative to equities fails:

R - R[:, 0]

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

The solution is to give that vector the shape of a column, (1302, 1), with None (or np.newaxis) in the index:

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]

The first column is zero (RVAR minus itself). When broadcasting fails, print the shapes (.shape) of the two operands and compare them from the right: almost always a dimension is missing or extra.

Linear algebra: covariances and portfolio volatility

The covariance matrix summarizes how the assets move together. With the returns already centered (X), it is a matrix product: $\Sigma = X^\top X / (n - 1)$. The @ operator is NumPy's matrix product, and np.cov does the same calculation:

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 tells np.cov that the variables (the assets) are in the columns. The diagonal of the matrix holds each asset's variance; annualized (×252 sessions) and square-rooted, the volatilities. The correlation of 0.85 between TECH and RVAR and the negative one of RFIJ with both are what make combining the assets interesting.

The variance of a portfolio with weights $w$ is the quadratic form $w^\top \Sigma w$. With equal weights:

w = np.full(4, 0.25)
vol_cartera = np.sqrt(w @ S_anual @ w)
media_vols = w @ volatilidades
print(f"Weighted average of volatilities: {media_vols:.2%}")
print(f"Portfolio volatility:             {vol_cartera:.2%}")

Weighted average of volatilities: 16.29%
Portfolio volatility:             11.68%

The portfolio is considerably less volatile than the average of its components (the figure at the top compares the six numbers): low or negative correlations offset part of the moves. That difference is diversification, and it only shows up if the full matrix is used; the weighted average of volatilities ignores it.

np.linalg solves systems of linear equations, inverts matrices and computes eigenvalues. An example with practical consequences: the (fully invested) minimum-variance portfolio has weights proportional to $\Sigma^{-1}\mathbf{1}$. Instead of inverting the matrix, you solve the system $\Sigma x = \mathbf{1}$ with np.linalg.solve, which is numerically more stable:

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"Volatility: {np.sqrt(w_mv @ S_anual @ w_mv):.2%}")

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

The result concentrates the portfolio in fixed income and puts a negative weight (a short position) on TECH: it is what an optimization that only looks at variance asks for. Why these portfolios are so sensitive to the estimate of $\Sigma$ and how they are constrained is covered in lesson 0.7.3; here it is enough to know that @ and np.linalg do the work.

Reproducible random numbers

Monte Carlo simulations, return bootstraps, validation splits: quantitative trading uses randomness everywhere and needs it to be reproducible. NumPy offers two interfaces:

• The legacy one: np.random.seed(42) followed by np.random.normal(...). It sets the seed of a global generator shared by the whole program, including the libraries you import.
• The current one: rng = np.random.default_rng(42) creates a Generator of your own (PCG64 algorithm) that is passed as an argument to whoever needs it. NumPy's policy (NEP 19) recommends it for all new code; the legacy one is kept only for compatibility.

The problem with global state shows up with two functions that share that generator:

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

def otra_libreria():
    np.random.normal(size=5)  # consumes numbers from the same generator

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]

Same seed, different results: between the seed and the simulation, another part of the program used up numbers. In a real project, that otra_libreria can be a cross-validation function or a package you do not even know uses np.random. With an explicit Generator, each component consumes only its own:

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

# independent streams for parallel processes
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 derives statistically independent child generators from a parent one: it is the correct way to distribute randomness across processes (lesson 0.4.6) without two of them repeating the same sequence. float(...) converts the NumPy scalar to a Python float so that the list prints without np.float64.

A typical use in trading is the bootstrap: resampling the observed returns with replacement to see how much a statistic could vary. For example, the annualized return of 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"Observed: {r.mean() * 252:.1%}")
print(f"90% interval: from {p5:.1%} to {p95:.1%}")

Observed: 8.6%
90% interval: from -4.1% to 21.7%

idx is a 5,000 × 1,302 matrix of random indices; r[idx] uses that matrix as an index (fancy indexing) and returns 5,000 samples at once, without a single loop. The interval is very wide: five years of data say little about the expected return, as lesson 0.1.2 anticipated.

Missing values: NaN

Real data has gaps: a session with no quote, a provider failure, an asset that did not exist yet. In NumPy, a gap in a float array is np.nan (not a number), and it has two dangerous properties:

con_hueco = cierres[:6].copy()
con_hueco[3] = np.nan          # session with no quote
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(), "returns lost")
print(np.nan == np.nan, r_hueco > 0)

[ 0.0011 -0.0188     nan     nan -0.0039]
nan -0.00722
2 returns lost
False [ True False False False False]

First, it propagates: a single missing close leaves two returns as nan and the mean of the whole array as nan. np.nanmean, np.nansum or np.nanstd ignore the missing values, but you have to ask for it. Second, every comparison with nan is false, even nan == nan; that is why it is detected with np.isnan, not with ==. For a signal, this means that r > 0 returns False on a day with no data, and the rule treats it as "do not buy" without warning.

Exercises

1. Compute the annualized volatility of each of the four assets without using np.cov: only with R, broadcasting and np.sqrt. Check that it matches 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. With weights of 60% in RVAR and 40% in RFIJ (and zero in the others), compute the portfolio volatility and compare it with the weighted average of the two volatilities. Why is the difference larger than in the equal-weight portfolio?

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

Portfolio:        9.92%
Weighted average: 12.08%

3. With default_rng(11), simulate 1,000 one-year paths (252 sessions) of normal returns with SINT's daily mean and standard deviation. In what percentage of the paths does the year close at a loss?

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"Losing years: {(final < 0).mean():.1%}")
print(f"Median: {np.median(final):.1%}")

Losing years: 36.7%
Median: 6.7%