## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4
)

## ----simulate-----------------------------------------------------------------
library(pvarife)

set.seed(1)
sim <- sim_pvarife(
  n_units   = 50,   # I
  n_time    = 30,   # T
  n_vars    = 2,    # K
  n_lags    = 1,    # lag order
  n_factors = 1,    # number of interactive fixed effects
  seed      = 42
)

dim(sim$y)          # I x T x K
sim$beta_true       # true coefficient vector
sim$sigma_true      # true reduced-form covariance

## ----estimate-----------------------------------------------------------------
fit <- pvarife(
  y         = sim$y,
  n_lags    = 1,
  n_factors = 1,
  n_out     = 20,   # outer GLS iterations (paper default: 50)
  n_in      = 5    # inner PCA/EM iterations (paper default: 10)
)

print(fit)

## ----coef---------------------------------------------------------------------
coef(fit)

## ----irf-point----------------------------------------------------------------
ir <- compute_irf(
  fit,
  n_periods = 8,
  shock     = 1L     # shock to first variable (Cholesky ordering)
)

dim(ir)    # K x n_periods

plot(ir, var_names = c("Variable 1", "Variable 2"))

## ----irf-bands, eval=FALSE----------------------------------------------------
# bands <- irf_bands(
#   fit,
#   n_periods = 8,
#   n_draw    = 200,    # paper default: 500
#   level     = 0.95,
#   seed      = 1
# )
# 
# plot(bands, var_names = c("Variable 1", "Variable 2"))

## ----bootstrap, eval=FALSE----------------------------------------------------
# bsb <- bootstrap_irf_bands(
#   fit,
#   n_periods = 8,
#   n_boot    = 100,   # computationally expensive
#   level     = 0.95,
#   seed      = 2
# )
# plot(bsb, var_names = c("Variable 1", "Variable 2"))

## ----lr-sim, eval=FALSE-------------------------------------------------------
# sim_lr <- sim_pvarife(
#   n_units        = 50,
#   n_time         = 30,
#   identification = "long_run",   # Blanchard-Quah DGP
#   seed           = 42
# )
# sim_lr$identification          # "long_run"
# sim_lr$diff_vars_suggested     # 1L  (display cumulative IRF for variable 1)

## ----lr-irf, eval=FALSE-------------------------------------------------------
# fit_lr <- pvarife(sim_lr$y, n_lags = 1, n_factors = 1, n_out = 20, n_in = 5)
# 
# # Point estimate at estimated beta (uncorrected)
# ir_lr <- compute_irf(
#   fit_lr,
#   n_periods      = 8,
#   identification = "long_run",
#   diff_vars      = sim_lr$diff_vars_suggested   # cumulate variable 1
# )
# plot(ir_lr, var_names = c("Variable 1", "Variable 2"))
# 
# # Bias-corrected point estimate
# ir_lr_bc <- compute_irf(
#   fit_lr,
#   n_periods      = 8,
#   identification = "long_run",
#   diff_vars      = 1L,
#   bias_correct   = TRUE    # apply Theorem 2.3 bias correction
# )
# 
# # Full confidence bands (median is already bias-corrected)
# bands_lr <- irf_bands(
#   fit_lr,
#   n_periods      = 8,
#   identification = "long_run",
#   diff_vars      = 1L,
#   n_draw         = 200,
#   seed           = 1
# )
# plot(bands_lr, var_names = c("Variable 1", "Variable 2"))

## ----bias-correct, eval=FALSE-------------------------------------------------
# ir_unc <- compute_irf(fit, n_periods = 8)
# ir_bc  <- compute_irf(fit, n_periods = 8, bias_correct = TRUE)
# # ir_bc is close to the median returned by irf_bands()

## ----avar, eval=FALSE---------------------------------------------------------
# avar <- asymptotic_var(fit)
# se   <- sqrt(diag(avar$variance))
# 
# # Bias-corrected 95% confidence intervals for each element of beta
# beta_hat   <- as.numeric(fit$beta)
# beta_biasC <- beta_hat - avar$bias
# ci_lower   <- beta_biasC - 1.96 * se
# ci_upper   <- beta_biasC + 1.96 * se
# 
# data.frame(
#   parameter = names(coef(fit)),
#   estimate  = round(beta_hat, 4),
#   std_err   = round(se, 4),
#   ci_lower  = round(ci_lower, 4),
#   ci_upper  = round(ci_upper, 4)
# )

## ----unbalanced, eval=FALSE---------------------------------------------------
# y_unbal <- sim$y
# # Randomly set 15% of unit-time observations to NA
# set.seed(10)
# mask <- array(runif(prod(dim(y_unbal))) < 0.15, dim = dim(y_unbal))
# y_unbal[mask] <- NA
# 
# fit_unbal <- pvarife(y_unbal, n_lags = 1, n_factors = 1, n_out = 10, n_in = 5)
# print(fit_unbal)