What SPI and SPEI measure

The Standardized Precipitation Index (SPI) describes how unusual an accumulated precipitation amount is relative to the precipitation climatology for the same location and calendar month. It was introduced as a probability-based, multiscalar index: only precipitation is required, and the accumulation period is selected to represent the memory of the system being studied (McKee et al., 1993; Guttman, 1998; Hayes et al., 1999).

The Standardized Precipitation-Evapotranspiration Index (SPEI) applies the same standardisation principle to climatic water balance,

\[D_t = P_t - PET_t,\]

where \(P_t\) is precipitation and \(PET_t\) is potential evapotranspiration. Consequently, SPEI responds both to water supply and to atmospheric evaporative demand. This makes it useful when temperature variability or long-term warming affects drought severity, but also makes the result sensitive to the selected PET method and the quality of its meteorological inputs (Vicente-Serrano et al., 2010; Beguería et al., 2014).

For either index, let \(X_t\) denote monthly precipitation for SPI or monthly water balance for SPEI. A backward-looking value at scale \(k\) is

\[Y_{t,k} = \sum_{j=0}^{k-1} w_j X_{t-j},\]

where \(w_j\) are the kernel weights. A separate probability distribution \(F_m\) is fitted to \(Y_{t,k}\) for each calendar month \(m\), because a January water balance should be compared with other Januaries rather than with months from a different season. The fitted non-exceedance probability is then mapped to a standard-normal variate:

\[I_{t,k} = \Phi^{-1}\{H_m(Y_{t,k})\},\]

where \(\Phi^{-1}\) is the standard-normal quantile function and \(H_m\) is either the fitted continuous distribution or, for zero precipitation, a mixed discrete-continuous distribution. If the probability model is well calibrated, the resulting index has an approximately standard-normal interpretation: zero is near the climatological median, negative values are drier than normal, and positive values are wetter than normal. Standard-normal thresholds of -1, -1.5, and -2 correspond approximately to lower-tail probabilities of 15.9%, 6.7%, and 2.3%, respectively.

SPI and SPEI describe statistical rarity, not drought impact by themselves. The same index value can have different consequences depending on season, exposure, water storage, crop stage, and antecedent conditions.

Temporal scales and system memory

The scale is the number of current and preceding months contributing to each index value. It is not the duration of an identified drought event. A drought event may persist for many index values, and neighbouring values are correlated because their accumulation windows overlap.

Common interpretations are:

• 1 month: immediate precipitation or water-balance anomaly and rapidly developing meteorological drought;
• 3 months: short seasonal conditions, soil-moisture stress, and rain-fed agricultural response;
• 6 months: seasonal water availability and slower agricultural or surface-water response;
• 12 months: annual water balance, streamflow, reservoir, and groundwater signals; and
• 24 months or longer: persistent hydrological storage deficits.

These are guidelines, not universal definitions. The appropriate scale should be selected against the response time of the variable or impact of interest. Multiscalar interpretation is one of the principal advantages of both indices (Guttman, 1998; Vicente-Serrano et al., 2010). Comparing several scales is often more informative than choosing a single scale in advance.

Why distribution choice matters

The probability model determines the percentile assigned to each accumulated value and therefore directly controls the final index, especially in the tails. A poor distribution can exaggerate or suppress drought severity even when the fitted mean appears reasonable (Stagge et al., 2015). Distribution choice should respect both the support of the data and the empirical shape:

• SPI precipitation is non-negative, usually right-skewed at short scales, and may contain an exact probability mass at zero. Skewness commonly decreases with longer accumulation.
• SPEI water balance can be negative or positive, so distributions with a location parameter and support extending below zero are generally required.

Stagge et al. (2015) evaluated seven SPI candidates (Gamma, Gumbel, logistic, log-logistic, lognormal, normal, and Weibull) and four SPEI candidates (generalized logistic, GEV, normal, and Pearson type III). They found the two-parameter Gamma to be a strong general SPI choice in their European experiment and recommended GEV for SPEI there. The original SPEI formulation used the three-parameter log-logistic distribution, which is represented here by the equivalent generalized-logistic formulation (Vicente-Serrano et al., 2010). These findings should not be treated as universal winners: performance can change with climate, record length, calendar month, accumulation scale, and observational or modelled data (Pieper et al., 2020).

A defensible selection checks every calendar month and scale, with particular attention to the dry tail. Useful diagnostics include Q-Q or P-P plots, empirical versus theoretical drought-category frequencies, tail-sensitive goodness-of-fit statistics, and validation outside the calibration period. Goodness-of-fit critical values must account for estimated parameters; applying an ordinary Kolmogorov-Smirnov reference distribution after fitting parameters to the same sample is not valid without adjustment or resampling.

