// docs: stats#realized-variance-from-intraday-returns
//! Backs the realized-variance example on the statistics catalog page.

use ndarray::Array1;
use stochastic_rs::stats::realized::variance::realized_variance;

#[test]
fn realized_variance_from_log_returns() {
  // Stand-in for "1-min log-returns over a trading day": a small fixed
  // series keeps this example free of a live data feed.
  let log_returns: Array1<f64> = Array1::from(vec![
    0.001, -0.0007, 0.0012, -0.0003, 0.0009, -0.0011, 0.0004, -0.0002, 0.0006, -0.0008,
  ]);

  let rv = realized_variance(log_returns.view());
  assert!(rv > 0.0);

  let annualised_vol = (rv * 252.0).sqrt();
  assert!(annualised_vol.is_finite());
}

import stochastic_rs as srs
import numpy as np

log_returns = ...   # 1-min log-returns over a trading day
# annualisation is a constructor argument, not a post-hoc multiplier — it
# scales straight into `.volatility`. There is no standard-error field:
# RealizedMoments reports variance/volatility/skewness/kurtosis/quarticity,
# not a point estimate with a `se` — see the ADF section below for the
# same "no se/p_value field" pattern on the stationarity tests.
res = srs.RealizedMoments(log_returns, annualisation=252.0)
print("RV     =", res.variance)
print("annVol =", res.volatility)

// docs: stats#ewma-riskmetrics-variance
//! Backs the EWMA / RiskMetrics example on the stats page.

use ndarray::Array1;
use stochastic_rs::stats::realized::ewma::RISKMETRICS_DAILY_LAMBDA;
use stochastic_rs::stats::realized::ewma::ewma_variance;
use stochastic_rs::stats::realized::ewma::riskmetrics_variance;

#[test]
fn ewma_variance_tracks_a_return_series() {
  let returns = Array1::from(vec![0.012_f64, -0.007, 0.021, 0.0, -0.015, 0.004]);

  // RiskMetrics daily decay (λ = 0.94) is the default convention…
  let rm = riskmetrics_variance(returns.view());
  assert_eq!(rm.lambda, RISKMETRICS_DAILY_LAMBDA);
  assert_eq!(rm.variance.len(), returns.len());

  // …and any λ ∈ (0, 1) is accepted explicitly.
  let fast = ewma_variance(returns.view(), 0.80);
  assert!(fast.forecast > 0.0);
  // A faster decay weights the newest shock harder, so the two forecasts
  // genuinely differ.
  assert!((fast.forecast - rm.forecast).abs() > 1e-9);
}

import stochastic_rs as srs

# lambda_ defaults to the RiskMetrics daily 0.94
res = srs.EwmaVariance(returns, lambda_=0.94)
print("sigma2 series:", res.variance())
print("1-step forecast:", res.forecast)

// docs: stats#lee-mykland-jump-test
//! Backs the Lee–Mykland example on the stats page.

use ndarray::Array1;
use stochastic_rs::distributions::normal::SimdNormal;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::realized::lee_mykland::lee_mykland_test;
use stochastic_rs::stats::realized::lee_mykland::lee_mykland_window;

#[test]
fn lee_mykland_locates_a_planted_jump() {
  // Stand-in for a day of 5-minute log-returns with one news jump in it.
  let dist = SimdNormal::<f64>::new(0.0, 0.001, &Deterministic::new(42));
  let mut returns = Array1::<f64>::zeros(2_000);
  dist.fill_slice(returns.as_slice_mut().unwrap());
  returns[1_200] = 0.02;

  // The paper's window rule for 78 five-minute returns a day is K = 141.
  let window = lee_mykland_window(78);
  assert_eq!(window, 141);

  let test = lee_mykland_test(returns.view(), window, 0.01);
  assert!(test.jump_indices.contains(&1_200));
  // Every flagged return sits beyond the Gumbel threshold on |L(i)|.
  for &i in &test.jump_indices {
    assert!(test.statistics[i].abs() > test.threshold);
  }
}

