chelis
GitHub Get started
Menu

Examples

Programs from the Chelis repositories and shells, grouped by the job they do. Each listing comes from a named source file and shows its recorded result.

  • Pricing and risk

    Price a call and get its Greeks from grad

    Delta and vega are grad of the pricing function with respect to spot and volatility, so the Greeks come from the same code that sets the price.

    call_price = 10.450584
delta      = 0.6368296
vega       = 37.524048

  • Time series and dataframes

    Rolling mean, rolling std and EWM

    Coral's window functions return one value per input, so a five-point series stays tensor[5, f32] through every transform.

    mean3()    = [NaN, NaN, 2.0, 3.0, 4.0]
std3()     = [NaN, NaN, 1.0, 1.0, 1.0]

  • Proven properties

    Prove the Gordon growth model behaves

    cvc5 proves over the reals that the present value D / (r - g) is positive, rises with the dividend and falls as the required return rises.

    gordon_positive         passed  smt
gordon_increasing_in_d  passed  smt

Numerics and statistics

Root finding, integration and descriptive statistics from Nautilus.

Find the square root of 2 three ways

hello-chelis/src/nautilus/roots.chcomplete

Bisection, Brent and Newton each take the function itself as an argument and converge on the square root of 2.

module Hello.Nautilus.Roots
import Nautilus.Roots (bisection, newton, brent)
export (sqrt2_bisect, sqrt2_brent, sqrt2_newton)
def quad_minus_2(x: f32) -> f32 = sub(mul(x, x), 2.0f32)
def dquad_minus_2(x: f32) -> f32 = mul(2.0f32, x)
def sqrt2_bisect() -> f32 = bisection(quad_minus_2, 1.0f32, 2.0f32, 1e-8f32, 100i64)
def sqrt2_brent() -> f32 = brent(quad_minus_2, 1.0f32, 2.0f32, 1e-10f32, 100i64)
def sqrt2_newton() -> f32 = newton(quad_minus_2, dquad_minus_2, 1.5f32, 1e-10f32, 50i64)

Result

sqrt2_bisect() = 1.4142137
sqrt2_brent()  = 1.4142132
sqrt2_newton() = 1.4142135

Integrate sin over [0, pi]

hello-chelis/src/nautilus/integration.chcomplete

Trapezoid, Simpson and five-point Gauss-Legendre rules integrate the same function to its exact value of 2.

module Hello.Nautilus.Integration
import Nautilus.Integrate (trapezoidal, simpsons, gauss_legendre_5)
export (sin_integral_trap, sin_integral_simpson, sin_integral_gl5)
def sin_fn(x: f32) -> f32 = sin(x)
def pi_f() -> f32 = 3.1415927f32
def sin_integral_trap() -> f32 = trapezoidal(sin_fn, 0.0f32, pi_f(), 100i64)
def sin_integral_simpson() -> f32 = simpsons(sin_fn, 0.0f32, pi_f(), 100i64)
def sin_integral_gl5() -> f32 = gauss_legendre_5(sin_fn, 0.0f32, pi_f(), 5i64)

Result

sin_integral_trap()    = 1.9998353
sin_integral_simpson() = 2.0
sin_integral_gl5()     = 2.0000002

Describe a sample

hello-chelis/src/nautilus/stats.chcomplete

Every statistic borrows a vector of length n, and sample_correlation only accepts two vectors that share the same n.

module Hello.Nautilus.Stats
import Nautilus.Stats (mean_vec, variance_vec, std_vec, median_vec, quantile_vec, correlation_scalar)
export (sample_mean, sample_var, sample_std, sample_median, sample_quartile_low, sample_correlation)
def sample_mean[n](v: &tensor[n, f32]) -> f32 = mean_vec(v)
def sample_var[n](v: &tensor[n, f32]) -> f32 = variance_vec(v, 1i64)
def sample_std[n](v: &tensor[n, f32]) -> f32 = std_vec(v, 1i64)
def sample_median[n](v: &tensor[n, f32]) -> f32 = median_vec(v)
def sample_quartile_low[n](v: &tensor[n, f32]) -> f32 = quantile_vec(v, 0.25f32)
def sample_correlation[n](a: &tensor[n, f32], b: &tensor[n, f32]) -> f32 = correlation_scalar(a, b)

Result · xs = [2, 4, 4, 5, 7, 9], ys = [1, 3, 2, 5, 6, 8]

sample_mean(xs)          = 5.1666665
sample_std(xs)           = 2.4832773
sample_median(xs)        = 4.5
sample_quartile_low(xs)  = 4.0
sample_correlation(xs, ys) = 0.97134686

Time series and dataframes

Rolling windows, group-by and joins, with shapes carried in the types.

Rolling mean, rolling std and EWM

hello-chelis/src/coral/windowrolling.chcomplete

Coral's window functions return one value per input, so a five-point series stays tensor[5, f32] through every transform.

