#####################################################################
#####################################################################
##                                                                 ##
## R CODE FOR GLUEVAR RISK MEASURES APPLICATIONS                   ##
## JAUME BELLES SAMPERA - MONTSERRAT GUILLÉN - MIGUEL SANTOLINO    ##
## RISKCENTER - UNIVERSITAT DE BARCELONA                           ##
## http://www.ub.edu/riskcenter                                    ##
## DECEMBER 2014                                                   ##
## Version 1.1                                                     ## 
##                                                                 ##
#####################################################################
#####################################################################


################## LOADING PACKAGES

library('Matrix')


################## EMPIRICAL PROBABILITY FUNCTIONS

################## Empirical Cumulative Distribution Function
CDFemp<- function(data)
{
	df<-dfemp(data)
	cdf<-matrix(0,ncol=ncol(df), nrow=nrow(df))
	cdf[1,1]<-df[1,1]
	cdf[1,2]<-df[1,2]
	for(i in 2:nrow(df))
		{
		cdf[i,1]<-df[i,1]
		cdf[i,2]<-df[i,2]+cdf[i-1,2]
		}
	return(cdf)
}

################## Empirical Cumulative Survival Function 
CSFemp<- function(data)
{		
	cdf<-CDFemp(data)
	csf<-matrix(0,ncol=ncol(cdf), nrow=nrow(cdf))
	
	for(i in 1:nrow(cdf))
		{
		csf[i,1]<-cdf[i,1]
		csf[i,2]<-1-cdf[i,2]
		}
	return(csf)

}


################## Empirical density function
dfemp<-function(data)
{
	###### Controlling N/A
	data_dep<-0
	for(i in 1:length(data))
	{
          if(is.na(data[i]) != TRUE) data_dep<-c(data_dep, data[i])
	}
	data_dep<-data_dep[-1]
	######

	n<-length(data_dep)
	inivec<-cbind(sort(data_dep),1/n)

	counter<-1
	aux<-inivec[1]
	for(i in 1:n)
	{
	 if(inivec[i]!=aux)
		{
		counter<-counter+1
		aux<-inivec[i]
		}
	}	

	dens<-matrix(0,ncol=2,nrow=n)
	for(i in 1:n)
	{
	  auxdens<-inivec[i]
	  dens[i,1]<-auxdens
	  dens[i,2]<-sum((inivec[,2]) [inivec[,1] == auxdens])
	}

	densres<-matrix(0,ncol=2,nrow=counter)
	aux<-inivec[1]
	auxcont<-1
	for(i in 1:n)
	{
	  if(inivec[i]!=aux)
		{
		densres[auxcont,1]<-aux
		densres[auxcont,2]<-dens[i-1,2]
		auxcont<-auxcont+1
		aux<-inivec[i]
		}
	}
	densres[counter,1]<-dens[n,1]
	densres[counter,2]<-dens[n,2]
	
	return(densres)
}


################# GLUEVAR PARAMETER MATRICES

gluevar_H<-function(alpha, beta)
{	
	H<-matrix(c(1,1,(1-beta)/(1-alpha),1), ncol=2, nrow=2)
	return(H)
}

gluevar_Hinv<-function(alpha, beta)
{	
	Hinv<-matrix(c((1-alpha)/(1-beta), (alpha-1)/(beta-alpha), (beta-1)/(beta-alpha), (1-alpha)/(beta-alpha)), ncol=2, nrow=2)
	return(Hinv)
}

gluevar_T<-function(alpha, beta)
{	
	T<-matrix(c(1-beta,0,beta-1,1-alpha), ncol=2, nrow=2)
	return(T)
}

gluevar_Tinv<-function(alpha, beta)
{
	Tinv<-matrix(c(1/(1-beta), 0, 1/(1-alpha), 1/(1-alpha)), ncol=2, nrow=2)
	return(Tinv)
}

gluevar_M<-function(alpha, beta)
{
	M<-matrix(c(1-beta,1-beta,0,beta-alpha), ncol=2, nrow=2)
	return(M)
}

gluevar_Minv<-function(alpha, beta)
{
	Minv<-matrix(c(1/(1-beta),-1/(beta-alpha),0,1/(beta-alpha)), ncol=2, nrow=2)
	return(Minv)
}

gluevarmatrix<-function(alpha, beta)
{
	H<-gluevar_H(alpha,beta)
	Hinv<-gluevar_Hinv(alpha,beta)

	T<-gluevar_T(alpha,beta)
	Tinv<-gluevar_Tinv(alpha,beta)

	M<-gluevar_M(alpha,beta)
	Minv<-gluevar_Minv(alpha,beta)

	print("-----------------------------------", quote=FALSE)
	print("-------- GLUEVAR PARAMETERS -------", quote=FALSE)
	print("-----------------------------------", quote=FALSE)
	print("", quote=FALSE)
	print("Weights to heights. Matrix H:", quote=FALSE)
	print(H, quote=FALSE)
	print("", quote=FALSE)
	print("Heights to weights. Matrix H inverse:", quote=FALSE)
	print(Hinv, quote=FALSE)
	print("", quote=FALSE)
	print("Slopes to weights. Matrix T:", quote=FALSE)
	print(T, quote=FALSE)
	print("", quote=FALSE)
	print("Weights to slopes. Matrix T inverse:", quote=FALSE)
	print(Tinv, quote=FALSE)
	print("", quote=FALSE)
	print("Slopes to heights. Matrix M:", quote=FALSE)
	print(M, quote=FALSE)
	print("", quote=FALSE)
	print("Heights to slopes. Matrix M inverse:", quote=FALSE)
	print(Minv, quote=FALSE)
	print("", quote=FALSE)
	print("-----------------------------------", quote=FALSE)
	print("-----------------------------------", quote=FALSE)

}


################# EMPIRICAL RISK CALCULATIONS

gluevar_emp_risk<-function(filename, header, obs, alpha, beta, w1, w2, w3)
{
	datafile<-read.table(filename, sep=";", header=TRUE)

	risks<-length(header)	
	data<-matrix(0, ncol=risks, nrow=obs)
	for(j in 1:risks)
	{
	 for(i in 1:obs)
		 data[i,j]<-datafile[i,j]
	}
			
	print("-----------------------------------", quote=FALSE)
	print("--- EMPIRICAL RISK CALCULATIONS ---", quote=FALSE)
	print("-----------------------------------", quote=FALSE)

	print("", quote=FALSE)
	print(paste("Alpha confidence level: ", alpha), quote=FALSE)
	print(paste("Beta confidence level: ", beta), quote=FALSE)
	print(paste("VaR alpha weight: ", round(w3,4)), quote=FALSE)
	print(paste("TVaR alpha weight: ", round(w2,4)), quote=FALSE)
	print(paste("TVaR beta weight: ", round(w1,4)), quote=FALSE)

	rnames<-c("VaR_alpha","TVaR_alpha","TVaR_beta","GlueVaR")
	
	taula<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnames,header))
	
	for(j in 1:risks)
	{
		distrib<-CDFemp(data[,j])
		dens<-dfemp(data[,j])
		nobs<-length(distrib[,1])
		
		rank.alpha<-0
		rank.beta<-0
		ES.alpha<-0
		ES.beta<-0
		VaR.alpha<-0
		VaR.beta<-0
		TVaR.alpha<-0
		TVaR.beta<-0
	
		#VaR alpha
		rank.alpha<-min(which(distrib[,2]>=alpha))
		VaR.alpha<-distrib[rank.alpha,1]
		taula[1,j]<-round(VaR.alpha,4)
		if(rank.alpha<nobs)
		{
		  for(i in (rank.alpha+1):nobs) 
			{
			ES.alpha<-ES.alpha + (distrib[i,1] - VaR.alpha)*dens[i,2]
			}
		}

		#TVaR alpha
		TVaR.alpha<-VaR.alpha+(1/(1-alpha))*ES.alpha
		taula[2,j]<-round(TVaR.alpha,4)
		
		#TVaR beta
		rank.beta<-min(which(distrib[,2]>=beta))
		VaR.beta<-distrib[rank.beta,1]
		if(rank.beta<nobs)
		{
		  for(i in (rank.beta+1):nobs) 
			{
			ES.beta<-ES.beta + (distrib[i,1] - VaR.beta)*dens[i,2]
			}
		}
		TVaR.beta<-VaR.beta+(1/(1-beta))*ES.beta
		taula[3,j]<-round(TVaR.beta,4)

		#GlueVaR(beta,alpha,w1,w2)
		taula[4,j]<-round(w1*TVaR.beta+w2*TVaR.alpha+w3*VaR.alpha ,4)

	}

	
	print("", quote=FALSE)
	print(taula, quote=FALSE)
	
	return(taula)
}


################# COMMON STATISTICAL DISTRIBUTIONS CALCULATIONS

