---
title: "0.5.1 NumPy 2: arrays, vectorization and linear algebra"
description: "NumPy 2 arrays for prices and returns: vectorization, broadcasting, a covariance matrix with np.linalg, seeded generators and NaN."
url: "https://pythonparatrading.com/en/python-trading-course/level-0/numpy-arrays-and-vectorization/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-0/numpy-arrays-and-vectorization.png"
---
# 0.5.1 · NumPy 2: arrays, broadcasting, vectorization, linear algebra and random generators
> The whole Python quant stack —pandas, Polars, scikit-learn, statsmodels— rests on NumPy arrays. This lesson teaches you to think in columns of numbers instead of loops, with the SINT prices and a four-asset portfolio.
- Course: [Level 0: Python from scratch and a data mindset](https://pythonparatrading.com/en/python-trading-course/level-0/)
- Module: 0.5 · Numerical computing and data
- Duration: 50 min
- Prerequisites: [0.3.5 · Functions: arguments, *args and **kwargs, lambdas and scope](https://pythonparatrading.com/en/python-trading-course/level-0/python-functions-arguments-scope/)
## By the end of this lesson
- Turn a series of closes into an `ndarray` and compute simple and log returns without loops.
- Apply the broadcasting rules to subtract means, normalize and compare the assets of a return matrix.
- Compute a portfolio's covariance matrix and its volatility with matrix products and `np.linalg`.
- Generate reproducible random numbers with `np.random.default_rng` and handle missing values (`NaN`) without distorting the results.
## 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:
| Property | Python list | NumPy `ndarray` |
|---|---|---|
| Element type | Any, mixed | A single one (`dtype`): `float64`, `int64`, `bool`… |
| Memory | Pointers to scattered objects | One contiguous block of numbers |
| Operations | Element by element, in Python | On the whole array, in compiled code |
| Dimensions | Lists of lists | Explicit `shape`: `(n,)`, `(n, k)`… |
The closes of SINT, the course's fictitious ETF, are read with pandas (lesson 0.5.2) and converted to NumPy with
`to_numpy()`:
```python
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]
```
```output
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_t = \frac{P_t}{P_{t-1}} - 1
$$
```python
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))
```
```output
(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:
```python skip
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:
```python
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%}")
```
```output
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`:
```python
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))
```
```output
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:
```python
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))
```
```output
['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.
| Shape A | Shape B | Result | Why |
|---|---|---|---|
| `(1302, 4)` | `(4,)` | `(1302, 4)` | 4 with 4; B has no rows: counts as 1 |
| `(1302, 4)` | `(1302, 1)` | `(1302, 4)` | 4 with 1; 1302 with 1302 |
| `(1302, 4)` | `(1302,)` | Error | 4 with 1302: neither equal nor 1 |
```python
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))
```
```output
(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:
```python raises
R - R[:, 0]
```
```output
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:
```python
relativa = R - R[:, 0][:, None]
print(relativa.shape)
print((relativa.mean(axis=0) * 252).round(4))
```
```output
(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:
```python
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))
```
```output
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:
```python
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%}")
```
```output
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:
```python
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%}")
```
```output
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:
```python
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())
```
```output
[ 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:
```python
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])
```
```output
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:
```python
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%}")
```
```output
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:
```python
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)
```
```output
[ 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`.
**Solution.**
```python
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))
```
```output
[0.1683 0.0496 0.1511 0.2827]
True
```
The variance is the mean of the squared deviations (with $n - 1$ in the denominator, as `np.cov` does by
default); it is the diagonal of the covariance matrix.
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?
**Solution.**
```python
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%}")
```
```output
Portfolio: 9.92%
Weighted average: 12.08%
```
The correlation between `RVAR` and `RFIJ` is negative (−0.19): some of the drops in one coincide with rises in the
other. The lower the correlation between the components, the further the portfolio volatility moves away from the
average of the volatilities.
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?
**Solution.**
```python
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%}")
```
```output
Losing years: 36.7%
Median: 6.7%
```
Each row is a simulated year; `np.prod(..., axis=1)` compounds the 252 returns of each row. The simulation assumes
normal, independent returns, a simplification that level 1 revisits with real data (fat tails,
volatility clustering).
## Self-assessment
1. `R` is a return matrix of shape `(1302, 4)` and `m = R.mean(axis=0)`. What shape does `R - m` have?
- A. `(4,)`
- B. `(1302,)`
- C. `(1302, 4)`: `m` is stretched along the rows
- D. Error: the shapes do not match
*Answer C: `(1302, 4)`: `m` is stretched along the rows.* Broadcasting compares the dimensions from the right: 4 with 4 matches and the dimension missing in `m` is treated as 1, so each column's mean is subtracted from every row.
2. Why does the NumPy documentation discourage `np.random.seed(42)` in new code?
- A. Because it is slower than `default_rng`
- B. Because it sets a global state shared by the whole program: any other call to `np.random` changes the sequence and the results stop being reproducible piece by piece
- C. Because it no longer exists in NumPy 2
- D. Because it only accepts small seeds
*Answer B: Because it sets a global state shared by the whole program: any other call to `np.random` changes the sequence and the results stop being reproducible piece by piece.* `np.random.seed` still exists, but it drives a legacy global generator. A `Generator` of your own, created with `default_rng(seed)`, is passed as an argument and nobody else consumes it.
3. An array of closes has a `NaN` for a session with no quote. What does `np.mean` of its returns return?
- A. The mean without that day
- B. Zero
- C. `nan`: the missing value propagates to any operation that touches it
- D. An error
*Answer C: `nan`: the missing value propagates to any operation that touches it.* A `NaN` contaminates the whole operation. To ignore it you have to ask for it (`np.nanmean`) or filter it out with a mask; never fill it with zero without thinking, because a zero is a data point.
4. A portfolio has weights `w` and annualized covariance matrix `S`. Which expression gives its annual volatility?
- A. `np.sqrt(w @ S @ w)`
- B. `(w * np.sqrt(np.diag(S))).sum()`
- C. `S.sum() * w.sum()`
- D. `np.linalg.solve(S, w)`
*Answer A: `np.sqrt(w @ S @ w)`.* The portfolio variance is the quadratic form w′Σw. The second option is the weighted average of the volatilities, which ignores correlations and overstates the risk of a diversified portfolio.
## Key points
- An `ndarray` stores numbers of a single type in contiguous memory; operations apply to the whole array in compiled code, without Python loops.
- Broadcasting aligns shapes from the right: a per-column mean `(4,)` is subtracted from an `(n, 4)` matrix without copying it.
- A portfolio's volatility is `np.sqrt(w @ S @ w)`: it depends on the correlations and is usually lower than the average of the volatilities.
- Reproducible randomness: a `Generator` of your own with `default_rng(seed)`, passed as an argument; never the global state of `np.random.seed`.
- `NaN` propagates and comparisons with it are always false: count the missing values before computing and decide what to do with them.
## Resources
### Documentation
- [NumPy: broadcasting (official documentation)](https://numpy.org/doc/stable/user/basics.broadcasting.html) · numpy.org
- [NumPy: Random Generator](https://numpy.org/doc/stable/reference/random/generator.html) · numpy.org
- [NEP 19: Random number generator policy](https://numpy.org/neps/nep-0019-rng-policy.html) · numpy.org
- [NumPy: linear algebra (numpy.linalg)](https://numpy.org/doc/stable/reference/routines.linalg.html) · numpy.org
### Books
- [Wes McKinney, Python for Data Analysis, 3rd ed. (O'Reilly, 2022; open access edition)](https://wesmckinney.com/book/) · wesmckinney.com
### Data
- [Level 0 data: closes of the fictitious four-asset portfolio (CSV)](https://pythonparatrading.com/curso/datos/nivel-0/cartera_sint.csv) · Download · cartera_sint.csv
Previous lesson: [0.4.6 · Concurrency: threading, multiprocessing, asyncio and free-threading](https://pythonparatrading.com/en/python-trading-course/level-0/python-concurrency-threading-multiprocessing-asyncio/)
Next lesson: [0.5.2 · pandas 3.0: Series, DataFrame, time indexes, Copy-on-Write, the str type and pd.col()](https://pythonparatrading.com/en/python-trading-course/level-0/pandas-3-series-and-dataframe/)