Step1<-function(x)
{
  if(!is.data.frame(x)) stop("The input data is not a data frame object") 
  x_list<-vector("list",length = nrow(x))
  for(i in 1:nrow(x))
  {
    x_list[[i]]<-x[,order(x[i,],decreasing=TRUE)] 
  }
  return(x_list)
}

Step2<-function(x,R,time_id) 
{
  for (j in 1:(ncol(x)-1)) 
  {
    cs<-ConvergenceSpeed(x[,j:(j+1)],R,time_id)
    if (cs>=0) 
    {
      k<-1
      core_group_indexes<-integer(2)
      core_group_indexes[1]<-j
      core_group_indexes[2]<-j+1
      while (cs>=0 & (j+1+k)<=ncol(x))
      {
        cs<-ConvergenceSpeed(x[,j:(j+1+k)],R,time_id)
        if (cs>=0) core_group_indexes[2]<-core_group_indexes[2]+1
        k<-k+1
      }
      break
    }
    if (cs<0 & j==(ncol(x)-1)) 
    {
      return(cs) 
    }
  }
  
  club<-logical(ncol(x)) 
  club[core_group_indexes[1]:core_group_indexes[2]]<-TRUE
  if (all(club)) 
  {
    attr(club,"units")<-names(x)[club]
    fcs<-ConvergenceSpeed(x[,club],R,time_id,boot=FALSE)
    attr(club,"cs")<-fcs
    return(club)
  }
  
  club1<-club 
  for (i in 1:length(club))
  {
    if (!club[i]==TRUE)
    {
      club[i]<-TRUE
      cs<-ConvergenceSpeed(x[,club],R,time_id)
      if (cs>=0) club1[i]<-TRUE 
      cs1<-ConvergenceSpeed(x[,club1],R,time_id) 
      if (cs1<0) club1[i]<-FALSE
      club[i]<-FALSE 
    }
  }
  attr(club1,"units")<-names(x)[club1]
  fcs<-ConvergenceSpeed(x[,club1],R,time_id,boot=FALSE)
  attr(club1,"cs")<-fcs
  return(club1)
}

ConvergenceSpeed<-function(x,R,time_id,boot=TRUE) 
{
  library(locfit)
  a<-0.4
  col_n<-ncol(x)
  x<-cbind(x, hdem = rowSums(x)/ncol(x))
  for (i in 1:(ncol(x)-1))
  {
    x<-cbind(x, h=((x[,i]/x$hdem)-1)^2)
  }
  x<-cbind(x, ht=rowSums(x[,(col_n+2):(2*col_n+1)])/col_n)
  X<-seq(from=R+1,to=nrow(x))
  Y<-log(x[1,ncol(x)]/x[(R+1):nrow(x),ncol(x)])-2*log(log(X))
  
  if (boot==FALSE)
  {
    model<-locfit(Y~log(X),alpha=a,kern="tcub",deg=1,deriv=1)
    local_slopes<-fitted(model)
    return(local_slopes[time_id])
  }
  else
  {
  
  model<-locfit(Y~log(X),alpha=a,kern="tcub",deg=1,deriv=1)
  local_slopes<-fitted(model)
  model_b<-locfit(Y~log(X),alpha=a,kern="tcub",deg=1)
  y_fit<-fitted(model_b) 
  l<-round(length(X)^(1/3))
  res<-residuals(model_b)
  res_blocks<-matrix(nrow =l,ncol=length(X)) 
  for (i in 1:length(X)) 
  {
    if (i<=(length(X)-l+1)) 
    {
      res_blocks[,i]<-res[i:(i+l-1)]
    }
    else
    {
      diff<-l-length(i:length(X))
      res_blocks[,i]<-res[c(i:length(X),1:diff)]
    }
  }
  local_slope_time<-numeric(1000)
  for (i in 1:1000) #ustawić liczbę replikacji
  {
    if ((length(X)%%l)==0) s<-sample(1:(length(X)/l),replace=TRUE)
    if ((length(X)%%l)!=0) s<-sample(1:((length(X)/l)+1),replace=TRUE)
    res_new<-c(res_blocks[,s])
    res_new<-res_new[1:length(X)]
    Y_new<-y_fit+res_new
    model_bot<-locfit(Y_new~log(X),alpha=a,kern="tcub",deg=1,deriv=1)
    local_slopes_b<-fitted(model_bot)
    local_slope_time[i]<-local_slopes_b[time_id]
  }
  local_slope_time<-sort(local_slope_time)
  return(local_slope_time[25])
  }
}