module Hello.Coral.WindowRolling
import Coral.Window (rolling_mean, rolling_std, ewm)
export (sample_series, mean3, std3, ewm_half)
def sample_series() -> tensor[5, f32] = to_tensor([1.0f32, 2.0f32, 3.0f32, 4.0f32, 5.0f32])
def mean3() -> tensor[5, f32] = rolling_mean(sample_series(), 3i64)
def std3() -> tensor[5, f32] = rolling_std(sample_series(), 3i64)
def ewm_half() -> tensor[5, f32] = ewm(sample_series(), 0.5f32)

Result

mean3()    = [NaN, NaN, 2.0, 3.0, 4.0]
std3()     = [NaN, NaN, 1.0, 1.0, 1.0]
ewm_half() = [1.0, 1.5, 2.25, 3.125, 4.0625]

Window reductions with the output shape in the type

chelis/examples/checked_window_geometry.chcomplete

Width-2, stride-2 windows over two rows of five give a tensor[2, 2, f32], and the annotation on each result is checked against that geometry.

-- Valid-padding windows preserve the leading axis and use the declared stride.
x: tensor[2, 5, f32] = [[1.0f32, 2.0f32, 3.0f32, 4.0f32, 5.0f32], [6.0f32, 7.0f32, 8.0f32, 9.0f32, 10.0f32]]
totals: tensor[2, 2, f32] = reduce_window_sum(&x, [2i64], [2i64])
means: tensor[2, 2, f32] = reduce_window_mean(&x, [2i64], [2i64])
maxima: tensor[2, 2, f32] = reduce_window_max(&x, [2i64], [2i64])
minima: tensor[2, 2, f32] = reduce_window_min(&x, [2i64], [2i64])

Result · each result has shape [2, 2]

totals = [3, 7, 13, 17]
means  = [1.5, 3.5, 6.5, 8.5]
maxima = [2, 4, 7, 9]
minima = [1, 3, 6, 8]

Group by city, then sum and average

hello-chelis/src/coral/groupbyagg.chcomplete

A four-row frame grouped on city becomes a Frame[2], one row per city.

module Hello.Coral.GroupByAgg
import Coral.Frame (Frame, Column, FloatCol, StringCol, from_pairs, nrows, ncols, get_float_col)
import Coral.GroupBy (group_by, agg_sum, agg_mean)
export (sales_frame, total_per_city, mean_per_city)
def sales_frame() -> Frame[4] = from_pairs([("city", StringCol(["a", "b", "a", "b"])), ("qty", FloatCol(to_tensor([1.0f32, 10.0f32, 2.0f32, 20.0f32])))])
def total_per_city() -> Frame[2] = agg_sum(group_by(sales_frame(), "city"), "qty")
def mean_per_city() -> Frame[2] = agg_mean(group_by(sales_frame(), "city"), "qty")

Result

city     = [a, b]
qty_sum  = [3.0, 30.0]
qty_mean = [1.5, 15.0]

Inner, left and outer joins

hello-chelis/src/coral/joins.chcomplete

Joining two three-row frames on customer gives 2, 3 and 4 rows, and each result type states its row count.

module Hello.Coral.Joins
import Coral.Frame (Frame, Column, FloatCol, StringCol, from_pairs, nrows)
import Coral.Join (inner_join, left_join, outer_join)
export (left_frame, right_frame, inner, left, outer)
def left_frame() -> Frame[3] = from_pairs([("customer", StringCol(["a", "b", "c"])), ("qty", FloatCol(to_tensor([1.0f32, 2.0f32, 3.0f32])))])
def right_frame() -> Frame[3] = from_pairs([("customer", StringCol(["b", "c", "d"])), ("score", FloatCol(to_tensor([10.0f32, 20.0f32, 30.0f32])))])
def inner() -> Frame[2] = inner_join(left_frame(), right_frame(), "customer")
def left() -> Frame[3] = left_join(left_frame(), right_frame(), "customer")
def outer() -> Frame[4] = outer_join(left_frame(), right_frame(), "customer")

Result

nrows(inner()) = 2
nrows(left())  = 3
nrows(outer()) = 4

Pricing and risk

Option prices, Greeks and value at risk.

Price a call and get its Greeks from grad

hello-chelis/src/capstone/blackscholes.chcomplete

Delta and vega are grad of the pricing function with respect to spot and volatility, so the Greeks come from the same code that sets the price.

