fitSCI: Standardized Climate Index (SCI)
In SCI: Standardized Climate Indices Such as SPI, SRI or SPEI

Description Usage Arguments Details Value Note Author(s) References See Also Examples

fitSCI identifies parameters for the Standardized Climate Index (SCI) transformation. transformSCI applies the transformation

fitSCI(x, ...)
## Default S3 method:
fitSCI(x, first.mon, time.scale, distr, p0, 
       p0.center.mass=FALSE, scaling=c("no","max","sd"),mledist.par =  list(),
       start.fun = dist.start, start.fun.fix = FALSE, warn = TRUE, ...)

transformSCI(x, ...)
## Default S3 method:
transformSCI(x, first.mon, obj, sci.limit = Inf, warn=TRUE, ...)

`x`	`numeric` vector, representing a monthly univariate time series.
`first.mon`	value in 1:12 indicating the month of the first element of x
`time.scale`	The time scale (`integer`) of the SCI calculation. The time scale is the window length of an backward looking running mean.
`distr`	A character string `"name"` naming a distribution for which the corresponding density function (`dname`), the corresponding distribution function (`pname`) and the quantile function (`qname`) must be defined (see for example `GammaDist`)
`p0`	if TRUE, model Probability of zero (precipitation) months is modeled with a mixed distribution as D(x) = p0 + (1-p0)G(x), where G(x) > 0 is the reference distribution (e.g. Gamma) p0 is the probability of a zero (precipitation) month.
`p0.center.mass`	If TRUE, the Probability of zero (precipitation) is estimated using a "center of mass" estimate based on the Weibull plotting position function (see details). Only applies if `p0=TRUE`.
`scaling`	Indicates whether to do some scaling of `x` prior to parameter identification. `"no"` (the default) indicates no scaling. `"max"` indicates scaling by the maximum of `x`, such that `x <- x/max(x,na.rm=TRUE)`. `"sd"` stands for scaling by the standard deviation. Scaling can stabilize parameter estimation.
`mledist.par`	named `list` that can be used to pass parameters to `mledist` in package fitdistrplus.
`start.fun`	Function with arguments `x` and `distr` estimating initial parameters of the function `distr` for each month. The function should return a named list corresponding to the parameters of `distr`. (See also `dist.start`)
`start.fun.fix`	`logical` argument, indicating if parameter estimates of `start.fun` should be used if maximum likelihood estimation breaks down. This stabilizes the implementation but can introduce biases in the resulting SCI.
`obj`	an object of class `fitSCI`, output from `fitSCI`.
`sci.limit`	Truncate absolute values of SCI that are lage than `sci.limit`. See details.
`warn`	Issue warnings if problems in parameter estimation occur.
`...`	further arguments passed to methods

fitSCI estimates the parameters for transforming a meteorological and environmental time series to a Standardized Climate Index (SCI). transformSCI applies the standardisation. Typical SCI are the Standardized Precipitation Index (SPI), the Standardized Runoff Index (SRI) or the Standardized Precipitation Evapotranspiration Index (SPEI).

To reduce biases in the presence of many zero (precipitation) events, the probability of these events (p0) can be estimated using a "center of mass" estimate based on the Weibull plotting position function (p0.center.mass=TRUE). Following Stagge et al. (2014) the probability of zero events is then estimated as p0 = (n_p)/(n + 1), where np refers to the number of zero events and n is the sample size. The resulting mixed distribution used fro SCI transformation is then

g(x) = if(x > 0) p0 + (1 - p0)G(x) else if(x == 0) (np + 1)/(2(n + 1))

where G(x)>0 is a model (e.g. gamma) distribution.

Uncertainty in distribution parameters can cause unrealistically large (small) SCI values if values in x exceed the values used for parameter estimation (see fitSCI). Therefore transformSCI allows for a truncation of the SCI series such that abs(sci)<=sci.limit. The truncation can be disabled by setting sci.limit=Inf.

fitSCI returns an object of class "fitSCI" with the following components:

`dist.para`	A column `matrix` containing the parameters of distribution `distr` for each month. Row names correspond to the distribution parameters. If `p0=TUE` an additional row named `P0` is introduced, indicating the probability of zero (precipitation) events.
`dist.para.flag`	an vector indicating possible issues occurring throughout parameter estimation. Possible values are: 0. no problems occurred; 1. starting values could not be estimated; 2. `mledist` crashed with unknown error; 3. `mledist` did not converge; 4. all values in this month are `NA`; 5. all values in this month are constant, distribution not defined.
`time.scale`	The time scale (`integer`) of the SCI calculation.
`distr`	A character string `"name"` naming a distribution used
`p0`	`logical` indicating whether probability of zero (precipitation) events is estimated separately.
`p0.center.mass`	`logical` indicating whether probability of zero (precipitation) events is estimated using the "centre of mass" estimator (see Stagge et al. (2014) for details).
`scaling`	`numeric` value that has been used to scale `x` (see argument `scaling`). A value of 1 results from `scaling="no"`, other values are the maximum value or the standard deviation of `x`, depending on the choice of the parameter `scaling`.
`call`	the function call

