use greeners::VARMA;
use ndarray::Array2;
use ndarray_rand::rand_distr::Normal;
use rand::prelude::*;
fn main() -> Result<(), Box<dyn std::error::Error>> {
let t_obs = 2000;
let mut rng = rand::thread_rng();
let norm = Normal::new(0.0, 1.0).unwrap();
let mut data = Array2::<f64>::zeros((t_obs, 2));
let mut errors = Array2::<f64>::zeros((t_obs, 2));
for t in 1..t_obs {
let u1 = norm.sample(&mut rng);
let u2 = norm.sample(&mut rng);
errors[[t, 0]] = u1;
errors[[t, 1]] = u2;
let y1_prev = data[[t - 1, 0]];
let y2_prev = data[[t - 1, 1]];
let u1_prev = errors[[t - 1, 0]];
let u2_prev = errors[[t - 1, 1]];
data[[t, 0]] = (0.6 * y1_prev) + u1 + (0.4 * u2_prev);
data[[t, 1]] = (0.6 * y2_prev) + u2 + (0.4 * u1_prev);
}
println!("--- VARMA(1, 1) Simulation ---");
println!("True AR (A) Diagonal: 0.6");
println!("True MA (M) Off-Diagonal: 0.4\n");
let model = VARMA::fit(&data, 1, 1)?;
println!("{}", model);
println!("Estimated AR Matrix (Should be close to 0.6 on diag):");
println!(
"Eq 1: {:.4} * Y1(-1) + {:.4} * Y2(-1)",
model.ar_params[[1, 0]],
model.ar_params[[1, 1]]
);
println!(
"Eq 2: {:.4} * Y1(-1) + {:.4} * Y2(-1)",
model.ar_params[[1, 0]],
model.ar_params[[1, 1]]
);
println!("\nCheck Coefficients:");
println!("Eq 1 (Y1):");
println!(" AR(Y1_L1): {:.4} (Target 0.6)", model.ar_params[[1, 0]]);
println!(" MA(U2_L1): {:.4} (Target 0.4)", model.ma_params[[1, 0]]);
println!("Eq 2 (Y2):");
println!(" AR(Y2_L1): {:.4} (Target 0.6)", model.ar_params[[2, 1]]);
println!(" MA(U1_L1): {:.4} (Target 0.4)", model.ma_params[[0, 1]]);
Ok(())
}