import stochastic_rs as srs

window = srs.LeeMyklandJumpTest.recommended_window(78)  # 5-minute data → 141
test = srs.LeeMyklandJumpTest(returns, window, alpha=0.01)
print("jumps at:", test.jump_indices)
print("threshold on |L(i)|:", test.threshold)

// docs: stats#johansen-cointegration-and-vecm
//! Backs the Johansen / VECM example on the stats page.

use ndarray::Array2;
use stochastic_rs::distributions::normal::SimdNormal;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::econometrics::johansen_test;
use stochastic_rs::stats::econometrics::vecm_fit;

#[test]
fn johansen_rank_then_vecm() {
  // Two prices sharing one random walk: y2 ≈ 0.7·y1, so the pair has
  // exactly one cointegrating relation.
  let steps = SimdNormal::<f64>::new(0.0, 1.0, &Deterministic::new(7));
  let noise = SimdNormal::<f64>::new(0.0, 0.1, &Deterministic::new(11));
  let mut dw = vec![0.0_f64; 500];
  let mut eps = vec![0.0_f64; 500];
  steps.fill_slice(&mut dw);
  noise.fill_slice(&mut eps);
  let mut y = Array2::<f64>::zeros((500, 2));
  let mut w = 0.0;
  for t in 0..500 {
    w += dw[t];
    y[[t, 0]] = w;
    y[[t, 1]] = 0.7 * w + eps[t];
  }

  // Rank tests at VAR order 2: both sequential procedures stop at r = 1.
  let test = johansen_test(y.view(), 2);
  assert_eq!(test.rank_trace, 1);
  assert_eq!(test.rank_max_eig, 1);
  assert!(test.max_eig_statistics[0] > test.max_eig_critical_5pct[0]);

  // The VECM at that rank recovers the relation y2 − 0.7·y1 up to scale.
  let fit = vecm_fit(y.view(), 2, test.rank_trace);
  let ratio = fit.beta[[1, 0]] / fit.beta[[0, 0]];
  assert!((ratio + 1.0 / 0.7).abs() < 0.1, "beta ratio {ratio}");
  assert_eq!(fit.gamma.len(), 1);
  assert_eq!(fit.residuals.dim(), (498, 2));
}

import stochastic_rs as srs

test = srs.Johansen(prices, lags=2)          # prices: (t, k) array
print("rank (trace):", test.rank_trace, "rank (max-eig):", test.rank_max_eig)
print("max-eig stats:", test.max_eig_statistics())

fit = srs.Vecm(prices, lags=2, rank=test.rank_trace)
print("pi:", fit.pi())
print("gamma_1:", fit.gamma()[0])

// docs: stats#garch-family-qmle
//! Backs the GARCH QMLE example on the stats page.

use ndarray::Array1;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::garch::GarchSpec;
use stochastic_rs::stats::garch::MeanSpec;
use stochastic_rs::stats::garch::garch_fit;
use stochastic_rs::stochastic::autoregressive::garch::Garch;
use stochastic_rs::traits::ProcessExt;

#[test]
fn garch_qmle_recovers_the_simulated_parameters() {
  // Two thousand daily returns from GARCH(1,1) with ω = 0.05, α = 0.10,
  // β = 0.85 (percent units, unconditional variance 1.0).
  let process = Garch::<f64, _>::new(
    0.05,
    Array1::from(vec![0.10]),
    Array1::from(vec![0.85]),
    2_000,
    Deterministic::new(42),
  );
  let returns = process.sample();

  // The simulator has no drift, so a zero-mean fit is the matching model;
  // `GarchSpec::garch(1, 1)` alone would estimate a constant mean as well.
  let fit = garch_fit(
    returns.view(),
    GarchSpec::garch(1, 1).with_mean(MeanSpec::Zero),
  );
  assert!(fit.converged);
  assert!((fit.alpha[0] - 0.10).abs() < 3.0 * fit.robust_std_errors[1]);
  assert!((fit.beta[0] - 0.85).abs() < 3.0 * fit.robust_std_errors[2]);
  assert!(fit.persistence < 1.0);
  assert_eq!(fit.conditional_variance.len(), returns.len());
}

