2Python and pandas

Answer to Question 2.1.

It prints [1], [1, 2], [3] and [1, 2, 4]. A default value is evaluated once, when the def statement runs, so every call that omits the list shares one list, which outlives the call.

The third call passes its own list, so it does not touch the shared one, which is why the fourth continues from [1, 2].

The fix is the usual sentinel: write items=None in the signature, and inside the function say if items is None: items = []. Then a new list is made on each call. Test is None and not if not items, because a caller who passes an empty list expects that same list back.

Answer to Question 2.2.

It prints {'A': [1], 'B': [1]} and [[9, 0], [9, 0]]. In both cases one list object is shared, so a change made through one name is seen through all.

dict.fromkeys evaluates its second argument once and stores that one object under every key. [[0] * 2] * 2 repeats a reference to the same inner row. Python does not copy on repetition, and it never copies on assignment.

The fixes make one list per key or row: {k: [] for k in names}, or a defaultdict(list), and [[0] * 2 for _ in range(2)], where the expression runs once per row.

The 0 is safe because an integer cannot change. d[k] += 1 rebinds that one key to a new integer. It does not alter the shared 0. So sharing does no harm for immutable values, and the inner [0] * 2 in the row above is fine for the same reason.

Answer to Question 2.3.

It prints True True, then False True. == compares values, and is compares identity: whether the two names hold the very same object.

CPython keeps one shared object for each integer from $-5$ to 256, so int(…) returns the cached object for 256, and a new object each time for 257. That is an implementation detail, and other implementations or versions need not follow it.

The listing goes through int on text on purpose. With the literals 257 is 257 in one block, the compiler can merge the constants, and you get True, which adds to the confusion.

Use is for singletons, where identity is the meaning: x is None, or a sentinel object of your own. Use == for everything else, including True and False from a library, where a NumPy boolean is not the object True.

Answer to Question 2.4.

It prints tuple, then a TypeError for the list, a TypeError for the tuple holding a list, and absent for the frozen set. A key must be hashable, and a list is not.

A dictionary stores a key by its hash, and it finds the key again by recomputing that hash. If a key could change after insertion, its hash would change, and the entry would sit in the wrong place, unreachable. So mutable built-ins refuse to be hashed. The check comes on any lookup, as get must hash its argument too.

A tuple’s hash is built from its elements, so a tuple is hashable only if all its elements are. (1, [2]) fails for that reason, and it has no check on being “a tuple”.

For memory, a tuple is smaller, as it needs no spare room for growth. On CPython 3.10, 64-bit, sys.getsizeof gives 56 for a tuple of two and 72 for a list made of two by a list display. It makes little difference for one object, and a lot for a million small records.

Answer to Question 2.5.

It prints {'A': [1, 2, 3]} and {'A': [1, 2]}. The shallow copy is a new dictionary holding the same lists, so it does not see the new key B, but it does see the append. The deep copy has new lists all the way down.

For the snapshot, deepcopy walks the whole structure with a memo table, at a Python-level cost for every object, so taken many times a second on a large book it will show up on a profile.

Better options depend on what is read. If a snapshot is only for one consumer, send a compact form: a tuple of (price, size) pairs for each side. Tuples of immutable values can be shared safely, so copying just the top level is then enough. Or hold levels as immutable values, replace one on each change, and take the shallow copy at $O(\text{levels})$.

Answer to Question 2.6.

It prints False True False, then nan 1.0, then [2.0, nan, 1.0]. NaN compares false with everything, itself included, and that one rule explains all three lines.

nan == nan is false. nan in [nan] is true, since in first checks identity and the object is the same, while a fresh float('nan') is another object, and equality fails.

max keeps its first argument unless a later one compares greater. Nothing compares greater than NaN, so max(nan, 1.0) stays at NaN; and NaN does not compare greater than 1.0, so max(1.0, nan) stays at 1.0. So the result depends on where the NaN falls.

sorted sees every comparison with NaN as false, so it takes the list as already sorted, and leaves it unsorted, with no error. The list then gives a wrong minimum, a wrong maximum and a wrong median.

So for the largest price, filter first: max(x for x in prices if not math.isnan(x)), or np.nanmax. Do not call max and hope.

Answer to Question 2.7.

It prints [(1, 3), (2, 4)], [], 5 0 and 3 3. zip and a generator are iterators, which can be walked once. A range is a sequence, which starts a new iterator each time.

So on a generator of 1, 2 and 3, the function counts 3, which uses up the stream, then sums nothing and returns $0 / 3 = 0.0$, a plausible mean. With the sum first, it gets the right sum and a length of 0, and raises ZeroDivisionError, which is the kinder failure.

The fix depends on the size. If the data fits, xs = list(xs) at the top. If not, do a single pass that keeps a running count and a running total. Testing iter(xs) is xs tells you whether you hold a one-shot iterator, if a function must refuse one.

Answer to Question 2.8.

It prints [-2147483648], [255] and 2147483648. NumPy arrays have fixed-width types and wrap around with no warning, while Python’s own integers grow as needed.

An int32 array stops at $2^{31} - 1$, and adding 1 wraps to the lowest value. A uint8 array cannot hold a negative value, so $1 - 2$ wraps to 255. Array arithmetic gives no warning and no exception, since a check per element would cost too much. A NumPy scalar is different: np.int32(2**31 - 1) + 1 does give a RuntimeWarning, so a check on a single value can make you think all is well.

This bites where a column’s type was chosen for one size and then accumulated. A running volume in int32 wraps to a large negative number as soon as it passes about 2.1 billion. Timestamps in nanoseconds need 64 bits. A difference of unsigned sequence numbers gives a huge positive value in place of a small negative one.

Answer to Question 2.9.

Only C works, and it returns the row with qty 2. A and B both raise ValueError: “The truth value of a Series is ambiguous”.

Python’s and has to reduce its left side to one true or false, so it calls bool on a Series, and pandas refuses, as one answer per row cannot make one truth value. There is no way for a library to override and. So pandas uses the bitwise operator &, which it can overload to work row by row.

B fails through precedence. & binds tighter than > and <, so the line means df.qty > (1 & df.px) < 6. That is a chained comparison, which Python expands with an and, so the same error appears. If the columns held floats, 1 & df.px would raise a TypeError first.

C has the brackets, so each comparison gives a Series of booleans, and & combines them. The row with qty 2 has px 5, which passes. The row with qty 3 fails on px.

