API.md

Classes

HurstifyError ⇐ Error

hurstify-specific error with a stable code field.

Hurstify

Randomized Kolmogorov-Smirnov Analysis of Volatility Roughness estimator.

Wraps the configuration (scales, sample size, sampler, optimizer, KS objective, h bounds) and exposes the same public API as v1.x — but every pluggable concern now lives behind a strategy.

Lifecycle:

1. The constructor stores the configuration and resolves the sampler / KS-objective / optimizer strategies once. The estimator is therefore stateful across calls — each call to estimate draws a fresh independent PRNG state via prng.js.
2. estimate averages iterations independent estimateSingle results for variance reduction.
3. rolling and rollingMultiScale are convenience wrappers around a sliding window of these estimates.
4. estimateBatch performs non-overlapping window estimation for parallel processing pipelines.

Optimizer

Strategy base class. Subclasses implement minimize polymorphically.

Forecaster

Abstract base class for H-series forecasters.

ArfimaForecaster

ARFIMA(p, d, q) forecaster.

HoltWintersForecaster

Holt-Winters (level + trend) forecaster.

LstmForecaster

Stateless LSTM-like recurrent cell forecaster.

AttentionForecaster

Stateless single-head self-attention block forecaster.

HypothesisTest

KsSignificanceTest

KS-distance significance test on the minimized statistic returned by Hurstify.estimateSingle.

Under the null of self-similarity at the estimated H the minimized KS distance should be near the asymptotic critical value; rejecting the null suggests the estimator should be treated with caution.

ConstancyTest

Likelihood-ratio constancy test for a series of H estimates under a 1D Kalman-filter state-space model.

CusumBreakTest

One-sided CUSUM structural-break detector on standardized residuals.

BootstrapConfidenceInterval

Percentile bootstrap CI for an arbitrary estimator function.

Kernel

Abstract base class for fractional-integration kernels.

RiemannLiouvilleKernel

Riemann–Liouville kernel K(t) = sqrt(2 H) * t^{H - 0.5} for H in (0, 1) and t > 0.

This is the kernel of choice for rBergomi and the exact-OU / mPRE simulators.

TimeVaryingKernel

Time-varying kernel whose local exponent H(t) is sampled at every step. The kernel evaluator picks the exponent based on t:

K(t) = sqrt(2 * H(t)) * t^{H(t) - 0.5}

The caller supplies a hPath array whose i-th entry is the local Hurst exponent at time i * dt. When H(t) is constant the kernel collapses to the Riemann–Liouville form.

KsObjective

Abstract base class for KS-distance objectives.

PairwiseKsObjective

Two-scale pairwise KS objective.

Rescales the two sorted samples by scales[i]^{-H} and returns the KS distance between them.

MultiScaleKsObjective

Multi-scale unweighted KS objective.

Returns the arithmetic mean of the pairwise KS distances over every unordered pair (i, j) with i < j. Equivalent to the paper's recommended extension to K > 2 scales.

WeightedMultiScaleKsObjective

Multi-scale weighted KS objective.

Weights the pair (i, j) by weights[i] * weights[j] and normalizes by the sum of those weights so the result stays in [0, 1] regardless of the absolute weight magnitudes.

StochasticModel

Abstract base class for stochastic-process simulators.

RoughBergomiModel

Rough Bergomi model.

dV_t / V_t = eta * dW^perp_t
I_t = int_0^t sqrt(2H) (t - s)^{H - 0.5} dW^perp_s
V_t = xi * exp(eta I_t - (eta^2 / 2) t^{2H})

RoughFsvModel

Rough Fractional Stochastic Volatility model.

dV_t = theta (mu - V_t) dt + nu V_t^alpha dW^V_t + roughComp

FractionalOuModel

Abstract base class for the Fractional Ornstein-Uhlenbeck model.

dX_t = theta (mu - X_t) dt + sigma dB^H_t

Concrete subclasses pick the discretization scheme:

• EulerMaruyamaFractionalOuModel — default, cheap
• ExactFractionalOuModel — Riemann-Liouville integral, accurate

EulerMaruyamaFractionalOuModel

Euler-Maruyama discretization of the fOU model.

Cheap O(n) integration; first-order accurate.

ExactFractionalOuModel

Exact Riemann-Liouville discretization of the fOU model.

O(n^2) per path; higher-order accurate.

MultifractionalPreModel

Abstract base class for the Multifractional Process with Random Exponent.

X_t = B_{H(t)}(t)

where H(t) itself is a stochastic Ornstein-Uhlenbeck process bounded between hMin and hMax. Two concrete subclasses pick the discretization scheme:

• LocalHolderMultifractionalPreModel — default, cheap
• ExactMultifractionalPreModel — time-varying kernel

LocalHolderMultifractionalPreModel

Local-Holder approximation of the mPRE model.

Cheap O(n) integration via cumulative sqrt(dt^{2 * H_avg}) scaling.

ExactMultifractionalPreModel

Exact time-varying-kernel discretization of the mPRE model.

O(n^2) per path; uses a time-varying Riemann-Liouville kernel.

Sampler

Abstract base class for sampling strategies.

A Sampler is a strategy: callers obtain a fresh instance and invoke draw(inc, n) once per variance-reduction iteration. The estimator never holds sampler state across calls so strategies can be safely shared across estimator instances.

ReservoirSampler

Floyd's Algorithm R reservoir sampler wrapped as a Sampler strategy.

Complexity: O(inc.length) time, O(n) extra memory. Reproducible when prng.setRandomSeed() has been called.

BlockPermutationSampler

Block random permutation followed by a reservoir draw.

This is the paper-faithful RK-SAVR pipeline: the increments are first sliced into blocks of length blockSize (optionally with a random phase offset) and the blocks are shuffled, then the desired number of increments is drawn without replacement from the permuted array.

IdentitySampler

Identity sampler — returns the input unchanged.

Useful when the caller has already prepared an array of exactly n elements (e.g. in deterministic unit tests).

Constants

NORMAL_QUANTILE_COEFFS

Coefficients for the Beasley-Springer-Malkin rational approximation of the inverse standard normal CDF. Used piecewise for p in [pLow, 1 - pLow] (central region) and tail rational functions for the extremes.

The standard deviation may be c/d constants at the tails is adapted from Peter Acklam's algorithm.

modelRegistry : Registry.<StochasticModel>

Strategy registry for stochastic models. The default fOU key resolves to the Euler-Maruyama discretization; consumers who want the exact Riemann-Liouville variant look up fOU-exact. Same convention for mPRE / mPRE-exact.

forecasterRegistry : Registry.<Forecaster>

Strategy registry for forecasters.

optimizerRegistry : Registry.<Optimizer>

Global optimizer registry.

Functions

parseCsv(csv, opts) ⇒ Array.<Object>

Parses a CSV string into an array of plain objects.

Expected input shape:

• The first non-empty line is the header row.
• Each subsequent line is a record with the same column count as the header.
• Fields can be optionally wrapped in double quotes; quotes may embed commas but not other escapes.

Type coercion:

• opts.dateField (default "date") is parsed via new Date(...).
• Any field listed in opts.numericFields is parsed via parseFloat.
• All other fields are kept as trimmed strings.

Error handling:

• Empty input returns [].
• Mismatched column counts throw with a descriptive message.
• Non-numeric values in declared numeric columns throw.

splitCSVLine(line) ⇒ Array.<string>

Splits a single CSV line respecting double-quoted regions.

States:

• Outside quotes: a comma terminates the current field.
• Inside quotes: a quote toggles back to "outside", all other chars are kept verbatim.

extractSeries(rows, field, opts) ⇒ Array.<{date: Date, value: number}>

Extracts a {date, value} series from a parsed CSV array.

Rows that are missing field are skipped; the resulting series is optionally sorted by dateField when the caller asks. Sorting uses the standard JS Date arithmetic, so the dates must be real Date instances.

parseJson(json) ⇒ Array.<Object>

Parses a JSON string that must encode an array of objects.

The function deliberately refuses non-array JSON to keep the loader simple. Empty or whitespace-only input returns [].

validateNoGaps(series, maxGapMs) ⇒ Object

Validates that a time series does not contain temporal gaps larger than maxGapMs.

Returns the maximum observed gap, the full list of pairwise gap lengths, and a valid flag for the threshold check. Series with fewer than two points are deemed valid by definition.

downsampleSeries(series, intervalMs) ⇒ Array.<{date: Date, value: number}>

Downsamples a time series by averaging values that fall into fixed intervalMs-wide buckets.

The bucket index is computed as floor(date.getTime() / intervalMs), so all buckets share the same left edge (0, intervalMs, 2 * intervalMs, ...). The output is sorted by date and every returned point carries the bucket start (not the average timestamp) as its date value.

preaverageReturns(prices, [windowSize]) ⇒ Array.<number>

Preaveraging of log-returns.

Implementation of the Jacod et al. (2009) preaveraging estimator (simplified single-bar variant):

1. Compute log-returns r_t = log(P_t / P_{t-1}).
2. For each i, average the windowSize consecutive returns ending at i (g_avg[i] = mean(r_{i - windowSize + 1}, ..., r_i)).
3. The "preaveraged return" is the first-difference sequence g_avg[i] - g_avg[i - 1]. This cancellation attenuates microstructure noise by 1/sqrt(windowSize) while preserving the drift and diffusion up to O(1 / windowSize).

Note: the result has length prices.length - windowSize - 1; for very short series the function throws rather than returning a few noisy points.

computeRealizedKernel(returns, [kernelType], [bandwidth]) ⇒ number

Realized-kernel variance estimator with pluggable kernels.

Given n returns, the estimator forms the autocorrelation sequence

gamma_k = sum_{i=k+1}^{n} r_i * r_{i - k},  k = 0..h

and combines them through a weighted sum

RV_K = gamma_0 + 2 * sum_{k=1..h} w_k * gamma_k

with weights w_k provided by the chosen kernel. The default bandwidth is floor(n^0.6), a rule-of-thumb that matches the optimal scaling under i.i.d. microstructure noise.

Kernels shipped:

• bartlett: w_k = 1 - k / h (default).
• parzen: the standard piecewise-cubic Parzen kernel.
• tukey-hanning: 0.5 (1 + cos(pi k / h)).

Any unknown kernel name falls back to Bartlett.

