##### This file gives the R code for the main simulation study in paper "serra son" #####

#########################################################################################
### Helper packages and functions #######################################################
#########################################################################################

library(ismev)     ## ismev package, just in case
library(evd)       ## evd package, just in case

## Jonathan and Dave's function for quick (non-buggy) simulation of GPD series with logistic dependence 
#First, Gumbel margins::
gumbelsim<-function(n,alpha,seed){      #n=length of series, alpha=logistic dependence, seed=anything
        f<-function (z,x,v){
                ((1+exp(-z/alpha))**(alpha-1))*exp(-exp(-x)*((1+exp(-z/alpha))**alpha-1))-v
                }                       
        sim<-vector("numeric",n)
        x<-seed
        for (i in 1:n){
                        v<-runif(1)
                        y<-uniroot(f,c(-100,100),tol=0.0001,x,v)
                        sim[i]<-x+y$root
                        x<-sim[i]
                }
        return(sim)     
}
#Now transform to GPD margins::
GPDsim<-function(n,alpha,sigma,xi,seed){     
    test<-gumbelsim(n,alpha,seed)
    test2<-(sigma/xi)*((1-exp(-exp(-test)))**(-xi)-1)
    return(test2)
  }

#Return level function::
lee.q<-function(sigma,xi,theta,thresh,lambda,period,obs)
  {
    q=thresh+(sigma/xi)*(((1/lambda)*((1-(1-(1/(period*obs)))**(1/theta))))**(-xi)-1)
    return(q)
  }

##MCMC function for bivariate threshold excess model: 
#bayes6(n,dataset,U,sigmastart,xistart,alphastart,erreta,errxi,erralpha,sdeta,sdxi,lik.tol)
source("logistic_mcmc.R")

#########################################################################################
### Simulation study ####################################################################
#########################################################################################

#First, for the logistic model with dependence alpha - see "evd" manual for model codes 
MODEL = "log"

#######################################################################################
######## Main settings - sample sizes, no. of iterations and no. of replications::#####
#######################################################################################
N = 10^4                ### Length of reference series (years)  (Aim for 10^4)
Ny = 365.25*1           ### No. of obs. per year (hourly here) (Aim for 365.25*24)
n = 10                  ### Length of each simulated series (years)  (Aim for 10)
ny = Ny                 ### No. of obs. per year in each simulated series (hourly here)

L = 1000                ### No. of replications in the simulation study (Aim for 1000)

S = 5000 	        ### No. of MCMC iterations (Aim for 5000 - MUST BE > 1000 GIVEN THE BURN-IN DISCARD!)

#######################################################################################
######## Sub-settings - marginal and dependence parameters, thresholds::###############
#######################################################################################
scale = 1                        ### Hold marginal scale unit constant.  Why not?!
xi = c(-0.4, -0.1, 0, 0.1, 0.3)  ### Different marginal shapes to use

alpha = seq(0.9, 0.1, -0.1)                                    ### Logistic dependence parameters
theta = 0.013 - 0.092*alpha +1.833*(alpha^2) - 0.756*(alpha^3) ### 'True' values for the extremal index

p = c(0.9, 0.95, 0.99)           ### 1 - (Threshold exceedance probability)

r = c(10, 50, 200, 1000, 10000)  ### Specific return periods of interest

########################################################################################
######## Specific settings, choosing elements from the vectors above::##################
########################################################################################
dep = alpha[1]
shape = xi[1]
exindex = theta[1]
excess_prob = p[2]


##################################
####### Reference dataset::#######
##################################
#y_Ref = GPDsim(n = (1000+N*Ny), alpha = dep, sigma = scale, xi = shape, seed = 2016)
#y_Ref = y_Ref[1001: length(y_Ref)]

##################################
### Results vectors::#############
##################################
scale_star.post.mean = vector("numeric", length = L)    ### Posterior means of sigma_star
shape_star.post.mean = vector("numeric", length = L)       ### Posterior means of xi_star
dep.post.mean = vector("numeric", length = L)           ### Posterior means of dependence parameter
extremal.index.post.mean = vector("numeric",length = L)      ### Posterior means of extremal index (via polynomials in case of asympt. dep.)
# Return levels::
cn = c(r[1], r[2], r[3], r[4], r[5])
z_post.mean = matrix(ncol = length(r), nrow = L)        ### Posterior means of return levels of interest
z_post.mode = matrix(ncol = length(r), nrow = L)        ### Posterior modes of return levels of interest
z_post.upper = matrix(ncol = length(r), nrow = L)       ### Posterior 95% upper bounds of return levels of interest
z_pred = matrix(ncol = length(r), nrow = L)             ### Posterior predictive return levels of interest
colnames(z_post.mean) = cn
colnames(z_post.mode) = cn
colnames(z_post.upper) = cn
colnames(z_pred) = cn

