μ

plotting stories in μ

Plotting Stories

Each card is a self-contained narrative: a mathematical object is constructed, mapped, or simulated, then rendered as an SVG chart. μ provides two plotting APIs: the original verb-based system (plot, plothist, plothline, heatmap) for quick single-purpose charts, and the new layered chart system (chart + layer constructors) for composing multiple visual layers — line, scatter, band, bar, hline, vline, diag, step — into a single chart with shared axes. Both systems work in the CLI (saving SVG to the current directory) and in the playground (rendering to canvas).

Basic

Simple Plot

plot 1 4 9 16 25
// → renders SVG plot

Monadic plot takes a vector y ∈ ℝⁿ and synthesises the abscissa as {0, 1, …, n−1}, rendering a polyline over (i, yᵢ). This is the fastest path when only the function values matter and domain scaling is irrelevant — for instance when visualising index sequences, filter outputs, or sorted statistics. When an explicit, non-uniform domain is required (log-strike grids, irregular time series), use the dyadic form plot[x; y] shown below.

X vs Y Plot

plot[1 2 3 4 5;1 4 9 16 25]
// → renders SVG plot with custom x-axis

Dyadic plot[x; y] maps the ordered pairs (xᵢ, yᵢ) to a polyline, requiring |x| = |y|. The semicolon inside brackets is μ's vector-argument separator, equivalent to passing a 2-tuple without the overhead of constructing an intermediate list. This matters whenever the domain is non-uniform: log-moneyness grids for volatility surfaces, irregular time stamps, or parameter sweeps with non-constant step size. An optional symbol third argument labels the chart title: plot[x; y; `MyChart].

Multi-Series Line Chart

// Time points (shared across all series)
t: 0.0 0.25 0.5 0.75 1.0

// Three different convergence paths
linear:     0.0 0.25 0.50 0.75 1.00
exponential:1.0 0.78 0.61 0.47 0.37
constant:   0.5 0.50 0.50 0.50 0.50

// Combine into list of lists using cut
paths: cut 5 (linear, exponential, constant)

// Multi-series line chart with legend
plot[t; paths; `Convergence]
// → three colored lines with legend showing 0, 1, 2

Multi-series plot renders multiple line series when the y argument is a list of lists — each inner list becomes a series, plotted with distinct colors from the palette. Up to 8 series are supported. A legend appears at top-right showing series indices (0, 1, 2…). This is ideal for comparing trajectories: "convergence rates," "sensitivity scenarios," or "model vs actual" over a shared x-axis (time, strike price, etc.). Use cut n list to reshape flat data into a list-of-lists: cut 5 (series1, series2, series3) creates 3 series of 5 points each. The x-axis values are shared across all series.

Multi-Series Grouped Bar Chart

// Exit multiples: 3x, 4x, 5x, 6x, 7x, 8x
exit_mults: 3.0 4.0 5.0 6.0 7.0 8.0

// P&L without contract protection (sensitive to full exit value)
without_contract: 15.0 20.0 25.0 30.0 35.0 40.0

// P&L with contract protection (capped at certain levels)
with_contract: 12.0 18.0 24.0 30.0 36.0 42.0

// Combine into list of lists using cut
series_data: cut 6 (without_contract, with_contract)

// Grouped bar chart comparing both scenarios
plothist[exit_mults; series_data; `ExitValueComparison]
// → side-by-side bars at each exit multiple, legend shows 0 and 1

Multi-series plothist renders grouped (side-by-side) bars when the counts argument is a list of lists — each inner list becomes a series, with bars grouped at each bin position. The bin positions (first argument) are shared across all series. Up to 8 series are supported, each rendered with a distinct color from the palette. A legend appears at top-right showing series indices (0, 1, 2…). This is ideal for comparing scenarios: "with vs without contract protection," "actual vs target," or "current vs simulated" at each exit multiple, strike level, or time bucket. The cut n list verb is the standard way to create a list-of-lists from flat data: cut 6 (series1, series2) reshapes 12 elements into 2 series of 6 values each. All series must have equal length matching the bins vector.

