#' Suboptimal multi-step System-GMM estimator for AR(1) panel data model
#'
#' @author Michal Oleszak
#'
#' @description Computes an enhanced version of the Blundell-Bond (System-GMM)
#' estimator for panel data by means of replacing the standard
#' GMM-weighting matrix by its sub-optimal version, thus increasing
#' estimator's efficiency.
#'
#' @references Youssef, A. and Abonazel, M. (2015). Alternative GMM estimators
#' for first-order autoregressive panel model: An improving
#' efficiency approach. MPRA Paper No. 68674; Forthcoming in:
#' Communications in Statistics - Simulation and Computation,
#' \url{https://mpra.ub.uni-muenchen.de/68674/1/MPRA_paper_68674.pdf}
#'
#' @param Y matrix of size (T x N) with the dependent variable
#' @param model one of: onestep, twosteps, threesteps; more steps should
#' increase efficiency, but might be computationally infeasible
#' (a singular matrix needs to be inverted); if this is the case,
#' generalised inverse is used
#'
#' @return The estimated value of the auto-regressive parameter.
#'
#' @importFrom plm pgmm
#' @importFrom Matrix bdiag solve
#' @importFrom pracma ones
#' @importFrom MASS ginv
#' @importFrom utils head tail
#' @export
ssys_gmm <- function(Y, model = c("onestep", "twosteps", "threesteps")) {
N <- ncol(Y)
T <- nrow(Y)
# tranform data to a format accepted by pgmm()
unit <- rep(1:N, each = T)
time <- rep(1:T, N)
yplm <- c(Y)
plmdata <- cbind(unit, time, yplm)
# obtain variance ratio 'r'
DIF <- pgmm(yplm ~ lag(yplm, 1) | lag(yplm, 2:99),
data = as.data.frame(plmdata), effect = "individual",
method = "onestep", transformation = "d"
)
SYS <- pgmm(yplm ~ lag(yplm, 1) | lag(yplm, 2:99),
data = as.data.frame(plmdata), effect = "individual",
method = "onestep", transformation = "ld"
)
delta_u_i_hat <- DIF$residuals
sysres <- SYS$residuals
delta_u_i_tilde <- lapply(sysres, function(x) head(x, (T - 2)))
u_i_tilde <- lapply(sysres, function(x) tail(x, (T - 1)))
sigma2eps <- (sum(unlist(lapply(delta_u_i_hat, function(x) t(x) %*% x)))) /
(2 * N * (T - 2))
sigma2gamma <- (sum(unlist(lapply(u_i_tilde, function(x) t(x) %*% x)) -
(unlist(lapply(delta_u_i_tilde, function(x) t(x) %*% x)) / 2))) /
(N * (T - 2))
r <- sigma2gamma / sigma2eps
# construct moment condition matrices
L_i <- list()
D_i <- list()
D_i_list <- rep(list(vector("list", T - 2)), N)
S_i <- list()
Y_diff <- matrix(NA, ncol = N, nrow = T - 2)
for (t in 2:(T - 1)) {
Y_diff[t - 1, ] <- Y[t, ] - Y[t - 1, ]
}
for (i in 1:N) {
L_i[[i]] <- diag(Y_diff[, i])
for (t in 1:(T - 2)) {
D_i_list[[i]][[t]] <- matrix(c(Y[1:t, i]), nrow = 1)
}
D_i[[i]] <- bdiag(D_i_list[[i]])
S_i[[i]] <- bdiag(D_i[[i]], L_i[[i]])
}
F <- matrix(0, nrow = T - 2, ncol = T - 1)
for (i in 1:nrow(F)) {
for (j in 1:ncol(F)) {
if (i == j) {
F[i, j] <- -1
}
if (j == i + 1) {
F[i, j] <- 1
}
}
}
A <- diag(T - 2) + r * ones(T - 2)
G <- bdiag(F %*% t(F), A)
Z1 <- matrix(0, nrow = ncol(S_i[[1]]), ncol = ncol(S_i[[1]]))
for (i in 1:N) {
Z1 <- Z1 + t(as.matrix(S_i[[i]])) %*% G %*% S_i[[i]]
}
Z1 <- solve(Z1)
S <- S_i[[1]]
for (i in 2:N) {
S <- rbind(S, S_i[[i]])
}
Y_diff_temp <- matrix(NA, ncol = N, nrow = T - 2)
for (t in 3:(T)) {
Y_diff_temp[t - 2, ] <- Y[t, ] - Y[t - 1, ]
}
y_stacked <- c(c(Y_diff_temp), c(Y[-c(1:2), ]))
y_stacked_lag <- c(c(Y_diff), c(Y[-c(1, T), ]))
rho1 <- solve(t(y_stacked_lag) %*% S %*% Z1 %*% t(as.matrix(S)) %*% y_stacked_lag) %*%
(t(y_stacked_lag) %*% S %*% Z1 %*% t(as.matrix(S)) %*% y_stacked)
if (model == "onestep") {
return(as.numeric(rho1))
}
# step 2
fitted1lev <- matrix(NA, nrow = T - 2, ncol = N)
fitted1diff <- matrix(NA, nrow = T - 2, ncol = N)
for (t in 1:(T - 2)) {
fitted1lev[t, ] <- Y[t + 1, ] * as.numeric(rho1)
fitted1diff[t, ] <- Y_diff_temp[t, ] * as.numeric(rho1)
}
res1lev <- Y[-(1:2), ] - fitted1lev
res1diff <- Y_diff_temp - fitted1diff
res1 <- rbind(res1diff, res1lev)
Z2 <- matrix(0, nrow = ncol(Z1), ncol = ncol(Z1))
for (i in 1:N) {
Z2 <- Z2 + t(as.matrix(S_i[[i]])) %*% (res1[, i] %*% t(res1[, i])) %*% S_i[[i]]
}
Z2try <- try(Z2 <- solve(Z2))
if ("try-error" %in% class(Z2try)) {
print("Second-step matrix is singular. Generalised inverse was used.")
Z2 <- as.matrix(Z2)
Z2 <- ginv(Z2)
}
rho2 <- solve(t(y_stacked_lag) %*% S %*% Z2 %*% t(as.matrix(S)) %*% y_stacked_lag) %*%
(t(y_stacked_lag) %*% S %*% Z2 %*% t(as.matrix(S)) %*% y_stacked)
if (model == "twosteps") {
return(as.numeric(rho2))
}
# step 3
fitted2lev <- matrix(NA, nrow = T - 2, ncol = N)
fitted2diff <- matrix(NA, nrow = T - 2, ncol = N)
for (t in 1:(T - 2)) {
fitted2lev[t, ] <- Y[t + 1, ] * as.numeric(rho2)
fitted2diff[t, ] <- Y_diff_temp[t, ] * as.numeric(rho2)
}
res2lev <- Y[-(1:2), ] - fitted2lev
res2diff <- Y_diff_temp - fitted2diff
res2 <- rbind(res2diff, res2lev)
Z3 <- matrix(0, nrow = ncol(Z1), ncol = ncol(Z1))
for (i in 1:N) {
Z3 <- Z3 + t(as.matrix(S_i[[i]])) %*% (res2[, i] %*% t(res2[, i])) %*% S_i[[i]]
}
Z3try <- try(Z3 <- solve(Z3))
if ("try-error" %in% class(Z3try)) {
print("Third-step matrix is singular. Generalised inverse was used.")
Z3 <- as.matrix(Z3)
Z3 <- ginv(Z3)
}
rho3 <- solve(t(y_stacked_lag) %*% S %*% Z3 %*% t(as.matrix(S)) %*% y_stacked_lag) %*%
(t(y_stacked_lag) %*% S %*% Z3 %*% t(as.matrix(S)) %*% y_stacked)
if (model == "threesteps") {
return(as.numeric(rho3))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.