


#R FUNCTIONS REQUIRED FOR the simulation table
#R functions documented:
#DRAWsi
#DRAWsisi2d
#DRAWsisiccs
#ESTratio 
#ESTtotal 
#ESTVARlinratio2d
#ESTVARlinratioccs
#ESTVARtotal2d
#ESTVARtotalccs
#LINratio
#MEASURESratio2d
#MEASURESratioccs
#POPratio 
#POPtotal 
#SAMPLEsisiccs
#SIMULATIONratio2d
#SIMULATIONratioccs
#VARIANCE_MCratio2d
#VARIANCE_MCratioccs


#Note replicability: the seed was fixed using the command set.seed() in the 
#functions SIMULATIONratio2d, SIMULATIONratioccs, VARIANCE_MCratio2d and VARIANCE_MCratioccs.


################################################################################################
#DESCRIPION: draws (and fixes) simple sample (equal probabilities) without replacement
#USAGE: DRAWsi(n,N) 
#ARGUMENTS: 
#		n 	size of the sample
#		N	size of the population
#VALUE: the function returns a vector of lenght n

DRAWsi <- function(n,N) 
	{ s=c(); s=sample(1:N,n); return(s) }


################################################################################################
#DESCRIPION: draws (and fixes) two-stage sample (equal probabilities at the two degrees) without replacement and with same sample size in each selected PSU
#USAGE: DRAWsisi2d(n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a vector of lenght n_m*n_d

DRAWsisi2d <- function(n_m,n_d,N_m,N_d)
	{
	Vs_d <- c()
	s_m=sort(DRAWsi(n_m,N_m)) 
	for (i in 1:n_m)
		{	
		s_d=sort(DRAWsi(n_d,N_d)) + N_d*(s_m[i]-1)
		Vs_d=c(Vs_d,s_d)
		}
	Dummy=rep(0,N_m*N_d) ; Dummy[Vs_d] <-1 ; return(Dummy)
	}


################################################################################################
#DESCRIPION: draws (and fixes) two simple samples (equal probabilities) without replacement in two populations
#USAGE: DRAWsisiccs(n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of two vectors of respective lenghts n_m and n_d

DRAWsisiccs <- function(n_m,n_d,N_m,N_d)
	{
	sd=DRAWsi(n_d,N_d); sm=DRAWsi(n_m,N_m)
	S=list(sm,sd) ; return(S) 
	}


################################################################################################
#DESCRIPION: estimates a population total
#USAGE: ESTtotal(ECH,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECH	matrix of the interest variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the estimated total

ESTtotal <- function(ECH,n_m,n_d,N_m,N_d)
	{w=(N_m/n_m)*(N_d/n_d); Tpi=sum(w*c(ECH)) ; return(Tpi) }


################################################################################################
#DESCRIPION: estimates a population ratio Y/X
#USAGE: ESTratio(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the estimated ratio

ESTratio <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{Rpi=ESTtotal(ECHY,n_m,n_d,N_m,N_d)/ESTtotal(ECHX,n_m,n_d,N_m,N_d) ; return(Rpi) }


################################################################################################
#DESCRIPION: estimates the variance of the estimated total of the estimated linearized variable of a ratio Y/X
#USAGE: ESTVARlinratio2d(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of two values (the unbiased estimator V2d and its approximations V2d,a)

ESTVARlinratio2d <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{
	Rlin=LINratio(ECHY,ECHX,n_m,n_d,N_m,N_d)
	vl=ESTVARtotal2d(Rlin,n_m,n_d,N_m,N_d) ; return(vl)
	}


################################################################################################
#DESCRIPION: estimates the variance of the estimated total of the estimated linearized variable of a ratio Y/X
#USAGE: ESTVARlinratioccs(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of three values (the unbiased estimator Vccs and its approximations V2d and V2d,a)

ESTVARlinratioccs <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{
	Rlin=LINratio(ECHY,ECHX,n_m,n_d,N_m,N_d)
	vl=ESTVARtotalccs(Rlin,n_m,n_d,N_m,N_d) ; return(vl)
	}


