quant/src/arbitrage.rs
The geometry behind the no-arbitrage chapter's proof of the fundamental theorem.
//! The geometry behind the no-arbitrage chapter's proof of the fundamental theorem.//!//! One period, two states, one risky asset. Small enough to draw, which is the//! point: the separation argument is easy to follow and hard to picture, and in//! two dimensions the picture is a line and a segment.//!//! Discounted, the payoffs reachable from nothing form the set//!//! ```text//! K = { theta (S1 - S0) : theta real }//! ```//!//! which is a line through the origin in the plane whose axes are the two//! states. Arbitrage is a point of `K` with both coordinates non-negative and at//! least one positive — so, after scaling, a point of `K` on the unit simplex.//! No arbitrage is exactly the statement that the line misses the segment. /// A one-period, two-state market, in discounted prices.#[derive(Clone, Copy, Debug)]pub struct TwoState { /// Today's discounted price of the risky asset. pub spot: f64, /// Its discounted price if the up state happens. pub up: f64, /// And if the down state happens. pub down: f64,} impl TwoState { /// The direction of the line of attainable payoffs. /// /// Holding one unit of the asset, funded from the money market, pays /// `S1 - S0` in each state. Every other zero-cost portfolio is a multiple. pub fn direction(&self) -> (f64, f64) { (self.up - self.spot, self.down - self.spot) } /// Whether the market admits arbitrage. /// /// The line meets the positive quadrant precisely when its direction, or its /// negative, has both coordinates non-negative — that is, when the spot sits /// outside the range of tomorrow's prices. Buy the asset when it cannot fall, /// sell it when it cannot rise. pub fn has_arbitrage(&self) -> bool { let (a, b) = self.direction(); let same_way_up = a >= 0.0 && b >= 0.0 && (a > 0.0 || b > 0.0); let same_way_down = a <= 0.0 && b <= 0.0 && (a < 0.0 || b < 0.0); same_way_up || same_way_down } /// The risk neutral probabilities, when they exist. /// /// The no-arbitrage chapter obtains these as a vector orthogonal to `K` with strictly /// positive entries, normalised to sum to one. In two dimensions the /// orthogonal complement of a line is a line, so there is exactly one /// candidate and the only question is whether it is positive — which is the /// no-arbitrage condition again, now visible as a sign. pub fn risk_neutral(&self) -> Option<(f64, f64)> { // Orthogonal to the direction (a, b), so proportional to (-b, a) — and // the pair that prices the asset correctly is (spot - down, up - spot) // over the range. let span = self.up - self.down; if span.abs() < 1e-12 { return None; } let q_up = (self.spot - self.down) / span; let q_down = 1.0 - q_up; if q_up > 0.0 && q_down > 0.0 { Some((q_up, q_down)) } else { None } } /// The line of attainable payoffs, as a function of the up-state payoff. /// /// Returned as a function so the figure can plot it against the simplex on a /// shared axis. Vertical lines — which happen when the asset cannot move in /// the up state — have no such representation and come back as `None`. pub fn payoff_line(&self, x_up: f64) -> Option<f64> { let (a, b) = self.direction(); if a.abs() < 1e-12 { return None; } Some(b / a * x_up) }} #[cfg(test)]mod tests { use super::*; fn market(spot: f64) -> TwoState { TwoState { spot, up: 120.0, down: 90.0 } } #[test] fn arbitrage_appears_exactly_when_the_spot_leaves_the_range() { assert!(market(85.0).has_arbitrage(), "a spot below every outcome is free money"); assert!(market(125.0).has_arbitrage(), "and above every outcome, sold short"); for spot in [90.5, 100.0, 105.0, 119.5] { assert!(!market(spot).has_arbitrage(), "spot {spot} should be arbitrage-free"); } } #[test] fn the_risk_neutral_measure_exists_exactly_when_there_is_no_arbitrage() { // The fundamental theorem, in the smallest case where it has content. for spot in [80.0, 89.0, 90.0, 95.0, 110.0, 120.0, 130.0] { let m = market(spot); assert_eq!( m.risk_neutral().is_some(), !m.has_arbitrage(), "the two disagreed at spot {spot}" ); } } #[test] fn the_risk_neutral_measure_reprices_the_asset() { let m = market(102.0); let (q_up, q_down) = m.risk_neutral().unwrap(); let priced = q_up * m.up + q_down * m.down; assert!((priced - m.spot).abs() < 1e-12, "priced at {priced}, spot is {}", m.spot); } #[test] fn the_measure_is_orthogonal_to_every_attainable_payoff() { // Which is what "separating" meant: the plane through the origin with // normal q contains the whole of K and misses the simplex. let m = market(102.0); let (q_up, q_down) = m.risk_neutral().unwrap(); let (a, b) = m.direction(); assert!((q_up * a + q_down * b).abs() < 1e-12); }} // ---------------------------------------------------------------------------// The separating lemma itself, rather than the theorem that uses it.// --------------------------------------------------------------------------- /// The construction inside Lemma "Separating a subspace from the simplex",/// drawn in the plane.////// The lemma's proof does not separate `K` from `C` directly. It forms the/// difference set `D = C - K`, observes that `K` missing `C` is exactly/// `0` not in `D`, and takes `lambda` to be the point of `D` closest to the/// origin. Every conclusion then falls out of that one choice.////// In two dimensions the whole thing can be seen at once. `C` is the segment/// from `(1,0)` to `(0,1)`; `K` is a line through the origin; and `D`, being a/// segment swept along a line, is an infinite strip parallel to `K`. So:////// * `lambda` is perpendicular to `K` because the nearest point of a strip is/// reached by going straight at it, and the strip runs parallel to `K` ---/// which is the lemma's conclusion `lambda . k = 0`;/// * the supporting hyperplane at `lambda` is not a separate object, it is the/// near edge of the strip;/// * `lambda` has positive coordinates because the strip lies on the far side/// of the origin from them, which is the conclusion that makes it a measure.pub struct Separation { /// Unit vector spanning `K`. pub direction: (f64, f64), /// Unit normal to `K`. The strip's edges are offsets along this. pub normal: (f64, f64), /// Whether `K` misses the simplex, which is the lemma's hypothesis. pub separated: bool, /// The point of `D` nearest the origin. The zero vector when the hypothesis /// fails, since then the origin is itself in `D`. pub lambda: (f64, f64), /// `lambda` scaled to sum to one: the risk neutral measure the lemma /// delivers. `None` when the hypothesis fails. pub measure: Option<(f64, f64)>, /// Offsets along `normal` of the strip's two edges, nearest first. pub near_offset: f64, pub far_offset: f64,} /// Build the picture for a subspace `K` at angle `theta` to the first axis.pub fn separation(theta: f64) -> Separation { let direction = (theta.cos(), theta.sin()); let normal = (-theta.sin(), theta.cos()); // Every point of the strip has the same range of offsets along the normal as // the simplex does, because sliding along K does not change that component. // The simplex is the segment between the two axis vectors, so the range is // the interval between their offsets. let at_e1 = normal.0; let at_e2 = normal.1; let (lo, hi) = (at_e1.min(at_e2), at_e1.max(at_e2)); // The origin is inside the strip exactly when zero is inside that interval, // which is exactly when K meets the simplex. let separated = lo > 0.0 || hi < 0.0; let (near, far) = if lo > 0.0 { (lo, hi) } else { (hi, lo) }; let offset = if separated { near } else { 0.0 }; let lambda = (offset * normal.0, offset * normal.1); let total = lambda.0 + lambda.1; let measure = if separated && total > 0.0 { Some((lambda.0 / total, lambda.1 / total)) } else { None }; Separation { direction, normal, separated, lambda, measure, near_offset: if separated { near } else { 0.0 }, far_offset: far, }} #[cfg(test)]mod separation_tests { use super::*; /// Brute force: the nearest point of the strip to the origin, found by /// scanning rather than reasoned about. Shares no algebra with `separation`. fn nearest_by_search(theta: f64) -> (f64, f64) { let (u, _n) = ((theta.cos(), theta.sin()), (-theta.sin(), theta.cos())); let mut best = (f64::INFINITY, 0.0, 0.0); for i in 0..=2000 { // A point of the simplex. let s = i as f64 / 2000.0; let c = (1.0 - s, s); for j in -4000..=4000 { let t = j as f64 * 0.005; let z = (c.0 + t * u.0, c.1 + t * u.1); let d2 = z.0 * z.0 + z.1 * z.1; if d2 < best.0 { best = (d2, z.0, z.1); } } } (best.1, best.2) } #[test] fn the_hypothesis_holds_exactly_when_the_line_misses_the_simplex() { // K is a line through the origin and the simplex sits in the positive // quadrant, so the line meets it precisely when its angle points into // that quadrant. for degrees in [-80.0, -30.0, 100.0, 135.0, 179.0] { assert!(separation(degrees * std::f64::consts::PI / 180.0).separated, "expected separation at {degrees} degrees"); } for degrees in [1.0, 30.0, 45.0, 60.0, 89.0] { assert!(!separation(degrees * std::f64::consts::PI / 180.0).separated, "expected no separation at {degrees} degrees"); } } #[test] fn lambda_is_perpendicular_to_the_subspace() { // The lemma's second conclusion, and the reason the strip picture works: // the nearest point of a strip is reached at right angles to it. for degrees in [-80.0, -30.0, 100.0, 135.0, 170.0] { let s = separation(degrees * std::f64::consts::PI / 180.0); let dot = s.lambda.0 * s.direction.0 + s.lambda.1 * s.direction.1; assert!(dot.abs() < 1e-12, "at {degrees} degrees lambda.k was {dot}"); } } #[test] fn lambda_has_strictly_positive_coordinates() { // The conclusion that turns the lemma into a pricing measure. Without it // there would be a separating vector but no equivalent measure. for degrees in [-89.0, -45.0, 91.0, 135.0, 179.0] { let s = separation(degrees * std::f64::consts::PI / 180.0); assert!(s.separated); assert!(s.lambda.0 > 0.0 && s.lambda.1 > 0.0, "at {degrees} degrees lambda was {:?}", s.lambda); let m = s.measure.expect("a measure"); assert!((m.0 + m.1 - 1.0).abs() < 1e-12); assert!(m.0 > 0.0 && m.1 > 0.0); } } #[test] fn lambda_really_is_the_nearest_point_of_the_difference_set() { // The claim the proof turns on, checked against a search that knows // nothing about normals, offsets or which case it is in. The angles // cover both: where the line misses the simplex the nearest point is // lambda, and where it meets it the nearest point is the origin itself, // which `separation` reports as the zero vector. // // The tolerance is the search grid's own resolution rather than a fudge: // it steps along the strip in units of 0.005, so it cannot locate any // point more precisely than that. for degrees in [-60.0, 45.0, 89.0, 120.0, 150.0] { let theta = degrees * std::f64::consts::PI / 180.0; let s = separation(theta); let (x, y) = nearest_by_search(theta); assert!( (s.lambda.0 - x).hypot(s.lambda.1 - y) < 6e-3, "at {degrees} degrees: {:?} against search ({x}, {y})", s.lambda ); // And the two agree on which case it is. assert_eq!( s.separated, x.hypot(y) > 1e-2, "at {degrees} degrees the two disagreed about separation" ); } }} /// How deep a doubling strategy must be prepared to go before it wins.////// The no-arbitrage chapter's doubling strategy bets on a Brownian motion until/// its winnings first reach one. That happens almost surely, which is why the/// naive definition of arbitrage fails in continuous time and why admissibility/// has to require the wealth to stay above a fixed level.////// This measures what that requirement excludes. Run the path until it reaches/// either `+1` or `-depth`; the gambler's ruin identity for Brownian motion says////// ```text/// P(reach -a before +1) = 1 / (1 + a)./// ```////// Returns the fraction of paths that hit the lower barrier first, and the mean/// drawdown *capped at* `depth` --- which by the same identity is////// ```text/// E[min(-M, a)] = integral_0^a P(-M > u) du = ln(1 + a),/// ```////// growing logarithmically and without limit. That is the divergence made/// simulable: every path terminates at one barrier or the other, so nothing is/// censored, and the growth in `a` is the thing being measured.////// # Why the expected drawdown is not computed here////// It is infinite --- a tail decaying like `1/a` has no mean --- and that is/// exactly what cannot be established by simulation. Any run has to stop/// somewhere, and stopping censors the tail that carries the divergence, so the/// sample mean converges to a finite number and converges to the wrong one. The/// divergence follows from the barrier law above by an integral, and that is/// where it should be established; the simulation's job is to confirm the law.pub fn ruin_probability( depth: f64, paths: usize, steps_per_unit: usize, seed: u64,) -> (f64, f64) { use crate::pathwise::Rng; let mut rng = Rng::new(seed); let dt = 1.0 / steps_per_unit as f64; let root_dt = dt.sqrt(); let (mut ruined, mut capped_total) = (0usize, 0.0); for _ in 0..paths { let (mut w, mut lowest) = (0.0f64, 0.0f64); loop { w += root_dt * rng.next_normal(); lowest = lowest.min(w); if w <= -depth { ruined += 1; capped_total += depth; break; } if w >= 1.0 { capped_total += -lowest; break; } } } (ruined as f64 / paths as f64, capped_total / paths as f64)} #[cfg(test)]mod doubling_tests { use super::*; #[test] fn the_drawdown_has_the_gamblers_ruin_distribution() { // The exact law: the chance of falling to -a before rising to +1 is // 1/(1+a). Two barriers, so every path terminates and nothing is // censored. // // The grid has to be fine, and for a reason worth knowing: a discretised // path overshoots a barrier it crosses, by about half a step on average, // which is the same as moving the barrier outwards. That biases the ruin // probability upwards, and the bias grows with the depth because the // overshoot is a larger fraction of the near barrier's distance. At // sixteen hundred steps a unit it is under one per cent out to a depth of // three and about two and a half per cent by nine. for depth in [0.5, 1.0, 3.0] { let (measured, _) = ruin_probability(depth, 40_000, 1600, 20260808); let exact = 1.0 / (1.0 + depth); assert!( (measured / exact - 1.0).abs() < 0.015, "depth={depth}: {measured:.5} against 1/(1+a) = {exact:.5}" ); } } #[test] fn the_capital_required_grows_without_limit() { // The consequence, and the honest way to see it. The expected drawdown is // infinite, which no simulation can show: any run stops somewhere, and // stopping censors the tail carrying the divergence, so a sample mean // converges to a finite number and converges to the wrong one. // // What is simulable is the divergence's rate. Capping the drawdown at a // makes every path terminate at one barrier or the other, and // // E[min(-M, a)] = integral_0^a du/(1+u) = ln(1 + a), // // which grows without limit as the cap is lifted. The float a doubling // strategy needs to survive a drawdown of a is therefore about ln(1+a), // and no finite float suffices -- which is what admissibility excludes, // and it excludes it with a bound of any size whatever. // // The capped mean tolerates a coarser grid than the ruin probability // does, since an overshoot at the lower barrier changes the recorded // drawdown by a fraction of a step rather than reclassifying the path. let mut previous = 0.0; for depth in [1.0, 4.0, 15.0] { let (_, mean) = ruin_probability(depth, 40_000, 400, 20260808); let exact = (1.0 + depth).ln(); assert!( (mean / exact - 1.0).abs() < 0.025, "depth={depth}: mean capped drawdown {mean:.4} against ln(1+a) = {exact:.4}" ); assert!(mean > previous, "lifting the cap must raise the mean"); previous = mean; } }} /// Whether the doubling strategy's time compression can actually be carried out.////// The no-arbitrage chapter needs a strategy on `[0, 1)` whose wealth reaches/// one almost surely, and the first passage of a Brownian motion to a level is/// almost surely finite but has no bound. The construction that reconciles them/// is a deterministic position,////// ```text/// Delta_t = 1 / (1 - t),/// ```////// whose wealth `pi_t = integral_0^t dW / (1 - s)` is a continuous local/// martingale with quadratic variation////// ```text/// <pi>_t = integral_0^t ds / (1 - s)^2 = t / (1 - t),/// ```////// which is unbounded on `[0, 1)`. By Dambis-Dubins-Schwarz `pi_t = B_{<pi>_t}`/// for some Brownian motion `B`, so the whole of `B` on `[0, infinity)` is/// traversed before calendar time one, and `B` reaches any level.////// That chain of reasoning is what this checks, and it checks it where it could/// fail. Simulate `pi` on a grid in the clock `u = t / (1 - t)`, where the steps/// are of equal variance and the exploding integrand is resolved rather than/// stepped over, and count the paths reaching one by clock time `U`. If the/// identification with a Brownian motion is right, that fraction is the/// reflection principle's////// ```text/// P(max_{u <= U} B_u >= 1) = 2 (1 - Phi(1 / sqrt(U))),/// ```////// which tends to one as `U` grows --- and every one of those clock times is a/// calendar time `u / (1 + u)`, strictly inside `[0, 1)`.////// # Monitoring on a grid moves the barrier////// A path is only inspected at grid points, so a crossing that happens and/// reverses within a step is missed. The bias is downwards and it is not small/// enough to ignore: the expected overshoot of a discretely monitored maximum/// puts the effective barrier at `1 + 0.5826 sqrt(du)`, the Broadie-Glasserman-Kou/// shift, which at sixteen hundred steps a clock unit is a barrier one and a half/// per cent too far away and a probability one per cent too low. Both the/// continuous law and the law with the shifted barrier are returned, since the/// first is what the chapter states and the second is what a grid can measure.pub struct Compression { /// Fraction of paths whose wealth reached one by clock time `U`. pub reached: f64, /// The continuous-monitoring law, `2 (1 - Phi(1 / sqrt(U)))`. pub exact: f64, /// The same law with the barrier moved out by the discretisation shift. pub exact_monitored: f64, /// Mean calendar time at which the winners get there. Below one, which is /// the whole point. pub mean_calendar_time: f64,} /// Run the compressed-clock strategy. See [`Compression`].pub fn compressed_clock( horizon: f64, paths: usize, steps_per_unit: usize, seed: u64,) -> Compression { use crate::black::norm_cdf; use crate::pathwise::Rng; let mut rng = Rng::new(seed); // Fixed step in clock time, so the monitoring is equally fine at every // horizon and the bias does not grow with it. let du = 1.0 / steps_per_unit as f64; let steps = (horizon * steps_per_unit as f64).ceil() as usize; let root_du = du.sqrt(); let (mut hits, mut calendar_total) = (0usize, 0.0); for _ in 0..paths { let mut pi = 0.0f64; for i in 0..steps { pi += root_du * rng.next_normal(); if pi >= 1.0 { // Clock time u maps back to calendar time u / (1 + u), the // inverse of the quadratic variation above. This is where the // compression is: an unbounded clock inside a bounded calendar. let u = du * (i + 1) as f64; hits += 1; calendar_total += u / (1.0 + u); break; } } } let root_u = horizon.sqrt(); Compression { reached: hits as f64 / paths as f64, exact: 2.0 * (1.0 - norm_cdf(1.0 / root_u)), exact_monitored: 2.0 * (1.0 - norm_cdf((1.0 + 0.5826 * root_du) / root_u)), mean_calendar_time: if hits == 0 { f64::NAN } else { calendar_total / hits as f64 }, }} #[cfg(test)]mod compression_tests { use super::*; #[test] fn the_clock_can_be_compressed_into_the_unit_interval() { // The construction is only an arbitrage if the target is reached before // calendar time one, so there are two things to check: that the target // is reached with the probability the reflection principle gives, which // is what confirms the wealth really is a Brownian motion on the // stretched clock, and that the calendar time it takes stays inside the // interval. for horizon in [1.0, 4.0, 25.0] { let c = compressed_clock(horizon, 20_000, 1600, 20260810); assert!( (c.reached - c.exact_monitored).abs() < 0.006, "U={horizon}: {:.4} against {:.4}", c.reached, c.exact_monitored ); // The shift is the only reason the two differ, and it differs one way. assert!( c.reached < c.exact, "U={horizon}: grid monitoring cannot exceed the continuous law" ); assert!( c.mean_calendar_time < horizon / (1.0 + horizon), "U={horizon}: mean calendar hit {:.6} is not inside the interval", c.mean_calendar_time ); } } #[test] fn lifting_the_clock_horizon_drives_the_probability_to_one() { // The almost-sure claim, seen as a limit. A simulation cannot exhibit an // almost-sure event; what it can do is confirm the law whose limit is // one, at horizons where the grid is fine enough to measure it. let mut previous = 0.0; for horizon in [1.0, 4.0, 25.0] { let c = compressed_clock(horizon, 8_000, 1600, 20260810); assert!( (c.reached - c.exact_monitored).abs() < 0.012, "U={horizon}: {:.4} against {:.4}", c.reached, c.exact_monitored ); assert!( c.reached > previous, "U={horizon}: {:.4} did not exceed {previous:.4}", c.reached ); previous = c.reached; } assert!(previous > 0.83, "at U=25 the target is already usual, got {previous:.4}"); }}