library(GeoModels)
library(geoR)
library(fields)
library(RandomFields)
#### SIMULATIONS CASE 1--4
#### DEFINE GRID AND REAL PARAMETERS
x <- seq(-3/2, 3/2, length.out = 20)
y <- seq(-3/2, 3/2, length.out = 20)
grid=expand.grid(x,y)
param=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,sill_1=1,sill_2=1,pcol=0.3,power2_1=4, power2_12=4, power2_2=4,scale_1=1.8 ,scale_12=1.4 ,scale_2=1.5 )

# SIMULATIONS  OF WENDLAND GAUSSIAN FIELD:
data=list()
for (i in 1:500){
  data[[i]] <- GeoSim(coordx=grid, corrmodel="Bi_Wend1", param=param,n=1)$data
}

# FIXED AND START PARAMETERS
fixed=list(mean_1=0,mean_2=0,power2_1=4, power2_12=4, power2_2=4,nugget_1=0,nugget_2=0)
start=list(sill_1=1,sill_2=1,pcol=0.3,scale_1=1.8 ,scale_12=1.4 ,scale_2=1.5  )
# FITING THE MODEL
fitwend=list()
for (i in 1:500){
  fitwend[[i]]<- GeoFit(data=data[[i]], coordx=grid, corrmodel="Bi_Wend1",maxdist = max(dist(grid)),n=1, start=start,fixed=fixed)
}


suma=matrix(rep(0,500*6),ncol=6)
for (i in 1:500){
  suma[i,]=fitwend[[i]]$param
}
mean_caso3=c(mean(suma[,1]),mean(suma[,2]),mean(suma[,3]),mean(suma[,4]),mean(suma[,5]),mean(suma[,6]))


############# CASE 1
data_caso1=list()
param_caso1=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,sill_1=1,sill_2=1,pcol=-0.15,power2_1=4, power2_12=4, power2_2=4,scale_1=0.5 ,scale_12=0.35 ,scale_2=0.4 )
start_caso1=list(sill_1=1,sill_2=1,pcol=-0.15,scale_1=0.5 ,scale_12=0.35 ,scale_2=0.4  )

for (i in 1:500){
  data_caso1[[i]] <- GeoSim(coordx=grid, corrmodel="Bi_Wend1", param=param_caso1,n=1)$data
}

############## CASE 2
data_caso2=list()
param_caso2=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,sill_1=1,sill_2=1,pcol=0.25,power2_1=4, power2_12=4, power2_2=4,scale_1=1.2 ,scale_12=1 ,scale_2=0.9 )
start_caso2=list(sill_1=1,sill_2=1,pcol=0.25,scale_1=1.2 ,scale_12=1 ,scale_2=0.9  )

for (i in 1:500){
  data_caso2[[i]] <- GeoSim(coordx=grid, corrmodel="Bi_Wend1", param=param_caso2,n=1)$data
}

fitwend_caso1=list()
fitwend_caso2=list()
for (i in 1:500){
  fitwend_caso1[[i]]<- GeoFit(data=data_caso1[[i]], coordx=grid, corrmodel="Bi_Wend1",maxdist = max(dist(grid)),n=1, start=start_caso1,fixed=fixed)
  fitwend_caso2[[i]]<- GeoFit(data=data_caso2[[i]], coordx=grid, corrmodel="Bi_Wend1",maxdist = max(dist(grid)),n=1, start=start_caso2,fixed=fixed)
}


suma_caso1=matrix(rep(0,500*6),ncol=6)
suma_caso2=matrix(rep(0,500*6),ncol=6)
for (i in 1:500){
  suma_caso1[i,]=fitwend_caso1[[i]]$param
  suma_caso2[i,]=fitwend_caso2[[i]]$param
}
mean_caso1=c(mean(suma_caso1[,1]),mean(suma_caso1[,2]),mean(suma_caso1[,3]),mean(suma_caso1[,4]),mean(suma_caso1[,5]),mean(suma_caso1[,6]))
mean_caso2=c(mean(suma_caso2[,1]),mean(suma_caso2[,2]),mean(suma_caso2[,3]),mean(suma_caso2[,4]),mean(suma_caso2[,5]),mean(suma_caso2[,6]))
image(fitwend[[1]]$data)

##### CASE 4
data_caso4=list()
param_caso4=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,sill_1=2,sill_2=0.5,pcol=0.45,power2_1=4, power2_12=4, power2_2=4,scale_1=1.1 ,scale_12=0.9 ,scale_2=1.1 )
start_caso4=list(sill_1=2,sill_2=0.5,pcol=0.45,scale_1=1.1 ,scale_12=0.9 ,scale_2=1  )
for (i in 1:500){
  data_caso4[[i]] <- GeoSim(coordx=grid, corrmodel="Bi_Wend1", param=param_caso4,n=1)$data
}


fitwend_caso4=list()
for (i in 1:500){
  fitwend_caso4[[i]]<- GeoFit(data=data_caso4[[i]], coordx=grid, corrmodel="Bi_Wend1",maxdist = max(dist(grid)),n=1, start=start_caso4,fixed=fixed)
}

suma_caso4=matrix(rep(0,500*6),ncol=6)
for (i in 1:500){
  suma_caso4[i,]=fitwend_caso4[[i]]$param
}
mean_caso4=c(mean(suma_caso4[,1]),mean(suma_caso4[,2]),mean(suma_caso4[,3]),mean(suma_caso4[,4]),mean(suma_caso4[,5]),mean(suma_caso4[,6]))
mean_caso4