import stochastic_rs as srs

fit = srs.GarchFit(returns, kind="gjr", p=1, q=1, mean="constant")
for name, value, se in zip(fit.param_names(), fit.params(), fit.robust_std_errors()):
    print(f"{name:9s} {value: .5f} ({se:.5f})")
print("persistence:", fit.persistence, "log-likelihood:", fit.log_likelihood)

// docs: stats#extreme-value-theory-hill-peaks-over-threshold-and-block-maxima
//! Backs the extreme-value example on the stats page.

use ndarray::Array1;
use stochastic_rs::distributions::pareto::SimdPareto;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::evt::block_maxima;
use stochastic_rs::stats::evt::gev_fit;
use stochastic_rs::stats::evt::hill_estimator;
use stochastic_rs::stats::evt::pot_fit;

#[test]
fn evt_tail_estimates_agree_on_a_pareto_tail() {
  // Losses with an exact Pareto(α = 3) tail, so the true tail index is
  // ξ = 1/3.
  let dist = SimdPareto::<f64>::new(1.0, 3.0, &Deterministic::new(7));
  let mut losses = vec![0.0; 20_000];
  dist.fill_slice(&mut losses);
  let losses = Array1::from(losses);

  let hill = hill_estimator(losses.view(), 500);
  assert!((hill.xi - 1.0 / 3.0).abs() < 3.0 * hill.std_error);

  // Peaks over a threshold of 3: the GPD shape estimates the same ξ, and
  // the tail quantiles follow.
  let pot = pot_fit(losses.view(), 3.0);
  assert!(pot.gpd.converged);
  assert!((pot.gpd.xi - 1.0 / 3.0).abs() < 3.0 * pot.gpd.std_errors[1]);
  let var99 = pot.quantile(0.99);
  assert!(pot.expected_shortfall(0.99) > var99 && var99 > pot.threshold);

  // Block maxima of 100 draws: GEV with ξ ≈ 1/3 again.
  let maxima = block_maxima(losses.view(), 100);
  let gev = gev_fit(maxima.view());
  assert!(gev.converged);
  assert!((gev.xi - 1.0 / 3.0).abs() < 3.0 * gev.std_errors[2]);
  assert!(gev.return_level(50.0) > gev.mu);
}

import numpy as np
import stochastic_rs as srs

losses = -returns                                   # right tail = losses
hill = srs.HillEstimator(losses, k=100)
print("xi:", hill.xi, "+/-", hill.std_error)

pot = srs.PotFit(losses, threshold=np.quantile(losses, 0.95))
print("VaR99:", pot.quantile(0.99), "ES99:", pot.expected_shortfall(0.99))

gev = srs.GevFit(srs.block_maxima(losses, 21))    # monthly maxima
print("100-block return level:", gev.return_level(100.0))

// docs: stats#fitting-the-skewed-return-distributions
//! Backs the distribution-fitting example on the stats page.

use ndarray::Array1;
use stochastic_rs::distributions::skew_t::SimdSkewT;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::distfit::johnson_su_fit;
use stochastic_rs::stats::distfit::skew_t_fit;

#[test]
fn skew_t_fit_recovers_the_simulated_shape() {
  // Returns from Hansen's skew-t with η = 6, λ = −0.3, scaled to 1.5%
  // daily volatility around a 0.05% drift.
  let dist = SimdSkewT::<f64>::new(6.0, -0.3, &Deterministic::new(11));
  let mut z = vec![0.0; 4_000];
  dist.fill_slice(&mut z);
  let returns = Array1::from_iter(z.iter().map(|v| 0.0005 + 0.015 * v));

  let fit = skew_t_fit(returns.view());
  assert!(fit.converged);
  assert!((fit.lambda + 0.3).abs() < 3.0 * fit.std_errors[3]);
  assert!((fit.sigma - 0.015).abs() < 3.0 * fit.std_errors[1]);

  // A Johnson SU fit of the same sample ranks below it on AIC.
  let jsu = johnson_su_fit(returns.view());
  assert!(jsu.converged && jsu.aic > fit.aic);
}

