load("r250")
load("Beta250")
load("SectorTable")

source('RFstepArchetypesModB.R')
source('RFFunctions_to_calculate_real_archetypes_with_swap2B.R')
source('RFstepArchetypoidsB.R')
source("auxiliaryfunctions.R")

library(dplyr)
library(data.table)
library(fda)
library(network)
library(GGally)
library(RColorBrewer)
library(colorspace)
library(colorRamps)
library(ggplot2)
library(zoo)
library(scales)


#### Remove with missings > 20%

u<-0.2
r250_filtrado<-del.Umbr.nas(r250,u)
Beta250_filtrado<-del.Umbr.nas(Beta250,u)


# Generate [BI]variate [F]unctional [D]ata object.
FD.Bivariable<-gen.FD.object.BI(N = 13,Betas = Beta250_filtrado,retornos = r250_filtrado)

#

extraer.coefs.norm<-function(tempfd){
  X=array(0,dim=c(dim(t(tempfd$coef[,,1])),2))
  X[,,1]=t(tempfd$coef[,,1])
  X[,,2]=t(tempfd$coef[,,2])
  
  Xs=X
  Xs[,,1]=scale(X[,,1])
  Xs[,,2]=scale(X[,,2])
  
  return(Xs)
}
Xs<-extraer.coefs.norm(FD.Bivariable)

x11 <- t(Xs[,,1]) 
x12 <- t(Xs[,,2])
x1=rbind(x11,x12)   

data=t(x1)      



### Calc achetypes with K=4 ###

dim(Xs)
data.length<-nrow(r250_filtrado)
fbasis<-create.bspline.basis(rangeval = c(1,data.length),nbasis = 13)
PM<-eval.penalty(fbasis)


set.seed(1)
norep=10

K<-4

lass10 <- RFstepArchetypesModB(data=Xs,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F,PO=0.5,CT=6) 
ai<- bestModel(lass10[[1]]) #archetypes
ad=stepRFArchetypoids3B(data=Xs, k=K, norep=norep,PM,PO=0.5,CT=6)
K_5<-ad[[1]]
Aqt<-ad[[1]][[1]]
Aqt.names<-colnames(FD.Bivariable$coefs)[Aqt]


#Aqt
# 251 24 301 180
#243  23 292 174

#####Calc alphas

huge=200
options<-K_5
K<-4
# errors<-unlist(lapply(options,FUN="[[",2))
# pos.min.error<-min(which(errors==min(errors))) #if ties: min in is used.

archs_ini <- K_5[[1]]
# data<-data

n <- ncol(t(data))
x_gvv <- rbind(t(data), rep(huge, n))

zs <- x_gvv[, archs_ini] 
zs <- as.matrix(zs)

alphas <- matrix(0, nrow = K, ncol = n)
for (j in 1 : n){
  alphas[, j] = coef(nnls(zs, x_gvv[,j]))
}

alphas_K5<-alphas

aqtp_index_K5<-K_5[[1]]
stock_names<-toupper(sub("\\w+_","",names(r250_filtrado)))
colnames(alphas_K5)<-stock_names
rownames(alphas_K5)<-stock_names[aqtp_index_K5]



####################### Graphic outputs ##########################

##################
#     Figure 1
#################


n<-Aqt
evaluated<-eval.fd(1:3422,FD.Bivariable[n])
archs_rets<-evaluated[,,1]
archs_betas<-evaluated[,,2]      

names<-names(r250_filtrado)[Aqt]
names<-toupper(gsub("\\w+_","",names))
colnames(evaluated)<-names
rownames(evaluated)<-rownames(Beta250_filtrado)
mdf<-melt(evaluated)
colnames(mdf)<-c("Date","stock","variable","value")
mdf$Date<-as.Date(mdf$Date)
mdf$variable[mdf$variable==1]<-"r250"
mdf$variable[mdf$variable==2]<-"Beta250"
mdf$var_f<-factor(mdf$variable,levels = c("r250","Beta250"))


### Overlapped archetypoids ###
p6<-ggplot(mdf,aes(x=Date,y=value,group=stock,colour=stock,linetype=stock))+geom_line()+facet_grid(var_f~.)+ggtitle("Functional Archetypoids")

p6<-ggplot(mdf,aes(x=Date,y=value,group=stock,colour=stock,linetype=stock))+geom_line()+facet_grid(var_f~.)

p6


##################
#     Figure 2
#################


####################### Vary  'umbral' for different networks.

K5.DT<-as.data.table(t(alphas_K5))
K5.DT<-cbind(ID=1:440,stock=colnames(alphas_K5),K5.DT)



umbral<-0.8 #Here we set the U for the network

n_stocks<-ncol(r250_filtrado)
netMat<-matrix(0, nrow = n_stocks, ncol = n_stocks)
colnames(netMat)<-K5.DT$stock
rownames(netMat)<-colnames(netMat)
K5.melt<-melt(K5.DT,id = c("stock","ID"),variable.name = "Arquetipoid")
K5.aqtpd.names<-colnames(K5.DT)[c(-1,-2)]

#Pure clusters

CP<-vector("list",length = K)
names(CP)<-K5.aqtpd.names
for (i in 1:K){
  
  CP[[i]]<-K5.melt[value>umbral & Arquetipoid==K5.aqtpd.names[i]][,ID]
}


#Mixture of 2