kernelWeight(type, k, h) ⇒ number

Kernel weight function used by realizedKernel.

debiasLogVolatility(rawHEstimates, sigmaObs, sigmaLatent) ⇒ Array.<number>

Heuristic de-biasing of log-volatility H estimates.

Microstructure noise inflates the variance of the log-volatility proxy relative to the latent signal, which in turn attenuates the observed roughness. This routine adds a small correction

h_debias = h + 0.01 * log(sigmaObs / sigmaLatent)

and clamps the result to [0.01, 0.99]. It is intentionally conservative — the user is expected to validate the calibration against a trust sample before relying on it for production.

computeRealizedVariance(prices, [interval]) ⇒ Array.<number>

Computes per-bucket realized variance from a price series.

The realized variance is the sum of squared log-returns within each non-overlapping bucket of interval observations:

RV_k = sum_{i in bucket k} (log P_i - log P_{i-1})^2

With interval = 1 the function emits one RV per log-return directly, which is the canonical "5-minute RV" form when prices are already sampled at 5-minute intervals.

computeRealizedVarianceParkinson(bars) ⇒ Array.<number>

Parkinson (1980) high-low RV estimator from OHLC bars.

For each bar the within-period variance is approximated by

sigma^2 ~= (log(H/L))^2 / (4 * ln 2)

which is 1/(4 ln 2) ~ 0.36 of the log-range-squared. Parkinson is strictly less efficient than tick-based RV but only requires four numbers per bar.

aggregateDailyRealizedVariance(intradayRVs) ⇒ number

Aggregates intraday (5-minute) realized variances into a single daily value via plain summation.

This is the standard "sum of squared returns" daily RV used in financial econometrics. It assumes the input is already free of overnight gaps.

applyLogTransform(rv) ⇒ Array.<number>

Maps realized variance to the log-volatility series consumed by hurstify.

The transformation is

X_t = 0.5 * log(RV_t)

i.e. log(sqrt(RV)). This converts multiplicative variance dynamics into a roughly additive (and therefore more stationary) signal, on top of which the self-similarity property exploited by the RK-SAVR algorithm is expressed.

centerSeries(series) ⇒ Array.<number>

Subtracts the arithmetic mean from every element.

Useful as a final step in the preprocessing pipeline when the user wants the series to mean-zero (which can stabilize variance-reducing permutations inside Hurstify).

standardizeSeries(series) ⇒ Array.<number>

Standardizes a time series to zero mean and unit variance.

Divides each centered value by the population standard deviation. A constant series has zero variance and triggers an explicit error rather than silently producing NaNs.

applyPreprocessingPipeline(prices, opts) ⇒ Array.<number>

Bundled preprocessing pipeline: prices -> RV -> log-vol -> (optional) centering.

Equivalent to running computeRV + logTransform + (optionally) centerSeries, but more compact for callers who want the canonical transformation.

splitTrainTest(series, [trainRatio]) ⇒ Object

Splits a series into contiguous training and test arrays.

The split point is floor(series.length * trainRatio) so the training set is the leftmost prefix of the series; this preserves temporal ordering, which is what hurstify forecasters and validation scripts typically need.

createSlidingWindows(series, windowSize, [step]) ⇒ Array.<Array.<number>>

Builds overlapping windows from a single time series.

The i-th window is series.slice(i, i + windowSize) for i = 0, step, 2*step, ... until no full window fits. Used by offline batch evaluation pipelines that want to score the estimator on every available segment of the series.

generateVixLogVolatility(nDays, h, opts) ⇒ Array.<number>

Synthetic VIX-style daily log-volatility.

Generates an fBM with the requested h and maps it to a log-volatility level around 2.0 (i.e. sqrt(RV) ~ 20%) by adding a small drift term and Gaussian observation noise:

X_t = 2.0 + drift * (fbm[t] / sqrt(n)) + 0.5 * fbm[t] + noise

Default tuning matches the empirical VIX roughness (h ~ 0.1) and annualized log-vol mean.

generateSpxLogVolatility(nDays, h, opts) ⇒ Array.<number>

Synthetic S&P 500 realized-volatility style daily log-volatility.

Same construction as generateVIXLogVol but with a smoother default Hurst (h = 0.14), a smaller drift, and a less volatile observation-noise level. Empirically these choices match the rough regime typically reported for SPX RV.

generateIntradayPrices([nIntraday], [nDays], h, opts) ⇒ Array.<Array.<number>>

Generates synthetic intraday 5-minute prices useful for testing realized-variance pipelines.

For every (re-)sampled day the generator draws an fBM with the requested h, exponentiates it into a volatility factor, and steps a log-return process

S_{i+1} = S_i * exp(drift + vol_i * z_i * sqrt(dt))

with drift set to the per-5-minute-bar annualized drift. The result is a nDays x nIntraday array of prices suitable for feeding into computeRV.

seriesToCsv(series, [dateHeader], [valueHeader]) ⇒ string

Serializes a {date, value} series as a CSV string.

Dates that are Date instances are formatted as their ISO yyyy-mm-dd prefix; everything else is stringified verbatim. Empty series produces a header-only CSV.

buildScaleProfile(sortedSamples, scales, H) ⇒ Array.<number>

Builds a flat "profile" of all pairwise KS distances at a fixed H.

Given K sorted samples, the profile has K * (K - 1) / 2 entries corresponding to every unordered scale pair. Useful for diagnostics.

getAsymptoticVariance(scaleA1, scaleA2, n, m) ⇒ number

Asymptotic variance of the Hurstify estimator.

Implements

Var(H_hat) = (2 * pi * e) / (ln(a2/a1))^2 * (1/sqrt(n) + 1/sqrt(m))^2.

When a1 == a2 (log ratio zero) the variance is degenerate and the function returns Infinity rather than dividing by zero; callers that intend to compute a SE/CI should reject equal scales up-front.

getStandardError(scaleA1, scaleA2, n, m) ⇒ number

Asymptotic standard error: square root of the asymptotic variance.

Thin convenience wrapper. The standard error has units of "Hurst" and can be read against the hMin/hMax bounds the estimator was configured with.

getConfidenceInterval(hEstimate, scaleA1, scaleA2, n, m, alpha) ⇒ Object

Two-sided asymptotic confidence interval for H.

Combines the asymptotic standard error with the standard-normal critical value z_{1 - alpha/2} (computed by the internal normalQuantile) to produce

CI = H_hat +/- z * SE.

Note: this CI is not clipped to [0, 1]. For practical reporting users may want to clamp to [hMin, hMax].

runKalmanFilter(observations, opts) ⇒ Object

One-dimensional Kalman filter for H(t) smoothing.

State: x_t = H_t. Transition: H_t = H_{t-1} + w_t, w_t ~ N(0, q). Observation: z_t = H_t + v_t, v_t ~ N(0, r).

The filter is seeded with the first observation (x_0 = z_0) and a unit prior covariance. Each subsequent step performs:

1. Predict: xPred = x, pPred = p + q.
2. Update: K = pPred / (pPred + r), x = xPred + K * (z - xPred), p = (1 - K) * pPred.

The result captures both the one-step-ahead predictions (before incorporating the observation) and the filtered states (after).

normalQuantile(p) ⇒ number

Inverse standard normal CDF (quantile function).

Implementation: piecewise rational approximation due to Beasley & Springer (1977) / Acklam (2010). The central region p in [pLow, 1 - pLow] uses a degree-5/4 rational function of r2 = (p - 0.5)^2; the tails use a degree-3/3 rational function of q = sqrt(-2 ln p) (or q = sqrt(-2 ln (1 - p)) for the upper tail).

• p <= 0 returns -Infinity.
• p >= 1 returns Infinity.
• p === 0.5 returns exactly 0.

Numerical accuracy is ~1e-9 across the open interval (0, 1).

normalCdf(x) ⇒ number

Standard normal CDF via the Abramowitz & Stegun rational approximation (7.1.26).

Numerical accuracy is ~7.5e-8 over the whole real line. This is the inverse-of-complement of normalQuantile and is shared by every inference routine that needs a closed-form normal tail probability (currently the constancy likelihood-ratio test in inference/filtering.js).

setLogLevel(level)

Sets the current log level.

getLogLevel() ⇒ number

Reads the current log level.

log(level, label, args)

Internal dispatcher: drops the message if it falls below the configured cut-off, otherwise forwards to the appropriate console.* channel.

debug(...args)

Emits a message at DEBUG level.

info(...args)

Emits a message at INFO level.

warn(...args)

Emits a message at WARN level (visible by default).

error(...args)

Emits a message at ERROR level (visible by default).

getModel(name) ⇒ StochasticModel | undefined

Retrieves a registered model strategy by name.

registerModel(name, factory)

Registers a new model strategy under the supplied name.

listModels() ⇒ Array.<string>

Lists every registered model strategy identifier.

getForecaster(name) ⇒ Forecaster | undefined

Retrieves a registered forecaster by name.

registerForecaster(name, factory)

Registers a new forecaster strategy under the supplied name.

listForecasters() ⇒ Array.<string>

Lists every registered forecaster identifier.

runAdaptiveGridSearch(f, min, max, opts) ⇒ Object

Adaptive grid search with Brent refinement for 1D minimization.

Algorithm:

1. Initialize with the midpoint of [min, max].
2. Repeat refineIters times:
   • Sample gridSize evenly spaced points across [a, b].
   • Track the best point.
   • Shrink [a, b] to [best - 2*step, best + 2*step] clamped to the original interval.
   • Stop early if [a, b] shrinks below tol.
3. Polish the local minimum with Brent's method using bestX as the initial guess.

The Brent refinement makes the function value at the returned x accurate to machine epsilon in nearly all cases.

runBrent(f, ax, bx, cx, tol) ⇒ Object

Minimizes f(x) on the interval [ax, cx] using Brent's method.

The algorithm tracks the best point x, the second-best w, and the third-best v; it uses a parabolic fit whenever the parabolic step is safe, otherwise falls back to a golden-section step. Convergence is declared when |x - midpoint| <= 2 * tol * |x| + EPS or when the iteration cap of 100 is reached.

Invariants:

• The bracket [a, b] always contains the minimum.
• f(x) <= f(w) <= f(v) at every iteration.

runDifferentialEvolution(f, x0, opts) ⇒ Object

Differential-evolution minimization over an arbitrary-dimensional space.

