plugin_bandwidth <- function(W, phi_U, sd_X, kernel_type) {
kernel_list <- kernel(kernel_type)
phiK <- kernel_list$phik
mu_K2 <- kernel_list$muk2
RK <- kernel_list$rk
tt <- kernel_list$tt
deltat <- tt[2] - tt[1]
n <- length(W)
# Grid of h on which to search for a solution
maxh <- (max(W) - min(W)) / 10
hnaive <- ((8 * sqrt(pi) * RK/3/mu_K2^2)^0.2) * sqrt(stats::var(W)) *
n^(-1/5)
lh <- 101
hgrid <- seq(hnaive / 3, maxh, length.out = lh)
# Quantities and functions that will be needed multiple times
toverh <- tt %*% t(1 / hgrid)
phiKsq <- phiK(tt)^2
phiUsq <- phi_U(toverh)^2
calculate_indh <- function(rr, th) {
term1 <- -hgrid^2 * mu_K2 * th
term2 <- kronecker(matrix(1, 1, lh), tt^(2*rr) * phiKsq) / phiUsq
term2 <- colSums(term2) * deltat / (2 * pi * n * hgrid^(2 * rr + 1))
ABias2 <- (term1 + term2)^2
indh <- which.min(ABias2)
if (indh == 1) {
warning(paste("Minimum of Abias2 for rr =", as.character(rr), "is
the first element of the grid of bandwidths. Consider enlarging
the grid."))
}
if(indh == lh){
warning(paste("Minimum of Abias2 for rr =", as.character(rr), "is
the last element of the grid of bandwidths. Consider enlarging
the grid."))
}
indh
}
calculate_th <- function(h, indh) {
phi_W <- ComputePhiEmp(W, tt/h)
th <- sum(tt^(2 * rr) * phi_W$norm^2 * phiKsq / phiUsq[, indh])
th * deltat / (2 * pi * h^(2 * rr + 1))
}
#---------------------------------------------------------------------------
# Start with theta4 and iterate down to get h1
#---------------------------------------------------------------------------
# theta 4
th4 <- (sd_X^(-9) * factorial(8) ) / (2^9 * factorial(4) * sqrt(pi))
# h3
rr <- 3
indh3 <- calculate_indh(rr, th4)
h3 <- hgrid[indh3]
# theta 3
th3 <- calculate_th(h3, indh3)
# h2
rr <- 2
indh2 <- calculate_indh(rr, th3)
h2 <- hgrid[indh2]
# theta 2
th2 <- calculate_th(h2, indh2)
# h1
term1 <- hgrid^4 * mu_K2^2 * th2 / 4
term2 <- kronecker(matrix(1,1,lh), phiKsq) / phiUsq
term2 <- colSums(term2) * deltat / (2 * pi * n * hgrid)
AMISE <- term1 + term2
indh1 <- which.min(AMISE)
hgrid[indh1]
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.