import stochastic_rs as srs

fit = srs.SkewTFit(returns)
print("eta:", fit.eta, "lambda:", fit.lambda_, "sigma:", fit.sigma)
print("SEs:", fit.std_errors(), "AIC:", fit.aic)

vg = srs.VarianceGammaFit(returns)
print("sigma, nu, theta, mu:", vg.sigma, vg.nu, vg.theta, vg.mu)

// docs: stats#adf-stationarity-test
//! Backs the ADF example on the statistics catalog page. `adf_test` lives
//! under the `stationarity` module (OLS regression via
//! `ndarray-linalg`), so the whole file is gated on that feature.

use ndarray::ArrayView1;
use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::stationarity::adf::AdfConfig;
use stochastic_rs::stats::stationarity::adf::adf_test;
use stochastic_rs::stochastic::diffusion::ou::Ou;
use stochastic_rs::traits::ProcessExt;

#[test]
fn adf_rejects_unit_root_on_a_mean_reverting_series() {
  // An OU path is mean-reverting by construction, so the null (unit
  // root) should be rejected — a stand-in for "prices or log-returns"
  // that avoids a live data feed.
  let series = Ou::<f64, _>::new(
    2.0,
    0.0,
    0.5,
    500,
    Some(0.0),
    Some(10.0),
    Deterministic::new(7),
  )
  .sample();

  let result = adf_test(
    ArrayView1::from(series.as_slice().unwrap()),
    AdfConfig::default(),
  );
  assert!(
    result.reject_unit_root,
    "expected the unit-root null to be rejected"
  );
}

import stochastic_rs as srs

# ADFTest is a class, not a free function; the fitted
# object exposes .statistic/.reject_unit_root/.critical_values directly —
# no .p_value, matching the "no separate p_value field" note above.
res = srs.ADFTest(series, max_lags=5)
print(f"t-stat={res.statistic:.3f}, reject_unit_root={res.reject_unit_root}")

// docs: stats#hurst-exponent-estimation
//! Backs the Hurst-estimator example on the statistics catalog page.
//!
//! Uses the rescaled-range estimator (`RescaledRange`) rather than the
//! Fukasawa/Whittle one (`stats::hurst::whittle`): Fukasawa estimates the
//! roughness of a *latent volatility* process from a realized-variance
//! series, not the Hurst exponent of a path sampled directly from `Fgn` —
//! feeding path samples straight into it silently recovers nonsense. See
//! `stochastic-rs-stats/src/hurst/whittle.rs`'s own `simulate_log_rv` test
//! helper for the realized-variance construction Fukasawa actually expects.

use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stats::hurst::HurstEstimator;
use stochastic_rs::stats::hurst::rs::RescaledRange;
use stochastic_rs::stochastic::noise::fgn::Fgn;
use stochastic_rs::traits::ProcessExt;

#[test]
fn rescaled_range_recovers_true_h_from_fgn_increments() {
  let true_h = 0.3;
  let path = Fgn::<f64, _>::new(true_h, 4096, Some(1.0), Deterministic::new(42)).sample();

  // `path` is already the stationary fGn increment series, so disable the
  // estimator's default first-differencing (meant for an fBM-like walk).
  let estimator = RescaledRange {
    take_differences: false,
    ..RescaledRange::default()
  };
  let res = estimator.estimate(path.view()).unwrap();
  assert!(
    (res.hurst - true_h).abs() < 0.1,
    "H = {:.3} (true = {:.3})",
    res.hurst,
    true_h
  );
}

import stochastic_rs as srs

x = ...   # a price / level series, e.g. np.cumsum(fgn.sample())
res = srs.RescaledRange().estimate(x)
print(f"H = {res.hurst:.3f}")