num_aqtp<-length(K5.aqtpd.names)
combs<-(num_aqtp^2-num_aqtp)/2
names_CM<-vector(length = combs)
k=1
for (i in 1:num_aqtp){
  for(j in (i+1):num_aqtp){
    if(i==4&j==5) break
    
    names_CM[[k]]<-paste(K5.aqtpd.names[i],K5.aqtpd.names[j],sep="-")
    k=k+1
  }
}
#The same but with codes
num_aqtp<-length(K5.aqtpd.names)
combs<-(num_aqtp^2-num_aqtp)/2
index_CM<-vector(length = combs)
k=1
for (i in 1:num_aqtp){
  for(j in (i+1):num_aqtp){
    if(i==5&j==6) break
    
    index_CM[[k]]<-paste(aqtp_index_K5[i],aqtp_index_K5[j],sep="-")
    k=k+1
  }
}

A1<-strsplit(names_CM[1],"-")[[1]][[1]]
A2<-strsplit(names_CM[1],"-")[[1]][[2]]
CM<-vector("list",length = combs)
for(i in 1:combs){
  
  A1<-strsplit(names_CM[i],"-")[[1]][[1]]
  A2<-strsplit(names_CM[i],"-")[[1]][[2]]
  
  CM[[i]]<-K5.melt[Arquetipoid%in%c(A1,A2)&value<umbral,.(suma=sum(value)),by=.(ID,stock)][suma>umbral,ID]
  
}
names(CM)<-names_CM

#####Adj Matrix

n_stocks<-ncol(r250_filtrado)
netMat<-matrix(0, nrow = n_stocks, ncol = n_stocks)
colnames(netMat)<-K5.DT$stock
rownames(netMat)<-colnames(netMat)
aqtp_index_K5
my_vector<-CP
names(my_vector)<-aqtp_index_K5
as.numeric(names(my_vector[1]))
vector_tipos<-vector(length = nrow(netMat))

#For pure
L<-length(aqtp_index_K5)
for( i in 1:L){
  aqtpd<-aqtp_index_K5[i]
  friends<-CP[[i]]
  netMat[aqtpd,friends]<-1
  vector_tipos[friends]<-names(CP[i])
}
#For mixture
L<-length(CM)

for (i in 1:L){
  aqtpds<-as.numeric(unlist(strsplit(index_CM[i],"-")))
  friends<-CM[[i]]
  netMat[aqtpds,friends]<-1
  
  vector_tipos[friends]<-names(CM[i])
  
}




unconected<-which(apply(netMat,2,sum)==0)


net<-network(netMat[-unconected,-unconected],directed = F)


net %v% "tipo"<- vector_tipos[-unconected]

nam<-colnames(netMat)

sectores<-sapply(nam,function(x) as.character(Sector.Table[Sector.Table$Symbol==x,"Sector"]))

net %v%"sector"=sectores[-unconected]





cols<-c("black","darkgreen","gold3","navy","green","turquoise3","yellow","greenyellow","slategrey","blue","red","darkred","yellow3","orange","mediumorchid","white")
y<-colorRampPalette(cols,bias = 1)(16)

names(y) = levels(as.factor(vector_tipos))



set.seed(1234)
net_plot<-ggnet2(net,size=6,label.size = 4,label=T,color="tipo",shape="sector",label.color="darkblue",palette = y, edge.color = c("color","grey50"),edge.alpha = 0.2,mode = "kamadakawai", layout.par = list(niter = 10000),legend.position="bottom",legend.size = 8)+scale_shape_manual(values=c(1:11)) +guides(col = guide_legend(nrow = 5))



plot2<-net_plot+geom_point(data=subset(net_plot$data,label%in%c("EBAY","LNC","ED","NTAP","EOG")),size=8,alpha=0.9,color="grey",shape=15)+geom_point(size=4,alpha=0.5,aes(color=color))

# net_plot$data

plot2$layers<-plot2$layers[c(5,4,1,2,3)]
plot2$labels$shape<-"Sector"


plot2
#uncoment to save plot
#ggsave("red80.pdf",width = 60,height = 30, units = "cm")
#ggsave("red80.pdf",width = 594,height = 420, units = "mm") # A4


##################
#     Figure 3
#################

check_sector<-function(x){
  return(as.character(Sector.Table[Sector.Table$Symbol==x,"Sector"]))
}


stock_names<-toupper(sub("\\w+_","",names(r250_filtrado)))

M1<-t(alphas_K5)

secs<-sapply(rownames(M1),check_sector)

M1<-as.data.table(M1)
M1<-cbind(M1,secs)
meltM1<-melt(M1,id.vars="secs")
meltM1[,secs:= sub("\\w+_","",secs)]
meltM1[,secs:= sub("character\\(0\\)","-",secs)]
sumado<-meltM1[,sum(value),by=.(variable,secs)]

sumado[secs!="-",]
str(sumado)
suma_sec<-sumado[,sum(V1),by=secs]


vector<-suma_sec[match(sumado$secs,suma_sec$secs),V1]
sumado[,tot_sec:=vector]

tabla_pesos<-sumado[,.(V1/tot_sec,by=secs,variable)]

tabla_pesos[order(by)]


ggplot(sumado[secs!="-",],aes(x = secs, y = V1,fill = variable)) + 
  geom_bar(position = "fill",stat = "identity",colour="black") +
  # or:
  # geom_bar(position = position_fill(), stat = "identity") 
  scale_y_continuous(labels = percent_format()) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  ggtitle("Sector composition")+ylab("weight")+xlab("Sectors") +
  scale_fill_grey()