Answer to Question 2.10.

You get no speed-up, and often a small slowdown. CPython 3.10 has a global interpreter lock: only one thread runs Python bytecode at a time. The lock does not make count += 1 safe, because that statement is several bytecode steps.

The disassembly of count += 1 shows a load, an add and a store, and nothing in the language promises that a switch cannot fall between them. Four threads running the bare statement lost none of 8 million increments, in three runs here on CPython 3.10.12, which is the trap: the test passes, and that is no guarantee. Make the threads wait for each other between the read and the write, and our test loses three updates in every four, every time. So use a threading.Lock, or give each thread its own counter and add them at the end.

Threads help where the lock is released: waiting on a socket or a disk, and inside NumPy operations and many C extensions, which run in parallel on large arrays. For pure Python computation, use processes, at the price of sending data between them.

Answer to Question 2.11.

Listing 2.18 takes the calls of a recursive Fibonacci for $n = 30$ from 2 692 537 to 31. A list argument raises TypeError, since it cannot be hashed. f(1, 2) and f(1, b=2) are two entries: right answers, computed twice.

from functools import wraps


def memo(f):
    cache = {}

    @wraps(f)
    def inner(*args, **kwargs):
        key = (args, tuple(sorted(kwargs.items())))
        if key not in cache:
            cache[key] = f(*args, **kwargs)
        return cache[key]

    return inner

Without a cache, the number of calls satisfies $c(n) = c(n - 1) + c(n - 2) + 1$, which comes to $2F_{n+1} - 1$, so 2 692 537 for $n = 30$. With it, each of the 31 values from 0 to 30 is computed once. The recursive calls go through the decorated name, which is why it works at all.

The test key not in cache matters because a result can be 0, None or an empty list, and the form cache.get(key) or f(…) recomputes on each of those, Fibonacci of 0 included.

wraps copies the name and the docstring, so that a stack trace and a profile show the real function, and not inner. To accept lists, convert them to tuples in the key, or leave the caching to the caller.

Answer to Question 2.12.

Listing 2.19 reads the file one row at a time and keeps one running total per symbol, so memory depends on the number of distinct symbols, not on the size of the file.

import csv
from collections import defaultdict


def rows(path):
    with open(path, newline="") as f:
        yield from csv.DictReader(f)


def notional_by_symbol(records):
    totals = defaultdict(int)
    for r in records:
        totals[r["symbol"]] += int(r["price_ticks"]) * int(r["qty"])
    return dict(totals)

rows is a generator, so nothing is read until notional_by_symbol asks for it, and each row can be freed once it has been added. In our test, on a file of 200 000 rows, the peak allocation stays below 1 MB, while list(rows(path)) uses more than 10 MB, and that gap grows with the file. Prices stay as integer ticks, so the totals are exact.

If the caller stops halfway, the file stays open while the generator is alive. It closes when the generator is closed, by close() or by CPython’s reference counting when the last reference goes. Both raise GeneratorExit at the yield, which runs the with. Other Python implementations may delay that, so a long-lived caller should call close() or wrap it in contextlib.closing.

Answer to Question 2.13.

It prints [1], then [9] [1]. orders is one list on the class, shared by all instances, so a.add(1) appends to it and b sees that. The assignment a.orders = [9] creates a new attribute on a, which hides the class’s list for a alone.

So reading and appending go to the class’s object, while assignment always writes on the instance. That is the whole of the inconsistency. For a plain class attribute, reading checks the instance first and then the class, and assignment writes on the instance. A property or other data descriptor on the class takes precedence over both. The fix is to set self.orders = [] in __init__.

With __slots__ = ('name', 'orders'), instances have no __dict__, just fixed fields. In our test, a slotted instance takes less than half the memory of an ordinary one, on CPython 3.10. A typo such as self.qyt = 1 raises AttributeError at once.

A name in __slots__ cannot also be a class attribute with a default: Python raises ValueError when it builds the class, so the original bug cannot be carried across.

Answer to Question 2.14.

Listing 2.20 shows both forms, and each logs “release” when the body raises. If __exit__ returns a true value, the exception is swallowed, and execution continues after the with block as though nothing went wrong.

from contextlib import contextmanager


class Held:
    def __init__(self, log):
        self.log = log

    def __enter__(self):
        self.log.append("acquire")
        return self

    def __exit__(self, exc_type, exc, tb):
        self.log.append("release")
        return False                  # True would swallow the exception


@contextmanager
def held(log):
    log.append("acquire")
    try:
        yield
    finally:
        log.append("release")

Python calls __exit__ however the block ends, passing the exception’s details if there is one. Returning nothing, which is None, lets the exception go on, and so does False. So the safe rule is to return true only when this manager means to handle that one exception, and to test its type.

In the contextlib version, an exception in the body is raised inside the generator at the yield. Without try and finally, it leaves the generator at that point, and the line after the yield never runs. So a lock stays held, or a file stays open, in exactly the failure case the manager was written for. The normal path shows nothing, so a test must raise on purpose to catch it.

Answer to Question 2.15.

A field tags: list = [] raises ValueError when the class is defined, not when an order is made. Listing 2.21 fixes it with field(default_factory=list). Objects of a plain dataclass cannot go in a set, because it has no hash: @dataclass sets __hash__ to None.

from dataclasses import dataclass, field


@dataclass
class Order:
    id: int
    qty: int
    tags: list = field(default_factory=list)


@dataclass(frozen=True)
class Key:
    symbol: str
    venue: str

The error at definition is Python protecting you from the shared default in a plain function: without it, all orders would share one list. A factory runs once per instance.

The hash goes because @dataclass writes an __eq__ by value. Two objects that compare equal must have the same hash, and a mutable value could change after it was put in a set, leaving it in the wrong place. So Python’s rule is: define equality without saying the object is immutable, and you get no hash.

frozen=True makes assignment to a field raise, and it generates __hash__ from the fields, so Key works in a set or as a dictionary key. It does not protect what the fields hold. A frozen object with a list field can still have that list appended to, and hashing it then raises TypeError, as hashing the list does.

Answer to Question 2.16.

Sort twice, least important key first, as in Listing 2.22: first by notional descending, then by symbol. Python’s sort is stable, so the second sort keeps the first sort’s order within each symbol.

from operator import itemgetter


