'Compute Monte Carlo integrated error bands assuming normally
distributed for impulse response functions
Parameters
orth: bool, default False
Compute orthoganalized impulse response error bands
repl: int
number of Monte Carlo replications to perform
T: int, default 10
number of impulse response periods
signif: float (0 < signif <1)
Significance level for error bars, defaults to 95% CI
seed: int
np.random.seed for replications
burn: int
number of initial observations to discard for simulation
cum: bool, default False
produce cumulative irf error bands
Notes
Lutkepohl (2005) Appendix D
Returns
Tuple of lower and upper arrays of ma_rep monte carlo standard errors'

def irf_errband_mc(self, orth=False, repl=1000, T=10, signif=0.05, seed=None, burn=100, cum=False):

    neqs = self.neqs
    k_ar = self.k_ar
    coefs = self.coefs
    sigma_u = self.sigma_u
    intercept = self.intercept
    nobs = self.nobs
    ma_coll = np.zeros((repl, (T + 1), neqs, neqs))
    def fill_coll(sim):
        ret = VAR(sim, exog=self.exog).fit(maxlags=k_ar, trend=self.trend)
        ret = (ret.orth_ma_rep(maxn=T) if orth else ret.ma_rep(maxn=T))
        return (ret.cumsum(axis=0) if cum else ret)
    for i in range(repl):
        sim = util.varsim(coefs, intercept, sigma_u, seed=seed, steps=(nobs + burn))
        sim = sim[burn:]
        ma_coll[i, :, :, :] = fill_coll(sim)
    ma_sort = np.sort(ma_coll, axis=0)
    index = ((round(((signif / 2) * repl)) - 1), (round(((1 - (signif / 2)) * repl)) - 1))
    lower = ma_sort[index[0], :, :, :]
    upper = ma_sort[index[1], :, :, :]
    return (lower, upper)
