Course 1 · Lesson 5 Vectorized Numerical Computing with NumPy
Python for Actuarial Students · Course 1
Lesson 5

Vectorized Numerical Computing with NumPy

Whole-array arithmetic: faster code, fewer loops, and formulas that look like the math on the page.

Lesson 5 · OverviewBy the end of this lesson you can…

  • Explain why NumPy arrays are faster than Python lists
  • Create arrays and reason about dtype, shape, and ndim
  • Index, slice, and mask arrays — and know when you have a view or a copy
  • Predict the result of broadcasting two arrays of different shapes
  • Use dot products, matrix multiplication, and transposes
  • Vectorize bond pricing, cohort mortality, and portfolio returns
75 min
Concepts & worked examples
75 min
Lab: three vectorized models

Companion notebook: lesson-05-numpy.ipynb — worked examples and self-checking exercises.

Part 1

Why NumPy?

What a Python list really stores — and what an array stores instead.

Under the hoodLists of objects vs. contiguous arrays

Python list of floats

A list of pointers, each to a separate float object (24 bytes each) scattered around memory. Every + is looked up and dispatched at run time.

NumPy ndarray

One contiguous block of raw 8-byte numbers (C-contiguous by default). Operations run in compiled C loops over the whole block.

import numpy as np

xs = list(range(1_000_000))
arr = np.arange(1_000_000, dtype=float)

sum(x * x for x in xs)    # ≈ 38 ms
(arr * arr).sum()         # ≈ 1 ms

arr.itemsize, arr.nbytes  # (8, 8000000)

Roughly 30–50× faster here (timings vary by machine) — and the code is shorter. Vectorized code is usually both quicker and clearer.

Part 2

Arrays: Creation, Shape & Slicing

And the difference between a view and a copy.

CreationMaking arrays

q = np.array([0.0010, 0.0012, 0.0015, 0.0019, 0.0024, 0.0031])
q.dtype, q.shape, q.ndim      # (dtype('float64'), (6,), 1)

np.arange(1, 6)               # [1 2 3 4 5]
np.linspace(0, 1, 5)          # [0.   0.25 0.5  0.75 1.  ]
np.zeros(3)                   # [0. 0. 0.]
np.ones((2, 3))               # 2 rows × 3 columns of 1.
np.full(4, 100.0)             # [100. 100. 100. 100.]

m = np.arange(12).reshape(3, 4)
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]
dtypeelement type — one per array (float64, int64, bool)
shapesize along each axis — a tuple
ndimnumber of axes

A shape (6,) array is a 1-D vector — not a row (1, 6) or a column (6, 1). That distinction matters for broadcasting.

np.arange excludes the stop, like range; np.linspace includes both ends.

IndexingSlices, rows, columns, and masks

m[1]          # [4 5 6 7]      a row
m[:, 2]       # [ 2  6 10]     a column
m[1, 2]       # 6              one element: [row, col]
m[0:2, 1:3]   # [[1 2]
              #  [5 6]]

# Boolean mask: keep where the condition is True
q[q > 0.0015]      # [0.0019 0.0024 0.0031]

# Fancy indexing: pick by position list
q[[0, 2, 4]]       # [0.001  0.0015 0.0024]
  • One index per axis, separated by commas
  • : alone means “everything along this axis”
  • Comparisons produce boolean arrays; indexing with one filters
  • Combine masks with &, |, ~ — and parentheses: q[(q > 0.001) & (q < 0.003)]

Masks replace filter loops: claims[claims > deductible] - deductible in one line.

GotchaViews vs. copies

q = np.array([0.0010, 0.0012, 0.0015, 0.0019, 0.0024, 0.0031])

window = q[1:4]        # a VIEW: same memory
window[0] = 99
q
# [1.0e-03 9.9e+01 1.5e-03 1.9e-03 2.4e-03 3.1e-03]