def by_symbol_then_notional(trades):
    """Symbol ascending, notional descending; equal trades keep their order."""
    out = sorted(trades, key=itemgetter(1), reverse=True)
    out.sort(key=itemgetter(0))
    return out

A stable sort never swaps two elements with equal keys. So after sorting by notional and then by symbol, trades with one symbol stay ordered by notional. reverse=True also preserves the order of ties, by the language’s definition, so that first pass reverses nothing that ties. Writing sorted(trades, key=…)[::-1] does not have that property, as it flips ties too.

Where the second key is a number, one pass is enough: key=lambda t: (t[0], -t[1]). That fails for a string or a None, neither of which can be negated. The two-pass form works for any key that orders, so a column mixing None with numbers still needs a key that places the None.

Each pass costs $O(n \log n)$ at worst.

Answer to Question 2.17.

The first print gives [[-0.5, -2.5], [1.5, 1.5]] and the second [[-0.5, 0.5], [-1.5, 1.5]]. The second is right. Each row should add up to zero after you take its mean out, and in the first the rows add up to $-3$ and $3$.

r.mean(axis=1) has shape (2,). NumPy’s broadcasting lines up shapes from the right, so (2,) matches the columns of (2, 2), and it subtracts the first day’s mean from the first stock on all days, and the second day’s mean from the second stock. It is a meaningless result, with no error, because the matrix happens to be square.

keepdims=True gives shape (2, 1), which lines up with the rows and repeats across the columns.

With three stocks, the shapes are (2, 3) and (2,): the last dimensions are 3 and 2, which do not match, so NumPy raises ValueError. So the bug hides when the days equal the stocks, which can happen with a small test matrix or a short window, or in one window of a rolling calculation. It also hides with a single stock, where (n, 1) and (n,) broadcast to (n, n).

Answer to Question 2.18.

It prints [ 0 99 2 3 4]. Only b is a view: a slice shares memory with a, while indexing with a list or a boolean mask makes a copy, so the writes through c and d do not reach a.

The rule follows from what can be expressed. A slice is a start, a stop and a step, so a new header over the same buffer can describe it, at $O(1)$ cost. Arbitrary positions cannot be described that way, so NumPy copies them. np.shares_memory(a, b) tells you which you have.

small = huge[:10] keeps all 10 GB alive, since small.base points at the original buffer. Use huge[:10].copy() and then let huge go.

df[df.x > 0]['y'] = 1 on pandas 2.3 leaves df as it was, with a SettingWithCopyWarning. The mask gives a new object, and the assignment goes to that object. It is then thrown away. Write df.loc[df.x > 0, 'y'] = 1, which is one operation on df.

Answer to Question 2.19.

(A, 1) gives 10.75, (B, 1) gives 20 and (A, 2) gives NaN, as in Listing 2.23. groupby(…).price.mean() gives 10.5 for (A, 1), which treats a trade of 100 shares and a trade of 300 shares alike.

import pandas as pd


def vwap_by_symbol_day(df: pd.DataFrame) -> pd.Series:
    money = df["price"] * df["qty"]
    keys = [df["symbol"], df["day"]]
    return money.groupby(keys).sum() / df["qty"].groupby(keys).sum()

The weighted average is $(10 \times 100 + 11 \times 300) / 400 = 4300 / 400 = 10.75$. Compute the two sums per group and divide. The two results have the same group index, so the division lines up by group with no further work.

The NaN for (A, 2) is the honest answer: no shares traded, so no average exists, and 12 would be a price at which nothing was done. Pandas gives NaN for $0 / 0$ here, and it needs no guard. A plain Python version would raise, and a C++ one would divide by zero.

The usual alternative is groupby(…).apply(lambda g: …), which is correct and runs a Python function once per group, while the listing’s two sums run in compiled code.

Answer to Question 2.20.

By default the trades at 1, 5 and 10 get no quote, the quote at 5, and the quote at 9. Listing 2.24 is the usual call. With allow_exact_matches=False, the trade at 5 gets the quote at 2. With direction='forward' they get 2, 5 and no quote, which is the backtest bug.

import pandas as pd


def quotes_at_trades(trades, quotes):
    return pd.merge_asof(
        trades.sort_values("t"), quotes.sort_values("t"),
        on="t", by="symbol", direction="backward",
    )

A forward match gives each trade a quote from its future. A signal built on that looks excellent and cannot be earned, since at time 1 nobody yet knows the quote of time 2.

An exact match is a choice about what the timestamps mean. If a quote at 5 is published before a trade at 5, keep it. If your clocks cannot say which came first, then allow_exact_matches=False is the careful choice, and costs you little.

Both inputs must be sorted on the key, and pandas raises ValueError (“left keys must be sorted”) if they are not, which is why the listing sorts first. With by='symbol', matches stay inside a symbol.

Add tolerance, in practice. Otherwise, a trade after a long gap gets a quote that is hours old, and nothing says so.

Answer to Question 2.21.

The 09:30 bar holds the ticks from 09:30:00 up to but not including 09:31:00: open 10, high 12, low 10, close 12. By default pandas labels it 09:30, the start of the bar. A join on that label treats 12, a price seen at 09:30:50, as known at 09:30:00. Listing 2.25 labels each bar at its end.

import pandas as pd


def minute_bars(prices: pd.Series) -> pd.DataFrame:
    bars = prices.resample("1min", label="right", closed="left").ohlc()
    bars["close"] = bars["close"].ffill()
    return bars

So a bar is known only once it closes, and label='right' says so: the same bar becomes 09:31. A strategy reading the bar labelled 09:31 at 09:31 then knows only the past. closed is a separate matter, which side of the minute a tick at exactly 09:31:00 falls on, and here it stays left.

The bar for 09:32 has no ticks, and so all four columns are NaN. A close for an empty bar can sensibly carry the last price, as the listing does, since that is still the last known price, which gives 11. Open, high and low should stay NaN or be set equal to it, and any volume column should be 0. A strategy must also know that the bar was empty, as a filled price is not a market that was tradable then.

Answer to Question 2.22.

It prints [nan, 12.0, 23.0, nan], [11, 22, 33] and [1.0, 12.0, 23.0, 30.0]. The second line is the bug that raises no error. The first is pandas working as designed: it lines up by index, over the union of the two.

a + b matches labels: at 2 it adds 2 and 10, at 3 it adds 3 and 20, and labels in only one have no partner, so NaN. a.values + b.values throws the labels away and adds by position, so label 1 of one is added to label 2 of the other.

