######################### Librerias #########################

library(fda)
library(ggplot2)
library(MASS)
library(dplyr)
library(purrr)
library(ggpubr)
library(tidyr)
library(funData)
library(MFPCA)
library(cluster)
library(factoextra)
library(imputeTS)
library(readxl)

old.dir=getwd() 



######################### lectura de datos #########################



Indice_tv <- read_excel("Indice_tv.xlsx")
II=Indice_tv 
II=II[-c(1:4),-1]
IT=apply(as.matrix(II), 2, as.numeric);IT=round(IT*100,2)

sum(is.na(IT))

IT=na_kalman(IT,model="auto.arima")




######################### Se suavisan las funciones #########################

BSpl     <- create.bspline.basis(norder=3, breaks=seq(0,40,length=3))
Ftemp    <- Data2fd(IT, basisobj=BSpl)
plot(Ftemp, main="Datos Suavizados usando la funcion Data2fd con una base de BSpl") 



######################### Se crean los datos funcionales #########################

dir.create("Colombia")
setwd("Colombia")

fundD = fd2funData(Ftemp,argvals = seq(0,40,length.out = 100))

png("Curvas.png",width=1200,height=820)
autoplot.funData(fundD,obs = 1:20) +
  geom_line(aes(colour = obs),show.legend = FALSE)
dev.off()



######################### Clasificacion a nivel Colombia #########################


km.res <- kmeans(t(Ftemp$coefs), 3 )

png("Clusters.png",width=1200,height=820)
fviz_cluster(km.res, data = t(Ftemp$coefs), frame.type = "convex")
dev.off()


d <- dist(t(Ftemp$coefs), method = "euclidean")
res.hc <- hclust(d, method = "ward.D2")

png("Dendograma.png",width=1200,height=820)
plot(res.hc,cex = 0.8)
rect.hclust(res.hc, k = 4, border = 2:5)
dev.off()


######################### Clasificacion a nivel Departamento #########################


BSpl  <- create.bspline.basis(norder=3, breaks=seq(0,40,length=3))
total <- read.delim("C:/Users/57313/OneDrive/Escritorio/Data Jamp/Television/departamentos.txt", header=FALSE)
total <- total[-3:-6,]


deps=unique(t(total[1,]))[-1]
Ndeps=length(deps)

 

i=33



setwd(old.dir)
dir.create(deps[i])
setwd(deps[i])

Departamento=total[,total[1,]==deps[i]][-3:-6,]
colnames(Departamento)=Departamento[2,]
Departamento=apply(Departamento[-1:-2,],2,as.numeric)

sum(is.na(Departamento))
Departamento=na_kalman(Departamento,model="auto.arima")


Fdep<- Data2fd(Departamento, basisobj=BSpl)
fundDep = fd2funData(Fdep,argvals = seq(0,40,length.out = 100))

{png("Curvas.png",width=1200,height=820)
  autoplot.funData(fundDep,obs = 1:20) +
    geom_line(aes(colour = obs),show.legend = FALSE)}
dev.off()

km.resdep <- kmeans(t(Fdep$coefs), 3 )


{png("Clusters.png",width=1200,height=820)
  fviz_cluster(km.resdep, data = t(Fdep$coefs), frame.type = "convex")}
dev.off()


ddep <- dist(t(Fdep$coefs), method = "euclidean")
res.hcdep <- hclust(ddep, method = "ward.D2")

{png("Dendograma.png",width=1200,height=820)
  plot(res.hcdep,cex = 0.8)
  rect.hclust(res.hcdep, k = 3, border = 2:5)}
dev.off()