########################################################################
######### SIMULATIONS ON 6 DIFFERENT GRIDS (4x4,8x8,12x12,16x16,20x20,24x24)
simu4=list()
simu8=list()
simu12=list()
simu16=list()
simu20=list()
simu24=list()
######## SIMULATING THE DATA
for (i in 1:500){
  simu4[[i]]=GeoSim(coordx=expand.grid(1:4,1:4), corrmodel="Bi_Wend1", param=param)$data
  simu8[[i]]=GeoSim(coordx=expand.grid(1:8,1:8), corrmodel="Bi_Wend1", param=param)$data
  simu12[[i]]=GeoSim(coordx=expand.grid(1:12,1:12), corrmodel="Bi_Wend1", param=param)$data
  simu16[[i]]=GeoSim(coordx=expand.grid(1:16,1:16), corrmodel="Bi_Wend1", param=param)$data
  simu20[[i]]=GeoSim(coordx=expand.grid(1:20,1:20), corrmodel="Bi_Wend1", param=param)$data
  simu24[[i]]=GeoSim(coordx=expand.grid(1:24,1:24), corrmodel="Bi_Wend1", param=param)$data
}
#### SETTING PARAMETERS
fixed=list(mean_1=0,mean_2=0,power2_1=4, power2_12=4, power2_2=4,nugget_1=0,nugget_2=0)
start=list(sill_1=1,sill_2=1,pcol=0.3,scale_1=1.8 ,scale_12=1.4 ,scale_2=1.5  )
# FITTING MODELS
fitsimu4=list()
fitsimu8=list()
fitsimu12=list()
fitsimu16=list()
fitsimu20=list()
fitsimu24=list()
for (i in 1:500){
  fitsimu4[[i]]<- GeoFit(data=simu4[[i]], coordx=expand.grid(1:4,1:4), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:4,1:4))), start=start,fixed=fixed)
  fitsimu8[[i]]<- GeoFit(data=simu8[[i]], coordx=expand.grid(1:8,1:8), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:8,1:8))), start=start,fixed=fixed)
  fitsimu12[[i]]<- GeoFit(data=simu12[[i]], coordx=expand.grid(1:12,1:12), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:12,1:12))), start=start,fixed=fixed)
  fitsimu16[[i]]<- GeoFit(data=simu16[[i]], coordx=expand.grid(1:16,1:16), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:16,1:16))), start=start,fixed=fixed)
  fitsimu20[[i]]<- GeoFit(data=simu20[[i]], coordx=expand.grid(1:20,1:20), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:20,1:20))), start=start,fixed=fixed)
  fitsimu24[[i]]<- GeoFit(data=simu24[[i]], coordx=expand.grid(1:24,1:24), corrmodel="Bi_Wend1",maxdist = max(dist(expand.grid(1:24,1:24))), start=start,fixed=fixed)
}

#EXTRACTING PARAMETERS OF THE SIMULATIONS
suma_4=matrix(rep(0,500*6),ncol=6)
suma_8=matrix(rep(0,500*6),ncol=6)
suma_12=matrix(rep(0,500*6),ncol=6)
suma_16=matrix(rep(0,500*6),ncol=6)
suma_20=matrix(rep(0,500*6),ncol=6)
suma_24=matrix(rep(0,500*6),ncol=6)
for (i in 1:500){
  suma_4[i,]=fitsimu4[[i]]$param
  suma_8[i,]=fitsimu8[[i]]$param
  suma_12[i,]=fitsimu12[[i]]$param
  suma_16[i,]=fitsimu16[[i]]$param
  suma_20[i,]=fitsimu20[[i]]$param
  suma_24[i,]=fitsimu24[[i]]$param
}

pars4=c(mean(suma_4[,1]),mean(suma_4[,3]),mean(suma_4[,5]),mean(suma_4[,6]),as.numeric(quantile(suma_4[,1],c(0.05,0.95))),as.numeric(quantile(suma_4[,3],c(0.05,0.95))),as.numeric(quantile(suma_4[,5],c(0.05,0.95))),as.numeric(quantile(suma_4[,6],c(0.05,0.95))))
pars8=c(mean(suma_8[,1]),mean(suma_8[,3]),mean(suma_8[,5]),mean(suma_8[,6]),as.numeric(quantile(suma_8[,1],c(0.05,0.95))),as.numeric(quantile(suma_8[,3],c(0.05,0.95))),as.numeric(quantile(suma_8[,5],c(0.05,0.95))),as.numeric(quantile(suma_8[,6],c(0.05,0.95))))
pars12=c(mean(suma_12[,1]),mean(suma_12[,3]),mean(suma_12[,5]),mean(suma_12[,6]),as.numeric(quantile(suma_12[,1],c(0.05,0.95))),as.numeric(quantile(suma_12[,3],c(0.05,0.95))),as.numeric(quantile(suma_12[,5],c(0.05,0.95))),as.numeric(quantile(suma_12[,6],c(0.05,0.95))))
pars16=c(mean(suma_16[,1]),mean(suma_16[,3]),mean(suma_16[,5]),mean(suma_16[,6]),as.numeric(quantile(suma_16[,1],c(0.05,0.95))),as.numeric(quantile(suma_16[,3],c(0.05,0.95))),as.numeric(quantile(suma_16[,5],c(0.05,0.95))),as.numeric(quantile(suma_16[,6],c(0.05,0.95))))
pars20=c(mean(suma_20[,1]),mean(suma_20[,3]),mean(suma_20[,5]),mean(suma_20[,6]),as.numeric(quantile(suma_20[,1],c(0.05,0.95))),as.numeric(quantile(suma_20[,3],c(0.05,0.95))),as.numeric(quantile(suma_20[,5],c(0.05,0.95))),as.numeric(quantile(suma_20[,6],c(0.05,0.95))))
pars24=c(mean(suma_24[,1]),mean(suma_24[,3]),mean(suma_24[,5]),mean(suma_24[,6]),as.numeric(quantile(suma_24[,1],c(0.05,0.95))),as.numeric(quantile(suma_24[,3],c(0.05,0.95))),as.numeric(quantile(suma_24[,5],c(0.05,0.95))),as.numeric(quantile(suma_24[,6],c(0.05,0.95))))

#### CREATING A MATRIX WITH ALL THE DATA
parssimus=matrix(rep(0,6*12),ncol=12)
parssimus[1,]=pars4
parssimus[2,]=pars8
parssimus[3,]=pars12
parssimus[4,]=pars16
parssimus[5,]=pars20
parssimus[6,]=pars24


#save(parssimus,file="Pars_simus_6casos.RData")  

###### SIMULATIONS LOCAL APPROACH

