From 57d873d3a17ae147c47f0252be0149cd780ce12f Mon Sep 17 00:00:00 2001
From: Young Geun Kim <ygeunkimstat@gmail.com>
Date: Mon, 23 Dec 2024 19:53:38 +0900
Subject: [PATCH 1/2] require r 4.1 following tidyverse schedule

---
 DESCRIPTION               |   2 +-
 NAMESPACE                 |   2 -
 R/criteria.R              |  66 +++----
 R/forecast.R              |  94 +++++-----
 R/generate-process.R      |   4 +-
 R/irf.R                   |  16 +-
 R/misc-r.R                | 126 ++++++++-----
 R/plot-spillover.R        |  14 +-
 R/plot.R                  | 113 ++++++------
 R/print-spillover.R       |  18 +-
 R/print-varlse.R          |   4 +-
 R/print-vharlse.R         |   8 +-
 R/spillover.R             |  40 ++---
 R/stable-process.R        |  40 ++---
 R/summary-bayes.R         |  12 +-
 R/summary-forecast.R      | 360 ++++++++++----------------------------
 R/summary-sparse.R        |  14 +-
 R/summary-varlse.R        |   8 +-
 R/summary-vharlse.R       |   8 +-
 R/utils-pipe.R            |  14 --
 data-raw/etf_vix.R        |  10 +-
 man/pipe.Rd               |  20 ---
 vignettes/bvhar.Rmd       |   2 +-
 vignettes/forecasting.Rmd |  24 +--
 24 files changed, 423 insertions(+), 596 deletions(-)
 delete mode 100644 R/utils-pipe.R
 delete mode 100644 man/pipe.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 7d5a84d2..4786fd3e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -50,6 +50,6 @@ LinkingTo:
     RcppThread
 VignetteBuilder: knitr
 Depends: 
-    R (>= 4.0.0)
+    R (>= 4.1.0)
 Roxygen: list(markdown = TRUE)
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index 3156db07..c786e426 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -180,7 +180,6 @@ S3method(summary,normaliw)
 S3method(summary,ssvsmod)
 S3method(summary,varlse)
 S3method(summary,vharlse)
-export("%>%")
 export(FPE)
 export(HQ)
 export(VARtoVMA)
@@ -320,7 +319,6 @@ importFrom(lifecycle,deprecate_soft)
 importFrom(lifecycle,deprecate_warn)
 importFrom(lifecycle,deprecated)
 importFrom(lifecycle,is_present)
-importFrom(magrittr,"%>%")
 importFrom(optimParallel,optimParallel)
 importFrom(posterior,as_draws_df)
 importFrom(posterior,as_draws_matrix)
diff --git a/R/criteria.R b/R/criteria.R
index b66e8800..d7402db8 100644
--- a/R/criteria.R
+++ b/R/criteria.R
@@ -168,8 +168,8 @@ logLik.bvharmn <- function(object, ...) {
 #' @importFrom stats AIC
 #' @export
 AIC.varlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -178,8 +178,8 @@ AIC.varlse <- function(object, ...) {
 #' @param ... not used
 #' @export
 AIC.vharlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -188,8 +188,8 @@ AIC.vharlse <- function(object, ...) {
 #' @param ... not used
 #' @export
 AIC.bvarmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -198,8 +198,8 @@ AIC.bvarmn <- function(object, ...) {
 #' @param ... not used
 #' @export
 AIC.bvarflat <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -208,8 +208,8 @@ AIC.bvarflat <- function(object, ...) {
 #' @param ... not used
 #' @export
 AIC.bvharmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -265,40 +265,40 @@ FPE.vharlse <- function(object, ...) {
 #' @importFrom stats BIC
 #' @export
 BIC.varlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     BIC()
 }
 
 #' @rdname vhar_lm
 #' @export
 BIC.vharlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     BIC()
 }
 
 #' @rdname bvar_minnesota
 #' @export
 BIC.bvarmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     BIC()
 }
 
 #' @rdname bvar_flat
 #' @export
 BIC.bvarflat <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     BIC()
 }
 
 #' @rdname bvhar_minnesota
 #' @export
 BIC.bvharmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     BIC()
 }
 
@@ -351,40 +351,40 @@ HQ.logLik <- function(object, ...) {
 #' Quinn, B.G. (1980). *Order Determination for a Multivariate Autoregression*. Journal of the Royal Statistical Society: Series B (Methodological), 42: 182-185.
 #' @export
 HQ.varlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     HQ()
 }
 
 #' @rdname HQ
 #' @export
 HQ.vharlse <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     HQ()
 }
 
 #' @rdname HQ
 #' @export
 HQ.bvarmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     HQ()
 }
 
 #' @rdname HQ
 #' @export
 HQ.bvarflat <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     HQ()
 }
 
 #' @rdname HQ
 #' @export
 HQ.bvharmn <- function(object, ...) {
-  object %>% 
-    logLik() %>% 
+  object |> 
+    logLik() |> 
     AIC()
 }
 
@@ -427,8 +427,8 @@ compute_dic.bvarmn <- function(object, n_iter = 100L, ...) {
   bmat_gen <- rand_gen$coefficients
   covmat_gen <- rand_gen$covmat
   log_lik <- 
-    object %>% 
-    logLik() %>% 
+    object |> 
+    logLik() |> 
     as.numeric()
   obs <- object$obs
   m <- object$m
@@ -446,7 +446,7 @@ compute_dic.bvarmn <- function(object, n_iter = 100L, ...) {
             )
           ) / 2
       }
-    ) %>% 
+    ) |> 
     unlist()
   eff_num <- 2 * (log_lik - mean(post_mean))
   -2 * log_lik + 2 * eff_num