For a pair, that means comparing Monday’s price of one stock with Tuesday’s price of the other, from the missing day onward, for every later day. Spreads stay plausible, and the correlation falls, or an apparent lead-lag appears. There is no NaN and no warning, and shapes match, because the second series has a later day to make up the count.

So never leave the index when there are two series. Reindex both to one calendar, decide what a missing day means, and keep the NaN until you have decided. add(b, fill_value=0) is for quantities, where 0 means none. For prices it makes a nonsense value, and a 100% move.

Answer to Question 2.23.

You get an annual Sharpe of about 21, on pure noise. The fix is shift: signal.shift(1) * ret, which gives about 0. Today’s return was in the signal, so the backtest trades on tomorrow’s newspaper. Listing 2.26 runs it.

import numpy as np
import pandas as pd


def sharpes(n=25_200, seed=23):
    rng = np.random.default_rng(seed)
    ret = pd.Series(rng.normal(0, 0.01, n))     # pure noise: nothing to find
    signal = np.sign(ret)
    def annual(pnl):
        return pnl.mean() / pnl.std() * np.sqrt(252)
    return annual(signal * ret), annual(signal.shift(1) * ret)

With $r \sim N(0, \sigma^2)$, the profit is $\text{sign}(r) \times r = |r|$. Its mean is $\sigma\sqrt{2/\pi}$, and its variance is $\sigma^2(1 - 2/\pi)$. The daily ratio is $\sqrt{2/\pi} / \sqrt{1 - 2/\pi} = 1.3236$, and multiplying by $\sqrt{252}$ gives 21.0. The result is the same for any $\sigma$, since a Sharpe ratio has no units.

Our simulation over 25 200 days gives 21.1 and then, with the shift, $-0.03$. Noise run through the shifted version should show nothing, and that is the test to keep.

A clean bug-finding habit follows. Run every strategy on shuffled or simulated noise. Anything well above zero there has a data leak.

To make it hard to write, give the research library a single place that joins signals and returns, one that applies the lag itself, so a strategy can only return positions decided with data up to each date. Then a look-ahead has to be written on purpose.

Answer to Question 2.24.

It prints 4.0 2.0 2 3 1.33, then 0.0 nan, then {'A': 4}, then False.

sum, mean and count all skip NaN by default, so the mean is $4 / 2$, and count is 2 while len is 3. fillna(0) first changes the meaning: the mean is now $4 / 3$. The sum of only NaN is 0, by the same skipping, which can make a missing day look like a real zero; min_count=1 makes it NaN.

groupby drops rows with a missing key by default, so the row with value 2 disappears from the total, with no warning. Use dropna=False to keep it. And s == np.nan is false everywhere, so use s.isna().

prices.fillna(0) says that the price fell to zero. The return on that day is $-100\%$, and on the next day it is a division by zero, which pandas gives as inf. Sums and means of those columns are then useless. For prices, leave the NaN, or carry the last price forward with ffill if a stale price is acceptable, and say so. Compute returns after that step, and do not fill the returns with the mean, as that invents data.

Answer to Question 2.25.

Use category for symbol, keep price as float64, and hold qty as int64 or int32 depending on its range. Prices take 800 MB, quantities 800 MB or 400 MB, and symbols 200 MB, against about 1 GB for a plain string column from read_csv and about 6 GB with a string per row.

A plain string column holds 8 bytes per row for the pointer, plus the string objects. Built from Python’s csv module, each row gets its own string, of about 53 bytes for a short one on CPython 3.10, so that comes to about 6 GB. The C parser of pd.read_csv shares equal strings, so there it is nearer 1 GB, though memory_usage(deep=True) still reports about 6 GB, as it counts each element. category stores each distinct string once, and a small integer code per row. For 4 000 values the code is int16: 2 bytes a row, or 200 MB. In our test on 200 000 rows with a string per row it is more than ten times smaller, and against read_csv’s shared strings it saves about four times.

A float32 has 24 bits of mantissa. A price of 101.27 reads back as 101.26999664…, so an equality with a tick grid fails, and a float price near 10 000 has a step of about 0.001. For totals it is much worse. Above $2^{24} = 16\,777\,216$, not every integer exists, so adding 1 to that value leaves it unchanged: a running volume simply stops growing, with no error.

Answer to Question 2.26.

Listing 2.27 is such a test file, and it passes against the function. 0.1 + 0.2 is 0.30000000000000004, since none of the three values has an exact binary form. pytest.approx allows a relative tolerance of $10^{-6}$ or an absolute one of $10^{-12}$, whichever is larger. So 1e-13 == approx(0) passes, and 1e-9 == approx(0) fails.

import math

import pandas as pd
import pytest
from iv_vwap import vwap_by_symbol_day


def table(*rows):
    return pd.DataFrame(rows, columns=["symbol", "day", "price", "qty"])


@pytest.mark.parametrize("rows, want", [
    ([("A", 1, 10.0, 100)], 10.0),                       # one trade
    ([("A", 1, 10.0, 100), ("A", 1, 11.0, 300)], 10.75), # weighted, not 10.5
    ([("A", 1, 0.1, 1), ("A", 1, 0.2, 1)], 0.15),        # needs approx
])
def test_one_group(rows, want):
    got = vwap_by_symbol_day(table(*rows))
    assert got.loc[("A", 1)] == pytest.approx(want)


def test_zero_quantity_is_not_a_number():
    assert math.isnan(vwap_by_symbol_day(table(("A", 1, 10.0, 0))).loc[("A", 1)])

parametrize runs one body over a table of cases, and it reports each as a separate test, so one failing row does not hide the others. The cases are chosen to fail a wrong model: 10.75 fails the plain mean, as that gives 10.5.

A relative tolerance is useless near zero, because $10^{-6}$ of zero is zero, and this is why the absolute limit exists. Against a target of 0, anything up to $10^{-12}$ passes, and a value like $10^{-9}$ fails even if it is only rounding noise, so pass abs= when the scale calls for it.

Add the cases that a table of ordinary rows misses: zero quantity, one row, an empty table, several groups at once, and a group that comes in two separated blocks. Also check a limit with a known answer: if all prices are equal, the average is that price whatever the sizes.

Answer to Question 2.27.

It prints 2 4 2.67, then the exact double, 2.67499999999999982…, then False, and then 2.68. Python’s round sends halves to the even neighbour, and 2.675 as a double is not a half at all.

