Commit 9f8b0952 authored by Reinhard Furrer's avatar Reinhard Furrer

From PRSim/PRSim

parent e38d6238
......@@ -17,7 +17,9 @@ Description: Provides a simulation framework to simulate streamflow time series
streamflow values and their temporal correlation as expressed by short- and long-range
dependence. The approach is based on the randomization of the phases of the Fourier
transform. We further use the flexible four-parameter Kappa distribution, which allows
for the extrapolation to yet unobserved low and high flows.
for the extrapolation to yet unobserved low and high flows. A detailed description of
the simulation approach and an application example can be found
in https://www.hydrol-earth-syst-sci-discuss.net/hess-2019-142/.
URL: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
BugReports: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
License: GPL-3
......
prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15, marginal="kappa",
verbose=TRUE, kappapar=TRUE, suppWarn=FALSE, KStest=FALSE){
ifft <- function (x) fft(x, inverse = TRUE)/length(x)
......@@ -79,66 +79,67 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### daily fitting of Kappa distribution
### fit the parameters of the Kappa distribution for each day separately.
### To enable a sufficient sample size by using daily values in moving window around day i (i.e., reduce uncertainty due to fitting)
### try Kappa distribution
p_vals <- numeric(365)
kap_par_day <- matrix(0, nrow=365, ncol=4)
density_kap <- list()
### define window length
win_length <- c(1:win_h_length)
for(d in c(1:365)){
### define start and end of window
before <- data$index[d+365-win_length]
after <- data$index[d+365+win_length-1]
### define days within window
ids <- c(before, after)
### determine values in window around day i
data_window <- data$Qobs[which(data$index%in%ids)]
# par.kappa(data_monthly)
ll<- Lmoments(data_window)
### test whether Kappa distribution can be fit
if (suppWarn) {
suppressWarnings( test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5])) )
} else {
test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5]))
}
if(length(test)>1){
kap_par <- test
kap_par_day[d,] <- unlist(kap_par)
### define vector of quantiles
quant <- sort(data_window)
thresh <- kap_par$xi + kap_par$alfa*(1 - kap_par$h^(-kap_par$k))/kap_par$k
if(!is.na(thresh)){
## min(quant)>thresh
### only use quantiles larger than threshold (as in f.kappa function documentation)
quant <- quant[which(quant>thresh)]
if(marginal=="kappa"){
p_vals <- numeric(365)
kap_par_day <- matrix(0, nrow=365, ncol=4)
density_kap <- list()
### define window length
win_length <- c(1:win_h_length)
for(d in c(1:365)){
### define start and end of window
before <- data$index[d+365-win_length]
after <- data$index[d+365+win_length-1]
### define days within window
ids <- c(before, after)
### determine values in window around day i
data_window <- data$Qobs[which(data$index%in%ids)]
# par.kappa(data_monthly)
ll<- Lmoments(data_window)
### test whether Kappa distribution can be fit
if (suppWarn) {
suppressWarnings( test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5])) )
} else {
test <- try(par.kappa(ll[1],ll[2],ll[4],ll[5]))
}
# kappa_density <- f.kappa(x=quant,xi=kap_par$xi,alfa=kap_par$alfa,k=kap_par$k,h=kap_par$h)
# plot(kappa_density)
data_kap <- rand.kappa(length(data_window), xi=kap_par$xi,alfa=kap_par$alfa, k=kap_par$k, h=kap_par$h)
# density_kap[[d]] <- density(data_kap)
# hist(data_window)
# hist(data_kap,add=T,col="red")
if (KStest)
p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
} else{
if(d==1){
p_vals[d] <- NA
kap_par_day[d,] <- NA
}else{
p_vals[d] <- p_vals[d-1]
kap_par_day[d,] <- kap_par_day[d-1,]
if(length(test)>1){
kap_par <- test
kap_par_day[d,] <- unlist(kap_par)
### define vector of quantiles
quant <- sort(data_window)
thresh <- kap_par$xi + kap_par$alfa*(1 - kap_par$h^(-kap_par$k))/kap_par$k
if(!is.na(thresh)){
## min(quant)>thresh
### only use quantiles larger than threshold (as in f.kappa function documentation)
quant <- quant[which(quant>thresh)]
}
# kappa_density <- f.kappa(x=quant,xi=kap_par$xi,alfa=kap_par$alfa,k=kap_par$k,h=kap_par$h)
# plot(kappa_density)
data_kap <- rand.kappa(length(data_window), xi=kap_par$xi,alfa=kap_par$alfa, k=kap_par$k, h=kap_par$h)
# density_kap[[d]] <- density(data_kap)
# hist(data_window)
# hist(data_kap,add=T,col="red")
if (KStest)
p_vals[d] <- ks.test(data_window, data_kap)$p.value ### kappa distribution not rejected at alpha=0.05
} else{
if(d==1){
p_vals[d] <- NA
kap_par_day[d,] <- NA
}else{
p_vals[d] <- p_vals[d-1]
kap_par_day[d,] <- kap_par_day[d-1,]
}
}
}
}
###qu: check if subsequent NAs work properly.
### replace NA entries by values of subsequent day
if(length(which(is.na(kap_par_day[,1])))>0){
indices <- rev(which(is.na(kap_par_day[,1])))
for(i in 1:length(indices)){
kap_par_day[indices[i],] <- kap_par_day[indices[i]+1,]
###qu: check if subsequent NAs work properly.
### replace NA entries by values of subsequent day
if(length(which(is.na(kap_par_day[,1])))>0){
indices <- rev(which(is.na(kap_par_day[,1])))
for(i in 1:length(indices)){
kap_par_day[indices[i],] <- kap_par_day[indices[i]+1,]
}
}
}
......@@ -232,47 +233,67 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### e) backtransform from normal to actual distribution
### apply daily backtransformation
colnames( kap_par_day) <- names(kap_par)
data_new$simulated_seasonal <- NA
for(d in c(1:365)){
# ### use empirical distribution for backtransformation
data_day <- data[which(data$index%in%c(d)),]
# data_month$rank <- rank(data_month$Qobs)
# data_new$rank[which(data$MM%in%c(m))] <- rank(data_new[which(data$MM%in%c(m)),]$seasonal)
# ### derive corresponding values from the empirical distribution
# ### identify value corresponding to rank in the original time series
# data_ordered <- data_month[order(data_month$rank),]
# data_new$simulated_seasonal[which(data_new$MM%in%c(m))] <- data_ordered$Qobs[data_new$rank[which(data$MM%in%c(m))]]
# }
# if(marginal=="kappa"){
### use monthly Kappa distribution for backtransformation
### simulate random sample of size n from Kappa disribution
data_day$kappa <- rand.kappa(length(data_day$Qobs),
xi=kap_par_day[d,"xi"],alfa=kap_par_day[d,"alfa"],
k=kap_par_day[d,"k"],h=kap_par_day[d,"h"])
# }
# if(marginal=="wakeby"){
# ### use Wakeby distribution for backtrasformation
# ### simulate random sample of size n from Wakeby distribution
# data_day$kappa <- rlmomco(length(data_day$Qobs),wak_par_day[[d]])
# }
data_day$rank <- rank(data_day$kappa)
##QUESTION: had to recopy the line here: possibly more not needed below..
data_new$rank <- rank(data_new$seasonal)
data_new$rank[ which(data$index%in%c(d)) ] <- rank(data_new[ which(data$index%in%c(d)), ]$seasonal)
### derive corresponding values from the kappa distribution
### identify value corresponding to rank in the kappa time series
data_ordered <- data_day[order(data_day$rank),]
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$kappa[data_new$rank[which(data$index%in%c(d))]]
### if error was applied, replace negative values by 0 values
### in any case, replace negative values by 0. Corresponds to a bounded Kappa distribution
# if(error==TRUE){
if(length(which(data_new$simulated_seasonal<0))>0){
data_new$simulated_seasonal[which(data_new$simulated_seasonal<0)] <- 0
}
# }
data_new$simulated_seasonal <- NA
for(d in c(1:365)){
data_day <- data[which(data$index%in%c(d)),]
# data_month$rank <- rank(data_month$Qobs)
# data_new$rank[which(data$MM%in%c(m))] <- rank(data_new[which(data$MM%in%c(m)),]$seasonal)
# ### derive corresponding values from the empirical distribution
# ### identify value corresponding to rank in the original time series
# data_ordered <- data_month[order(data_month$rank),]
# data_new$simulated_seasonal[which(data_new$MM%in%c(m))] <- data_ordered$Qobs[data_new$rank[which(data$MM%in%c(m))]]
# }
### use kappa distribution for backtransformation
if(marginal=="kappa"){
colnames(kap_par_day) <- names(kap_par)
### use monthly Kappa distribution for backtransformation
### simulate random sample of size n from Kappa disribution
data_day$kappa <- rand.kappa(length(data_day$Qobs),
xi=kap_par_day[d,"xi"],alfa=kap_par_day[d,"alfa"],
k=kap_par_day[d,"k"],h=kap_par_day[d,"h"])
# }
# if(marginal=="wakeby"){
# ### use Wakeby distribution for backtrasformation
# ### simulate random sample of size n from Wakeby distribution
# data_day$kappa <- rlmomco(length(data_day$Qobs),wak_par_day[[d]])
# }
data_day$rank <- rank(data_day$kappa)
##QUESTION: had to recopy the line here: possibly more not needed below..
data_new$rank <- rank(data_new$seasonal)
data_new$rank[ which(data$index%in%c(d)) ] <- rank(data_new[ which(data$index%in%c(d)), ]$seasonal)
### derive corresponding values from the kappa distribution
### identify value corresponding to rank in the kappa time series
data_ordered <- data_day[order(data_day$rank),]
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$kappa[data_new$rank[which(data$index%in%c(d))]]
### if error was applied, replace negative values by 0 values
### in any case, replace negative values by 0. Corresponds to a bounded Kappa distribution
# if(error==TRUE){
if(length(which(data_new$simulated_seasonal<0))>0){
### do not use 0 as a replacement value directly
# data_new$simulated_seasonal[which(data_new$simulated_seasonal<0)] <- 0
### sample value from a uniform distribution limited by 0 and the minimum observed value
### determine replacement value
rep_value <- runif(n=1,min=0,max=min(data_day$Qobs))
data_new$simulated_seasonal[which(data_new$simulated_seasonal<0)]<-rep_value
}
# }
}
### use empirical distribution for backtransformation
if(marginal=="empirical"){
data_day$rank <- rank(data_day$Qobs)
data_new$rank <- rank(data_new$seasonal)
data_new$rank[which(data$index%in%c(d))] <- rank(data_new[which(data$index%in%c(d)),]$seasonal)
### derive corresponding values from the empirical distribution
### identify value corresponding to rank in the original time series
data_ordered <- data_day[order(data_day$rank),]
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$Qobs[data_new$rank[which(data$index%in%c(d))]]
# }
}
}
data_sim[[r]] <- data_new$simulated_seasonal
......@@ -292,9 +313,13 @@ prsim <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
p_vals <- NULL
}
if (kappapar) { # also return intermediate results
return(list( simulation=data_stoch, kappa_pars=kap_par_day, p_val=p_vals))
} else {
return(list( simulation=data_stoch, kappa_pars=NULL, p_val=p_vals))
if(marginal=="kappa"){
if (kappapar) { # also return intermediate results
return(list( simulation=data_stoch, kappa_pars=kap_par_day, p_val=p_vals))
} else {
return(list( simulation=data_stoch, kappa_pars=NULL, p_val=p_vals))
}
}else{
return(list( simulation=data_stoch) )
}
}
# short demo
data(runoff)
out <- prsim(runoff, "Qobs", 10)
out <- prsim(data=runoff, number_sim=1, marginal="empirical")
out <- prsim(data=runoff, number_sim=1, marginal="kappa")
sim <- out$simulation
par(mai=c(.9,.9,.1,.1))
plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
......
......@@ -6,7 +6,7 @@
Applies the algorithm to a single station
}
\usage{
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,marginal="kappa",
verbose=TRUE, kappapar=TRUE, suppWarn=FALSE, KStest=FALSE)
}
\arguments{
......@@ -14,6 +14,7 @@ prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,
\item{station_id}{identifies the station in case several runoffs are present in \code{data}. See \sQuote{Details}.}
\item{win_h_length}{(half-)length of moving window size.}
\item{number_sim}{number of simulations to be carried out.}
\item{marginal}{marginal distribution to be used for the backtransformation.}
\item{verbose}{logical. Should progress be reported?}
\item{kappapar}{logical. Should the estimated parameters of the kappa distribution be returned?}
\item{suppWarn}{logical. See \sQuote{Details}.}
......@@ -37,7 +38,7 @@ The function \code{homtest::par.kappa} might issue quite a few warnings of type
\item{p_val}{A vector containing the p-values of \code{ks.test} applied to the daily detrended data (if \code{KStest}).}
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
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/.
}
\author{
Manuela Brunner
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment