rm(list=ls(all=TRUE))
setwd("H:/CoVaR05")
# setwd("E:/Researches/20130212 CCP/Code/CoVaR05")
#install.packages("quantreg")
#install.packages("KernSmooth")
library(quantreg)
library(KernSmooth)
library(foreign)
################ Source code ###############

# source("MD_tests.r")
source("spacin_lprq.r")

################ Data #####################

RET<-read.csv("CDS_G14.csv") # CDS returns
 MV<-read.table("Market_Data_Total_WithoutCDS.txt") # Market variables
var<-read.csv("VaR05_Gesamt_Lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")
  BE<-array(list(0),c(13,1))  # Lasso betas
 BE[[1,1]]<-read.csv("02_Citi_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")     # Citi      # ,sep=";",dec=","
 BE[[2,1]]<-read.csv("03_BoA_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[3,1]]<-read.csv("04_Barclays_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[4,1]]<-read.csv("05_BNP_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[5,1]]<-read.csv("06_CreditSuisse_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[6,1]]<-read.csv("07_DB_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[7,1]]<-read.csv("08_GS_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[8,1]]<-read.csv("09_HSBC_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[9,1]]<-read.csv("10_JPM_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[10,1]]<-read.csv("11_MS_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[11,1]]<-read.csv("12_RBS_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[12,1]]<-read.csv("13_SocGen_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")  
 BE[[13,1]]<-read.csv("14_UBS_lasso-qr_para33_ohneCDS.csv",sep=";",dec=",")   

#### Remove the first "date" column
 RET<-RET[-2:-1,-1]
 MV <-as.matrix(MV[,-1]) 
 var<-as.matrix(var[,-1])
 for(l in 1:13){
 BE[[l,1]]<-as.matrix(BE[[l,1]][,-1]) 
 } 
### # Transfer from bp to decimals
# RET<-RET*0.01
# MV<-MV*0.01
############### Product ########################

BETA<-array(list(0),c(13,13))
COVAR<-array(list(0),c(13,13))

############## Parameters #####################

   q  <- 0.01
n_all <- length(RET[,1])        # Length of the sample
winsize <- 250                     # no. of samples in estimation
winlimit<- n_all-winsize
n_var <- 33                      # no. of regressors xcl. other asset returns (incl. intercept, market or firm specific variables)

########### Main code #########################

for(i in 1:13){
for(j in 1:13){
if(i != j){

# betatemp <- matrix(0,nrow=n_test,ncol=n_var) # temp storage for coefficients
     covar <- numeric(0)                       # temp storage for estimates
         T <- winsize
#      vio <- numeric(0)                       # temp storage for violations


for(seuil in 1:winlimit){

    VAR2 <- RET[(seuil+1):T,j]                   # Regressand
    VAR1 <- RET[(seuil+1):T,i]                   # Regrdssor
   selec <- which(BE[[j,1]][seuil,]!=0) 
  MV_est <- as.matrix(MV[seuil:(T-1),selec])          # Transfer from data to matrix  
 n_selec <- length(selec) 
   #### plm filtering procedure ########################
    orderVAR2<-VAR2[order(VAR1)]
    orderMV_est<-as.matrix(MV_est[order(VAR1),])
    leng <- 1/round((winsize-1)^(1/3)/2)                                   # length of each partition
    bbeta<- matrix(0,(1/leng),n_selec+1) 					   # +1 for interception term
    for(k in 1: (1/leng)){                                           # go through each parition for the parametric quantile estimate
      	
		uu   <- orderMV_est[(round((k-1)*leng*(winsize-1))+1): (round(k*leng*(winsize-1))),]
      	yy   <- orderVAR2[(round((k-1)*leng*(winsize-1))+1): (round((k)*leng*(winsize-1)))]
      	fit1 <- rq(yy~uu,tau=q)
      	bbeta[k,] <- coef(fit1)[1:(n_selec+1)]                             # beta estimate in each small interval is just the coefficient
    }


      beta <- as.matrix(colMeans(as.matrix(bbeta[,2:(n_selec+1)])))# c(mean(bbeta[,2]),mean(bbeta[,3]),mean(bbeta[,4]),mean(bbeta[,5]),mean(bbeta[,6]),mean(bbeta[,7]),mean(bbeta[,8])) 
   filVAR2 <- VAR2 - mean(bbeta[,1]) - MV_est%*%beta              
   ######### end of plm filtering procedure ################

#   betatemp[seuil,] <- beta
   
   fit_plm <- spacin_lprq(VAR1, VAR2, q, gridn = 60)
     l_tau <- approxfun(fit_plm$x0, fit_plm$est, method="linear")

   Mt <- t(as.matrix(MV[T,selec]))
    
         # (turn into "vector" type for inner product) Get current MV for prediction

  ### if the present asset return i is below the support of VAR1, use linear approximation #### 
     if(is.na(l_tau(var[seuil,i]))==FALSE){
		covar[seuil]<-l_tau(var[seuil,i])+Mt%*%beta 
	}
else{
		fit_out   <-lm(fit_plm$est[1:10]~fit_plm$x0[1:10])
		covar[seuil]<-coef(fit_out)[1]+coef(fit_out)[2]*var[seuil,i]+Mt%*%beta
	}
  #############################################################################

T<-T+1
}
BETA[[i,j]] <- beta
COVAR[[i,j]]<-covar 
}}}


COVAR1<-array(list(0),c(13,1))
for(h in 1:13){
COVAR1[[h,1]]<-cbind(COVAR[[1,h]],COVAR[[2,h]],COVAR[[3,h]],COVAR[[4,h]],COVAR[[5,h]],COVAR[[6,h]],COVAR[[7,h]],COVAR[[8,h]],COVAR[[9,h]],COVAR[[10,h]],COVAR[[11,h]],COVAR[[12,h]],COVAR[[13,h]])
}

for(n_y in 2:14){
file_path <- c("H:/")
if (n_y == 2){
  file_name <-  paste (file_path, "02_Citi_", sep="")
} else if (n_y == 3){
  file_name <-  paste (file_path, "03_BoA_", sep="")
} else if (n_y == 4){
  file_name <-  paste (file_path, "04_Barclays_", sep="")
} else if (n_y == 5){
  file_name <-  paste (file_path, "05_BNP_", sep="")
} else if (n_y == 6){
  file_name <-  paste (file_path, "06_CreditSuisse_", sep="")
} else if (n_y == 7){
  file_name <-  paste (file_path, "07_DB_", sep="")
} else if (n_y == 8){
  file_name <-  paste (file_path, "08_GS_", sep="")
} else if (n_y == 9){
  file_name <-  paste (file_path, "09_HSBC_", sep="")
} else if (n_y == 10){
  file_name <-  paste (file_path, "10_JPM_", sep="")
} else if (n_y == 11){
  file_name <-  paste (file_path, "11_MS_", sep="")
} else if (n_y == 12){
  file_name <-  paste (file_path, "12_RBS_", sep="")
} else if (n_y == 13){
  file_name <-  paste (file_path, "13_SocGen_", sep="")
} else if (n_y == 14){
  file_name <-  paste (file_path, "14_UBS_", sep="")
}
file_name_lasso <- paste( file_name,"covar05_adapt.csv",sep="")
write.csv2(COVAR1[[(n_y-1),1]], file = file_name_lasso)
}