Round half to even keeps long columns of rounded values from drifting upward, so round(2.5) is 2 and round(3.5) is 4. NumPy’s round also sends exact halves to the even neighbour, so “it rounds up” is wrong on both. With decimals it scales first, so on NumPy 2.2 np.round(2.675, 2) is 2.68, where Python gives 2.67.

round(2.675, 2) would be a tie only if the input really were 2.675. The nearest double lies just below it, as the second line shows, so the correct rounding of that double is 2.67. Decimal(2.675) converts the double exactly and shows you the problem, which is why Decimal(0.1) does not equal Decimal('0.1').

So the fee must come in as text, or as an integer number of smaller units. Decimal('2.675') holds the exact value, and quantize with ROUND_HALF_UP gives 2.68. That is the fourth line. State the rounding rule in the code, since a fee schedule or a rule book sets it, and not the language.

Answer to Question 2.28.

sum returns 0.0, and math.fsum returns 1.0, the true answer. Doubles near $10^{16}$ are 2 apart, so $10^{16} + 1$ rounds back to $10^{16}$, and the 1 is gone before the subtraction. Listing 2.28 is Neumaier’s form of compensated summation, which also returns 1.0. Classic Kahan returns 0.0 on this list.

def neumaier(xs):
    total = lost = 0.0
    for x in xs:
        t = total + x
        if abs(total) >= abs(x):
            lost += (total - t) + x   # the low bits of x that the addition dropped
        else:
            lost += (x - t) + total
        total = t
    return total + lost

Each addition loses some low bits, and those bits can be computed exactly: for two doubles, the larger minus the rounded total, plus the smaller, is exactly what was dropped. The listing adds those to a separate small total, and only combines the two at the end.

Kahan folds the correction into the next input, so the next value, $-10^{16}$, with its correction of $+1$ sits at $-10^{16} + 1$, which rounds back. The loss is real: its correction assumes that the running total is larger than each new term, and here the new term is as large as the total.

Neither is exact in all cases, and fsum is: it keeps enough partial sums to give the correctly rounded result, at some cost in speed. On a million copies of 0.1, plain sum is off by about $1.3 \times 10^{-6}$, while fsum and the listing both give 100 000.0.

Answer to Question 2.29.

The heartbeat stops for the whole 200 ms. An event loop is cooperative: on one thread, one coroutine runs until it reaches an await that yields, and blocking work has none. Listing 2.29 shows this as an order of events: with blocking, all the work runs before the next beat.

import asyncio


async def heartbeat(log):
    for _ in range(3):
        log.append("beat")
        await asyncio.sleep(0)


def blocking_work(log):
    log.append("work")


async def handler_blocks(log):
    for _ in range(3):
        blocking_work(log)            # no await: nothing else can run during this


async def handler_yields(log):
    for _ in range(3):
        await asyncio.to_thread(blocking_work, log)


async def main(handler):
    log = []
    await asyncio.gather(heartbeat(log), handler(log))
    return log

On the test’s log, handler_blocks produces one beat, three pieces of work and then two beats. No clock is needed to show it. An exchange that sees no heartbeat for that long may disconnect the session, and every other coroutine has stalled as well, including the ones reading market data.

For blocking input or output, use an asynchronous library, or asyncio.to_thread, which runs the call on a worker thread and yields to the loop while waiting. Since the lock is released during a real wait, the loop keeps running. For pure Python computation that does not help for long, because of the interpreter lock: the thread still competes with the loop. So send it to a process with run_in_executor and a ProcessPoolExecutor.

Answer to Question 2.30.

The right slope is 2, using the three complete rows, and Listing 2.30 computes it. Demeaning each column on its own gives 0.00056. A symbol with one row returns NaN, as does one with a constant x.

import pandas as pd


def group_betas(df: pd.DataFrame) -> pd.Series:
    d = df.dropna(subset=["x", "y"])          # a row counts only if it has both
    means = d.groupby("symbol")[["x", "y"]].transform("mean")
    dx, dy = d["x"] - means["x"], d["y"] - means["y"]
    by = d["symbol"]
    sxx = (dx * dx).groupby(by).sum()
    flat = sxx <= 1e-24 * (d["x"] * d["x"]).groupby(by).sum()    # constant x, up to rounding
    return ((dx * dy).groupby(by).sum() / sxx).mask(flat)

The slope is

with both means taken over the rows used. transform('mean') gives each row its group’s mean, and the rest is arithmetic on columns plus two group sums, all in compiled code.

The error comes from the row (100, NaN). Column by column, pandas takes $\bar x = 26.5$ over all four rows and adds the 100 into the denominator, while the numerator skips that row because its product is NaN. So the two sums cover different rows, and nothing warns you. Dropping incomplete rows first makes every sum cover the same rows.

A slope from one row or from a constant x has no meaning, and NaN says so. Both would be $0 / 0$ in exact arithmetic, but in floats a group’s mean can be off by one bit, and a constant x then gives a finite slope of pure rounding. So the listing also returns NaN where the sum of $dx^2$ is too small, relative to the sum of $x^2$, to be anything else. A filter on a minimum row count belongs after this step, as a slope from 3 rows prints just like a slope from 3 000.

Answer to Question 2.31.

Remove impossible rows first, then sort, then remove duplicates, then test outliers against past values only, and resample within each symbol and day, as in Listing 2.31. With a limit of 5% on a median of 20, a true jump of 6% loses its first 10 ticks, and the rest pass.

import numpy as np


def clean(q, window=20, limit=0.05):
    q = q[(q.bid > 0) & (q.ask > q.bid)]                   # bad sides, locked, crossed
    q = q.sort_values(["symbol", "t"], kind="stable")
    q = q.drop_duplicates(["symbol", "t"], keep="last")
    q = q.assign(mid=(q.bid + q.ask) / 2)
    past = q.groupby("symbol").mid.transform(              # the past only: no centring
        lambda s: s.rolling(window, min_periods=5).median().shift(1))
    return q[~((q.mid / past - 1).abs() > limit)]          # NaN compares false: kept


def minute_returns(q, stale_minutes=5):
    q = q.assign(day=q.t.dt.normalize()).set_index("t")
    keys = ["symbol", "day"]
    last = q.groupby(keys).mid.resample("1min", label="right", closed="left").last()
    filled = last.groupby(level=keys).ffill(limit=stale_minutes)
    return np.log(filled).groupby(level=keys).diff()       # never across a night

