library(np)
library(KernSmooth)

n = 100
x = rnorm(n)

m = function(x) {
  ifelse(x < 5, sin(x), sin(x) + 1)
}
sig = function(x) {
  0.1 + 0.05 * x
}

y = m(x) + sig(x) * rnorm(n)

x0 = seq(-3, 3, length.out = 100)
plot(x, y, pch = 20)
lines(x0, m(x0), lwd = 1.5)

h_RT = function(x, y) {
  R_K = 0.5 / sqrt(pi)
  # raw = T : obicni polinomi, ne ortogonalni
  mQ = lm(y ~ poly(x, raw = TRUE, degree = 4))
  mQ2 = 2 * mQ$coefficients[3] + 6 * mQ$coefficients[4] * x +
    12 * mQ$coefficients[5] * x^2
  theta_22_hat = mean((mQ2)^2)
  sigma2_hat = sum(mQ$residuals^2) / (length(x) - 5)
 
  ((R_K * (max(x) - min(x)) * sigma2_hat) / (theta_22_hat * length(x)))^(1 / 5)
  
}
(h_rt = h_RT(x,y))

# lc - local-constant (p=0), ll local-linear (p=1)
h_CV0 = npregbw(xdat = x, ydat = y, regtype = "lc")
(h_cv0 = h_CV0$bandwidth$x)
h_CV1 = npregbw(xdat = x, ydat = y, regtype = "ll")
(h_cv1 = h_CV1$bandwidth$x)

m0_cv = locpoly(x, y, degree = 0, bandwidth = h_cv0,
                gridsize = length(x0),
                range.x = range(x0))

m0_rt = locpoly(x, y, degree = 0, bandwidth = h_rt,
                   gridsize = length(x0),
                   range.x = range(x0))

m1_cv = locpoly(x, y, degree = 1, bandwidth = h_cv1,
                gridsize = length(x0),
                range.x = range(x0))

m1_rt = locpoly(x, y, degree = 1, bandwidth = h_rt,
                   gridsize = length(x0),
                   range.x = range(x0))

lines(x0, m0_rt$y, col = 2, lwd = 1.5)
lines(x0, m0_cv$y, col = 3, lwd = 1.5)
lines(x0, m1_rt$y, col = 4, lwd = 1.5)
lines(x0, m1_cv$y, col = 5, lwd = 1.5)
legend("bottomright", legend = c("true", "p=0, rt", "p=0, cv", "p=1, rt", "p=1, cv"),
       col=1:5, lty = 1, cex = 0.5)