q = np.array([0.0010, 0.0012, 0.0015, 0.0019, 0.0024, 0.0031])
safe = q[1:4].copy()   # a COPY: independent
safe[0] = 99
q                      # unchanged

np.shares_memory(q, q[1:4])        # True
np.shares_memory(q, q[[1, 2, 3]])  # False
OperationReturns
Basic slice a[1:4]view
reshape, .Tview (usually)
Mask a[a > 0]copy
Fancy a[[1, 2]]copy
Arithmetic a * 2new array

The opposite of Python lists, where a slice is always a copy. Views are why NumPy slicing is cheap — and why stressing a “slice” of base assumptions can change the base.

Part 3

Broadcasting & Element-wise Math

How arrays of different shapes combine.

BroadcastingThe rules

Compare shapes from the right. Two dimensions are compatible when they are equal, or one of them is 1. Size-1 dimensions are stretched to match.

ABResult
(6,)scalar(6,)
(3, 1)(5,)(3, 5)
(3, 1)(4,)(3, 4)
(3,)(4,)error
 
ValueError: operands could not be broadcast together
with shapes (3,) (4,)

A whole table of discount factors — three rates × five years — with no loop:

rates = np.array([0.03, 0.05, 0.07])[:, None]  # (3, 1)
t = np.arange(1, 6)                            # (5,)

D = (1 + rates) ** -t                          # (3, 5)
D.round(4)
# [[0.9709 0.9426 0.9151 0.8885 0.8626]
#  [0.9524 0.907  0.8638 0.8227 0.7835]
#  [0.9346 0.8734 0.8163 0.7629 0.713 ]]

[:, None] adds an axis, turning a (3,) vector into a (3, 1) column.

Element-wiseUniversal functions and aggregates

x = np.array([0.02, 0.05, 0.10])

np.exp(x), np.log1p(x), np.sqrt(x)   # element-wise "ufuncs"
np.maximum(claims - 500, 0)          # vectorized max(loss - d, 0)
np.minimum(claims, 10_000)           # cap at a limit
np.where(x > 0.04, "high", "low")   # vectorized if/else

D.sum()            # everything
D.sum(axis=0)      # down the rows → one value per column
D.sum(axis=1)      # across columns → one value per row
D.mean(), D.max(), D.cumsum(), D.cumprod()
  • ufuncs apply element by element in C
  • Use np.maximum/np.minimum — not Python’s max/min — on arrays
  • axis=0 collapses rows; axis=1 collapses columns

D.sum(axis=1) on the discount table gives the 5-year annuity factor at each of the three rates. Lesson 1’s deductible-and-limit becomes np.minimum(np.maximum(loss - d, 0), u) for a million losses at once.

Linear algebraDot products, matrix products, transposes

x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
x @ y                  # 32.0      dot product x · y

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

A @ B      # [[19 22]    matrix product
           #  [43 50]]
A * B      # [[ 5 12]    element-wise — different!
           #  [21 32]]
A.T        # [[1 3]      transpose
           #  [2 4]]

Shapes must line up: (m, k) @ (k, n) → (m, n)

(4, 5)4 bonds × 5 periods of cash flows
@
(5,)discount factors
=
(4,)4 prices

* multiplies element by element; @ is matrix multiplication. Mixing them up gives a result of the wrong shape — or, worse, the right shape and wrong numbers.

Actuarial Lab · 75 minThree vectorized models

1 · Multi-period pricing

Price a whole portfolio of bonds with a single matrix product.

2 · Cohort mortality

Build survival probabilities with cumprod, then value term insurance and an annuity with dot products.

3 · Portfolio returns

Combine 10,000 scenarios of three asset returns into a portfolio distribution.

The rule for today: no Python for loops over cash flows, ages, or scenarios. If you write one, find the array operation that replaces it.

Lab · Part 1Pricing a bond portfolio in one line