gluevar_distributions<-function(filename, header, obs, alpha, beta, w1, w2, w3)
{
	datafile<-read.table(filename, sep=";", header=TRUE)

	risks<-length(header)	
	data<-matrix(0, ncol=risks, nrow=obs)
	for(j in 1:risks)
	{
	 for(i in 1:obs)
		 data[i,j]<-datafile[i,j]
	}
			
	print("-----------------------------------", quote=FALSE)
	print("---- STATISTICAL DISTRIBUTIONS ----", quote=FALSE)
	print("-----------------------------------", quote=FALSE)

	print("", quote=FALSE)
	print(paste("Alpha confidence level: ", alpha), quote=FALSE)
	print(paste("Beta confidence level: ", beta), quote=FALSE)
	print(paste("VaR alpha weight: ", round(w3,4)), quote=FALSE)
	print(paste("TVaR alpha weight: ", round(w2,4)), quote=FALSE)
	print(paste("TVaR beta weight: ", round(w1,4)), quote=FALSE)

	rnamesT2<-c("VaR_alpha","TVaR_alpha","TVaR_beta","GlueVaR")
	T2<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnamesT2,header))
	
	##### NORMAL RISK ESTIMATION
	print("--------------------------------------------", quote=FALSE)
	print("---- RISK ASSUMING NORMAL DISTRIBUTIONS ----", quote=FALSE)
	print("--------------------------------------------", quote=FALSE)

	rnamesT1<-c("mu.est","(sigma^2).est")
	T1<-matrix(0, ncol=risks, nrow=2, dimnames=list(rnamesT1,header))	

	for(i in 1:risks)
	{
		###### Controlling N/A
		data_dep<-0
		for(j in 1:length(data[,i]))		
		{
          	  if(is.na(data[j,i]) != TRUE) data_dep<-c(data_dep, data[j,i])		
		}
		data_dep<-data_dep[-1]
		######		

		mu.est<-mean(data_dep)
		sigma2.est<-var(data_dep)
			
		VaR.alpha<-mu.est+sqrt(sigma2.est)*qnorm(alpha,0,1)
		TVaR.alpha<-mu.est+sqrt(sigma2.est)*dnorm(qnorm(alpha,0,1),0,1)/(1-alpha)
		TVaR.beta<-mu.est+sqrt(sigma2.est)*dnorm(qnorm(beta,0,1),0,1)/(1-beta)
		
		T1[1,i]<-round(mu.est,4)
		T1[2,i]<-round(sigma2.est,4)
		
		T2[1,i]<-round(VaR.alpha,4)
		T2[2,i]<-round(TVaR.alpha,4)
		T2[3,i]<-round(TVaR.beta,4)
		T2[4,i]<-round(w1*TVaR.beta+w2*TVaR.alpha+w3*VaR.alpha,4)
	}

	print("", quote=FALSE)

	print("---- Parameter estimates - Normal ----", quote=FALSE)
	print("", quote=FALSE)
	print(T1, quote=FALSE)
	print("", quote=FALSE)

	print("---- Risk estimation - Normal ----", quote=FALSE)
	print("", quote=FALSE)
	print(T2, quote=FALSE)
	print("", quote=FALSE)
		

	##### LOG-NORMAL RISK ESTIMATION
	print("------------------------------------------------", quote=FALSE)
	print("---- RISK ASSUMING LOG-NORMAL DISTRIBUTIONS ----", quote=FALSE)
	print("------------------------------------------------", quote=FALSE)
	
	print("", quote=FALSE)
	print("Note: if data take negative values, results are not reliable.", quote=FALSE)

	rnamesT1<-c("mu.est","(sigma^2).est")
	T1<-matrix(0, ncol=risks, nrow=2, dimnames=list(rnamesT1,header))	

	for(i in 1:risks)
	{
		###### Controlling N/A
		data_dep<-0
		for(j in 1:length(data[,i]))		
		{
          	  if(is.na(data[j,i]) != TRUE) data_dep<-c(data_dep, data[j,i])		
		}
		data_dep<-data_dep[-1]
		######		
		
		est.mean<-mean(data_dep)
		est.var<-var(data_dep)
		mu.est<-log(est.mean)-0.5*log(est.var/(est.mean^2)+1)
		sigma2.est<-log(est.var/(est.mean^2)+1)


		VaR.alpha<-exp(mu.est+sqrt(sigma2.est)*qnorm(alpha,0,1))
		TVaR.alpha<-exp(mu.est+sigma2.est/2)*pnorm(sqrt(sigma2.est)-qnorm(alpha,0,1),0,1)/(1-alpha)
		TVaR.beta<-exp(mu.est+sigma2.est/2)*pnorm(sqrt(sigma2.est)-qnorm(beta,0,1),0,1)/(1-beta)
		
		T1[1,i]<-round(mu.est,4)
		T1[2,i]<-round(sigma2.est,4)

		T2[1,i]<-round(VaR.alpha,4)
		T2[2,i]<-round(TVaR.alpha,4)
		T2[3,i]<-round(TVaR.beta,4)
		T2[4,i]<-round(w1*TVaR.beta+w2*TVaR.alpha+w3*VaR.alpha,4)
	}


	print("", quote=FALSE)

	print("---- Parameter estimates - Log-normal ----", quote=FALSE)
	print("", quote=FALSE)
	print(T1, quote=FALSE)
	print("", quote=FALSE)

	print("---- Risk estimation - Log-normal ----", quote=FALSE)
	print("", quote=FALSE)
	print(T2, quote=FALSE)
	print("", quote=FALSE)	


	##### STUDENT'S-T RISK ESTIMATION
	print("-------------------------------------------------", quote=FALSE)
	print("---- RISK ASSUMING STUDENT'S-T DISTRIBUTIONS ----", quote=FALSE)
	print("-------------------------------------------------", quote=FALSE)
	
	rnamesT1<-c("mu.est","(sigma^2).est")
	T1<-matrix(0, ncol=risks, nrow=2, dimnames=list(rnamesT1,header))	

	###### Look for different degrees of freedom: 3,4,7,10
	T2<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnamesT2,header))
	T3<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnamesT2,header))
	T4<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnamesT2,header))
	T5<-matrix(0, ncol=risks, nrow=4, dimnames=list(rnamesT2,header))

	for(i in 1:risks)
	{
		###### Controlling N/A
		data_dep<-0
		for(j in 1:length(data[,i]))		
		{
          	  if(is.na(data[j,i]) != TRUE) data_dep<-c(data_dep, data[j,i])		
		}
		data_dep<-data_dep[-1]
		######		

		mu.est<-mean(data_dep)
		sigma2.est<-var(data_dep)

		###### 3 df
		VaR.alpha.3df<-mu.est+sqrt(sigma2.est)*qt(alpha,df=3)
		TVaR.alpha.3df<-mu.est+sqrt(sigma2.est)*dt(qt(alpha,df=3),df=3)*(3+qt(alpha,df=3)^2)/((1-alpha)*2)
		TVaR.beta.3df<-mu.est+sqrt(sigma2.est)*dt(qt(beta,df=3),df=3)*(3+qt(beta,df=3)^2)/((1-beta)*2)

		###### 4 df
		VaR.alpha.4df<-mu.est+sqrt(sigma2.est)*qt(alpha,df=4)
		TVaR.alpha.4df<-mu.est+sqrt(sigma2.est)*dt(qt(alpha,df=4),df=4)*(4+qt(alpha,df=4)^2)/((1-alpha)*3)
		TVaR.beta.4df<-mu.est+sqrt(sigma2.est)*dt(qt(beta,df=4),df=4)*(4+qt(beta,df=4)^2)/((1-beta)*3)
		
		###### 7 df
		VaR.alpha.7df<-mu.est+sqrt(sigma2.est)*qt(alpha,df=7)
		TVaR.alpha.7df<-mu.est+sqrt(sigma2.est)*dt(qt(alpha,df=7),df=7)*(7+qt(alpha,df=7)^2)/((1-alpha)*6)
		TVaR.beta.7df<-mu.est+sqrt(sigma2.est)*dt(qt(beta,df=7),df=7)*(7+qt(beta,df=7)^2)/((1-beta)*6)
		
		###### 10 df
		VaR.alpha.10df<-mu.est+sqrt(sigma2.est)*qt(alpha,df=10)
		TVaR.alpha.10df<-mu.est+sqrt(sigma2.est)*dt(qt(alpha,df=10),df=10)*(10+qt(alpha,df=10)^2)/((1-alpha)*9)
		TVaR.beta.10df<-mu.est+sqrt(sigma2.est)*dt(qt(beta,df=10),df=10)*(10+qt(beta,df=10)^2)/((1-beta)*9)
		
		
		T1[1,i]<-round(mu.est,4)
		T1[2,i]<-round(sigma2.est,4)

		T2[1,i]<-round(VaR.alpha.3df,4)
		T2[2,i]<-round(TVaR.alpha.3df,4)
		T2[3,i]<-round(TVaR.beta.3df,4)
		T2[4,i]<-round(w1*TVaR.beta.3df+w2*TVaR.alpha.3df+w3*VaR.alpha.3df,4)

		T3[1,i]<-round(VaR.alpha.4df,4)
		T3[2,i]<-round(TVaR.alpha.4df,4)
		T3[3,i]<-round(TVaR.beta.4df,4)
		T3[4,i]<-round(w1*TVaR.beta.4df+w2*TVaR.alpha.4df+w3*VaR.alpha.4df,4)

		T4[1,i]<-round(VaR.alpha.7df,4)
		T4[2,i]<-round(TVaR.alpha.7df,4)
		T4[3,i]<-round(TVaR.beta.7df,4)
		T4[4,i]<-round(w1*TVaR.beta.7df+w2*TVaR.alpha.7df+w3*VaR.alpha.7df,4)

		T5[1,i]<-round(VaR.alpha.10df,4)
		T5[2,i]<-round(TVaR.alpha.10df,4)
		T5[3,i]<-round(TVaR.beta.10df,4)
		T5[4,i]<-round(w1*TVaR.beta.10df+w2*TVaR.alpha.10df+w3*VaR.alpha.10df,4)


	}


	print("", quote=FALSE)

	print("---- Parameter estimates - Student's-t ----", quote=FALSE)
	print("", quote=FALSE)
	print(T1, quote=FALSE)
	print("", quote=FALSE)

	print("---- Risk estimation - Student's-t: 3df ----", quote=FALSE)
	print("", quote=FALSE)
	print(T2, quote=FALSE)
	print("", quote=FALSE)	

	print("---- Risk estimation - Student's-t: 4df ----", quote=FALSE)
	print("", quote=FALSE)
	print(T3, quote=FALSE)
	print("", quote=FALSE)	

	print("---- Risk estimation - Student's-t: 7df ----", quote=FALSE)
	print("", quote=FALSE)
	print(T4, quote=FALSE)
	print("", quote=FALSE)	

	print("---- Risk estimation - Student's-t: 10df ----", quote=FALSE)
	print("", quote=FALSE)
	print(T5, quote=FALSE)
	print("", quote=FALSE)	


	##### GENERALIZED PARETO RISK ESTIMATION
	print("--------------------------------------------------------", quote=FALSE)
	print("---- RISK ASSUMING GENERALIZED PARETO DISTRIBUTIONS ----", quote=FALSE)
	print("--------------------------------------------------------", quote=FALSE)
	
	print("", quote=FALSE)
	print("Note: if data take negative values, results are not reliable.", quote=FALSE)

	rnamesT1<-c("k.est","sigma.est")
	T1<-matrix(0, ncol=risks, nrow=2, dimnames=list(rnamesT1,header))	

	for(i in 1:risks)
	{
		###### Controlling N/A
		data_dep<-0
		for(j in 1:length(data[,i]))		
		{
          	  if(is.na(data[j,i]) != TRUE) data_dep<-c(data_dep, data[j,i])		
		}
		data_dep<-data_dep[-1]
		######		
		

		est.mean<-mean(data_dep)
		est.var<-var(data_dep)
		
		k.est<-0.5*((est.mean^2)/est.var-1)
		sigma.est<-0.5*est.mean*((est.mean^2)/est.var+1)

		if(round(k.est,4) != 0) #### If Exponential distribution is discarded
		{
		  if(round(k.est,4)>0) #### Type II Pareto distribution
		  {
			VaR.alpha<- (sigma.est/k.est)*(1-(1-alpha)^k.est)
		  	TVaR.alpha<- (sigma.est/k.est)*(1-(1-alpha)^k.est) + (sigma.est/k.est)*(k.est*((1-alpha)^k.est)/(k.est+1))
		  	TVaR.beta<- (sigma.est/k.est)*(1-(1-beta)^k.est) + (sigma.est/k.est)*(k.est*((1-beta)^k.est)/(k.est+1))
		   }
		  else
		    {
			if(round(k.est,4)>-1) #### Pareto distribution with finite TVaR
			{
				VaR.alpha<- (sigma.est/k.est)*(1-(1-alpha)^k.est)
			  	TVaR.alpha<- (sigma.est/k.est)*(1-(1-alpha)^k.est) + (sigma.est/k.est)*(k.est*((1-alpha)^k.est)/(k.est+1))
			  	TVaR.beta<- (sigma.est/k.est)*(1-(1-beta)^k.est)+ (sigma.est/k.est)*(k.est*((1-beta)^k.est)/(k.est+1))
			}

			else #### Pareto distribution with infinite TVaR
			{
				VaR.alpha<- (sigma.est/k.est)*(1-(1-alpha)^k.est)
			  	TVaR.alpha<- NA
			  	TVaR.beta<- NA				
			}
		    }
		}
		else #### Exponential distribution
		{
		
		  VaR.alpha<- -sigma.est*log(1-alpha)
		  TVaR.alpha<- sigma.est*(1-log(1-alpha))
		  TVaR.beta<- sigma.est*(1-log(1-beta))

		}
		
		T1[1,i]<-round(k.est,4)
		T1[2,i]<-round(sigma.est,4)

		T2[1,i]<-round(VaR.alpha,4)
		T2[2,i]<-round(TVaR.alpha,4)
		T2[3,i]<-round(TVaR.beta,4)
		T2[4,i]<-round(w1*TVaR.beta+w2*TVaR.alpha+w3*VaR.alpha,4)
	}


	print("", quote=FALSE)

	print("---- Parameter estimates - GPD ----", quote=FALSE)
	print("", quote=FALSE)
	print(T1, quote=FALSE)
	print("", quote=FALSE)

	print("---- Risk estimation - GDP ----", quote=FALSE)
	print("", quote=FALSE)
	print(T2, quote=FALSE)
	print("", quote=FALSE)	



}