The initial population is drawn uniformly inside [lb, ub]. Each member produces one trial per generation; the trial survives to the next generation only when its objective is strictly better.

runNelderMead(f, x0, opts) ⇒ Object

Nelder-Mead minimization over a multidimensional space.

Builds an initial simplex by perturbing each axis of x0 by 1e-4 and then iterates the standard reflection / expansion / contraction / shrink move until either the spread of function values is below tol or maxIter iterations have been performed.

runSimulatedAnnealing(f, x0, opts) ⇒ Object

Simulated-annealing minimization over an arbitrary-dimensional space.

The neighbor for each iteration is generated by perturbing every coordinate by a uniform offset in [-stepSize, stepSize]. The acceptance temperature decays geometrically: temp *= coolingRate. The loop terminates once either maxIter iterations are performed or the temperature drops below finalTemp.

mulberry32(seed) ⇒ function

Constructs a mulberry32 generator with the given 32-bit seed.

The algorithm packs the state into a single unsigned 32-bit integer a. Each call applies two well-known integer mixing steps (Math.imul & bitwise shift) and returns the result divided by 2^32 so the output is in [0, 1).

setRandomSeed(seed)

Sets a global seed for reproducible simulations.

Passing null or undefined clears the seed and reverts to Math.random(). Calling setRandomSeed twice restarts the deterministic sequence from scratch.

resetRandomSeed()

Resets the PRNG to use Math.random() for all subsequent draws.

Equivalent to setRandomSeed(null). Use this at the end of a deterministic experiment to restore nondeterministic behavior.

nextRandom() ⇒ number

Returns a uniform random number in [0, 1).

Uses the seeded generator when one has been installed via setRandomSeed, otherwise falls through to Math.random(). Because this dispatcher is called from every stochastic primitive in the library, the entire computation tree is reproducible from a single seed.

computeKsDistance(sample1, sample2, isSorted) ⇒ number

Computes the two-sample Kolmogorov-Smirnov distance.

Algorithm: a linear merged-pointer walk over the sorted order statistics. As we walk through the sorted union we maintain the empirical CDF values F_n(x) = (i + 1) / n and G_m(x) = j / m at the current position and record the absolute difference. Sorting first dominates the cost; the walk itself is O(n + m) where n = sample1.length and m = sample2.length.

Input validation:

• Both samples must be non-empty arrays or Float64Arrays.
• All values must be finite (no NaN, +Infinity, -Infinity).

Ties: when values are equal the walk advances both pointers and uses (i + 1) / n vs. (j + 1) / m for the distance — this matches the standard two-sided statistic.

computeKsDistanceRescaled(sortedA, sortedB, factorA, factorB) ⇒ number

Kolmogorov-Smirnov distance for already sorted samples that need rescaling.

Equivalent to ksDistance(a, b, true) but applies the rescaling factors during the merged-pointer walk so no auxiliary allocation is needed. Multiplication by a positive scalar is order-preserving, so the pre-sorting of the inputs is unaffected by the choice of factorA and factorB.

This is the hot path of the Hurstify estimator's inner loop: O(n + m) per evaluation, no allocations beyond the locals below.

shuffleArray(array) ⇒ Array.<*>

Unbiased Fisher-Yates shuffle.

Returns a new array; the input is never mutated. Uses the seeded PRNG exposed by prng.js, so the result is reproducible when a seed is set.

Complexity: O(n) time, O(n) extra memory.

permuteBlocks(data, blockSize, randomPhase) ⇒ Array.<*>

Block random permutation for decorrelating serial dependence.

Conceptually this is the paper's "preserves marginals, kills short-range autocorrelation" operation:

1. (Optional) shift the starting index by a uniform [-0, blockSize) offset so two calls with the same seed still produce different alignments.
2. Slice the resulting series into blocks of length blockSize (the first block may be shorter than blockSize when a phase offset was applied).
3. Apply a Fisher-Yates shuffle to the block list.
4. Concatenate the shuffled blocks back into a single sequence.

Picking blockSize is the user's responsibility: it should be larger than the dominant autocorrelation length in data. Too small and serial dependence survives; too large and the number of blocks — and therefore the effective randomization — shrinks.

getRandomSample(array, n) ⇒ Array.<*>

Floyd's Algorithm R reservoir sampler.

Streams over the input producing a uniformly random sample of size n without replacement. Equivalent to shuffle(array).slice(0, n) but uses only O(n) auxiliary memory and a single pass through array, which matters when sampling from very large arrays (e.g. millions of increments).

Edge cases:

• n >= array.length: returns a shuffled full copy of array.
• n <= 0: returns an empty array.

nextGaussian() ⇒ number

Draws a single standard normal via Box-Muller.

The polar variant is implemented by guarding against degenerate u === 0 draws from nextRandom(). One Box-Muller pair yields two independent standard normals; this routine keeps the cosine component and discards the sine. Use generateCorrelatedGaussian if you need both halves, or call nextGaussian twice with distinct nextRandom() outputs.

generateGaussianBatch(n) ⇒ Float64Array

Pre-allocates a Float64Array of standard normals.

Useful when an inner loop needs a contiguous buffer of normals; the allocation is amortized across a single batch draw, whereas repeated nextGaussian calls would each allocate internally.

generateCorrelatedGaussian(n, rho) ⇒ Array.<Float64Array>

Generates two correlated standard-normal streams via Cholesky.

Mathematically the model is (Z1, Z2) with unit marginals and Corr(Z1, Z2) = rho. Implementation: draw an i.i.d. Box-Muller pair (z1, z2); set Z1 = z1; set Z2 = rho * z1 + sqrt(1 - rho^2) * z2. Both Z1 and Z2 have unit variance and exactly correlation rho.

Important: rho must be strictly in (-1, 1); the implementation silently clamps 1 - rho^2 to zero via Math.max(0, ...) so the endpoints collapse to the trivial deterministic case.

generateFractionalNoise(n, H) ⇒ Float64Array

Fractional Gaussian Noise via Hosking's method.

Hosking's method is an exact O(n^2) Cholesky-style recursion that generates samples from the autocovariance gamma(k) = 0.5 (|k+1|^{2H} - 2|k|^{2H} + |k-1|^{2H}).

It uses O(n) recursion updates to compute the conditional mean and variance (phi, v) incrementally, so the per-step cost is O(k) and the total O(n^2). This is fine for the scales used in the paper (a few hundred to a few thousand samples) but dominates for n >> 1e4.

Assumptions:

• n > 0 and H in (0, 1).
• The result is mean-zero (the recursion conditions on x_0 ~ N(0, 1)).

generateFractionalBrownianMotion(n, H) ⇒ Float64Array

Fractional Brownian Motion by cumulative summation of fGN.

The implementation delegates the heavy lifting to generateFractionalNoise and then performs a single O(n) cumulative-sum pass. The first sample is fixed at 0 (the standard convention for fBM(0) = 0), so paths always start at the origin.

For non-zero means, simply add a constant afterwards — fGn is mean-zero by construction.

computeFractionalKernel(H, nSteps, dt) ⇒ Float64Array

Precomputes the Riemann-Liouville fractional kernel used by the rough-volatility simulators.

Mathematically K(t) = sqrt(2 H) * t^{H - 0.5} for t > 0. The result is a length-nSteps array where entry i corresponds to t = (i + 1) * dt.

Reusing a precomputed kernel for every path avoids the O(n^2) cost of re-evaluating the power function per integration step.

computeFractionalIntegral(dW, kernel, t) ⇒ number

Computes a single time-step of the Riemann-Liouville fractional integral.

Given precomputed Brownian increments dW and a kernel from computeFractionalKernel, returns I_t = sum_{j=0}^{t-1} K(t - j) * dW_j.

Used inside the rBergomi path generator and the exact fOU driver.

Complexity: O(t) per call, so building a full path is O(n^2). This is acceptable for paths up to a few hundred steps; for long simulations switch to a circulant-embedding FFT approximation (not implemented here).

xavierInit(rows, cols) ⇒ Array.<Array.<number>>

Xavier (Glorot-uniform) weight initialization.

Produces a rows x cols matrix where each entry is sampled uniformly in [-scale, scale] with scale = sqrt(2 / (rows + cols)). This is the standard initializer for tanh/sigmoid-activated layers (Glorot & Bengio, 2010).

getBinomialCoeffs(d, lag) ⇒ Float64Array

Returns the binomial coefficient sequence [C(d, 0), ..., C(d, lag)].

Uses a tiny FIFO cache keyed by ${d}:${lag} so that identical lookups within a rolling ARFIMA run are O(1). When the cache is full the oldest entry is evicted.

fractionalDifference(data, d, [lag]) ⇒ Array.<number>

Computes the (truncated) fractional difference of a series for a given d and lag cap. The truncation to lag keeps the per-step cost O(lag) rather than O(t), which is essential for long-history forecasting.

ksCriticalValue(n, m, alpha) ⇒ number

Two-sample Kolmogorov–Smirnov asymptotic critical value.

D_alpha = sqrt(-0.5 * ln(alpha / 2)) * sqrt((n + m) / (n * m))

ksPvalue(D, n, m) ⇒ number

Approximate two-sample KS p-value via the asymptotic Kolmogorov distribution.

Q(lambda) ~ 2 * sum_{j=1..3} (-1)^{j-1} * exp(-2 j^2 lambda^2)

with the standard lambda correction.

kalmanLogLikelihood(observations, q, r) ⇒ number

Log-likelihood of the observations under a 1D Kalman filter.

detectCusumBreakpoints(hHistory, windowSize, threshold) ⇒ Array.<{index: number, H_before: number, H_after: number}>

Detects breakpoints in a series of H estimates via a sliding-window CUSUM.

chooseKsObjective([scales], [weights]) ⇒ KsObjective

Selects the right KsObjective for a configuration.

defaultSampler([blockSize]) ⇒ Sampler

Convenience: selects the default sampler based on blockSize.

• When blockSize is a positive number a BlockPermutationSampler is returned.
• Otherwise a ReservoirSampler is returned.

Typedefs

KsSignificanceResult : Object

ConstancyResult : Object

CusumBreakResult : Object

BootstrapCiResult : Object

SimulationResult : Object

PriceResult : Object

Optimizer

Strategy base class. Subclasses implement minimize polymorphically.

Kind: global class