################################################################################################
#DESCRIPION: estimates the variance of an estimated total for a two-stage sampling {SI,SI}
#USAGE: ESTVARtotal2d(ECH,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECH	matrix of the interest variable, of size  [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of two values (the unbiased estimator V2d and its approximation V2d,a)

ESTVARtotal2d <- function(ECH,n_m,n_d,N_m,N_d)
	{
	w_m=N_m/n_m; w_d=N_d/n_d
	VarM=var(apply(ECH*w_d,2,sum));  V2da=(N_m^2)*(1/n_m - 1/N_m)*VarM	
	Vard=(N_d^2)*(1/n_d - 1/N_d)*apply(ECH,2,var); V2db=w_m*sum(Vard)
	V2d=V2da+V2db
	v=list(V2d,V2da) ; return(v)
	}


################################################################################################
#DESCRIPION: estimates the variance of an estimated total for a CCS SIxSI
#USAGE: ESTVARtotalccs(ECH,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECH	matrix of the interest variable, of size  [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of three values (the unbiased estimator Vccs and its approximations V2d and V2d,a)

ESTVARtotalccs <- function(ECH,n_m,n_d,N_m,N_d)
	{
	w_m=N_m/n_m; w_d=N_d/n_d
	VarM=var(apply(ECH*w_d,2,sum));  V2da=(N_m^2)*(1/n_m - 1/N_m)*VarM	
	VarD=var(apply(ECH*w_m,1,sum));  Vsimp=(N_d^2)*(1/n_d - 1/N_d)*VarD
	d=apply(ECH,1,mean); e=apply(ECH,2,mean); f=mean(ECH)
	X = t(t(ECH - d) - e) + f; X2=X*X
	Vco=(1/(n_d - 1))*(1/(n_m - 1))*sum(X2)
	VCORR=(N_d^2)*(1/n_d - 1/N_d)*(N_m^2)*(1/n_m - 1/N_m)*Vco	
	VCCS=V2da+Vsimp-VCORR
	Vard=(N_d^2)*(1/n_d - 1/N_d)*apply(ECH,2,var); V2db=w_m*sum(Vard)
	V2d=V2da+V2db
	v=list(VCCS,V2d,V2da) ; return(v)
	}


################################################################################################
#DESCRIPION: estimates the linearized variable of a ratio Y/X
#USAGE: LINratio(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a matrix of length [n_d,n_m]

LINratio <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{
	rpi=ESTratio(ECHY,ECHX,n_m,n_d,N_m,N_d)
	txpi=ESTtotal(ECHX,n_m,n_d,N_m,N_d)
	upi = (ECHY - ECHX*rpi) / txpi ; return(upi)
	}


################################################################################################
#DESCRIPION: measures (with simulations in Population) the relative biases of an estimated ratio Y/X and of its variance estimators for two-stage sampling
#USAGE: MEASURESratio2d(NBech,NBechMC,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBech		number of simulations to calculate the estimated ratio and its variance estimators
#		NBechMC	number of simulations to calculate the 'true' variance
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of 12 values: the sizes n_m, N_m, n_d, N_d, 
#the ratio population, its estimate under NBech simulations, its relative bias, 
#the variance monteCarlo under NBechMC simulation, the unbiased estimated variance under NBech simulations, its relative bias, 
#the simplified variance estimator V2d,a under NBech simulations, its relative bias

