###############################################################################################
#                                                                                             #
#                 Trabajo 1: COEFICIENTE DE CORRELACIÓN (BOOTSTRAP)                           #
#                                                                                             #
###############################################################################################

#-----------------------------------------------------------------------------------------------#
#--1.-Caso 2: Y, Z variables indepenientes con distribución N(0,4)y N(2,1), respectivamente.    #
#Para n=20 y B=500, obtén una estimación de la distribución del coeficiente de correlación.     #
#Genera una muestra (X,Y) y aproxima el p-valor asociado al coeficiente de correlación muestral.#
#-----------------------------------------------------------------------------------------------#

B<-500
n<-20

y<-rnorm(n,0,4)
z<-rnorm(n,2,1)
coef_corrB<-numeric(B)

for(b in 1:B){
    yb<-sample(y,n,replace=TRUE)
    zb<-sample(z,n,replace=TRUE)
    covb<-mean(yb*zb)-mean(yb)*mean(zb)
    coef_corrB[b]<-covb/sqrt(var(yb)*var(zb))
}

Fn<-ecdf(coef_corrB)#Estimación de la función de distribución del coeficiente de correlación muestral.
plot(Fn,main="Función de distribución empírica del coeficiente de correlación muestral")
y0<-rnorm(n,0,4)
z0<-rnorm(n,2,1)
coef_corr_obs<-cov(y0,z0)/sqrt(var(y0)*var(z0))
p_valor=1-Fn(coef_corr_obs)

p_valor  #p-valor asociado al coeficiente de correlación muestral.


#---------------------------------------------------------------------------------------------#
#--2.-Caso 2: Y,Z variables independientes con distribución N(0,4) y N(0,1),respectivamente.  #
#Para n=20 y B=500, obtén una estimación de la distribución del coeficiente de correlación    # 
#muestral y del p-valor saociado a una muestra de (X,Y) con distribución N(0,4) y N(0,1),     #
#respectivamente,con correlación 0.25.                                                        #
#---------------------------------------------------------------------------------------------#

B<-500
n<-20
y<-rnorm(n,0,4)
z<-rnorm(n,0,1)
coef_corrB2<-numeric(B)

for(b in 1:B){
    yb<-sample(y,n,replace=TRUE)
    zb<-sample(z,n,replace=TRUE)
    covb<-mean(yb*zb)-mean(yb)*mean(zb)
    coef_corrB2[b]<-covb/sqrt(var(yb)*var(zb))

}

Fn<-ecdf(coef_corrB2)#Estimación de la función de distribución del coeficiente de correlación muestral.
plot(Fn,main="Función de distribución empírica del coeficiente de correlación muestral")
cov<-0.25*2
YZ<-rmvnorm(2*n,mean=c(0,0),sigma=matrix(c(4,cov,cov,1),ncol=2))
coef_corr_obs2<-cov(YZ[,1],YZ[,2])/sqrt(var(YZ[,1])*var(YZ[,2]))
p_valor=1-Fn(coef_corr_obs2)
p_valor  #p-valor asociado al coeficiente de correlación muestral.


#------------------------------------------------------------------------------------------------#
#---3.-En el caso 1, obten un intervalo de confianza Boostrap para el coeficiente de correlación.#
#con alpha=0.05 y evalúa su cobertura.                                                           #
#------------------------------------------------------------------------------------------------#

#a) Calcula el intervalo por Boostrap, con B=500 y alpha=0.05
#------------------------------------------------------------

n<-20
B<-500
alpha<-0.05

y<-rnorm(n,0,4)
z<-rnorm(n,2,1)
cov<-mean(y*z)-mean(y)*mean(z)
coef_corr<-cov/sqrt(var(y)*var(z))

coef_corrB<-numeric(B)
for(b in 1:B){
     yb<-sample(y,n,replace=TRUE)
     zb<-sample(z,n,replace=TRUE)
     covb<-mean(yb*zb)-mean(yb)*mean(zb)
     coef_corrB[b]<-covb/sqrt(var(yb)*var(zb))
}

scoef_corrB<-sort(coef_corrB)
clb<-scoef_corrB[floor((alpha/2)*B)]
cub<-scoef_corrB[floor((1-alpha/2)*B)]


IC<-c(coef_corr-cub,coef_corr-clb)
IC

par(mfrow=c(1,2))
hist(coef_corrB,freq=F,main="coef_corrB")
abline(v=clb,col=2,lwd=2)
abline(v=cub,col=2,lwd=2)
boxplot(coef_corrB,freq=F,main="coef_corrB")

#b)Covertura del intervalo, con 100 réplicas.  
#--------------------------------------------


n<-20
B<-500
M<-500
alpha<-0.05

coef_corrB<-numeric(B)
ICM<-matrix(0,nrow=M,ncol=2)


for(m in 1:M){
y<-rnorm(n,0,4)
z<-rnorm(n,0,1)
cov<-mean(y*z)-mean(y)*mean(z)
coef_corr<-cov/sqrt(var(y)*var(z))

for(b in 1:B){
     yb<-sample(y,n,replace=TRUE)
     zb<-sample(z,n,replace=TRUE)
     covb<-mean(yb*zb)-mean(yb)*mean(zb)
     coef_corrB[b]<-covb/sqrt(var(yb)*var(zb))
}
scoef_corrB<-sort(coef_corrB)
clb<-scoef_corrB[floor((alpha/2)*B)]
cub<-scoef_corrB[floor((1-alpha/2)*B)]

ICM[m,]<-c(coef_corr-cub,coef_corr-clb)
}
coverage<-sum((ICM[,1]<=0)*(ICM[,2]>=0))
porcentaje<-coverage/M

ind<-1:M
plot(ICM[,2]~ind,type="n",xlab="",ylab="",main="IC Bootstrap",ylim=c(min(ICM[,1]),max(ICM[,2])))
for(i in 1:M){
segments(ind[i],ICM[i,1],ind[i],ICM[i,2])
}
abline(h=0,col=2,lwd=2)