The order matters. A zero side left in would pull the median, and a late row left unsorted would be compared with the wrong neighbours. keep='last' after a stable sort keeps the last to arrive at one timestamp, which is a choice worth stating.

A centred window compares a tick with ticks that came after it. In a backtest, that removes bad prints using knowledge you did not have then, and it makes the cleaned history smoother than anything you could trade.

The day key stops diff from crossing a night. The carry forward has a limit of 5 minutes, so a stale mid is not presented as a quiet market for hours.

The first five ticks of a symbol have no past median and pass unchecked, so a bad print there survives. Test them against a reference price, such as the previous close.

On the jump, the past median stays at the old level until 10 of the last 20 mids are new. Those ticks look like errors, so a real move and a bad print start out alike. So count rejections per symbol and alert on them, and never let the filter run silently, since a long run of rejects is news.

Answer to Question 2.32.

With $a = 0.1$, the default gives 2.0701 at the third value and adjust=False gives 1.29. Listing 2.32 keeps two numbers per stream and matches pandas in the tests, including on NaN. The two forms agree once $(1 - a)^n$ is negligible.

class Ewma:
    """Same numbers as pandas ewm(alpha=alpha, adjust=True, ignore_na=False)."""

    def __init__(self, alpha):
        self.keep = 1.0 - alpha
        self.num = self.den = 0.0
        self.mean = math.nan

    def update(self, x):
        self.num *= self.keep
        self.den *= self.keep
        if not math.isnan(x):
            self.num += x
            self.den += 1.0
            self.mean = self.num / self.den
        return self.mean


class TimeEwma:
    """Weights fall by half every halflife, whatever the spacing of the ticks."""

    def __init__(self, halflife):
        self.rate = math.log(2.0) / halflife
        self.num = self.den = 0.0
        self.t, self.mean = None, math.nan

    def update(self, t, x):
        if self.t is not None:
            keep = math.exp(-self.rate * (t - self.t))
            self.num *= keep
            self.den *= keep
        self.t = t
        if not math.isnan(x):                 # a NaN ages the weights, adds nothing
            self.num += x
            self.den += 1.0
            self.mean = self.num / self.den
        return self.mean

The default is a true weighted mean, $\sum_i (1 - a)^i x_{n-i} / \sum_i (1 - a)^i$. Here that is $(3 + 0.9 \times 2 + 0.81 \times 1) / 2.71 = 5.61 / 2.71$. Numerator and denominator each have a one-line update, which is the whole trick. adjust=False is the recursion $y_n = (1 - a) y_{n-1} + a x_n$ started at $x_0$, so the first value keeps a weight of $(1 - a)^n$ and not the small share it would get in a long history. In the test the two differ by less than $10^{-6}$ after 200 values.

For a NaN, pandas ages the old weights and adds nothing, which the first class copies. The mean itself does not move at that step, but the next value counts for more.

For uneven times, use a decay of $e^{-\lambda \, \Delta t}$ with $\lambda = \ln 2 / h$. The second class matches ewm called with halflife and times, in the test, NaN included. At the extremes it behaves sensibly. A shared timestamp has $\Delta t = 0$, so the two ticks get equal weight. After a gap much longer than the half-life, the old weights have almost vanished, so the mean jumps to the new tick.

Answer to Question 2.33.

z @ L has covariance $L^\top L$, which for 0.8 is variances 1.64 and 0.36 with covariance 0.48. The correct line is z @ L.T, as in Listing 2.33. The matrix with 0.9, 0.9 and $-0.9$ makes cholesky raise LinAlgError, and so does a correlation of exactly 1.

import numpy as np


def correlated_normals(corr, n, rng):
    lower = np.linalg.cholesky(corr)               # corr = lower @ lower.T
    z = rng.standard_normal((n, len(corr)))        # one row per draw
    return z @ lower.T


def nearest_correlation(corr):
    """Clip negative eigenvalues and rescale: valid, close, not the true nearest."""
    w, v = np.linalg.eigh(corr)
    fixed = v @ np.diag(np.clip(w, 1e-8, None)) @ v.T
    d = np.sqrt(np.diag(fixed))
    return fixed / np.outer(d, d)

Take a row $z$ with covariance $I$. Then $x = z L^\top$ has $\text{Cov}(x) = L\, I\, L^\top = C$. The usual statement $x = L z$ is the same thing with $z$ as a column, which is how the mistake arises. For 0.8, $L$ has rows $(1, 0)$ and $(0.8, 0.6)$. Both versions give a positive correlation, 0.8 against 0.625, and the wrong one also breaks the unit variances, so testing the sample covariance catches it in one line.

If A moves closely with B and with C, then B and C cannot move oppositely. That shows in the eigenvalues, which here are $-0.8$, 1.9 and 1.9. A variance of a portfolio cannot be negative, so no data can give this matrix if all pairs use the same dates. Such matrices arise from pairs estimated over different windows, or from hand edits. So look at the source first. A quick repair clips the negative eigenvalues and rescales the diagonal to 1, which here gives 0.5, 0.5 and $-0.5$. That is a big change, and a warning about the inputs.

At a correlation of 1, the matrix is singular and the factor has a zero pivot. Then use one normal for both variables, or draw from the eigenvectors with $\sqrt{\lambda}$.

Answer to Question 2.34.

Listing 2.34 gives 10.4585 with a standard error of 0.0073, on a million pairs: 1.1 errors from the true value. Treating the two members as independent gives 0.0104 here, too large. On a nearly symmetric payoff, such as a straddle, it is too small, which is the dangerous direction.

import numpy as np


def mc_call(spot, strike, rate, vol, t, pairs, rng, chunk=500_000):
    """Antithetic estimate of a European call: returns (price, standard error)."""
    drift = (rate - 0.5 * vol * vol) * t
    root = vol * np.sqrt(t)
    n, total, squares = 0, 0.0, 0.0
    while n < pairs:                               # memory is bounded by chunk
        m = min(chunk, pairs - n)
        z = rng.standard_normal(m)
        up = np.maximum(spot * np.exp(drift + root * z) - strike, 0.0)
        down = np.maximum(spot * np.exp(drift - root * z) - strike, 0.0)
        one = 0.5 * (up + down)                    # one sample per pair
        total += one.sum()
        squares += (one * one).sum()
        n += m
    mean = total / pairs
    var = max(squares / pairs - mean * mean, 0.0)
    disc = np.exp(-rate * t)
    return disc * mean, disc * np.sqrt(var / pairs)

