Rama cont model implemented in Rust

Cargo.toml

[package]
name = "agent_based_trading_julia"
version = "0.1.0"
authors = ["cgmossa <cgmossa@gmail.com>"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
rand = {version = "0.7.3", features = ["small_rng"]}
rand_distr = "0.2.2"

[dev-dependencies]
criterion = "0.3"

[[bench]]
name = "cont_run"
harness = false

cont_run.rs

// resides in benches/cont_run.rs
//

use agent_based_trading_julia::agent_simulation::cont_run;
use criterion::{black_box, criterion_group, criterion_main, Criterion};

pub fn criterion_benchmark(c: &mut Criterion) {
    c.bench_function("cont_run default", |b| {
        b.iter(|| {
            black_box(cont_run(
                black_box(10_000),
                black_box(10_000),
                black_box(0.05),
                black_box(0.1),
            ))
        })
    });
}

criterion_group!(benches, criterion_benchmark);
criterion_main!(benches);

julia_impl.jl

using Random
using StatsBase

function cont_run(time=10000, n=10000, λ=0.05, q=0.1)
    r = zeros(time)
    θ = zeros(n)
    pchange = zeros(n)
    for t = 1:time
        ε = randn()
        if ε > 0
            r[t] =  sum(<(ε), θ) / (λ * n)
        else
            r[t] =  -sum(<(-ε), θ) / (λ * n)
        end
        θ .= ifelse.(rand!(pchange) .< q, abs(r[t]), θ)
    end
    return kurtosis(r)
end

lib.rs

// resides in src/lib.rs
use rand::prelude::SmallRng;
use rand::{thread_rng, Rng, SeedableRng};

pub fn cont_run(time: usize, n: usize, lambda: f64, q: f64) -> f64 {
    let mut theta = vec![0.; n];
    let n = n as f64;
    let mut eps_sampler = SmallRng::from_rng(thread_rng())
        .unwrap()
        .sample_iter(rand_distr::StandardNormal);
    let mut pchange_sampler = SmallRng::from_rng(thread_rng())
        .unwrap()
        .sample_iter(rand::distributions::Uniform::new_inclusive(0., 1.));
    let r = std::iter::repeat_with(move || {
        let eps: f64 = eps_sampler.next().unwrap();
        let r_t = if eps > 0. {
            theta.iter().filter(|&&x| eps > x).count() as f64 / (lambda * n)
        } else {
            -(theta.iter().filter(|&&x| -eps > x).count() as f64) / (lambda * n)
        };
        theta
            .iter_mut()
            .filter(|_| pchange_sampler.next().unwrap() < q)
            .for_each(|x| {
                *x = r_t.abs();
            });
        r_t
    });
    let r = r.take(time).collect::<Vec<_>>();
    kurtosis(r)
}

fn kurtosis(x: Vec<f64>) -> f64 {
    let n = x.len() as f64;
    let mean_x = x.iter().sum::<f64>() / n;
    let x = x.iter().copied().map(|x| x - mean_x);
    let r: f64 = n * x.clone().map(|x| x.clone().powi(4)).sum::<f64>()
        / (x.map(|x| x.powi(2)).sum::<f64>().powi(2));

    r * (1. - 1. / n).powi(2) - 3.
}

Right, I always set CC=clang, sorry I didn't communicate that.

That sounds better, but I don't think it has SIMD support?

Mersenne Twister typically comes with a very large state size, and still has some statistical deficiencies. There's a paper summarizing a few arguments against it https://arxiv.org/pdf/1910.06437.

Just for completion sake: Julia 1.5v didn't improve the benchmark results at all.