pmethod <- function(theta=2){
    result <- c()
    for (j in seq(1,1000)){
        data <- c()
        for(i in seq(1,1000)){
            sample <- 2.07*rt(100,df=5) # sample from distribution where we estimated $theta$ which is 2.07
            est  <- sqrt((3/500)*sum(sample^2))
            data <- c(data,est)
        }
        qq <- quantile(data,probs=c(0.025,0.975))
        if(theta >= qq[1] & theta <= qq[2]){
            t <- 1
        }else{
            t <-0
        }
        result <- c(result,t)
    }
    return(mean(result))
}


pmethod()