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

## Description

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

## Usage

 ```1 2 3 4 5 6 7 8 9``` ```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, ...) ```

## Arguments

 `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

## Details

`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`.

## Value

`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.

## Note

This function is intended to be used together with `transformSCI`.

## Author(s)

Lukas Gudmundsson & James Stagge

## References

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`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120``` ```## ## 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) ```

### Example output

```Loading required package: fitdistrplus
\$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"

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
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
```

SCI documentation built on May 2, 2019, 3:04 a.m.