optimizer.minimize(_objective, _lower, _upper, _initial) ⇒ number

Minimizes the given objective on the closed interval [lower, upper] starting from initial. Subclasses implement the algorithm-specific search.

Kind: instance method of Optimizer
Returns: number - Argmin h inside [_lower, _upper].

Param	Type	Description
_objective	function	Scalar objective function.
_lower	number	Lower bound of the search interval.
_upper	number	Upper bound of the search interval.
_initial	number	Initial guess inside [_lower, _upper].

Forecaster

Abstract base class for H-series forecasters.

Kind: global abstract class

• Forecaster
  • .predict(history) ⇒ number
  • .forecast(history) ⇒ number

forecaster.predict(history) ⇒ number

Predicts the next H value from the supplied history.

Kind: instance abstract method of Forecaster
Returns: number - Predicted H.

Param	Type	Description
history	Array.<number>	Time-ordered H estimates.

forecaster.forecast(history) ⇒ number

Alias for predict.

Kind: instance method of Forecaster
Returns: number - Predicted H.

Param	Type	Description
history	Array.<number>	Time-ordered H estimates.

HypothesisTest

Kind: global abstract class

hypothesisTest.run(data, opts) ⇒ Result

Runs the test on data with the supplied options.

Kind: instance abstract method of HypothesisTest
Returns: Result - Test result bundle.

Param	Type	Description
data	*	Test input (depends on the concrete test).
opts	Object	Test options.

Kernel

Abstract base class for fractional-integration kernels.

Kind: global abstract class

• Kernel
  • .evaluate(t) ⇒ number
  • .precompute(nSteps, dt) ⇒ Float64Array

kernel.evaluate(t) ⇒ number

Evaluates the kernel at lag t > 0.

Kind: instance abstract method of Kernel
Returns: number - Kernel weight at t.

Param	Type	Description
t	number	Positive lag (units of dt).

kernel.precompute(nSteps, dt) ⇒ Float64Array

Precomputes the kernel over nSteps time steps at stride dt.

Kind: instance method of Kernel
Returns: Float64Array - Cached kernel values.

Param	Type	Description
nSteps	number	Number of time steps.
dt	number	Per-step time increment.

KsObjective

Abstract base class for KS-distance objectives.

Kind: global abstract class

ksObjective.evaluate(sortedSamples, scales, H) ⇒ number

Evaluates the objective at the trial H.

Kind: instance abstract method of KsObjective
Returns: number - Non-negative objective value.

Param	Type	Description
sortedSamples	Array.<Float64Array>	One pre-sorted sample per scale (each of equal length).
scales	Array.<number>	Scale values matching sortedSamples.
H	number	Trial Hurst parameter.

MultiScaleKsObjective

Multi-scale unweighted KS objective.

Returns the arithmetic mean of the pairwise KS distances over every unordered pair (i, j) with i < j. Equivalent to the paper's recommended extension to K > 2 scales.

Kind: global class

multiScaleKsObjective.evaluate()

Kind: instance method of MultiScaleKsObjective

StochasticModel

Abstract base class for stochastic-process simulators.

Kind: global abstract class

• StochasticModel
  • .simulate(opts) ⇒ SimulationResult
  • .price(sim, opts) ⇒ PriceResult

stochasticModel.simulate(opts) ⇒ SimulationResult

Simulates nPaths paths of the underlying stochastic process.

Kind: instance abstract method of StochasticModel
Returns: SimulationResult - Simulated paths + time grid.

Param	Type	Description
opts	Object	Model-specific options.

stochasticModel.price(sim, opts) ⇒ PriceResult

Optionally drives a price SDE using the simulator's noise realization. Default: throws — only models with a price SDE implement this method.

Kind: instance method of StochasticModel
Returns: PriceResult - Simulated prices.

Param	Type	Description
sim	SimulationResult	Output of simulate.
opts	Object	Price-SDE options.

RoughBergomiModel

Rough Bergomi model.

dV_t / V_t = eta * dW^perp_t
I_t = int_0^t sqrt(2H) (t - s)^{H - 0.5} dW^perp_s
V_t = xi * exp(eta I_t - (eta^2 / 2) t^{2H})

Kind: global class

RoughFsvModel

Rough Fractional Stochastic Volatility model.

dV_t = theta (mu - V_t) dt + nu V_t^alpha dW^V_t + roughComp

Kind: global class

FractionalOuModel

Abstract base class for the Fractional Ornstein-Uhlenbeck model.

dX_t = theta (mu - X_t) dt + sigma dB^H_t

Concrete subclasses pick the discretization scheme:

• EulerMaruyamaFractionalOuModel — default, cheap
• ExactFractionalOuModel — Riemann-Liouville integral, accurate

Kind: global abstract class

• FractionalOuModel
  • instance
    • .vasicekPath(opts) ⇒ Object
  • static
    • .parseOpts([opts]) ⇒ Object
    • .buildTimes(nSteps, dt) ⇒ Array.<number>

fractionalOuModel.vasicekPath(opts) ⇒ Object

Special-case fast path: exact Vasicek recursion when H = 0.5.

Kind: instance method of FractionalOuModel
Returns: Object - Generated path and matching time grid.

Param	Type	Description
opts	Object	Model options.

FractionalOuModel.parseOpts([opts]) ⇒ Object

Shared parameter parsing for the fOU family.

Kind: static method of FractionalOuModel
Returns: Object - Normalized parameter bundle.

Param	Type	Description
[opts]	Object	Model options.

FractionalOuModel.buildTimes(nSteps, dt) ⇒ Array.<number>

Shared time-grid construction for the fOU family.

Kind: static method of FractionalOuModel
Returns: Array.<number> - Array of length nSteps + 1 of cumulative times.

Param	Type	Description
nSteps	number	Number of discrete steps.
dt	number	Time step size.

EulerMaruyamaFractionalOuModel

Euler-Maruyama discretization of the fOU model.

Cheap O(n) integration; first-order accurate.

Kind: global class

ExactFractionalOuModel

Exact Riemann-Liouville discretization of the fOU model.

O(n^2) per path; higher-order accurate.

Kind: global class

MultifractionalPreModel

Abstract base class for the Multifractional Process with Random Exponent.

X_t = B_{H(t)}(t)

where H(t) itself is a stochastic Ornstein-Uhlenbeck process bounded between hMin and hMax. Two concrete subclasses pick the discretization scheme:

• LocalHolderMultifractionalPreModel — default, cheap
• ExactMultifractionalPreModel — time-varying kernel

Kind: global abstract class

MultifractionalPreModel.generateHPath(opts) ⇒ Array.<number>

Shared H(t) path generation under an OU bridge.

Kind: static method of MultifractionalPreModel
Returns: Array.<number> - Generated H(t) path of length nSteps + 1.

Param	Type	Description
opts	Object	Model options.

LocalHolderMultifractionalPreModel

Local-Holder approximation of the mPRE model.

Cheap O(n) integration via cumulative sqrt(dt^{2 * H_avg}) scaling.

Kind: global class

ExactMultifractionalPreModel

Exact time-varying-kernel discretization of the mPRE model.

O(n^2) per path; uses a time-varying Riemann-Liouville kernel.

Kind: global class

Sampler

Abstract base class for sampling strategies.

A Sampler is a strategy: callers obtain a fresh instance and invoke draw(inc, n) once per variance-reduction iteration. The estimator never holds sampler state across calls so strategies can be safely shared across estimator instances.

Kind: global abstract class

sampler.draw(inc, n) ⇒ Array.<number>

Draws a sub-sample of size n from inc.

Kind: instance abstract method of Sampler
Returns: Array.<number> - Sampled sub-array.

Param	Type	Description
inc	Array.<number> | Float64Array	Increment array at a single scale.
n	number	Desired sample size.

HurstifyErrorCode : enum

Stable error codes. Treat the string values as part of the public API — renaming them is a breaking change.

Kind: global enum

LogLevel : enum

Severity levels, numerically ordered from most to least verbose.

• DEBUG (0): per-step diagnostics, only useful for tracing algorithm internals.
• INFO (1): high-level progress messages.
• WARN (2): recoverable issues (default cut-off).
• ERROR (3): unhandled failures during processing.
• SILENT (4): disables all logging; convenience for tests.

Kind: global enum

NORMAL_QUANTILE_COEFFS

Coefficients for the Beasley-Springer-Malkin rational approximation of the inverse standard normal CDF. Used piecewise for p in [pLow, 1 - pLow] (central region) and tail rational functions for the extremes.

The standard deviation may be c/d constants at the tails is adapted from Peter Acklam's algorithm.

Kind: global constant

modelRegistry : Registry.<StochasticModel>

Strategy registry for stochastic models. The default fOU key resolves to the Euler-Maruyama discretization; consumers who want the exact Riemann-Liouville variant look up fOU-exact. Same convention for mPRE / mPRE-exact.

Kind: global constant

forecasterRegistry : Registry.<Forecaster>

Strategy registry for forecasters.

Kind: global constant

optimizerRegistry : Registry.<Optimizer>

Global optimizer registry.

Kind: global constant

parseCsv(csv, opts) ⇒ Array.<Object>

Parses a CSV string into an array of plain objects.

Expected input shape:

• The first non-empty line is the header row.
• Each subsequent line is a record with the same column count as the header.
• Fields can be optionally wrapped in double quotes; quotes may embed commas but not other escapes.

Type coercion:

• opts.dateField (default "date") is parsed via new Date(...).
• Any field listed in opts.numericFields is parsed via parseFloat.
• All other fields are kept as trimmed strings.

Error handling:

• Empty input returns [].
• Mismatched column counts throw with a descriptive message.
• Non-numeric values in declared numeric columns throw.

Kind: global function
Returns: Array.<Object> - Parsed rows, one per non-empty CSV line.
Throws:

• Error When the input is malformed.

Param	Type	Description
csv	string	Raw CSV content.
opts	Object	Parser options.
[opts.dateField]	string	Date column name (default "date").
[opts.numericFields]	Array.<string>	Columns to coerce to numbers.

splitCSVLine(line) ⇒ Array.<string>

Splits a single CSV line respecting double-quoted regions.

States:

• Outside quotes: a comma terminates the current field.
• Inside quotes: a quote toggles back to "outside", all other chars are kept verbatim.