##################################
####### Now the replications::####
##################################
for(l in 1:L){
	print(l)
	#Simulated series, typical in length of "real" datasets::
	y_l = GPDsim(1000+(n*ny), alpha = dep, sigma = scale, xi = shape, seed = 2016)
	y_l = y_l[1001: length(y_l)]	

	#For temporal bivariate modelling, split into two columns of neighbouring values::
	y_l.m = matrix(ncol = 2, nrow = (n*ny-1))
	y_l.m[,1] = y_l[1: (length(y_l)-1)];   y_l.m[,2] = y_l[2: length(y_l)]

	#Threshold - needs to be of length 2 for bivariate modelling::
	U = rep(as.numeric(quantile(y_l, excess_prob)), 2)
	
	#ML estimation, just to check (and to obtain decent statting values in the MCMC)::
	fit.ML = fbvpot(y_l.m, U, model = MODEL, std.err = FALSE)

	################################
	###Bayesian fitting via MCMC::##
	################################
	#Starting values::
	sigma.start = c(fit.ML$estimate[1], fit.ML$estimate[3]); xi.start = c(fit.ML$estimate[2], fit.ML$estimate[4]); alpha.start = fit.ML$estimate[5]
	# Random walk innovation variances for tuning::
	err.eta = c(0.15, 0.15)
	err.xi = c(0.15, 0.15)
	err.alpha = 0.15
	# St. dev. in the Normal priors for log(sigma) and xi::
	sd.eta = c(1000, 1000)
	sd.xi = c(100, 100)
	# Rounding tolerance in likelihood function (probably don't touch)::
	T = 10
	# Now the MCMC::	
	fit.bayes = bayes6(S, y_l.m, U, sigmastart=sigma.start, xistart=xi.start, alphastart=alpha.start, erreta=err.eta, errxi=err.xi, erralpha=err.alpha, sdeta=sd.eta, sdxi=sd.xi, lik.tol = T)

	#################################################################
	###Extract the posterior means for each parameter of interest::##
	#################################################################
	# Allow for some burn-in; should be minimal/non-existent, given MLE starting values::	
	
	# Marginal parameters::
	mcmc.sigma = exp(fit.bayes[1001:S, 1])
	mcmc.xi = fit.bayes[1001:S, 3]
	mcmc.alpha = fit.bayes[1001:S, 5]
	mcmc.theta = 0.013 - 0.092*mcmc.alpha + 1.833*(mcmc.alpha^2) - 0.756*(mcmc.alpha^3)

	scale_star.post.mean[l] = mean(mcmc.sigma)     #Doesn't matter, really, if we use column 1 or 2         
	shape_star.post.mean[l] = mean(mcmc.xi)        #Doesn't matter, really, if we use column 3 or 4
	dep.post.mean[l] = mean(mcmc.alpha)
	extremal.index.post.mean[l] = mean(mcmc.theta)
	
	# Return levels::
	LAMBDA = length(y_l[y_l > U[1]])/length(y_l)
	for(i in 1:ncol(z_pred)){

		z = lee.q(mcmc.sigma, mcmc.xi, mcmc.theta, U[1], LAMBDA, r[i], obs = ny)

		z_post.mean[l, i] =  mean(z)
		z_post.upper[l, i] = quantile(z, p=0.975)
 
		A = density(z)
		z_post.mode[l, i] = A$x[which(A$y == max(A$y))]

	#Now for the predictive return level:
	f = function(zp){
		inner = vector("numeric", length(mcmc.sigma))
		core = vector("numeric", length(inner))
		for(j in 1:length(inner)){
			inner[j] = max(0, (1 + (mcmc.xi[j]/mcmc.sigma[j])*(zp - U[1])))
			core[j] = (1-LAMBDA*((inner[j])**(-1/mcmc.xi[j])))**(mcmc.theta[j])}
		answer = mean(core) - (1-1/(ny*r[i]))
		return(answer)}
	B = uniroot(f, lower = 2, upper = 6)
	
		z_pred[l,i] = B$root }

	if(l==L){print("See 'scale_star.post.mean','shape_star.post.mean', 'dep.post.mean', 'extremal.index.post.mean', 'z_post.mean', 'z_post.mode', 'z_post.upper'")}

write.table(z_post.mean,"z-mean.txt",quote=FALSE, col.names=TRUE, row.names=FALSE)
write.table(z_post.mode,"z-mode.txt",quote=FALSE, col.names=TRUE, row.names=FALSE)
write.table(z_post.upper,"z-upper.txt",quote=FALSE, col.names=TRUE, row.names=FALSE)
write.table(z_pred,"z-pred.txt",quote=FALSE, col.names=TRUE, row.names=FALSE)
write.table(scale_star.post.mean,"scale-mean.txt",quote=FALSE, col.names=FALSE, row.names=FALSE)
write.table(shape_star.post.mean,"shape-mean.txt",quote=FALSE, col.names=FALSE, row.names=FALSE)
write.table(dep.post.mean,"alpha-mean.txt",quote=FALSE, col.names=FALSE, row.names=FALSE)
write.table(extremal.index.post.mean,"theta-mean.txt",quote=FALSE, col.names=FALSE, row.names=FALSE)
	}