## load packages required

rm(list = ls())

library(tidyverse)

library(data.table)

library(SPEI)

## define drought indices

dIndex <- list()

dIndex$spi <-

  function(ds, scale)

    SPEI::spi(ds, scale = scale) %>% getElement("fitted") %>% as.numeric

dIndex$spei <-

  function(rainfall, ET0, scale = 1)

    SPEI::spei(rainfall - ET0, scale = scale) %>% getElement("fitted") %>% as.numeric

dIndex$CZI <- function(ds, scale = 1) {

  innerCZI <- function(xx) {

    idx <- !is.na(xx)

    txx <- xx[idx]

    mean.X <- mean(txx)

    sigma.X <- sqrt(sum((txx - mean.X) ^ 2) / length(txx))

    Cs <- sum((txx - mean.X) ^ 3) / (length(txx) * sigma.X ^ 3)

    Phi <- (txx - mean.X) / sigma.X

    xx[idx] <-

      6 / Cs * abs(Cs / 2 * Phi + 1) ^ (1 / 3) - 6 / Cs + Cs / 6

    xx

  }

  .tmp <- ts(ds, frequency = 12)

  .tmproll <- frollsum(.tmp, n = scale)

  .tmpds <-

    data.frame(ID = seq_along(.tmproll),

               cycle = cycle(.tmp),

               tmproll = .tmproll) %>% arrange(ID, cycle) %>%

    group_by(cycle) %>% mutate(Z = innerCZI(tmproll))

  .tmpds$Z

}

## create a dataset for testing

ET0 <- runif(12*50,min = 0,max = 10) # monthly ET0

rainfall <- runif(12*50,min = 0,max = 300) # monthly Rainfall

ds <- expand.grid(Y=1971:2020,M=1:12) %>% arrange(Y,M) %>% 

  dplyr::mutate(ET0 =ET0,rainfall=rainfall)

##  Drought indices (Z, SPI and SPEI) calculated at different time scales.

ds %>% dplyr::mutate(

  SPI_1=dIndex$spi(ds =rainfall,scale = 1),    ## scale =1

  SPEI_1=dIndex$spei(rainfall = rainfall,ET0 = ET0,scale = 1),

  Z_1=dIndex$CZI(ds =rainfall,scale = 1),

  SPI_3=dIndex$spi(ds =rainfall,scale = 3),   ## scale =3

  SPEI_3=dIndex$spei(rainfall = rainfall,ET0 = ET0,scale = 3),

  Z_3=dIndex$CZI(ds =rainfall,scale = 3),

  SPI_12=dIndex$spi(ds =rainfall,scale =12),   ## scale =12

  SPEI_12=dIndex$spei(rainfall = rainfall,ET0 = ET0,scale = 12),

  Z_12=dIndex$CZI(ds =rainfall,scale = 12)

)

计算结果（前20行）如下：