Histogram with Custom Bins

\S 42
// Generate 1000 normal random variates
bm:{[n] sqrt(-2*log(1.0-?n)) * cos(2*acos(-1)*?n)}
samples: bm 1000

// Bin into 20 buckets and plot distribution
bucket_counts: hist[20; samples]
plothist[bucket_counts; `NormalDist]
// → bell-shaped bar chart centered at 0

The hist[n; data] verb bins numeric data into n equal-width buckets spanning the data range [min, max], returning a count vector. When passed to plothist, these counts render as a histogram bar chart. This is the standard workflow for visualizing empirical distributions: bin with hist, then plot with plothist. For custom bin centers (e.g., log-spaced strikes), pass explicit bin positions as the first argument and your pre-computed counts as the second.

Layered Chart Composition

The chart Verb

The chart verb composes multiple visual layers into a single chart with shared axes, automatic range computation, and palette color assignment. Its general form is:

chart[layer1; layer2; ...; `Title]

Each layer is constructed with a layer constructor verb that returns a tagged list. The chart verb parses these tagged lists, computes shared axis extents, and renders all layers to SVG (CLI) or to the Canvas binary protocol (playground). Up to 16 layers are supported.

Layer Constructors

// Layer constructors return tagged lists:
line[0 1 2 3; 1 4 9 16]
// → (`line; 0 1 2 3; 1 4 9 16)

scatter[0 1 2; 3 6 9]
// → (`scatter; 0 1 2; 3 6 9)

hline 2.5
// → (`hline; 2.5)

diag[]
// → (`diag)

Layer constructors are composable values — they return tagged lists that can be stored in variables, passed to functions, or composed programmatically before being passed to chart. This decouples data preparation from rendering.

Line + Confidence Band

// x-axis: 5 time steps
x: 0 1 2 3 4

// Mean and ± 1 standard deviation
mean: 10 12 11 14 13
lo:    8 10  9 11 10
hi:   12 14 13 17 16

// Compose: band for uncertainty, line for mean
chart[band[x; lo; hi]; line[x; mean]; `Forecast]
// → shaded band with line overlay

The band layer renders a semi-transparent filled polygon between ylo and yhi, ideal for confidence intervals, prediction bands, or Bollinger bands. Layers are drawn in order, so place band before line to keep the mean visible above the fill.

Scatter + Horizontal Reference

x: 1 2 3 4 5 6 7 8
y: 2.1 3.5 2.8 4.9 3.2 5.1 4.4 6.0

chart[scatter[x; y]; hline 4.0; `Threshold]
// → circles with dashed line at y=4

Combining scatter with hline is useful for visualizing threshold decisions, breakeven levels, or target values. The horizontal line is rendered with a dashed stroke to distinguish it from data layers.

QQ Plot (Normal Probability)

\S 42
bm:{[n] sqrt(-2*log(1.0-?n)) * cos(2*acos(-1)*?n)}

// Generate 100 normal samples
data: bm 100

// QQ plot: scatter of (theoretical quantile, sorted sample)
// plus y=x reference diagonal
qqplot[data; `NormalityCheck]
// → points hug the diagonal if data is normal

The qqplot verb is a convenience that computes theoretical normal quantiles Φ⁻¹((i + 0.5)/n) for each sorted sample value and renders a scatter plot against the y = x diagonal. If the sample is truly normal, the points lie on the diagonal; systematic curvature indicates skewness or heavy tails. This uses the Winitzki approximation for the inverse error function, accurate to ~10⁻⁵. The verb accepts a title as an optional second argument: qqplot[data; `Title].

Step Function (Empirical CDF)

// Sorted sample quantiles
x: 0 1 2 3 4 5

