diff --git a/src/primary/correlation.rs b/src/primary/correlation.rs
new file mode 100644
index 0000000..004e7d9
--- /dev/null
+++ b/src/primary/correlation.rs
@@ -0,0 +1,38 @@
+/// Returns the correlation of a sequence of `f64` values.
+///
+/// Empty input returns `None`.
+///
+/// # Example
+/// ```
+/// use rustats::primary::correlation;
+///
+/// let values_x = [1.0, 2.0, 3.0];
+/// let values_y = [4.0, 5.0, 6.0];
+/// let r = correlation(&values_x, &values_y);
+/// assert_eq!(r, Some(1.0));
+/// ```
+use super::{covariance, std_deviation};
+
+pub fn correlation(x: &[f64], y: &[f64]) -> Option<f64> {
+
+    if x.len() != y.len(){
+        return None;
+    }
+
+    if x.len() < 2 {
+            return None;
+    }
+
+    let cov = covariance(x, y)?;
+    let std_x = std_deviation(x)?;
+    let std_y = std_deviation(y)?;
+
+    if std_x == 0.0 || std_y == 0.0 {
+        return None;
+    }
+
+    Some(cov / (std_x*std_y))
+
+}
+    
+
diff --git a/src/primary/mod.rs b/src/primary/mod.rs
index 2ea7492..9f12e79 100644
--- a/src/primary/mod.rs
+++ b/src/primary/mod.rs
@@ -7,6 +7,7 @@ pub mod std_deviation;
 pub mod interquartile_range;
 pub mod range;
 pub mod covariance;
+pub mod correlation;
 
 
 pub use mean::mean;
@@ -18,3 +19,4 @@ pub use std_deviation::std_deviation;
 pub use covariance::covariance;
 pub use range::range;
 pub use interquartile_range::interquartile_range;
+pub use correlation::correlation;
\ No newline at end of file
diff --git a/tests/correlation.rs b/tests/correlation.rs
new file mode 100644
index 0000000..edc2392
--- /dev/null
+++ b/tests/correlation.rs
@@ -0,0 +1,114 @@
+use rustats::primary::correlation;
+
+fn approx_eq(left: f64, right: f64, epsilon: f64) -> bool {
+    (left - right).abs() <= epsilon
+}
+
+#[test]
+fn correlation_empty_returns_none() {
+    // Correlation requires at least 2 values in each input (denominator is n-1)
+    assert_eq!(correlation(&[], &[]), None);
+}
+
+#[test]
+fn correlation_different_lengths(){
+    //Correlational requires both input arrays to have the same length
+    let values_x = [1.0, 2.0, 3.0];
+    let values_y = [2.0, 3.0];
+    assert_eq!(correlation(&values_x, &values_y),None);
+}
+
+#[test]
+fn correlation_single_element_returns_none() {
+// Correlation requires at least two observations.
+assert_eq!(correlation(&[5.0], &[9.0]), None);
+}
+
+#[test]
+fn correlation_zero_std_x_returns_none() {
+    // Correlation requires the std deviation of x to be non zero
+    assert_eq!(correlation(&[4.0, 4.0, 4.0], &[1.0, 2.0, 3.0]), None);
+}
+
+#[test]
+fn correlation_zero_std_y_returns_none() {
+    // Correlation requires the std deviation of y to be non zero
+    assert_eq!(correlation(&[1.0, 2.0 , 3.0], &[4.0, 4.0, 4.0]), None);
+}
+
+#[test]
+fn correlation_two_identical_sample_arrays() {
+    // Identical arrays have perfect positive linear correlation
+    let values_x = [1.0, 2.0, 3.0, 4.0];
+    let values_y = [1.0, 2.0, 3.0, 4.0];
+    let r = correlation(&values_x, &values_y).expect("two or more values should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+}
+
+#[test]
+fn correlation_smallest_valid_input() {
+    // covariance = 0.5, std_x = sqrt(0.5), std_y = sqrt(0.5), correlation = 0.5 / (sqrt(0.5)*sqrt(0.5)) 
+    let values_x = [1.0, 2.0];
+    let values_y = [3.0, 4.0];
+    let r = correlation(&values_x, &values_y).expect("two or more values should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+}
+
+#[test]
+fn correlation_multiple_values() {
+    //correlation should be exactly 1.0
+    let values_x = [1.0, 2.0, 3.0, 4.0, 5.0];
+    let values_y = [5.0, 6.0, 7.0, 8.0, 9.0];
+    let r = correlation(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+}
+
+#[test]
+fn correlation_with_negative_values() {
+    // Both variables increase together, so the correlation is 1.0
+    let values_x = [-2.0, -1.0, 1.0, 2.0];
+    let values_y = [-4.0, -2.0, 2.0, 4.0];
+    let r = correlation(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+}
+
+#[test]
+fn correlation_with_decimals() {
+    // Decimal values preserve a perfect positive linear relationship
+    let values_x = [0.5, 1.5, 2.5, 3.5];
+    let values_y = [0.5, 1.5, 2.5, 3.5];
+    let r = correlation(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+}
+
+#[test]
+fn correlation_negative() {
+    // When one variable increases while the other decreases linearly, the correlation is 1.0
+    let values_x = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let values_y = [10.0, 0.0, -10.0, -5.0, 5.0];
+    let r = correlation(&values_x, &values_y).expect("non-empty input should return Some");
+    dbg!(r);
+    assert!(approx_eq(r, -1.0, 1e-12));
+}
+
+#[test]
+fn correlation_positive() {
+    // 
+    let values_x = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let values_y = [-10.0, 0.0, 10.0, 5.0, -5.0];
+    let r = correlation(&values_x, &values_y).expect("non-empty input should return Some");
+    assert!(approx_eq(r, 1.0, 1e-12));
+
+}
+
+#[test]
+fn correlation_zero() {
+// These inputs have zero linear relationship.
+let values_x = [-1.0, 1.0, -1.0, 1.0];
+let values_y = [-1.0, -1.0, 1.0, 1.0];
+
+let r = correlation(&values_x, &values_y)
+.expect("valid input should return Some");
+
+assert!(approx_eq(r, 0.0, 1e-12));
+}
\ No newline at end of file