feat: add auto regressive processes

src/stochastic.rs

@@ -15,6 +15,7 @@
 //! | **volatility**   | Focuses on modeling stochastic volatility, including processes like the Heston model, which are used to simulate changes in volatility over time in financial markets.                                                    |
 //!
 
+pub mod autoregressive;
 pub mod diffusion;
 pub mod interest;
 pub mod isonormal;

src/stochastic/autoregressive.rs

@@ -0,0 +1,9 @@
+pub mod agrach;
+pub mod ar;
+pub mod arch;
+pub mod arima;
+pub mod egarch;
+pub mod garch;
+pub mod ma;
+pub mod sarima;
+pub mod tgarch;

src/stochastic/autoregressive/agrach.rs

@@ -0,0 +1,115 @@
+use impl_new_derive::ImplNew;
+use ndarray::Array1;
+use ndarray_rand::RandomExt;
+use rand_distr::Normal;
+
+use crate::stochastic::Sampling;
+
+/// A generic Asymmetric GARCH(p,q) model (A-GARCH),
+/// allowing a separate "delta" term for negative-lag effects.
+///
+/// One possible form:
+/// \[
+///   \sigma_t^2
+///     = \omega
+///       + \sum_{i=1}^p \Bigl[\alpha_i X_{t-i}^2
+///                             + \delta_i X_{t-i}^2 \mathbf{1}_{\{X_{t-i}<0\}}\Bigr]
+///       + \sum_{j=1}^q \beta_j \sigma_{t-j}^2,
+///   \quad X_t = \sigma_t \cdot z_t, \quad z_t \sim \mathcal{N}(0,1).
+/// \]
+///
+/// # Parameters
+/// - `omega`: Constant term \(\omega\).
+/// - `alpha`: Array \(\{\alpha_1, \ldots, \alpha_p\}\) for positive squared terms.
+/// - `delta`: Array \(\{\delta_1, \ldots, \delta_p\}\) for negative-lag extra effect.
+/// - `beta`:  Array \(\{\beta_1, \ldots, \beta_q\}\).
+/// - `n`:     Length of the time series.
+/// - `m`:     Optional batch size (unused by default).
+///
+/// # Notes
+/// - This is essentially a T-GARCH-like structure but with different naming (`delta`).
+/// - Stationarity constraints typically require \(\sum \alpha_i + \tfrac{1}{2}\sum \delta_i + \sum \beta_j < 1\).
+#[derive(ImplNew)]
+pub struct AGARCH {
+  pub omega: f64,
+  pub alpha: Array1<f64>,
+  pub delta: Array1<f64>,
+  pub beta: Array1<f64>,
+  pub n: usize,
+  pub m: Option<usize>,
+}
+
+impl Sampling<f64> for AGARCH {
+  fn sample(&self) -> Array1<f64> {
+    let p = self.alpha.len();
+    let q = self.beta.len();
+
+    // Generate white noise
+    let z = Array1::random(self.n, Normal::new(0.0, 1.0).unwrap());
+
+    // Arrays for X_t and sigma_t^2
+    let mut x = Array1::<f64>::zeros(self.n);
+    let mut sigma2 = Array1::<f64>::zeros(self.n);
+
+    // Summation for unconditional variance init
+    let sum_alpha: f64 = self.alpha.iter().sum();
+    let sum_delta_half: f64 = self.delta.iter().sum::<f64>() * 0.5;
+    let sum_beta: f64 = self.beta.iter().sum();
+    let denom = (1.0 - sum_alpha - sum_delta_half - sum_beta).max(1e-8);
+
+    for t in 0..self.n {
+      if t == 0 {
+        sigma2[t] = self.omega / denom;
+      } else {
+        let mut var_t = self.omega;
+        // p-lag terms
+        for i in 1..=p {
+          if t >= i {
+            let x_lag = x[t - i];
+            let indicator = if x_lag < 0.0 { 1.0 } else { 0.0 };
+
+            var_t +=
+              self.alpha[i - 1] * x_lag.powi(2) + self.delta[i - 1] * x_lag.powi(2) * indicator;
+          }
+        }
+        // q-lag terms
+        for j in 1..=q {
+          if t >= j {
+            var_t += self.beta[j - 1] * sigma2[t - j];
+          }
+        }
+        sigma2[t] = var_t;
+      }
+      x[t] = sigma2[t].sqrt() * z[t];
+    }
+
+    x
+  }
+
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use ndarray::arr1;
+
+  use crate::{
+    plot_1d,
+    stochastic::{autoregressive::agrach::AGARCH, Sampling},
+  };
+
+  #[test]
+  fn agarch_plot() {
+    let alpha = arr1(&[0.05, 0.01]); // p=2
+    let delta = arr1(&[0.03, 0.01]); // p=2
+    let beta = arr1(&[0.8]); // q=1
+    let agarchpq = AGARCH::new(0.1, alpha, delta, beta, 100, None);
+    plot_1d!(agarchpq.sample(), "A-GARCH(p,q) process");
+  }
+}

