################################################################################
# Title:     Visual Tools for Identifying Influential Observations in Bivariate
#            Geostatistical Data
# Authors:   Jerfson B. N. Honório, Fernanda De Bastiani, Isabel S. D. de Oliveira,
#            Manuel J. Galea Rojas
# Journal:   (Target) Austrian Journal of Statistics
# Date:      [Insira a data]
# Version:   1.0
# 
# Description:
# This script contains the full implementation of the methods and analysis 
# presented in the manuscript. It includes:
#   - Simulation of bivariate spatial data on a regular grid
#   - Additive perturbation to simulate influential observations
#   - Cross-semivariogram computation using method of moments
#   - Visualization tools including Hair-plots, Boundary Curves, and 
#     Andrews Curves for outlier detection
#   - Application to real dataset on soil carbon and nitrogen
#
# Reproducibility:
# All code runs in R version 4.3.2 with the packages listed below.
# A fixed random seed ensures reproducible results.
#
# Required Packages:
#   - tidyverse       # data manipulation and plotting
#   - gstat           # spatial modeling and semivariogram
#   - sp              # spatial data structures
#   - ggrepel         # non-overlapping labels in ggplot2
#   - numbers         # utility functions (e.g., prime numbers, divisors)
#   - reshape2        # data reshaping
#   - gridExtra       # combining ggplot2 plots
#   - MASS            # multivariate normal generation
#   - viridis         # color palettes for plots
#   - RColorBrewer    # additional color palettes
################################################################################


################################################################################
# 1. LOAD REQUIRED PACKAGES AND SET GLOBAL OPTIONS
################################################################################


library(tidyverse)
library(gstat)
library(sp)
library(ggrepel)
library(numbers)
library(reshape2)
library(gridExtra)
library(MASS)
library(viridis)
library(RColorBrewer)
library(MASS)

theme_set(theme_classic() + theme(text = element_text(size = 16), 
                                  axis.text = element_text(size = 14),  
                                  axis.title = element_text(size = 16) 
))


################################################################################
# 2. SIMULATION OF BIVARIATE SPATIAL DATA ON A REGULAR GRID
################################################################################

set.seed(21) # Semente original [cite: 192, 219]
grid <- expand.grid(x = 1:9, y = 1:9)
coordinates(grid) <- ~x+y

v_model <- vgm(psill = 2.0, model = "Exp", range = 5, nugget = 0.5)

g <- gstat(formula = z ~ 1, locations = grid, dummy = TRUE, beta = 3, model = v_model)
z_sim <- predict(g, newdata = grid, nsim = 2) # nsim=2 para bivariada [cite: 193]


dados_bivariados <- as.data.frame(z_sim)
names(dados_bivariados) <- c("x", "y", "Z1", "Z2")

dados_bivariados[21, "Z1"] <- dados_bivariados[21, "Z1"] + 2 * sd(dados_bivariados$Z1)


################################################################################
# 3. CONTAMINATION SCENARIO 1: PERTURBATION OF A SINGLE VARIABLE
################################################################################

c1.1<-c2.1<-dados_bivariados
c1.2<-c2.2<-dados_bivariados
c1.3<-c2.3<-dados_bivariados

c1.1[21,3]<-dados_bivariados[21, "Z1"] + sd(dados_bivariados$Z1)
c1.2[21,3]<-dados_bivariados[21, "Z1"] + 2 * sd(dados_bivariados$Z1)
c1.3[21,3]<-dados_bivariados[21, "Z1"] + 3 * sd(dados_bivariados$Z1)

c2.1[21,3]<-dados_bivariados[21, "Z1"] +  sd(dados_bivariados$Z1)
c2.2[21,3]<-dados_bivariados[21, "Z1"] + 2 * sd(dados_bivariados$Z1)
c2.3[21,3]<-dados_bivariados[21, "Z1"] + 3 * sd(dados_bivariados$Z1)

c2.1[21,4]<-dados_bivariados[21, "Z2"] + sd(dados_bivariados$Z1)
c2.2[21,4]<-dados_bivariados[21, "Z2"] + 2 * sd(dados_bivariados$Z1)
c2.3[21,4]<-dados_bivariados[21, "Z2"] + 3 * sd(dados_bivariados$Z1)