x1=y1=1:50
grid1=expand.grid(x1,y1)
param=list(mean_1=0,mean_2=0,nugget_1=0,nugget_2=0,sill_1=1,sill_2=1,pcol=0.3,power2_1=4, power2_12=4, power2_2=4,scale_1=1.8 ,scale_12=1.5 ,scale_2=1.4 )
data_loc=list()
for (i in 1:100){
  data_loc[[i]] <- GeoSim(coordx=grid1, corrmodel="Bi_Wend1", param=param,n=1)$data
}
grid_local=list()
grid=expand.grid(1:50,1:50)
c1=c(1:10,51:60,101:110,151:160,201:210,251:260,301:310,351:360,401:410,451:460)
c2=c(11:20,61:70,111:120,161:170,211:220,261:270,311:320,361:370,411:420,461:470)
c3=c(21:30,71:80,121:130,171:180,221:230,271:280,321:330,371:380,421:430,471:480)
c4=c(31:40,81:90,131:140,181:190,231:240,281:290,331:340,381:390,431:440,481:490)
c5=c(41:50,91:100,141:150,191:200,241:250,291:300,341:350,391:400,441:450,491:500)

c6=c(501:510,551:560,601:610,651:660,701:710,751:760,801:810,851:860,901:910,951:960)
c7=c(511:520,561:570,611:620,661:670,711:720,761:770,811:820,861:870,911:920,961:970)
c8=c(521:530,571:580,621:630,671:680,721:730,771:780,821:830,871:880,921:930,971:980)
c9=c(531:540,581:590,631:640,681:690,731:740,781:790,831:840,881:890,931:940,981:990)
c10=c(541:550,591:600,641:650,691:700,741:750,791:800,841:850,891:900,941:950,991:1000)

c11=c(1001:1010,1051:1060,1101:1110,1151:1160,1201:1210,1251:1260,1301:1310,1351:1360,1401:1410,1451:1460)
c12=c(1011:1020,1061:1070,1111:1120,1161:1170,1211:1220,1261:1270,1311:1320,1361:1370,1411:1420,1461:1470)
c13=c(1021:1030,1071:1080,1121:1130,1171:1180,1221:1230,1271:1280,1321:1330,1371:1380,1421:1430,1471:1480)
c14=c(1031:1040,1081:1090,1131:1140,1181:1190,1231:1240,1281:1290,1331:1340,1381:1390,1431:1440,1481:1490)
c15=c(1041:1050,1091:1100,1141:1150,1191:1200,1241:1250,1291:1300,1341:1350,1391:1400,1441:1450,1491:1500)

c16=c(1501:1510,1551:1560,1601:1610,1651:1660,1701:1710,1751:1760,1801:1810,1851:1860,1901:1910,1951:1960)
c17=c(1511:1520,1561:1570,1611:1620,1661:1670,1711:1720,1761:1770,1811:1820,1861:1870,1911:1920,1961:1970)
c18=c(1521:1530,1571:1580,1621:1630,1671:1680,1721:1730,1771:1780,1821:1830,1871:1880,1921:1930,1971:1980)
c19=c(1531:1540,1581:1590,1631:1640,1681:1690,1731:1740,1781:1790,1831:1840,1881:1890,1931:1940,1981:1990)
c20=c(1541:1550,1591:1600,1641:1650,1691:1700,1741:1750,1791:1800,1841:1850,1891:1900,1941:1950,1991:2000)

c21=c(2001:2010,2051:2060,2101:2110,2151:2160,2201:2210,2251:2260,2301:2310,2351:2360,2401:2410,2451:2460)
c22=c(2011:2020,2061:2070,2111:2120,2161:2170,2211:2220,2261:2270,2311:2320,2361:2370,2411:2420,2461:2470)
c23=c(2021:2030,2071:2080,2121:2130,2171:2180,2221:2230,2271:2280,2321:2330,2371:2380,2421:2430,2471:2480)
c24=c(2031:2040,2081:2090,2131:2140,2181:2190,2231:2240,2281:2290,2331:2340,2381:2390,2431:2440,2481:2490)
c25=c(2041:2050,2091:2100,2141:2150,2191:2200,2241:2250,2291:2300,2341:2350,2391:2400,2441:2450,2491:2500)

grid_local[[1]]=grid[c1,]
grid_local[[2]]=grid[c2,]
grid_local[[3]]=grid[c3,]
grid_local[[4]]=grid[c4,]
grid_local[[5]]=grid[c5,]
grid_local[[6]]=grid[c6,]
grid_local[[7]]=grid[c7,]
grid_local[[8]]=grid[c8,]
grid_local[[9]]=grid[c9,]
grid_local[[10]]=grid[c10,]
grid_local[[11]]=grid[c11,]
grid_local[[12]]=grid[c12,]
grid_local[[13]]=grid[c13,]
grid_local[[14]]=grid[c14,]
grid_local[[15]]=grid[c15,]
grid_local[[16]]=grid[c16,]
grid_local[[17]]=grid[c17,]
grid_local[[18]]=grid[c18,]
grid_local[[19]]=grid[c19,]
grid_local[[20]]=grid[c20,]
grid_local[[21]]=grid[c21,]
grid_local[[22]]=grid[c22,]
grid_local[[23]]=grid[c23,]
grid_local[[24]]=grid[c24,]
grid_local[[25]]=grid[c25,]


fixed=list(mean_1=0,mean_2=0,power2_1=4, power2_12=4, power2_2=4,nugget_1=0,nugget_2=0)
start=list(sill_1=1,sill_2=1,pcol=0.3,scale_1=1.8 ,scale_12=1.5 ,scale_2=1.4  )