transformSCI returns a numeric vector containing the SCI, having values of the standard normal distribution.

This function is intended to be used together with transformSCI.

Lukas Gudmundsson & James Stagge

Stagee, J.H. ; Tallaksen, L.M.; Gudmundsson, L.; van Loon, A.; Stahl, K.: Candidate Distributions for Climatological Drought Indices (SPI and SPEI), 2015, International Journal of Climatology, 35, 4027-4040, doi:10.1002/joc.4267.

Stagee, J.H. ; Tallaksen, L.M.; Gudmundsson, L.; van Loon, A.; Stahl, K.: Response to comment on "Candidate Distributions for Climatological Drought Indices (SPI and SPEI)", 2016, International Journal of Climatology, 36, 2132-2138, doi:10.1002/joc.4564.

McKee, T.; Doesken, N. & Kleist, J.: The relationship of drought frequency and duration to time scales Preprints, 8th Conference on Applied Climatology, 1993, 179-184.

Shukla, S. & Wood, A. W.: Use of a standardized runoff index for characterizing hydrologic drought Geophysical Research Letters, 2008, 35, L02405.

Vicente-Serrano, S. M.; Begueria, S. & Lopez-Moreno, J. I.: A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index J. Climate, Journal of Climate, American Meteorological Society, 2009, 23, 1696-1718.

dist.start

##
## generate artificial data
##
set.seed(101)
n.years <- 60
date <- rep(1:n.years,each=12) + 1950 + rep((0:11)/12,times=n.years)
## Precipitation
PRECIP <- (0.25*sin( 2 * pi * date) + 0.3)*rgamma(n.years*12, shape = 3, scale = 1)
PRECIP[PRECIP<0.1] <- 0
## Potential Evapotranspiration
PET <- 0.5*sin( 2 * pi * date)+1.2+rnorm(n.years*12,0,0.2)
## display test data
matplot(date,cbind(PRECIP,PET),t=c("h","l"),col=c("blue","red"),lty=1)
legend("topright",legend=c("PRECIPitation","temperature"),fill=c("blue","red"))

##
## example SPI
##
spi.para <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
spi.para
spi <- transformSCI(PRECIP,first.mon=1,obj=spi.para)
plot(date,spi,t="l")

##
## effect of time.scale on SPI
##
spi.1.para <- fitSCI(PRECIP,first.mon=1,time.scale=1,distr="gamma",p0=TRUE)
spi.12.para <- fitSCI(PRECIP,first.mon=1,time.scale=12,distr="gamma",p0=TRUE)
spi.1 <- transformSCI(PRECIP,first.mon=1,obj=spi.1.para)
spi.12 <- transformSCI(PRECIP,first.mon=1,obj=spi.12.para)
matplot(date,cbind(spi.1,spi.12),t="l",lty=1,col=c("red","blue"),lwd=c(1,2))
legend("topright",legend=c("time.scale=1","time.scale=12"),fill=c("red","blue"))

##
## example SPEI
##
if(require(evd)){
    spei.para <- fitSCI(PRECIP-PET,first.mon=1,time.scale=6,distr="gev",p0=FALSE)
    spei <- transformSCI(PRECIP-PET,first.mon=1,obj=spei.para)
    plot(date,spei,t="l")
}

##
## effect of changing different distribution for SPEI computation
##
spei.genlog.para <- fitSCI(PRECIP-PET,first.mon=1,time.scale=6,distr="genlog",p0=FALSE)
spei.genlog <- transformSCI(PRECIP-PET,first.mon=1,obj=spei.genlog.para)
if(require(evd)){lines(date,spei.genlog, col="red")} else {plot(date,spei.genlog,t="l")}
## in this case: only limited effect.
## generally: optimal choice of distribution: user responsibility.

##
## use a 30 year reference period for SPI parameter estimation
##
sel.date <- date>=1970 & date<2000
spi.ref.para <- fitSCI(PRECIP[sel.date],first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
## apply the the parameters of the reference period to all data
## also outside the reference period
spi.ref <- transformSCI(PRECIP,first.mon=1,obj=spi.ref.para)
plot(date,spi.ref,t="l",col="blue",ylim=c(-5,5),lwd=2)
lines(date[sel.date],spi.ref[sel.date],col="red",lwd=3)
legend("bottom",legend=c("reference period","extrapolation period"),fill=c("red","blue"),
       horiz=TRUE)

##
## use "start.fun.fix" in instances where maximum likelyhood estimation fails
##
## force failure of maximum likelyhood estimation by adding "strange" value
## a warning should be issued
xx <- PRECIP-PET; xx[300] <- 1000
spei.para <- fitSCI(xx,first.mon=2,time.scale=1,p0=FALSE,distr="gev")
spei.para$dist.para
## use start.fun, usually ment for estimating inital values for
## parameter optimisation if maximum likelihood estimation fails
spei.para <- fitSCI(xx,first.mon=2,time.scale=1,p0=FALSE,distr="gev",
                    start.fun.fix=TRUE)
spei.para$dist.para

##
## usage of sci.limit to truncate unrealistic SCI values
## 
PRECIP.mod <- PRECIP
PRECIP.mod[300] <- 100 ## introduce spuriously large value
spi.mod.para <- fitSCI(PRECIP.mod,first.mon=1,time.scale=3,p0=TRUE,distr="gamma")
plot(transformSCI(PRECIP.mod,first.mon=1,obj=spi.mod.para,sci.limit=Inf),
     t="l",col="blue",lwd=2)
lines(transformSCI(PRECIP.mod,first.mon=1,obj=spi.mod.para,sci.limit=4),col="red")

##
## how to modify settings of function "mledist" used for parameter identification
## 
## identify parameters with standard settings
spi.para <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
## add lower and upper limits for parameter identification
lower.lim <- apply(spi.para$dist.para,1,min) - 0.5*apply(spi.para$dist.para,1,sd)
upper.lim <- apply(spi.para$dist.para,1,max) + 0.5*apply(spi.para$dist.para,1,sd)
spi.para.limit <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE,
                         mledist.par=list(lower=lower.lim, upper=upper.lim))