Kind: global function
Returns: Array.<string> - Split fields. Empty trailing field is preserved.

Param	Type	Description
line	string	Raw CSV line (no trailing newline).

extractSeries(rows, field, opts) ⇒ Array.<{date: Date, value: number}>

Extracts a {date, value} series from a parsed CSV array.

Rows that are missing field are skipped; the resulting series is optionally sorted by dateField when the caller asks. Sorting uses the standard JS Date arithmetic, so the dates must be real Date instances.

Kind: global function
Returns: Array.<{date: Date, value: number}> - Series of {date, value} points.
Throws:

• Error When rows is not an array or field is not a string.

Param	Type	Description
rows	Array.<Object>	Parsed CSV rows.
field	string	Numeric field name to extract.
opts	Object	Extraction options.
[opts.sortByDate]	boolean	When true, sort by the date field before extraction (default false).
[opts.dateField]	string	Date field name (default "date").

parseJson(json) ⇒ Array.<Object>

Parses a JSON string that must encode an array of objects.

The function deliberately refuses non-array JSON to keep the loader simple. Empty or whitespace-only input returns [].

Kind: global function
Returns: Array.<Object> - Parsed objects (empty if the input is empty).
Throws:

• Error When the input is not valid JSON or does not decode to an array.

Param	Type	Description
json	string	Raw JSON string.

validateNoGaps(series, maxGapMs) ⇒ Object

Validates that a time series does not contain temporal gaps larger than maxGapMs.

Returns the maximum observed gap, the full list of pairwise gap lengths, and a valid flag for the threshold check. Series with fewer than two points are deemed valid by definition.

Kind: global function
Returns: Object - Validation result.

Param	Type	Description
series	Array.<{date: Date}>	Time series with Date fields.
maxGapMs	number	Maximum allowed gap in milliseconds.

downsampleSeries(series, intervalMs) ⇒ Array.<{date: Date, value: number}>

Downsamples a time series by averaging values that fall into fixed intervalMs-wide buckets.

The bucket index is computed as floor(date.getTime() / intervalMs), so all buckets share the same left edge (0, intervalMs, 2 * intervalMs, ...). The output is sorted by date and every returned point carries the bucket start (not the average timestamp) as its date value.

Kind: global function
Returns: Array.<{date: Date, value: number}> - One entry per non-empty bucket, sorted chronologically.
Throws:

• Error When series is not an array or intervalMs <= 0.

Param	Type	Description
series	Array.<{date: Date, value: number}>	Input series.
intervalMs	number	Bucket length in milliseconds.

preaverageReturns(prices, [windowSize]) ⇒ Array.<number>

Preaveraging of log-returns.

Implementation of the Jacod et al. (2009) preaveraging estimator (simplified single-bar variant):

1. Compute log-returns r_t = log(P_t / P_{t-1}).
2. For each i, average the windowSize consecutive returns ending at i (g_avg[i] = mean(r_{i - windowSize + 1}, ..., r_i)).
3. The "preaveraged return" is the first-difference sequence g_avg[i] - g_avg[i - 1]. This cancellation attenuates microstructure noise by 1/sqrt(windowSize) while preserving the drift and diffusion up to O(1 / windowSize).

Note: the result has length prices.length - windowSize - 1; for very short series the function throws rather than returning a few noisy points.

Kind: global function
Returns: Array.<number> - Preaveraged-returns series.
Throws:

• Error When prices has fewer than windowSize + 1 elements.

Param	Type	Description
prices	Array.<number>	Price series.
[windowSize]	number	Preaveraging window (default 2).

computeRealizedKernel(returns, [kernelType], [bandwidth]) ⇒ number

Realized-kernel variance estimator with pluggable kernels.

Given n returns, the estimator forms the autocorrelation sequence

gamma_k = sum_{i=k+1}^{n} r_i * r_{i - k},  k = 0..h

and combines them through a weighted sum

RV_K = gamma_0 + 2 * sum_{k=1..h} w_k * gamma_k

with weights w_k provided by the chosen kernel. The default bandwidth is floor(n^0.6), a rule-of-thumb that matches the optimal scaling under i.i.d. microstructure noise.

Kernels shipped:

• bartlett: w_k = 1 - k / h (default).
• parzen: the standard piecewise-cubic Parzen kernel.
• tukey-hanning: 0.5 (1 + cos(pi k / h)).

Any unknown kernel name falls back to Bartlett.

Kind: global function
Returns: number - Realized-kernel variance (clamped to be non-negative).
Throws:

• Error When returns is empty.

Param	Type	Description
returns	Array.<number>	Log-return series.
[kernelType]	string	One of "bartlett", "parzen", "tukey-hanning" (default "bartlett").
[bandwidth]	number	Optional explicit bandwidth; defaults to floor(n^0.6).

kernelWeight(type, k, h) ⇒ number

Kernel weight function used by realizedKernel.

Kind: global function
Returns: number - Weight for the k-th autocorrelation lag.

Param	Type	Description
type	string	Kernel identifier ("bartlett", "parzen", "tukey-hanning").
k	number	Lag index (k >= 0).
h	number	Bandwidth (h > 0).

debiasLogVolatility(rawHEstimates, sigmaObs, sigmaLatent) ⇒ Array.<number>

Heuristic de-biasing of log-volatility H estimates.

Microstructure noise inflates the variance of the log-volatility proxy relative to the latent signal, which in turn attenuates the observed roughness. This routine adds a small correction

h_debias = h + 0.01 * log(sigmaObs / sigmaLatent)

and clamps the result to [0.01, 0.99]. It is intentionally conservative — the user is expected to validate the calibration against a trust sample before relying on it for production.

Kind: global function
Returns: Array.<number> - De-biased H estimates.
Throws:

• Error When sigmaLatent <= 0.

Param	Type	Description
rawHEstimates	Array.<number>	Raw H estimates from Hurstify.estimate or rolling.
sigmaObs	number	Standard deviation of the observed log-vol series.
sigmaLatent	number	Standard deviation of the latent (denoised) log-vol series.

computeRealizedVariance(prices, [interval]) ⇒ Array.<number>

Computes per-bucket realized variance from a price series.

The realized variance is the sum of squared log-returns within each non-overlapping bucket of interval observations:

RV_k = sum_{i in bucket k} (log P_i - log P_{i-1})^2

With interval = 1 the function emits one RV per log-return directly, which is the canonical "5-minute RV" form when prices are already sampled at 5-minute intervals.

Kind: global function
Returns: Array.<number> - Realized-variance series.
Throws:

• Error When prices is missing, has fewer than two elements, contains non-finite or non-positive values, or interval is not a positive integer.

Param	Type	Description
prices	Array.<number>	Chronological price series (strictly positive, finite).
[interval]	number	Bucket size (default 1; must be a positive integer).

computeRealizedVarianceParkinson(bars) ⇒ Array.<number>

Parkinson (1980) high-low RV estimator from OHLC bars.

For each bar the within-period variance is approximated by

sigma^2 ~= (log(H/L))^2 / (4 * ln 2)

which is 1/(4 ln 2) ~ 0.36 of the log-range-squared. Parkinson is strictly less efficient than tick-based RV but only requires four numbers per bar.

Kind: global function
Returns: Array.<number> - One Parkinson variance estimate per bar.
Throws:

• Error When bars is not an array or any bar has non-positive/non-finite high/low values.

Param	Type	Description
bars	Array.<{open: number, high: number, low: number, close: number}>	OHLC bars.

aggregateDailyRealizedVariance(intradayRVs) ⇒ number

Aggregates intraday (5-minute) realized variances into a single daily value via plain summation.

This is the standard "sum of squared returns" daily RV used in financial econometrics. It assumes the input is already free of overnight gaps.

Kind: global function
Returns: number - Sum of the intraday RVs (zero for an empty input).
Throws:

• Error When intradayRVs is not an array.

Param	Type	Description
intradayRVs	Array.<number>	Sequence of 5-minute RVs.

applyLogTransform(rv) ⇒ Array.<number>

Maps realized variance to the log-volatility series consumed by hurstify.

The transformation is

X_t = 0.5 * log(RV_t)

i.e. log(sqrt(RV)). This converts multiplicative variance dynamics into a roughly additive (and therefore more stationary) signal, on top of which the self-similarity property exploited by the RK-SAVR algorithm is expressed.

Kind: global function
Returns: Array.<number> - Log-volatility series.
Throws:

• Error When rv is not an array or contains non-positive / non-finite values.

Param	Type	Description
rv	Array.<number>	Realized-variance series.

centerSeries(series) ⇒ Array.<number>

Subtracts the arithmetic mean from every element.

Useful as a final step in the preprocessing pipeline when the user wants the series to mean-zero (which can stabilize variance-reducing permutations inside Hurstify).

Kind: global function
Returns: Array.<number> - New array of length series.length with the mean subtracted. Empty input yields [].

Param	Type	Description
series	Array.<number>	Input series.

standardizeSeries(series) ⇒ Array.<number>

Standardizes a time series to zero mean and unit variance.

Divides each centered value by the population standard deviation. A constant series has zero variance and triggers an explicit error rather than silently producing NaNs.

Kind: global function
Returns: Array.<number> - Standardized copy of series.
Throws:

• Error When series has fewer than two elements or population variance zero.

Param	Type	Description
series	Array.<number>	Input series (needs at least two points).

applyPreprocessingPipeline(prices, opts) ⇒ Array.<number>

Bundled preprocessing pipeline: prices -> RV -> log-vol -> (optional) centering.

Equivalent to running computeRV + logTransform + (optionally) centerSeries, but more compact for callers who want the canonical transformation.

Kind: global function
Returns: Array.<number> - Preprocessed log-volatility series.

Param	Type	Description
prices	Array.<number>	Chronological price series.
opts	Object	Pipeline options.
[opts.interval]	number	RV aggregation interval (default 1).
[opts.center]	boolean	When true, subtract the mean from the log-volatility series at the end (default false).

splitTrainTest(series, [trainRatio]) ⇒ Object

Splits a series into contiguous training and test arrays.

The split point is floor(series.length * trainRatio) so the training set is the leftmost prefix of the series; this preserves temporal ordering, which is what hurstify forecasters and validation scripts typically need.

Kind: global function
Returns: Object - Train/test arrays.
Throws:

• Error When series is not an array or trainRatio is out of range.