src/stochastic/autoregressive/ar.rs

@@ -0,0 +1,91 @@
+use impl_new_derive::ImplNew;
+use ndarray::Array1;
+use ndarray_rand::RandomExt;
+use rand_distr::Normal;
+
+use crate::stochastic::Sampling;
+
+/// Implements an AR(p) model:
+///
+/// \[
+///   X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \dots + \phi_p X_{t-p}
+///         + \epsilon_t,
+///   \quad \epsilon_t \sim \mathcal{N}(0, \sigma^2).
+/// \]
+///
+/// # Fields
+/// - `phi`: Vector of AR coefficients (\(\phi_1, \ldots, \phi_p\)).
+/// - `sigma`: Standard deviation of the noise \(\epsilon_t\).
+/// - `n`: Length of the time series.
+/// - `m`: Optional batch size (for parallel sampling).
+/// - `x0`: Optional array of initial values. If provided, should have length at least `phi.len()`.
+#[derive(ImplNew)]
+pub struct ARp {
+  /// AR coefficients
+  pub phi: Array1<f64>,
+  /// Noise std dev
+  pub sigma: f64,
+  /// Number of observations
+  pub n: usize,
+  /// Optional batch size
+  pub m: Option<usize>,
+  /// Optional initial conditions
+  pub x0: Option<Array1<f64>>,
+}
+
+impl Sampling<f64> for ARp {
+  fn sample(&self) -> Array1<f64> {
+    let p = self.phi.len();
+    let noise = Array1::random(self.n, Normal::new(0.0, self.sigma).unwrap());
+    let mut series = Array1::<f64>::zeros(self.n);
+
+    // Fill initial conditions if provided
+    if let Some(init) = &self.x0 {
+      // Copy up to min(p, n)
+      for i in 0..p.min(self.n) {
+        series[i] = init[i];
+      }
+    }
+
+    // AR recursion
+    for t in 0..self.n {
+      let mut val = 0.0;
+      // Sum over AR lags
+      for k in 1..=p {
+        if t >= k {
+          val += self.phi[k - 1] * series[t - k];
+        }
+      }
+      // Add noise
+      series[t] += val + noise[t];
+    }
+
+    series
+  }
+
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use ndarray::arr1;
+
+  use crate::{
+    plot_1d,
+    stochastic::{autoregressive::ar::ARp, Sampling},
+  };
+
+  #[test]
+  fn ar_plot() {
+    // Suppose p=2 with user-defined coefficients
+    let phi = arr1(&[0.5, -0.25]);
+    let ar_model = ARp::new(phi, 1.0, 100, None, None);
+    plot_1d!(ar_model.sample(), "AR(p) process");
+  }
+}

src/stochastic/autoregressive/arch.rs

@@ -0,0 +1,78 @@
+use impl_new_derive::ImplNew;
+use ndarray::Array1;
+use ndarray_rand::RandomExt;
+use rand_distr::Normal;
+
+use crate::stochastic::Sampling;
+
+/// Implements an ARCH(m) model:
+///
+/// \[
+///   \sigma_t^2 = \omega + \sum_{i=1}^m \alpha_i X_{t-i}^2,
+///   \quad X_t = \sigma_t \cdot z_t, \quad z_t \sim \mathcal{N}(0,1).
+/// \]
+///
+/// # Fields
+/// - `omega`: Constant term.
+/// - `alpha`: Array of ARCH coefficients.
+/// - `n`: Number of observations.
+/// - `m`: Optional batch size.
+#[derive(ImplNew)]
+pub struct ARCH {
+  /// Omega (constant term in variance)
+  pub omega: f64,
+  /// Coefficients alpha_i
+  pub alpha: Array1<f64>,
+  /// Length of series
+  pub n: usize,
+  /// Optional batch size
+  pub m: Option<usize>,
+}
+
+impl Sampling<f64> for ARCH {
+  fn sample(&self) -> Array1<f64> {
+    let m = self.alpha.len();
+    let z = Array1::random(self.n, Normal::new(0.0, 1.0).unwrap());
+    let mut x = Array1::<f64>::zeros(self.n);
+
+    for t in 0..self.n {
+      // compute sigma_t^2
+      let mut var_t = self.omega;
+      for i in 1..=m {
+        if t >= i {
+          let x_lag = x[t - i];
+          var_t += self.alpha[i - 1] * x_lag.powi(2);
+        }
+      }
+      let sigma_t = var_t.sqrt();
+      x[t] = sigma_t * z[t];
+    }
+
+    x
+  }
+
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use ndarray::arr1;
+
+  use crate::{
+    plot_1d,
+    stochastic::{autoregressive::arch::ARCH, Sampling},
+  };
+
+  #[test]
+  fn arch_plot() {
+    let alpha = arr1(&[0.2, 0.1]);
+    let arch_model = ARCH::new(0.1, alpha, 100, None);
+    plot_1d!(arch_model.sample(), "ARCH(m) process");
+  }
+}