# rows = bonds, columns = years 1..5
CF = np.array([
    [5, 5,   5, 5, 105],   # 5y, 5% coupon
    [3, 3, 103, 0,   0],   # 3y, 3% coupon
    [0, 0,   0, 0, 100],   # 5y zero-coupon
    [8, 8,   8, 8, 108],   # 5y, 8% coupon
], dtype=float)

d = 1.04 ** -np.arange(1, 6)     # flat 4% curve

prices = CF @ d
prices.round(4)
# [104.4518  97.2249  82.1927 117.8073]

Tasks

  1. Reproduce the prices above
  2. Check row 0 against your Lesson 4 CouponBond(100, 0.05, 5, freq=1).price(0.04)
  3. Price all four bonds at 3%, 4%, and 5% at once: build a (3, 5) discount matrix by broadcasting, then CF @ D.T → shape (4, 3)
  4. Macaulay durations of all four: (CF * d * t).sum(axis=1) / prices

Real portfolios have thousands of bonds and hundreds of periods. The code doesn’t change.

Lab · Part 2Cohort mortality with dot products

The SOA’s Standard Ultimate Survival Model (Makeham): A = 0.00022, B = 2.7×10−6, c = 1.124.

$$p_x = \exp\!\Big(-A - \frac{B c^x (c-1)}{\ln c}\Big)$$

A, B, c = 0.00022, 2.7e-6, 1.124
ages = np.arange(60, 70)
px = np.exp(-A - B * c**ages * (c - 1) / np.log(c))
qx = 1 - px
# q_60 = 0.003398 ... q_69 = 0.009294

kpx = np.concatenate([[1.0], np.cumprod(px)[:-1]])  # _k p_60
v = 1.05 ** -np.arange(1, 11)

A_term = v @ (kpx * qx)          # 10-year term, age 60
a_due = (1.05 ** -np.arange(10)) @ kpx
Quantity (age 60, 10 years, 5%)Value
$A^{1}_{60:\overline{10}|}$0.042521
$\ddot{a}_{60:\overline{10}|}$7.955548
EPV of 100,000 death benefit4,252.09
Level annual net premium534.48

Tasks

  1. Reproduce the table
  2. A cohort of l60 = 100,000: build lx with cumprod and deaths dx with np.diff — 5,745.08 deaths over 10 years
  3. Premiums for ages 30–70 in one go (broadcast ages against durations)

Lab · Part 3Multi-asset portfolio returns

rng = np.random.default_rng(2026)     # next lesson's topic!

mu = np.array([0.06, 0.03, 0.08])     # equities, bonds, property
sd = np.array([0.15, 0.05, 0.20])
corr = np.array([[1.0, 0.2, 0.7],
                 [0.2, 1.0, 0.1],
                 [0.7, 0.1, 1.0]])
cov = np.outer(sd, sd) * corr

R = rng.multivariate_normal(mu, cov, size=10_000)  # (10000, 3)
w = np.array([0.5, 0.3, 0.2])                      # weights

port = R @ w                                       # (10000,)
SimulatedTheory
Mean return0.0564w · μ = 0.0550
Std. deviation0.1108√(wᵀΣw) = 0.1105
5th percentile−0.1253—

Tasks

  1. Reproduce the table (same seed → same numbers)
  2. Compute the theory column with @ and np.sqrt
  3. Try 1,000 random weight vectors at once: R @ W.T with W of shape (1000, 3), rows summing to 1

RecapLesson 5 checklist

  • I can explain why arrays beat lists for numbers
  • I check shape before I combine arrays
  • I can slice rows, columns, and use boolean masks
  • I know which operations return views and use .copy() when needed
  • I can predict a broadcast result shape
  • I use axis= deliberately in aggregates
  • I know @ from *
  • I priced bonds, a cohort, and a portfolio without a single loop
Up next
Lesson 6

Discrete Distributions & Stochastic Simulation

Reproducible random numbers, claim-frequency models, and simulating 10,000 lives to price stop-loss cover.

← Course home