Param	Type	Description
series	Array.<number>	Input series.
[trainRatio]	number	Training fraction in (0, 1) (default 0.8).

createSlidingWindows(series, windowSize, [step]) ⇒ Array.<Array.<number>>

Builds overlapping windows from a single time series.

The i-th window is series.slice(i, i + windowSize) for i = 0, step, 2*step, ... until no full window fits. Used by offline batch evaluation pipelines that want to score the estimator on every available segment of the series.

Kind: global function
Returns: Array.<Array.<number>> - One entry per non-truncated window.
Throws:

• Error When series is not an array or windowSize/step are non-positive.

Param	Type	Description
series	Array.<number>	Input series.
windowSize	number	Window length (positive integer).
[step]	number	Stride between consecutive windows (default 1).

generateVixLogVolatility(nDays, h, opts) ⇒ Array.<number>

Synthetic VIX-style daily log-volatility.

Generates an fBM with the requested h and maps it to a log-volatility level around 2.0 (i.e. sqrt(RV) ~ 20%) by adding a small drift term and Gaussian observation noise:

X_t = 2.0 + drift * (fbm[t] / sqrt(n)) + 0.5 * fbm[t] + noise

Default tuning matches the empirical VIX roughness (h ~ 0.1) and annualized log-vol mean.

Kind: global function
Returns: Array.<number> - Daily log-volatility series. Empty when nDays <= 0.
Throws:

• Error When h is out of (0, 1).

Param	Type	Description
nDays	number	Number of trading days.
h	number	Hurst parameter (default 0.1).
opts	Object	Generation options.
[opts.seed]	number	PRNG seed for reproducibility.
[opts.noiseStd]	number	Observation-noise standard deviation (default 0.05).
[opts.drift]	number	Log-volatility drift (default 0.02).

generateSpxLogVolatility(nDays, h, opts) ⇒ Array.<number>

Synthetic S&P 500 realized-volatility style daily log-volatility.

Same construction as generateVIXLogVol but with a smoother default Hurst (h = 0.14), a smaller drift, and a less volatile observation-noise level. Empirically these choices match the rough regime typically reported for SPX RV.

Kind: global function
Returns: Array.<number> - Daily log-volatility series. Empty when nDays <= 0.
Throws:

• Error When h is out of (0, 1).

Param	Type	Description
nDays	number	Number of trading days.
h	number	Hurst parameter (default 0.14).
opts	Object	Generation options.
[opts.seed]	number	PRNG seed.
[opts.noiseStd]	number	Observation-noise standard deviation (default 0.03).
[opts.drift]	number	Log-volatility drift (default 0.015).

generateIntradayPrices([nIntraday], [nDays], h, opts) ⇒ Array.<Array.<number>>

Generates synthetic intraday 5-minute prices useful for testing realized-variance pipelines.

For every (re-)sampled day the generator draws an fBM with the requested h, exponentiates it into a volatility factor, and steps a log-return process

S_{i+1} = S_i * exp(drift + vol_i * z_i * sqrt(dt))

with drift set to the per-5-minute-bar annualized drift. The result is a nDays x nIntraday array of prices suitable for feeding into computeRV.

Kind: global function
Returns: Array.<Array.<number>> - Array of daily price arrays.
Throws:

• Error When nIntraday <= 0, nDays <= 0, or h is out of (0, 1).

Param	Type	Description
[nIntraday]	number	Number of 5-minute bars per day (default 78, the typical US-equities count).
[nDays]	number	Number of days to simulate (default 1).
h	number	Hurst parameter.
opts	Object	Generation options.
[opts.seed]	number	PRNG seed for reproducibility.
[opts.drift]	number	Annualized drift (default 0.05).

seriesToCsv(series, [dateHeader], [valueHeader]) ⇒ string

Serializes a {date, value} series as a CSV string.

Dates that are Date instances are formatted as their ISO yyyy-mm-dd prefix; everything else is stringified verbatim. Empty series produces a header-only CSV.

Kind: global function
Returns: string - CSV-encoded content joined with \n.
Throws:

• Error When series is not an array.

Param	Type	Description
series	Array.<Object>	Time series with date and value fields.
[dateHeader]	string	Header for the date column (default "date").
[valueHeader]	string	Header for the value column (default "value").

buildScaleProfile(sortedSamples, scales, H) ⇒ Array.<number>

Builds a flat "profile" of all pairwise KS distances at a fixed H.

Given K sorted samples, the profile has K * (K - 1) / 2 entries corresponding to every unordered scale pair. Useful for diagnostics.

Kind: global function
Returns: Array.<number> - Flat array of pairwise KS distances.

Param	Type	Description
sortedSamples	Array.<Float64Array>	Pre-sorted samples.
scales	Array.<number>	Scale values.
H	number	Hurst parameter.

getAsymptoticVariance(scaleA1, scaleA2, n, m) ⇒ number

Asymptotic variance of the Hurstify estimator.

Implements

Var(H_hat) = (2 * pi * e) / (ln(a2/a1))^2 * (1/sqrt(n) + 1/sqrt(m))^2.

When a1 == a2 (log ratio zero) the variance is degenerate and the function returns Infinity rather than dividing by zero; callers that intend to compute a SE/CI should reject equal scales up-front.

Kind: global function
Returns: number - Non-negative asymptotic variance (Infinity if the scales coincide).

Param	Type	Description
scaleA1	number	Lower scale a_1.
scaleA2	number	Upper scale a_2.
n	number	Sample size at a_1.
m	number	Sample size at a_2.

getStandardError(scaleA1, scaleA2, n, m) ⇒ number

Asymptotic standard error: square root of the asymptotic variance.

Thin convenience wrapper. The standard error has units of "Hurst" and can be read against the hMin/hMax bounds the estimator was configured with.

Kind: global function
Returns: number - Non-negative standard error (Infinity for degenerate scale choices).

Param	Type	Description
scaleA1	number	Lower scale a_1.
scaleA2	number	Upper scale a_2.
n	number	Sample size at a_1.
m	number	Sample size at a_2.

getConfidenceInterval(hEstimate, scaleA1, scaleA2, n, m, alpha) ⇒ Object

Two-sided asymptotic confidence interval for H.

Combines the asymptotic standard error with the standard-normal critical value z_{1 - alpha/2} (computed by the internal normalQuantile) to produce

CI = H_hat +/- z * SE.

Note: this CI is not clipped to [0, 1]. For practical reporting users may want to clamp to [hMin, hMax].

Kind: global function
Returns: Object - Confidence interval bounds.

Param	Type	Description
hEstimate	number	Point estimate of H.
scaleA1	number	Lower scale a_1.
scaleA2	number	Upper scale a_2.
n	number	Sample size at a_1.
m	number	Sample size at a_2.
alpha	number	Significance level (default 0.05).

runKalmanFilter(observations, opts) ⇒ Object

One-dimensional Kalman filter for H(t) smoothing.

State: x_t = H_t. Transition: H_t = H_{t-1} + w_t, w_t ~ N(0, q). Observation: z_t = H_t + v_t, v_t ~ N(0, r).

The filter is seeded with the first observation (x_0 = z_0) and a unit prior covariance. Each subsequent step performs:

1. Predict: xPred = x, pPred = p + q.
2. Update: K = pPred / (pPred + r), x = xPred + K * (z - xPred), p = (1 - K) * pPred.

The result captures both the one-step-ahead predictions (before incorporating the observation) and the filtered states (after).

Kind: global function
Returns: Object - Filtered and one-step-predicted states, each of length n.

Param	Type	Description
observations	Array.<number>	Time-ordered H estimates.
opts	Object	Filter options.
[opts.q]	number	Process noise variance (default 0.01).
[opts.r]	number	Measurement noise variance (default 0.1).

normalQuantile(p) ⇒ number

Inverse standard normal CDF (quantile function).

Implementation: piecewise rational approximation due to Beasley & Springer (1977) / Acklam (2010). The central region p in [pLow, 1 - pLow] uses a degree-5/4 rational function of r2 = (p - 0.5)^2; the tails use a degree-3/3 rational function of q = sqrt(-2 ln p) (or q = sqrt(-2 ln (1 - p)) for the upper tail).

• p <= 0 returns -Infinity.
• p >= 1 returns Infinity.
• p === 0.5 returns exactly 0.

Numerical accuracy is ~1e-9 across the open interval (0, 1).

Kind: global function
Returns: number - Quantile Phi^{-1}(p).

Param	Type	Description
p	number	Probability in [0, 1].

normalCdf(x) ⇒ number

Standard normal CDF via the Abramowitz & Stegun rational approximation (7.1.26).

Numerical accuracy is ~7.5e-8 over the whole real line. This is the inverse-of-complement of normalQuantile and is shared by every inference routine that needs a closed-form normal tail probability (currently the constancy likelihood-ratio test in inference/filtering.js).

Kind: global function
Returns: number - P(Z <= x) for Z ~ N(0, 1), in [0, 1].

Param	Type	Description
x	number	Input value (any real number).

setLogLevel(level)

Sets the current log level.

Kind: global function

Param	Type	Description
level	number	One of the LogLevel numeric constants.

getLogLevel() ⇒ number

Reads the current log level.

Kind: global function
Returns: number - Active LogLevel value.

log(level, label, args)

Internal dispatcher: drops the message if it falls below the configured cut-off, otherwise forwards to the appropriate console.* channel.

Kind: global function

Param	Type	Description
level	number	Log level (one of LogLevel.*).
label	string	Short human label (DEBUG, INFO, ...).
args	Array.<*>	Arguments to forward to the underlying console.

debug(...args)

Emits a message at DEBUG level.

Kind: global function

Param	Type	Description
...args	*	Values forwarded to console.debug.

info(...args)

Emits a message at INFO level.

Kind: global function

Param	Type	Description
...args	*	Values forwarded to console.info.

warn(...args)

Emits a message at WARN level (visible by default).

Kind: global function

Param	Type	Description
...args	*	Values forwarded to console.warn.

error(...args)

Emits a message at ERROR level (visible by default).

Kind: global function

Param	Type	Description
...args	*	Values forwarded to console.error.

getModel(name) ⇒ StochasticModel | undefined

Retrieves a registered model strategy by name.

Kind: global function
Returns: StochasticModel | undefined - The strategy instance, or undefined when the name is unknown.

Param	Type	Description
name	string	Model identifier.