DynamicConvergenceClubs<-function(x,end_date,fr,trend=TRUE)
{
  if (trend==TRUE) 
  {
    library(bHP)
    y<-lapply(x,ts)
    y<-lapply(y,BoostedHP,lambda=6.25)
    for (i in 1:length(y)) 
    {
      y[[i]]<-y[[i]]$trend
    }
    y<-lapply(y,as.vector)
    y<-do.call(cbind, lapply(y, as.data.frame))
    names(y)<-names(x)
    x<-y
  }
  unit_names<-names(x)
  t<-nrow(x)
  if (t<=50) r<-0.3
  if (t>=100) r<-0.2
  if ((t<100)&(t>50)) r<-0.3-((t-50)*(0.1/50))
  R<-round(r*t) 
  
  x<-Step1(x) 
  clubs_number<-integer(nrow(x[[1]])-R)
  FinalClubs<-vector(mode="list",length =(nrow(x[[1]])-R))
  ConSpeedCoeff<-numeric()
  for (i in 1:(nrow(x[[1]])-R))
  {
    y<-x[[i+R]]
    while (ncol(y)>1)
    {
      club<-Step2(y,R,i) 
      if (is.logical(club)) 
      {
        if (is.null(FinalClubs[[i]]))
        {
          FinalClubs[i]<-list(as.character(attr(club,"units")))
          ConSpeedCoeff<-c(ConSpeedCoeff,attr(club,"cs"))
          print(ConSpeedCoeff)
        }
        else
        {
          FinalClubs[i]<-list(c(FinalClubs[i],list(as.character(attr(club,"units")))))
          ConSpeedCoeff<-c(ConSpeedCoeff,attr(club,"cs"))
          print(ConSpeedCoeff)
        }
        clubs_number[i]<-clubs_number[i]+1
        y<-y[,!club]
        if (is.null(ncol(y))) break
      }
      else break
    }
  }
 

  for (i in 1:(nrow(x[[1]])-R))
  {
    FinalClubs[[i]]<-lapply(rapply(FinalClubs[i], enquote, how="unlist"), eval)
  }
  
  DivergentUnit<-list(nrow(x[[1]])-R)
  for (i in 1:(nrow(x[[1]])-R)) 
  {
    DivergentUnit[[i]]<-c(setdiff(unlist(FinalClubs[[i]]),unit_names), setdiff(unit_names,unlist(FinalClubs[[i]])))
  }
  
  clubs_number<-ts(clubs_number,frequency = fr,end = end_date)
  results<-list(ClubNumber=clubs_number,ClubStructure=FinalClubs,DivergentUnit=DivergentUnit,UnitNames=unit_names,TimeStamp=nrow(x[[1]])-R,ConvergenceCoefficients=ConSpeedCoeff)
  attr(results,"dcc")<-TRUE
  return(results)
  
}

ClubSimilarity<-function(x,end_date,fr) 
{
  library(zoo)
  if(is.null(attr(x,"dcc"))) stop ("This is not an object from the DynamicConvergenceClubs function")
  Similarity<-matrix(nrow=x[[5]],ncol=x[[5]])
  D<-numeric(x[[5]])
  P<-numeric(x[[5]]) 
  
  for (i in 1:x[[5]])
  {
    for (j in 1:(length(x[[2]][[i]]))) 
    {
        D[i]<-D[i]+(length(x[[2]][[i]][[j]])*(length(x[[2]][[i]][[j]])-1))
    } 
    P[i]<-sqrt(D[i]/(length(x[[4]])*(length(x[[4]])-1)))
  }
  
   for (i in 1:x[[5]]) 
   {
     for (j in 1:x[[5]]) 
     {
       if (i==j) 
       {
         Similarity[i,j]<-1
         break
       }
       tmp<-0
       d<-0
       for (k in 1:length(x[[2]][[i]]))
       {
         for (l in 1:length(x[[2]][[j]]))
          {
            tmp<-length(intersect(x[[2]][[i]][[k]],x[[2]][[j]][[l]]))
            tmp<-tmp*(tmp-1)
            d<-tmp+d
          }
       }
       Similarity[i,j]<-sqrt(d/(sqrt(D[i])*sqrt(D[j])))
       Similarity[j,i]<-sqrt(d/(sqrt(D[i])*sqrt(D[j])))
     }
   }
  P<-ts(P,frequency = fr,end = end_date)
  if (fr==1)
  {
    colnames(Similarity)<-as.character(time(P))
    rownames(Similarity)<-as.character(time(P)) 
  }
  if (fr==4)
  {
    colnames(Similarity)<-as.character(as.yearqtr(time(P))) 
    rownames(Similarity)<-as.character(as.yearqtr(time(P))) 
  }
  if (fr==12)
  {
    colnames(Similarity)<-as.character(as.yearmon(time(P))) 
    rownames(Similarity)<-as.character(as.yearmon(time(P)))
  }
  Results<-list(ConvergenceDegree=P,ClubSimilarity=Similarity)
  return(Results)
}


#example
library(openxlsx)
library(locfit)
library(bHP)
library(zoo)
library(corrplot)
x<-read.xlsx("C:/Users/Mateusz/Documents/DynamicConvergenceClubs/data.xlsx")
model<-DynamicConvergenceClubs(x,end_date = 2021,fr=1)
similarity<-ClubSimilarity(model,end_date = 2021,fr=1)
corrplot(similarity[[2]],type="lower",tl.col="black",method="shade", number.cex=0.5,cl.cex=1, addCoef.col="white")