## ----setup, include = FALSE-----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE,
                      fig.width = 6.5, fig.height = 3.6, fig.align = "center")
old <- options(width = 88)
set.seed(20260810)

## -------------------------------------------------------------------------------------
library(hdcce)

## -------------------------------------------------------------------------------------
data("data_estimation")
data("data_inference")

obs_N <- 20
obs_T <- 20
p     <- ncol(data_estimation$x)

dim(data_estimation$x)

## -------------------------------------------------------------------------------------
gsize     <- (p - 1) / 3
beta_true <- c(1, rep(c(1, 1, 1, rep(0, gsize - 3)), 3))
which(beta_true != 0)

## -------------------------------------------------------------------------------------
make_dict <- function(X) {
  Phi <- do.call(cbind, lapply(seq_len(ncol(X)), function(j) cbind(X[, j], X[, j]^2)))
  list(Phi = Phi, group = rep(seq_len(ncol(X)), each = 2))
}
dict_est <- make_dict(data_estimation$x)
dict_inf <- make_dict(data_inference$x)
dim(dict_est$Phi)

## -------------------------------------------------------------------------------------
fit <- hdcce_estimator(data_estimation, obs_N = obs_N, obs_T = obs_T,
                       NFOLDS = 5)
fit$K_hat

## ----scree----------------------------------------------------------------------------
plot(fit$eigenvalues[1:15], type = "b", pch = 16, ylim = c(0, 1),
     xlab = "index", ylab = "normalised eigenvalue",
     main = sprintf("K_hat = %d", fit$K_hat))
abline(h = 0.01, col = "red", lty = 2)
legend("topright", "TRUNC", lty = 2, col = "red", bty = "n")

## -------------------------------------------------------------------------------------
est <- as.numeric(fit$coefs)
round(head(est, 6), 3)

c(selected     = sum(est != 0),
  true_nonzero = sum(beta_true != 0),
  found        = sum(est[beta_true != 0] != 0))

## -------------------------------------------------------------------------------------
fit_d <- hdcce_estimator(data_estimation, obs_N = obs_N, obs_T = obs_T,
                         dictionaries = dict_est, NFOLDS = 5)
c(K_hat = fit_d$K_hat, n_coef = length(as.numeric(fit_d$coefs)))
round(head(fit_d$coefs, 6), 3)

## -------------------------------------------------------------------------------------
dat0 <- list(x = data_inference$x, y = data_inference$y[, 1])   # beta_1 = 0
inf0 <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T,
                        COEF_INDEX_VEC = 1, NFOLDS = 5)

r <- inf0$results[["1"]]
c(estimate = r$coef_despar, se = r$se, p_value = r$p_value)
r$confidence_band

## -------------------------------------------------------------------------------------
dat2 <- list(x = data_inference$x, y = data_inference$y[, 3])   # beta_1 = 0.2
inf2 <- hdcce_inference(dat2, obs_N = obs_N, obs_T = obs_T,
                        COEF_INDEX_VEC = 1, NFOLDS = 5)
inf2$results[["1"]]$confidence_band

## -------------------------------------------------------------------------------------
cb <- r$confidence_band
cbind(alpha         = inf0$alpha,
      p_le_alpha    = r$p_value <= inf0$alpha,
      zero_excluded = !(cb[, 1] <= 0 & 0 <= cb[, 2]))

## -------------------------------------------------------------------------------------
inf_multi <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T,
                             COEF_INDEX_VEC = c(1, 2, 5), NFOLDS = 5)
t(sapply(inf_multi$results,
         function(z) c(estimate = z$coef_despar, se = z$se, p = z$p_value)))

## -------------------------------------------------------------------------------------
sapply(1:3, function(k)
  hdcce_inference(dat2, obs_N, obs_T, COEF_INDEX_VEC = 1,
                  NFOLDS = 5, HAC = k)$results[["1"]]$se)

## -------------------------------------------------------------------------------------
tst1 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]),
                        obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 1,
                        dictionaries = dict_inf, NFOLDS = 5, B = 500)
q1 <- tst1$results[["1"]]
c(statistic = q1$statistic, p_value = q1$p_value)
q1$critical_values

## -------------------------------------------------------------------------------------
tst2 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]),
                        obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 2,
                        dictionaries = dict_inf, NFOLDS = 5, B = 500)
q2 <- tst2$results[["2"]]
c(statistic = q2$statistic, p_value = q2$p_value)
q2$critical_values

## ----profiles, fig.height = 3.2-------------------------------------------------------
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.4, 1))
for (z in list(list(q = q1, main = "covariate 1  (no effect)"),
               list(q = q2, main = "covariate 2  (effect)"))) {
  pr <- z$q$profile
  plot(pr$w, abs(pr$Psi_w), type = "b", pch = 16, main = z$main,
       xlab = "location w", ylab = expression(group("|", Psi[list(w, h)], "|")),
       ylim = range(0, abs(pr$Psi_w), z$q$critical_values))
  abline(h = z$q$critical_values[2], col = "red", lty = 2)
}
par(op)

## -------------------------------------------------------------------------------------
q2$profile

## -------------------------------------------------------------------------------------
rbind(covariate_1 = c(min = min(q1$profile$n_eff), max = max(q1$profile$n_eff)),
      covariate_2 = c(min = min(q2$profile$n_eff), max = max(q2$profile$n_eff)))
c(nT = obs_N * obs_T)

## ----cleanup, include=FALSE---------------------------------------------------
options(old)