How to calculate Kernel Density for Bootstrap Likelihood

I am attempting to write R code to generate bootstrap likelihood as described in section 3 of this paper https://arxiv.org/pdf/1510.07287.pdf.

I am confident that I performed the bootstraps correct, but I am new to KDE and I am confused on how to aggregate the KDE densities to bootstrap likelihoods. I have attached a completely reproducible example of my attempt at bootstrap likelihood.

The steps they use in section 3 of the paper are as follows.

1. Generate K = 100 bootstraps of the data.

2. Calculate a theta (parameter of interest) for each of the K bootstraps

3. Generate L = 1000 second level bootstraps for each K boostraps

4. Calculate theta for each of the L bootstraps

5. Pass the theta calculated above to a kernel function to estimate the distribution

6. Log the KDE distribution to get the log likelihood of the K thetas

Using the mtcars data and the model $$mpg = \theta*cyl$$. I wanted to build a bootstrap likelihood for $$\theta$$. The point estimate from lm(mpg ~ cyl, data = mtcars)is about -2.9.

My likelihood function seems reasonable and is shown below: # install.packages("strapgod")
# install.packages("tidyverse")
# install.packages("KernSmooth")
library(strapgod)
library(tidyverse)
library(KernSmooth)
data(mtcars)

# check point estimate for cyl coeff
f <- lm(mpg ~ cyl, data = mtcars)
summary(f)

###########################
# Set up vars and functions
###########################
K = 100  # firest level boot
L = 1000 # second level boot

# model and pull out cyl coeff
lm_fun <- function(data){
fit <- lm(mpg ~ cyl, data = data)
summary(fit)\$coefficients
}

# kernel density functions
kern_funy <- function(params){
p <- params$$model est <- bkde(p, bandwidth=0.25) est$$y
}

kern_funx <- function(params){
p <- params$$model est <- bkde(p, bandwidth=0.25) est$$x
}

# Execute the method
# ------------------
set.seed(42)

# K bootstrap samples
P <- bootstrapify(mtcars, K) %>%
collect()

# generate theta_i estimates
thetai <- P %>%
group_by(.bootstrap) %>%
nest() %>%
mutate(theta = map_dbl(data, lm_fun))

# generate theta_ij estimates
thetaij <- P %>%
group_by(.bootstrap) %>%
bootstrapify(L) %>%
collect() %>%
group_by(.bootstrap, .bootstrap1) %>%
nest() %>%
mutate(model = map_dbl(data, lm_fun)) %>%
unnest(model)

# get kernel density from thetaij
kd <- thetaij %>%
select(.bootstrap, model) %>%
group_by(.bootstrap) %>%
nest() %>%
mutate(kernx = map(data, kern_funx)) %>%
mutate(kerny = map(data, kern_funy)) %>%
unnest(kernx, kerny) %>%
group_by(kernx) %>%
summarise(avg_p = mean(kerny))

# look up fun to find likelihood for a theta
get_p <- function(x){
kd %>%
filter(abs(kernx - x) == min(abs(kernx - x))) %>%
pull(avg_p)
}

# join params to kern
ll <- thetai %>%
rowwise() %>%
mutate(kerny = get_p(theta)) %>%
mutate(ll = log(kerny))

# Plot
ll %>%
ggplot(aes(x = theta, y = ll)) +
geom_point(shape = "o") +
stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1, se = F) +
ggtitle("log bootstrap likelihood") +
theme_classic()