library(devtools)
library(magrittr)
#require(rgdal)
library(sf)


pkg <- "/Users/andresvallone/Dropbox/tesis_AndresVallone/PAPERS/Rpackage/03-estdaR/estdaR"

load_all(pkg)

Year<-c(1930,1940,1952,1960,1970,1982,1992,2002)
per <- c("1930-1940","1940-1952","1952-1960","1960-1970","1970-1982","1982-1992","1992-2002")

set.seed(2010)


#base<-readOGR("./SHP","c5000")
base <- st_read(dsn="./SHP/c5000.shp",quiet=TRUE)


base$P70[180] <- 1525 #fix data problems in La islita 

b <- as.data.frame(base)
b <- as.data.frame(b[,17:24])
b <- apply(b,2,as.numeric)
colnames(b)<-Year

#Regime def.

Regime <- rep("South", dim(base)[1])

Regime[which(base$FIELD2=="I"| base$FIELD2=="XV"| base$FIELD2=="II" | base$FIELD2=="III" | base$FIELD2 =="IV")] <- "North"

Regime[which(base$FIELD2 == "RM" | base$FIELD2 =="V" | base$FIELD2 =="VI" | base$FIELD2 =="VII"  | base$FIELD2 =="VIII")] <- "Center"

###################################
#Kernel
###################################


den<-function(x){ #log calc and mean normalization
	y<-log(x)
	z<-mean(y)
	out<-y/z
	out
}

D1 <- apply(b,2,den)

D <- D1

geom.text.size = 7
theme.size = (14/5) * geom.text.size

#full sample
a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig3.a <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig3.a

#Center

