R Under development (unstable) (2023-09-17 r85158 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # some simple testing commands...
> #   testthat would be an overshoot...
> 
> require(PRSim)
Loading required package: PRSim
> 
> # raw demos
> demo( "PRSim", ask=FALSE)


	demo(PRSim)
	---- ~~~~~

> # short demo
> data(runoff)

> out <- prsim(data=runoff, number_sim=1, marginal="empirical")
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.

> out <- prsim(data=runoff, number_sim=1, marginal="kappa", GoFtest = "KS")
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.

> ### GEV distribution
> require("evd")
Loading required package: evd

> require("ismev")
Loading required package: ismev
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.9-0. For overview type 'help("mgcv-package")'.

> 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, show=FALSE, ...)$mle

> ### GEV
> out <- prsim(data=runoff, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.

> sim <- out$simulation

> # p_val <- out$p_val
> # par(mai=c(.9,.9,.1,.1))
> plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l", 
+      xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))

> matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
+          lty=1, col="gray")
> demo( "PRSim-validate", ask=FALSE)


	demo(PRSim-validate)
	---- ~~~~~~~~~~~~~~

> # short demo
> data(simulations)

> sim <- simulations$simulation

> # periodogram of deseasonalized
> kern <- kernel("modified.daniell",c(10,10))

> sp1 <- spec.pgram(sim$Qobs, k=kern, taper=0, log="no", plot=FALSE)

> sp2 <- spec.pgram(sim$des, k=kern, taper=0, log="no", plot=FALSE)

> plot(sp1, xlim=c(0,.05))     

> plot( sp2, add=TRUE, col=2)

> # Peaks correspond to the following cycles:
> 1/sp1$freq[head(order(sp1$spec, decreasing=TRUE))]
[1] 355.2632 375.0000 337.5000 321.4286 397.0588 306.8182

> # compare periodogram of simulated series
> plot(sp1, xlim=c(0,.05))     # would be nice to identify the peaks...

> for (i in grep("r",names(sim))) {
+   spi <- spec.pgram(sim[,i], k=kern, taper=0, log="no", plot=FALSE)
+   plot( spi, add=TRUE, col="gray")
+ }

> sp3 <- spec.pgram(sim$Qobs, taper=0, log="no", plot=FALSE)

> 1/sp3$freq[head(order(sp3$spec, decreasing=TRUE))]
[1] 375.0000 355.2632 182.4324 122.7273 337.5000 421.8750

> # Annual, 6 months and 4 months
> 
> 
> ### plot mean regime for each simulation run and compare to observed regime
> ### define plotting colors
> col_sim <- adjustcolor("#fd8d3c",alpha=0.8)

> col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)

> col_obs <- adjustcolor( "black", alpha.f = 0.2)

> year <- unique(sim$YYYY)

> ### compute mean runoff hydrograph
> sim$day_id <- rep(seq(1:365),times=length(year))

> mean_hydrograph_obs <- aggregate(sim$Qobs, by=list(sim$day_id), FUN=mean,simplify=FALSE)

> plot(unlist(mean_hydrograph_obs[,2]), lty=1, lwd=1, col="black", ylab=expression(paste("Discharge [m"^3,"/s]")),
+      xlab="Time [d]", main="Mean hydrographs", ylim=c(0,max(unlist(mean_hydrograph_obs[,2]))*1.5),type="l")

> ### add mean runoff hydrographs
> for(r in 7:(length(names(sim))-1)){
+   mean_hydrograph <- aggregate(sim[,r], by=list(sim$day_id), FUN=mean,simplify=FALSE)
+   lines(mean_hydrograph, lty=1, lwd=1, col=col_sim)
+ }

> ### redo observed mean
> lines(mean_hydrograph_obs, lty=1, lwd=1, col="black")

> ### autocorrelation
> acf_mare <- list()

> acf_obs <- acf(sim$Qobs, plot=FALSE)