fitwend_local1=list()
fitwend_local2=list()
fitwend_local3=list()
fitwend_local4=list()
fitwend_local5=list()
fitwend_local6=list()
fitwend_local7=list()
fitwend_local8=list()
fitwend_local9=list()
fitwend_local10=list()
fitwend_local11=list()
fitwend_local12=list()
fitwend_local13=list()
fitwend_local14=list()
fitwend_local15=list()
fitwend_local16=list()
fitwend_local17=list()
fitwend_local18=list()
fitwend_local19=list()
fitwend_local20=list()
fitwend_local21=list()
fitwend_local22=list()
fitwend_local23=list()
fitwend_local24=list()
fitwend_local25=list()
for (i in 1:100){
  fitwend_local1[[i]]<- GeoFit(data=data_loc[[i]][,c1], coordx=grid_local[[1]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[1]])),n=1, start=start,fixed=fixed)
  fitwend_local2[[i]]<- GeoFit(data=data_loc[[i]][,c2], coordx=grid_local[[2]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[2]])),n=1, start=start,fixed=fixed)
  fitwend_local3[[i]]<- GeoFit(data=data_loc[[i]][,c3], coordx=grid_local[[3]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[3]])),n=1, start=start,fixed=fixed)
  fitwend_local4[[i]]<- GeoFit(data=data_loc[[i]][,c4], coordx=grid_local[[4]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[4]])),n=1, start=start,fixed=fixed)
  fitwend_local5[[i]]<- GeoFit(data=data_loc[[i]][,c5], coordx=grid_local[[5]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[5]])),n=1, start=start,fixed=fixed)
  fitwend_local6[[i]]<- GeoFit(data=data_loc[[i]][,c6], coordx=grid_local[[6]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[6]])),n=1, start=start,fixed=fixed)
  fitwend_local7[[i]]<- GeoFit(data=data_loc[[i]][,c7], coordx=grid_local[[7]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[7]])),n=1, start=start,fixed=fixed)
  fitwend_local8[[i]]<- GeoFit(data=data_loc[[i]][,c8], coordx=grid_local[[8]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[8]])),n=1, start=start,fixed=fixed)
  fitwend_local9[[i]]<- GeoFit(data=data_loc[[i]][,c9], coordx=grid_local[[9]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[9]])),n=1, start=start,fixed=fixed)
  fitwend_local10[[i]]<- GeoFit(data=data_loc[[i]][,c10], coordx=grid_local[[10]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[10]])),n=1, start=start,fixed=fixed)
  fitwend_local11[[i]]<- GeoFit(data=data_loc[[i]][,c11], coordx=grid_local[[11]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[11]])),n=1, start=start,fixed=fixed)
  fitwend_local12[[i]]<- GeoFit(data=data_loc[[i]][,c12], coordx=grid_local[[12]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[12]])),n=1, start=start,fixed=fixed)
  fitwend_local13[[i]]<- GeoFit(data=data_loc[[i]][,c13], coordx=grid_local[[13]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[13]])),n=1, start=start,fixed=fixed)
  fitwend_local14[[i]]<- GeoFit(data=data_loc[[i]][,c14], coordx=grid_local[[14]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[14]])),n=1, start=start,fixed=fixed)
  fitwend_local15[[i]]<- GeoFit(data=data_loc[[i]][,c15], coordx=grid_local[[15]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[15]])),n=1, start=start,fixed=fixed)
  fitwend_local16[[i]]<- GeoFit(data=data_loc[[i]][,c16], coordx=grid_local[[16]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[16]])),n=1, start=start,fixed=fixed)
  fitwend_local17[[i]]<- GeoFit(data=data_loc[[i]][,c17], coordx=grid_local[[17]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[17]])),n=1, start=start,fixed=fixed)
  fitwend_local18[[i]]<- GeoFit(data=data_loc[[i]][,c18], coordx=grid_local[[18]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[18]])),n=1, start=start,fixed=fixed)
  fitwend_local19[[i]]<- GeoFit(data=data_loc[[i]][,c19], coordx=grid_local[[19]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[19]])),n=1, start=start,fixed=fixed)
  fitwend_local20[[i]]<- GeoFit(data=data_loc[[i]][,c20], coordx=grid_local[[20]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[20]])),n=1, start=start,fixed=fixed)
  fitwend_local21[[i]]<- GeoFit(data=data_loc[[i]][,c21], coordx=grid_local[[21]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[21]])),n=1, start=start,fixed=fixed)
  fitwend_local22[[i]]<- GeoFit(data=data_loc[[i]][,c22], coordx=grid_local[[22]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[22]])),n=1, start=start,fixed=fixed)
  fitwend_local23[[i]]<- GeoFit(data=data_loc[[i]][,c23], coordx=grid_local[[23]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[23]])),n=1, start=start,fixed=fixed)
  fitwend_local24[[i]]<- GeoFit(data=data_loc[[i]][,c24], coordx=grid_local[[24]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[24]])),n=1, start=start,fixed=fixed)
  fitwend_local25[[i]]<- GeoFit(data=data_loc[[i]][,c25], coordx=grid_local[[25]], corrmodel="Bi_Wend1",maxdist = max(dist(grid_local[[25]])),n=1, start=start,fixed=fixed)
}



grid[cbind(1:10, 1:10)]


suma_local1=matrix(rep(0,100*6),ncol=6)
suma_local2=matrix(rep(0,100*6),ncol=6)
suma_local3=matrix(rep(0,100*6),ncol=6)
suma_local4=matrix(rep(0,100*6),ncol=6)
suma_local5=matrix(rep(0,100*6),ncol=6)
suma_local6=matrix(rep(0,100*6),ncol=6)
suma_local7=matrix(rep(0,100*6),ncol=6)
suma_local8=matrix(rep(0,100*6),ncol=6)
suma_local9=matrix(rep(0,100*6),ncol=6)
suma_local10=matrix(rep(0,100*6),ncol=6)
suma_local11=matrix(rep(0,100*6),ncol=6)
suma_local12=matrix(rep(0,100*6),ncol=6)
suma_local13=matrix(rep(0,100*6),ncol=6)
suma_local14=matrix(rep(0,100*6),ncol=6)
suma_local15=matrix(rep(0,100*6),ncol=6)
suma_local16=matrix(rep(0,100*6),ncol=6)
suma_local17=matrix(rep(0,100*6),ncol=6)
suma_local18=matrix(rep(0,100*6),ncol=6)
suma_local19=matrix(rep(0,100*6),ncol=6)
suma_local20=matrix(rep(0,100*6),ncol=6)
suma_local21=matrix(rep(0,100*6),ncol=6)
suma_local22=matrix(rep(0,100*6),ncol=6)
suma_local23=matrix(rep(0,100*6),ncol=6)
suma_local24=matrix(rep(0,100*6),ncol=6)
suma_local25=matrix(rep(0,100*6),ncol=6)