registerModel(name, factory)

Registers a new model strategy under the supplied name.

Kind: global function

Param	Type	Description
name	string	Unique identifier.
factory	function	Factory returning a fresh instance.

listModels() ⇒ Array.<string>

Lists every registered model strategy identifier.

Kind: global function
Returns: Array.<string> - Snapshot of registered model keys.

getForecaster(name) ⇒ Forecaster | undefined

Retrieves a registered forecaster by name.

Kind: global function
Returns: Forecaster | undefined - The strategy instance.

Param	Type	Description
name	string	Forecaster identifier.

registerForecaster(name, factory)

Registers a new forecaster strategy under the supplied name.

Kind: global function

Param	Type	Description
name	string	Unique identifier.
factory	function	Factory returning a fresh instance.

listForecasters() ⇒ Array.<string>

Lists every registered forecaster identifier.

Kind: global function
Returns: Array.<string> - Snapshot of registered forecaster keys.

runAdaptiveGridSearch(f, min, max, opts) ⇒ Object

Adaptive grid search with Brent refinement for 1D minimization.

Algorithm:

1. Initialize with the midpoint of [min, max].
2. Repeat refineIters times:
   • Sample gridSize evenly spaced points across [a, b].
   • Track the best point.
   • Shrink [a, b] to [best - 2*step, best + 2*step] clamped to the original interval.
   • Stop early if [a, b] shrinks below tol.
3. Polish the local minimum with Brent's method using bestX as the initial guess.

The Brent refinement makes the function value at the returned x accurate to machine epsilon in nearly all cases.

Kind: global function
Returns: Object - Best point and its objective value.
Throws:

• Error When gridSize <= 1.

Param	Type	Description
f	function	Objective function (1D).
min	number	Lower bound.
max	number	Upper bound.
opts	Object	Algorithm options.
[opts.gridSize]	number	Number of coarse-grid points per refinement (default 50).
[opts.refineIters]	number	Number of refinement rounds (default 3).
[opts.tol]	number	Convergence tolerance (default 1e-7).

runBrent(f, ax, bx, cx, tol) ⇒ Object

Minimizes f(x) on the interval [ax, cx] using Brent's method.

The algorithm tracks the best point x, the second-best w, and the third-best v; it uses a parabolic fit whenever the parabolic step is safe, otherwise falls back to a golden-section step. Convergence is declared when |x - midpoint| <= 2 * tol * |x| + EPS or when the iteration cap of 100 is reached.

Invariants:

• The bracket [a, b] always contains the minimum.
• f(x) <= f(w) <= f(v) at every iteration.

Kind: global function
Returns: Object - The argmin x and the value f(x).
Throws:

• Error When the bounds are equal or do not bracket bx.

Param	Type	Description
f	function	The function to minimize.
ax	number	Lower bound of the search interval.
bx	number	Initial guess within [ax, cx].
cx	number	Upper bound of the search interval.
tol	number	Convergence tolerance (default 1e-6).

runDifferentialEvolution(f, x0, opts) ⇒ Object

Differential-evolution minimization over an arbitrary-dimensional space.

The initial population is drawn uniformly inside [lb, ub]. Each member produces one trial per generation; the trial survives to the next generation only when its objective is strictly better.

Kind: global function
Returns: Object - Best point found and its objective value.

Param	Type	Description
f	function	Objective function.
x0	Array.<number>	Initial guess; used only to size the search space and the lower-bound default (x0[i] is ignored otherwise).
opts	Object	Algorithm options.
[opts.maxIter]	number	Maximum generations (default 500).
[opts.popSize]	number	Population size (default max(20, 10 * dim)).
[opts.cr]	number	Per-coordinate crossover probability (default 0.7).
[opts.f]	number	Differential scale factor F (default 0.8).
[opts.lb]	Array.<number>	Per-dimension lower bounds (default -5 for every dimension).
[opts.ub]	Array.<number>	Per-dimension upper bounds (default 5 for every dimension).

runNelderMead(f, x0, opts) ⇒ Object

Nelder-Mead minimization over a multidimensional space.

Builds an initial simplex by perturbing each axis of x0 by 1e-4 and then iterates the standard reflection / expansion / contraction / shrink move until either the spread of function values is below tol or maxIter iterations have been performed.

Kind: global function
Returns: Object - Best point, its function value, and the iteration count at termination.

Param	Type	Description
f	function	Objective function.
x0	Array.<number>	Initial guess (length determines dimension).
opts	Object	Algorithm options.
[opts.maxIter]	number	Maximum iterations (default 1000).
[opts.tol]	number	Convergence tolerance on the spread of f values across the simplex (default 1e-6).
[opts.alpha]	number	Reflection coefficient (default 1.0).
[opts.gamma]	number	Expansion coefficient (default 2.0).
[opts.rho]	number	Contraction coefficient (default 0.5).
[opts.sigma]	number	Shrink coefficient (default 0.5).

runSimulatedAnnealing(f, x0, opts) ⇒ Object

Simulated-annealing minimization over an arbitrary-dimensional space.

The neighbor for each iteration is generated by perturbing every coordinate by a uniform offset in [-stepSize, stepSize]. The acceptance temperature decays geometrically: temp *= coolingRate. The loop terminates once either maxIter iterations are performed or the temperature drops below finalTemp.

Kind: global function
Returns: Object - The best point found and its function value.

Param	Type	Description
f	function	Objective function.
x0	Array.<number>	Initial guess.
opts	Object	Algorithm options.
[opts.maxIter]	number	Maximum iterations (default 5000).
[opts.initialTemp]	number	Initial temperature (default 100).
[opts.finalTemp]	number	Temperature cut-off (default 0.001).
[opts.coolingRate]	number	Per-iteration multiplier (default 0.995).
[opts.stepSize]	number	Half-width of the uniform proposal distribution (default 0.1).

mulberry32(seed) ⇒ function

Constructs a mulberry32 generator with the given 32-bit seed.

The algorithm packs the state into a single unsigned 32-bit integer a. Each call applies two well-known integer mixing steps (Math.imul & bitwise shift) and returns the result divided by 2^32 so the output is in [0, 1).

Kind: global function
Returns: function - A function that returns the next uniform sample on every call.

Param	Type	Description
seed	number	PRNG seed (will be coerced to a 32-bit unsigned integer; >>> 0 performs the conversion).

setRandomSeed(seed)

Sets a global seed for reproducible simulations.

Passing null or undefined clears the seed and reverts to Math.random(). Calling setRandomSeed twice restarts the deterministic sequence from scratch.

Kind: global function

Param	Type	Description
seed	number | null | undefined	Integer seed (coerced to 32-bit). null/undefined clears the seed.

resetRandomSeed()

Resets the PRNG to use Math.random() for all subsequent draws.

Equivalent to setRandomSeed(null). Use this at the end of a deterministic experiment to restore nondeterministic behavior.

Kind: global function

nextRandom() ⇒ number

Returns a uniform random number in [0, 1).

Uses the seeded generator when one has been installed via setRandomSeed, otherwise falls through to Math.random(). Because this dispatcher is called from every stochastic primitive in the library, the entire computation tree is reproducible from a single seed.

Kind: global function
Returns: number - A pseudo-random number in [0, 1).

computeKsDistance(sample1, sample2, isSorted) ⇒ number

Computes the two-sample Kolmogorov-Smirnov distance.

Algorithm: a linear merged-pointer walk over the sorted order statistics. As we walk through the sorted union we maintain the empirical CDF values F_n(x) = (i + 1) / n and G_m(x) = j / m at the current position and record the absolute difference. Sorting first dominates the cost; the walk itself is O(n + m) where n = sample1.length and m = sample2.length.

Input validation:

• Both samples must be non-empty arrays or Float64Arrays.
• All values must be finite (no NaN, +Infinity, -Infinity).

Ties: when values are equal the walk advances both pointers and uses (i + 1) / n vs. (j + 1) / m for the distance — this matches the standard two-sided statistic.

Kind: global function
Returns: number - KS distance sup_x |F_n(x) - G_m(x)| in [0, 1].
Throws:

• Error When either input is not an array/typed array, is empty, or contains non-finite values.

Param	Type	Description
sample1	Array.<number> | Float64Array	First empirical sample.
sample2	Array.<number> | Float64Array	Second empirical sample.
isSorted	boolean	If true, skip sorting both samples. Off by default; setting this to true is the user's responsibility and is the hot path used inside rkSAVR's prepared-samples loop.

computeKsDistanceRescaled(sortedA, sortedB, factorA, factorB) ⇒ number

Kolmogorov-Smirnov distance for already sorted samples that need rescaling.

Equivalent to ksDistance(a, b, true) but applies the rescaling factors during the merged-pointer walk so no auxiliary allocation is needed. Multiplication by a positive scalar is order-preserving, so the pre-sorting of the inputs is unaffected by the choice of factorA and factorB.

This is the hot path of the Hurstify estimator's inner loop: O(n + m) per evaluation, no allocations beyond the locals below.

Kind: global function
Returns: number - KS distance between the rescaled samples in [0, 1].
Throws:

• Error When either input is not an array/typed array or is empty.

Param	Type	Description
sortedA	Array.<number> | Float64Array	Pre-sorted sample A.
sortedB	Array.<number> | Float64Array	Pre-sorted sample B.
factorA	number	Positive rescaling factor for A (typically a^{-H}).
factorB	number	Positive rescaling factor for B.

shuffleArray(array) ⇒ Array.<*>

Unbiased Fisher-Yates shuffle.

Returns a new array; the input is never mutated. Uses the seeded PRNG exposed by prng.js, so the result is reproducible when a seed is set.

Complexity: O(n) time, O(n) extra memory.

Kind: global function
Returns: Array.<*> - Shuffled copy of array.

Param	Type	Description
array	Array.<*>	Input array (not modified).

permuteBlocks(data, blockSize, randomPhase) ⇒ Array.<*>

Block random permutation for decorrelating serial dependence.

Conceptually this is the paper's "preserves marginals, kills short-range autocorrelation" operation:

1. (Optional) shift the starting index by a uniform [-0, blockSize) offset so two calls with the same seed still produce different alignments.
2. Slice the resulting series into blocks of length blockSize (the first block may be shorter than blockSize when a phase offset was applied).
3. Apply a Fisher-Yates shuffle to the block list.
4. Concatenate the shuffled blocks back into a single sequence.