#### Empirical VaR as Choquet Integral 
emp_VaR_asCI<-function(filename,obs,cols1,alpha)
{
	datafile<-read.table(filename, sep=";", header=TRUE)

	data<-matrix(0, ncol=1, nrow=obs)
	
	data[,1]<-sort(datafile[,cols1])
	dens<-dfemp(data[,1])
	n<-length(dens[,1])
	rm<-0

	## We first calculate g() for VaR_alpha
	gPA<-matrix(0,nrow=n+1,ncol=1)

	for(i in 1:n)
	{ 
		PA<-0
		for(k in i:n) PA<-PA+dens[k,2]

		if(PA>=1-alpha){
		 	gPA[i,1]<-1
		 	}
	}	
	#### Vector of weights to calculate DRM
	for(i in 1:n) rm<-rm+dens[i,1]*(gPA[i,]-gPA[i+1,])
	return(rm)
}

#### Empirical TVaR as Choquet Integral 
emp_TVaR_asCI<-function(filename,obs,cols1,alpha)
{
	datafile<-read.table(filename, sep=";", header=TRUE)

	data<-matrix(0, ncol=1, nrow=obs)
	
	data[,1]<-sort(datafile[,cols1])
	dens<-dfemp(data[,1])
	n<-length(dens[,1])
	rm<-0

	## We first calculate g() for TVaR_alpha
	gPA<-matrix(0,nrow=n+1,ncol=1)

	for(i in 1:n)
	{ 
		PA<-0
		for(k in i:n) PA<-PA+dens[k,2]

		if(PA>=1-alpha){
		 	gPA[i,1]<-1
		 	}
		else gPA[i,1]<-PA/(1-alpha)
	}	
	#### Vector of weights to calculate DRM
	for(i in 1:n) rm<-rm+dens[i,1]*(gPA[i,]-gPA[i+1,])
	return(rm)
}


################################################################################
################################################################################

############ EXAMPLE: HOW TO USE THESE FUNCTIONS ###############################

# Data management

# Load selected data
data <- read.csv2("IBEX.csv")
data <- read.csv2("CAC.csv") 
data <- read.csv2("DAX.csv") 

val <- as.matrix(data[,2]) 
n <- nrow(val)

val.ln <- -diff(log(val)) # Factor risk changes

write.table(val.ln, file="data.csv", dec=".", sep=";")

filenamed<-"data.csv"
maxnumdata<-length(data[,1])
headtitle<-c("Index") 

alpha<-0.95
beta<-0.995

w1<-0.33333
w2<-0.33333
w3<-1-w1-w2

gluevarmatrix(alpha, beta)

# Empirical GlueVaR calculation
emprisk<-gluevar_emp_risk(filenamed, headtitle, maxnumdata, alpha, beta, w1, w2, w3)