for (i in 1:100){
  suma_local1[i,]=fitwend_local1[[i]]$param
  suma_local2[i,]=fitwend_local2[[i]]$param
  suma_local3[i,]=fitwend_local3[[i]]$param
  suma_local4[i,]=fitwend_local4[[i]]$param
  suma_local5[i,]=fitwend_local5[[i]]$param
  suma_local6[i,]=fitwend_local6[[i]]$param
  suma_local7[i,]=fitwend_local7[[i]]$param
  suma_local8[i,]=fitwend_local8[[i]]$param
  suma_local9[i,]=fitwend_local9[[i]]$param
  suma_local10[i,]=fitwend_local10[[i]]$param
  suma_local11[i,]=fitwend_local11[[i]]$param
  suma_local12[i,]=fitwend_local12[[i]]$param
  suma_local13[i,]=fitwend_local13[[i]]$param
  suma_local14[i,]=fitwend_local14[[i]]$param
  suma_local15[i,]=fitwend_local15[[i]]$param
  suma_local16[i,]=fitwend_local16[[i]]$param
  suma_local17[i,]=fitwend_local17[[i]]$param
  suma_local18[i,]=fitwend_local18[[i]]$param
  suma_local19[i,]=fitwend_local19[[i]]$param
  suma_local20[i,]=fitwend_local20[[i]]$param
  suma_local21[i,]=fitwend_local21[[i]]$param
  suma_local22[i,]=fitwend_local22[[i]]$param
  suma_local23[i,]=fitwend_local23[[i]]$param
  suma_local24[i,]=fitwend_local24[[i]]$param
  suma_local25[i,]=fitwend_local25[[i]]$param
  
}



### BOX PLOTS OF ESTIMATIONS
box_rho12=boxplot(suma_local1[,1],suma_local2[,1],suma_local3[,1],suma_local4[,1],suma_local5[,1],suma_local6[,1],suma_local7[,1],suma_local8[,1],suma_local9[,1],suma_local10[,1],suma_local11[,1],suma_local12[,1],suma_local13[,1],suma_local14[,1],suma_local15[,1],suma_local16[,1],suma_local17[,1],suma_local18[,1],suma_local19[,1],suma_local20[,1],suma_local21[,1],suma_local22[,1],suma_local23[,1],suma_local24[,1],suma_local25[,1])
box_b1=boxplot(suma_local1[,2],suma_local2[,2],suma_local3[,2],suma_local4[,2],suma_local5[,2],suma_local6[,2],suma_local7[,2],suma_local8[,2],suma_local9[,2],suma_local10[,2],suma_local11[,2],suma_local12[,2],suma_local13[,2],suma_local14[,2],suma_local15[,2],suma_local16[,2],suma_local17[,2],suma_local18[,2],suma_local19[,2],suma_local20[,2],suma_local21[,2],suma_local22[,2],suma_local23[,2],suma_local24[,2],suma_local25[,2])
box_b12=boxplot(suma_local1[,3],suma_local2[,3],suma_local3[,3],suma_local4[,3],suma_local5[,3],suma_local6[,3],suma_local7[,3],suma_local8[,3],suma_local9[,3],suma_local10[,3],suma_local11[,3],suma_local12[,3],suma_local13[,3],suma_local14[,3],suma_local15[,3],suma_local16[,3],suma_local17[,3],suma_local18[,3],suma_local19[,3],suma_local20[,3],suma_local21[,3],suma_local22[,3],suma_local23[,3],suma_local24[,3],suma_local25[,3])
box_b2=boxplot(suma_local1[,4],suma_local2[,4],suma_local3[,4],suma_local4[,4],suma_local5[,4],suma_local6[,4],suma_local7[,4],suma_local8[,4],suma_local9[,4],suma_local10[,4],suma_local11[,4],suma_local12[,4],suma_local13[,4],suma_local14[,4],suma_local15[,4],suma_local16[,4],suma_local17[,4],suma_local18[,4],suma_local19[,4],suma_local20[,4],suma_local21[,4],suma_local22[,4],suma_local23[,4],suma_local24[,4],suma_local25[,4])
box_sigma1=boxplot(suma_local1[,5],suma_local2[,5],suma_local3[,5],suma_local4[,5],suma_local5[,5],suma_local6[,5],suma_local7[,5],suma_local8[,5],suma_local9[,5],suma_local10[,5],suma_local11[,5],suma_local12[,5],suma_local13[,5],suma_local14[,5],suma_local15[,5],suma_local16[,5],suma_local17[,5],suma_local18[,5],suma_local19[,5],suma_local20[,5],suma_local21[,5],suma_local22[,5],suma_local23[,5],suma_local24[,5],suma_local25[,5])
box_sigma2=boxplot(suma_local1[,6],suma_local2[,6],suma_local3[,6],suma_local4[,6],suma_local5[,6],suma_local6[,6],suma_local7[,6],suma_local8[,6],suma_local9[,6],suma_local10[,6],suma_local11[,6],suma_local12[,6],suma_local13[,6],suma_local14[,6],suma_local15[,6],suma_local16[,6],suma_local17[,6],suma_local18[,6],suma_local19[,6],suma_local20[,6],suma_local21[,6],suma_local22[,6],suma_local23[,6],suma_local24[,6],suma_local25[,6])

sim_1=im(data_caso1[[4]][1,],1:20,1:20)
sim_1_2=im(data_caso1[[4]][2,],1:20,1:20)
sim1=list(sim_1,sim_1_2)
plot(as.listof(sim1),main="",ribside="bottom",ribsep=0.05,main.panel = c("X(s)","Y(s)"),mar.panel = 0.01,hsep=0,vsep=0)

sim_2=im(data_caso2[[4]][1,],1:20,1:20)
sim_2_2=im(data_caso2[[4]][2,],1:20,1:20)
sim2=list(sim_2,sim_2_2)
plot(as.listof(sim2),main="",ribside="bottom",ribsep=0.05,main.panel = c("X(s)","Y(s)"),mar.panel = 0.01,hsep=0,vsep=0)

sim_3=im(data[[101]][1,],1:20,1:20)
sim_3_2=im(data[[101]][2,],1:20,1:20)
sim3=list(sim_3,sim_3_2)
plot(as.listof(sim3),main="",ribside="bottom",ribsep=0.05,main.panel = c("X(s)","Y(s)"),mar.panel = 0.01,hsep=0,vsep=0)

sim_4=im(data_caso4[[4]][1,],1:20,1:20)
sim_4_2=im(data_caso4[[4]][2,],1:20,1:20)
sim4=list(sim_4,sim_4_2)
plot(as.listof(sim4),main="",ribside="bottom",ribsep=0.05,main.panel = c("X(s)","Y(s)"),mar.panel = 0.01,hsep=0,vsep=0)