diff --git a/R/forecast.R b/R/forecast.R
index 8e7ff27c..ad29536b 100644
--- a/R/forecast.R
+++ b/R/forecast.R
@@ -56,9 +56,9 @@ predict.varlse <- function(object, n_ahead, level = .05, ...) {
   pred_res <- forecast_var(object, n_ahead)
   colnames(pred_res) <- colnames(object$y0)
   SE <- 
-    compute_covmse(object, n_ahead) %>% # concatenated matrix
-    split.data.frame(gl(n_ahead, object$m)) %>% # list of forecast MSE covariance matrix
-    sapply(diag) %>% 
+    compute_covmse(object, n_ahead) |> # concatenated matrix
+    split.data.frame(gl(n_ahead, object$m)) |> # list of forecast MSE covariance matrix
+    sapply(diag) |> 
     t() # extract only diagonal element to compute CIs
   SE <- sqrt(SE)
   colnames(SE) <- colnames(object$y0)
@@ -109,9 +109,9 @@ predict.vharlse <- function(object, n_ahead, level = .05, ...) {
   pred_res <- forecast_vhar(object, n_ahead)
   colnames(pred_res) <- colnames(object$y0)
   SE <- 
-    compute_covmse_har(object, n_ahead) %>% # concatenated matrix
-    split.data.frame(gl(n_ahead, object$m)) %>% # list of forecast MSE covariance matrix
-    sapply(diag) %>% 
+    compute_covmse_har(object, n_ahead) |> # concatenated matrix
+    split.data.frame(gl(n_ahead, object$m)) |> # list of forecast MSE covariance matrix
+    sapply(diag) |> 
     t() # extract only diagonal element to compute CIs
   SE <- sqrt(SE)
   colnames(SE) <- colnames(object$y0)
@@ -190,12 +190,12 @@ predict.bvarmn <- function(object, n_ahead, n_iter = 100L, level = .05, num_thre
   # Predictive distribution-------------------------
   dim_data <- ncol(pred_mean)
   y_distn <-
-    pred_res$predictive %>%
+    pred_res$predictive |>
     array(dim = c(n_ahead, dim_data, n_iter)) # 3d array: h x m x B
   # num_draw <- nrow(alpha_record) # concatenate multiple chains
   # y_distn <-
-  #   pred_res %>%
-  #   unlist() %>%
+  #   pred_res |>
+  #   unlist() |>
   #   array(dim = c(n_ahead, dim_data, num_draw))
   # pred_mean <- apply(y_distn, c(1, 2), mean)
   lower_quantile <- apply(y_distn, c(1, 2), quantile, probs = level / 2)
@@ -271,12 +271,12 @@ predict.bvharmn <- function(object, n_ahead, n_iter = 100L, level = .05, num_thr
   # Predictive distribution-------------------------
   dim_data <- ncol(pred_mean)
   y_distn <- 
-    pred_res$predictive %>% 
+    pred_res$predictive |> 
     array(dim = c(n_ahead, dim_data, n_iter)) # 3d array: h x m x B
   # num_draw <- nrow(phi_record) # concatenate multiple chains
   # y_distn <-
-  #   pred_res %>%
-  #   unlist() %>%
+  #   pred_res |>
+  #   unlist() |>
   #   array(dim = c(n_ahead, dim_data, num_draw))
   # pred_mean <- apply(y_distn, c(1, 2), mean)
   lower_quantile <- apply(y_distn, c(1, 2), quantile, probs = level / 2)
@@ -343,12 +343,12 @@ predict.bvarflat <- function(object, n_ahead, n_iter = 100L, level = .05, num_th
   # Predictive distribution-------------------------
   dim_data <- ncol(pred_mean)
   y_distn <-
-    pred_res$predictive %>%
+    pred_res$predictive |>
     array(dim = c(n_ahead, dim_data, n_iter)) # 3d array: h x m x B
   # num_draw <- nrow(alpha_record) # concatenate multiple chains
   # y_distn <-
-  #   pred_res %>%
-  #   unlist() %>%
+  #   pred_res |>
+  #   unlist() |>
   #   array(dim = c(n_ahead, dim_data, num_draw))
   # pred_mean <- apply(y_distn, c(1, 2), mean)
   lower_quantile <- apply(y_distn, c(1, 2), quantile, probs = level / 2)
@@ -399,13 +399,11 @@ predict.bvarldlt <- function(object, n_ahead, level = .05, stable = FALSE, num_t
       alpha_record,
       1,
       function(x) {
-        all(
-          matrix(x, ncol = object$m) %>%
-            compute_stablemat() %>%
-            eigen() %>%
-            .$values %>%
-            Mod() < 1
-        )
+        eigen_vals <-
+          matrix(x, ncol = object$m) |>
+          compute_stablemat() |>
+          eigen()
+        all(Mod(eigen_vals$values) < 1)
       }
     )
     if (any(!is_stable)) {
@@ -458,8 +456,8 @@ predict.bvarldlt <- function(object, n_ahead, level = .05, stable = FALSE, num_t
   # Predictive distribution------------------------------------
   num_draw <- nrow(alpha_record) # concatenate multiple chains
   y_distn <-
-    pred_res %>%
-    unlist() %>%
+    pred_res |>
+    unlist() |>
     array(dim = c(n_ahead, dim_data, num_draw))
   if (med) {
     pred_mean <- apply(y_distn, c(1, 2), median)
@@ -512,13 +510,11 @@ predict.bvharldlt <- function(object, n_ahead, level = .05, stable = FALSE, num_
       1,
       function(x) {
         coef <- t(object$HARtrans[1:(object$p * dim_data), 1:(object$month * dim_data)]) %*% matrix(x, ncol = object$m)
-        all(
-          coef %>%
-            compute_stablemat() %>%
-            eigen() %>%
-            .$values %>%
-            Mod() < 1
-        )
+        eigen_vals <-
+          coef |>
+          compute_stablemat() |>
+          eigen()
+        all(Mod(eigen_vals$values) < 1)
       }
     )
     if (any(!is_stable)) {
@@ -571,8 +567,8 @@ predict.bvharldlt <- function(object, n_ahead, level = .05, stable = FALSE, num_
   # Predictive distribution------------------------------------
   num_draw <- nrow(phi_record) # concatenate multiple chains
   y_distn <-
-    pred_res %>% 
-    unlist() %>% 
+    pred_res |> 
+    unlist() |> 
     array(dim = c(n_ahead, dim_data, num_draw))
   if (med) {
     pred_mean <- apply(y_distn, c(1, 2), median)
@@ -629,13 +625,11 @@ predict.bvarsv <- function(object, n_ahead, level = .05, stable = FALSE, num_thr
       alpha_record,
       1,
       function(x) {
-        all(
-          matrix(x, ncol = object$m) %>%
-            compute_stablemat() %>%
-            eigen() %>%
-            .$values %>%
-            Mod() < 1
-        )
+        eigen_vals <-
+          matrix(x, ncol = object$m) |>
+          compute_stablemat() |>
+          eigen()
+        all(Mod(eigen_vals$values) < 1)
       }
     )
     if (any(!is_stable)) {
@@ -689,8 +683,8 @@ predict.bvarsv <- function(object, n_ahead, level = .05, stable = FALSE, num_thr
   # Predictive distribution------------------------------------
   num_draw <- nrow(alpha_record) # concatenate multiple chains
   y_distn <-
-    pred_res %>% 
-    unlist() %>% 
+    pred_res |> 
+    unlist() |> 
     array(dim = c(n_ahead, dim_data, num_draw))
   if (med) {
     pred_mean <- apply(y_distn, c(1, 2), median)
@@ -744,13 +738,11 @@ predict.bvharsv <- function(object, n_ahead, level = .05, stable = FALSE, num_th
       1,
       function(x) {
         coef <- t(object$HARtrans[1:(object$p * dim_data), 1:(object$month * dim_data)]) %*% matrix(x, ncol = object$m)
-        all(
-          coef %>%
-            compute_stablemat() %>%
-            eigen() %>%
-            .$values %>%
-            Mod() < 1
-        )
+        eigen_vals <-
+          coef |>
+          compute_stablemat() |>
+          eigen()
+        all(Mod(eigen_vals$values) < 1)
       }
     )
     if (any(!is_stable)) {
@@ -805,8 +797,8 @@ predict.bvharsv <- function(object, n_ahead, level = .05, stable = FALSE, num_th
   # Predictive distribution------------------------------------
   num_draw <- nrow(phi_record) # concatenate multiple chains
   y_distn <-
-    pred_res %>% 
-    unlist() %>% 
+    pred_res |> 
+    unlist() |> 
     array(dim = c(n_ahead, dim_data, num_draw))
   if (med) {
     pred_mean <- apply(y_distn, c(1, 2), median)
diff --git a/R/generate-process.R b/R/generate-process.R
index b1c31219..8f0ed2ce 100644
--- a/R/generate-process.R
+++ b/R/generate-process.R
@@ -132,8 +132,8 @@ sim_var <- function(num_sim,
 #' @export
 sim_mniw <- function(num_sim, mat_mean, mat_scale_u, mat_scale, shape, u_prec = FALSE) {
   res <-
-    sim_mniw_export(num_sim, mat_mean, mat_scale_u, mat_scale, shape, u_prec) %>%
-    simplify2array() %>%
+    sim_mniw_export(num_sim, mat_mean, mat_scale_u, mat_scale, shape, u_prec) |>
+    simplify2array() |>
     apply(1, function(x) x)
   names(res) <- c("mn", "iw")
   res
diff --git a/R/irf.R b/R/irf.R
index e89e0722..6bd40db8 100644
--- a/R/irf.R
+++ b/R/irf.R
@@ -75,17 +75,17 @@ irf.varlse <- function(object,
   )
   res <- list(coefficients = mat_irf)
   res$df_long <- 
-    mat_irf %>% 
-    as.data.frame() %>% 
+    mat_irf |> 
+    as.data.frame() |> 
     mutate(
       impulse = impulse_name,
       period = period_name
-    ) %>% 
+    ) |> 
     pivot_longer(
       -c(period, impulse),
       names_to = "response",
       values_to = "value"
-    ) %>% 
+    ) |> 
     filter(impulse %in% impulse_var, response %in% response_var)
   # return----------------------
   res$lag_max <- lag_max
@@ -144,17 +144,17 @@ irf.vharlse <- function(object,
   )
   res <- list(coefficients = mat_irf)
   res$df_long <- 
-    mat_irf %>% 
-    as.data.frame() %>% 
+    mat_irf |> 
+    as.data.frame() |> 
     mutate(
       impulse = impulse_name,
       period = period_name
-    ) %>% 
+    ) |> 
     pivot_longer(
       -c(period, impulse),
       names_to = "response",
       values_to = "value"
-    ) %>% 
+    ) |> 
     filter(impulse %in% impulse_var, response %in% response_var)
   # return----------------------
   res$lag_max <- lag_max
diff --git a/R/misc-r.R b/R/misc-r.R
index 5e567d7f..a1f2c739 100644
--- a/R/misc-r.R
+++ b/R/misc-r.R
@@ -4,7 +4,7 @@ concatenate_colnames <- function(var_name, prefix, include_mean = TRUE) {
     lapply(
       prefix,
       function(lag) paste(var_name, lag, sep = "_")
-    ) %>% 
+    ) |> 
     unlist()
   if (!include_mean) {
     return(nm)
@@ -33,11 +33,11 @@ split_coef <- function(object, ...) {
     return(
       switch(object$type,
         "const" = {
-          split.data.frame(object$coefficients[-object$df, ], gl(object$p, object$m)) %>%
+          split.data.frame(object$coefficients[-object$df, ], gl(object$p, object$m)) |>
             lapply(t)
         },
         "none" = {
-          split.data.frame(object$coefficients, gl(object$p, object$m)) %>%
+          split.data.frame(object$coefficients, gl(object$p, object$m)) |>
             lapply(t)
         }
       )
@@ -68,9 +68,9 @@ split_coef <- function(object, ...) {
 split_paramarray <- function(x, chain, param_name) {
   num_var <- ncol(x) / chain
   res <- 
-    split.data.frame(t(x), gl(num_var, 1, ncol(x))) %>%
-    lapply(t) %>%
-    unlist() %>% 
+    split.data.frame(t(x), gl(num_var, 1, ncol(x))) |>
+    lapply(t) |>
+    unlist() |> 
     array(
       dim = c(nrow(x), chain, num_var),
       dimnames = list(
@@ -91,7 +91,7 @@ split_paramarray <- function(x, chain, param_name) {
 #' @noRd
 split_psirecord <- function(x, chain = 1, varname = "cholesky") {
   res <- 
-    x %>% 
+    x |> 
     split.data.frame(gl(nrow(x) / ncol(x), ncol(x)))
   if (chain == 1) {
     return(res)
@@ -100,9 +100,9 @@ split_psirecord <- function(x, chain = 1, varname = "cholesky") {
       res,
       function(y) {
         num_var <- ncol(y) / chain
-        split.data.frame(t(y), gl(num_var, 1, ncol(y))) %>% 
-          lapply(t) %>% 
-          unlist() %>% 
+        split.data.frame(t(y), gl(num_var, 1, ncol(y))) |> 
+          lapply(t) |> 
+          unlist() |> 
           array(
             dim = c(nrow(y), chain, num_var),
             dimnames = list(
@@ -132,10 +132,10 @@ split_chain <- function(x, chain = 1, varname = "alpha") {
     # num_var <- ncol(x) / chain
     num_row <- nrow(x) / chain
     res <-
-      # split.data.frame(t(x), gl(num_var, 1, ncol(x))) %>%
-      # lapply(t) %>%
-      split.data.frame(x, gl(chain, num_row)) %>%
-      unlist(x) %>%
+      # split.data.frame(t(x), gl(num_var, 1, ncol(x))) |>
+      # lapply(t) |>
+      split.data.frame(x, gl(chain, num_row)) |>
+      unlist(x) |>
       array(
         # dim = c(nrow(x), chain, num_var),
         dim = c(num_row, chain, ncol(x)),
@@ -180,21 +180,67 @@ get_gammaparam <- function(mode, sd) {
 #' @param var_names Variable names
 #' @param level level for lower and upper quantiles
 #' @param med Get median instead of mean?
+#' @param roll Is the `draws` the result of rolling or expanding windows?
 #' 
 #' @noRd 
-process_forecast_draws <- function(draws, n_ahead, dim_data, num_draw, var_names, level = .05, med = FALSE) {
-  mcmc_distn <-
-    draws %>%
-    unlist() %>% 
-    array(dim = c(n_ahead, dim_data, num_draw))
-  if (med) {
-    pred_mean <- apply(mcmc_distn, c(1, 2), median)
+process_forecast_draws <- function(draws, n_ahead, dim_data, num_draw, var_names, level = .05, roll = FALSE, med = FALSE) {
+  if (roll) {
+    if (med) {
+      pred_mean <-
+        draws |>
+        lapply(function(res) {
+          unlist(res) |>
+            array(dim = c(n_ahead, dim_data, num_draw)) |>
+            apply(c(1, 2), median)
+        })
+    } else {
+      pred_mean <-
+        draws |>
+        lapply(function(res) {
+          unlist(res) |>
+            array(dim = c(n_ahead, dim_data, num_draw)) |>
+            apply(c(1, 2), mean)
+        })
+    }
+    pred_mean <- do.call(rbind, pred_mean)
+    pred_se <-
+      draws |>
+      lapply(function(res) {
+        unlist(res) |>
+          array(dim = c(n_ahead, dim_data, num_draw)) |>
+          apply(c(1, 2), sd)
+      })
+    pred_se <- do.call(pred_se, pred_se)
+    pred_lower <-
+      draws |> 
+      lapply(function(res) {
+        unlist(res) |> 
+          array(dim = c(n_ahead, dim_data, num_draw)) |> 
+          apply(c(1, 2), quantile, probs = level / 2)
+      })
+    pred_lower <- do.call(rbind, pred_lower)
+    pred_upper <-
+      draws |>
+      lapply(function(res) {
+        unlist(res) |>
+          array(dim = c(n_ahead, dim_data, num_draw)) |>
+          apply(c(1, 2), quantile, probs = 1 - level / 2)
+      })
+    pred_upper <- do.call(rbind, pred_upper)
   } else {
-    pred_mean <- apply(mcmc_distn, c(1, 2), mean)
+    mcmc_distn <-
+      draws |>
+      unlist() |>
+      array(dim = c(n_ahead, dim_data, num_draw))
+    if (med) {
+      pred_mean <- apply(mcmc_distn, c(1, 2), median)
+    } else {
+      pred_mean <- apply(mcmc_distn, c(1, 2), mean)
+    }
+    pred_se <- apply(mcmc_distn, c(1, 2), sd)
+    pred_lower <- apply(mcmc_distn, c(1, 2), quantile, probs = level / 2)
+    pred_upper <- apply(mcmc_distn, c(1, 2), quantile, probs = 1 - level / 2)
   }
-  pred_se <- apply(mcmc_distn, c(1, 2), sd)
-  pred_lower <- apply(mcmc_distn, c(1, 2), quantile, probs = level / 2)
-  pred_upper <- apply(mcmc_distn, c(1, 2), quantile, probs = 1 - level / 2)
   colnames(pred_mean) <- var_names
   rownames(pred_mean) <- var_names
   colnames(pred_se) <- var_names
@@ -244,16 +290,16 @@ process_vector_draws <- function(draws, dim_data, level = .05, med = FALSE) {
 #' @importFrom dplyr mutate
 #' @noRd 
 process_dynamic_spdraws <- function(draws, dim_data, level = .05, med = FALSE, var_names) {
-  lapply(
+  sp_draws <- lapply(
     draws,
     function(x) {
-      process_vector_draws(unlist(x), dim_data = dim_data, level = level, med = med) %>%
-        do.call(cbind, .) %>%
-        as.data.frame() %>%
+      process_vector_draws(unlist(x), dim_data = dim_data, level = level, med = med) |>
+        do.call(cbind, .) |>
+        as.data.frame() |>
         mutate(series = var_names)
     }
-  ) %>%
-    do.call(rbind, .) %>%
+  )
+  do.call(rbind, sp_draws) |> 
     as_tibble()
 }
 
@@ -263,9 +309,9 @@ process_dynamic_spdraws <- function(draws, dim_data, level = .05, med = FALSE, v
 #' @param col_names Column name for value
 #' @noRd 
 gather_spillover <- function(connect, col_names = "spillover") {
-  connect %>%
-    as.data.frame() %>%
-    rownames_to_column(var = "series") %>%
+  connect |>
+    as.data.frame() |>
+    rownames_to_column(var = "series") |>
     pivot_longer(-"series", names_to = "shock", values_to = col_names)
 }
 
@@ -276,9 +322,9 @@ gather_spillover <- function(connect, col_names = "spillover") {
 #' 
 #' @noRd
 join_long_spillover <- function(connect, prefix = "spillover") {
-  gather_spillover(connect$mean, col_names = prefix) %>%
-    left_join(gather_spillover(connect$lower, col_names = paste(prefix, "lower", sep = "_")), by = c("series", "shock")) %>%
-    left_join(gather_spillover(connect$upper, col_names = paste(prefix, "upper", sep = "_")), by = c("series", "shock")) %>%
+  gather_spillover(connect$mean, col_names = prefix) |>
+    left_join(gather_spillover(connect$lower, col_names = paste(prefix, "lower", sep = "_")), by = c("series", "shock")) |>
+    left_join(gather_spillover(connect$upper, col_names = paste(prefix, "upper", sep = "_")), by = c("series", "shock")) |>
     left_join(gather_spillover(connect$sd, col_names = paste(prefix, "sd", sep = "_")), by = c("series", "shock"))
 }
 
@@ -354,10 +400,10 @@ get_records <- function(object, split_chain = TRUE) {
   lapply(
     object$param_names,
     function(x) {
-      subset_draws(object$param, variable = x) %>%
-        as_draws_matrix() %>%
+      subset_draws(object$param, variable = x) |>
+        as_draws_matrix() |>
         split.data.frame(gl(num_chains, nrow(object$param) / num_chains))
     }
-  ) %>%
+  ) |>
     setNames(paste(object$param_names, "record", sep = "_"))
 }
diff --git a/R/plot-spillover.R b/R/plot-spillover.R
index d5f211df..bb9ba511 100644
--- a/R/plot-spillover.R
+++ b/R/plot-spillover.R
@@ -25,29 +25,29 @@ autoplot.bvhardynsp <- function(object,
       data.frame(
         id = object$index,
         y = object$tot
-      ) %>%
+      ) |>
         ggplot(aes(x = id, y = y)) +
           geom_path()
     },
     "to" = {
-      cbind(id = object$index, object$to) %>%
-        pivot_longer(-id, names_to = "series", values_to = "value") %>%
+      cbind(id = object$index, object$to) |>
+        pivot_longer(-id, names_to = "series", values_to = "value") |>
         ggplot(aes(x = id)) +
         # geom_ribbon(aes(ymin = 0, ymax = value)) +
         geom_path(aes(y = value)) +
         facet_wrap(series ~ ., nrow = row_facet, ncol = col_facet)
     },
     "from" = {
-      cbind(id = object$index, object$from) %>%
-        pivot_longer(-id, names_to = "series", values_to = "value") %>%
+      cbind(id = object$index, object$from) |>
+        pivot_longer(-id, names_to = "series", values_to = "value") |>
         ggplot(aes(x = id)) +
         # geom_ribbon(aes(ymin = 0, ymax = value)) +
         geom_path(aes(y = value)) +
         facet_wrap(series ~ ., nrow = row_facet, ncol = col_facet)
     },
     "net" = {
-      cbind(id = object$index, object$net) %>%
-        pivot_longer(-id, names_to = "series", values_to = "value") %>%
+      cbind(id = object$index, object$net) |>
+        pivot_longer(-id, names_to = "series", values_to = "value") |>
         ggplot(aes(x = id)) +
         # geom_ribbon(aes(ymin = 0, ymax = value)) +
         geom_hline(yintercept = 0, col = hcol, size = hsize) +
diff --git a/R/plot.R b/R/plot.R
index c9e0f69a..2b6bfb94 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -34,12 +34,12 @@ autoplot.summary.normaliw <- function(object, type = c("trace", "dens", "area"),
 #' @importFrom tidyr pivot_longer
 #' @export
 autoplot.normaliw <- function(object, hcol = "grey", hsize = 1.5, ...) {
-  X <- object$residuals %>% as.data.frame()
+  X <- object$residuals |> as.data.frame()
   X[["id"]] <- 1:object$obs
   X <- 
-    X %>% 
+    X |> 
     pivot_longer(-id, names_to = "name", values_to = "value")
-  X %>% 
+  X |> 
     ggplot(aes(x = id, y = value)) +
     geom_hline(yintercept = 0, col = hcol, size = hsize) +
     geom_point(...) +
@@ -70,27 +70,27 @@ autoplot.normaliw <- function(object, hcol = "grey", hsize = 1.5, ...) {
 #' @noRd
 gather_predbvhar <- function(object) {
   Y <- 
-    object$y %>% 
-    as.data.frame() %>% 
+    object$y |> 
+    as.data.frame() |> 
     mutate(forecast = FALSE)
   PRED <- 
-    object$forecast %>% 
-    as.data.frame() %>% 
+    object$forecast |> 
+    as.data.frame() |> 
     mutate(forecast = TRUE)
   lapply(
     c("forecast", "lower_joint", "upper_joint"),
     function(comp) {
       PRED <- 
-        object[[comp]] %>% 
-        as.data.frame() %>% 
+        object[[comp]] |> 
+        as.data.frame() |> 
         mutate(forecast = TRUE)
-      Y %>% 
-        bind_rows(PRED) %>% 
-        mutate(id = 1:n()) %>% 
+      Y |> 
+        bind_rows(PRED) |> 
+        mutate(id = 1:n()) |> 
         pivot_longer(-c(id, forecast), names_to = "variable", values_to = paste("value", comp, sep = "_"))
     }
-  ) %>% 
-    reduce(left_join, by = c("id", "variable", "forecast")) %>% 
+  ) |> 
+    reduce(left_join, by = c("id", "variable", "forecast")) |> 
     mutate(Model = object$process)
 }
 
@@ -182,10 +182,10 @@ autoplot.predbvhar <- function(object,
                                NCOL = NULL, ...) {
   type <- match.arg(type)
   forecast_list <- 
-    gather_predbvhar(object) %>% 
+    gather_predbvhar(object) |> 
     filter(id >= x_cut)
   p <- 
-    forecast_list %>% 
+    forecast_list |> 
     ggplot(aes(x = id, y = value_forecast)) +
     geom_path(data = filter(forecast_list, forecast == FALSE))
   p <- 
@@ -239,17 +239,16 @@ autolayer.predbvhar <- function(object,
                                 ci_fill = "grey70", 
                                 ci_alpha = .5,
                                 alpha_scale = .3, ...) {
-  aes_data <- 
-    last_plot() %>% 
-    ggplot_build() %>% 
-    .$plot %>% 
-    .$data # same form as forecast_list above
+  aes_data <-
+    last_plot() |>
+    ggplot_build()
+  aes_data <- aes_data$plot$data # same form as forecast_list above
   x_cut <- aes_data[["id"]][1]
   NEW_list <- 
-    gather_predbvhar(object) %>% # new forecast_list
+    gather_predbvhar(object) |> # new forecast_list
     filter(id >= x_cut)
   geom_predbvhar(
-    data = NEW_list %>% filter(forecast == TRUE),
+    data = NEW_list |> filter(forecast == TRUE),
     alpha_scale = alpha_scale,
     ci_param = list(alpha = ci_alpha, colour = NA),
     line_param = list(...)
@@ -279,15 +278,15 @@ geom_eval <- function(data, colour = "red", ...) {
     stop("Every column must be numeric class.")
   }
   aes_id <-
-    last_plot() %>%
-    ggplot_build() %>%
-    .$plot %>%
-    .$data %>% # same form as forecast_list above
-    filter(forecast == FALSE) %>%
-    .$id # index
+    last_plot() |>
+    ggplot_build()
+  aes_id <-
+    aes_id$plot$data |>
+    filter(forecast == FALSE)
+  aes_id <- aes_id$id # index
   new_data <- 
-    data %>% 
-    mutate(id = 1:n() + max(aes_id)) %>% 
+    data |> 
+    mutate(id = 1:n() + max(aes_id)) |> 
     pivot_longer(-id, names_to = "variable", values_to = "value")
   geom_path(
     aes(x = id, y = value),
@@ -319,34 +318,34 @@ gather_loss <- function(object, y_test, loss = c("mse", "mae", "mape", "mase"))
   }
   # Model names-------------------------
   mod_name <- 
-    object %>% 
-    lapply(function(PRED) PRED$process) %>% 
+    object |> 
+    lapply(function(PRED) PRED$process) |> 
     unlist()
   # error for each model----------------
   score_dt <- 
     switch(
       loss,
       "mse" = {
-        object %>% 
+        object |> 
           lapply(mse, y_test)
       },
       "mae" = {
-        object %>% 
+        object |> 
           lapply(mae, y_test)
       },
       "mape" = {
-        object %>% 
+        object |> 
           lapply(mape, y_test)
       },
       "mase" = {
-        object %>% 
+        object |> 
           lapply(mase, y_test)
       }
     )
-  score_dt %>% 
-    bind_rows() %>% 
-    mutate(Model = mod_name) %>% 
-    pivot_longer(-Model, names_to = "name", values_to = "score") %>% 
+  score_dt |> 
+    bind_rows() |> 
+    mutate(Model = mod_name) |> 
+    pivot_longer(-Model, names_to = "name", values_to = "score") |> 
     mutate(Loss = toupper(loss))
 }
 
@@ -356,8 +355,8 @@ gather_loss <- function(object, y_test, loss = c("mse", "mae", "mape", "mase"))
 #' @importFrom dplyr group_by mutate
 #' @noRd
 summarise_loss <- function(data) {
-  data %>% 
-    group_by(Model, Loss) %>% 
+  data |> 
+    group_by(Model, Loss) |> 
     mutate(average = mean(score))
 }
 
@@ -482,13 +481,13 @@ gg_loss <- function(mod_list,
   # Input data for geom_loss----------------
   fct_arrange <- toupper(type)
   data <- 
-    type %>% 
+    type |> 
     lapply(
       function(loss) {
         gather_loss(mod_list, y, loss)
       }
-    ) %>% 
-    bind_rows() %>% 
+    ) |> 
+    bind_rows() |> 
     mutate(Loss = factor(Loss, levels = fct_arrange))
   # plot------------------------------------
   p <- 
@@ -534,8 +533,8 @@ gg_loss <- function(mod_list,
 #' @export
 autoplot.bvharirf <- function(object, ...) {
   irf_df <- object$df_long
-  irf_df %>% 
-    unite("term", impulse, response, c("impulse", "response"), sep = "->") %>% 
+  irf_df |> 
+    unite("term", impulse, response, c("impulse", "response"), sep = "->") |> 
     ggplot(aes(x = period, y = value)) +
     geom_path(...) +
     scale_x_continuous(breaks = 0:(object$lag_max)) +
@@ -584,15 +583,15 @@ autoplot.bvharsp <- function(object,
 #' @noRd
 gather_heat <- function(object) {
   heat_coef <-
-    object$coefficients %>%
-    as.data.frame() %>%
-    rownames_to_column("term") %>%
-    pivot_longer(-term, names_to = "x", values_to = "value") %>%
+    object$coefficients |>
+    as.data.frame() |>
+    rownames_to_column("term") |>
+    pivot_longer(-term, names_to = "x", values_to = "value") |>
     mutate(x = factor(x, levels = colnames(object$coefficients)))
   is_vhar <- gsub(pattern = "(?=\\_).*", replacement = "", object$process, perl = TRUE) == "VHAR"
   if (object$type == "const") {
     heat_coef <-
-      heat_coef %>% 
+      heat_coef |> 
       mutate(term = ifelse(
         term == "const",
         paste("const", term, sep = "_"),
@@ -600,16 +599,16 @@ gather_heat <- function(object) {
       ))
   }
   heat_coef <-
-    heat_coef %>% 
+    heat_coef |> 
     separate_wider_regex(
       term,
       patterns = c(y = ".*", "_", ord = ".*$")
-    ) %>% 
+    ) |> 
     mutate(y = factor(y, levels = c(rev(colnames(object$coefficients)), "const")))
   # VHAR model--------------------------------
   if (is_vhar) {
     heat_coef <- 
-      heat_coef %>% 
+      heat_coef |> 
       mutate(
         ord = case_when(
           ord == "day" ~ "Daily",
@@ -636,7 +635,7 @@ gather_heat <- function(object) {
 autoplot.summary.bvharsp <- function(object, point = FALSE, ...) {
   heat_coef <- gather_heat(object)
   p <- 
-    heat_coef %>% 
+    heat_coef |> 
     ggplot(aes(x = x, y = y))
   if (point) {
     p <- 
diff --git a/R/print-spillover.R b/R/print-spillover.R
index 7410d8e1..8412ccec 100644
--- a/R/print-spillover.R
+++ b/R/print-spillover.R
@@ -97,19 +97,19 @@ print.bvhardynsp <- function(x, digits = max(3L, getOption("digits") - 3L), ...)
   if (is_mcmc) {
     dim_data <- nrow(x$to) / length(x$index)
     to_distn <-
-      x$to %>%
-      select("series", "mean") %>%
-      mutate(id = rep(x$index, each = dim_data)) %>%
+      x$to |>
+      select("series", "mean") |>
+      mutate(id = rep(x$index, each = dim_data)) |>
       pivot_wider(names_from = "series", values_from = "mean")
     from_distn <-
-      x$from %>%
-      select("series", "mean") %>%
-      mutate(id = rep(x$index, each = dim_data)) %>%
+      x$from |>
+      select("series", "mean") |>
+      mutate(id = rep(x$index, each = dim_data)) |>
       pivot_wider(names_from = "series", values_from = "mean")
     net_distn <-
-      x$net %>%
-      select("series", "mean") %>%
-      mutate(id = rep(x$index, each = dim_data)) %>%
+      x$net |>
+      select("series", "mean") |>
+      mutate(id = rep(x$index, each = dim_data)) |>
       pivot_wider(names_from = "series", values_from = "mean")
   } else {
     to_distn <- x$to
diff --git a/R/print-varlse.R b/R/print-varlse.R
index eaed13df..10f37d53 100644
--- a/R/print-varlse.R
+++ b/R/print-varlse.R
@@ -93,8 +93,8 @@ print.summary.varlse <- function(x, digits = max(3L, getOption("digits") - 3L),
   dim_data <- ncol(x$covmat)
   dim_design <- nrow(coef_mat) / dim_data
   coef_mat <- 
-    coef_mat %>% 
-    separate(term, into = c("term", "variable"), sep = "\\.") %>% 
+    coef_mat |> 
+    separate(term, into = c("term", "variable"), sep = "\\.") |> 
     split.data.frame(f = gl(dim_data, dim_design))
   if (signif_code) {
     sig_star <- numeric(dim_design)
diff --git a/R/print-vharlse.R b/R/print-vharlse.R
index ffe74f86..8c14ab2d 100644
--- a/R/print-vharlse.R
+++ b/R/print-vharlse.R
@@ -13,11 +13,11 @@ print.vharlse <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
   # phihat_mat <- switch(
   #   x$type,
   #   "const" = {
-  #     split.data.frame(x$coefficients[-(3 * x$m + 1),], gl(3, x$m)) %>% 
+  #     split.data.frame(x$coefficients[-(3 * x$m + 1),], gl(3, x$m)) |> 
   #       lapply(t)
   #   },
   #   "none" = {
-  #     split.data.frame(x$coefficients, gl(3, x$m)) %>% 
+  #     split.data.frame(x$coefficients, gl(3, x$m)) |> 
   #       lapply(t)
   #   }
   # )
@@ -99,8 +99,8 @@ print.summary.vharlse <- function(x, digits = max(3L, getOption("digits") - 3L),
   dim_data <- ncol(x$covmat)
   dim_design <- nrow(coef_mat) / dim_data
   coef_mat <- 
-    coef_mat %>% 
-    separate(term, into = c("term", "variable"), sep = "\\.") %>% 
+    coef_mat |> 
+    separate(term, into = c("term", "variable"), sep = "\\.") |> 
     split.data.frame(f = gl(dim_data, dim_design))
   if (signif_code) {
     sig_star <- numeric(dim_design)
diff --git a/R/spillover.R b/R/spillover.R
index a3874881..62324a09 100644
--- a/R/spillover.R
+++ b/R/spillover.R
@@ -39,9 +39,9 @@ spillover.olsmod <- function(object, n_ahead = 10L, ...) {
   colnames(res$connect) <- colnames(object$coefficients)
   rownames(res$connect) <- colnames(object$coefficients)
   res$df_long <-
-    res$connect %>%
-    as.data.frame() %>%
-    rownames_to_column(var = "series") %>%
+    res$connect |>
+    as.data.frame() |>
+    rownames_to_column(var = "series") |>
     pivot_longer(-"series", names_to = "shock", values_to = "spillover")
   colnames(res$net_pairwise) <- colnames(res$connect)
   rownames(res$net_pairwise) <- rownames(res$connect)
@@ -108,9 +108,9 @@ spillover.normaliw <- function(object, n_ahead = 10L, num_iter = 5000L, num_burn
   colnames(res$connect) <- colnames(object$coefficients)
   rownames(res$connect) <- colnames(object$coefficients)
   res$df_long <-
-    res$connect %>%
-    as.data.frame() %>%
-    rownames_to_column(var = "series") %>%
+    res$connect |>
+    as.data.frame() |>
+    rownames_to_column(var = "series") |>
     pivot_longer(-"series", names_to = "shock", values_to = "spillover")
   colnames(res$net_pairwise) <- colnames(res$connect)
   rownames(res$net_pairwise) <- rownames(res$connect)
@@ -171,7 +171,7 @@ spillover.bvarldlt <- function(object, n_ahead = 10L, level = .05, sparse = FALS
   from_distn <- process_vector_draws(sp_res$from, dim_data = dim_data, level = level, med = FALSE)
   net_distn <- process_vector_draws(sp_res$net, dim_data = dim_data, level = level, med = FALSE)
   df_long <-
-    join_long_spillover(connect_distn, prefix = "spillover") %>%
+    join_long_spillover(connect_distn, prefix = "spillover") |>
     left_join(join_long_spillover(net_pairwise_distn, prefix = "net"), by = c("series", "shock"))
   res <- list(
     connect = connect_distn,
@@ -237,7 +237,7 @@ spillover.bvharldlt <- function(object, n_ahead = 10L, level = .05, sparse = FAL
   from_distn <- process_vector_draws(sp_res$from, dim_data = dim_data, level = level, med = FALSE)
   net_distn <- process_vector_draws(sp_res$net, dim_data = dim_data, level = level, med = FALSE)
   df_long <-
-    join_long_spillover(connect_distn, prefix = "spillover") %>%
+    join_long_spillover(connect_distn, prefix = "spillover") |>
     left_join(join_long_spillover(net_pairwise_distn, prefix = "net"), by = c("series", "shock"))
   res <- list(
     connect = connect_distn,
@@ -383,7 +383,7 @@ dynamic_spillover.normaliw <- function(object, n_ahead = 10L, window,
         # num_chains = num_chains, num_iter = object$iter, num_burn = object$burn, thin = object$thin,
         lag = object$p, bayes_spec = object$spec, include_mean = include_mean,
         seed_chain = sample.int(.Machine$integer.max, size = num_horizon),
-        # seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        # seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         nthreads = num_thread
       )
     },
@@ -394,7 +394,7 @@ dynamic_spillover.normaliw <- function(object, n_ahead = 10L, window,
         # num_chains = num_chains, num_iter = object$iter, num_burn = object$burn, thin = object$thin,
         week = object$week, month = object$month, bayes_spec = object$spec, include_mean = include_mean,
         seed_chain = sample.int(.Machine$integer.max, size = num_horizon),
-        # seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        # seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         nthreads = num_thread
       )
     },
@@ -509,7 +509,7 @@ dynamic_spillover.ldltmod <- function(object, n_ahead = 10L, window, level = .05
         grp_id = grp_id, own_id = own_id, cross_id = cross_id, grp_mat = object$group,
         include_mean = include_mean,
         # seed_chain = sample.int(.Machine$integer.max, size = num_horizon),
-        seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         nthreads = num_thread
       )
     },
@@ -534,7 +534,7 @@ dynamic_spillover.ldltmod <- function(object, n_ahead = 10L, window, level = .05
         ggl = object$ggl,
         grp_id = grp_id, own_id = own_id, cross_id = cross_id, grp_mat = object$group,
         include_mean = include_mean,
-        seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        seed_chain = sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         nthreads = num_thread
       )
     },
@@ -549,11 +549,11 @@ dynamic_spillover.ldltmod <- function(object, n_ahead = 10L, window, level = .05
     lapply(
       sp_list$tot,
       function(x) {
-        process_vector_draws(unlist(x), dim_data = 1, level = level, med = FALSE) %>%
-          do.call(cbind, .)
+        processed <- process_vector_draws(unlist(x), dim_data = 1, level = level, med = FALSE)
+        do.call(cbind, processed)
       }
-    ) %>% 
-    do.call(rbind, .)
+    )
+  tot_distn <- do.call(rbind, tot_distn)
   # sp_list <- lapply(sp_list, function(x) {
   #   if (is.matrix(x)) {
   #     return(apply(x, 1, mean))
@@ -632,11 +632,11 @@ dynamic_spillover.svmod <- function(object, n_ahead = 10L, level = .05, sparse =
     lapply(
       sp_list$tot,
       function(x) {
-        process_vector_draws(unlist(x), dim_data = 1, level = level, med = FALSE) %>%
-          do.call(cbind, .)
+        processed <- process_vector_draws(unlist(x), dim_data = 1, level = level, med = FALSE)
+        do.call(cbind, processed)
       }
-    ) %>%
-    do.call(rbind, .)
+    )
+  tot_distn <- do.call(rbind, tot_distn)
   # colnames(sp_list$to) <- paste(colnames(object$y), "to", sep = "_")
   # colnames(sp_list$from) <- paste(colnames(object$y), "from", sep = "_")
   # colnames(sp_list$to) <- colnames(object$y)
diff --git a/R/stable-process.R b/R/stable-process.R
index 77333837..bee1e32b 100644
--- a/R/stable-process.R
+++ b/R/stable-process.R
@@ -32,10 +32,10 @@ is.stable <- function(x, ...) {
 #' where \eqn{A} is VAR(1) coefficient matrix representation.
 #' @export
 stableroot.varlse <- function(x, ...) {
-  compute_var_stablemat(x) %>% 
-    eigen() %>% 
-    .$values %>% 
-    Mod()
+  eigen_vals <-
+    compute_var_stablemat(x) |>
+    eigen()
+  Mod(eigen_vals$values)
 }
 
 #' @rdname is.stable
@@ -55,10 +55,10 @@ is.stable.varlse <- function(x, ...) {
 #' @rdname stableroot
 #' @export
 stableroot.vharlse <- function(x, ...) {
-  compute_vhar_stablemat(x) %>% 
-    eigen() %>% 
-    .$values %>% 
-    Mod()
+  eigen_vals <-
+    compute_vhar_stablemat(x) |>
+    eigen()
+  Mod(eigen_vals$values)
 }
 
 #' @rdname is.stable
@@ -70,10 +70,10 @@ is.stable.vharlse <- function(x, ...) {
 #' @rdname stableroot
 #' @export
 stableroot.bvarmn <- function(x, ...) {
-  compute_var_stablemat(x) %>% 
-    eigen() %>% 
-    .$values %>% 
-    Mod()
+  eigen_vals <- 
+    compute_var_stablemat(x) |> 
+    eigen()
+  Mod(eigen_vals$values)
 }
 
 #' @rdname is.stable
@@ -85,10 +85,10 @@ is.stable.bvarmn <- function(x, ...) {
 #' @rdname stableroot
 #' @export
 stableroot.bvarflat <- function(x, ...) {
-  compute_var_stablemat(x) %>% 
-    eigen() %>% 
-    .$values %>% 
-    Mod()
+  eigen_vals <- 
+    compute_var_stablemat(x) |> 
+    eigen()
+  Mod(eigen_vals$values)
 }
 
 #' @rdname is.stable
@@ -100,10 +100,10 @@ is.stable.bvarflat <- function(x, ...) {
 #' @rdname stableroot
 #' @export
 stableroot.bvharmn <- function(x, ...) {
-  compute_vhar_stablemat(x) %>% 
-    eigen() %>% 
-    .$values %>% 
-    Mod()
+  eigen_vals <-
+    compute_vhar_stablemat(x) |>
+    eigen()
+  Mod(eigen_vals$values)
 }
 
 #' @rdname is.stable
diff --git a/R/summary-bayes.R b/R/summary-bayes.R
index cfe7988c..4ab700a8 100644
--- a/R/summary-bayes.R
+++ b/R/summary-bayes.R
@@ -95,7 +95,7 @@ summary.normaliw <- function(object, num_chains = 1, num_iter = 1000,  num_burn
   #   iw_scale, # scale of IW
   #   nu, # shape of IW
   #   TRUE
-  # ) %>%
+  # ) |>
   #   simplify2array()
   # # preprocess--------------------------------
   # dim_design <- object$df # k or h = 3m + 1 or 3m
@@ -122,16 +122,16 @@ summary.normaliw <- function(object, num_chains = 1, num_iter = 1000,  num_burn
   # len_res <- length(thin_id)
   mn_name <- ifelse(grepl(pattern = "^BVAR_", object$process), "alpha", "phi")
   # # coef_record <-
-  # #   coef_and_sig$mn %>%
-  # #   t() %>%
-  # #   split.data.frame(gl(num_iter, object$m)) %>%
+  # #   coef_and_sig$mn |>
+  # #   t() |>
+  # #   split.data.frame(gl(num_iter, object$m)) |>
   # #   lapply(function(x) c(t(x)))
   # coef_record <- lapply(coef_and_sig[1,], c)
   # coef_record <- coef_record[thin_id]
   # coef_record <- do.call(rbind, coef_record)
   # colnames(coef_record) <- paste0(mn_name, "[", seq_len(ncol(coef_record)), "]")
   # res$coefficients <-
-  #   colMeans(coef_record) %>%
+  #   colMeans(coef_record) |>
   #   matrix(ncol = object$m)
   # # coef_and_sig$iw <- split_psirecord(t(coef_and_sig$iw), chain = 1, varname = "psi")
   # coef_and_sig$iw <- coef_and_sig[2,]
@@ -148,7 +148,7 @@ summary.normaliw <- function(object, num_chains = 1, num_iter = 1000,  num_burn
   # prec_record <- lapply(coef_and_sig$iw, function(x) chol2inv(chol(x)))
   # res$covmat <- Reduce("+", coef_and_sig$iw) / length(coef_and_sig$iw)
   # res$omega_record <-
-  #   lapply(prec_record, diag) %>%
+  #   lapply(prec_record, diag) |>
   #   do.call(rbind, .)
   # colnames(res$omega_record) <- paste0("omega[", seq_len(ncol(res$omega_record)), "]")
   # res$omega_record <- as_draws_df(res$omega_record)
diff --git a/R/summary-forecast.R b/R/summary-forecast.R
index a69e627b..f55e51b8 100644
--- a/R/summary-forecast.R
+++ b/R/summary-forecast.R
@@ -11,10 +11,10 @@ divide_ts <- function(y, n_ahead) {
   num_ts <- nrow(y)
   fac_train <- rep(1, num_ts - n_ahead)
   fac_test <- rep(2, n_ahead)
-  y %>% 
+  y |> 
     split.data.frame(
       factor(c(fac_train, fac_test))
-    ) %>% 
+    ) |> 
     setNames(c("train", "test"))
 }
 
@@ -160,11 +160,11 @@ forecast_roll.normaliw <- function(object, n_ahead, y_test, num_thread = 1, use_
   #   fit_ls <- lapply(
   #     object$param_names,
   #     function(x) {
-  #       subset_draws(object$param, variable = x) %>%
-  #         as_draws_matrix() %>%
+  #       subset_draws(object$param, variable = x) |>
+  #         as_draws_matrix() |>
   #         split.data.frame(gl(num_chains, nrow(object$param) / num_chains))
   #     }
-  #   ) %>%
+  #   ) |>
   #     setNames(paste(object$param_names, "record", sep = "_"))
   # }
   res_mat <- switch(model_type,
@@ -174,7 +174,7 @@ forecast_roll.normaliw <- function(object, n_ahead, y_test, num_thread = 1, use_
       #   y, object$p, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     },
@@ -184,7 +184,7 @@ forecast_roll.normaliw <- function(object, n_ahead, y_test, num_thread = 1, use_
       #   y, object$p, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec$U, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     },
@@ -194,19 +194,19 @@ forecast_roll.normaliw <- function(object, n_ahead, y_test, num_thread = 1, use_
       #   y, object$week, object$month, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     }
   )
   # num_draw <- nrow(object$a_record) # concatenate multiple chains
   # res_mat <-
-  #   res_mat %>%
+  #   res_mat |>
   #   lapply(function(res) {
-  #     unlist(res) %>%
-  #       array(dim = c(1, object$m, num_draw)) %>%
+  #     unlist(res) |>
+  #       array(dim = c(1, object$m, num_draw)) |>
   #       apply(c(1, 2), mean)
-  #   }) %>%
+  #   }) |>
   #   do.call(rbind, .)
   colnames(res_mat) <- name_var
   res <- list(
@@ -305,7 +305,7 @@ forecast_roll.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, level
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -354,7 +354,7 @@ forecast_roll.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, level
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -366,61 +366,23 @@ forecast_roll.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, level
     lpl_val <- res_mat$lpl
     res_mat$lpl <- NULL
   }
-  if (med) {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), median)
-      }) %>%
-      do.call(rbind, .)
-  } else {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), mean)
-      }) %>%
-      do.call(rbind, .)
-  }
-  colnames(pred_mean) <- name_var
-  est_se <- 
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), sd)
-    }) %>%
-    do.call(rbind, .)
-  colnames(est_se) <- name_var
-  lower_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(lower_quantile) <- name_var
-  upper_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = 1 - level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(upper_quantile) <- name_var
+  y_distn <- process_forecast_draws(
+    draws = res_mat,
+    n_ahead = num_horizon,
+    dim_data = object$m,
+    num_draw = num_draw,
+    var_names = name_var,
+    roll = TRUE,
+    med = med
+  )
   res <- list(
     process = object$process,
-    forecast = pred_mean,
-    se = est_se,
-    lower = lower_quantile,
-    upper = upper_quantile,
-    lower_joint = lower_quantile,
-    upper_joint = upper_quantile,
+    forecast = y_distn$mean,
+    se = y_distn$sd,
+    lower = y_distn$lower,
+    upper = y_distn$upper,
+    lower_joint = y_distn$lower,
+    upper_joint = y_distn$upper,
     eval_id = n_ahead:nrow(y_test),
     y = y
   )
@@ -525,7 +487,7 @@ forecast_roll.svmod <- function(object, n_ahead, y_test, num_thread = 1, level =
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -574,7 +536,7 @@ forecast_roll.svmod <- function(object, n_ahead, y_test, num_thread = 1, level =
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -586,61 +548,23 @@ forecast_roll.svmod <- function(object, n_ahead, y_test, num_thread = 1, level =
     lpl_val <- res_mat$lpl
     res_mat$lpl <- NULL
   }
-  if (med) {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), median)
-      }) %>%
-      do.call(rbind, .)
-  } else {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), mean)
-      }) %>%
-      do.call(rbind, .)
-  }
-  colnames(pred_mean) <- name_var
-  est_se <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), sd)
-    }) %>%
-    do.call(rbind, .)
-  colnames(est_se) <- name_var
-  lower_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(lower_quantile) <- name_var
-  upper_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = 1 - level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(upper_quantile) <- name_var
+  y_distn <- process_forecast_draws(
+    draws = res_mat,
+    n_ahead = num_horizon,
+    dim_data = object$m,
+    num_draw = num_draw,
+    var_names = name_var,
+    roll = TRUE,
+    med = med
+  )
   res <- list(
     process = object$process,
-    forecast = pred_mean,
-    se = est_se,
-    lower = lower_quantile,
-    upper = upper_quantile,
-    lower_joint = lower_quantile,
-    upper_joint = upper_quantile,
+    forecast = y_distn$mean,
+    se = y_distn$sd,
+    lower = y_distn$lower,
+    upper = y_distn$upper,
+    lower_joint = y_distn$lower,
+    upper_joint = y_distn$upper,
     eval_id = n_ahead:nrow(y_test),
     y = y
   )
@@ -768,7 +692,7 @@ forecast_expand.normaliw <- function(object, n_ahead, y_test, num_thread = 1, us
       #   y, object$p, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     },
@@ -778,7 +702,7 @@ forecast_expand.normaliw <- function(object, n_ahead, y_test, num_thread = 1, us
       #   y, object$p, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec$U, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     },
@@ -788,19 +712,19 @@ forecast_expand.normaliw <- function(object, n_ahead, y_test, num_thread = 1, us
       #   y, object$week, object$month, num_chains, object$iter, object$burn, object$thin,
       #   fit_ls,
       #   object$spec, include_mean, n_ahead, y_test,
-      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+      #   sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
       #   num_thread, chunk_size
       # )
     }
   )
   # num_draw <- nrow(object$a_record) # concatenate multiple chains
   # res_mat <-
-  #   res_mat %>%
+  #   res_mat |>
   #   lapply(function(res) {
-  #     unlist(res) %>%
-  #       array(dim = c(1, object$m, num_draw)) %>%
+  #     unlist(res) |>
+  #       array(dim = c(1, object$m, num_draw)) |>
   #       apply(c(1, 2), mean)
-  #   }) %>%
+  #   }) |>
   #   do.call(rbind, .)
   colnames(res_mat) <- name_var
   res <- list(
@@ -900,7 +824,7 @@ forecast_expand.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, lev
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -949,7 +873,7 @@ forecast_expand.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, lev
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -961,70 +885,23 @@ forecast_expand.ldltmod <- function(object, n_ahead, y_test, num_thread = 1, lev
     lpl_val <- res_mat$lpl
     res_mat$lpl <- NULL
   }
-  # res_mat <-
-  #   res_mat %>%
-  #   lapply(function(res) {
-  #     unlist(res) %>%
-  #       array(dim = c(1, object$m, num_draw)) %>%
-  #       apply(c(1, 2), mean)
-  #   }) %>%
-  #   do.call(rbind, .)
-  # colnames(res_mat) <- name_var
-  if (med) {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), median)
-      }) %>%
-      do.call(rbind, .)
-  } else {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), mean)
-      }) %>%
-      do.call(rbind, .)
-  }
-  colnames(pred_mean) <- name_var
-  est_se <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), sd)
-    }) %>%
-    do.call(rbind, .)
-  colnames(est_se) <- name_var
-  lower_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(lower_quantile) <- name_var
-  upper_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = 1 - level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(upper_quantile) <- name_var
+  y_distn <- process_forecast_draws(
+    draws = res_mat,
+    n_ahead = num_horizon,
+    dim_data = object$m,
+    num_draw = num_draw,
+    var_names = name_var,
+    roll = TRUE,
+    med = med
+  )
   res <- list(
     process = object$process,
-    forecast = pred_mean,
-    se = est_se,
-    lower = lower_quantile,
-    upper = upper_quantile,
-    lower_joint = lower_quantile,
-    upper_joint = upper_quantile,
+    forecast = y_distn$mean,
+    se = y_distn$sd,
+    lower = y_distn$lower,
+    upper = y_distn$upper,
+    lower_joint = y_distn$lower,
+    upper_joint = y_distn$upper,
     eval_id = n_ahead:nrow(y_test),
     y = y
   )
@@ -1123,7 +1000,7 @@ forecast_expand.svmod <- function(object, n_ahead, y_test, num_thread = 1, level
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -1172,7 +1049,7 @@ forecast_expand.svmod <- function(object, n_ahead, y_test, num_thread = 1, level
         grp_id, own_id, cross_id, grp_mat,
         include_mean, stable, n_ahead, y_test,
         lpl,
-        sample.int(.Machine$integer.max, size = num_chains * num_horizon) %>% matrix(ncol = num_chains),
+        sample.int(.Machine$integer.max, size = num_chains * num_horizon) |> matrix(ncol = num_chains),
         sample.int(.Machine$integer.max, size = num_chains),
         verbose, num_thread
       )
@@ -1184,70 +1061,23 @@ forecast_expand.svmod <- function(object, n_ahead, y_test, num_thread = 1, level
     lpl_val <- res_mat$lpl
     res_mat$lpl <- NULL
   }
-  # res_mat <-
-  #   res_mat %>%
-  #   lapply(function(res) {
-  #     unlist(res) %>%
-  #       array(dim = c(1, object$m, num_draw)) %>%
-  #       apply(c(1, 2), mean)
-  #   }) %>%
-  #   do.call(rbind, .)
-  # colnames(res_mat) <- name_var
-  if (med) {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), median)
-      }) %>%
-      do.call(rbind, .)
-  } else {
-    pred_mean <-
-      res_mat %>%
-      lapply(function(res) {
-        unlist(res) %>%
-          array(dim = c(num_horizon, object$m, num_draw)) %>%
-          apply(c(1, 2), mean)
-      }) %>%
-      do.call(rbind, .)
-  }
-  colnames(pred_mean) <- name_var
-  est_se <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), sd)
-    }) %>%
-    do.call(rbind, .)
-  colnames(est_se) <- name_var
-  lower_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(lower_quantile) <- name_var
-  upper_quantile <-
-    res_mat %>%
-    lapply(function(res) {
-      unlist(res) %>%
-        array(dim = c(num_horizon, object$m, num_draw)) %>%
-        apply(c(1, 2), quantile, probs = 1 - level / 2)
-    }) %>%
-    do.call(rbind, .)
-  colnames(upper_quantile) <- name_var
+  y_distn <- process_forecast_draws(
+    draws = res_mat,
+    n_ahead = num_horizon,
+    dim_data = object$m,
+    num_draw = num_draw,
+    var_names = name_var,
+    roll = TRUE,
+    med = med
+  )
   res <- list(
     process = object$process,
-    forecast = pred_mean,
-    se = est_se,
-    lower = lower_quantile,
-    upper = upper_quantile,
-    lower_joint = lower_quantile,
-    upper_joint = upper_quantile,
+    forecast = y_distn$mean,
+    se = y_distn$sd,
+    lower = y_distn$lower,
+    upper = y_distn$upper,
+    lower_joint = y_distn$lower,
+    upper_joint = y_distn$upper,
     eval_id = n_ahead:nrow(y_test),
     y = y
   )
@@ -1282,7 +1112,7 @@ mse <- function(x, y, ...) {
 #' @references Hyndman, R. J., & Koehler, A. B. (2006). *Another look at measures of forecast accuracy*. International Journal of Forecasting, 22(4), 679-688.
 #' @export
 mse.predbvhar <- function(x, y, ...) {
-  (y - x$forecast)^2 %>% 
+  (y - x$forecast)^2 |> 
     colMeans()
 }
 
@@ -1293,7 +1123,7 @@ mse.predbvhar <- function(x, y, ...) {
 #' @export
 mse.bvharcv <- function(x, y, ...) {
   y_test <- y[x$eval_id,]
-  (y_test - x$forecast)^2 %>% 
+  (y_test - x$forecast)^2 |> 
     colMeans()
 }
 
@@ -1431,9 +1261,9 @@ mase <- function(x, y, ...) {
 #' @export
 mase.predbvhar <- function(x, y, ...) {
   scaled_err <- 
-    x$y %>% 
-    diff() %>% 
-    abs() %>% 
+    x$y |> 
+    diff() |> 
+    abs() |> 
     colMeans()
   apply(
     100 * (y - x$forecast) / scaled_err, 
@@ -1451,9 +1281,9 @@ mase.predbvhar <- function(x, y, ...) {
 #' @export
 mase.bvharcv <- function(x, y, ...) {
   scaled_err <- 
-    x$y %>% 
-    diff() %>% 
-    abs() %>% 
+    x$y |> 
+    diff() |> 
+    abs() |> 
     colMeans()
   y_test <- y[x$eval_id,]
   apply(
diff --git a/R/summary-sparse.R b/R/summary-sparse.R
index 4c81ec4e..0e8ec0ed 100644
--- a/R/summary-sparse.R
+++ b/R/summary-sparse.R
@@ -7,14 +7,10 @@
 compute_ci <- function(draws, level = .05) {
   # low_lev <- ifelse(correction, level / (2 * object$df * object$m), level / 2)
   low_lev <- level / 2
-  cred_int <-
-    draws %>% 
-    summarise_draws(
-      ~quantile(
-        .,
-        prob = c(low_lev, 1 - low_lev)
-      )
-    )
+  cred_int <- summarise_draws(
+    draws,
+    ~quantile(., prob = c(low_lev, 1 - low_lev))
+  )
   colnames(cred_int) <- c("term", "conf.low", "conf.high")
   cred_int
 }
@@ -286,7 +282,7 @@ confusion <- function(x, y, ...) {
 #' @export
 confusion.summary.bvharsp <- function(x, y, truth_thr = 0, ...) {
   est <- factor(c(x$choose_coef * 1), levels = c(0L, 1L))
-  truth <- ifelse(c(abs(y)) <= truth_thr, 0L, 1L) %>% factor(levels = c(0L, 1L))
+  truth <- ifelse(c(abs(y)) <= truth_thr, 0L, 1L) |> factor(levels = c(0L, 1L))
   table(truth = truth, estimation = est)
 }
 
diff --git a/R/summary-varlse.R b/R/summary-varlse.R
index 36039245..af75a471 100644
--- a/R/summary-varlse.R
+++ b/R/summary-varlse.R
@@ -45,15 +45,15 @@ summary.varlse <- function(object, ...) {
   term_name <- lapply(
     var_name,
     function(x) paste(rownames(coef_mat), x, sep = ".")
-  ) %>% 
+  ) |> 
     unlist()
   var_coef <- 
-    var_coef %>% 
-    as.data.frame() %>% 
+    var_coef |> 
+    as.data.frame() |> 
     add_column(
       term = term_name,
       .before = 1
-    ) %>% 
+    ) |> 
     mutate(p.value = 2 * (1 - pt(abs(statistic), df = var_stat$df)))
   log_lik <- logLik(object)
   res <- list(
diff --git a/R/summary-vharlse.R b/R/summary-vharlse.R
index ecf734e4..795e4847 100644
--- a/R/summary-vharlse.R
+++ b/R/summary-vharlse.R
@@ -50,15 +50,15 @@ summary.vharlse <- function(object, ...) {
   term_name <- lapply(
     vhar_name,
     function(x) paste(rownames(coef_mat), x, sep = ".")
-  ) %>% 
+  ) |> 
     unlist()
   vhar_coef <- 
-    vhar_coef %>% 
-    as.data.frame() %>% 
+    vhar_coef |> 
+    as.data.frame() |> 
     add_column(
       term = term_name,
       .before = 1
-    ) %>% 
+    ) |> 
     mutate(p.value = 2 * (1 - pt(abs(statistic), df = vhar_stat$df)))
   log_lik <- logLik(object)
   res <- list(
diff --git a/R/utils-pipe.R b/R/utils-pipe.R
deleted file mode 100644
index fd0b1d13..00000000
--- a/R/utils-pipe.R
+++ /dev/null
@@ -1,14 +0,0 @@
-#' Pipe operator
-#'
-#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-#'
-#' @name %>%
-#' @rdname pipe
-#' @keywords internal
-#' @export
-#' @importFrom magrittr %>%
-#' @usage lhs \%>\% rhs
-#' @param lhs A value or the magrittr placeholder.
-#' @param rhs A function call using the magrittr semantics.
-#' @return The result of calling `rhs(lhs)`.
-NULL
diff --git a/data-raw/etf_vix.R b/data-raw/etf_vix.R
index c610c70e..c9dcbf70 100644
--- a/data-raw/etf_vix.R
+++ b/data-raw/etf_vix.R
@@ -11,17 +11,17 @@ etf_vix_long <- purrr::map_dfr(
   fredr,
   observation_start = as.Date("2012-01-09"), # after Italian debt crisis
   observation_end = as.Date("2015-06-27") # before Grexit
-) %>% 
+) |> 
   select(date, series_id, value)
 # date, variables-----------------------------
 etf_vix_raw <- 
-  etf_vix_long %>% 
+  etf_vix_long |> 
   tidyr::pivot_wider(names_from = "series_id", values_from = "value")
 # only variables and impute missing-----------
 etf_vix <- 
-  etf_vix_raw %>% 
-  select(-date) %>% 
-  apply(2, imputeTS::na_interpolation) %>% 
+  etf_vix_raw |> 
+  select(-date) |> 
+  apply(2, imputeTS::na_interpolation) |> 
   as_tibble()
 # only date-----------------------------------
 trading_day <- etf_vix_raw$date
diff --git a/man/pipe.Rd b/man/pipe.Rd
deleted file mode 100644
index a648c296..00000000
--- a/man/pipe.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils-pipe.R
-\name{\%>\%}
-\alias{\%>\%}
-\title{Pipe operator}
-\usage{
-lhs \%>\% rhs
-}
-\arguments{
-\item{lhs}{A value or the magrittr placeholder.}
-
-\item{rhs}{A function call using the magrittr semantics.}
-}
-\value{
-The result of calling \code{rhs(lhs)}.
-}
-\description{
-See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-}
-\keyword{internal}
diff --git a/vignettes/bvhar.Rmd b/vignettes/bvhar.Rmd
index 5f286f84..963f73d7 100644
--- a/vignettes/bvhar.Rmd
+++ b/vignettes/bvhar.Rmd
@@ -55,7 +55,7 @@ Since this is just an example, we arbitrarily extract a small number of variable
 ```{r etfdat}
 var_idx <- c("GVZCLS", "OVXCLS", "EVZCLS", "VXFXICLS")
 etf <- 
-  etf_vix %>% 
+  etf_vix |> 
   dplyr::select(dplyr::all_of(var_idx))
 etf
 ```
diff --git a/vignettes/forecasting.Rmd b/vignettes/forecasting.Rmd
index b22374a1..ab80e264 100644
--- a/vignettes/forecasting.Rmd
+++ b/vignettes/forecasting.Rmd
@@ -54,8 +54,8 @@ We use coefficient matrix estimated by VAR(5) in introduction vignette.
 
 ```{r evalcoef, echo=FALSE}
 etf_eval <- 
-  etf_vix %>% 
-  dplyr::select(GVZCLS, OVXCLS, EVZCLS, VXFXICLS) %>% 
+  etf_vix |> 
+  dplyr::select(GVZCLS, OVXCLS, EVZCLS, VXFXICLS) |> 
   divide_ts(20)
 etf_train <- etf_eval$train
 etf_test <- etf_eval$test
@@ -258,9 +258,9 @@ list(
   BVAR = mse_bvar,
   BVHAR1 = mse_bvhar_v1,
   BVHAR2 = mse_bvhar_v2
-) %>% 
-  lapply(mean) %>% 
-  unlist() %>% 
+) |> 
+  lapply(mean) |> 
+  unlist() |> 
   sort()
 ```
 
@@ -273,7 +273,7 @@ list(
   pred_bvar,
   pred_bvhar_v1,
   pred_bvhar_v2
-) %>% 
+) |> 
   gg_loss(y = y_test, "mse")
 ```
 
@@ -286,8 +286,8 @@ list(
   BVAR = pred_bvar,
   BVHAR1 = pred_bvhar_v1,
   BVHAR2 = pred_bvhar_v2
-) %>% 
-  lapply(rmape, pred_bench = pred_var, y = y_test) %>% 
+) |> 
+  lapply(rmape, pred_bench = pred_var, y = y_test) |> 
   unlist()
 ```
 
@@ -345,8 +345,8 @@ list(
   BVAR = bvar_roll,
   BVHAR1 = bvhar_roll_v1,
   BVHAR2 = bvhar_roll_v2
-) %>% 
-  lapply(rmape, pred_bench = var_roll, y = y_test) %>% 
+) |> 
+  lapply(rmape, pred_bench = var_roll, y = y_test) |> 
   unlist()
 ```
 
@@ -383,8 +383,8 @@ list(
   BVAR = bvar_expand,
   BVHAR1 = bvhar_expand_v1,
   BVHAR2 = bvhar_expand_v2
-) %>% 
-  lapply(rmape, pred_bench = var_expand, y = y_test) %>% 
+) |> 
+  lapply(rmape, pred_bench = var_expand, y = y_test) |> 
   unlist()
 ```
 

From eaa04d4ad3af0a1f1dfbbe1adb119a46e3959e0c Mon Sep 17 00:00:00 2001
From: Young Geun Kim <ygeunkimstat@gmail.com>
Date: Mon, 23 Dec 2024 19:57:15 +0900
Subject: [PATCH 2/2] Increment version number to 2.1.2.9015

---
 DESCRIPTION | 2 +-
 NEWS.md     | 2 ++
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 4786fd3e..e7a0092e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: bvhar
 Type: Package
 Title: Bayesian Vector Heterogeneous Autoregressive Modeling
-Version: 2.1.2.9014
+Version: 2.1.2.9015
 Authors@R: 
     c(person(given = "Young Geun",
              family = "Kim",
diff --git a/NEWS.md b/NEWS.md
index c316f3bd..87fc3b23 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # bvhar (development version)
 
+* Requires `R >= 4.1` following [tidyverse R version support schedule](https://www.tidyverse.org/blog/2019/04/r-version-support/)
+
 * Use `spdlog` (using `RcppSpdlog`) logger instead of custom progress bar (`bvharprogress`).
 
 * Use `RcppThread` to make the logger thread-safe ([eddelbuettel/rcppspdlog#22](https://github.com/eddelbuettel/rcppspdlog/issues/22))