names(c1.1)<-c("x","y","cadmium", "copper")
names(c1.2)<-c("x","y","cadmium", "copper")
names(c1.3)<-c("x","y","cadmium", "copper")

meuse<-c1.1
coordinates(meuse)<-~x+y

meuse_p<-meuse
zeta1<- -3:3*sd(meuse_p@data[,1])
zeta2<- -3:3*sd(meuse_p@data[,2])  

base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
base<-as.data.frame(base)
names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
valores<-NULL

for (z in 1:length(zeta1)) {
  for (p in 1:length(meuse@data[,1])) {
    meuse_p<-meuse
    meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
    v <- variogram(cadmium~copper, data=meuse_p)
    valores[p]<-v$gamma[1]
  }
  base[,z]<-valores
}

################################################################################
# 4. VISUALIZATION FOR SCENARIO 1: HAIR-PLOTS AND ANDREWS CURVES
################################################################################


obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*1.5
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

out<-obsv %>% filter(Row==21)

as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) + 
  geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(1))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(color="red",label=21,x="3sd",y=1.62,size = 5, family = "sans", check_overlap = TRUE) +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")

ggsave("c1.1.pdf", height = 4.13, width = 5.83, units = "in", dpi = 800)

meuse<-c1.2
coordinates(meuse)<-~x+y

meuse_p<-meuse
zeta1<- -3:3*sd(meuse_p@data[,1])
zeta2<- -3:3*sd(meuse_p@data[,2])  

base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
base<-as.data.frame(base)
names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
valores<-NULL

for (z in 1:length(zeta1)) {
  for (p in 1:length(meuse@data[,1])) {
    meuse_p<-meuse
    meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
    v <- variogram(cadmium~copper, data=meuse_p)
    valores[p]<-v$gamma[1]
  }
  base[,z]<-valores
}


obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*1.5
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")



out<-obsv %>% filter(Row==21)


as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(2))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(color="red",size = 5,label=21,x="3sd",y=1.9, family = "sans", check_overlap = TRUE) +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none") 

ggsave("c1.2.pdf",, height = 4.13, width = 5.83, units = "in", dpi = 800)


meuse<-c1.3
coordinates(meuse)<-~x+y

meuse_p<-meuse
zeta1<- -3:3*sd(meuse_p@data[,1])
zeta2<- -3:3*sd(meuse_p@data[,2])  

base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
base<-as.data.frame(base)
names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
valores<-NULL

for (z in 1:length(zeta1)) {
  for (p in 1:length(meuse@data[,1])) {
    meuse_p<-meuse
    meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
    v <- variogram(cadmium~copper, data=meuse_p)
    valores[p]<-v$gamma[1]
  }
  base[,z]<-valores
}



obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*1.5
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")


out<-obsv %>% filter(Row==21)

