Nothing
#' Extract PVARFEOLS(p) Model Coefficients
#' @param object object
#' @param ... further arguments
#' @export
coef.pvarfeols <- function(object, ...) {
object$OLS$coef
}
#' @rdname se
#' @export
se.pvarfeols <- function(object, ...) {
object$OLS$se
}
#' @rdname pvalue
#' @export
pvalue.pvarfeols <- function(object, ...) {
object$OLS$pvalues
}
#' @rdname extract
#' @export
extract.pvarfeols <- function(model, ...) {
co.m <- coef(model)
modelnames <- row.names(co.m)
se.m <- se(model)
pval.m <- pvalue(model)
equationList <- list()
for (eq in 1:length(model$dependent_vars)) {
tr <- texreg::createTexreg(
coef.names = colnames(co.m),
coef = as.vector(co.m[eq,]),
se = as.vector(se.m[eq,]),
pvalues = as.vector(pval.m[eq,]),
model.name = modelnames[eq]
)
equationList[[eq]] <- tr
}
names(equationList) <- modelnames
equationList
}
#' S3 Print Method for pvarfeols
#' @param x object
#' @param ... further arguments
#' @method print pvarfeols
#' @export
print.pvarfeols <- function(x, ...) {
res <- extract(x)
print(texreg::screenreg(res, custom.model.names = names(res), digits = 4))
}
#' Knit Print Method for pvarfeols
#' @param x object
#' @param ... further arguments
#' @export
knit_print.pvarfeols <- function(x, ...) {
res <- extract(x)
knitr::asis_output(texreg::htmlreg(res, custom.model.names = names(res), digits = 4))
}
#' S3 Summary Method for pvarfeols
#' @param object object
#' @param ... further arguments
#' @method summary pvarfeols
#' @export
summary.pvarfeols <- function(object, ...)
{
x <- object
sumry <- list()
sumry$model_name <- "Fixed Effects OLS Panel VAR estimation"
sumry$transformation <- x$transformation
sumry$groupvariable <- x$panel_identifier[1]
sumry$timevariable <- x$panel_identifier[2]
sumry$nof_observations <- x$nof_observations
sumry$obs_per_group_min <- x$obs_per_group_min
sumry$obs_per_group_avg <- x$obs_per_group_avg
sumry$obs_per_group_max <- x$obs_per_group_max
sumry$nof_groups = x$nof_groups
sumry$results <- extract(x)
class(sumry) <- "summary.pvarfeols"
sumry
}
#' S3 Print Method for summary.pvarfeols
#' @param x object
#' @param ... further arguments
#' @method print summary.pvarfeols
#' @export
print.summary.pvarfeols <- function(x, ...) {
cat("---------------------------------------------------\n")
cat(x$model_name, "\n")
cat("---------------------------------------------------\n")
cat("Transformation:", x$transformation,"\n")
cat("Group variable:", x$groupvariable,"\n")
cat("Time variable:",x$timevariable,"\n")
cat("Number of observations =", x$nof_observations,"\n")
cat("Number of groups =", x$nof_groups, "\n")
cat("Obs per group: min =",x$obs_per_group_min,"\n")
cat(" avg =",x$obs_per_group_avg,"\n")
cat(" max =",x$obs_per_group_max,"\n")
#cat("---------------------------------------------------\n")
print(texreg::screenreg(x$results, custom.model.names = names(x$results), digits = 4))
cat("\n")
}
#' Knit Print summary Method
#' @param x object
#' @param ... further arguments
#' @export
knit_print.summary.pvarfeols <- function(x, ...) {
res <- paste0(c("<p><b>",x$model_name,"</b></p>",
"<p>Transformation: <em>",x$transformation,"</em><br>",
"Group variable: <em>",x$groupvariable,"</em><br>",
"Time variable: <em>",x$timevariable,"</em><br>",
"Number of observations = <em>",x$nof_observations,"</em><br>",
"Number of groups = <em>",x$nof_groups,"</em><br>",
"Obs per group: min = <em>",x$obs_per_group_min,"</em><br>",
"Obs per group: avg = <em>",x$obs_per_group_avg,"</em><br>",
"Obs per group: max = <em>",x$obs_per_group_max,"</em><br>",
"</p>",
texreg::htmlreg(x$results, custom.model.names = names(x$results), digits = 4, caption = "", center = FALSE),
"</p>"
)
)
knitr::asis_output(res)
}
#' @export
vcov.pvarfeols <- function(object, ...){
residuals <- na.exclude( as.matrix(object$residuals) )
Sigma_Hat1 <- ( t(residuals) %*% residuals ) / (nrow(residuals) - ncol(object$OLS$coef))
return(Sigma_Hat1)
}
#' @export
oirf.pvarfeols <- function(model, n.ahead){
Phi = coef(model)
# P <- t(chol(covm))
MA_Phi <- ma_phi_representation(Phi = Phi,
ma_approx_steps = n.ahead,
lags = model$lags)
Sigma_Hat1 <- vcov(model)
P <- chol(Sigma_Hat1)
irf_output <- list()
MA_Phi_P <- list()
# Calculate 2.3.27 in Luetkepohl p. 58.
for (i0 in 1:n.ahead){
MA_Phi_P[[i0]] <- MA_Phi[[i0]] %*% t(P)
}
# Redistribute the MA_Phi_P into the correct list form.
for (i0 in 1:nrow(Phi)){
zwischen <- matrix(NA, nrow = length(MA_Phi_P), ncol = nrow(Phi))
for (i1 in 1:length(MA_Phi_P)){
zwischen[i1,] <- MA_Phi_P[[i1]][,i0]
}
irf_output[[i0]] <- zwischen
colnames(irf_output[[i0]]) <- model$dependent_vars
}
names(irf_output) <- model$dependent_vars
class(irf_output) <- "pvarfeolsoirf"
return(irf_output)
}
#' @export
plot.pvarfeolsoirf <- function(x, cibootstrap_results, ...) {
dependent_vars <- names(x)
steps <- nrow(x[[1]])
plotdata <- do.call(rbind,lapply(x, reshape2::melt))
names(plotdata) <- c("periods", "impulse", "value")
plotdata$impulse <- paste(rep(dependent_vars, each = steps*length(dependent_vars)), "on", plotdata$impulse)
p <-
ggplot2::ggplot(plotdata, ggplot2::aes_(x = ~periods, y = ~value)) +
ggplot2::geom_line(colour = "orangered4") + ggplot2::facet_wrap(~impulse) +
ggplot2::labs(title = "Orthogonalized impulse response function") +
ggplot2::xlab("steps") + ggplot2::ylab("") +
ggplot2::scale_x_continuous(breaks = pretty(plotdata$periods)) +
ggplot2::theme_bw()
if (!missing(cibootstrap_results)) {
ci_data_upper <-
data.frame(do.call(rbind,lapply(cibootstrap_results$Upper, reshape2::melt)), Bound = "Upper")
names(ci_data_upper)[1:3] <- c("periods", "impulse", "value")
row.names(ci_data_upper) <- NULL
ci_data_upper$impulse <- paste(rep(dependent_vars,
each = steps*length(dependent_vars)), "on", ci_data_upper$impulse)
ci_data_lower <-
data.frame(do.call(rbind,lapply(cibootstrap_results$Lower, reshape2::melt)), Bound = "Lower")
names(ci_data_lower)[1:3] <- c("periods", "impulse", "value")
row.names(ci_data_lower) <- NULL
ci_data_lower$impulse <- paste(rep(dependent_vars,
each = steps*length(dependent_vars)), "on", ci_data_lower$impulse)
ci_data <- reshape2::dcast(rbind(ci_data_upper, ci_data_lower), periods + impulse ~ Bound, value.var = "value")
plotdata <-
merge(plotdata, ci_data, all.x = TRUE, by = c("periods", "impulse"), sort = FALSE)
p <-
ggplot2::ggplot(plotdata, ggplot2::aes_(x = ~periods, y = ~value)) +
ggplot2::geom_line(colour = "orangered4") +
ggplot2::geom_ribbon(ggplot2::aes_(x = ~periods,ymin = ~Lower, ymax = ~Upper), alpha = 0.3,fill = "steelblue2") +
ggplot2::facet_wrap(~impulse) +
ggplot2::labs(title = "Orthogonalized impulse response function") +
ggplot2::xlab("steps") + ggplot2::ylab("") +
ggplot2::scale_x_continuous(breaks = pretty(plotdata$periods)) +
ggplot2::theme_bw() +
ggplot2::labs(subtitle = paste0("OIRF and ",cibootstrap_results$CI*100,"%", " confidence bands"))
}
p
}
#' @export
girf.pvarfeols <- function(model, n.ahead, ma_approx_steps) {
Phi = coef(model)
# Phi is a function for the MA representaton in the package vars.
MA_Phi <- ma_phi_representation(Phi = Phi,
ma_approx_steps = ma_approx_steps,
lags = model$lags)
Sigma_Hat1 <- vcov(model)
sigmas <- sqrt(diag(Sigma_Hat1))
girf_output_time <- list()
for (i0 in 1:n.ahead){
girf_output_time[[i0]] <- t( MA_Phi[[i0]] %*% Sigma_Hat1) / sigmas
}
girf_output_shock <- list()
for (i0 in 1:nrow(Sigma_Hat1)){
zwischen <- matrix(NA, nrow = n.ahead, ncol = nrow(Sigma_Hat1))
for (i1 in 1:length(girf_output_time)){
zwischen[i1,] <- girf_output_time[[i1]][i0,]
}
girf_output_shock[[i0]] <- zwischen
colnames(girf_output_shock[[i0]]) <- model$dependent_vars
}
names(girf_output_shock) <- model$dependent_vars
class(girf_output_shock) <- "pvarfeolsgirf"
return(girf_output_shock)
}
#' @export
plot.pvarfeolsgirf <- function(x, cibootstrap_results, ...) {
dependent_vars <- names(x)
steps <- nrow(x[[1]])
plotdata <- do.call(rbind,lapply(x, reshape2::melt))
names(plotdata) <- c("periods", "impulse", "value")
plotdata$impulse <- paste(rep(dependent_vars, each = steps*length(dependent_vars)), "on", plotdata$impulse)
p <-
ggplot2::ggplot(plotdata, ggplot2::aes_(x = ~periods, y = ~value)) +
ggplot2::geom_line(colour = "orangered4") + ggplot2::facet_wrap(~impulse) +
ggplot2::labs(title = "Generalized impulse response function") +
ggplot2::xlab("steps") + ggplot2::ylab("") +
ggplot2::scale_x_continuous(breaks = pretty(plotdata$periods)) +
ggplot2::theme_bw()
if (!missing(cibootstrap_results)) {
ci_data_upper <-
data.frame(do.call(rbind,lapply(cibootstrap_results$Upper, reshape2::melt)), Bound = "Upper")
names(ci_data_upper)[1:3] <- c("periods", "impulse", "value")
row.names(ci_data_upper) <- NULL
ci_data_upper$impulse <- paste(rep(dependent_vars,
each = steps*length(dependent_vars)), "on", ci_data_upper$impulse)
ci_data_lower <-
data.frame(do.call(rbind,lapply(cibootstrap_results$Lower, reshape2::melt)), Bound = "Lower")
names(ci_data_lower)[1:3] <- c("periods", "impulse", "value")
row.names(ci_data_lower) <- NULL
ci_data_lower$impulse <- paste(rep(dependent_vars,
each = steps*length(dependent_vars)), "on", ci_data_lower$impulse)
ci_data <- reshape2::dcast(rbind(ci_data_upper, ci_data_lower), periods + impulse ~ Bound, value.var = "value")
plotdata <-
merge(plotdata, ci_data, all.x = TRUE, by = c("periods", "impulse"), sort = FALSE)
p <-
ggplot2::ggplot(plotdata, ggplot2::aes_(x = ~periods, y = ~value)) +
ggplot2::geom_line(colour = "orangered4") +
ggplot2::geom_ribbon(ggplot2::aes_(x = ~periods,ymin = ~Lower, ymax = ~Upper), alpha = 0.3,fill = "steelblue2") +
ggplot2::facet_wrap(~impulse) +
ggplot2::labs(title = "Generalized impulse response function") +
ggplot2::xlab("steps") + ggplot2::ylab("") +
ggplot2::scale_x_continuous(breaks = pretty(plotdata$periods)) +
ggplot2::theme_bw() +
ggplot2::labs(subtitle = paste0("GIRF and ",cibootstrap_results$CI*100,"%", " confidence bands"))
}
p
}
#' @rdname stability
#' @export
stability.pvarfeols <- function(model, ...) {
Phi <- coef(model)
# Get the PVAR(1) represenation
A_full <- pvar1_phi_representation(Phi = Phi, lags = model$lags)
# Calculate the eigenvalues of the matrix A_full.
A_full_eigen <- eigen(A_full, only.values = TRUE)
res <- data.frame(Eigenvalue = A_full_eigen$values, Modulus = abs(A_full_eigen$values))
class(res) <- c("pvarstability", "data.frame")
return(res)
}
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.