diff --git a/NAMESPACE b/NAMESPACE
index dc312044..3e783711 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -25,6 +25,7 @@ S3method(pcor_to_cor,matrix)
 S3method(plot,easycor_test)
 S3method(plot,easycormatrix)
 S3method(plot,easycorrelation)
+S3method(print,cor_diff)
 S3method(print,easycormatrix)
 S3method(print,easycorrelation)
 S3method(print,easymatrixlist)
@@ -38,6 +39,7 @@ S3method(summary,easycorrelation)
 S3method(visualisation_recipe,easycor_test)
 S3method(visualisation_recipe,easycormatrix)
 S3method(visualisation_recipe,easycorrelation)
+export(cor_diff)
 export(cor_lower)
 export(cor_smooth)
 export(cor_sort)
diff --git a/R/cor_diff.R b/R/cor_diff.R
new file mode 100644
index 00000000..108f8385
--- /dev/null
+++ b/R/cor_diff.R
@@ -0,0 +1,121 @@
+#' Test differences between correlations
+#'
+#' @description
+#' Tests whether the correlation between two variables `x` and `y` is different
+#' from the correlation between `x2` and `y2`.
+#'
+#' `cor_diff()` returns a table containing an index of difference precision (i.e.,
+#' the estimated difference divided by its standard error) and an associated p-value.
+#' A significant p-value indicates that the correlation between `x` and `y` is
+#' different from the correlation between `x2` and `y2`.
+#'
+#' @param data A data frame of observations.
+#' @param x,y,x2,y2 The variable names in `data` to be used. `x` and `y` can also
+#' be pairs of variables, in which case the second variable is used as `x2` and `y2`.
+#' @param method Can be `"parametric"` or `"bootstrapping"`. If `"parametric"`,
+#' the [psych::r.test()] function is used. If `"bootstrapping"`, a bootstrapping
+#' procedure is used.
+#' @param ... Other arguments to be passed, for instance `iterations` (default: 1000)
+#' if method is bootstrapping.
+#'
+#' @examples
+#' cor_diff(iris, c("Sepal.Length", "Sepal.Width"), c("Sepal.Length", "Petal.Width"))
+#' cor_diff(iris,
+#'   c("Sepal.Length", "Sepal.Width"),
+#'   c("Sepal.Length", "Petal.Width"),
+#'   method = "bootstrapping", iterations = 100
+#' )
+#' @export
+cor_diff <- function(data, x, y, x2 = NULL, y2 = NULL, method = "parametric", ...) {
+  # If pairs are passed
+  if (length(x) == 2 && length(y) == 2) {
+    x2 <- y[1]
+    y2 <- y[2]
+    y <- x[2]
+    x <- x[1]
+  }
+
+  # Compute
+  if (method %in% c("bootstrapping")) {
+    out <- .cor_diff_bootstrapping(data, x, y, x2, y2, ...)
+  } else {
+    out <- .cor_diff_parametric(data, x, y, x2, y2, ...)
+  }
+  class(out) <- c("cor_diff", class(out))
+  out
+}
+
+
+
+# Methods -----------------------------------------------------------------
+
+
+
+#' @keywords internal
+.cor_diff_parametric <- function(data, x, y, x2, y2, ...) {
+  insight::check_if_installed("psych", "for 'parametric' correlation difference method")
+
+  args <- list(n = nrow(data), r12 = stats::cor(data[[x]], data[[y]]))
+  if (x == x2 && y != y2) {
+    args$r13 <- stats::cor(data[[x]], data[[y2]])
+    args$r23 <- stats::cor(data[[y]], data[[y2]])
+  } else if (y == y2 && x != x2) {
+    args$r13 <- stats::cor(data[[y]], data[[x2]])
+    args$r23 <- stats::cor(data[[x]], data[[x2]])
+  } else {
+    args$r34 <- stats::cor(data[[x2]], data[[y2]])
+  }
+  test <- do.call(psych::r.test, args)
+
+  out <- data.frame(
+    Method = "parametric"
+  )
+  if ("t" %in% names(test)) {
+    out$t <- test$t
+  } else {
+    out$z <- test$z
+  }
+  out$p <- test$p
+  out
+}
+
+#' @keywords internal
+.cor_diff_bootstrapping <- function(data, x, y, x2, y2, iterations = 1000, robust = FALSE, ...) {
+  diff <- rep(NA, iterations) # Initialize vector
+
+  # Bootstrap
+  for (i in 1:iterations) {
+    # Take random sample of data
+    dat <- data[sample(nrow(data), nrow(data), replace = TRUE), ]
+    # Compute diff
+    diff[i] <- stats::cor(dat[[x]], dat[[y]]) - stats::cor(dat[[x2]], dat[[y2]])
+  }
+
+  # Summarize
+  if (robust == FALSE) {
+    out <- data.frame(
+      Method = "bootstrapping",
+      z = mean(diff) / stats::sd(diff),
+      p = bayestestR::pd_to_p(as.numeric(bayestestR::p_direction(diff)))
+    )
+  } else {
+    out <- data.frame(
+      Method = "bootstrapping_robust",
+      z = stats::median(diff) / stats::mad(diff),
+      p = bayestestR::pd_to_p(as.numeric(bayestestR::p_direction(diff)))
+    )
+  }
+  out
+}
+
+
+
+# Printing ----------------------------------------------------------------
+
+#' @export
+print.cor_diff <- function(x, ...) {
+  insight::format_table(x, ...) |>
+    insight::export_table(title = "Correlation Difference Test") |>
+    print()
+  invisible(x)
+}
diff --git a/man/cor_diff.Rd b/man/cor_diff.Rd
new file mode 100644
index 00000000..71afb1fb
--- /dev/null
+++ b/man/cor_diff.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cor_diff.R
+\name{cor_diff}
+\alias{cor_diff}
+\title{Test differences between correlations}
+\usage{
+cor_diff(data, x, y, x2 = NULL, y2 = NULL, method = "parametric", ...)
+}
+\arguments{
+\item{data}{A data frame of observations.}
+
+\item{x, y, x2, y2}{The variable names in \code{data} to be used. \code{x} and \code{y} can also
+be pairs of variables, in which case the second variable is used as \code{x2} and \code{y2}.}
+
+\item{method}{Can be \code{"parametric"} or \code{"bootstrapping"}. If \code{"parametric"},
+the \code{\link[psych:r.test]{psych::r.test()}} function is used. If \code{"bootstrapping"}, a bootstrapping
+procedure is used.}
+
+\item{...}{Other arguments to be passed, for instance \code{iterations} (default: 1000)
+if method is bootstrapping.}
+}
+\description{
+Tests whether the correlation between two variables \code{x} and \code{y} is different
+from the correlation between \code{x2} and \code{y2}.
+
+\code{cor_diff()} returns a table containing an index of difference precision (i.e.,
+the estimated difference divided by its standard error) and an associated p-value.
+A significant p-value indicates that the correlation between \code{x} and \code{y} is
+different from the correlation between \code{x2} and \code{y2}.
+}
+\examples{
+cor_diff(iris, c("Sepal.Length", "Sepal.Width"), c("Sepal.Length", "Petal.Width"))
+cor_diff(iris,
+  c("Sepal.Length", "Sepal.Width"),
+  c("Sepal.Length", "Petal.Width"),
+  method = "bootstrapping", iterations = 100
+)
+}
diff --git a/tests/testthat/test-cor_diff.R b/tests/testthat/test-cor_diff.R
new file mode 100644
index 00000000..fea85889
--- /dev/null
+++ b/tests/testthat/test-cor_diff.R
@@ -0,0 +1,7 @@
+test_that("cor_diff", {
+  expect_equal(
+    cor_diff(iris, "Sepal.Length", "Sepal.Width", "Sepal.Length", "Petal.Width")$t,
+    -10,
+    tolerance = 0.001
+  )
+})