#====================================================================EXAMPLE 1 (lognormal distribution)
bootstr<-function(){     
x<- rlnorm(NN)
for (a in 1:M){

sigma<-vector(length=B)
for (i in 1:B){
sigma[i]<-sum(sample(x,nn, replace="TRUE")^a)/nn
}
s_N<-sum(x^a)/NN

m<-mean(as.numeric(sigma>xi*s_N))     #p-value

if (m<0.05) break  #Break when the minimal moment that null hypothesis is rejected is found
}

a
}

set.seed(123456)
xi<-0.9

NN<-200     #Number of observations
M<-20000    #Maximal moment (Never achieved)
MN<-10000   #Number of Monte-Carlo simulations

nn<-round(log(NN))  #Size of bootstrap subsamples
nn

B<-5000              

AA<-vector(length=MN)

for (u in 1:MN){

AA[u]<-bootstr2()

hist(AA, main=u)
}


hist(sort(AA)[1:9900], main="", xlab="k") #Histogram
mean(AA==1)
mean(AA==2)   # Some results reported in the text
mean(AA==3)




#============================================================================EXAMPLE 2 (log-Pareto distribution)


library(VGAM)              #VGAM package may need to be installed

bootstr2<-function(){      
x1<-rpareto(NN, 1, 1)     #simple Pareto
for (a in 100:0){
a1<-a/100                 

x2<-a1*x1              #Taking log-Pareto-sample into the power 
                       #in order to avoid infinities in the next line as much as possible 
                       #few infinities with large "a" is not a problem, 
                       # since the test rejects the null hypothesis, in this case as it sholud do theoretically

x<-exp(x2)             #transformation to log-Preto

sigma<-vector(length=B)
for (i in 1:B){
sigma[i]<-sum(sample(x,nn, replace=TRUE))/nn  # As the sample is already taken to the power a few lines ago, 
                                                # here powers are not needed
}
s_N<-sum(x)/NN

m<-mean(as.numeric(sigma>xi*s_N))     #p-value

if (m>0.05) break  #Break when the maximal moment that nul hypothesis is accepted is found
}

a1
}

set.seed(123456)
xi<-0.9

NN<-200   # Number of observations
MN<-10000 # Number of Monte-Carlo simulations

nn<-round(log(NN)) #Size of bootstrap subsamples
nn

B<-5000                #Number of bootstraps 5000 - standard

AA<-vector(length=MN)

for (u in 1:MN){

AA[u]<-bootstr2()

hist(AA, main=u)
}


hist(sort(AA), main="", xlab="k") # Histogram
mean(AA==1)
mean(AA==0.01) # Some results reported in the text
mean(AA<=0.05)