| Title: | Multi-Site Stochastic Models for Daily Precipitation and Temperature |
|---|---|
| Description: | Application of multi-site models for daily precipitation and temperature data. This package is designed for an application to 105 precipitation and 26 temperature gauges located in Switzerland. It applies fitting procedures and provides weather generators described in the following references: - Evin, G., A.-C. Favre, and B. Hingray. (2018) <doi:10.5194/hess-22-655-2018>. - Evin, G., A.-C. Favre, and B. Hingray. (2018) <doi:10.1007/s00704-018-2404-x>. |
| Authors: | Guillaume Evin [aut, cre] |
| Maintainer: | Guillaume Evin <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.5 |
| Built: | 2025-02-04 05:56:08 UTC |
| Source: | https://github.com/guillaumeevin/gwex |
Simple accumulation of a matrix of precipitation
agg.matrix(mat, k, average = F)agg.matrix(mat, k, average = F)
mat |
matrix nDates x nStations to be aggregated |
k |
number of days for the accumulation |
average |
logical: should we average over the different periods (default=F) |
aggregated matrix
Guillaume Evin
Finds empirical autocorrelations (lag-1) between intensities corresponding to a degree of autocorrelation of an AR(1) process
autocor.emp.int(rho, nChainFit, Xt, parMargin, typeMargin)autocor.emp.int(rho, nChainFit, Xt, parMargin, typeMargin)
rho |
autocorrelation of the AR(1) process |
nChainFit |
number of simulated variates |
Xt |
simulated occurrences, nChainFit x 2 matrix |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
scalar |
correlation between simulated intensities |
Guillaume Evin
Finds observed correlations between intensities corresponding to a degree of correlation of Gaussian multivariate random numbers
cor.emp.int(zeta, nChainFit, Xt, parMargin, typeMargin)cor.emp.int(zeta, nChainFit, Xt, parMargin, typeMargin)
zeta |
correlation of Gaussian multivariates |
nChainFit |
number of simulated variates |
Xt |
simulated occurrences, n x 2 matrix |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
scalar |
correlation between simulated intensities |
Guillaume Evin
Finds observed correlations between occurrences corresponding to a degree of correlation of Gaussian multivariate random numbers
cor.emp.occ(w, Qtrans.mat, mat.comb, nLag, nChainFit, myseed = 1)cor.emp.occ(w, Qtrans.mat, mat.comb, nLag, nChainFit, myseed = 1)
w |
correlation of Gaussian multivariates |
Qtrans.mat |
transition probabilities, 2 x ncomb matrix |
mat.comb |
matrix of logical: ncomb x nlag |
nLag |
order of the Markov chain |
nChainFit |
number of simulated variates |
myseed |
seed of random variates |
scalar |
correlation between occurrences |
Guillaume Evin
provide observed correlations between occurrences for all pairs of stations see Mhanna et al. (2012)
cor.obs.occ(pi00, pi0, pi1)cor.obs.occ(pi00, pi0, pi1)
pi00 |
joint probability of having dry states |
pi0 |
probability of having a dry state |
pi1 |
probability of having a wet state |
scalar |
matrix of observed correlations |
Guillaume Evin
Mhanna, Muamaraldin, and Willy Bauwens. “A Stochastic Space-Time Model for the Generation of Daily Rainfall in the Gaza Strip.” International Journal of Climatology 32, no. 7 (June 15, 2012): 1098–1112. doi:10.1002/joc.2305.
Example of daily observations of precipitation (mm) for three fictive stations, for a period of ten years.
data(dailyPrecipGWEX)data(dailyPrecipGWEX)
matrix of Observed precipitation: 3652 days x 3 stations
Guillaume Evin [email protected]
Evin, G., A.-C. Favre, and B. Hingray. 2018. “Stochastic Generation of Multi-Site Daily Precipitation Focusing on Extreme Events". Hydrol. Earth Syst. Sci. 22 (1): 655–672.
Example of daily observations of temperature (mm) for three fictive stations, for a period of ten years.
data(dailyTemperGWEX)data(dailyTemperGWEX)
matrix of Observed temperature: 3652 days x 3 stations
Guillaume Evin [email protected]
Evin G., A.C. Favre, and B. Hingray. 2018. Stochastic Generators of Multi Site Daily Temperature: Comparison of Performances in Various Applications. Theoretical and Applied Climatology.
disag.3D.to.1D
disag.3D.to.1D(Yobs, YObsAgg, mObsAgg, YSimAgg, mSimAgg, prob.class)disag.3D.to.1D(Yobs, YObsAgg, mObsAgg, YSimAgg, mSimAgg, prob.class)
Yobs |
matrix of observed intensities at 24h: (nTobs*3) x nStation |
YObsAgg |
matrix of observed 3-day intensities: nTobs x nStation |
mObsAgg |
vector of season corresponding to YobsAgg |
YSimAgg |
matrix of simulated intensities per 3-day period: nTsim x nStation |
mSimAgg |
vector of season corresponding to the period simulated |
prob.class |
vector of probabilities indicating class of "similar" mean intensities |
list |
Ysim matrix of disagregated daily precipitation, codeDisag matrix of disagregation codes |
Guillaume Evin
Density function, distribution function, quantile function, random generation for the unified EGPD distribution
dEGPD.GI(x, kappa, sig, xi) pEGPD.GI(x, kappa, sig, xi) qEGPD.GI(p, kappa, sig, xi) rEGPD.GI(n, kappa, sig, xi)dEGPD.GI(x, kappa, sig, xi) pEGPD.GI(x, kappa, sig, xi) qEGPD.GI(p, kappa, sig, xi) rEGPD.GI(n, kappa, sig, xi)
x |
Vector of quantiles |
kappa |
transformation parameter greater than 0 |
sig |
Scale parameter |
xi |
Shape parameter |
p |
Vector of probabilities |
n |
Number of observations |
dEGPD.GI gives the density function, pEGPD.GI gives the distribution function, qEGPD.GI gives the quantile function, and rEGPD.GI generates random deviates.
Guillaume Evin
Estimate the dry day frequency (proportion of dry days) for all stations
dry.day.frequency(mat.prec, th)dry.day.frequency(mat.prec, th)
mat.prec |
matrix of precipitation (possibly for one month/period) |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
vector of numeric |
dry day frequencies |
Guillaume Evin
Parameter estimation of the unified EGPD distribution with the PWM method. Numerical solver of the system of nonlinear equations
EGPD.GI.fit.PWM(x, xi = 0.05)EGPD.GI.fit.PWM(x, xi = 0.05)
x |
vector of parameters kappa,sig |
xi |
shape parameter |
estimated parameters kappa, sig, xi
Guillaume Evin
Parameter estimation of the unified EGPD distribution with the PWM method. Set of equations which have to be equal to zero
EGPD.GI.fPWM(par, pwm, xi)EGPD.GI.fPWM(par, pwm, xi)
par |
vector of parameters kappa,sig,xi |
pwm |
set of probability weighted moments of order 0, 1 and 2 |
xi |
shape parameter |
differences between expected and target weighted moments
Guillaume Evin
finds the autocorrelation leading to observed autocorrelation
find.autocor(autocor.emp, nChainFit, Xt, parMargin, typeMargin)find.autocor(autocor.emp, nChainFit, Xt, parMargin, typeMargin)
autocor.emp |
target correlation between intensities |
nChainFit |
number of simulations |
Xt |
simulated occurrences, nChainFit x 2 matrix |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
scalar |
needed correlation |
Guillaume Evin
finds the correlation between normal variates leading to correlation between occurrences
find.omega(rho.emp, Qtrans.mat, mat.comb, nLag, nChainFit)find.omega(rho.emp, Qtrans.mat, mat.comb, nLag, nChainFit)
rho.emp |
target correlation between occurences |
Qtrans.mat |
transition probabilities, 2 x ncomb matrix |
mat.comb |
matrix of logical: ncomb x nlag |
nLag |
order of the Markov chain |
nChainFit |
length of the simulated chains used during the fitting |
scalar |
needed correlation |
Guillaume Evin
finds the correlation between normal variates leading to correlation between intensities
find.zeta(eta.emp, nChainFit, Xt, parMargin, typeMargin)find.zeta(eta.emp, nChainFit, Xt, parMargin, typeMargin)
eta.emp |
target correlation between intensities |
nChainFit |
number of simulations |
Xt |
simulated occurrences, n x 2 matrix |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
scalar |
needed correlation |
Guillaume Evin
estimate parameters which control the spatial dependence between intensities using a copula
fit.copula.amount(P.mat, isPeriod, th, copulaInt, M0)fit.copula.amount(P.mat, isPeriod, th, copulaInt, M0)
P.mat |
precipitation matrix |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
copulaInt |
type of dependence between inter-site amounts: 'Gaussian' or 'Student' |
M0 |
covariance matrix of gaussianized prec. amounts for all pairs of stations |
list |
list of estimates (e.g., M0, dfStudent) |
Guillaume Evin
estimate all the parameters for the G-Wex model of precipitation
fit.GWex.prec(objGwexObs, parMargin, listOption = NULL)fit.GWex.prec(objGwexObs, parMargin, listOption = NULL)
objGwexObs |
object of class |
parMargin |
if not NULL, list where each element parMargin[[iM]] corresponds to a month iM=1...12 and contains a matrix nStation x 3 of estimated parameters of the marginal distributions (EGPD or mixture of exponentials) |
listOption |
list with the following fields:
|
a list containing the list of options listOption and the list of estimated parameters listPar.
The parameters of the occurrence process are contained in parOcc and the parameters related to the precipitation
amounts are contained in parInt. Each type of parameter is a list containing the estimates for each month. In parOcc, we find:
p01: For each station, the probability of transition from a dry state to a wet state.
p11: For each station, the probability of staying in a wet state.
list.pr.state: For each station, the probabilities of transitions for a Markov chain with lag p.
list.mat.omega: The spatial correlation matrix of occurrences (see Evin et al., 2018).
In parInt, we have:
parMargin: list of matrices nStation x nPar of parameters for the marginal distributions (one element per Class).
cor.int: Matrices nStation x nStation , , representing
the spatial and temporal correlations between all the stations (see Evin et al., 2018). For the
Student copula, dfStudent indicates the parameter.
Guillaume Evin
Evin, G., A.-C. Favre, and B. Hingray. 2018. 'Stochastic Generation of Multi-Site Daily Precipitation Focusing on Extreme Events.' Hydrol. Earth Syst. Sci. 22 (1): 655-672. doi.org/10.5194/hess-22-655-2018.
estimate parameters which control the dependence between intensities with a MAR(1) process
fit.MAR1.amount(P.mat, isPeriod, th, copulaInt, M0, A)fit.MAR1.amount(P.mat, isPeriod, th, copulaInt, M0, A)
P.mat |
precipitation matrix |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
copulaInt |
type of dependance between inter-site amounts: 'Gaussian' or 'Student' |
M0 |
covariance matrix of gaussianized prec. amounts for all pairs of stations |
A |
Matrix containing the autocorrelation (temporal) correlations |
list with the following items
M0 covariance matrix of gaussianized prec. amounts for all pairs of stations
A omega correlations for all pairs of stations
covZ covariance matrix of the MAR(1) process
sdZ standard deviation of the diagonal elements
corZ correlation matrix of the MAR(1) process
dfStudent degrees of freedom for the Student copula if CopulaInt is equal to "Student"
Guillaume Evin
Matalas, N. C. 1967. “Mathematical Assessment of Synthetic Hydrology.” Water Resources Research 3 (4): 937–45. https://doi.org/10.1029/WR003i004p00937.
Bárdossy, A., and G. G. S. Pegram. 2009. “Copula Based Multisite Model for Daily Precipitation Simulation.” Hydrology and Earth System Sciences 13 (12): 2299–2314. https://doi.org/10.5194/hess-13-2299-2009.
estimate parameters which control the marginal distribution of precipitation amounts
fit.margin.cdf(P.mat, isPeriod, th, type = c("EGPD", "mixExp"))fit.margin.cdf(P.mat, isPeriod, th, type = c("EGPD", "mixExp"))
P.mat |
precipitation matrix |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
type |
distribution: 'EGPD' or 'mixExp' |
matrix |
matrix of estimates p x 3 |
Guillaume Evin
fitGwexModel: fit a GWex model to observations.
fitGwexModel(objGwexObs, parMargin = NULL, listOption = NULL)fitGwexModel(objGwexObs, parMargin = NULL, listOption = NULL)
objGwexObs |
an object of class |
parMargin |
(not required for temperature) list parMargin where and each element corresponds to a month (1...12) and contains a matrix nStation x 3 of pre-estimated parameters of the marginal distributions (EGPD or Mixture of Exponentials) |
listOption |
for precipitation, a list with the following fields:
and for temperature, a list with the following fields:
|
Return an object of class GwexFit with:
p: The number of station,
version: package version,
variable: the type of variable,
fit: a list containing the list of options listOption
and the list of estimated parameters listPar.
Guillaume Evin
# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') ############################################################### # FIT THE PRECIPITATION MODEL ############################################################### # Format observations: create a Gwex object for one station only to show a quick # example. The syntax is similar for multi-site applications. myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX[,1,drop=FALSE]) # Fit precipitation model with a threshold of 0.5 mm to distinguish wet and dry # states (th) and keep default options otherwise, e.g. a Gaussian # copula for the spatial dependence (copulaInt) and a mixExp distribution for # marginal intensities ('typeMargin') myParPrec = fitGwexModel(myObsPrec,listOption=list(th=0.5)) myParPrec # print object ############################################################### # FIT THE TEMPERATURE MODEL, COND. TO PRECIPITATION ############################################################### # Format observations: create a G-Wex object myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # Fit temperature model with a long-term linear trend ('hasTrend'), Gaussian margins # ('typeMargin') and Gaussian spatial dependence ('depStation') myParTemp = fitGwexModel(myObsTemp,listOption=list(hasTrend=TRUE,typeMargin='Gaussian', depStation='Gaussian')) myParTemp # print object# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') ############################################################### # FIT THE PRECIPITATION MODEL ############################################################### # Format observations: create a Gwex object for one station only to show a quick # example. The syntax is similar for multi-site applications. myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX[,1,drop=FALSE]) # Fit precipitation model with a threshold of 0.5 mm to distinguish wet and dry # states (th) and keep default options otherwise, e.g. a Gaussian # copula for the spatial dependence (copulaInt) and a mixExp distribution for # marginal intensities ('typeMargin') myParPrec = fitGwexModel(myObsPrec,listOption=list(th=0.5)) myParPrec # print object ############################################################### # FIT THE TEMPERATURE MODEL, COND. TO PRECIPITATION ############################################################### # Format observations: create a G-Wex object myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # Fit temperature model with a long-term linear trend ('hasTrend'), Gaussian margins # ('typeMargin') and Gaussian spatial dependence ('depStation') myParTemp = fitGwexModel(myObsTemp,listOption=list(hasTrend=TRUE,typeMargin='Gaussian', depStation='Gaussian')) myParTemp # print object
First parametric family for G(v) = v^kappa: distribution, density and quantile function
EGPD.pGI(v, kappa) EGPD.dGI(v, kappa) EGPD.qGI(p, kappa)EGPD.pGI(v, kappa) EGPD.dGI(v, kappa) EGPD.qGI(p, kappa)
v |
probability |
kappa |
transformation parameter greater than 0 |
p |
probability |
distribution, density and quantile of EGPD
Guillaume Evin
Estimates the nu parameter (degrees of freedom) of the multivariate Student distribution when the correlation matrix Sig is given
get.df.Student(P, Sig, max.df = 20)get.df.Student(P, Sig, max.df = 20)
P |
matrix of non-zero precipitation (zero precipitation are set to NA) |
Sig |
correlation matrix |
max.df |
maximum degrees of freedom tested (default=20) |
nu estimate
Guillaume Evin
McNeil et al. (2005) "Quantitative Risk Management"
get the cdf values (empirical distribution) of positive precipitation
get.emp.cdf.matrix(X)get.emp.cdf.matrix(X)
X |
matrix of positive precipitation |
matrix with cdf values (NA if zero precipitation)
Guillaume Evin
return a vector of 3-char tags of the 12 months
get.list.month()get.list.month()
get the vector of the four seasons c('DJF','MAM','JJA','SON')
get.list.season()get.list.season()
Guillaume Evin
get default options and check values proposed by the user
get.listOption(listOption)get.listOption(listOption)
listOption |
list containing fields corr. to the different options. Can be NULL if no options are set |
listOption |
list of options |
Guillaume Evin
find matrix of correlations leading to estimates cor between intensities
get.M0( cor.obs, infer.mat.omega.out, nLag, parMargin, typeMargin, nChainFit, isParallel )get.M0( cor.obs, infer.mat.omega.out, nLag, parMargin, typeMargin, nChainFit, isParallel )
cor.obs |
matrix p x p of observed correlations between intensities for all pairs of stations |
infer.mat.omega.out |
output of |
nLag |
order of the Markov chain |
parMargin |
parameters of the margins p x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
nChainFit |
integer indicating the length of simulated chains |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
list with two items
Xt long simulation of the wet/dry states according to the model
M0 covariance matrix of gaussianized prec. amounts for all pairs of stations
Guillaume Evin
find omega correlation leading to estimates cor between occurrences
get.mat.omega(cor.obs, Qtrans.mat, mat.comb, nLag, nChainFit, isParallel)get.mat.omega(cor.obs, Qtrans.mat, mat.comb, nLag, nChainFit, isParallel)
cor.obs |
matrix p x p of observed correlations between occurrences for all pairs of stations |
Qtrans.mat |
transition probabilities, 2 x ncomb matrix |
mat.comb |
matrix of logical: ncomb x nlag |
nLag |
order of the Markov chain |
nChainFit |
length of the simulated chains used during the fitting |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
matrix |
omega correlations for all pairs of stations |
Guillaume Evin
get.period.fitting.month
get.period.fitting.month(m.char)get.period.fitting.month(m.char)
m.char |
3-letter name of a month (e.g. 'JAN') return the 3 indices corresponding to the 3-month period of a month ('JAN') |
find rho autocorrelation leading to empirical estimates
get.vec.autocor(vec.ar1.obs, Xt, parMargin, typeMargin, nChainFit, isParallel)get.vec.autocor(vec.ar1.obs, Xt, parMargin, typeMargin, nChainFit, isParallel)
vec.ar1.obs |
vector of observed autocorrelations for all stations |
Xt |
simulated occurrences given model parameters of wet/dry states |
parMargin |
parameters of the margins p x 3 |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
nChainFit |
integer indicating the length of the simulated chains |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
vector |
vector of rho parameters to simulate the MAR process |
Guillaume Evin
get object GwexFit derived from the parameters replicated for each month
getGwexFitPrec( listOption = NULL, p, condProbaWDstates, parMargin, vec.ar1 = NULL, M0 = NULL, mat.omega = NULL )getGwexFitPrec( listOption = NULL, p, condProbaWDstates, parMargin, vec.ar1 = NULL, M0 = NULL, mat.omega = NULL )
listOption |
list of options (see |
p |
number of stations |
condProbaWDstates |
vector of length nLag^2 of transition probabilities
corresponding to the nlag possible transitions between dry/wet states
|
parMargin |
parameters of the margins: vector of length 3 |
vec.ar1 |
vector of observed autocorrelations for all stations |
M0 |
M0: covariance matrix of gaussianized prec. amounts for all pairs of stations |
mat.omega |
mat.omega: The spatial correlation matrix of occurrences |
Return an object of class GwexFit with:
p: The number of station,
version: package version,
variable: the type of variable,
fit: a list containing the list of options listOption
and the list of estimated parameters listPar.
exFitGwexPrec = getGwexFitPrec(p=2,condProbaWDstates=c(0.7,0.3,0.2,0.1), parMargin=c(0.5,0.1,0.4),vec.ar1=rep(0.7,2),M0=rbind(c(1,0.6),c(0.6,1)), mat.omega=rbind(c(1,0.8),c(0.8,1)))exFitGwexPrec = getGwexFitPrec(p=2,condProbaWDstates=c(0.7,0.3,0.2,0.1), parMargin=c(0.5,0.1,0.4),vec.ar1=rep(0.7,2),M0=rbind(c(1,0.6),c(0.6,1)), mat.omega=rbind(c(1,0.8),c(0.8,1)))
Defines a generic Gwex object.
GWex objects contain two slots:
- the version ('vX.X.X')
- the type of variable ('Prec' or 'Temp')
Guillaume Evin
Defines a GwexFit object which is a Gwex object containing 'fit', a list containing the fitted parameters, and 'p', the number of stations.
See fitGwexModel for some examples.
Guillaume Evin
Constructor of class [GwexObs]
GwexObs(variable, date, obs)GwexObs(variable, date, obs)
variable |
'Prec' or 'Temp' |
date |
vector of class 'Date' |
obs |
matrix nTime x nStations of observations |
An object of class [GwexObs]
# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with precipitation data myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX) # print GwexObs object myObsPrec # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with precipitation data myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX) # print GwexObs object myObsPrec # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp
GwexObs
Defines a GwexObs object which is a Gwex object containing dates and a matrix of observations.
Guillaume Evin
# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with precipitation data myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX) # print GwexObs object myObsPrec # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with precipitation data myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX) # print GwexObs object myObsPrec # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp
GwexSim object which is a Gwex object containing 'sim', an array containing the simulations, and 'dates', a vector of dates.
See simGwexModel for some examples.Defines a GwexSim object which is a Gwex object containing 'sim', an array containing the simulations, and 'dates', a vector of dates.
See simGwexModel for some examples.
Guillaume Evin
special case of infer.dep.amount where there is only one station
infer.autocor.amount( P.mat, pr.state, isPeriod, nLag, th, parMargin, typeMargin, nChainFit, isMAR, isParallel )infer.autocor.amount( P.mat, pr.state, isPeriod, nLag, th, parMargin, typeMargin, nChainFit, isMAR, isParallel )
P.mat |
precipitation matrix |
pr.state |
probabilities of transitions for a Markov chain with lag |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
nLag |
order of he Markov chain for the transitions between dry and wet states (=2 by default) |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
'EGPD' (Extended GPD) or 'mixExp' (Mixture of Exponentials). 'EGPD' by default |
nChainFit |
integer, length of the runs used during the fitting procedure. =100000 by default |
isMAR |
logical value, do we apply a Autoregressive Multivariate Autoregressive model (order 1) =TRUE by default |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
list |
list of estimates (e.g., M0, dfStudent) |
Guillaume Evin
estimate parameters which control the spatial dependence between intensities using a copula
infer.dep.amount( P.mat, isPeriod, infer.mat.omega.out, nLag, th, parMargin, typeMargin, nChainFit, isMAR, copulaInt, isParallel )infer.dep.amount( P.mat, isPeriod, infer.mat.omega.out, nLag, th, parMargin, typeMargin, nChainFit, isMAR, copulaInt, isParallel )
P.mat |
precipitation matrix |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
infer.mat.omega.out |
output of |
nLag |
order of he Markov chain for the transitions between dry and wet states (=2 by default) |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
parMargin |
parameters of the margins 2 x 3 |
typeMargin |
'EGPD' (Extended GPD) or 'mixExp' (Mixture of Exponentials). 'EGPD' by default |
nChainFit |
integer, length of the runs used during the fitting procedure. =100000 by default |
isMAR |
logical value, do we apply a Autoregressive Multivariate Autoregressive model (order 1) =TRUE by default |
copulaInt |
'Gaussian' or 'Student': type of dependence for amounts (='Student' by default) |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
list |
list of estimates (e.g., M0, dfStudent) |
Guillaume Evin
find omega correlation leading to estimates cor between occurrences
infer.mat.omega(P.mat, isPeriod, th, nLag, pr.state, nChainFit, isParallel)infer.mat.omega(P.mat, isPeriod, th, nLag, pr.state, nChainFit, isParallel)
P.mat |
matrix of precipitation n x p |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
nLag |
order of the Markov chain |
pr.state |
output of function |
nChainFit |
length of the simulated chains used during the fitting |
isParallel |
logical: indicate computation in parallel or not (easier for debugging) |
A list with different objects
Qtrans.mat: matrix nStation x n.comb of transition probabilites
mat.comb: matrix of possible combination n.comb x nLag
mat.omega: The spatial correlation matrix of occurrences (see Evin et al., 2018).
Guillaume Evin
joint probabilities of occurrences for all pairs of stations
joint.proba.occ(P, th)joint.proba.occ(P, th)
P |
matrix of precipitation |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
list |
list of joint probabilities |
Guillaume Evin
Estimate the transition probabilities between wet and dry states, for nlag previous days, for all stations
lagTransProbaMatrix(mat.prec, isPeriod, th, nlag)lagTransProbaMatrix(mat.prec, isPeriod, th, nlag)
mat.prec |
matrix of precipitation |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
nlag |
number of lag days |
list |
list with one item per station, where each item is a matrix nLag^2 x (nLag+1) of transition probability between dry/wet state. The first nLag columns indicate the wet/dry states for the previous nLag days |
Guillaume Evin
Estimate the transition probabilities between wet and dry states, for nlag previous days, for one station
lagTransProbaVector(vec.prec, isPeriod, th, nlag)lagTransProbaVector(vec.prec, isPeriod, th, nlag)
vec.prec |
vector nx1 of precipitation for one station |
isPeriod |
vector of logical n x 1 indicating the days concerned by a 3-month period |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
nlag |
number of lag days |
matrix |
matrix nLag^2 x (nLag+1) of transition probability between dry/wet state. The first nLag columns indicate the wet/dry states for the previous nLag days |
Guillaume Evin
Mask intensities where there is no occurrence
mask.GWex.Yt(Xt, Yt)mask.GWex.Yt(Xt, Yt)
Xt |
simulated occurrences |
Yt |
simulated intensities |
matrix |
matrix n x p of simulated precipitations |
Guillaume Evin
Modify a non-positive definite correlation matrix in order to have a positive definite matrix
modify.cor.matrix(cor.matrix)modify.cor.matrix(cor.matrix)
cor.matrix |
possibly non-positive definite correlation matrix |
positive definite correlation matrix
Guillaume Evin
Rousseeuw, P. J. and G. Molenberghs. 1993. Transformation of non positive semidefinite correlation matrices. Communications in Statistics: Theory and Methods 22(4):965-984.
Rebonato, R., & Jackel, P. (2000). The most general methodology to create a valid correlation matrix for risk management and option pricing purposes. J. Risk, 2(2), 17-26.
transform vector of months to seasons
month2season(vecMonth)month2season(vecMonth)
vecMonth |
a vector of months given as integers 1:12 |
Guillaume Evin
print-methods: Create a method to print Gwex objects.
## S4 method for signature 'Gwex' print(x) ## S4 method for signature 'GwexObs' print(x) ## S4 method for signature 'GwexFit' print(x) ## S4 method for signature 'GwexSim' print(x)## S4 method for signature 'Gwex' print(x) ## S4 method for signature 'GwexObs' print(x) ## S4 method for signature 'GwexFit' print(x) ## S4 method for signature 'GwexSim' print(x)
x |
|
# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # print GwexObs object myObsTemp
Probability Weighted Moments of order 0, 1 and 2 of the unified EGPD distribution
EGPD.GI.mu0(kappa, sig, xi) EGPD.GI.mu1(kappa, sig, xi) EGPD.GI.mu2(kappa, sig, xi)EGPD.GI.mu0(kappa, sig, xi) EGPD.GI.mu1(kappa, sig, xi) EGPD.GI.mu2(kappa, sig, xi)
kappa |
transformation parameter greater than 0 |
sig |
Scale parameter |
xi |
Shape parameter |
Probability Weighted Moments
Guillaume Evin
reshape Qtrans.mat to an array
QtransMat2Array(n, p, Qtrans.mat)QtransMat2Array(n, p, Qtrans.mat)
n |
matrix of precipitation |
p |
number of stations |
Qtrans.mat |
transition probabilities, 2 x ncomb matrix |
array |
array of transition probabilities with dimension n x p x n.comb |
Guillaume Evin
show-methods: Create a method to show Gwex objects.
## S4 method for signature 'Gwex' show(object) ## S4 method for signature 'GwexObs' show(object) ## S4 method for signature 'GwexFit' show(object) ## S4 method for signature 'GwexSim' show(object)## S4 method for signature 'Gwex' show(object) ## S4 method for signature 'GwexObs' show(object) ## S4 method for signature 'GwexFit' show(object) ## S4 method for signature 'GwexSim' show(object)
object |
|
# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # show GwexObs object myObsTemp# Format dates corresponding to daily observations of precipitation and temperature vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') # build GwexObs object with temperature data myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # show GwexObs object myObsTemp
generate boolean variates which describe the dependence between intersite occurrence correlations and wet/dry persistence
sim.GWex.occ(objGwexFit, vecMonth)sim.GWex.occ(objGwexFit, vecMonth)
objGwexFit |
object of class GwexFit |
vecMonth |
vector n x 1 of integers indicating the months |
matrix of logical |
occurrences simulated |
Guillaume Evin
Simulate one scenario of precipitation from the GWex model
sim.GWex.prec.1it(objGwexFit, vecDates, myseed, objGwexObs, prob.class)sim.GWex.prec.1it(objGwexFit, vecDates, myseed, objGwexObs, prob.class)
objGwexFit |
object of class GwexFit |
vecDates |
vector of continuous dates |
myseed |
seed of the random generation, to be fixed if the results need to be replicated |
objGwexObs |
optional: necessary if we need observations to simulate (e.g. disaggregation of 3-day periods) |
prob.class |
vector of probabilities indicating class of "similar" mean intensities |
matrix |
Precipitation simulated for the dates contained in vec.Dates at the different stations |
Guillaume Evin
Inverse PIT: from the probability space to the precipitation space
sim.GWex.Yt(objGwexFit, vecMonth, Yt.Pr)sim.GWex.Yt(objGwexFit, vecMonth, Yt.Pr)
objGwexFit |
object of class GwexFit |
vecMonth |
vector of integer indicating the months |
Yt.Pr |
uniform variates describing dependence between inter-site amounts |
matrix |
matrix n x p of simulated non-zero precipitation intensities |
Guillaume Evin
generate uniform variates which describe the dependence between intersite amount correlations
sim.GWex.Yt.Pr(objGwexFit, vecMonth)sim.GWex.Yt.Pr(objGwexFit, vecMonth)
objGwexFit |
object of class GwexFit |
vecMonth |
vector n x 1 of integer indicating the months |
matrix |
matrix n x p of uniform dependent variates |
Guillaume Evin
get relevant parameters
sim.GWex.Yt.Pr.get.param(objGwexFit, iM)sim.GWex.Yt.Pr.get.param(objGwexFit, iM)
objGwexFit |
object of class GwexFit |
iM |
integer indicating the month |
list |
list of parameters |
Guillaume Evin
generate gaussian variates which describe the spatial and temporal dependence between the sites (MAR(1) process)
sim.Zt.MAR(PAR, copulaInt, Zprev, p)sim.Zt.MAR(PAR, copulaInt, Zprev, p)
PAR |
parameters for this class |
copulaInt |
'Gaussian' or 'Student' |
Zprev |
previous Gaussian variate |
p |
number of stations |
matrix |
matrix n x p of uniform dependent variates |
Guillaume Evin
generate gaussian variates which describe the spatial dependence between the sites
sim.Zt.Spatial(PAR, copulaInt, p)sim.Zt.Spatial(PAR, copulaInt, p)
PAR |
parameters for a class |
copulaInt |
'Gaussian' or 'Student' |
p |
number of stations |
matrix |
matrix n x p of uniform dependent variates |
Guillaume Evin
Simulate from a GWex model
simGwexModel( objGwexFit, nb.rep = 10, d.start = as.Date("01011900", "%d%m%Y"), d.end = as.Date("31121999", "%d%m%Y"), objGwexObs = NULL, prob.class = c(0.5, 0.75, 0.9, 0.99), objGwexSim = NULL, nCluster = 1, useSeed = FALSE )simGwexModel( objGwexFit, nb.rep = 10, d.start = as.Date("01011900", "%d%m%Y"), d.end = as.Date("31121999", "%d%m%Y"), objGwexObs = NULL, prob.class = c(0.5, 0.75, 0.9, 0.99), objGwexSim = NULL, nCluster = 1, useSeed = FALSE )
objGwexFit |
an object of class |
nb.rep |
number of repetitions of scenarios |
d.start |
a starting date for the simulation |
d.end |
an ending date for the simulation |
objGwexObs |
optional: an object of class |
prob.class |
vector of probabilities indicating class of "similar" mean intensities |
objGwexSim |
optional: an object of class |
nCluster |
optional, number of clusters which can be used for the parallel computation |
useSeed |
optional: a logical variable indicating if we control the seed of the simulations or not |
GwexSim |
an object of class |
Guillaume Evin
# vector of dates vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') ############################################################### # FIT AND SIMULATE FROM THE PRECIPITATION MODEL ############################################################### # Format observations: create a G-Wex object myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX[,1,drop=FALSE]) # Fit GWEX precipitation model, default options except for the threshold th myParPrec = fitGwexModel(myObsPrec,listOption=list(th=0.5)) # fit model # Generate 2 scenarios for one year, using the 'GwexFit' object mySimPrec = simGwexModel(objGwexFit=myParPrec, nb.rep=2, d.start=vecDates[1], d.end=vecDates[10]) mySimPrec # print object ############################################################### # FIT AND SIMULATE FROM THE TEMPERATURE MODEL ############################################################### # Format observations: create a G-Wex object myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # Fit GWEX temperature model myParTemp = fitGwexModel(myObsTemp,listOption=list(hasTrend=TRUE,typeMargin='Gaussian', depStation='Gaussian')) # Generate 2 scenarios for one year, using an existing 'GwexFit' object mySimTemp = simGwexModel(objGwexFit=myParTemp, nb.rep=2, d.start=vecDates[1], d.end=vecDates[365],objGwexObs=myObsPrec) mySimTemp # print object# vector of dates vecDates = seq(from=as.Date("01/01/2005",format="%d/%m/%Y"), to=as.Date("31/12/2014",format="%d/%m/%Y"),by='day') ############################################################### # FIT AND SIMULATE FROM THE PRECIPITATION MODEL ############################################################### # Format observations: create a G-Wex object myObsPrec = GwexObs(variable='Prec',date=vecDates,obs=dailyPrecipGWEX[,1,drop=FALSE]) # Fit GWEX precipitation model, default options except for the threshold th myParPrec = fitGwexModel(myObsPrec,listOption=list(th=0.5)) # fit model # Generate 2 scenarios for one year, using the 'GwexFit' object mySimPrec = simGwexModel(objGwexFit=myParPrec, nb.rep=2, d.start=vecDates[1], d.end=vecDates[10]) mySimPrec # print object ############################################################### # FIT AND SIMULATE FROM THE TEMPERATURE MODEL ############################################################### # Format observations: create a G-Wex object myObsTemp = GwexObs(variable='Temp',date=vecDates,obs=dailyTemperGWEX) # Fit GWEX temperature model myParTemp = fitGwexModel(myObsTemp,listOption=list(hasTrend=TRUE,typeMargin='Gaussian', depStation='Gaussian')) # Generate 2 scenarios for one year, using an existing 'GwexFit' object mySimTemp = simGwexModel(objGwexFit=myParTemp, nb.rep=2, d.start=vecDates[1], d.end=vecDates[365],objGwexObs=myObsPrec) mySimTemp # print object
find matrix of correlations leading to estimates cor between intensities
simPrecipOcc(nLag, n, pr)simPrecipOcc(nLag, n, pr)
nLag |
order of the Markov chain |
n |
integer indicating the length of simulated chains |
pr |
vector of probabilies corr. to the conditional transition probabilities |
a vector Xt of length n with values 0/1 corr. to dry/wet states
Guillaume Evin
from uniform variates to precipitation variates
unif.to.prec(pI, typeMargin, U)unif.to.prec(pI, typeMargin, U)
pI |
vector of three parameters of the marginal distributions |
typeMargin |
type of marginal distribution: 'EGPD' or 'mixExp' |
U |
vector of uniform variates |
matrix |
matrix of estimates p x 3 |
Guillaume Evin
Estimate the wet day frequency (proportion of wet days) for all stations
wet.day.frequency(mat.prec, th)wet.day.frequency(mat.prec, th)
mat.prec |
matrix of precipitation (possibly for one month/period) |
th |
threshold above which we consider that a day is wet (e.g. 0.2 mm) |
vector of numeric |
wet day frequencies |
Guillaume Evin