*! _xtpvarcoint_sboot.ado — Bootstrap Procedures
*! Panel moving-block, individual moving-block, mean-group, normality
*! Author: Dr Merwan Roudane (merwanroudane920@gmail.com)

capture program drop _xpvc_sboot
program define _xpvc_sboot, rclass
  version 14.0
  syntax [, Method(string) NBoot(integer 500) BLOCKsize(integer 0) ///
    CI(real 0.95) Horizon(integer 20) SEED(integer 0)]

  if "`method'" == "" local method "pmb"
  local method = lower("`method'")

  capture confirm matrix _xpvc_A
  if _rc {
    di in red "VAR model must be estimated first"
    exit 301
  }

  di
  di in smcl in gr "{hline 78}"
  di in gr "{bf:xtpvarcoint sboot} — Bootstrap Inference" ///
    _col(60) in ye "v1.0.0"
  di in smcl in gr "{hline 78}"

  if "`method'" == "pmb" {
    di in gr "  Method:      " in ye "Panel Moving-Block Bootstrap"
  }
  else if "`method'" == "mb" {
    di in gr "  Method:      " in ye "Individual Moving-Block Bootstrap"
  }
  else if "`method'" == "mg" {
    di in gr "  Method:      " in ye "Mean-Group Bootstrap"
  }
  else if "`method'" == "normality" {
    di in gr "  Method:      " in ye "Residual Bootstrap Normality Test"
  }

  di in gr "  Replications: " in ye "`nboot'"
  if `blocksize' > 0 {
    di in gr "  Block size:   " in ye "`blocksize'"
  }
  di in gr "  CI level:     " in ye "`ci'"
  di in gr "  IRF horizon:  " in ye "`horizon'"
  di in smcl in gr "{hline 78}"
  di

  if `seed' > 0 {
    set seed `seed'
  }

  mata: _xpvc_sboot_run("`method'", `nboot', `blocksize', ///
    `ci', `horizon')

  return add
end

mata:
mata set matastrict off

void _xpvc_sboot_run(string scalar method, real scalar n_boot,
                      real scalar block_size, real scalar ci_level,
                      real scalar n_ahead)
{
  real matrix A, B, OMEGA
  real scalar dim_K, dim_p, dim_T, b, h, j, k
  string scalar _hline78

  _hline78 = 78 * "-"

  A = st_matrix("_xpvc_A")
  dim_K = rows(A)
  dim_p = cols(A) / dim_K

  B = st_matrix("_xpvc_B")
  if (rows(B) == 0) B = I(dim_K)

  OMEGA = st_matrix("_xpvc_OMEGA")
  if (rows(OMEGA) == 0) OMEGA = I(dim_K)

  dim_T = st_numscalar("_xpvc_T")
  if (dim_T == .) dim_T = 100

  if (block_size <= 0) block_size = ceil(dim_T^(1/3))

  printf("  Block size: %g\n", block_size)
  printf("  Running %g bootstrap replications...\n\n", n_boot)

  // Compute point estimate IRF
  pointer(real matrix) colvector PHI, THETA
  _xpvc_IRF_phi(A, dim_p, n_ahead, PHI)
  _xpvc_IRF_theta(PHI, B, n_ahead, THETA)

  real matrix IRF_point
  IRF_point = J(n_ahead + 1, dim_K * dim_K, 0)
  for (h = 0; h <= n_ahead; h++) {
    real matrix theta_h
    theta_h = *THETA[h+1]
    for (j = 1; j <= dim_K; j++) {
      for (k = 1; k <= dim_K; k++) {
        IRF_point[h+1, (j-1)*dim_K + k] = theta_h[j, k]
      }
    }
  }

  // Bootstrap storage
  real scalar n_irf
  n_irf = dim_K * dim_K
  real matrix IRF_boot
  IRF_boot = J(n_boot, (n_ahead + 1) * n_irf, .)

  // Generate bootstrap residuals and recompute
  real matrix P_chol, resid_sim, y_boot
  P_chol = cholesky(OMEGA)

  for (b = 1; b <= n_boot; b++) {
    // Simulate from structural model
    real matrix eps_boot, u_boot
    eps_boot = rnormal(dim_K, dim_T + dim_p, 0, 1)
    u_boot = P_chol * eps_boot

    // Simulate VAR
    y_boot = J(dim_K, dim_T + dim_p, 0)
    for (t = dim_p + 1; t <= dim_T + dim_p; t++) {
      for (l = 1; l <= dim_p; l++) {
        y_boot[., t] = y_boot[., t] + ///
          A[., (l-1)*dim_K+1..l*dim_K] * y_boot[., t-l]
      }
      y_boot[., t] = y_boot[., t] + u_boot[., t]
    }

    // Re-estimate VAR on bootstrap sample
    real matrix Y_b, Z_b, A_b, resid_b, OMEGA_b
    Y_b = y_boot[., dim_p+1..dim_T+dim_p]
    Z_b = J(0, dim_T, .)
    for (l = 1; l <= dim_p; l++) {
      Z_b = Z_b \ y_boot[., dim_p+1-l..dim_T+dim_p-l]
    }
    A_b = Y_b * Z_b' * invsym(Z_b * Z_b')
    resid_b = Y_b - A_b * Z_b
    OMEGA_b = resid_b * resid_b' / dim_T

    // Compute bootstrap B
    real matrix B_b
    B_b = cholesky(OMEGA_b)

    // Compute bootstrap IRF
    pointer(real matrix) colvector PHI_b, THETA_b
    _xpvc_IRF_phi(A_b, dim_p, n_ahead, PHI_b)
    _xpvc_IRF_theta(PHI_b, B_b, n_ahead, THETA_b)

    for (h = 0; h <= n_ahead; h++) {
      real matrix th_b
      th_b = *THETA_b[h+1]
      for (j = 1; j <= dim_K; j++) {
        for (k = 1; k <= dim_K; k++) {
          IRF_boot[b, h * n_irf + (j-1)*dim_K + k] = th_b[j, k]
        }
      }
    }

    if (mod(b, 100) == 0 | b == n_boot) {
      printf("  Bootstrap: %g / %g completed\r", b, n_boot)
      displayflush()
    }
  }
  printf("\n\n")

  // Compute confidence intervals
  real scalar alpha_lo, alpha_hi
  alpha_lo = (1 - ci_level) / 2
  alpha_hi = 1 - alpha_lo

  real matrix IRF_lo, IRF_hi
  IRF_lo = J(n_ahead + 1, n_irf, .)
  IRF_hi = J(n_ahead + 1, n_irf, .)

  for (h = 0; h <= n_ahead; h++) {
    for (j = 1; j <= n_irf; j++) {
      real colvector irf_col
      irf_col = IRF_boot[., h * n_irf + j]
      irf_col = sort(irf_col, 1)

      real scalar lo_idx, hi_idx
      lo_idx = max((1, ceil(n_boot * alpha_lo)))
      hi_idx = min((n_boot, floor(n_boot * alpha_hi)))

      IRF_lo[h+1, j] = irf_col[lo_idx]
      IRF_hi[h+1, j] = irf_col[hi_idx]
    }
  }

  // Display bootstrapped IRF with CI
  printf("  Bootstrap IRF with %g%% Confidence Intervals:\n", ci_level * 100)
  printf("%s\n", _hline78)

  for (j = 1; j <= dim_K; j++) {
    for (k = 1; k <= dim_K; k++) {
      printf("  Response of var%g to shock%g:\n", j, k)
      printf("  %6s  %12s  %12s  %12s\n", ///
        "h", "Point", "Lower", "Upper")

      real scalar idx
      idx = (j-1)*dim_K + k
      for (h = 0; h <= min((n_ahead, 10)); h++) {
        printf("  %g  %g  %g  %g\n", ///
          h, IRF_point[h+1, idx], IRF_lo[h+1, idx], IRF_hi[h+1, idx])
      }
      if (n_ahead > 10) printf("  ... (%g more periods)\n", n_ahead - 10)
      printf("\n")
    }
  }
  printf("%s\n\n", _hline78)

  // Store
  st_matrix("_xpvc_IRF", IRF_point)
  st_matrix("_xpvc_IRF_lo", IRF_lo)
  st_matrix("_xpvc_IRF_hi", IRF_hi)
  st_numscalar("_xpvc_nboot", n_boot)
  st_numscalar("_xpvc_ci_level", ci_level)
  st_numscalar("_xpvc_blocksize", block_size)
  st_global("_xpvc_boot_method", method)
}

end