Nothing
R/densitytoj.R
In panelhetero: Panel Data Analysis with Heterogeneous Dynamics
Defines functions
tojkd
Documented in
tojkd
#' The TOJ bias-corrected kernel density estimation
#'
#' The `tojkd()` function enables to implement the TOJ bias-corrected kernel
#' density estimation for the heterogeneous mean, the autocovariance,
#' and the autocorrelation.
#' The method is developed by Okui and Yanagi (2020).
#' For more details, see the package vignette with `vignette("panelhetero")`.
#'
#' @param data A matrix of panel data.
#' Each row corresponds to individual time series.
#' @param acov_order A non-negative integer of the order of autocovariance.
#' Default is 0.
#' @param acor_order A positive integer of the order of autocorrelation.
#' Default is 1.
#' @param mean_bw A scalar of bandwidth used for the estimation of
#' the denisty of mean.
#' Default is NULL, and the plug-in bandwidth is used.
#' @param acov_bw A scalar of bandwidth used for the estimation of
#' the denisty of autocovariance.
#' Default is NULL, and the plug-in bandwidth is used.
#' @param acor_bw A scalar of bandwidth used for the estimation of
#' the denisty of autocorrelation.
#' Default is NULL, and the plug-in bandwidth is used.
#'
#' @returns A list that contains the following elements:
#' \item{mean}{A plot of the corresponding density}
#' \item{acov}{A plot of the corresponding density}
#' \item{acor}{A plot of the corresponding density}
#' \item{mean_func}{A function that returns the corresponding density}
#' \item{acov_func}{A function that returns the corresponding density}
#' \item{acor_func}{A function that returns the corresponding density}
#' \item{bandwidth}{A Vector of the bandwidths}
#' \item{quantity}{A matrix of the estimated heterogeneous quantities}
#' \item{acov_order}{The order of autocovariance}
#' \item{acor_order}{The order of autocorrelation}
#' \item{N}{The number of cross-sectional units}
#' \item{S}{The length of time series}
#'
#' @examples
#' data <- panelhetero::simulation(N = 300, S = 50)
#' panelhetero::tojkd(data = data)
#'
#' @references Okui, R. and Yanagi, T., 2020.
#' Kernel estimation for panel data with heterogeneous dynamics.
#' The Econometrics Journal, 23(1), pp.156-175.
#'
#' @export
#'
tojkd <- function(data,
acov_order = 0,
acor_order = 1,
mean_bw = NULL,
acov_bw = NULL,
acor_bw = NULL) {
# Error handling -------------------------------------------------------------
error3(data = data,
acov_order = acov_order,
acor_order = acor_order,
mean_bw = mean_bw,
acov_bw = acov_bw,
acor_bw = acor_bw)
# Variable definitions -------------------------------------------------------
# Initialization
x <- NULL
# Omit NA
data <- stats::na.omit(data)
# Sample size
N <- nrow(data)
S <- ncol(data)
# Estimated means, autocovariances, autocorrelations
mean_est <- rowMeans(data)
acov_est <- apply(X = data,
MARGIN = 1,
FUN = acov,
acov_order = acov_order)
acor_est <- apply(X = data,
MARGIN = 1,
FUN = acor,
acor_order = acor_order)
# Plug-in bandwidth
if (is.null(mean_bw)) {
mean_bw <- KernSmooth::dpik(x = mean_est,
scalest = "minim",
kernel = "normal")
}
if (is.null(acov_bw)) {
acov_bw <- KernSmooth::dpik(x = acov_est,
scalest = "minim",
kernel = "normal")
}
if (is.null(acor_bw)) {
acor_bw <- KernSmooth::dpik(x = acor_est,
scalest = "minim",
kernel = "normal")
}
# Limits used for ggplot2
mean_lim <- c(min(mean_est),
max(mean_est))
acov_lim <- c(min(acov_est),
max(acov_est))
acor_lim <- c(min(acor_est),
max(acor_est))
# TOJ bias-correction
if (S %% 6 == 0) {
# Split panel data for T equivalent to 0 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:(S / 3)]
data32 <- data[, (S / 3 + 1):(2*S / 3)]
data33 <- data[, (2 * S / 3 + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest0(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest0(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest0(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)
}
} else if (S %% 6 == 1) {
# Split panel data for T equivalent to 1 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3))]
data33 <- data[, (2 * floor(S / 3) + 1):S]
data34 <- data[, 1:floor(S / 3)]
data35 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * floor(S / 3) + 2):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data39 <- data[, (2 * floor(S / 3) + 2):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest1(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest1(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest1(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else if (S %% 6 == 2) {
# Split panel data for T equivalent to 2 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1) ]
data33 <- data[, (2 * ceiling(S / 3)):S]
data34 <- data[, 1:ceiling(S / 3)]
data35 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * ceiling(S / 3)):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * ceiling(S / 3))]
data39 <- data[, (2 * ceiling(S / 3) + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest2(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest2(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest2(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else if (S %% 6 == 3) {
# Split panel data for T equivalent to 3 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:(S / 3)]
data32 <- data[, (S / 3 + 1):(2*S / 3)]
data33 <- data[, (2 * S / 3 + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest3(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest3(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest3(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)
}
} else if (S %% 6 == 4) {
# Split panel data for T equivalent to 4 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3))]
data33 <- data[, (2 * floor(S / 3) + 1):S]
data34 <- data[, 1:floor(S / 3)]
data35 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * floor(S / 3) + 2):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data39 <- data[, (2 * floor(S / 3) + 2):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest4(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest4(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest4(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else {
# Split panel data for T equivalent to 5 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1) ]
data33 <- data[, (2 * ceiling(S / 3)):S]
data34 <- data[, 1:ceiling(S / 3)]
data35 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * ceiling(S / 3)):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * ceiling(S / 3))]
data39 <- data[, (2 * ceiling(S / 3) + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures by ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest5(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest5(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest5(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
}
# Results
bandwidth <- c(mean_bw,
acov_bw,
acor_bw)
quantity <- cbind(mean_est,
acov_est,
acor_est)
names(bandwidth) <- colnames(quantity) <-
c("mean", "autocovariance", "autocorrelation")
return(list(mean = mean_plot,
acov = acov_plot,
acor = acor_plot,
mean_func = mean_func,
acov_func = acov_func,
acor_func = acor_func,
bandwidth = bandwidth,
quantity = quantity,
acov_order = acov_order,
acor_order = acor_order,
N = N,
S = S)
)
}
#' Compute TOJ kernel density estimates for T equivalent to 0 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest0 <- Vectorize(FUN = function(x, X, X21, X22, X31, X32, X33, h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33) / 3
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 1 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest1 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 2 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest2 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 3 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest3 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
# TOJ estimate
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33) / 3
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 4 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest4 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 5 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest5 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
Try the panelhetero package in your browser
Any scripts or data that you put into this service are public.
panelhetero documentation built on May 31, 2023, 5:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions tojkd
Documented in tojkd
#' The TOJ bias-corrected kernel density estimation
#'
#' The `tojkd()` function enables to implement the TOJ bias-corrected kernel
#' density estimation for the heterogeneous mean, the autocovariance,
#' and the autocorrelation.
#' The method is developed by Okui and Yanagi (2020).
#' For more details, see the package vignette with `vignette("panelhetero")`.
#'
#' @param data A matrix of panel data.
#' Each row corresponds to individual time series.
#' @param acov_order A non-negative integer of the order of autocovariance.
#' Default is 0.
#' @param acor_order A positive integer of the order of autocorrelation.
#' Default is 1.
#' @param mean_bw A scalar of bandwidth used for the estimation of
#' the denisty of mean.
#' Default is NULL, and the plug-in bandwidth is used.
#' @param acov_bw A scalar of bandwidth used for the estimation of
#' the denisty of autocovariance.
#' Default is NULL, and the plug-in bandwidth is used.
#' @param acor_bw A scalar of bandwidth used for the estimation of
#' the denisty of autocorrelation.
#' Default is NULL, and the plug-in bandwidth is used.
#'
#' @returns A list that contains the following elements:
#' \item{mean}{A plot of the corresponding density}
#' \item{acov}{A plot of the corresponding density}
#' \item{acor}{A plot of the corresponding density}
#' \item{mean_func}{A function that returns the corresponding density}
#' \item{acov_func}{A function that returns the corresponding density}
#' \item{acor_func}{A function that returns the corresponding density}
#' \item{bandwidth}{A Vector of the bandwidths}
#' \item{quantity}{A matrix of the estimated heterogeneous quantities}
#' \item{acov_order}{The order of autocovariance}
#' \item{acor_order}{The order of autocorrelation}
#' \item{N}{The number of cross-sectional units}
#' \item{S}{The length of time series}
#'
#' @examples
#' data <- panelhetero::simulation(N = 300, S = 50)
#' panelhetero::tojkd(data = data)
#'
#' @references Okui, R. and Yanagi, T., 2020.
#' Kernel estimation for panel data with heterogeneous dynamics.
#' The Econometrics Journal, 23(1), pp.156-175.
#'
#' @export
#'
tojkd <- function(data,
acov_order = 0,
acor_order = 1,
mean_bw = NULL,
acov_bw = NULL,
acor_bw = NULL) {
# Error handling -------------------------------------------------------------
error3(data = data,
acov_order = acov_order,
acor_order = acor_order,
mean_bw = mean_bw,
acov_bw = acov_bw,
acor_bw = acor_bw)
# Variable definitions -------------------------------------------------------
# Initialization
x <- NULL
# Omit NA
data <- stats::na.omit(data)
# Sample size
N <- nrow(data)
S <- ncol(data)
# Estimated means, autocovariances, autocorrelations
mean_est <- rowMeans(data)
acov_est <- apply(X = data,
MARGIN = 1,
FUN = acov,
acov_order = acov_order)
acor_est <- apply(X = data,
MARGIN = 1,
FUN = acor,
acor_order = acor_order)
# Plug-in bandwidth
if (is.null(mean_bw)) {
mean_bw <- KernSmooth::dpik(x = mean_est,
scalest = "minim",
kernel = "normal")
}
if (is.null(acov_bw)) {
acov_bw <- KernSmooth::dpik(x = acov_est,
scalest = "minim",
kernel = "normal")
}
if (is.null(acor_bw)) {
acor_bw <- KernSmooth::dpik(x = acor_est,
scalest = "minim",
kernel = "normal")
}
# Limits used for ggplot2
mean_lim <- c(min(mean_est),
max(mean_est))
acov_lim <- c(min(acov_est),
max(acov_est))
acor_lim <- c(min(acor_est),
max(acor_est))
# TOJ bias-correction
if (S %% 6 == 0) {
# Split panel data for T equivalent to 0 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:(S / 3)]
data32 <- data[, (S / 3 + 1):(2*S / 3)]
data33 <- data[, (2 * S / 3 + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest0,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest0(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest0(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest0(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)
}
} else if (S %% 6 == 1) {
# Split panel data for T equivalent to 1 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3))]
data33 <- data[, (2 * floor(S / 3) + 1):S]
data34 <- data[, 1:floor(S / 3)]
data35 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * floor(S / 3) + 2):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data39 <- data[, (2 * floor(S / 3) + 2):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest1,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest1(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest1(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest1(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else if (S %% 6 == 2) {
# Split panel data for T equivalent to 2 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1) ]
data33 <- data[, (2 * ceiling(S / 3)):S]
data34 <- data[, 1:ceiling(S / 3)]
data35 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * ceiling(S / 3)):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * ceiling(S / 3))]
data39 <- data[, (2 * ceiling(S / 3) + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest2,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest2(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest2(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest2(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else if (S %% 6 == 3) {
# Split panel data for T equivalent to 3 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:(S / 3)]
data32 <- data[, (S / 3 + 1):(2*S / 3)]
data33 <- data[, (2 * S / 3 + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest3,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest3(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest3(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest3(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
h = acor_bw)
}
} else if (S %% 6 == 4) {
# Split panel data for T equivalent to 4 modulo 6
data21 <- data[, 1:(S / 2)]
data22 <- data[, (S / 2 + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3))]
data33 <- data[, (2 * floor(S / 3) + 1):S]
data34 <- data[, 1:floor(S / 3)]
data35 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * floor(S / 3) + 2):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data39 <- data[, (2 * floor(S / 3) + 2):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures using ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest4,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest4(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest4(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest4(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
} else {
# Split panel data for T equivalent to 5 modulo 6
data21 <- data[, 1:floor(S / 2)]
data22 <- data[, (floor(S / 2) + 1):S]
data23 <- data[, 1:ceiling(S / 2)]
data24 <- data[, (ceiling(S / 2) + 1):S]
data31 <- data[, 1:floor(S / 3)]
data32 <- data[, (floor(S / 3) + 1):(2 * floor(S / 3) + 1) ]
data33 <- data[, (2 * ceiling(S / 3)):S]
data34 <- data[, 1:ceiling(S / 3)]
data35 <- data[, (ceiling(S / 3) + 1):(2 * floor(S / 3) + 1)]
data36 <- data[, (2 * ceiling(S / 3)):S]
data37 <- data[, 1:ceiling(S / 3)]
data38 <- data[, (ceiling(S / 3) + 1):(2 * ceiling(S / 3))]
data39 <- data[, (2 * ceiling(S / 3) + 1):S]
# Estimated quantities for split panel data
mean_est21 <- rowMeans(data21)
mean_est22 <- rowMeans(data22)
mean_est23 <- rowMeans(data23)
mean_est24 <- rowMeans(data24)
mean_est31 <- rowMeans(data31)
mean_est32 <- rowMeans(data32)
mean_est33 <- rowMeans(data33)
mean_est34 <- rowMeans(data34)
mean_est35 <- rowMeans(data35)
mean_est36 <- rowMeans(data36)
mean_est37 <- rowMeans(data37)
mean_est38 <- rowMeans(data38)
mean_est39 <- rowMeans(data39)
acov_est21 <- apply(data21, MARGIN = 1, acov, acov_order = acov_order)
acov_est22 <- apply(data22, MARGIN = 1, acov, acov_order = acov_order)
acov_est23 <- apply(data23, MARGIN = 1, acov, acov_order = acov_order)
acov_est24 <- apply(data24, MARGIN = 1, acov, acov_order = acov_order)
acov_est31 <- apply(data31, MARGIN = 1, acov, acov_order = acov_order)
acov_est32 <- apply(data32, MARGIN = 1, acov, acov_order = acov_order)
acov_est33 <- apply(data33, MARGIN = 1, acov, acov_order = acov_order)
acov_est34 <- apply(data34, MARGIN = 1, acov, acov_order = acov_order)
acov_est35 <- apply(data35, MARGIN = 1, acov, acov_order = acov_order)
acov_est36 <- apply(data36, MARGIN = 1, acov, acov_order = acov_order)
acov_est37 <- apply(data37, MARGIN = 1, acov, acov_order = acov_order)
acov_est38 <- apply(data38, MARGIN = 1, acov, acov_order = acov_order)
acov_est39 <- apply(data39, MARGIN = 1, acov, acov_order = acov_order)
acor_est21 <- apply(data21, MARGIN = 1, acor, acor_order = acor_order)
acor_est22 <- apply(data22, MARGIN = 1, acor, acor_order = acor_order)
acor_est23 <- apply(data23, MARGIN = 1, acor, acor_order = acor_order)
acor_est24 <- apply(data24, MARGIN = 1, acor, acor_order = acor_order)
acor_est31 <- apply(data31, MARGIN = 1, acor, acor_order = acor_order)
acor_est32 <- apply(data32, MARGIN = 1, acor, acor_order = acor_order)
acor_est33 <- apply(data33, MARGIN = 1, acor, acor_order = acor_order)
acor_est34 <- apply(data34, MARGIN = 1, acor, acor_order = acor_order)
acor_est35 <- apply(data35, MARGIN = 1, acor, acor_order = acor_order)
acor_est36 <- apply(data36, MARGIN = 1, acor, acor_order = acor_order)
acor_est37 <- apply(data37, MARGIN = 1, acor, acor_order = acor_order)
acor_est38 <- apply(data38, MARGIN = 1, acor, acor_order = acor_order)
acor_est39 <- apply(data39, MARGIN = 1, acor, acor_order = acor_order)
# Make figures by ggplot2
mean_plot <- ggplot2::ggplot(data = data.frame(x = mean_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous mean") +
ggplot2::theme_bw()
acov_plot <- ggplot2::ggplot(data = data.frame(x = acov_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocovariance") +
ggplot2::theme_bw()
acor_plot <- ggplot2::ggplot(data = data.frame(x = acor_lim),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = tojkdest5,
args = list(X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)) +
ggplot2::labs(x = "x", y = "") +
ggplot2::ggtitle("The heterogeneous autocorrelation") +
ggplot2::theme_bw()
# Functions
mean_func <- function(x) {
tojkdest5(x = x,
X = mean_est,
X21 = mean_est21,
X22 = mean_est22,
X23 = mean_est23,
X24 = mean_est24,
X31 = mean_est31,
X32 = mean_est32,
X33 = mean_est33,
X34 = mean_est34,
X35 = mean_est35,
X36 = mean_est36,
X37 = mean_est37,
X38 = mean_est38,
X39 = mean_est39,
h = mean_bw)
}
acov_func <- function(x) {
tojkdest5(x = x,
X = acov_est,
X21 = acov_est21,
X22 = acov_est22,
X23 = acov_est23,
X24 = acov_est24,
X31 = acov_est31,
X32 = acov_est32,
X33 = acov_est33,
X34 = acov_est34,
X35 = acov_est35,
X36 = acov_est36,
X37 = acov_est37,
X38 = acov_est38,
X39 = acov_est39,
h = acov_bw)
}
acor_func <- function(x) {
tojkdest5(x = x,
X = acor_est,
X21 = acor_est21,
X22 = acor_est22,
X23 = acor_est23,
X24 = acor_est24,
X31 = acor_est31,
X32 = acor_est32,
X33 = acor_est33,
X34 = acor_est34,
X35 = acor_est35,
X36 = acor_est36,
X37 = acor_est37,
X38 = acor_est38,
X39 = acor_est39,
h = acor_bw)
}
}
# Results
bandwidth <- c(mean_bw,
acov_bw,
acor_bw)
quantity <- cbind(mean_est,
acov_est,
acor_est)
names(bandwidth) <- colnames(quantity) <-
c("mean", "autocovariance", "autocorrelation")
return(list(mean = mean_plot,
acov = acov_plot,
acor = acor_plot,
mean_func = mean_func,
acov_func = acov_func,
acor_func = acor_func,
bandwidth = bandwidth,
quantity = quantity,
acov_order = acov_order,
acor_order = acor_order,
N = N,
S = S)
)
}
#' Compute TOJ kernel density estimates for T equivalent to 0 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest0 <- Vectorize(FUN = function(x, X, X21, X22, X31, X32, X33, h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33) / 3
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 1 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest1 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 2 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest2 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 3 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest3 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
# TOJ estimate
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33) / 3
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 4 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest4 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22) / 2 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
#' Compute TOJ kernel density estimate for T equivalent to 5 modulo 6
#'
#' @param x An evaluation point
#' @param X A vector of cross-sectional data
#' @param X21 A vector of half-panel cross-sectional data 21
#' @param X22 A vector of half-panel cross-sectional data 22
#' @param X23 A vector of half-panel cross-sectional data 23
#' @param X24 A vector of half-panel cross-sectional data 24
#' @param X31 A vector of third-panel cross-sectional data 31
#' @param X32 A vector of third-panel cross-sectional data 32
#' @param X33 A vector of third-panel cross-sectional data 33
#' @param X34 A vector of third-panel cross-sectional data 34
#' @param X35 A vector of third-panel cross-sectional data 35
#' @param X36 A vector of third-panel cross-sectional data 36
#' @param X37 A vector of third-panel cross-sectional data 37
#' @param X38 A vector of third-panel cross-sectional data 38
#' @param X39 A vector of third-panel cross-sectional data 39
#' @param h A scalar of bandwidth
#'
#' @returns A vector of kernel density estimates
#'
#' @noRd
#'
tojkdest5 <- Vectorize(FUN = function(x,
X,
X21,
X22,
X23,
X24,
X31,
X32,
X33,
X34,
X35,
X36,
X37,
X38,
X39,
h) {
# Sample size
N <- length(X)
# Estimates
est <- sum(stats::dnorm((x - X) / h)) / (N * h)
est21 <- sum(stats::dnorm((x - X21) / h)) / (N * h)
est22 <- sum(stats::dnorm((x - X22) / h)) / (N * h)
est23 <- sum(stats::dnorm((x - X23) / h)) / (N * h)
est24 <- sum(stats::dnorm((x - X24) / h)) / (N * h)
est31 <- sum(stats::dnorm((x - X31) / h)) / (N * h)
est32 <- sum(stats::dnorm((x - X32) / h)) / (N * h)
est33 <- sum(stats::dnorm((x - X33) / h)) / (N * h)
est34 <- sum(stats::dnorm((x - X34) / h)) / (N * h)
est35 <- sum(stats::dnorm((x - X35) / h)) / (N * h)
est36 <- sum(stats::dnorm((x - X36) / h)) / (N * h)
est37 <- sum(stats::dnorm((x - X37) / h)) / (N * h)
est38 <- sum(stats::dnorm((x - X38) / h)) / (N * h)
est39 <- sum(stats::dnorm((x - X39) / h)) / (N * h)
# TOJ estimates
tojest <- 3.536 * est -
4.072 * (est21 + est22 + est23 + est24) / 4 +
1.536 * (est31 + est32 + est33 + est34 +
est35 + est36 + est37 + est38 + est39) / 9
# Ensure non-negative estimates
tojest <- ifelse(tojest >= 0, tojest, 0)
return(tojest)
}, vectorize.args = "x")
Try the panelhetero package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.