##
## how to write an own start.fun
## (required if distributions not mentioned in "dist.start" are used)
##
## function with same arguments as "dist.start" 
my.start <- function(x,distr="gamma"){
### code based on "mmedist" in package "fitdistrplus"
    ppar <- try({
        n <- length(x)
        m <- mean(x)
        v <- (n - 1)/n * var(x)
        shape <- m^2/v
        rate <- m/v
        list(shape = shape, rate = rate)},TRUE)
    if (class(ppar) == "try-error") ## function has to be able to return NA parameters
        ppar <- list(shape = NA, rate = NA)
    return(ppar)
}
my.start(PRECIP)
spi.para <- fitSCI(PRECIP,first.mon=1,time.scale=6,p0=TRUE,
                   distr="gamma",start.fun=my.start)

Loading required package: fitdistrplus
Loading required package: MASS
Loading required package: survival
Loading required package: lmomco
$dist.para
            M1       M2       M3       M4       M5       M6       M7       M8
shape 17.26002 12.06338 11.87662 14.43881 17.15217 16.96217 19.35115 18.24642
rate  40.58648 21.93490 15.56737 13.88936 13.67164 12.44925 14.01031 14.52150
P0     0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
            M9      M10      M11      M12
shape 13.14131 11.98158 12.62739 13.78922
rate  12.57110 15.58746 22.90729 31.16647
P0     0.00000  0.00000  0.00000  0.00000

$dist.para.flag
 M1  M2  M3  M4  M5  M6  M7  M8  M9 M10 M11 M12 
  0   0   0   0   0   0   0   0   0   0   0   0 

$time.scale
[1] 6

$distr
[1] "gamma"

$p0
[1] TRUE

$p0.center.mass
[1] FALSE

$scaling
[1] 1

$call
fitSCI.default(x = PRECIP, first.mon = 1, time.scale = 6, distr = "gamma", 
    p0 = TRUE)

attr(,"class")
[1] "fitSCI"
Loading required package: evd

Attaching package: 'evd'

The following object is masked from 'package:lmomco':

    pp

Warning message:
In fitSCI.default(xx, first.mon = 2, time.scale = 1, p0 = FALSE,  :
  maximum likelihood estimation failed for month 1
 parameters are set to NA.
      M1          M2          M3         M4           M5          M6
loc   NA -0.55849573 -0.50301068 -0.5835045 -0.409466891 -0.49832445
scale NA  0.45547909  0.59330109  0.6095776  0.790411104  0.67581457
shape NA -0.06179259 -0.05244057  0.1298105  0.008576031  0.04723189
               M7         M8         M9        M10        M11        M12
loc   -0.57548855 -0.4402906 -0.6148371 -0.5782434 -0.6970006 -0.6202785
scale  0.59705376  0.5540465  0.3030101  0.2469026  0.2805855  0.2545655
shape  0.03251138 -0.1631595 -0.1024841 -0.3290109 -0.2429432 -0.2780727
Warning message:
In fitSCI.default(xx, first.mon = 2, time.scale = 1, p0 = FALSE,  :
  maximum likelihood estimation failed for month 1
 use starting values instead.
              M1          M2          M3         M4           M5          M6
loc   -0.6356714 -0.55849573 -0.50301068 -0.5835045 -0.409466891 -0.49832445
scale  0.1814620  0.45547909  0.59330109  0.6095776  0.790411104  0.67581457
shape  0.9893455 -0.06179259 -0.05244057  0.1298105  0.008576031  0.04723189
               M7         M8         M9        M10        M11        M12
loc   -0.57548855 -0.4402906 -0.6148371 -0.5782434 -0.6970006 -0.6202785
scale  0.59705376  0.5540465  0.3030101  0.2469026  0.2805855  0.2545655
shape  0.03251138 -0.1631595 -0.1024841 -0.3290109 -0.2429432 -0.2780727
Warning message:
In fitSCI.default(PRECIP, first.mon = 1, distr = "gamma", time.scale = 6,  :
  maximum likelihood estimation failed for month 12
 parameters are set to NA.
$shape
[1] 1.222377

$rate
[1] 1.354323