diff --git a/NAMESPACE b/NAMESPACE
index 42ebf32aa133f32d15189d4a9f26abe9df0ba500..7cbec955260f4cd54997046277a3d9c886e599db 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -55,12 +55,14 @@ export(REOF)
 export(RMS)
 export(RMSSS)
 export(RandomWalkTest)
+export(RatioPredictableComponents)
 export(RatioRMS)
 export(RatioSDRMS)
 export(Regression)
 export(Reorder)
 export(SPOD)
 export(Season)
+export(SignalNoiseRatio)
 export(Smoothing)
 export(Spectrum)
 export(Spread)
diff --git a/R/RatioPredictableComponents.R b/R/RatioPredictableComponents.R
new file mode 100644
index 0000000000000000000000000000000000000000..270420ab272bc392caff74494c46f244844e2db2
--- /dev/null
+++ b/R/RatioPredictableComponents.R
@@ -0,0 +1,108 @@
+#'Calculate ratio of predictable components (RPC)
+#'
+#'This function computes the ratio of predictable components (RPC; Eade et al., 2014).
+#'
+#'@param exp A numerical array with, at least, 'time_dim' and 'member_dim' 
+#'  dimensions.
+#'@param obs A numerical array with the same dimensions than 'exp' except the
+#'  'member_dim' dimension.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'year'.
+#'@param member_dim A character string indicating the name of the member 
+#'  dimension. The default value is 'member'.
+#'@param na.rm A logical value indicating whether to remove NA values during
+#'  the computation. The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return An array of the ratio of the predictable components. it has the same
+#'  dimensions as 'exp' except 'time_dim' and 'member_dim' dimensions.
+#'
+#'@examples
+#'exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+#'obs <- array(data = runif(60), dim = c(year = 15, lat = 2, lon = 2))
+#'RatioPredictableComponents(exp, obs)
+#'
+#'@import multiApply
+#'@export
+RatioPredictableComponents <- function(exp, obs, time_dim = 'year', member_dim = 'member', na.rm = FALSE, ncores = NULL) {
+
+  C_cor <- stats:::C_cor
+  C_cov <- stats:::C_cov
+    
+  ## Checkings
+  if (is.null(exp)) {
+    stop("Parameter 'exp' cannot be NULL.")
+  }
+  if (!is.numeric(exp)) {
+    stop("Parameter 'exp' must be a numeric array.")
+  }
+  if (is.null(obs)) {
+    stop("Parameter 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(obs)) {
+    stop("Parameter 'obs' must be a numeric array.")
+  }
+  if (!(is.character(time_dim) & length(time_dim) == 1)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!(is.character(member_dim) & length(member_dim) == 1)) {
+    stop("Parameter 'member_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp))) {
+    stop("'exp' must have 'time_dim' dimension.")
+  }
+  if (!member_dim %in% names(dim(exp))) {
+    stop("'exp' must have 'member_dim' dimension.")
+  }
+  if (!time_dim %in% names(dim(obs))) {
+    stop("'obs' must have 'time_dim' dimension.")
+  }
+  if (!identical(dim(exp)[time_dim], dim(obs)[time_dim])) {
+    stop("'exp' and 'obs' must have the same length for 'time_dim'.")
+  }
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE.")
+  }
+  if (!is.null(ncores)){
+    if (!is.numeric(ncores) | ncores %% 1 != 0 |  ncores <= 0 | length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
+  RPC <- multiApply::Apply(data = list(exp, obs), 
+                           target_dims = list(exp = c(time_dim, member_dim),
+                                              obs = time_dim), 
+                           output_dims = NULL,
+                           fun = .RatioPredictableComponents,
+                           na.rm = na.rm,
+                           ncores = ncores)$output1
+  return(RPC)
+}
+
+.RatioPredictableComponents <- function(exp, obs, na.rm = na.rm) {
+  # exp: [time, member]
+  # obs: [time]
+
+  ## Ensemble mean and spread
+  ens_mean <- apply(exp, 1, mean, na.rm = na.rm)
+  ens_spread <- apply(exp, 2, "-", ens_mean)
+
+  ## Ensemble mean variance -> signal
+  var_signal <- var(ens_mean, na.rm = na.rm)
+  ## Variance of ensemble members about ensemble mean (= spread) -> noise
+  # var_noise <- var(as.vector(ens_spread), na.rm = na.rm)
+
+  ## Total variance
+  # var_total <- var_signal + var_noise
+  var_total <- mean(apply(exp, 2, var, na.rm = na.rm), na.rm = na.rm)
+
+  ## Correlation between observations and the ensemble mean
+  r <- cor(ens_mean, obs, method = 'pearson')
+
+  ## Ratio of predictable components
+  RPC <- r / sqrt(var_signal / var_total)
+  
+  return(RPC)
+}
+
diff --git a/R/SignalNoiseRatio.R b/R/SignalNoiseRatio.R
new file mode 100644
index 0000000000000000000000000000000000000000..a8f235637fd074990f0c573def1ee686ededc99b
--- /dev/null
+++ b/R/SignalNoiseRatio.R
@@ -0,0 +1,82 @@
+#'Calculate Signal-to-noise ratio
+#'
+#'This function computes the signal-to-noise ratio, where the signal is the 
+#'ensemble mean variance and the noise is the variance of the ensemble members 
+#'about the ensemble mean (Eade et al., 2014; Scaife and Smith, 2018).
+#'
+#'@param data A numerical array with, at least, 'time_dim' and 'member_dim' 
+#'  dimensions.
+#'@param time_dim A character string indicating the name of the time dimension
+#'  in 'data'. The default value is 'year'.
+#'@param member_dim A character string indicating the name of the member 
+#'  dimension in 'data'. The default value is 'member'.
+#'@param na.rm A logical value indicating whether to remove NA values during the
+#'  computation. The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return An array with of the signal-to-noise ratio. It has the same dimensions
+#'  as 'data' except 'time_dim' and 'member_dim' dimensions.
+#'
+#'@examples
+#' exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+#' SignalNoiseRatio(exp)
+#'
+#'@import multiApply
+#'@export
+SignalNoiseRatio <- function(data, time_dim = 'year', member_dim = 'member', na.rm = FALSE, ncores = NULL) {
+  
+  ## Input Check
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be a numeric array.")
+  }
+  if (!(is.character(time_dim) & length(time_dim) == 1)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!(is.character(member_dim) & length(member_dim) == 1)) {
+    stop("Parameter 'member_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("'data' must have 'time_dim' dimension.")
+  }
+  if (!member_dim %in% names(dim(data))) {
+    stop("'data' must have 'member_dim' dimension.")
+  }
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE.")
+  }
+  if (!is.null(ncores)){
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | length(ncores) > 1){
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
+  SNR <- multiApply::Apply(data = data, 
+                           target_dims = c(time_dim, member_dim), 
+                           output_dims = NULL,
+                           fun = .SignalNoiseRatio,
+                           na.rm = na.rm,
+                           ncores = ncores)$output1
+  return(SNR)
+}
+
+.SignalNoiseRatio <- function(data, na.rm = na.rm) {
+  # data: [time, member]
+
+  ## Ensemble mean and spread
+  ens_mean <- apply(data, 1, mean, na.rm = na.rm)
+  ens_spread <- apply(data, 2, "-", ens_mean)
+
+  ## Ensemble mean variance -> signal
+  var_signal <- var(ens_mean, na.rm = na.rm)
+  ## Variance of ensemble members about ensemble mean (= spread) -> noise
+  var_noise <- var(as.vector(ens_spread), na.rm = na.rm)
+
+  ## Signal to noise ratio
+  SNR <- var_signal / var_noise
+
+  return(SNR)
+}
diff --git a/man/RatioPredictableComponents.Rd b/man/RatioPredictableComponents.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..3e7fbad6c000b42fc6fb1bd689db8007d250b853
--- /dev/null
+++ b/man/RatioPredictableComponents.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RatioPredictableComponents.R
+\name{RatioPredictableComponents}
+\alias{RatioPredictableComponents}
+\title{Calculate ratio of predictable components (RPC)}
+\usage{
+RatioPredictableComponents(
+  exp,
+  obs,
+  time_dim = "year",
+  member_dim = "member",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A numerical array with, at least, 'time_dim' and 'member_dim' 
+dimensions.}
+
+\item{obs}{A numerical array with the same dimensions than 'exp' except the
+'member_dim' dimension.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'year'.}
+
+\item{member_dim}{A character string indicating the name of the member 
+dimension. The default value is 'member'.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values during
+the computation. The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+An array of the ratio of the predictable components. it has the same
+ dimensions as 'exp' except 'time_dim' and 'member_dim' dimensions.
+}
+\description{
+This function computes the ratio of predictable components (RPC; Eade et al., 2014).
+}
+\examples{
+exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+obs <- array(data = runif(60), dim = c(year = 15, lat = 2, lon = 2))
+RatioPredictableComponents(exp, obs)
+
+}
diff --git a/man/SignalNoiseRatio.Rd b/man/SignalNoiseRatio.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f92cf0a1bff6262e0e886c47f153c9f22cfa9989
--- /dev/null
+++ b/man/SignalNoiseRatio.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/SignalNoiseRatio.R
+\name{SignalNoiseRatio}
+\alias{SignalNoiseRatio}
+\title{Calculate Signal-to-noise ratio}
+\usage{
+SignalNoiseRatio(
+  data,
+  time_dim = "year",
+  member_dim = "member",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A numerical array with, at least, 'time_dim' and 'member_dim' 
+dimensions.}
+
+\item{time_dim}{A character string indicating the name of the time dimension
+in 'data'. The default value is 'year'.}
+
+\item{member_dim}{A character string indicating the name of the member 
+dimension in 'data'. The default value is 'member'.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values during the
+computation. The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+An array with of the signal-to-noise ratio. It has the same dimensions
+ as 'data' except 'time_dim' and 'member_dim' dimensions.
+}
+\description{
+This function computes the signal-to-noise ratio, where the signal is the 
+ensemble mean variance and the noise is the variance of the ensemble members 
+about the ensemble mean (Eade et al., 2014; Scaife and Smith, 2018).
+}
+\examples{
+exp <- array(data = runif(600), dim = c(year = 15, member = 10, lat = 2, lon = 2))
+SignalNoiseRatio(exp)
+
+}
diff --git a/tests/testthat/test-RatioPredictableComponents.R b/tests/testthat/test-RatioPredictableComponents.R
new file mode 100644
index 0000000000000000000000000000000000000000..b6a080835fa7fbfd2ccd51c756c2c1b9124cd1fb
--- /dev/null
+++ b/tests/testthat/test-RatioPredictableComponents.R
@@ -0,0 +1,99 @@
+context("RatioPredictableComponents test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  exp1 <- array(rnorm(12), c(dat = 1, member = 2, year = 3, time = 2)) 
+  set.seed(2)
+  obs1 <- array(rnorm(6), c(dat = 1, year = 3, time = 2))
+
+  # dat2
+  exp2 <- exp1
+  exp2[1, 1, 2, 1] <- NA
+  obs2 <- obs1
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+  RatioPredictableComponents(c()),
+  "Parameter 'exp' cannot be NULL."
+  )
+  expect_error(
+    RatioPredictableComponents(c(NA, NA)),
+    "Parameter 'exp' must be a numeric array."
+  )
+  expect_error(
+  RatioPredictableComponents(exp1, c()),
+  "Parameter 'obs' cannot be NULL."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, c(NA, NA)),
+    "Parameter 'obs' must be a numeric array."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, time_dim = 'a'),
+    "'exp' must have 'time_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, member_dim = 'ens'),
+    "'exp' must have 'member_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, array(rnorm(6), dim = c(sdate = 3, time = 2))),
+    "'obs' must have 'time_dim' dimension."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, array(rnorm(4), dim = c(year = 2, time = 2))),
+    "'exp' and 'obs' must have the same length for 'time_dim'."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, na.rm = 'TRUE'),
+    "Parameter 'na.rm' must be TRUE or FALSE."
+  )
+  expect_error(
+    RatioPredictableComponents(exp1, obs1, ncore = 0),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+dim(RatioPredictableComponents(exp1, obs1)),
+c(dat = 1, time = 2)
+)
+expect_equal(
+dim(RatioPredictableComponents(exp1, obs1, time_dim = 'time')),
+c(dat = 1, year = 3)
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp1, obs1)),
+c(-0.29462414, 0.07955797),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp1, obs1, time_dim = 'time')),
+c(-1.054665, 6.843041, -1.000069),
+tolerance = 0.0001
+)
+
+})
+
+test_that("2. Output checks: dat2", {
+
+expect_equal(
+as.vector(RatioPredictableComponents(exp2, obs2, na.rm = T)),
+c(-0.07981636, 0.07955797),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RatioPredictableComponents(exp2, obs2, na.rm = F)),
+c(NA, 0.07955797),
+tolerance = 0.0001
+)
+})
+
diff --git a/tests/testthat/test-SignalNoiseRatio.R b/tests/testthat/test-SignalNoiseRatio.R
new file mode 100644
index 0000000000000000000000000000000000000000..9fc5ce35478b9988d35523c4f19be5562f4185cf
--- /dev/null
+++ b/tests/testthat/test-SignalNoiseRatio.R
@@ -0,0 +1,73 @@
+
+context("SignalNoiseRatio test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  data1 <- array(rnorm(12), c(dat = 1, member = 2, year = 3, time = 2)) 
+
+  # dat2
+  data2 <- data1
+  data2[1, 1, 2, 1] <- NA
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+  SignalNoiseRatio(c()),
+  "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    SignalNoiseRatio(c(NA, NA)),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, time_dim = 'a'),
+    "'data' must have 'time_dim' dimension."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, member_dim = 'ens'),
+    "'data' must have 'member_dim' dimension."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, na.rm = 'TRUE'),
+    "Parameter 'na.rm' must be TRUE or FALSE."
+  )
+  expect_error(
+    SignalNoiseRatio(data1, ncore = 0),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+dim(SignalNoiseRatio(data1)),
+c(dat = 1, time = 2)
+)
+expect_equal(
+as.vector(SignalNoiseRatio(data1)),
+c(0.1591161, 0.8003926),
+tolerance = 0.0001
+)
+
+})
+
+test_that("2. Output checks: dat2", {
+expect_equal(
+as.vector(SignalNoiseRatio(data2, na.rm = T)),
+c(4.5075526, 0.8003926),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(SignalNoiseRatio(data2, na.rm = F)),
+c(NA, 0.8003926),
+tolerance = 0.0001
+)
+
+})
+