#====================================================================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)