Nothing
#' Hahn Kuehrsteiner Estimator for PVAR Model
#'
#' @param dependent_vars Dependent variables
#' @param exog_vars Exogenous variables
#' @param transformation Demeaning \code{"demean"}
#' @param data Data set
#' @param panel_identifier Vector of panel identifiers
#'
#' @export
#' @description
#' This function estimates a stationary PVAR with fixed effects.
#' @examples
#' data(Dahlberg)
#' ex1_hk <-
#' pvarhk(dependent_vars = c("expenditures", "revenues", "grants"),
#' transformation = "demean",
#' data = Dahlberg,
#' panel_identifier= c("id", "year"))
#'
#' summary(ex1_hk)
#'
#' @references
#' Hahn J., Kuehrsteiner G. (2002) Asymptotically Unbiased Inference for a Dynamic Panel Model with Fixed Effects When Both n and T Are Large, \emph{Econometrica}, \bold{70}(4), 1639--1657
pvarhk <- function(dependent_vars,
# lags = 1,
exog_vars,
transformation = c("demean"),
data,
panel_identifier = c(1, 2)){
# works only with lag 1
lags <- 1
# Set up Set_Vars ---------------------------------------------------------
# Drop factor levels from data
data <- droplevels(data)
# Select only required variables from data
required_vars <- c(dependent_vars)
if(!missing(exog_vars)) required_vars <- c(required_vars, exog_vars)
if (is.numeric(panel_identifier)) {
Set_Vars <-
data[, c(colnames(data)[panel_identifier], required_vars)]
} else {
Set_Vars <-
data[, c(panel_identifier, required_vars)]
}
nof_dependent_vars <- length(dependent_vars)
if(missing(exog_vars)) {
nof_exog_vars <- 0
} else {
nof_exog_vars <- length(exog_vars)
}
# sort categories and periods
categories <- sort(unique(Set_Vars[, 1]))
periods <- sort(unique(Set_Vars[, 2]))
nof_categories <- length(categories)
nof_periods <- length(periods)
# keep names of panel identifiers
name_category <- names(Set_Vars)[1]
name_period <- names(Set_Vars)[2]
names(Set_Vars)[1:2] <- c("category", "period")
Set_Vars <- Set_Vars[order(Set_Vars$category, Set_Vars$period,decreasing = FALSE),]
Set_Vars$category <- factor(Set_Vars$category)
Set_Vars$period <- factor(Set_Vars$period)
# Descriptive statistics --------------------------------------------------
nof_observations <- dim(data)[1]
data_panel_identifier <- data[panel_identifier]
obs_per_group_avg <- mean(table(data_panel_identifier[,1]), na.rm = TRUE)
obs_per_group_min <- min(table(data_panel_identifier[,1]), na.rm = TRUE)
obs_per_group_max <- max(table(data_panel_identifier[,1]), na.rm = TRUE)
nof_groups <- length(unique(data_panel_identifier[,1]))
# First demean and then lag the data --------------------------------------
# Add lags of dependent_vars
Set_Vars <- cbind(Set_Vars,
do.call(
rbind,
mapply(
function(i)
panel_lag(Set_Vars[Set_Vars$category == categories[i], c("period", dependent_vars)], lags),
1:length(categories),
SIMPLIFY = FALSE
)
))
# Add forwardlags of exogenous vars to do the trick suggested by Hahn Kuersteiner: page 1640
# Eq (2) by adding an additional equation for the exogenous variable, but restriciting its
# coefficients.
if (!missing(exog_vars)) {
Set_Vars <- cbind(Set_Vars,
do.call(
rbind,
mapply(
function(i)
panel_forward_lag(Set_Vars[Set_Vars$category == categories[i],
c("period", exog_vars)], lags),
1:length(categories),
SIMPLIFY = FALSE
)
))
}
# Add demeaned variables, but also from the lagged dependent variables above
Set_Vars <- cbind(Set_Vars,
do.call(
rbind,
mapply(
function(i)
panel_demean(Set_Vars[Set_Vars$category == categories[i],], demean),
1:length(categories),
SIMPLIFY = FALSE
)
))
# Handle exogenous variables (if necessary) -------------------------------
# If no exogenous variables then proceed in a standard way.
if (missing(exog_vars)) {
dependent.vars <- paste0("demeaned_", dependent_vars)
# Add all chosen lags to the lagged variables.
lagged.vars <- c()
for (l1 in 1:lags){
lagged.vars <- c(lagged.vars, paste0("demeaned_", "lag", l1,"_", dependent_vars))
}
# Generate Matrix R in order to estimate the restricted OLS-estimator
l <- length(dependent.vars) - 1
k <- lags
R <- matrix(0, nrow = l * k, ncol = (l + k) ^ 2)
R[1:(l * k), ((l + k) * l + 1):(l ^ 2 + 2 * l * k)] <-
diag(1, nrow = l * k)
}
# If exogenous variables then set up the helper additional equation with parameter restrictions.
if (!missing(exog_vars)) {
dependent.vars <- c(paste0("demeaned_", dependent_vars),
paste0("demeaned_forwardlag1_", exog_vars))
lagged.vars <- c(paste0("demeaned_", "lag1_", dependent_vars),
paste0("demeaned_", exog_vars))
# Generate Matrix R in order to estimate the restricted OLS-estimator
l <- length(dependent.vars) - length(exog_vars)
k <- nof_exog_vars
# R is the restriction matrix:
# ncol of R is the number of parameters in the model
# nrow is the number of restrictions. e.g. 3 endo + 1 exo. ==>
# 16 parameters in the extended model, 3 equations for the endos,
# 1 equation for the exogenous variables.
if (k == 1){
R <- matrix(0, nrow = l * k, ncol = (l + k) ^ 2)
R[1:(l * k), ((l + k) * l + 1):(l ^ 2 + 2 * l * k)] <-
diag(1, nrow = l * k)
}
# Only restrict the lagged endo parameters in the exo-equations to zero.
if (k > 1){
R <- matrix(0, nrow = l * k, ncol = (l + k) ^ 2)
for (i0 in 1:nof_exog_vars){
R[(i0*l-l+1):(l * i0), ((l + nof_exog_vars) * l + (1+(i0-1)*(l+nof_exog_vars) ) ):
(((l + nof_exog_vars) * l + (1+(i0-1)*(l+nof_exog_vars) ) )+l-1)] <-
diag(1, nrow = l)
}
}
}
Set_Vars <- na.exclude(Set_Vars)
n <-length(unique(Set_Vars$category))
catgnames<-unique(Set_Vars$category)
Matrix.Greek.T<- as.matrix(Set_Vars[,lagged.vars])
List.Te <- matrix(0, nrow = length(catgnames), ncol = 1)
# sth for unbalanced panels...to calculate an average T.
for (i1 in 1:length(catgnames)){
# Schreibe Vektor welcher f?r jede ILZ-Nummer die Zeitpunkte, an denen Daten erfasst wurden, ausgibt
Number.of.periods <- Set_Vars$Date_index[Set_Vars$category==catgnames[i1]]
Te <- sum(Set_Vars$category==catgnames[i1])
Matrix.Greek.T[Set_Vars$category == catgnames[i1],]<- 1/Te *
Matrix.Greek.T[Set_Vars$category == catgnames[i1],]
List.Te[i1,1] <- Te
}
############################################################################
# Build OLS Sum1, page 1641 Hahn-Kuersteiner 2002 --------------------------
M1 <- t(as.matrix(Matrix.Greek.T)) %*% as.matrix(Set_Vars[,lagged.vars])
M2 <- t(as.matrix(Set_Vars[,lagged.vars]))%*% as.matrix(Set_Vars[,lagged.vars])
M3 <- t(as.matrix(Set_Vars[, lagged.vars])) %*% as.matrix(Set_Vars[,dependent.vars])
M4 <- t(as.matrix(Matrix.Greek.T))%*% as.matrix(Set_Vars[,dependent.vars])
M5 <- t(as.matrix(Matrix.Greek.T))%*% as.matrix(Set_Vars[,lagged.vars])
Greek.T <- 1/n*M1
# Calculate OLS ------------------------------------------------------------
sum1.theta.hat <-M2
sum2.theta.hat<-M3
# here the equation are in the columns:
theta.hat.t <- solve(sum1.theta.hat)%*%as.matrix(sum2.theta.hat)
# Calculate restricted OLS (in case of exogenous variables) ----------------
if (!missing(exog_vars)) {
X <- diag(1, l + k) %x% as.matrix(Set_Vars[, lagged.vars])
XX_inv <- diag(1, l + nof_exog_vars) %x% solve(M2)
vec.theta.hat.R <-
matrix(theta.hat.t, ncol = 1) + XX_inv %*% t(R) %*% solve(R %*% XX_inv %*% t(R)) %*% (-R %*% matrix(theta.hat.t, ncol = 1))
theta.hat.R <- t(matrix(vec.theta.hat.R, ncol = l + k))
colnames(theta.hat.R) <- colnames(t(theta.hat.t))
rownames(theta.hat.R) <- rownames(t(theta.hat.t))
}
if (missing(exog_vars)){
theta.hat.R <- t(theta.hat.t)
}
# Standard Errors for OLS/restricted OLS estimation ------------------------
ols.res <- as.matrix(Set_Vars[, dependent.vars]) - as.matrix(Set_Vars[, lagged.vars]) %*% t(theta.hat.R)
ols.res.squ <- ( t(ols.res) %*% ols.res ) / ( dim(ols.res)[1] - nrow(theta.hat.t))
ols_variance <- ols.res.squ %x% solve(sum1.theta.hat)
ols_variance_sde <- sqrt(matrix(diag(ols_variance), nrow=length(lagged.vars),byrow=F))
z_values_ols <- t(theta.hat.R)/ols_variance_sde
p_values_ols <- 2*pnorm(-abs(z_values_ols))
# Calculate Theta-Double-Hat -----------------------------------------------
sum2.vec.theta.double.hat <- 1/n * M4
omega.hat <- Greek.T - theta.hat.R%*%Greek.T%*%t(theta.hat.R)
# Coefficients of an equation are in the rows of the matrix.
theta.double.hat <- t( solve(Greek.T) %*% (sum2.vec.theta.double.hat +
1/mean(List.Te[,1])*solve(diag(l+k)-theta.hat.R)%*%omega.hat) )
# Standard Error per Equation-Block for Hahn Kuehrsteiner ------------------
hk.res <- as.matrix(Set_Vars[, dependent.vars]) - as.matrix(Set_Vars[, lagged.vars]) %*% t(theta.double.hat)
hk.res.squ <- ( t(hk.res) %*% hk.res ) / ( dim(hk.res)[1] - nrow(theta.double.hat))
hk_variance <- hk.res.squ %x% solve(sum1.theta.hat)
hk_variance_sde <- sqrt(matrix(diag(hk_variance), nrow=length(lagged.vars),byrow=F))
z_values_hk <- theta.double.hat/hk_variance_sde
p_values_hk <- 2*pnorm(-abs(z_values_hk))
# Results ------------------------------------------------------------------
OLS.results <- list(coef = t(theta.hat.R), se = ols_variance_sde, pvalues = p_values_ols)
HK.results <- list(coef = theta.double.hat, se = hk_variance_sde, pvalues = p_values_hk)
inputargs <- list(dependent_vars = dependent_vars,
lags = lags,
transformation = transformation,
Set_Vars = Set_Vars,
panel_identifier = panel_identifier,
nof_observations = nof_observations,
obs_per_group_avg = obs_per_group_avg,
obs_per_group_min = obs_per_group_min,
obs_per_group_max = obs_per_group_max,
nof_groups = nof_groups
)
results <- append(list(OLS = OLS.results, HK = HK.results), c(inputargs))
class(results) <- "pvarhk"
results
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.