The pair average is the independent sample. So the variance must come from those averages, over the number of pairs, which is what the listing does. For the call, up and down are negatively correlated: a high $z$ pays and a low $z$ pays nothing. So averaging the pair removes more variance than the naive formula knows, and the naive error is too large. A straddle at the forward pays on both sides, so the members are positively correlated, and the pair carries less information than two draws. The naive error then promises more accuracy than you have, and the test confirms it on such a payoff.

Memory stays bounded because only running totals survive each chunk, so more paths cost time alone. With a path-dependent payoff, the same idea applies to the steps. A full table of a million paths and 252 steps takes about 2 GB, while a loop over steps keeps only the current state.

One exact check needs no option formula: at zero volatility, the price must be $S - K e^{-rT}$ when in the money, with an error of rounding size. A second check is statistical: at strike 0, the call is the stock, so the estimate must lie within a few errors of the spot, which tests the drift.

Answer to Question 2.35.

With orders 9, $-1$ and $L = 8$, the clipped sum gives 8, 8, and the true position is 8, 7. A single order can never show it, so two is the smallest. A cut changes every later position, while cumsum remembers the uncut total. Listing 2.35 has the correct loop.

import numpy as np


def capped_position(trades, limit):
    """Each order is cut to fit; the next one starts from the cut position."""
    pos, out = 0, []
    for q in trades:
        pos = max(-limit, min(limit, pos + q))
        out.append(pos)
    return out


def wrong_clip(trades, limit):
    return np.clip(np.cumsum(trades), -limit, limit).tolist()


def capped_above(trades, limit):
    """A cap on one side only has a closed form."""
    c = np.cumsum(trades)
    return (c - np.maximum(np.maximum.accumulate(c) - limit, 0)).tolist()

The first order is cut from 9 to 8, so one share is never bought. The running sum still counts it, so after $-1$ it says 8 when the position is 7. All it takes is one cut and one reversal.

The two-sided version is path-dependent, since each step needs the last result, so it does not reduce to an array operation in the usual way. Keep the loop, which is clear and tested, and if it shows up on a profile, compile it. That is a fine answer: vectorising is a means, and a wrong fast version is no gain.

With only an upper cap, there is a closed form. Let $c$ be the cumulative sum. All the amount cut so far is $\max(0, \max_{j \le i} c_j - L)$, so the position is $c_i$ minus that, which is capped_above in the listing. The cut total grows only when $c$ sets a new high above $L$. The test compares it with the loop on random orders, and at $L$ larger than all sums it reduces to $c$ itself.

Answer to Question 2.36.

The cache holds them. lru_cache on a method makes self part of the key, and with no bound it keeps every key forever. After 50 requests, cache_info reports 0 hits, 50 misses and a size of 50. Every request has a new self, so a stored result can never be found again.

gc.collect() frees nothing since these are not lost cycles. They are reachable from the function object on the class, which lives as long as the process. In our test a weakref to a finished Pricer is still alive afterwards.

To find growth, compare two tracemalloc snapshots taken some way apart. Here the comparison over 50 requests puts 4.0 MB on the single line that builds curve, and that points to what is kept. It shows where a block was allocated, not who holds it, so gc.get_referrers or a heap tool then answers that part.

It sees only memory requested through Python’s allocators, so a C library calling malloc directly stays invisible. Also, a falling total inside Python need not reduce the process’s resident size, as freed memory can stay with the allocator. So compare the traced total with the resident size: if only the second grows, look below Python.

The fix is to cache a function of plain values at module level, with a bound, so that self is not in the key. Or give each instance its own dictionary, which dies with it. In the test, either version lets the instance go.

maxsize=128 is a bound and not a fix. Memory stops growing at 128 pricers, each still held for nothing, and the hit count stays at 0.

Answer to Question 2.37.

String hashes are randomised per process, so the iteration order of a set of strings changes from run to run. The first line depends on order because all four symbols tie and max keeps the first it meets. The second depends on order because the floating-point sum absorbs the 1.0 or not, depending on when it comes. Integer hashes are not randomised, so codes would not show this.

In the tests, six values of PYTHONHASHSEED give four winners, and a sum of both 0.0 and 1.0. With a fixed seed, the output repeats. An integer set gives one order across all seeds.

The fix is to stop depending on set order wherever order matters: iterate over sorted(symbols), and state a tie-break in the key, such as (score, symbol). A dictionary is safe, since it keeps insertion order by language rule, as long as what fills it comes in a defined order.

Setting PYTHONHASHSEED makes runs repeat, and that hides the fault. The order is still arbitrary, and a different universe, a Python upgrade or one added symbol changes the result with no code change. It is fine as a diagnostic: if fixing it makes the runs agree, you have found the class of bug.

Check next: parallel work collected in order of completion, directory listings read without sorting, threaded linear algebra, whose rounding can differ with the thread count, and a database query with no ORDER BY. Find the first divergence by logging a hash of the state per day in two runs.

Answer to Question 2.38.

Listing 2.36 computes it with group ranks and group sums, and it matches scipy.stats.spearmanr date by date in the test. For jointly normal variables with correlation 0.1, expect $\frac{6}{\pi} \arcsin(0.05) = 0.0955$. One day’s value with 500 stocks has noise near $1 / \sqrt{500} = 0.045$.

def daily_ic(df):
    """df: date, symbol, price, signal. Rank correlation of signal and next return."""
    df = df.sort_values(["symbol", "date"])
    nxt = df.groupby("symbol").price.shift(-1)           # within a symbol only
    df = df.assign(fwd=nxt / df.price - 1).dropna(subset=["signal", "fwd"])
    r = df.groupby("date")[["signal", "fwd"]].rank()     # ranks among complete rows
    r = r - r.groupby(df["date"]).transform("mean")
    num = (r["signal"] * r["fwd"]).groupby(df["date"]).sum()
    sxx = (r["signal"] ** 2).groupby(df["date"]).sum()
    syy = (r["fwd"] ** 2).groupby(df["date"]).sum()
    return num / (sxx * syy) ** 0.5

A rank correlation is a Pearson correlation of ranks, so rank within each date, remove each date’s mean, and form the sum of products over the root of the product of the two sums of squares. Ranking comes after the dropna, so that both columns rank the same rows, for the same reason as in any pairwise statistic.

