PRSim (version 1.1)

pRsim: Simulate for one station

Description

Applies the algorithm to a single station

Usage

prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,
      marginal=c("kappa","empirical"), n_par=4, marginalpar=TRUE, 
      GoFtest=NULL, verbose=TRUE, suppWarn=FALSE, ...)

Arguments

data

data frame containing the time indications and runoff of at least one station. See ‘Details’.

station_id

identifies the station in case several runoffs are present in data. See ‘Details’.

number_sim

number of simulations to be carried out.

win_h_length

(half-)length of moving window size.

marginal

marginal distribution to be used for the backtransformation. Can be either "kappa", "empirical", or any type of CDF (see ‘Details’). "kappa" uses the four-parameter kappa distribution for backtransformation, "empirical" uses the empirical distribution. CDF allows for specifying any distribution ‘Examples’.

n_par

number of parameters of the marginal distribution used

GoFtest

If (non-null) a GoF test for daily data should be performed: "KS" performs a Kolmogorof-Smirnov test, and "AD" performs an Anderson-Darling test. see ‘Details’)

verbose

logical. Should progress be reported?

marginalpar

logical. Should the estimated parameters of the distribution used be returned?

suppWarn

logical. See ‘Details’.

...

any other argument passed to the sub-function specifying the cdf for fitting. See ‘Details’ and ‘Examples’.

Value

A list with elements

simulation

A data frame with time information, observations, deseaonalized observations and number_sim columns containing the simulated runoff.

pars

A matrix containing the estimated parameters of the marginal distribution (if marginalpar).

p_val

A vector containing the p-values of ks.test or ad.test applied to the daily detrended data (if GoFtest is not NULL)

Details

Time can be given with three columns named "YYYY", "MM", "DD", or as in POSIXct format YYYY-MM-DD. All leap days (Feb 29th) will be omitted from the analysis, but no missing observations are allowed.

Stations are identified by column name (default "Qobs"), or by column index.

The function homtest::par.kappa might issue quite a few warnings of type In fn(par, ...) : value out of range in 'gammafn'. The argument suppWarn allows to silence warnings for the specific function call via suppressWarnings(). Of course, a subsequent check via warnings() is recommended.

Alternative distributions can be specified by providing three functions: (1) a function fitting the parameters of a distributions and providing a vector of these parameters as output (CDF_fit), (2) a function simulating random numbers from this distribution (rCDF), and (3) a function specifying the distribution (pCDF). See ‘Examples’ for the generalized beta for the second kind and for the Generalized Extreme Values (GEV) distribution.

When using the kappa distribution, the AD test can for certain values of the parameter h not be performed.

References

Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Under review. https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/.

See Also

ks.test

Examples

# NOT RUN {
data( runoff)
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
# warnings() # as a follow-up to `suppWarn=TRUE`

## Specifying particular CDFs:
## (1) example with the Generalized Extreme Value (GEV) distribution
require("evd")
require("ismev")
rGEV <- function(n, theta)   rgev(n, theta[1], theta[2], theta[3])
pGEV <- function(x, theta)   pgev(x, theta[1], theta[2], theta[3])
GEV_fit <- function( xdat, ...)   gev.fit( xdat, ...)$mle
  
# }
# NOT RUN {
 # The following call requires 5 seconds to execute
out <- prsim( runoff[ runoff$YYYY<1978, ], "Qobs", 1, 
    marginal="GEV", n_par=3, verbose=FALSE, marginalpar=FALSE,
    show=FALSE)  # Supress 'gev.fit' output.
# }
# NOT RUN {
## (2) example with generalized Beta distribution of the second kind
require( "GB2")
rGB2 <- function(n, theta)  rgb2(n, theta[1], theta[2], theta[3], theta[4])
pGB2 <- function(x, theta)  pgb2(x, theta[1], theta[2], theta[3], theta[4])
GB2_fit <- function( xdat, ...)   ml.gb2( xdat, ...)$opt1$par

# }
# NOT RUN {
 # The following call requires half minute or so to execute. Some warnings are issued
out <- prsim( runoff[ runoff$YYYY<1987, ], "Qobs", 1, suppWarn=TRUE,
    marginal="GB2")

# }
# NOT RUN {
# }