as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(3))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(color="red",size = 5,label=21,x="3sd",y=2.3, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")

ggsave("c1.3.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)

################################################################################
# 5. CONTAMINATION SCENARIO 2: PERTURBATION OF BOTH VARIABLES
################################################################################

names(c2.1)<-c("x","y","cadmium", "copper")
names(c2.2)<-c("x","y","cadmium", "copper")
names(c2.3)<-c("x","y","cadmium", "copper")

meuse<-c2.1
coordinates(meuse)<-~x+y

meuse_p<-meuse
zeta1<- -3:3*sd(meuse_p@data[,1])
zeta2<- -3:3*sd(meuse_p@data[,2])  

base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
base<-as.data.frame(base)
names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
valores<-NULL

for (z in 1:length(zeta1)) {
  for (p in 1:length(meuse@data[,1])) {
    meuse_p<-meuse
    meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
    v <- variogram(cadmium~copper, data=meuse_p)
    valores[p]<-v$gamma[1]
  }
  base[,z]<-valores
}

################################################################################
# 6. VISUALIZATION FOR SCENARIO 2: HAIR-PLOTS AND ANDREWS CURVES
################################################################################

obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*1.5
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")


out<-obsv %>% filter(Row==21)


as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(1))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(color="red",size = 5,label=21,x="3sd",y=1.58, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")

ggsave("c2.1.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)

meuse<-c2.2
coordinates(meuse)<-~x+y

meuse_p<-meuse
zeta1<- -3:3*sd(meuse_p@data[,1])
zeta2<- -3:3*sd(meuse_p@data[,2])  

base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
base<-as.data.frame(base)
names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
valores<-NULL

for (z in 1:length(zeta1)) {
  for (p in 1:length(meuse@data[,1])) {
    meuse_p<-meuse
    meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
    v <- variogram(cadmium~copper, data=meuse_p)
    valores[p]<-v$gamma[1]
  }
  base[,z]<-valores
}

obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*1.5
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")


out<-obsv %>% filter(Row==21)


as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(2))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(color="red",size = 5,label=21,x="3sd",y=1.74, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")

ggsave("c2.2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


pertubacao<-function(meuse, lag){
  coordinates(meuse)<-~x+y
  meuse_p<-meuse
  zeta1<- -3:3*sd(meuse_p@data[,1])
  zeta2<- -3:3*sd(meuse_p@data[,2])  
  
  base<-matrix(NA, length(zeta1), nrow = length(meuse@data[,1]))
  base<-as.data.frame(base)
  names(base)<-c("-3sd","-2sd","-1sd", "sp","1sd","2sd","3sd")
  valores<-NULL
  
  for (z in 1:length(zeta1)) {
    for (p in 1:length(meuse@data[,1])) {
      meuse_p<-meuse
      meuse_p@data[p,] <- c(meuse_p@data[p,1] + zeta1[z], meuse_p@data[p,2] + zeta2[z])
      v <- variogram(cadmium~copper, data=meuse_p)
      valores[p]<-v$gamma[lag]
    }
    base[,z]<-valores
  }
  return(base)
}


base1<-pertubacao(c1.1,1)
base2<-pertubacao(c1.1,2)
base3<-pertubacao(c1.1,3)
base4<-pertubacao(c1.1,4)
base5<-pertubacao(c1.1,5)


base1$lag<-"1"
base2$lag<-"2"
base3$lag<-"3"
base4$lag<-"4"
base5$lag<-"5"


bf<-rbind(base1,base2)
bf<-rbind(bf,base3)
bf<-rbind(bf,base4)
bf<-rbind(bf,base5)

t=seq(-pi, pi, length.out=100)
root2 <- sqrt(2)


bf %>%  melt(id.vars = "lag") %>% filter( variable=="-3sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") 


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.5,y=2.8, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t")

ggsave("n3sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="-2sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.6,y=2.2, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t")
ggsave("n2sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="-1sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.5,y=3, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none") + xlab("t")

ggsave("n1sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="1sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.2,y=5.3, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none") + xlab("t")

ggsave("1sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="2sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) 


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.13,y=6.2, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t") 
ggsave("2sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="3sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1,2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.1,y=7.4, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")+ xlab("t")

ggsave("3sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


t=seq(-pi, pi, length.out=50)
root2 <- sqrt(2)


bf %>%  melt(id.vars = "lag") -> a 

a %>% add_column(Row=rep(1:81,nrow(a)/81)) %>% 
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb



obsv <- nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D")+
  geom_text(size = 5,label=21,x=0,y=2, color = "red", family = "sans", check_overlap = TRUE) +
  geom_line(data=out,aes(x=T, y=Andrews), color = "red")+
  theme(legend.position = "none")+ xlab("t")

ggsave("tlag.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


base1<-pertubacao(c2.1,1)
base2<-pertubacao(c2.1,2)
base3<-pertubacao(c2.1,3)
base4<-pertubacao(c2.1,4)
base5<-pertubacao(c2.1,5)


base1$lag<-"1"
base2$lag<-"2"
base3$lag<-"3"
base4$lag<-"4"
base5$lag<-"5"


bf<-rbind(base1,base2)
bf<-rbind(bf,base3)
bf<-rbind(bf,base4)
bf<-rbind(bf,base5)

t=seq(-pi, pi, length.out=100)
root2 <- sqrt(2)
q=3.5


bf %>%  melt(id.vars = "lag") %>% filter( variable=="-3sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.5,y=2.5, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t") 

ggsave("n3sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="-2sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.5,y=2.5, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t")
ggsave("n2sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="-1sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.5,y=3.4, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")+ xlab("t") 

ggsave("n1sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="1sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.1,y=5.1, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")+ xlab("t") 

ggsave("1sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="2sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) ;j=4

obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0.1,y=6.45, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) + theme(legend.position = "none")+ xlab("t") 
ggsave("2sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="3sd") %>% 
  mutate(Row=rep(1:81,5), .before=1) %>%
  filter(lag %in% c(1,2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)) 

out<-obsv %>% filter(Row==21)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  geom_text(color="red",size = 5,label=21,x=0,y=6.4, family = "sans", check_overlap = TRUE)  +
  geom_path(data=out,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8)+ theme(legend.position = "none")+ xlab("t")

ggsave("3sd2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


t=seq(-pi, pi, length.out=50)
root2 <- sqrt(2)


bf %>%  melt(id.vars = "lag") -> a 

a %>% add_column(Row=rep(1:81,nrow(a)/81)) %>% 
  filter(lag %in% c(1, 2,3,4,5)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`) -> nb

obsv<-nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))  %>% filter(Row==21)

nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D")+
  geom_text(size = 5,label=21,x=1.4,y=3, color = "red", family = "sans", check_overlap = TRUE) +
  geom_line(data=obsv,aes(x=T, y=Andrews), color = "red")+
  theme(legend.position = "none")+ xlab("t")

ggsave("tlag2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


################################################################################
# 7. REAL DATA ANALYSIS: SOIL CARBON AND NITROGEN IN SOUTHERN WISCONSIN
################################################################################

meuse<-read.csv(file ="soil_data.csv",sep = ",")
names(meuse)<-c("x","y","solo","depth", "mag")

ggplot(meuse, aes(x=x, y=y, color = depth, size= mag)) +
  geom_point() +geom_text_repel(aes(label = 1:nrow(meuse)),size = 4,vjust = -3) +
  scale_x_continuous(breaks = c(43, 43.4), labels = c("43", "43.4")) +
  facet_wrap(~solo, ncol =3) +
  labs(
    x = "Lat",
    y = "Lon",
    color = "SOC",
    size = "SON"
  ) +scale_color_viridis(discrete=F,  option = "D") + theme_light()+ theme(
    legend.text = element_text(size = 9),     
    legend.title = element_text(size = 9)) 

ggsave("maps.pdf",  height = 5.13, width = 6.83, units = "in", dpi = 800)


latlong<-meuse
latlong<-meuse[,1:2]

coordinates(latlong)<-~x+y

proj4string(latlong) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
utm <- spTransform(latlong, CRS("+proj=utm +zone=16 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

utm_coordinates <- coordinates(utm)
meuse<-cbind(utm_coordinates,meuse)
names(meuse)<-c("x","y","x1","y1","solo","cadmium", "copper")

meuse<-meuse[,-(3:5)]

base1<-pertubacao(meuse,1)
base2<-pertubacao(meuse,2)
base3<-pertubacao(meuse,3)
base4<-pertubacao(meuse,4)
base5<-pertubacao(meuse,5)
base6<-pertubacao(meuse,6)
base7<-pertubacao(meuse,7)

base = base1
obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>%
  melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*q
quantiles[1,] <- quantiles[1,] - RIC
quantiles[2,] <- quantiles[2,] + RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%
  mutate(Row=row_number(), .before=1) %>%
  melt(id.vars = "Row")

out<-obsv %>% filter(Row==10 | Row==85 | Row==124 | Row==88 | Row==6 )


as_tibble(base1) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +  
  geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +  
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(1))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(size = 5,label=10,x="3sd",y=2.65, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=85,x="3sd",y=2.53, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=124,x="3sd",y=2.39, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=88,x="-3sd",y=2.5, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=6,x="-3sd",y=2.40, family = "sans", check_overlap = TRUE,color = "red") +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.5) + theme(legend.position = "none") 


ggsave("r1.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


base = base2
obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*j
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv <- as_tibble(base2) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

out<-obsv %>% filter(Row==10 | Row==85|Row==110)

as_tibble(base2) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +  
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(2))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(size = 5,label=10,x="3sd",y=4.13, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=85,x="3sd",y=4.04, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=110,x="-3sd",y=4.04, family = "sans", check_overlap = TRUE,color = "red") +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.5) +
  theme(legend.position = "none") 

ggsave("r2.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


base = base3
obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*q
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")




obsv <- as_tibble(base3) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

out<-obsv %>% filter(Row==10 | Row==85)

as_tibble(base3) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +  
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(3))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(size = 5,label=85,x="3sd",y=4.61, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=10,x="3sd",y=4.92, family = "sans", check_overlap = TRUE,color = "red") +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.5) + theme(legend.position = "none") 

ggsave("r3.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


base = base4
obsv <- as_tibble(base) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

obsv$x=rep(-3:3, each = length(base[,1]))

quantiles <- apply(base, 2, function(x) quantile(x, c(0.25, 0.75)))
RIC=(quantiles[2,]-quantiles[1,])*q
quantiles[1,] <- quantiles[1,] -RIC
quantiles[2,] <- quantiles[2,] +RIC
quantiles <-as.data.frame(quantiles)[,-median(1:(ncol(base)))] %>%  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")


obsv <- as_tibble(base4) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row")

out<-obsv %>% filter(Row==10 | Row==85)

as_tibble(base4) %>%
  mutate(Row=row_number(), .before=1) %>% melt(id.vars = "Row") %>% 
  ggplot() +
  geom_line(aes(x=variable , y=value, group=Row, color= variable)) +   geom_line(data = quantiles, aes(x = variable, y = value,group = Row),linetype = "dashed",linewidth=0.7, color = "black") +  
  scale_x_discrete(breaks = c("-3sd", "-2sd", "-1sd", "sp", "1sd","2sd","3sd"), labels = c('-3', '-2', '-1', '0', '1',"2","3")) +
  xlab(expression(zeta)) +
  ylab(expression(2 * gamma(4))) +
  scale_color_viridis(discrete=T,  option = "D") +
  geom_text(size = 5,label=10,x="3sd",y=2.94, family = "sans", check_overlap = TRUE,color = "red") +
  geom_text(size = 5,label=85,x="3sd",y=2.69, family = "sans", check_overlap = TRUE,color = "red") +
  geom_path(data=out,aes(x=variable , y=value, group=Row), col='red', lwd = 0.5) + theme(legend.position = "none") 

ggsave("r4.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)




base1$lag<-"1"
base2$lag<-"2"
base3$lag<-"3"
base4$lag<-"4"
base5$lag<-"5"
base6$lag<-"6"
base7$lag<-"7"



bf<-rbind(base1,base2)
bf<-rbind(bf,base3)
bf<-rbind(bf,base4)
bf<-rbind(bf,base5)
bf<-rbind(bf,base6)
bf<-rbind(bf,base7)

t=seq(-pi, pi, length.out=100)
root2 <- sqrt(2)

bf %>%  melt(id.vars = "lag") %>% filter( variable=="-3sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))

out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)
out3<-obsv %>% filter(Row==124)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=2,y=2, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-1.3,y=-1.4, color = "red", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=124,x=1.1,y=9.5, color = "#F314D4", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8)+
  geom_path(data=out3,aes(x=T, y=Andrews, group=Row), col="#F314D4", lwd = 0.8) + xlab("t")

ggsave("rn3sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


bf %>%  melt(id.vars = "lag") %>% filter( variable=="-2sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() + geom_point(aes(x=T, y=Andrews, group=Row, color= factor(Row)))  + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))


out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)
out3<-obsv %>% filter(Row==124)

nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=2,y=2, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-1.2,y=-1.6, color = "red", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=124,x=1.1,y=9.5, color = "#F314D4", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8)+
  geom_path(data=out3,aes(x=T, y=Andrews, group=Row), col="#F314D4", lwd = 0.8)+ xlab("t")

ggsave("rn2sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)



bf %>%  melt(id.vars = "lag") %>% filter( variable=="-1sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))


out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)
nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=2,y=2, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-1.3,y=-1.4, color = "red", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8)+ xlab("t")  

ggsave("rn1sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


bf %>%  melt(id.vars = "lag") %>% filter( variable=="3sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))


out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)
out3<-obsv %>% filter(Row==124)

nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=1.6,y=3.32, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-0.7,y=-5, color = "red", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=124,x=0.4,y=9, color = "#F314D4", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8)+
  geom_path(data=out3,aes(x=T, y=Andrews, group=Row), col="#F314D4", lwd = 0.8)+ xlab("t")

ggsave("r3sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


bf %>%  melt(id.vars = "lag") %>% filter( variable=="2sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))


out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)
out3<-obsv %>% filter(Row==124)

nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=1.6,y=3.32, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-0.7,y=-5, color = "red", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=124,x=0.4,y=9, color = "#F314D4", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8) +
  geom_path(data=out3,aes(x=T, y=Andrews, group=Row), col="#F314D4", lwd = 0.8)+ xlab("t")
ggsave("r2sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


bf %>%  melt(id.vars = "lag") %>% filter( variable=="1sd") %>% 
  mutate(Row=rep(1:124,7), .before=1) %>%
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t) + lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= Row)) + theme(legend.position = "none")


obsv <- nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))


out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)


nb %>%  
  uncount(100) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t)+ lag6*sin(3*t) + lag7*cos(3*t))    %>% 
  ggplot() +
  geom_text(size = 5,label=85,x=1.5,y=3.32, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=10,x=-0.7,y=-4, color = "red", family = "sans", check_overlap = TRUE) +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D") +
  theme(legend.position = "none") +
  geom_path(data=out1,aes(x=T, y=Andrews, group=Row), col='red', lwd = 0.8) +
  geom_path(data=out2,aes(x=T, y=Andrews, group=Row), col="#7e0f7e", lwd = 0.8)+ xlab("t") 
ggsave("r1sd.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)


t=seq(-pi, pi, length.out=50)
root2 <- sqrt(2)


bf %>%  melt(id.vars = "lag") -> a 

a %>% add_column(Row=rep(1:124,nrow(a)/124)) %>% 
  filter(lag %in% c(1, 2,3,4,5,6,7)) %>%  # Filtra as observações com lag 1 e 2
  pivot_wider(names_from = lag, values_from = value) %>%
  rename(lag1 = `1`, lag2 = `2`, lag3 = `3`, lag4 = `4`, lag5 = `5`, lag6 = `6`, lag7 = `7`) -> nb


obsv<-nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))

out1<-obsv %>% filter(Row==10)
out2<-obsv %>% filter(Row==85)


nb %>%  
  uncount(50) %>% 
  add_column(T=rep(t, times=nrow(nb))) %>% 
  mutate(Andrews=lag1/root2 + lag2*sin(t) + lag3*cos(t) + lag4*sin(2*t) + lag5*cos(2*t))   %>% 
  ggplot() +
  geom_line(aes(x=T, y=Andrews, group=Row, color= T)) +
  scale_color_viridis(discrete=F,  option = "D")+
  geom_text(size = 5,label=10,x=1.5,y=7, color = "red", family = "sans", check_overlap = TRUE) +
  geom_text(size = 5,label=85,x=-0.5,y=7, color = "#7e0f7e", family = "sans", check_overlap = TRUE) +
  geom_line(data=out1,aes(x=T, y=Andrews), color = "red")+
  geom_line(data=out2,aes(x=T, y=Andrews), color = "#7e0f7e", alpha = 0.8)+
  theme(legend.position = "none")+ xlab("t")

ggsave("trl.pdf",  height = 4.13, width = 5.83, units = "in", dpi = 800)