// Empirical CDF values F(x) = (rank - 0.5) / n
y: 0.0 0.17 0.33 0.5 0.67 0.83

chart[step[x; y]; hline 0.5; `EmpiricalCDF]
// → staircase with median reference line

The step layer renders a piecewise-constant function (staircase) by drawing horizontal then vertical segments between consecutive points. This is the standard visual representation of empirical CDFs, survival functions, and binned rates. Combined with hline 0.5, it immediately shows the sample median.

Multi-Layer Composition

x: 0 1 2 3 4 5 6 7 8 9

// Two time series
actual:   3 5 4 7 6 8 9 7 10 11
forecast: 2 4 5 6 7 7 8 9 9 10

// Prediction interval
lo: 1 3 3 4 5 5 6 7 7 8
hi: 3 5 7 8 9 9 10 11 11 12

chart[band[x;lo;hi]; line[x;actual]; line[x;forecast]; hline 7; `ModelVsActual]
// → band + two lines + threshold: 4 layers

The layered system scales naturally: each layer gets an auto-assigned palette color (up to 8 distinct colors), and the chart computes axis extents across all layers simultaneously. This makes it trivial to overlay actuals, forecasts, confidence bands, and reference levels in a single composable expression. The chart title (a backtick-symbol) is always the last argument.

Finance — Plotting Stories

Black-Scholes: Call Price vs Strike

// --- normal CDF and bscall from the Black-Scholes example ---
erf0:{[t] p:1.0%(1.0+0.3275911*t);
  1.0-exp(-t*t)*(
    (0.254829592*p)+((-0.284496736)*p*p)
    +(1.421413741*p*p*p)
    +((-1.453152027)*p*p*p*p)
    +(1.061405429*p*p*p*p*p))}
sq2: sqrt 2.0
sign:{1-2*(0>x)}
ncdf:{[v] 0.5+0.5*(sign v)*(erf0((abs v)%sq2))}
bscall:{[S;K;r;sig;T]
  (S*(ncdf (log(S%K)+(r+0.5*sig*sig)*T)%(sig*sqrt T)))
  -(K*exp(-r*T)*(ncdf ((log(S%K)+(r+0.5*sig*sig)*T)%(sig*sqrt T))-sig*sqrt T))}

// build strike grid 60..140 in steps of 2
n:41
ks: 60.0 + 2.0 * !n
// map bscall over strikes (S=100, r=0.05, sig=0.2, T=1)
cs:(each[{bscall 100 x 0.05 0.2 1.0};ks])
plot[ks;cs;`BSCall]
// → smooth call-price curve, ~39 at K=60 down to ~3 at K=140

This traces the Black-Scholes call price C(K) = S Φ(d₁) − K e^−rT Φ(d₂) over a uniform strike grid K ∈ {60, 62, …, 140} with S = 100, r = 5%, σ = 20%, T = 1. each[f; list] is a first-class map: it applies the partially-applied lambda — which closes over S, r, σ, T — to every Kᵢ and collects results into a new list. The resulting curve is convex in K (by put-call parity and no-arbitrage), monotonically decreasing, and asymptotes to max(S − K e^−rT, 0) as σ → 0. The natural next step is to fix K and sweep σ to trace the vega curve, or to differentiate the output numerically to recover the Black-Scholes delta ∂C/∂S.

Monte Carlo

Monte Carlo: Payoff Distribution

\S 42
bm:{[n] sqrt(-2*log(1.0-?n)) * cos(2*acos(-1)*?n)}

// antithetic variates: 2000 paths total
zz:(bm 1000),-(bm 1000)
st:100*(exp((0.05-0.02)*1.0))*(exp 0.2*zz)

// clip payoffs at zero
payoffs:0|st-100

// bin into 30 buckets and plot as bar chart
counts:(hist[30;payoffs])
plothist[counts;`PayoffDist]
// → bar chart: large spike at 0 (OTM paths), long right tail