D <- D1[Regime=="Center",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig3.b <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig3.b

#North

D <- D1[Regime=="North",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig3.c <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig3.c

#South

#RM

D <- D1[Regime=="South",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig3.d <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig3.d


#Markov

full.markov <- mkv(b)
full.markov[[1]] %>% round(4)
st.st(full.markov[[1]])

#sp markov

## W matrix def
#coords<-coordinates(base)
coords <- st_coordinates(base)
k1 <- knn2nb(knearneigh(coords, k=1))
k1d <- nbdists(k1, coords)
max_k1 <- max(unlist(k1d))
dnn_nb <- dnearneigh(coords, 0, max_k1)
dlist<-nbdists(dnn_nb,coords)
idlist<-lapply(dlist, function(x) 1/x)
Wd<-nb2listw(dnn_nb,glist=idlist,style="W")

spatial.markov <- sp.mkv(b,Wd,fixed=FALSE)

sapply(1:5,function(x)st.st(spatial.markov[[1]][,,x]))
lapply(1:5,function(x)prais(spatial.markov[[1]][,,x]))

sp.homo.test(b,Wd,fixed=FALSE)

#LISA Markov

L.M <- lisamkv(b,Wd,geoda=FALSE)

L.M[[3]] %>% round(4)
st.st(L.M$p.lisamatrix)

join.d(b,Wd)

#Directional LISA



full.p8 <- d.LISA(D1[,1],D1[,8],Wd,Regime=Regime,nsim=999)
full.p4 <- d.LISA(D1[,1],D1[,8],Wd,Regime=Regime,nsim=999,k=4)


# k= 8

step <- (2*pi)/8
breaks <- seq(0,2*pi,step)
lmt <- (breaks * 180)/pi # rad to deg
symb<-c()
for (i in seq_len(length(lmt)-1)){
	symb <- c(symb, paste(lmt[i],lmt[i+1],sep="-"))
}
lim=c(symb[1],symb[length(symb):2])

dL <- lapply(seq_len((ncol(b)-1)),function(x)d.LISA(D[,x],D[,x+1],Wd,Regime=Regime,nsim=999))

p1 <- list()
for(i in 1:4){p1[[i]]<-dL[[i]]$data} 

p2 <- list()
for(i in 5:7){p2[[i]]<-dL[[i]]$data}

pt <- list()
for(i in 1:7){pt[[i]]<-dL[[i]]$data}


p1 <- do.call(rbind.data.frame,p1)
p2 <- do.call(rbind.data.frame,p2)
pt <- do.call(rbind.data.frame,pt)


# total
quartz()
rose.pt<-ggplot(pt,aes(angle,fill=Regime),xlab=" ",ylab=" ") +geom_bar(width=1,colour="black",size=0.1)+scale_x_discrete(name="",limit=lim)+scale_y_continuous(name="",breaks=seq(0,max(table(pt$angle)),5),labels=seq(0,max(table(pt$angle)),5))+coord_polar(start=pi/4)+theme(legend.position="bottom",panel.grid.major = element_line( color="gray50",linetype="solid"),panel.background = element_blank())+scale_fill_brewer(palette="Set1",direction=-1)
rose.pt




# Rank Decomposition

Reg <- as.numeric(as.factor(Regime))

th <- theta(b,Reg,nsim=999)

theta.table <- data.frame("Period"=per,"Theta"=th[,1],"p-value"=th[,2])

#GIMA

#block matrix

WR <- matrix(0,nrow(base),nrow(base))
for(i in seq_len(nrow(base))){
	WR[,i] <- 1*(Regime[i]==Regime)
}
diag(WR) <- 0

WR <- mat2listw(WR)

spTau <- lapply(1:(ncol(b)-1),function(i) sp.tau(b[,i],b[,(i+1)],WR,999))

GIMA <- do.call(rbind.data.frame,spTau)

colnames(GIMA)<-names(spTau[[1]])


p <- matrix(1,nrow(b),nrow(b))
diag(p)<-0

phi <- spweights.constants(WR)$S0/sum(p)

nonng.tau <- (GIMA[,5]-(phi*GIMA[,1]))/(1-phi)
ng.tau <- (GIMA[,1])
ng.M <- (ng.tau-1)/-2
nonng.M <- (nonng.tau-1)/-2
pvalue <- GIMA[,2]

tabla.GIMA <- data.frame("Period"=per,"MW"=ng.M,"MnonW"=nonng.M,"TW"=ng.tau,"TnonW"=nonng.tau,"p-value"=pvalue)

spTau_1 <- lapply(1:(ncol(b)-1),function(i) sp.tau(b[,i],b[,(i+1)],Wd,999))
GIMA_1 <- do.call(rbind.data.frame,spTau_1)

colnames(GIMA_1)<-names(spTau_1[[1]])

p <- matrix(1,nrow(b),nrow(b))
diag(p)<-0

phi <- spweights.constants(Wd)$S0/sum(p)

nonng.tau_1 <- (GIMA_1[,5]-(phi*GIMA_1[,1]))/(1-phi)
ng.tau_1 <- (GIMA_1[,1])
ng.M_1 <- (ng.tau_1-1)/-2
nonng.M_1 <- (nonng.tau_1-1)/-2
pvalue_1 <- GIMA_1[,2]

tabla.GIMA_1 <- data.frame("Period"=per,"MW"=ng.M_1,"MnonW"=nonng.M_1,"TW"=ng.tau_1,"TnonW"=nonng.tau_1,"p-value"=pvalue_1)



################################
#OTRA BASE
################################

base1<-st_read("./FMLA data/FLMAs.shp",quiet=TRUE)
data <- as.data.frame(base1)

corr <- openxlsx::read.xlsx("./FMLA data/Correspondence Mun to LMA.xlsx")

nn <- sapply(base$NAME,as.character)

fmla_cod <- rep(999, length(nn))

for (i in seq_along(nn)){
	a <- which(corr$Name_Mun == nn[i])
	if(length(a) == 0){
		fmla_cod[i] <- 0
	} else {
		fmla_cod[i] <- corr[a,3]
	}
}

fmla_cod[4] <- "5108"
fmla_cod[11] <- "3202"
fmla_cod[13] <- "13101"
fmla_cod[14] <- "7301"
fmla_cod[15] <- "13101"
fmla_cod[17] <- "6112"
fmla_cod[20] <- "5108" #Puchuncavi
fmla_cod[22] <- "6101"
fmla_cod[23] <- "13101"
fmla_cod[24] <- "8301"
fmla_cod[26] <- "2101" #Error de nombre en la base es Antofagasta
fmla_cod[28] <- "6101"
fmla_cod[30] <- "3202"
fmla_cod[31] <- "8401"
fmla_cod[33] <- "9210"
fmla_cod[34] <- "13101" #Lampa
fmla_cod[39] <- "4301"
fmla_cod[43] <- "5501"
fmla_cod[44] <- "10208"
fmla_cod[46] <- "8416" #Tucapel
fmla_cod[51] <- "4301" # Ovalle por cercania, no aparece en la lista a pesar de ser un muncipio importante de la IV
fmla_cod[55] <- "7103"
fmla_cod[58] <- "6101"
fmla_cod[61] <- "8203"
fmla_cod[64] <- "13101" #Buin
fmla_cod[65] <- "4301" #Está en Monte Paria que no está caputrado en la lista de municipios.
fmla_cod[67] <- "8301"
fmla_cod[79] <- "5701"
fmla_cod[84] <- "8401" #Yumbel por cercania
fmla_cod[85] <- "9117"
fmla_cod[88] <- "8101"
fmla_cod[89] <- "5108" #Region metropolitana de Valparaiso, no lo considera en el listado
fmla_cod[90] <- "9101"
fmla_cod[93] <- "10510"
fmla_cod[95] <- "2201" #Calama, campamento qu ese cierra en 2007 pero activo en el2002.
fmla_cod[96] <- "9210"
fmla_cod[99] <- "13101" #Calera de Tango
fmla_cod[104] <- "13101"
fmla_cod[105] <- "11101"
fmla_cod[108] <- "5108"
fmla_cod[110] <- "5501" #Nogales
fmla_cod[111] <- "10504"
fmla_cod[113] <- "8301"
fmla_cod[116] <- "8101" #Santa Juana Coronel por cercania
fmla_cod[117] <- "13101" #Paine
fmla_cod[121] <- "7401"
fmla_cod[130] <- "8401"  #Monte aguila, Yumbel por cercania.
fmla_cod[132] <- "2301"
fmla_cod[134] <- "3101"
fmla_cod[136] <- "4101" #vicuña, Coquimbo por cercania
fmla_cod[137] <- "10301"
fmla_cod[139] <- "5108"
fmla_cod[143] <- "5108"
fmla_cod[144] <- "8401" #Uso San Carlos, pero no es lo mismo.
fmla_cod[149] <- "10504" # La Union por cercania
fmla_cod[150] <- "8401" #Yumbel por cercania
fmla_cod[153] <- "8101"
fmla_cod[154] <- "9203"
fmla_cod[157] <- "7102"
fmla_cod[161] <- "12401" #Natales
fmla_cod[149] <- "10205"
fmla_cod[165] <- "6101" #comparte el mismo distrito electoral de Doñihue por eso el mismo codigo
fmla_cod[166] <- "6101" # Idem 165
fmla_cod[167] <- "6101" # Iden a 165
fmla_cod[168] <- "6307" #Por cercania a Placilla
fmla_cod[171] <- "6101" #Conurbacion Rancagua
fmla_cod[175] <- "13101" 
fmla_cod[180] <- "13101"
fmla_cod[181] <- "5108"
fmla_cod[183] <- "8301" #Uso Los Ageles por cercania, pero no es lo mismo.
fmla_cod[184] <- "10504"



##col def

n <- fmla_cod %>% unique %>% length
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))

library(dplyr)
data_fmla<-cbind(b,fmla_cod) %>% as.data.frame
data_fmla[,1:8] <- apply(data_fmla[,1:8],2,as.character)
data_fmla[,1:8] <- apply(data_fmla[,1:8],2,as.numeric)
fmla_data <-  data_fmla %>% group_by(fmla_cod) %>% summarise_all(funs(sum))

##Agregando datos
peralillo <- cbind("6304",721,970,1512,2064,2147,2570,3884,4439) %>% as.data.frame
chaiten <-  cbind("10401",300,461,480,	663,1375,2599,3258,4065) %>%as.data.frame #no hay datos reales de la poblacion en 1930, se usa una estimació historica.
c.names <- colnames(fmla_data)
colnames(peralillo) <- c.names
colnames(chaiten) <- c.names
fmla_data <- rbind(fmla_data,peralillo)
fmla_data <- rbind(fmla_data,chaiten)
fmla_data <- apply(fmla_data,2,as.character)
fmla_data <- apply(fmla_data,2,as.numeric)
fmla_data %<>% as.data.frame
colnames(fmla_data) <- c.names
colnames(fmla_data)[1] <- "code_flma"
fmla_data$code_fmla %>% unique -> u1
corr[,3] %>% unique -> u2
u2[!u2 %in% u1] -> f1

# 5103, 5201 y 12202 son islas alejadas del continente. 


# Construyendo la matriz W

id_base <- base1[,140]
data_f <- inner_join(id_base,fmla_data,by="code_flma")

cent <- data_f %>% st_geometry %>% st_centroid

coords1 <- st_coordinates(cent)
k1 <- knn2nb(knearneigh(coords1, k=1))
k1d <- nbdists(k1, coords1)
max_k1 <- max(unlist(k1d))
dnn_nb <- dnearneigh(coords1, 0, max_k1)
dlist<-nbdists(dnn_nb,coords1)
idlist<-lapply(dlist, function(x) 1/x)
Wd1<-nb2listw(dnn_nb,glist=idlist,style="W")

#Regime def.

Regime1 <- rep("Center",dim(fmla_data)[1])
Regime1[which( fmla_data$code_flma==9101 | fmla_data$code_flma==9102 | fmla_data$code_flma==9117 |  fmla_data$code_flma==9201 | fmla_data$code_flma==9203 | fmla_data$code_flma==9210 | fmla_data$code_flma==9211 | fmla_data$code_flma==10101 | fmla_data$code_flma==10201 | fmla_data$code_flma==10202 | fmla_data$code_flma==10205 | fmla_data$code_flma==10208 |fmla_data$code_flma==10301 | fmla_data$code_flma==10303 | fmla_data$code_flma==10401 |fmla_data$code_flma==10501 |fmla_data$code_flma==10504 |fmla_data$code_flma==10510 |fmla_data$code_flma==11101 |fmla_data$code_flma==12101 |fmla_data$code_flma==12401)] <- "South"
Regime1[which(fmla_data$code_flma==1101 | fmla_data$code_flma==1201 | fmla_data$code_flma==2101 |  fmla_data$code_flma==2201 | fmla_data$code_flma==2301 | fmla_data$code_flma==3101 | fmla_data$code_flma==3202 | fmla_data$code_flma==3301 | fmla_data$code_flma==4101 | fmla_data$code_flma==4201 | fmla_data$code_flma==4301)] <- "North"


#################
#Kernel
#################

D2 <- apply(fmla_data[,2:9],2,den)

D <- D2

geom.text.size = 7
theme.size = (14/5) * geom.text.size

#full sample
a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig4.a <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig4.a

#Center

D <- D2[Regime1=="Center",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig4.b <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig4.b

#North

D <- D2[Regime1=="North",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig4.c <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig4.c

#South


D <- D2[Regime1=="South",]

a30 <- density(D[,1])
x30 <- a30$x
y30 <- a30$y
l30 <- rep("1930",length(x30))

a52 <- density(D[,3])
x52 <- a52$x
y52 <- a52$y
l52 <- rep("1952",length(x52))

a70 <- density(D[,6])
x70 <- a70$x
y70 <- a70$y
l70 <- rep("1970",length(x70))

a02 <- density(D[,8])
x02 <- a02$x
y02 <- a02$y
l02 <- rep("2002",length(x02))

data <- data.frame("Population"=c(x30,x52,x70,x02),"Density"=c(y30,y52,y70,y02),"Year"=c(l30,l52,l70,l02))

fig4.d <- ggplot(data,aes(x=Population,y=Density,color=Year))+theme_bw()+geom_line(aes(linetype=Year),size=1)+theme(legend.title=element_blank(),legend.position=c(0.9,0.85),axis.text = element_text(size = theme.size, colour="black"),axis.title = element_text(size = 20),legend.text=element_text(size=24),legend.background=element_blank(),legend.key.size = unit(2,"line"))+scale_color_manual(values=c("black","red","blue","Darkgreen"))
fig4.d


#Markov

mkv(fmla_data[,2:9])[[1]] %>% round(4)


#SP Markov
sp.mkv(fmla_data[,2:9],Wd1,fixed=FALSE)

sp.homo.test(fmla_data[,2:9],Wd1,fixed=FALSE) # No se rechaza la nula de homogeneidad espacial

#El comportamiento mejora al reducir el numero de clases a 4

sp.mkv(fmla_data[,2:9],Wd1,fixed=FALSE,classes=4)

sp.homo.test(fmla_data[,2:9],Wd1,fixed=FALSE,classes=4)

# LISA Markov

lisamkv(fmla_data[,2:9],Wd1)[[3]] %>% round(4)

lisamkv(fmla_data[,2:9],Wd1)[[3]] %>% st.st %>% round(4) #algo pasa que la st.st se acumula solo en el 2, puede ser qeu los eigen vector no se encuentren, voy a ver paso por paso.

join.d(fmla_data[,2:9],Wd1) #al parecer la dinamica conjunta es dependiente, pero se contradice con el test anterior

#GIMA

WR2 <- matrix(0,nrow(base1),nrow(base1))
for(i in seq_len(nrow(base1))){
	WR2[,i] <- 1*(Regime1[i]==Regime1)
}
diag(WR2) <- 0

WR2 <- mat2listw(WR2)

spTau <- lapply(2:(ncol(fmla_data)-1),function(i) sp.tau(fmla_data[,i],fmla_data[,(i+1)],WR2,999))

GIMA <- do.call(rbind.data.frame,spTau)

colnames(GIMA)<-names(spTau[[1]])


p <- matrix(1,nrow(b),nrow(b))
diag(p)<-0

phi <- spweights.constants(WR)$S0/sum(p)

nonng.tau <- (GIMA[,5]-(phi*GIMA[,1]))/(1-phi)
ng.tau <- (GIMA[,1])
ng.M <- (ng.tau-1)/-2
nonng.M <- (nonng.tau-1)/-2
pvalue <- GIMA[,2]

tabla.GIMA_2 <- data.frame("Period"=per,"MW"=ng.M,"MnonW"=nonng.M,"TW"=ng.tau,"TnonW"=nonng.tau,"p-value"=pvalue)

# Rank Decomposition

Reg <- as.numeric(as.factor(Regime1))

th <- theta(fmla_data[,2:9],Reg,nsim=999)

theta.table_1 <- data.frame("Period"=per,"Theta"=th[,1],"p-value"=th[,2])

#Directional LISA



full.p8_1 <- d.LISA(D2[,1],D2[,8],Wd1,Regime=Regime1,nsim=999)
full.p4_1 <- d.LISA(D2[,1],D2[,8],Wd1,Regime=Regime1,nsim=999,k=4)


# k= 8


dL <- lapply(seq_len((ncol(D22)-1)),function(x)d.LISA(D2[,x],D2[,x+1],Wd1,Regime=Regime1,nsim=999))

p1 <- list()
for(i in 1:4){p1[[i]]<-dL[[i]]$data} 

p2 <- list()
for(i in 5:7){p2[[i]]<-dL[[i]]$data}

pt <- list()
for(i in 1:7){pt[[i]]<-dL[[i]]$data}


p1 <- do.call(rbind.data.frame,p1)
p2 <- do.call(rbind.data.frame,p2)
pt <- do.call(rbind.data.frame,pt)


# total
quartz()
rose.pt_1<-ggplot(pt,aes(angle,fill=Regime),xlab=" ",ylab=" ") +geom_bar(width=1,colour="black",size=0.1)+scale_x_discrete(name="",limit=lim)+scale_y_continuous(name="",breaks=seq(0,max(table(pt$angle)),5),labels=seq(0,max(table(pt$angle)),5))+coord_polar(start=pi/4)+theme(legend.position="bottom",panel.grid.major = element_line( color="gray50",linetype="solid"),panel.background = element_blank())+scale_fill_brewer(palette="Set1",direction=-1)
rose.pt_1