'Simulates impulse response function, returning an array of simulations.
Used for Sims-Zha error band calculation.
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
Sims, Christoper A., and Tao Zha. 1999. "Error Bands for Impulse Response." Econometrica 67: 1113-1155.
Returns
Array of simulated impulse response functions'

def irf_resim(self, orth=False, repl=1000, T=10, 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)
    return ma_coll