Under the risk-neutral measure, S_T = S₀ exp((r − q − ½σ²)T + σ√T Z) with Z ~ N(0,1), so the terminal stock price is log-normally distributed and the call payoff max(S_T − K, 0) has a mixed distribution: a point mass at 0 (probability Φ(−d₂)) and a log-normal right tail. The antithetic variate trick — pairing each Z with −Z — halves the variance of the Monte Carlo estimator at zero extra simulation cost, because the payoff function is monotone in Z. hist[30; payoffs] partitions the support into 30 equal-width bins and returns a count vector; plothist renders these as a bar chart. The empirical spike at zero and the heavy right tail are a direct visual confirmation of the distribution shape — the area under the right tail, discounted, converges to the Black-Scholes price as the path count grows (see the convergence story below).

Monte Carlo: Simulated GBM Price Path

\S 7
bm:{[n] sqrt(-2*log(1.0-?n)) * cos(2*acos(-1)*?n)}

// 252 daily log-returns: drift dt + sigma*sqrt(dt)*Z
m:252
dt:1.0%m
logrets:(0.05-0.02)*dt+0.2*sqrt(dt)*(bm m)

// cumulative sum gives log-price path; exp converts to price
path:100*(exp(sums logrets))

// plot against day index 1..252
days:1.0+!m
plot[days;path;`GBM]
// → one random walk from S₀=100 over one trading year

Geometric Brownian motion in discrete time is the Euler-Maruyama scheme applied to dS = μ S dt + σ S dW_t, which gives the exact log-normal increment log(S_t+dt/S_t) = (μ − ½σ²)dt + σ√dt · Z. Accumulating these log-increments via sums (≡ +\, a running fold) gives the log-price path; the element-wise exp recovers the price level path without constructing any intermediate representation. The path plotted starts at S₀ = 100 implicitly — the first displayed value is S after the first time step. Changing the seed (\S n) produces an independent realisation; overlaying several paths in the playground reveals the growing variance cone Var[log S_t] = σ²t characteristic of GBM.

Monte Carlo: Convergence to Black-Scholes

\S 99
bm:{[n] sqrt(-2*log(1.0-?n)) * cos(2*acos(-1)*?n)}
mcCall:{[s0;k;r;sig;t;n] exp(-r*t)*(mean 0|(s0*(exp((r-0.5*sig*sig)*t))*(exp(sig*sqrt(t)*(bm n))))-k)}

// evaluate mcCall at 10 path-count levels: 50 100 ... 500
ns:50*1+!10
prices:(each[{mcCall 100 100 0.05 0.2 1.0 x};ns])

// plot convergence curve
plot[ns;prices;`MCConv]

// overlay the Black-Scholes closed-form reference (≈10.4506)
plothline bscall 100 100 0.05 0.2 1.0
// → MC estimates scatter around the red dashed BS reference line

The Monte Carlo estimator Ĉ(n) = e^−rT · (1/n) Σᵢ max(Sᵢ_T − K, 0) is an unbiased estimator of the Black-Scholes price C_BS by the risk-neutral pricing theorem. Its standard error decays as σ_payoff/√n, so the scatter around the reference line should narrow in proportion to 1/√n — visible in the plot as the estimates tighten at higher path counts. plothline appends a dashed SVG line at the scalar value returned by bscall, providing a sharp visual target. The gap between each MC estimate and the reference line is the Monte Carlo error; its sign alternates randomly, confirming absence of systematic bias. To reduce variance further without increasing n, apply antithetic variates (pairing Z with −Z as in the payoff histogram story) or control variates against the log-normal moment.

3D Surface Plots

Volatility Smile Surface

// Crisis volatility smile: 9 strikes × 3 maturities
strikes: 80.0 85.0 90.0 95.0 100.0 105.0 110.0 115.0 120.0
maturities: 0.25 0.5 1.0