> plot(acf_obs$acf, type="l", xlab="Lag", main="Autocorrelation", ylab="ACF")

> for(r in 7:(length(names(sim))-2)){
+   acf_sim <- acf(sim[,r], plot=FALSE)
+   lines(acf_sim$acf, col=col_sim, type="l")
+   ### compute mean relative error in the acf
+   acf_mare[[r]]<- mean(abs((acf_obs$acf-acf_sim$acf)/acf_obs$acf))
+ }

> lines(acf_obs$acf)

> ### partial autocorrelation function
> pacf_obs <- pacf(sim$Qobs, plot=FALSE)

> pacf_mare <- list()

> plot(pacf_obs$acf, type="l", xlab="Lag", main="Partial autocorrelation", ylab="PACF")

> for(r in 7:(length(names(sim))-2)){
+   pacf_sim <- pacf(sim[,r], plot=FALSE)
+   lines(pacf_sim$acf, col=col_sim, type="l")
+   ### compute mean relative error in the acf
+   pacf_mare[[r]] <- mean(abs((pacf_obs$acf-pacf_sim$acf)/pacf_obs$acf))
+ }

> lines(pacf_obs$acf)

> ### compute seasonal statistics
> ### Q50,Q05,Q95, boxplots
> ### define seasons: Winter:12,1,2; spring:3,4,5; summer: 6,7,8; fall: 9,10,11
> sim$season <- "winter"

> sim$season[which(sim$MM%in%c(3,4,5))] <- "spring"

> sim$season[which(sim$MM%in%c(6,7,8))] <- "summer"

> sim$season[which(sim$MM%in%c(9,10,11))] <- "fall"

> ### all simulated series show the same seasonal statistics. plot only one 
> boxplot(sim$Qobs[which(sim$season=="winter")], sim$r1[which(sim$season=="winter")],
+         sim$Qobs[which(sim$season=="spring")], sim$r1[which(sim$season=="spring")],
+         sim$Qobs[which(sim$season=="summer")], sim$r1[which(sim$season=="summer")],
+         sim$Qobs[which(sim$season=="fall")], sim$r1[which(sim$season=="fall")],
+         border=c("black", col_sim, "black", col_sim, "black", col_sim, "black", col_sim), 
+         xaxt="n", main="Seasonal statistics", outline=FALSE)

> mtext(side=1, text=c("Winter", "Spring", "Summer", "Fall"), at=c(1.5,3.5,5.5,7.5))
> 
> 
> 
> # testing input
> data(runoff)
> unique(runoff$YYYY)
 [1] 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989
[16] 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
[31] 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
> 
> 
> try( prsim( runoff[1:130,] ))   #  At least one year of data required.
Error in prsim(runoff[1:130, ]) : At least one year of data required.
> try( prsim( runoff[1:730,] ))   #  No missing values allowed. Some days are missing.
Error in prsim(runoff[1:730, ]) : 
  No missing values allowed. Some days are missing.
> try( prsim( runoff[1:1445,] ))  #  No missing values allowed. Some days are missing.
Error in prsim(runoff[1:1445, ]) : 
  No missing values allowed. Some days are missing.
> try( prsim( runoff[runoff$YYYY<1976,] )) # At least one year of data required.
Error in prsim(runoff[runoff$YYYY < 1976, ]) : 
  At least one year of data required.
> 
> 
> 
> suppressWarnings( out <- prsim( runoff[runoff$YYYY<1977,] ) )
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> 
> runof <- runoff[runoff$YYYY<1980,]
> 
> set.seed(1)
> str(out1 <- prsim( runof, marginalpar=FALSE, suppWarn=TRUE))
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
List of 3
 $ simulation:'data.frame':	1825 obs. of  7 variables:
  ..$ YYYY         : int [1:1825] 1975 1975 1975 1975 1975 1975 1975 1975 1975 1975 ...
  ..$ MM           : int [1:1825] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ DD           : int [1:1825] 1 2 3 4 5 6 7 8 9 10 ...
  ..$ timestamp    : POSIXct[1:1825], format: "1975-01-01" "1975-01-02" ...
  ..$ Qobs         : num [1:1825] 2.05 1.75 1.62 1.58 1.47 ...
  ..$ deseaonalized: num [1:1825] 1.595 0.738 1.512 0.944 0.919 ...
  ..$ r1           : num [1:1825] 0.812 0.665 0.449 0.842 0.928 ...
 $ pars      : NULL
 $ p_val     : NULL