df.price.shift(-1) on a table sorted by symbol gives the last row of each symbol the first price of the next symbol. That makes a return of something like $-40\%$, with no NaN to warn you, on one row per symbol. A group shift gives NaN there, and the dropna removes it.

A simulation of 500 stocks and 199 days gives a mean of 0.0955 and a daily spread of 0.041, near that estimate, for a t-statistic of 33, namely the mean divided by the spread over $\sqrt{199}$, that is $0.0955 / (0.041 / \sqrt{199})$. So a signal that is small on any one day is plain over a long sample. At the extremes, a signal equal to the return gives exactly 1, its negative gives $-1$, and a constant signal gives NaN.

With a five-day target, neighbouring days share four of their five returns, so the daily values overlap and the plain t-statistic is too high, by about $\sqrt{5}$ if the signal changes slowly. So sample every fifth day, or use an error that allows for that overlap.

Answer to Question 2.39.

Listing 2.37 replaces times by their ranks across both arrays, builds one integer key per symbol and rank, then does a single binary search per trade. It matches pd.merge_asof and a slow loop on random data. A multiplier large enough for the times overflows int64, whose limit is about $9.2 \times 10^{18}$.

import numpy as np


def asof(trade_sym, trade_t, quote_sym, quote_t):
    """Row of the latest quote at or before each trade, same symbol; -1 if none."""
    if len(quote_t) == 0:
        return np.full(len(trade_t), -1)
    times, ranks = np.unique(np.concatenate([trade_t, quote_t]), return_inverse=True)
    span = len(times)                      # int64 codes: sym * span cannot overflow
    trade_key = np.asarray(trade_sym, dtype=np.int64) * span + ranks[: len(trade_t)]
    quote_key = np.asarray(quote_sym, dtype=np.int64) * span + ranks[len(trade_t) :]
    order = np.argsort(quote_key, kind="stable")     # equal keys keep arrival order
    i = np.searchsorted(quote_key[order], trade_key, side="right") - 1
    found = order[np.maximum(i, 0)]
    return np.where((i >= 0) & (quote_sym[found] == trade_sym), found, -1)

The idea is to reduce two sort keys to one. If the key orders by symbol and then time, a binary search for a trade’s key finds the last quote at or below it, and only one check remains: that this quote has the same symbol. If it has another symbol, then this symbol has no quote at or before the trade, since such a quote would have a key between the two.

With $10^{19}$ itself, NumPy 2.2 raises OverflowError, since that does not fit in int64. A multiplier that does fit must still be above $1.8 \times 10^{18}$ for the times, so symbol 5 plus a time, about $1.08 \times 10^{19}$, already passes the limit, and NumPy wraps with no error, making the order wrong. np.unique with return_inverse fixes that. There are at most $m + n$ distinct times, so the key stays below $5000 (m + n)$, and the ranks keep the order, which is all the search needs. The stem gives no type for the codes, so the listing casts them to int64 first: int32 codes would wrap with no error once $m + n$ reaches about a million, and int16 codes raise.

side='right' finds quotes at or before the trade, which is the stem’s rule. With left it would be strictly before. The stable argsort keeps arrival order among equal keys, so the search lands on the last to arrive.

The sorts dominate, at $O((m + n) \log (m + n))$, and the searches take $O(m \log n)$. If both inputs arrive sorted, then a merge walk per symbol is linear and it also needs less memory. Here the rank step needs an array of $m + n$ times, and one of ranks.

Answer to Question 2.40.

The cumulative-sum method gives a large negative number at row 999, $-77$ on our draw, where the true slope is 2.00. Listing 2.38 avoids this in $O(n)$, because each window’s sum uses only values inside that window. Pandas 2.3.3’s rolling variance shows a whole multiple of $128/49 = 2.61$ there, $-2.6$ on our draw and positive or zero on others, where the true value is about $10^{-4}$. Its rolling.corr shows NaN, a tiny number, or a value above 1.

import numpy as np


def window_sum(v, w):
    """Sum of the last w values at each row, using only values inside that window."""
    n = len(v)
    a = np.concatenate([v, np.zeros(-n % w)]).reshape(-1, w)
    prefix = np.cumsum(a, axis=1).ravel()
    suffix = np.cumsum(a[:, ::-1], axis=1)[:, ::-1].ravel()
    out = prefix[:n].copy()                     # rows up to w - 1, and ends of blocks
    i = np.arange(w, n)
    i = i[i % w != w - 1]                       # a window spans exactly two blocks
    out[i] = suffix[i - w + 1] + prefix[i]
    return out


def rolling_beta(x, y, w, min_obs):
    """Beta of y on x over the last w rows; a row counts only if both are present."""
    ok = ~(np.isnan(x) | np.isnan(y))
    x0, y0 = np.where(ok, x, 0.0), np.where(ok, y, 0.0)
    n = window_sum(ok.astype(float), w)
    sx, sy = window_sum(x0, w), window_sum(y0, w)
    sxx, sxy = window_sum(x0 * x0, w), window_sum(x0 * y0, w)
    with np.errstate(divide="ignore", invalid="ignore"):
        beta = (n * sxy - sx * sy) / (n * sxx - sx * sx)
    return np.where(n >= min_obs, beta, np.nan)

Masking both columns first makes all five window sums cover the same rows, and the slope is $(n S_{xy} - S_x S_y) / (n S_{xx} - S_x^2)$.

The damage comes from $x^2 = 10^{18}$. Doubles near that size are 128 apart, so once it enters the running total, every later small term is rounded away. Subtracting two such totals gives exactly 0, while the true $S_{xx}$ is near $5 \times 10^{-3}$. The denominator is then $-S_x^2$, which makes the slope large and negative: from $-16$ to $-44\,000$ over ten draws. The fault stays in every later row, long after the bad tick has left the window.

The listing cuts the series into blocks of $w$. Any window of length $w$ is then a suffix of one block plus a prefix of the next, and both are plain cumulative sums within a block. So a sum never contains a value from outside its window, and each row costs one addition. In the test it matches a slow fit on every window clear of the bad tick. Windows that contain the tick are still dominated by it, as arithmetic demands, so cleaning remains a separate step.

The pandas figures are what version 2.3.3 does here. Its rolling variance adds and removes values as the window moves, and it leaves a residue of that step on this input, while rolling.cov on these columns is still right. Run a check of this sort on your own library version, because the behaviour may change.