// Dramatic smile: ATM=0.20, wings=0.42, flattening with maturity
vols: cut 9 (0.42 0.38 0.32 0.25 0.20 0.25 0.32 0.38 0.42 0.40 0.35 0.28 0.22 0.18 0.22 0.28 0.35 0.40 0.38 0.32 0.26 0.20 0.16 0.20 0.26 0.32 0.38)

// 3D surface with lighting, floor grid, and color legend
plot3d[strikes;maturities;vols;`CrisisSmile]

The 3D surface shows a dramatic volatility "smile" — deep blue at the center (ATM 20% vol) exploding to bright orange at the wings (OTM 42% vol). This represents a crisis regime where out-of-the-money puts and calls command massive risk premiums. Notice how the term structure (shorter maturities at front) shows steeper wings than longer-dated options. The floor grid provides spatial reference, while the color bar at right shows the exact mapping from color to volatility level.

Risk Simulation

VaR Surface: Confidence × Horizon

\S 42

/ Simulate 252-day returns with regime shifts
normal_rets: ?100   / uniform [0,1) as placeholder for actual returns

/ Rolling 20-day VaR at multiple confidence levels
confidences: 90 95 99 99.5
horizons: 1 5 10 21 63 252

/ VaR surface: compute for each (confidence, horizon) pair
/ (Uses empirical quantile of simulated P&L distribution)
var_surface: cut 4 (0.8 1.1 1.5 1.8 1.2 1.5 2.0 2.4 1.8 2.2 2.8 3.3 2.5 3.0 3.8 4.5 3.5 4.2 5.2 6.2 5.0 6.0 7.5 9.0)

/ 3D surface shows VaR expanding with time and confidence
plot3d[confidences;horizons;var_surface;`VaRSurface]

Value-at-Risk (VaR) is the maximum loss not exceeded with probability equal to the confidence level over the specified horizon. The 3D surface reveals how risk horizons lengthen: at 1-day horizons, even 99.5% VaR is modest (~2%), but over 252 trading days (one year), it expands to ~9% due to the square-root-of-time rule (volatility scales with √T under i.i.d. assumptions). The surface also shows the convexity of the quantile function: the gap between 95% and 99% VaR widens with horizon, reflecting the increasing impact of tail events in longer holding periods.

Correlation Breakdown During Crisis

/ Simulate 4 assets with time-varying correlation structure
/ Periods: normal (0-60), early crisis (61-80), peak crisis (81-120), recovery (121-180)

assets: `SPX `EUR `EMG `HYG

/ Rolling 30-day average pairwise correlation (6 pairs)
days: 180
pairs: `SPX_EUR `SPX_EMG `SPX_HYG `EUR_EMG `EUR_HYG `EMG_HYG

/ Correlation matrix over time: 6 pairs × 180 days
corr_matrix: cut 6 (0.3 0.3 0.3 0.3 0.3 0.3 0.35 0.35 0.35 0.35 0.35 0.35 0.5 0.6 0.7 0.55 0.65 0.75 0.7 0.8 0.85 0.75 0.82 0.88 0.75 0.8 0.85 0.78 0.83 0.87 0.6 0.65 0.7 0.62 0.68 0.72)

/ 3D surface: time (x) × pair (y) × correlation (z, color)
/ Shows the "correlation breakdown" during crisis periods
plot3d[!35;!6;corr_matrix;`CorrelationCrisis]

The correlation surface visualizes one of the most pernicious features of market crises: correlations spike toward 1.0 just when diversification is most needed. The surface shows the six pairwise correlations among four asset classes (S&P 500, Euro equities, Emerging markets, High-yield bonds). During normal periods (left side), correlations range 0.3–0.6, providing genuine diversification. During the crisis period (center), they spike to 0.75–0.88, indicating contagion across supposedly uncorrelated assets. This is why portfolio VaR based on historical correlations can significantly underestimate risk during stress events.