MatthieuStigler / tsDyn

tsDyn
tsdyn.googlecode.com
GNU General Public License v2.0
34 stars 20 forks source link

tidy for irf #52

Open MatthieuStigler opened 1 year ago

MatthieuStigler commented 1 year ago

Maybe tsDyn could provide tidy.varirf? Admittedly, it would be better though if in vars!?

Here's some code:

suppressPackageStartupMessages(library(vars))
suppressPackageStartupMessages(library(ggplot2))
data(Canada)
var.2c <- VAR(Canada, p = 2, type = "const")
i <- irf(var.2c)

get_as_df <- function(x, long=FALSE, long_name = "irf") {

  N_T <- unique(sapply(x, nrow))
  if(length(N_T)>1) warning("UUps, problem")
  DF <- do.call("rbind", lapply(x, as.data.frame))
  DF$time <- rep(1:N_T, times= length(names(x)))
  DF$impulse <- rep(names(x), each=N_T)
  DF <- DF[,c(ncol(DF) + c(0, -1),1:(ncol(DF)-2))]

  if(long){
    DF <- DF |> 
      reshape(varying = 3:ncol(DF),
              times = colnames(DF)[3:ncol(DF)],
              v.names = long_name,
              timevar = "response",
              direction="long") 
    rownames(DF) <- NULL
    DF <- DF[,c("impulse", "response", "time", long_name)]
  }
  DF
}

# get_as_df(x=i$irf)
# get_as_df(x=i$irf, long=TRUE)
# get_as_df(x=i$Lower)
# get_as_df(x=i$Upper)

tidy.varirf <- function(x){
  irf <- get_as_df(x=i$irf, long=TRUE, long_name = "irf")
  lower <- get_as_df(x=i$Lower, long=TRUE, long_name = "conf.low")
  upper <- get_as_df(x=i$Upper, long=TRUE, long_name = "conf.high")
  Reduce(function(x, y) merge(x, y, all=TRUE, by = c("impulse", "response", "time")), list(irf, lower, upper))
}

tidy.varirf(i) |> head()
#>   impulse response time         irf   conf.low conf.high
#> 1       e        e    1  0.36281502  0.2836301 0.3927026
#> 2       e        e   10  0.04490000 -0.4518087 0.3914526
#> 3       e        e   11 -0.03568143 -0.5401913 0.3596537
#> 4       e        e    2  0.54753375  0.3792261 0.5956236
#> 5       e        e    3  0.61791814  0.3563412 0.6881031
#> 6       e        e    4  0.61135633  0.2631357 0.7109548

plot <- tidy.varirf(i) |>
  tibble::as_tibble() |>
  ggplot2::ggplot(aes(x = time, y=irf))+
  ggplot2::geom_hline(yintercept = 0, linetype=2)+
  ggplot2::geom_ribbon(aes(ymin = conf.low, ymax= conf.high),
              fill ="grey", alpha= 0.6)+
  ggplot2::geom_line()+
  ggplot2::facet_grid(impulse~response, switch = "y")

Created on 2023-10-05 with reprex v2.0.2