par(pty="s",cex.axis=2,cex=1)
b1=boxplot(suma_local1[,1],suma_local2[,1],suma_local3[,1],suma_local4[,1],suma_local5[,1],suma_local6[,1],suma_local7[,1],suma_local8[,1],suma_local9[,1],suma_local10[,1],suma_local11[,1],suma_local12[,1],suma_local13[,1],suma_local14[,1],suma_local15[,1],suma_local16[,1],suma_local17[,1],suma_local18[,1],suma_local19[,1],suma_local20[,1],suma_local21[,1],suma_local22[,1],suma_local23[,1],suma_local24[,1],suma_local25[,1],col="gray",border = "darkblue")
b2=boxplot(suma_local1[,2],suma_local2[,2],suma_local3[,2],suma_local4[,2],suma_local5[,2],suma_local6[,2],suma_local7[,2],suma_local8[,2],suma_local9[,2],suma_local10[,2],suma_local11[,2],suma_local12[,2],suma_local13[,2],suma_local14[,2],suma_local15[,2],suma_local16[,2],suma_local17[,2],suma_local18[,2],suma_local19[,2],suma_local20[,2],suma_local21[,2],suma_local22[,2],suma_local23[,2],suma_local24[,2],suma_local25[,2])
b3=boxplot(suma_local1[,3],suma_local2[,3],suma_local3[,3],suma_local4[,3],suma_local5[,3],suma_local6[,3],suma_local7[,3],suma_local8[,3],suma_local9[,3],suma_local10[,3],suma_local11[,3],suma_local12[,3],suma_local13[,3],suma_local14[,3],suma_local15[,3],suma_local16[,3],suma_local17[,3],suma_local18[,3],suma_local19[,3],suma_local20[,3],suma_local21[,3],suma_local22[,3],suma_local23[,3],suma_local24[,3],suma_local25[,3],col="gray",border = "darkblue")
b4=boxplot(suma_local1[,4],suma_local2[,4],suma_local3[,4],suma_local4[,4],suma_local5[,4],suma_local6[,4],suma_local7[,4],suma_local8[,4],suma_local9[,4],suma_local10[,4],suma_local11[,4],suma_local12[,4],suma_local13[,4],suma_local14[,4],suma_local15[,4],suma_local16[,4],suma_local17[,4],suma_local18[,4],suma_local19[,4],suma_local20[,4],suma_local21[,4],suma_local22[,4],suma_local23[,4],suma_local24[,4],suma_local25[,4])
b5=boxplot(suma_local1[,5],suma_local2[,5],suma_local3[,5],suma_local4[,5],suma_local5[,5],suma_local6[,5],suma_local7[,5],suma_local8[,5],suma_local9[,5],suma_local10[,5],suma_local11[,5],suma_local12[,5],suma_local13[,5],suma_local14[,5],suma_local15[,5],suma_local16[,5],suma_local17[,5],suma_local18[,5],suma_local19[,5],suma_local20[,5],suma_local21[,5],suma_local22[,5],suma_local23[,5],suma_local24[,5],suma_local25[,5],col="gray",border = "darkblue")
b6=boxplot(suma_local1[,6],suma_local2[,6],suma_local3[,6],suma_local4[,6],suma_local5[,6],suma_local6[,6],suma_local7[,6],suma_local8[,6],suma_local9[,6],suma_local10[,6],suma_local11[,6],suma_local12[,6],suma_local13[,6],suma_local14[,6],suma_local15[,6],suma_local16[,6],suma_local17[,6],suma_local18[,6],suma_local19[,6],suma_local20[,6],suma_local21[,6],suma_local22[,6],suma_local23[,6],suma_local24[,6],suma_local25[,6],col="gray",border = "darkblue")



par1=mean(c(suma_local1[,1],suma_local2[,1],suma_local3[,1],suma_local4[,1],suma_local5[,1],suma_local6[,1],suma_local7[,1],suma_local8[,1],suma_local9[,1],suma_local10[,1],suma_local11[,1],suma_local12[,1],suma_local13[,1],suma_local14[,1],suma_local15[,1],suma_local16[,1],suma_local17[,1],suma_local18[,1],suma_local19[,1],suma_local20[,1],suma_local21[,1],suma_local22[,1],suma_local23[,1],suma_local24[,1],suma_local25[,1]))
par3=mean(c(suma_local1[,3],suma_local2[,3],suma_local3[,3],suma_local4[,3],suma_local5[,3],suma_local6[,3],suma_local7[,3],suma_local8[,3],suma_local9[,3],suma_local10[,3],suma_local11[,3],suma_local12[,3],suma_local13[,3],suma_local14[,3],suma_local15[,3],suma_local16[,3],suma_local17[,3],suma_local18[,3],suma_local19[,3],suma_local20[,3],suma_local21[,3],suma_local22[,3],suma_local23[,3],suma_local24[,3],suma_local25[,3]))
par5=mean(c(suma_local1[,5],suma_local2[,5],suma_local3[,5],suma_local4[,5],suma_local5[,5],suma_local6[,5],suma_local7[,5],suma_local8[,5],suma_local9[,5],suma_local10[,5],suma_local11[,5],suma_local12[,5],suma_local13[,5],suma_local14[,5],suma_local15[,5],suma_local16[,5],suma_local17[,5],suma_local18[,5],suma_local19[,5],suma_local20[,5],suma_local21[,5],suma_local22[,5],suma_local23[,5],suma_local24[,5],suma_local25[,5]))
par6=mean(c(suma_local1[,6],suma_local2[,6],suma_local3[,6],suma_local4[,6],suma_local5[,6],suma_local6[,6],suma_local7[,6],suma_local8[,6],suma_local9[,6],suma_local10[,6],suma_local11[,6],suma_local12[,6],suma_local13[,6],suma_local14[,6],suma_local15[,6],suma_local16[,6],suma_local17[,6],suma_local18[,6],suma_local19[,6],suma_local20[,6],suma_local21[,6],suma_local22[,6],suma_local23[,6],suma_local24[,6],suma_local25[,6]))
rho_mean=c(par1,par3,par5,par6)

