#Notation 

EhZ <- 0.379 #The expectation of h(Z), where Z follows the standard normal distribution and h(x)=1/(x^2+2)

hsup <- (3*sqrt(1.5))/16 # ||h'|| for h(x) as above and with ||.|| denoting the supremum norm.

boundexp <- function(s){(hsup/sqrt(s))*(4.41456) + 4*(s+2)/((s-1)*(s-2)) + 8*hsup*sqrt(s)*(s+2)/((s-1)*(s-2))}
#The upper bound for the exponential distribution in the canonical parameterisation.

boundexpnon <- function(s){(hsup/sqrt(s))*(4.41456) + 4/s + 2*hsup/sqrt(s) + 80*sqrt(3)*sqrt(2/s+1)*hsup/sqrt(s)}
#The upper bound for the exponential distribution in the non-canonical case.

boundexpnonL <- function(s){(hsup/sqrt(s))*(4.41456)}
#The upper bound for Exp(1/theta) using directly Lemma 1.1; see also (3) of Remark 3.3.

#Simulation results for the exponential distribution in the canonical case
## here the choice is rate = 1 and n=10 
n=10
theta=1
ydataexp10 <- matrix(rep(NA,10000*n),n,10000) #creating the matrix that has all the simulated data. 10,000 trials of n
for (i in 1:10000) {                          #observations each time.      
  ydataexp10[,i] <- rexp(n,theta)}
MLE <- rep(0,10000)  
for(i in 1:10000){
  MLE[i] <- 1/mean(ydataexp10[,i])} #Creating the 10,000 MLEs. For this example, the MLe is the reciprocal of the sample mean.
expstandardized10 <- (MLE-theta)*sqrt(n)/theta # The standardised MLE.
exphstandardized10 <- 1/((expstandardized10)^2+2) #Apply function h on the standardised MLE.
result <- mean(exphstandardized10)
result2 <- abs(result-EhZ) #The value of the quantity of interest for which an upper bound is desired.
result3 <- boundexp(n) #The value of the upper bound for the specific value of the sample size.
finalresult <- list(hdtand = result, absolutediff=result2, bound=result3)
finalresult

#Other values of theta and other sample sizes are treated in the same way.

## As the code for the exponential distribution in the non-canonical case is similar.

## scale = 2 n=10
n=10
theta=2
ydataexpnon10 <- matrix(rep(NA,10000*n),n,10000)
for (i in 1:10000) {
  ydataexpnon10[,i] <- rexp(n,1/theta)}
MLE <- rep(0,10000)
for(i in 1:10000){
  MLE[i] <- mean(ydataexpnon10[,i])}
expnonstandardized10 <- (MLE-theta)*sqrt(n)/theta
expnonhstandardized10 <- 1/((expnonstandardized10)^2+2)
result <- mean(expnonhstandardized10)
result2 <- abs(result-EhZ)
result3 <- boundexpnon(n)
result4 <- boundexpnonL(10)
finalresult <- list(hdtand = result, absolutediff=result2, bound=result3, boundLemma=result4)
finalresult


# Code for the simulations related to the MSE of the Beta distribution. There is a slight
# calculation error for the Psi functions due to R precision.

## For beta=1
sample_beta <- function(theta){
  psi1theta = sum((1/(theta+0:2^27))^2) #Calculating psi_1(theta)
  psi1theta1 = sum((1/(theta+1+0:2^27))^2) #Calculating psi_1(theta+beta) when beta=1
  psid = (1/theta)^2 #The difference between psi_1(theta) and psi_1(theta+1) is (1/theta)^2
  B2 = (96 + 6.6*theta^4)/(theta)^4 #The value B2 as defined in the paper.
  min_s = ceiling((B2*theta/2 + sqrt((B2*theta/2)^2 + 8*psid^2))^2/(theta^2*(psid)^3))
  min_s #this is the minimum value of the sample size such that assumption (Fur.3) in the paper is satisfied.
}

ms<-sample_beta(1.5) ##for theta=1.5 then n>=7460
B1 <- function(theta){
  psi1theta = sum((1/(theta+0:2^27))^2)
  psi1theta1 = sum((1/(theta+1+0:2^27))^2)
  psi3theta = 6*sum((1/(theta+0:2^27))^4)
  psi3theta1 = psi_3 = 6*sum((1/(theta+1+0:2^27))^4)
  result = 8*(psi3theta + psi3theta1 + 3*(psi1theta^2 + psi1theta1^2))
  result
} ##Here is the function B1(theta) as given in Corollary 5.1.

b1=B1(1.5) ##B1 for theta=1.5
final <- function(theta){
  thetahat = matrix(NA,10000,1000) #Creating the matrix that will include the data for the MLE
  squarerootnMSEhat = rep(NA,1000) #Creating the vector that contains the 1000 values of the sample MSE.
  bound = rep(NA,1000) #Creating the vector for the bounds
  diff=rep(NA,1000) #Creating the vector for the difference between the bound and the sample MSE value.
  psi1theta = sum((1/(theta+0:2^27))^2)
  psi1theta1 = sum((1/(theta+1+0:2^27))^2)
  psid = (1/theta)^2
  psi3theta = 6*sum((1/(theta+0:2^27))^4) #Calculating psi_3(theta)
  psi3theta1 = 6*sum((1/(theta+1+0:2^27))^4) #Calculating psi_3(theta+beta) for beta=1
  B2 = (96 + 6.6*theta^4)/((theta)^4) #B2 as defined in Corollary 5.1.
  for (i in ms:(ms+999)){
    for(j in 1:10000){
      sample = rbeta(i,theta,1) 
      thetahat[j, i-ms+1] = -i/sum(log(sample))} # This 'for' loop creates 10000 values of the MLE for each sample size
                                                 # from 7460 up to 8459
    
    squarerootnMSEhat[i-ms+1] = (1/10000)*sum((thetahat[,i-ms+1]-theta)^2) #The 1000 sample MSEs corresponding to each
                                                                           #sample size value.
    
    bound[i-ms+1] = (1/i)*((-sqrt(-4*(1 + (2/sqrt(i))*(2+((b1^(0.75))/(psid^(3/2)))))
                                  *((8/(i*theta^2*psid))+(B2/(sqrt(i)*(psid)^(3/2)))-1)))
                           /(2*((8/(i*theta^2*sqrt(psid))) + (B2/(sqrt(i)*psid)) -sqrt(psid))))^2 #Calculating the bound
                                                                                      #as defined in (5.5) of the paper.
    
    diff[i-ms+1] = bound[i-ms+1]-squarerootnMSEhat[i-ms+1] #The difference between the bound and the sample MSE.
  }
  result=list(thetahat=thetahat, nMSEhat = squarerootnMSEhat, bound=bound,difference=diff)
  result
}
#The results when theta=1.5
result = final(1.5)
result$nMSEhat
result$bound
result$difference
########################