Resampling

In this notebook we will sample random variables from a Pareto distribution with parameters \(\alpha\) and \(\beta\): \(Y_i \sim \mathcal{P}(\alpha,\beta)\) for \(i=1,\ldots,n\). Pareto random variables are continuous and are characterized by the probability density function: \[\begin{align} f(y|\alpha,\beta) & = \frac{\beta \alpha^{\beta}}{y^{\beta+1}}. \end{align}\] In all cases \(\alpha > 0\), \(\beta > 0\) and \(Y_i > \alpha\). Since there is no function in R to sample directly from a Pareto density, we can draw random variable \(Y\) via the inverse probability transform: \[\begin{align} U & \sim \mathcal{U}(0,1) \\ Y & \sim \left(\frac{\beta \alpha^{\beta}}{U}\right)^{\frac{1}{\beta+1}}. \end{align}\]

That is, a random draw, \(Y\), from the Pareto distribution is obtained by drawing a uniform random variable, \(U\), and applying the inverse of the Pareto density function to that uniform value. The following code accomplishes this procedure for 30 values.

# True parameters
alpha = 2
beta = 5
# Number of draws
n = 30#100000
# Draw a bunch of uniforms
U = runif(n)
# Apply to inverse of the Pareto 
y = alpha/((1-U)^(1/beta))

We can get a sense of the shape of the Pareto by plotting a density estimate of our random sample.

We will assume that \(\alpha\) = 2 is a known parameter but that \(\beta\) is unknown and must be estimated with data (that we simulate). Using the derivations for the Pareto likelihood in the lecture slides, we can compute and plot the likelihood for a variety of values of \(\beta\) and we can also compute the MLE and asymptotic variance (inverse of Fisher information) of \(\beta\). The vertical dotted line in the plot depicts the value of the MLE.

betaGrid = seq(0.1, 10, length=10000)
logLike = n*log(betaGrid) + n*betaGrid*log(alpha) - (1+betaGrid)*sum(log(y))
betaHat = n/(sum(log(y)) - n*(log(alpha)))
sdBetaHat = betaHat/sqrt(n)
plot(betaGrid, logLike, type='l')
abline(v=betaHat,lty=3)

We can also compute the MLE variance (or standard deviation) via the bootstrap, using 1000 bootstrap samples of size 30.

# Number of resamples
B = 1000
# Draw B resamples and store each sample in column of yMat
yMat = matrix(0, nrow=n, ncol=B)
for(i in 1:B) yMat[,i] = sample(y, n, replace=TRUE)
# Compute MLE for each resample (i.e. using each column of yMat)
betaHatBoot = rep(0, length=B)
for(i in 1:B){
    betaHatBoot[i] = n/(sum(log(yMat[,i])) - n*(log(alpha)))
}
muBetaHatBoot = mean(betaHatBoot)
sdBetaHatBoot = sd(betaHatBoot)

To finish, we compare the bootstrap standard deviation with asymptotic standard deviation of \(\hat{\beta}\):

print(sdBetaHat)

[1] 0.8934158

print(sdBetaHatBoot)

[1] 0.6328841

We can also compute and compare the three confidence intervals that we studied in lecture:

print(c(betaHat + qnorm(0.025)*sdBetaHat, betaHat + qnorm(0.975)*sdBetaHat))

[1] 3.142377 6.644503

print(c(betaHat + qnorm(0.025)*sdBetaHatBoot, betaHat + qnorm(0.975)*sdBetaHatBoot))

[1] 3.65301 6.13387

print(c(as.numeric(quantile(betaHatBoot,0.025)), as.numeric(quantile(betaHatBoot,0.975))))

[1] 3.949741 6.354060