list1=list(suma_local1[,1],suma_local2[,1],suma_local3[,1],suma_local4[,1],suma_local5[,1],suma_local6[,1],suma_local7[,1],suma_local8[,1],suma_local9[,1],suma_local10[,1],suma_local11[,1],suma_local12[,1],suma_local13[,1],suma_local14[,1],suma_local15[,1],suma_local16[,1],suma_local17[,1],suma_local18[,1],suma_local19[,1],suma_local20[,1],suma_local21[,1],suma_local22[,1],suma_local23[,1],suma_local24[,1],suma_local25[,1])
list3=list(suma_local1[,3],suma_local2[,3],suma_local3[,3],suma_local4[,3],suma_local5[,3],suma_local6[,3],suma_local7[,3],suma_local8[,3],suma_local9[,3],suma_local10[,3],suma_local11[,3],suma_local12[,3],suma_local13[,3],suma_local14[,3],suma_local15[,3],suma_local16[,3],suma_local17[,3],suma_local18[,3],suma_local19[,3],suma_local20[,3],suma_local21[,3],suma_local22[,3],suma_local23[,3],suma_local24[,3],suma_local25[,3])
list5=list(suma_local1[,5],suma_local2[,5],suma_local3[,5],suma_local4[,5],suma_local5[,5],suma_local6[,5],suma_local7[,5],suma_local8[,5],suma_local9[,5],suma_local10[,5],suma_local11[,5],suma_local12[,5],suma_local13[,5],suma_local14[,5],suma_local15[,5],suma_local16[,5],suma_local17[,5],suma_local18[,5],suma_local19[,5],suma_local20[,5],suma_local21[,5],suma_local22[,5],suma_local23[,5],suma_local24[,5],suma_local25[,5])
list6=list(suma_local1[,6],suma_local2[,6],suma_local3[,6],suma_local4[,6],suma_local5[,6],suma_local6[,6],suma_local7[,6],suma_local8[,6],suma_local9[,6],suma_local10[,6],suma_local11[,6],suma_local12[,6],suma_local13[,6],suma_local14[,6],suma_local15[,6],suma_local16[,6],suma_local17[,6],suma_local18[,6],suma_local19[,6],suma_local20[,6],suma_local21[,6],suma_local22[,6],suma_local23[,6],suma_local24[,6],suma_local25[,6])

mean1=rep(0,25)
mean3=rep(0,25)
mean5=rep(0,25)
mean6=rep(0,25)

for (i in 1:25){
  mean1[i]=mean(list1[[i]])
  mean3[i]=mean(list3[[i]])
  mean5[i]=mean(list5[[i]])
  mean6[i]=mean(list6[[i]])
}

wend_local=matrix(rep(0,25*81),ncol=81)
for (i in 1:25){
  wend_local[i,]=rho_wend(h,4,mean3[i],mean5[i],mean6[i],mean1[i])
}
prom_wend=rep(0,81)
for (i in 1:81){
  prom_wend[i]=mean(wend_local[,i])
  
}

###### prom_wend es el SCCC promedio

############## CREATING THE WENDLAND FUNCTION AND PLOTS
#####SETTING VARIANCES
sig1=sig2=1
rho12=0.3
sig12=sqrt(sig1*sig2)*rho12

sig1est=rho_mean[3]^2
sig2est=rho_mean[4]^2
rho12est=rho_mean[1]
sig12est=sqrt(sig1est*sig2est)*rho12est

##################### WENDLAND-GNEITING##############################
##### SETTING DISTANCE VECTOR#############################
h=seq(0,2,0.025)
#### SETTING PARAMETERS#####################################
kappa=1
nu=4
mu=5/2
b12=1.5
b12est=rho_mean[2]
########### FUNCION WENDLAND
indicator<-function(h,b12) ifelse(h>=0 & h<=b12 ,1,0)
funwend=function(h,nu,b12){
  (1+(nu+1)*h/b12)*(1-h/b12)^(nu+1)*indicator(h,b12)
}
#########WEND

rho_wend=function(h,nu,b12,sig1,sig2,rho12){
  w1=funwend(h,nu,b12)
  sig12=(sig1*sig2)*rho12
  rhow1=2*sig12*w1/(sig1^2+sig2^2)
  return (rhow1)
}


par(pty="s")
plot(h,rho_wend(h,4,1.5,1,1,0.3),type="l",xlab=expression(group("|",group("|",h,"|"),"|")),ylab="",cex.axis=2,cex.lab=1.5,ylim=c(-0.2,0.45))
lines(h,prom_wend,type="p",pch=1,col="blue")
lines(h,rho_wend(h,4,par3,par5,par6,par1),type="p",pch=4,col="red")
legend(0.4,0.5, 
       legend = expression(rho^c ~ (h)),
       lty=1, lwd=3,bty="n",cex = 2,pt.cex = 1, col="black",seg.len=1,text.col = "black")
legend(0.7,0.4, 
       legend = expression(hat(rho)[1] ~~ (h)),
       lty=0,pch=1, lwd=2,bty="n",cex = 2,pt.cex = 1, col="blue",seg.len = 1,text.col = "blue")
legend(0.7,0.3, 
       legend = expression(hat(rho)[2] ~~ (h)),
       lty=0,pch=4, lwd=2,bty="n",cex = 2,pt.cex = 1, col="red",seg.len = 1,text.col = "red")
text(1.74,0.23, expression(c),col="blue",cex=1.5)
text(1.74,0.13, expression(c),col="red",cex=1.5)
abline(h=0)



percen1=as.numeric(c(quantile(suma_caso1[,1],c(0.05,0.95)),quantile(suma_caso1[,3],c(0.05,0.95)),quantile(suma_caso1[,5],c(0.05,0.95)),quantile(suma_caso1[,6],c(0.05,0.95))))
percen2=as.numeric(c(quantile(suma_caso2[,1],c(0.05,0.95)),quantile(suma_caso2[,3],c(0.05,0.95)),quantile(suma_caso2[,5],c(0.05,0.95)),quantile(suma_caso2[,6],c(0.05,0.95))))
percen3=as.numeric(c(quantile(suma[,1],c(0.05,0.95)),quantile(suma[,3],c(0.05,0.95)),quantile(suma[,5],c(0.05,0.95)),quantile(suma[,6],c(0.05,0.95))))
percen4=as.numeric(c(quantile(suma_caso4[,1],c(0.05,0.95)),quantile(suma_caso4[,3],c(0.05,0.95)),quantile(suma_caso4[,5],c(0.05,0.95)),quantile(suma_caso4[,6],c(0.05,0.95))))
mean_caso1=c(mean(suma_caso1[,1]),mean(suma_caso1[,3]),mean(suma_caso1[,5]),mean(suma_caso1[,6]))
mean_caso2=c(mean(suma_caso2[,1]),mean(suma_caso2[,3]),mean(suma_caso2[,5]),mean(suma_caso2[,6]))
mean_caso3=c(mean(suma[,1]),mean(suma[,3]),mean(suma[,5]),mean(suma[,6]))
mean_caso4=c(mean(suma_caso4[,1]),mean(suma_caso4[,3]),mean(suma_caso4[,5]),mean(suma_caso4[,6]))