module Hello.Capstone.BlackScholes
import Nautilus.Distributions (normal_cdf)
export (call_price, delta, vega)
def d1(s: f32, k: f32, r: f32, sigma: f32, t: f32) -> f32 = {
  num = add(log(div(s, k)), mul(add(r, mul(0.5f32, mul(sigma, sigma))), t))
  den = mul(sigma, sqrt(t))
  div(num, den)
}
def d2(s: f32, k: f32, r: f32, sigma: f32, t: f32) -> f32 = sub(d1(s, k, r, sigma, t), mul(sigma, sqrt(t)))
def call_price(s: f32, k: f32, r: f32, sigma: f32, t: f32) -> f32 = {
  d_1 = d1(s, k, r, sigma, t)
  d_2 = d2(s, k, r, sigma, t)
  discount = exp(neg(mul(r, t)))
  sub(mul(s, normal_cdf(d_1, 0.0f32, 1.0f32)), mul(mul(k, discount), normal_cdf(d_2, 0.0f32, 1.0f32)))
}
def delta(s: f32, k: f32, r: f32, sigma: f32, t: f32) -> f32 = grad(call_price, wrt=s)(s, k, r, sigma, t)
def vega(s: f32, k: f32, r: f32, sigma: f32, t: f32) -> f32 = grad(call_price, wrt=sigma)(s, k, r, sigma, t)

Result · S = K = 100, r = 0.05, sigma = 0.2, T = 1

call_price = 10.450584
delta      = 0.6368296
vega       = 37.524048

Parametric and historical value at risk

shoals/src/risk.chexcerpt

Both VaR measures take a loss vector of any length n and a confidence level, one from the normal quantile and one from the empirical quantile.

module Shoals.Risk
import Nautilus.Distributions (normal_inv_cdf)
import Nautilus.Stats (mean_vec, std_vec, quantile_vec)
def parametric_var[n](losses: tensor[n, f32], confidence: f32) -> f32 = {
  mu = mean_vec(&losses)
  sigma_loss = std_vec(&losses, 1i64)
  z = normal_inv_cdf(confidence, 0.0f32, 1.0f32)
  add(mu, mul(sigma_loss, z))
}
def historical_var[n](losses: tensor[n, f32], confidence: f32) -> f32 = quantile_vec(losses, confidence)

Result · 10 losses, confidence 0.95

parametric_var = 4.489891
historical_var = 4.13

Proven properties

Economic and statistical facts checked by an SMT solver over the reals.

Prove the Gordon growth model behaves

economoist/properties/growth.chexcerpt

cvc5 proves over the reals that the present value D / (r - g) is positive, rises with the dividend and falls as the required return rises.

module Economoist.Properties.Growth
import Economoist.Growth (gordon_pv, gordon_pv_strict)
@property gordon_positive forall(d: f32, r: f32, g: f32) where d > 0.0, r > g:
  (gordon_pv(d, r, g) > 0.0)
@property gordon_increasing_in_d forall(d1: f32, d2: f32, r: f32, g: f32) where r > g, d2 > d1:
  (gordon_pv(d2, r, g) > gordon_pv(d1, r, g))
@property gordon_decreasing_in_r forall(d: f32, r1: f32, r2: f32, g: f32) where d > 0.0, r1 > g, r2 > r1:
  (gordon_pv(d, r2, g) < gordon_pv(d, r1, g))

chelis prove

gordon_positive         passed  smt
gordon_increasing_in_d  passed  smt
gordon_decreasing_in_r  passed  smt

Quantiles are monotone

chelis/examples/nautilus_quantile_contract/src/proofs.chcomplete

The prover shows that for any p <= q the p-quantile of the sample is at most the q-quantile, using the monotonicity contract Nautilus declares for quantile_vec.

module QuantileContract.Proofs
import Nautilus.Stats (quantile_vec)
def observations() -> tensor[5, f32] = (to_tensor([8.0, 1.0, 5.0, 3.0, 2.0], f32) : tensor[5, f32])
@property example_quantiles_are_monotone forall(p: f32, q: f32) where 0.0 <= p, p <= q, q <= 1.0:
  (quantile_vec(observations(), p) <= quantile_vec(observations(), q))
  with contract = "std.quantile.monotonicity"

chelis prove

example_quantiles_are_monotone  passed  smt

Generic code

One function for f32 and f64, and a moving sum over any numeric type.

One contract for f32 and f64

chelis/examples/declared_operation_bounds.chexcerpt

average3 declares that it works for any float type, and runs on f32 and f64 data; moving_sum states the same kind of contract for a two-point window over any numeric type.

-- A generic definition states the dtype family its operations require.
sig average3[p: Float]: tensor[3, p] -> tensor[p]
def average3(values) = mean(values, 0i32)
-- A wrapper preserves that contract, even without naming the primitive.
def summarize[p: Float](values: tensor[3, p]) -> tensor[p] = average3(values)
float_result = to_tensor([1.0f32, 2.0f32, 6.0f32]) |> summarize
double_result = to_tensor([1.0f64, 2.0f64, 6.0f64]) |> summarize
-- Window reductions carry the same operation contract through their
-- specialized shape checker.
def moving_sum[p: Numeric](values: tensor[3, p]) -> tensor[2, p] = reduce_window_sum(values, [2i64], [1i64])
window_result = to_tensor([1.0f32, 2.0f32, 4.0f32]) |> moving_sum

Result

float_result  = 3.0
double_result = 3.0
window_result = [3.0, 6.0]

The libraries these programs import