MEASURESratio2d <- function(NBech,NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	esp=SIMULATIONratio2d(NBech,POPY,POPX,n_m,n_d,N_m,N_d)
	varmc=VARIANCE_MCratio2d(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	t=POPratio(POPY,POPX)
	RBratio=((esp[[1]]-t)/t) *100
	RBvar=((esp[[2]]-varmc)/varmc) *100
	RBvar1=((esp[[3]]-varmc)/varmc) *100
	m=list(n_m,N_m,n_d,N_d,t,esp[[1]],RBratio,varmc,esp[[2]],RBvar,esp[[3]],RBvar1)
	return(m)
	}


################################################################################################
#DESCRIPION: measures (with simulations in Population) the relative biases of an estimated ratio Y/X and of its variance estimators for CCS 
#USAGE: MEASURESratioccs(NBech,NBechMC,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBech		number of simulations to calculate the estimated ratio and its variance estimators
#		NBechMC	number of simulations to calculate the 'true' variance
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of 12 values: the sizes n_m, N_m, n_d, N_d, 
#the ratio population, its estimate under NBech simulations, its relative bias, 
#the variance monteCarlo under NBechMC simulation, the unbiased estimated variance under NBech simulations, its relative bias, 
#the simplified variance estimator V2d under NBech simulations, its relative bias,
#the simplified variance estimator V2d,a under NBech simulations, its relative bias

MEASURESratioccs <- function(NBech,NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	esp=SIMULATIONratioccs(NBech,POPY,POPX,n_m,n_d,N_m,N_d)
	varmc=VARIANCE_MCratioccs(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	t=POPratio(POPY,POPX)
	RBratio=((esp[[1]]-t)/t) *100
	RBvar=((esp[[2]]-varmc)/varmc) *100
	RBvar1=((esp[[3]]-varmc)/varmc) *100;RBvar2=((esp[[4]]-varmc)/varmc) *100
	m=list(n_m,N_m,n_d,N_d,t,esp[[1]],RBratio,varmc,esp[[2]],RBvar,esp[[3]],RBvar1,esp[[4]],RBvar2)
	return(m)
	}


################################################################################################
#DESCRIPION: calculates the population ratio Y/X
#USAGE: POPratio(POPY,POPX) 
#ARGUMENTS: 
#		POPY 	matrix of the numerator variable
#		POPX 	matrix of the denominator variable
#VALUE: the function returns a value for the population ratio

POPratio <- function(POPY,POPX) 
	{ return(sum(POPY)/sum(POPX)) }


################################################################################################
#DESCRIPION: calculates the population total
#USAGE: POPtotal(POP) 
#ARGUMENTS: 
#		POP 	matrix of the interest variable
#VALUE: the function returns a value for the population total

POPtotal <- function(POP) 
	{ return(sum(POP)) }


################################################################################################
#DESCRIPION: calculates the population ratio Y/X
#USAGE: SAMPLEsisiccs(POP,s_m,s_d)	 
#ARGUMENTS: 
#		POP 	matrix of the interest variable
#		s_m 	vector of lenght<dim(POP)[2]
#		s_d 	vector of lenght<dim(POP)[1]
#VALUE: the function returns a matrix of length [length(s_m),length(s_d)]

SAMPLEsisiccs <- function(POP,s_m,s_d)	
	{s=POP[s_d,s_m] ; return(s) }


################################################################################################
#DESCRIPION: simulations in Population for an estimated ratio Y/X and its variance estimators for a two-stage sampling {SI,SI}
#USAGE: SIMULATIONratio2d(NBech,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBech		number of simulations to calculate the estimated ratio and its variance estimators
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of 3 values:  
#the ratio estimator under NBech simulations, its unbiased variance estimator under NBech simulations,
#the simplified variance estimator V2d,a under NBech simulations 

SIMULATIONratio2d <- function(NBech,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	Est<-c(); EstVar<-c(); EstVar1<-c(); set.seed(12345)
	for (i in 1:NBech)		
		{
		echY<-c(); echX<-c(); Dummy_2=DRAWsisi2d(n_m,n_d,N_m,N_d)
		echY=matrix(POPY[Dummy_2==1],nrow=n_d)
		echX=matrix(POPX[Dummy_2==1],nrow=n_d)
		Est<-c(Est, ESTratio(echY,echX,n_m,n_d,N_m,N_d))
		v=ESTVARlinratio2d (echY,echX,n_m,n_d,N_m,N_d)
		EstVar<-c(EstVar,v[[1]]) 
		EstVar1<-c(EstVar1,v[[2]])
		}
	ESPest=mean(Est);ESPvar=mean(EstVar);ESPvar1=mean(EstVar1)
	E=list(ESPest,ESPvar,ESPvar1) ; return(E)
	}


################################################################################################
#DESCRIPION: simulations in Population for an estimated ratio Y/X and its variance estimators for a CCS SIxSI
#USAGE: SIMULATIONratioccs(NBech,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBech		number of simulations to calculate the estimated ratio and its variance estimators
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of 4 values:  
#the ratio estimator under NBech simulations, its unbiased variance estimator under NBech simulations,
#the simplified variance estimator V2d under NBech simulations 
#the simplified variance estimator V2d,a under NBech simulations

SIMULATIONratioccs <- function(NBech,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	Est<-c(); EstVar<-c(); EstVar1<-c(); EstVar2<-c() ; set.seed(1235)
	for (i in 1:NBech)		
		{
		s=DRAWsisiccs(n_m,n_d,N_m,N_d)
		POPYmat=matrix(POPY,nrow=N_d)
		POPXmat=matrix(POPX,nrow=N_d)
		echY=SAMPLEsisiccs(POPYmat,s[[1]],s[[2]])
		echX=SAMPLEsisiccs(POPXmat,s[[1]],s[[2]])
		Est<-c(Est, ESTratio(echY,echX,n_m,n_d,N_m,N_d))
		v=ESTVARlinratioccs (echY,echX,n_m,n_d,N_m,N_d)
		EstVar<-c(EstVar,v[[1]]) 
		EstVar1<-c(EstVar1,v[[2]]); EstVar2<-c(EstVar2,v[[3]])
		}
	ESPest=mean(Est);ESPvar=mean(EstVar);ESPvar1=mean(EstVar1);ESPvar2=mean(EstVar2)
	E=list(ESPest,ESPvar,ESPvar1,ESPvar2) ; return(E)
	}


################################################################################################
#DESCRIPION: calculates independently the MonteCarlo variance of an estimated ratio Y/X for a two-stage sampling {SI,SI}
#USAGE: VARIANCE_MCratio2d(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBechMC		number of simulations to calculate the MC variance
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the MC variance

VARIANCE_MCratio2d <- function(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	Est<-c() ; set.seed(13456)
	for (i in 1:NBechMC)		
		{
		echY<-c(); echX<-c() ; Dummy_2=DRAWsisi2d(n_m,n_d,N_m,N_d)
		echY=matrix(POPY[Dummy_2==1],nrow=n_d)
		echX=matrix(POPX[Dummy_2==1],nrow=n_d)
		Est<-c(Est, ESTratio(echY,echX,n_m,n_d,N_m,N_d))
		}
	return(var(Est))
	}


################################################################################################
#DESCRIPION: calculates independently the MonteCarlo variance of an estimated ratio Y/X for CCS SIxSI
#USAGE: VARIANCE_MCratioccs(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		NBechMC		number of simulations to calculate the MC variance
#		POPY	matrix of the numerator variable, of size [N_d,N_m]
#		POPX	matrix of the denominator variable, of size [N_d,N_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the MC variance

VARIANCE_MCratioccs <- function(NBechMC,POPY,POPX,n_m,n_d,N_m,N_d)
	{
	Est<-c() ; set.seed(1256)
	for (i in 1:NBechMC)		
		{
		s=DRAWsisiccs(n_m,n_d,N_m,N_d)
		POPYmat=matrix(POPY,nrow=N_d)
		POPXmat=matrix(POPX,nrow=N_d)
		echY=SAMPLEsisiccs(POPYmat,s[[1]],s[[2]])
		echX=SAMPLEsisiccs(POPXmat,s[[1]],s[[2]])
		Est<-c(Est, ESTratio(echY,echX,n_m,n_d,N_m,N_d))
		}
	return(var(Est))
	}