> 
> runo <- runof
> names( runo) <- tolower( names(runof))
> try( prsim( runo, marginalpar=FALSE, suppWarn=TRUE)) #  Wrong column for observations selected.
Error in prsim(runo, marginalpar = FALSE, suppWarn = TRUE) : 
  Wrong column (name) for observations selected.
> 
> runo <- runof[,4:1]
> set.seed(1)
> out3 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) #  ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> identical(out1,out3) 
[1] TRUE
> 
> runo <- runof[,4:1]
> set.seed(1)
> out4 <- prsim( runo, station_id=1, marginalpar=FALSE, suppWarn=TRUE) #  ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> identical(out1,out4) 
[1] TRUE
> 
> tmp <- paste(runof$YYYY, runof$MM, runof$DD,sep=" ")
> runo <- data.frame(time=as.POSIXct(strptime(tmp, format="%Y %m %d", tz="GMT")), Qobs=runof$Qobs)
> set.seed(1)
> out5 <- prsim( runo, marginalpar=FALSE, suppWarn=TRUE) #  ok
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> identical(out1,out5) 
[1] TRUE
> 
> # 
> 
> ######################
> # Test 'kappa' distribution with manual construction:
> rKappa <- function(n, theta) rand.kappa(n, theta[1], theta[2], theta[3], theta[4])
> Kappa_fit <- function(xdat, ...) {
+   ll <- Lmoments(xdat)  
+   unlist(par.kappa(ll[1],ll[2],ll[4],ll[5]))
+ }
> set.seed(1)
> out6a <- prsim( runo, marginalpar=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> set.seed(1)
> out6b <- prsim( runo, marginal="kappa", marginalpar=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> identical(out6a$pars, out6b$pars)   # columns are differently named...
[1] TRUE
> colSums( (as.matrix(out6a$pars)-out6b$pars)^2)
  xi alfa    k    h 
   0    0    0    0 
> summary(out6a$simulation-out6b$simulation)
      YYYY         MM          DD     timestamp             Qobs  
 Min.   :0   Min.   :0   Min.   :0   Length:1825       Min.   :0  
 1st Qu.:0   1st Qu.:0   1st Qu.:0   Class :difftime   1st Qu.:0  
 Median :0   Median :0   Median :0   Mode  :numeric    Median :0  
 Mean   :0   Mean   :0   Mean   :0                     Mean   :0  
 3rd Qu.:0   3rd Qu.:0   3rd Qu.:0                     3rd Qu.:0  
 Max.   :0   Max.   :0   Max.   :0                     Max.   :0  
 deseaonalized       r1   
 Min.   :0     Min.   :0  
 1st Qu.:0     1st Qu.:0  
 Median :0     Median :0  
 Mean   :0     Mean   :0  
 3rd Qu.:0     3rd Qu.:0  
 Max.   :0     Max.   :0  
> 
> plot(out6a$simulation$r1, type='l')
> lines.default(out6b$simulation$r1, col=3)
> rug( which(out6b$simulation$r1 != out6a$simulation$r1))
> 
> days_diff <- matrix(out6a$simulation$r1 != out6b$simulation$r1,nrow=365)
> kap_par <- data.frame(out6a$pars)
> thresh <- kap_par$xi + kap_par$alfa*(1 - kap_par$h^(-kap_par$k))/kap_par$k
> image( cbind( is.na(thresh), days_diff))
> 
> 
> 
> proc.time()
   user  system elapsed 
  47.48    6.76   54.23