###CASO1
par(pty="s")


rho1_teo=rho_wend(h,4,0.35,1,1,-0.15)
rho1_fit=rho_wend(h,4,mean_caso1[2],mean_caso1[3],mean_caso1[4],mean_caso1[1])
rho1_lower=rho_wend(h,4,percen1[3],percen1[5],percen1[7],percen1[1])
rho1_upper=rho_wend(h,4,percen1[4],percen1[6],percen1[8],percen1[2])

plot(h,rho1_teo,type="l",xlab=expression(group("|",group("|",h,"|"),"|")),ylab="",cex.axis=2,cex.lab=1.5,ylim=c(-0.2,0.45))
abline(h=0)
polygon(c(h,rev(h)),c(rho1_lower,rev(rho1_upper)),col = "grey75", border = FALSE)
lines(h,rho1_teo,lwd=2)
lines(h,rho1_upper,lty=2,col="red")
lines(h,rho1_lower,lty=2,col="red")
lines(h,rho1_fit,type="p",pch=1,col="blue")
legend(0.4,0.5, 
       legend = expression(rho^c ~ (h)),
       lty=1, lwd=3,bty="n",cex = 2,pt.cex = 1, col="black",seg.len=1,text.col = "black")
legend(0.7,0.4, 
       legend = expression(hat(rho) ~~ (h)),
       lty=0,pch=1, lwd=2,bty="n",cex = 2,pt.cex = 1, col="blue",seg.len = 1,text.col = "blue")
text(1.74,0.23, expression(c),col="blue",cex=1.5)


###########CASO2

rho2_teo=rho_wend(h,4,0.9,1,1,0.25)
rho2_fit=rho_wend(h,4,mean_caso2[2],mean_caso2[3],mean_caso2[4],mean_caso2[1])
rho2_lower=rho_wend(h,4,percen2[3],percen2[5],percen2[7],percen2[1])
rho2_upper=rho_wend(h,4,percen2[4],percen2[6],percen2[8],percen2[2])

plot(h,rho2_teo,type="l",xlab=expression(group("|",group("|",h,"|"),"|")),ylab="",cex.axis=2,cex.lab=1.5,ylim=c(-0.2,0.45))
abline(h=0)
polygon(c(h,rev(h)),c(rho2_lower,rev(rho2_upper)),col = "grey75", border = FALSE)
lines(h,rho2_teo,lwd=2)
lines(h,rho2_upper,lty=2,col="red")
lines(h,rho2_lower,lty=2,col="red")
lines(h,rho2_fit,type="p",pch=1,col="blue")
legend(0.4,0.5, 
       legend = expression(rho^c ~ (h)),
       lty=1, lwd=3,bty="n",cex = 2,pt.cex = 1, col="black",seg.len=1,text.col = "black")
legend(0.7,0.4, 
       legend = expression(hat(rho) ~~ (h)),
       lty=0,pch=1, lwd=2,bty="n",cex = 2,pt.cex = 1, col="blue",seg.len = 1,text.col = "blue")
text(1.74,0.23, expression(c),col="blue",cex=1.5)

###########CASO3

rho3_teo=rho_wend(h,4,1.4,1,1,0.3)
rho3_fit=rho_wend(h,4,mean_caso3[2],mean_caso3[3],mean_caso3[4],mean_caso3[1])
rho3_lower=rho_wend(h,4,percen3[3],percen3[5],percen3[7],percen3[1])
rho3_upper=rho_wend(h,4,percen3[4],percen3[6],percen3[8],percen3[2])

plot(h,rho3_teo,type="l",xlab=expression(group("|",group("|",h,"|"),"|")),ylab="",cex.axis=2,cex.lab=1.5,ylim=c(-0.2,0.45))
abline(h=0)
polygon(c(h,rev(h)),c(rho3_lower,rev(rho3_upper)),col = "grey75", border = FALSE)
lines(h,rho3_teo,lwd=2)
lines(h,rho3_upper,lty=2,col="red")
lines(h,rho3_lower,lty=2,col="red")
lines(h,rho3_fit,type="p",pch=1,col="blue")
legend(0.4,0.5, 
       legend = expression(rho^c ~ (h)),
       lty=1, lwd=3,bty="n",cex = 2,pt.cex = 1, col="black",seg.len=1,text.col = "black")
legend(0.7,0.4, 
       legend = expression(hat(rho) ~~ (h)),
       lty=0,pch=1, lwd=2,bty="n",cex = 2,pt.cex = 1, col="blue",seg.len = 1,text.col = "blue")
text(1.74,0.23, expression(c),col="blue",cex=1.5)

###########CASO4

rho4_teo=rho_wend(h,4,0.9,2,0.5,0.45)
rho4_fit=rho_wend(h,4,mean_caso4[2],mean_caso4[3],mean_caso4[4],mean_caso4[1])
rho4_lower=rho_wend(h,4,percen4[3],percen4[5],percen4[7],percen4[1])
rho4_upper=rho_wend(h,4,percen4[4],percen4[6],percen4[8],percen4[2])

plot(h,rho4_teo,type="l",xlab=expression(group("|",group("|",h,"|"),"|")),ylab="",cex.axis=2,cex.lab=1.5,ylim=c(-0.2,0.45))
abline(h=0)
polygon(c(h,rev(h)),c(rho4_lower,rev(rho4_upper)),col = "grey75", border = FALSE)
lines(h,rho4_teo,lwd=2)
lines(h,rho4_upper,lty=2,col="red")
lines(h,rho4_lower,lty=2,col="red")
lines(h,rho4_fit,type="p",pch=1,col="blue")
legend(0.4,0.5, 
       legend = expression(rho^c ~ (h)),
       lty=1, lwd=3,bty="n",cex = 2,pt.cex = 1, col="black",seg.len=1,text.col = "black")
legend(0.7,0.4, 
       legend = expression(hat(rho) ~~ (h)),
       lty=0,pch=1, lwd=2,bty="n",cex = 2,pt.cex = 1, col="blue",seg.len = 1,text.col = "blue")
text(1.74,0.23, expression(c),col="blue",cex=1.5)