Picking blockSize is the user's responsibility: it should be larger than the dominant autocorrelation length in data. Too small and serial dependence survives; too large and the number of blocks — and therefore the effective randomization — shrinks.

Kind: global function
Returns: Array.<*> - Permuted array containing exactly the same elements as data.
Throws:

• Error When data is not array-like or blockSize is out of range.

Param	Type	Description
data	Array.<*>	Input array (not modified).
blockSize	number	Block length; must satisfy 0 < blockSize <= data.length.
randomPhase	boolean	Whether to apply a random starting phase offset.

getRandomSample(array, n) ⇒ Array.<*>

Floyd's Algorithm R reservoir sampler.

Streams over the input producing a uniformly random sample of size n without replacement. Equivalent to shuffle(array).slice(0, n) but uses only O(n) auxiliary memory and a single pass through array, which matters when sampling from very large arrays (e.g. millions of increments).

Edge cases:

• n >= array.length: returns a shuffled full copy of array.
• n <= 0: returns an empty array.

Kind: global function
Returns: Array.<*> - Random sample of size min(n, array.length).

Param	Type	Description
array	Array.<*>	Input array.
n	number	Number of elements to sample.

nextGaussian() ⇒ number

Draws a single standard normal via Box-Muller.

The polar variant is implemented by guarding against degenerate u === 0 draws from nextRandom(). One Box-Muller pair yields two independent standard normals; this routine keeps the cosine component and discards the sine. Use generateCorrelatedGaussian if you need both halves, or call nextGaussian twice with distinct nextRandom() outputs.

Kind: global function
Returns: number - A standard normal random variable.

generateGaussianBatch(n) ⇒ Float64Array

Pre-allocates a Float64Array of standard normals.

Useful when an inner loop needs a contiguous buffer of normals; the allocation is amortized across a single batch draw, whereas repeated nextGaussian calls would each allocate internally.

Kind: global function
Returns: Float64Array - Buffer of n independent standard normals.

Param	Type	Description
n	number	Number of samples (n >= 0).

generateCorrelatedGaussian(n, rho) ⇒ Array.<Float64Array>

Generates two correlated standard-normal streams via Cholesky.

Mathematically the model is (Z1, Z2) with unit marginals and Corr(Z1, Z2) = rho. Implementation: draw an i.i.d. Box-Muller pair (z1, z2); set Z1 = z1; set Z2 = rho * z1 + sqrt(1 - rho^2) * z2. Both Z1 and Z2 have unit variance and exactly correlation rho.

Important: rho must be strictly in (-1, 1); the implementation silently clamps 1 - rho^2 to zero via Math.max(0, ...) so the endpoints collapse to the trivial deterministic case.

Kind: global function
Returns: Array.<Float64Array> - [Z1, Z2] of length n.

Param	Type	Description
n	number	Number of samples.
rho	number	Target correlation in (-1, 1).

generateFractionalNoise(n, H) ⇒ Float64Array

Fractional Gaussian Noise via Hosking's method.

Hosking's method is an exact O(n^2) Cholesky-style recursion that generates samples from the autocovariance gamma(k) = 0.5 (|k+1|^{2H} - 2|k|^{2H} + |k-1|^{2H}).

It uses O(n) recursion updates to compute the conditional mean and variance (phi, v) incrementally, so the per-step cost is O(k) and the total O(n^2). This is fine for the scales used in the paper (a few hundred to a few thousand samples) but dominates for n >> 1e4.

Assumptions:

• n > 0 and H in (0, 1).
• The result is mean-zero (the recursion conditions on x_0 ~ N(0, 1)).

Kind: global function
Returns: Float64Array - A contiguous fGN sample of length n.
Throws:

• Error When n is not a positive finite integer or H is out of range.

Param	Type	Description
n	number	Length of the desired sample.
H	number	Hurst parameter; must satisfy 0 < H < 1.

generateFractionalBrownianMotion(n, H) ⇒ Float64Array

Fractional Brownian Motion by cumulative summation of fGN.

The implementation delegates the heavy lifting to generateFractionalNoise and then performs a single O(n) cumulative-sum pass. The first sample is fixed at 0 (the standard convention for fBM(0) = 0), so paths always start at the origin.

For non-zero means, simply add a constant afterwards — fGn is mean-zero by construction.

Kind: global function
Returns: Float64Array - fBm path of length n (Float64Array(0) when n <= 0).
Throws:

• Error When H is out of range (propagated from generateFractionalNoise).

Param	Type	Description
n	number	Length of the path.
H	number	Hurst parameter; must satisfy 0 < H < 1.

computeFractionalKernel(H, nSteps, dt) ⇒ Float64Array

Precomputes the Riemann-Liouville fractional kernel used by the rough-volatility simulators.

Mathematically K(t) = sqrt(2 H) * t^{H - 0.5} for t > 0. The result is a length-nSteps array where entry i corresponds to t = (i + 1) * dt.

Reusing a precomputed kernel for every path avoids the O(n^2) cost of re-evaluating the power function per integration step.

Kind: global function
Returns: Float64Array - Kernel values of length nSteps.

Param	Type	Description
H	number	Hurst parameter.
nSteps	number	Number of time steps covered by the kernel.
dt	number	Per-step time increment.

computeFractionalIntegral(dW, kernel, t) ⇒ number

Computes a single time-step of the Riemann-Liouville fractional integral.

Given precomputed Brownian increments dW and a kernel from computeFractionalKernel, returns I_t = sum_{j=0}^{t-1} K(t - j) * dW_j.

Used inside the rBergomi path generator and the exact fOU driver.

Complexity: O(t) per call, so building a full path is O(n^2). This is acceptable for paths up to a few hundred steps; for long simulations switch to a circulant-embedding FFT approximation (not implemented here).

Kind: global function
Returns: number - Fractional integral value at time t.

Param	Type	Description
dW	Float64Array	Brownian increments.
kernel	Float64Array	Precomputed kernel of length >= t.
t	number	Current time index (exclusive upper bound).

xavierInit(rows, cols) ⇒ Array.<Array.<number>>

Xavier (Glorot-uniform) weight initialization.

Produces a rows x cols matrix where each entry is sampled uniformly in [-scale, scale] with scale = sqrt(2 / (rows + cols)). This is the standard initializer for tanh/sigmoid-activated layers (Glorot & Bengio, 2010).

Kind: global function
Returns: Array.<Array.<number>> - Initialized weight matrix.

Param	Type	Description
rows	number	Number of rows.
cols	number	Number of columns.

getBinomialCoeffs(d, lag) ⇒ Float64Array

Returns the binomial coefficient sequence [C(d, 0), ..., C(d, lag)].

Uses a tiny FIFO cache keyed by ${d}:${lag} so that identical lookups within a rolling ARFIMA run are O(1). When the cache is full the oldest entry is evicted.

Kind: global function
Returns: Float64Array - Coefficient vector of length lag + 1.

Param	Type	Description
d	number	Differencing parameter.
lag	number	Maximum lag (inclusive).

fractionalDifference(data, d, [lag]) ⇒ Array.<number>

Computes the (truncated) fractional difference of a series for a given d and lag cap. The truncation to lag keeps the per-step cost O(lag) rather than O(t), which is essential for long-history forecasting.

Kind: global function
Returns: Array.<number> - Fractionally differenced series.

Param	Type	Default	Description
data	Array.<number>		Input series.
d	number		Differencing parameter.
[lag]	number	50	Maximum lag for the binomial expansion (default 50).

ksCriticalValue(n, m, alpha) ⇒ number

Two-sample Kolmogorov–Smirnov asymptotic critical value.

D_alpha = sqrt(-0.5 * ln(alpha / 2)) * sqrt((n + m) / (n * m))

Kind: global function
Returns: number - Critical value D_alpha.

Param	Type	Description
n	number	First sample size.
m	number	Second sample size.
alpha	number	Significance level (default 0.05).

ksPvalue(D, n, m) ⇒ number

Approximate two-sample KS p-value via the asymptotic Kolmogorov distribution.

Q(lambda) ~ 2 * sum_{j=1..3} (-1)^{j-1} * exp(-2 j^2 lambda^2)

with the standard lambda correction.

Kind: global function
Returns: number - Approximate p-value in [0, 1].

Param	Type	Description
D	number	Observed KS distance.
n	number	First sample size.
m	number	Second sample size.

kalmanLogLikelihood(observations, q, r) ⇒ number

Log-likelihood of the observations under a 1D Kalman filter.

Kind: global function
Returns: number - Total log-likelihood (or -Infinity for empty input).

Param	Type	Description
observations	Array.<number>	Time-ordered H estimates.
q	number	Process-noise variance.
r	number	Measurement-noise variance.

detectCusumBreakpoints(hHistory, windowSize, threshold) ⇒ Array.<{index: number, H_before: number, H_after: number}>

Detects breakpoints in a series of H estimates via a sliding-window CUSUM.

Kind: global function
Returns: Array.<{index: number, H_before: number, H_after: number}> - Detected breakpoints in chronological order.

Param	Type	Description
hHistory	Array.<number>	Time-ordered series of H estimates.
windowSize	number	Sliding window size (default 50).
threshold	number	CUSUM threshold (default 3.0).

chooseKsObjective([scales], [weights]) ⇒ KsObjective

Selects the right KsObjective for a configuration.

Kind: global function
Returns: KsObjective - The matching strategy.

Param	Type	Description
[scales]	Array.<number>	Optional scale array.
[weights]	Array.<number>	Optional weights.

defaultSampler([blockSize]) ⇒ Sampler

Convenience: selects the default sampler based on blockSize.

• When blockSize is a positive number a BlockPermutationSampler is returned.
• Otherwise a ReservoirSampler is returned.

Kind: global function
Returns: Sampler - Either a BlockPermutationSampler or a ReservoirSampler.

Param	Type	Description
[blockSize]	number	Block length for the permutation sampler; omit (or pass 0/negative) to get the reservoir sampler instead.

KsSignificanceResult : Object

Kind: global typedef

ConstancyResult : Object

Kind: global typedef

CusumBreakResult : Object

Kind: global typedef

BootstrapCiResult : Object

Kind: global typedef

SimulationResult : Object

Kind: global typedef

PriceResult : Object

Kind: global typedef