Why parameter-estimation method matters

Distribution family and parameter estimator are separate decisions. Two implementations can use the same named distribution but produce different indices because their parameter estimates differ.

• Maximum likelihood (fit="max-lik") selects parameters that maximise the probability of the calibration sample. Under a correct model and standard regularity conditions it is statistically efficient, but it can be sensitive to starting values, outliers, parameter boundaries, and small samples.
• Unbiased probability-weighted moments (fit="ub-pwm") estimate moments that combine observations with their non-exceedance probabilities using unbiased order-statistic weights. PWMs and the related L-moments provide stable estimators for many skewed hydrological distributions (Greenwood et al., 1979; Hosking, 1990). Beguería et al. (2014) recommend unbiased PWM for the SPEI log-logistic model.
• Plotting-position PWM (fit="pp-pwm") replaces the unknown probabilities by empirical plotting positions. hydroTSM uses \(F_i=(i-0.35)/n\), matching the constants used by the SPEI package. This was used in the original SPEI formulation, but it is not identical to unbiased PWM and can differ appreciably in the short monthly calibration samples.

There is no estimator that compensates for a badly chosen distribution. Selection should consider convergence, plausibility of parameters, fit in the tails, and sensitivity of the resulting drought classes. Comparisons between packages are fair only when distribution, estimator, accumulation, reference period, zero treatment, and scaling conventions are equivalent.

Theoretical guide to the function arguments

The arguments below separate scientific choices from numerical and output controls. Defaults are intentionally stated because hydroTSM defaults do not always match those of other R packages.

Differences from the SPEI and SCI packages

At the time of this benchmark, the CRAN versions were SPEI 1.8.1 and SCI 1.0-3. The comparison below concerns their documented public interfaces and the hydroTSM 0.8-8 implementation, rather than suggesting that one design is preferable for every workflow.

The most important comparison consequence is the estimator default: hydroTSM uses maximum likelihood, whereas SPEI defaults to unbiased PWM and SCI fits by maximum likelihood. Equivalent numerical comparisons must therefore set fit explicitly. The generalized-logistic model used by hydroTSM and SCI for SPEI is functionally equivalent to the three-parameter log-logistic formulation used by the SPEI package, but parameter names and signs must be mapped consistently when parameters are supplied.

The current CRAN package descriptions and manuals are available at https://CRAN.R-project.org/package=SPEI and https://CRAN.R-project.org/package=SCI.

Basic usage

library(hydroTSM)

The main argument is a numeric monthly zoo object with a Date, POSIXt, or yearmon index. The observations must cover consecutive months.

The following example uses monthly precipitation from one station in the EbroPPtsMonthly dataset.

data(EbroPPtsMonthly)

pcp <- zoo::zoo(EbroPPtsMonthly$P9001, EbroPPtsMonthly$Date)
head(pcp)
#> 1941-01-01 1941-02-01 1941-03-01 1941-04-01 1941-05-01 1941-06-01 
#>      311.6      158.9       91.2       39.5      135.3       41.7

spi12 <- spi(pcp, scale=12, warn=FALSE)
head(spi12, 15)
#> 1941-01-01 1941-02-01 1941-03-01 1941-04-01 1941-05-01 1941-06-01 1941-07-01 1941-08-01 
#>         NA         NA         NA         NA         NA         NA         NA         NA 
#> 1941-09-01 1941-10-01 1941-11-01 1941-12-01 1942-01-01 1942-02-01 1942-03-01 
#>         NA         NA         NA  1.4114490  0.7628638  0.7692549  1.0063587

spi1 <- spi(pcp, scale=1, zero.threshold=0.1, warn=FALSE)
spi3 <- spi(pcp, scale=3, zero.threshold=0.1, warn=FALSE)
spi6 <- spi(pcp, scale=6, zero.threshold=0.1, warn=FALSE)

spi3.bounded <- spi(
  pcp, scale=3, zero.threshold=0.1,
  p0=TRUE, p0.center.mass=TRUE,
  sci.limit=3, warn=FALSE
)
range(spi3.bounded, na.rm=TRUE)
#> [1] -2.367121  2.035581

SPEI

spei() expects monthly climatic water balance, normally precipitation minus potential evapotranspiration. Its defaults select the generalized-logistic distribution and do not use a zero-probability mixture.

For illustration, a deterministic seasonal PET series is constructed below:

month.number <- as.integer(format(zoo::index(pcp), "%m"))
pet <- zoo::zoo(
  60 + 35 * sin(2 * pi * (month.number - 1) / 12),
  zoo::index(pcp)
)
water.balance <- pcp - pet

spei3 <- spei(water.balance, scale=3, warn=FALSE)
head(spei3)
#> 1941-01-01 1941-02-01 1941-03-01 1941-04-01 1941-05-01 1941-06-01 
#>         NA         NA  1.6907612  0.3953445  0.6236039  0.4519379

SPI and SPEI retain the input time index. Initial values are NA when a complete accumulation window is not yet available.

indices <- zoo::merge.zoo(SPI=spi3, SPEI=spei3)
plot(
  indices, plot.type="single", col=c("#1b9e77", "#d95f02"), lwd=1.2,
  xlab="Date", ylab="Standardized drought index"
)
abline(h=c(-2, -1, 0, 1, 2), col="grey85", lty=c(3, 3, 1, 3, 3))
legend(
  "topright", legend=colnames(indices),
  col=c("#1b9e77", "#d95f02"), lty=1, lwd=1.2, bty="n"
)

Figure 1. Example three-month SPI and SPEI series.

Distributions and parameter estimation

The supported distributions are those evaluated for SPI and SPEI by Stagge et al. (2015). The parameter column shows the required order when params is supplied without names.

Three estimation methods are available:

• fit="max-lik": maximum likelihood;
• fit="ub-pwm": unbiased probability-weighted moments; and
• fit="pp-pwm": plotting-position probability-weighted moments.

For example:

spi3.mle <- spi(pcp, scale=3, fit="max-lik", warn=FALSE)
spi3.ub  <- spi(pcp, scale=3, fit="ub-pwm", warn=FALSE)
spi3.pp  <- spi(pcp, scale=3, fit="pp-pwm", warn=FALSE)

spi3.reference <- spi(
  pcp, scale=3,
  ref.start="1943-01", ref.end="1948-12",
  warn=FALSE
)

spei.fixed.vector <- spei(
  water.balance, scale=1, distribution="norm",
  params=c(mean=0, sd=40), warn=FALSE
)
head(spei.fixed.vector)
#>  1941-01-01  1941-02-01  1941-03-01  1941-04-01  1941-05-01  1941-06-01 
#>  6.28999996  2.03500000  0.02222777 -1.38750000  1.12472777 -0.89500000

monthly.normal.params <- rbind(
  mean=seq(-10, 10, length.out=12),
  sd=rep(40, 12)
)
spei.fixed.monthly <- spei(
  water.balance, scale=1, distribution="norm",
  params=monthly.normal.params, warn=FALSE
)
head(spei.fixed.monthly)
#> 1941-01-01 1941-02-01 1941-03-01 1941-04-01 1941-05-01 1941-06-01 
#>  6.5400001  2.2395455  0.1813187 -1.2738636  1.1929096 -0.8722727

water.balance.two <- zoo::merge.zoo(
  stationA=water.balance,
  stationB=water.balance + 5
)
array.normal.params <- array(
  NA_real_, dim=c(2, 2, 12),
  dimnames=list(c("mean", "sd"), colnames(water.balance.two), month.abb)
)
array.normal.params["mean", "stationA", ] <- 0
array.normal.params["mean", "stationB", ] <- 5
array.normal.params["sd", , ] <- 40

spei.fixed.array <- spei(
  water.balance.two, scale=1, distribution="norm",
  params=array.normal.params, warn=FALSE
)
head(spei.fixed.array)
#>               stationA    stationB
#> 1941-01-01  6.28999996  6.28999996
#> 1941-02-01  2.03500000  2.03500000
#> 1941-03-01  0.02222777  0.02222777
#> 1941-04-01 -1.38750000 -1.38750000
#> 1941-05-01  1.12472777  1.12472777
#> 1941-06-01 -0.89500000 -0.89500000

normal.start <- function(x, distr) {
  c(mean=mean(x), sd=stats::sd(x))
}

spei3.custom <- spei(
  water.balance, scale=3, distribution="norm",
  fit="max-lik", start.fun=normal.start,
  start.fun.fix=FALSE, warn=FALSE
)

invisible(spi(pcp, scale=3, fit="ub-pwm", verbose=TRUE, warn=FALSE))
#> [SPI: scale=3; distribution=gamma; fit=ub-pwm; kernel=rectangular; shift=0]
#> [SPI: reference period 1941-01 to 1950-12]
#> [SPI: processing series 1 of 1 (1)]
#> [SPI: computation completed]

Comparison design

The objective is to compare the SPI and SPEI series computed by the current hydroTSM implementation with series computed by the SPEI and SCI packages. The comparisons are made independently for accumulation scales 1, 3, 6, and 12 and for hydroTSM fit="max-lik", fit="pp-pwm", and fit="ub-pwm".

The benchmark uses the Wichita monthly dataset distributed with SPEI. Climatic water balance was calculated with Thornthwaite PET at latitude 37.6475 degrees north. All implementations used a rectangular kernel and the full record.

The package versions and data period used to create the bundled benchmark are:

All tables and plots use hydroTSM as the first series. The benchmark uses:

• Gamma with mixed zero probability for SPI;
• generalized logistic/log-logistic for SPEI;
• scaling="sd" for hydroTSM and SCI;
• scales 1, 3, 6, and 12; and
• all three hydroTSM fitting methods.

SCI supplies the maximum-likelihood reference. It has no fit argument corresponding to hydroTSM’s PWM alternatives, so comparisons of SCI with ub-pwm or pp-pwm are descriptive rather than equivalent.

For SPEI 1.8.1, the same fit argument was requested from the package. ub-pwm produced finite results. In this dataset, the max-lik and pp-pwm runs produced no finite index values, so those cases are reported as not computable rather than assigned an artificial error statistic.

At scale 1, Wichita contains exact zero precipitation. The two packages treat those zero values differently in this configuration, so the scale-1 SPI comparison with SPEI is not fully equivalent. At scales 3, 6, and 12 there are no zero accumulated totals and the ub-pwm Gamma comparison is equivalent.

Independent comparison tables

Each table fixes one hydroTSM index and one hydroTSM fitting method, then compares that hydroTSM series separately with SCI and SPEI at scales 1, 3, 6, and 12. The hydroTSM result column identifies the series produced by the current hydroTSM function. The reference result column identifies the external package series used as the comparator. Both columns state the package, index, distribution, and fitting method.

The status column identifies the conditions behind each row. Fair means the distribution, fitting method, accumulation, reference period, scaling convention, and zero treatment are equivalent for the comparison being made. Other status values identify why the row is descriptive rather than a strict reproduction test. The nonfinite hydroTSM/reference column reports the number of non-finite values in each input series before pairwise metrics are calculated.

Graphical summaries

Figures 2 to 5 show direct value-to-value agreement for the main fair comparison families. A dashed one-to-one line represents exact agreement.

Figure 2. SPI: hydroTSM maximum likelihood versus SCI maximum likelihood.

Figure 3. SPEI: hydroTSM maximum likelihood versus SCI maximum likelihood.

Figure 4. SPI: hydroTSM unbiased PWM versus the SPEI package. The scale-1 panel has different zero treatment.

Figure 5. SPEI: hydroTSM unbiased PWM versus the SPEI package.

Interpretation

Under equivalent maximum-likelihood conditions, the largest RMSE against SCI was 0.001231 for SPI and 0.0001518 for SPEI. All correlations exceeded 0.99999926.

Under equivalent unbiased-PWM conditions, the largest RMSE against SPEI was 8.286e-09 for SPI at scales 3, 6, and 12, and 3.535e-15 for SPEI over all four scales. These differences are at or close to floating-point precision.

The direct scale-1 SPI comparison with SPEI had correlation 0.9978624 and RMSE 0.07103, but it is not a fair zero-treatment comparison. It must not be combined with the equivalent scale-3, scale-6, and scale-12 results when making a parity claim.

These results support two bounded conclusions for the tested data:

1. hydroTSM closely reproduces SCI when both use maximum likelihood, the same distributions, scaling, accumulation, and zero treatment; and
2. hydroTSM reproduces SPEI’s unbiased-PWM results when the distribution and zero treatment are equivalent.

They do not establish universal equality for other datasets, fitting failures, different reference periods, or intentionally different estimation methods.

Reproducing the external-package benchmark

The benchmark values are bundled so this vignette can be built without making SCI or SPEI package dependencies. A standalone script containing the complete external-package computation is installed with hydroTSM:

comparison.script <- system.file(
  "extdata", "compare_spi_spei.R.txt", package="hydroTSM"
)
source(comparison.script)

SCI, SPEI, and zoo must be installed to run that validation script. By default, it writes new compressed comparison and metadata CSV files to tempdir(). From the hydroTSM source directory, explicit output files can be supplied with:

Rscript inst/extdata/compare_spi_spei.R.txt values.csv.gz metadata.csv.gz

Because SCI and SPEI are comparison oracles rather than production dependencies, updating either package should be followed by regenerating the snapshot and reviewing all changes in finite-value coverage and accuracy.