code: """
  Impulse reponse-related code
  """
  
  from __future__ import division
  
  import numpy as np
  import numpy.linalg as la
  import scipy.linalg as L
  
  from scipy import stats
  
  from statsmodels.tools.decorators import cache_readonly
  from statsmodels.tools.tools import chain_dot
  #from statsmodels.tsa.api import VAR
  
  import statsmodels.tsa.tsatools as tsa
  import statsmodels.tsa.vector_ar.plotting as plotting
  import statsmodels.tsa.vector_ar.util as util
  
  mat = np.array
  
  class BaseIRAnalysis(object):
      """
      Base class for plotting and computing IRF-related statistics, want to be
      able to handle known and estimated processes
      """
  
      def __init__(self, model, P=None, periods=10, order=None, svar=False):
          self.model = model
          self.periods = periods
          self.neqs, self.lags, self.T = model.neqs, model.k_ar, model.nobs
  
          self.order = order
  
          if P is None:
              sigma = model.sigma_u
  
              # TODO, may be difficult at the moment
              # if order is not None:
              #     indexer = [model.get_eq_index(name) for name in order]
              #     sigma = sigma[:, indexer][indexer, :]
  
              #     if sigma.shape != model.sigma_u.shape:
              #         raise ValueError('variable order is wrong length')
  
              P = la.cholesky(sigma)
  
          self.P = P
  
          self.svar = svar
  
          self.irfs = model.ma_rep(periods)
          if svar:
              self.svar_irfs = model.svar_ma_rep(periods, P=P)
          else:
              self.orth_irfs = model.orth_ma_rep(periods)
  
          self.cum_effects = self.irfs.cumsum(axis=0)
          if svar:
              self.svar_cum_effects = self.svar_irfs.cumsum(axis=0)
          else:
              self.orth_cum_effects = self.orth_irfs.cumsum(axis=0)
  
          self.lr_effects = model.long_run_effects()
          if svar:
              self.svar_lr_effects = np.dot(model.long_run_effects(), P)
          else:
              self.orth_lr_effects = np.dot(model.long_run_effects(), P)
  
  
          # auxiliary stuff
          self._A = util.comp_matrix(model.coefs)
  
      def cov(self, *args, **kwargs):
          raise NotImplementedError
  
      def cum_effect_cov(self, *args, **kwargs):
          raise NotImplementedError
  
      def plot(self, orth=False, impulse=None, response=None,
               signif=0.05, plot_params=None, subplot_params=None,
               plot_stderr=True, stderr_type='asym', repl=1000,
               seed=None, component=None):
          """
          Plot impulse responses
  
          Parameters
          ----------
          orth : bool, default False
              Compute orthogonalized impulse responses
          impulse : string or int
              variable providing the impulse
          response : string or int
              variable affected by the impulse
          signif : float (0 < signif < 1)
              Significance level for error bars, defaults to 95% CI
          subplot_params : dict
              To pass to subplot plotting funcions. Example: if fonts are too big,
              pass {'fontsize' : 8} or some number to your taste.
          plot_params : dict
  
          plot_stderr: bool, default True
              Plot standard impulse response error bands
          stderr_type: string
              'asym': default, computes asymptotic standard errors
              'mc': monte carlo standard errors (use rpl)
          repl: int, default 1000
              Number of replications for Monte Carlo and Sims-Zha standard errors
          seed: int
              np.random.seed for Monte Carlo replications
          component: array or vector of principal component indices
          """
          periods = self.periods
          model = self.model
          svar = self.svar
  
          if orth and svar:
              raise ValueError("For SVAR system, set orth=False")
  
          if orth:
              title = 'Impulse responses (orthogonalized)'
              irfs = self.orth_irfs
          elif svar:
              title = 'Impulse responses (structural)'
              irfs = self.svar_irfs
          else:
              title = 'Impulse responses'
              irfs = self.irfs
  
          if plot_stderr == False:
              stderr = None
  
          elif stderr_type not in ['asym', 'mc', 'sz1', 'sz2','sz3']:
              raise ValueError("Error type must be either 'asym', 'mc','sz1','sz2', or 'sz3'")
          else:
              if stderr_type == 'asym':
                  stderr = self.cov(orth=orth)
              if stderr_type == 'mc':
                  stderr = self.errband_mc(orth=orth, svar=svar,
                                           repl=repl, signif=signif,
                                           seed=seed)
              if stderr_type == 'sz1':
                  stderr = self.err_band_sz1(orth=orth, svar=svar,
                                             repl=repl, signif=signif,
                                             seed=seed,
                                             component=component)
              if stderr_type == 'sz2':
                  stderr = self.err_band_sz2(orth=orth, svar=svar,
                                             repl=repl, signif=signif,
                                             seed=seed,
                                             component=component)
              if stderr_type == 'sz3':
                  stderr = self.err_band_sz3(orth=orth, svar=svar,
                                             repl=repl, signif=signif,
                                             seed=seed,
                                             component=component)
  
          plotting.irf_grid_plot(irfs, stderr, impulse, response,
                                 self.model.names, title, signif=signif,
                                 subplot_params=subplot_params,
                                 plot_params=plot_params, stderr_type=stderr_type)
  
      def plot_cum_effects(self, orth=False, impulse=None, response=None,
                           signif=0.05, plot_params=None,
                           subplot_params=None, plot_stderr=True,
                           stderr_type='asym', repl=1000, seed=None):
          """
          Plot cumulative impulse response functions
  
          Parameters
          ----------
          orth : bool, default False
              Compute orthogonalized impulse responses
          impulse : string or int
              variable providing the impulse
          response : string or int
              variable affected by the impulse
          signif : float (0 < signif < 1)
              Significance level for error bars, defaults to 95% CI
          subplot_params : dict
              To pass to subplot plotting funcions. Example: if fonts are too big,
              pass {'fontsize' : 8} or some number to your taste.
          plot_params : dict
  
          plot_stderr: bool, default True
              Plot standard impulse response error bands
          stderr_type: string
              'asym': default, computes asymptotic standard errors
              'mc': monte carlo standard errors (use rpl)
          repl: int, default 1000
              Number of replications for monte carlo standard errors
          seed: int
              np.random.seed for Monte Carlo replications
  
          """
  
          if orth:
              title = 'Cumulative responses responses (orthogonalized)'
              cum_effects = self.orth_cum_effects
              lr_effects = self.orth_lr_effects
          else:
              title = 'Cumulative responses'
              cum_effects = self.cum_effects
              lr_effects = self.lr_effects
  
          if stderr_type not in ['asym', 'mc']:
              raise TypeError
          else:
              if stderr_type == 'asym':
                  stderr = self.cum_effect_cov(orth=orth)
              if stderr_type == 'mc':
                  stderr = self.cum_errband_mc(orth=orth, repl=repl,
                                                  signif=signif, seed=seed)
          if not plot_stderr:
              stderr = None
  
          plotting.irf_grid_plot(cum_effects, stderr, impulse, response,
                                 self.model.names, title, signif=signif,
                                 hlines=lr_effects, subplot_params=subplot_params,
                                 plot_params=plot_params, stderr_type=stderr_type)
  
  class IRAnalysis(BaseIRAnalysis):
      """
      Impulse response analysis class. Computes impulse responses, asymptotic
      standard errors, and produces relevant plots
  
      Parameters
      ----------
      model : VAR instance
  
      Notes
      -----
      Using Lutkepohl (2005) notation
      """
      def __init__(self, model, P=None, periods=10, order=None, svar=False):
          BaseIRAnalysis.__init__(self, model, P=P, periods=periods,
                                  order=order, svar=svar)
  
          self.cov_a = model._cov_alpha
          self.cov_sig = model._cov_sigma
  
          # memoize dict for G matrix function
          self._g_memo = {}
  
      def cov(self, orth=False):
          """
          Compute asymptotic standard errors for impulse response coefficients
  
          Notes
          -----
          Lutkepohl eq 3.7.5
  
          Returns
          -------
          """
          if orth:
              return self._orth_cov()
  
          covs = self._empty_covm(self.periods + 1)
          covs[0] = np.zeros((self.neqs ** 2, self.neqs ** 2))
          for i in range(1, self.periods + 1):
              Gi = self.G[i - 1]
              covs[i] = chain_dot(Gi, self.cov_a, Gi.T)
  
          return covs
  
      def errband_mc(self, orth=False, svar=False, repl=1000,
                     signif=0.05, seed=None, burn=100):
          """
          IRF Monte Carlo integrated error bands
          """
          model = self.model
          periods = self.periods
          if svar == True:
              return model.sirf_errband_mc(orth=orth, repl=repl, T=periods,
                                          signif=signif, seed=seed,
                                          burn=burn, cum=False)
          else:
              return model.irf_errband_mc(orth=orth, repl=repl, T=periods,
                                          signif=signif, seed=seed,
                                          burn=burn, cum=False)
      def err_band_sz1(self, orth=False, svar=False, repl=1000,
                       signif=0.05, seed=None, burn=100, component=None):
          """
          IRF Sims-Zha error band method 1. Assumes symmetric error bands around
          mean.
  
          Parameters
          ----------
          orth : bool, default False
              Compute orthogonalized impulse responses
          repl : int, default 1000
              Number of MC replications
          signif : float (0 < signif < 1)
              Significance level for error bars, defaults to 95% CI
          seed : int, default None
              np.random seed
          burn : int, default 100
              Number of initial simulated obs to discard
          component : neqs x neqs array, default to largest for each
              Index of column of eigenvector/value to use for each error band
              Note: period of impulse (t=0) is not included when computing
              principle component
  
          References
          ----------
          Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
          Response". Econometrica 67: 1113-1155.
          """
  
          model = self.model
          periods = self.periods
          if orth:
              irfs = self.orth_irfs
          elif svar:
              irfs = self.svar_irfs
          else:
              irfs = self.irfs
          neqs = self.neqs
          irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
                                     burn=100)
          q = util.norm_signif_level(signif)
  
          W, eigva, k =self._eigval_decomp_SZ(irf_resim)
  
          if component != None:
              if np.shape(component) != (neqs,neqs):
                  raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs))
              if np.argmax(component) >= neqs*periods:
                  raise ValueError("Atleast one of the components does not exist")
              else:
                  k = component
  
          # here take the kth column of W, which we determine by finding the largest eigenvalue of the covaraince matrix
          lower = np.copy(irfs)
          upper = np.copy(irfs)
          for i in xrange(neqs):
              for j in xrange(neqs):
                  lower[1:,i,j] = irfs[1:,i,j] + W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]])
                  upper[1:,i,j] = irfs[1:,i,j] - W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]])
  
  
          return lower, upper
  
      def err_band_sz2(self, orth=False, repl=1000, signif=0.05,
                       seed=None, burn=100, component=None):
          """
          IRF Sims-Zha error band method 2.
  
          This method Does not assume symmetric error bands around mean.
  
          Parameters
          ----------
          orth : bool, default False
              Compute orthogonalized impulse responses
          repl : int, default 1000
              Number of MC replications
          signif : float (0 < signif < 1)
              Significance level for error bars, defaults to 95% CI
          seed : int, default None
              np.random seed
          burn : int, default 100
              Number of initial simulated obs to discard
          component : neqs x neqs array, default to largest for each
              Index of column of eigenvector/value to use for each error band
              Note: period of impulse (t=0) is not included when computing
              principle component
  
          References
          ----------
          Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
          Response". Econometrica 67: 1113-1155.
          """
          model = self.model
          periods = self.periods
          if orth:
              irfs = self.orth_irfs
          elif svar:
              irfs = self.svar_irfs
          else:
              irfs = self.irfs
          neqs = self.neqs
          irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
                                     burn=100)
  
          W, eigva, k = self._eigval_decomp_SZ(irf_resim)
  
          if component != None:
              if np.shape(component) != (neqs,neqs):
                  raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs))
              if np.argmax(component) >= neqs*periods:
                  raise ValueError("Atleast one of the components does not exist")
              else:
                  k = component
  
          gamma = np.zeros((repl, periods+1, neqs, neqs))
          for p in xrange(repl):
              for i in xrange(neqs):
                  for j in xrange(neqs):
                      gamma[p,1:,i,j] = W[i,j,k[i,j],:] * irf_resim[p,1:,i,j]
  
          gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles
          indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1
  
          lower = np.copy(irfs)
          upper = np.copy(irfs)
          for i in xrange(neqs):
              for j in xrange(neqs):
                  lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j]
                  upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j]
  
          return lower, upper
  
      def err_band_sz3(self, orth=False, repl=1000, signif=0.05,
                       seed=None, burn=100, component=None):
          """
          IRF Sims-Zha error band method 3. Does not assume symmetric error bands around mean.
  
          Parameters
          ----------
          orth : bool, default False
              Compute orthogonalized impulse responses
          repl : int, default 1000
              Number of MC replications
          signif : float (0 < signif < 1)
              Significance level for error bars, defaults to 95% CI
          seed : int, default None
              np.random seed
          burn : int, default 100
              Number of initial simulated obs to discard
          component : vector length neqs, default to largest for each
              Index of column of eigenvector/value to use for each error band
              Note: period of impulse (t=0) is not included when computing
              principle component
  
          References
          ----------
          Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
          Response". Econometrica 67: 1113-1155.
          """
  
          model = self.model
          periods = self.periods
          if orth:
              irfs = self.orth_irfs
          elif svar:
              irfs = self.svar_irfs
          else:
              irfs = self.irfs
          neqs = self.neqs
          irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
                                     burn=100)
          stack = np.zeros((neqs, repl, periods*neqs))
  
          #stack left to right, up and down
  
          for p in xrange(repl):
              for i in xrange(neqs):
                  stack[i, p,:] = np.ravel(irf_resim[p,1:,:,i].T)
  
          stack_cov=np.zeros((neqs, periods*neqs, periods*neqs))
          W = np.zeros((neqs, periods*neqs, periods*neqs))
          eigva = np.zeros((neqs, periods*neqs))
          k = np.zeros((neqs))
  
          if component != None:
              if np.size(component) != (neqs):
                  raise ValueError("Component array must be of length " + str(neqs))
              if np.argmax(component) >= neqs*periods:
                  raise ValueError("Atleast one of the components does not exist")
              else:
                  k = component
  
          #compute for eigen decomp for each stack
          for i in xrange(neqs):
              stack_cov[i] = np.cov(stack[i],rowvar=0)
              W[i], eigva[i], k[i] = util.eigval_decomp(stack_cov[i])
  
          gamma = np.zeros((repl, periods+1, neqs, neqs))
          for p in xrange(repl):
              c=0
              for j in xrange(neqs):
                  for i in xrange(neqs):
                          gamma[p,1:,i,j] = W[j,k[j],i*periods:(i+1)*periods] * irf_resim[p,1:,i,j]
                          if i == neqs-1:
                              gamma[p,1:,i,j] = W[j,k[j],i*periods:] * irf_resim[p,1:,i,j]
  
          gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles
          indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1
  
          lower = np.copy(irfs)
          upper = np.copy(irfs)
          for i in xrange(neqs):
              for j in xrange(neqs):
                  lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j]
                  upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j]
  
          return lower, upper
  
      def _eigval_decomp_SZ(self, irf_resim):
          """
          Returns
          -------
          W: array of eigenvectors
          eigva: list of eigenvalues
          k: matrix indicating column # of largest eigenvalue for each c_i,j
  
          """
          neqs = self.neqs
          periods = self.periods
  
          cov_hold = np.zeros((neqs, neqs, periods, periods))
          for i in xrange(neqs):
              for j in xrange(neqs):
                  cov_hold[i,j,:,:] = np.cov(irf_resim[:,1:,i,j],rowvar=0)
  
          W = np.zeros((neqs, neqs, periods, periods))
          eigva = np.zeros((neqs, neqs, periods, 1))
          k = np.zeros((neqs, neqs))
  
          for i in xrange(neqs):
              for j in xrange(neqs):
                  W[i,j,:,:], eigva[i,j,:,0], k[i,j] = util.eigval_decomp(cov_hold[i,j,:,:])
          return W, eigva, k
  
      @cache_readonly
      def G(self):
          # Gi matrices as defined on p. 111
  
          K = self.neqs
  
          # nlags = self.model.p
          # J = np.hstack((np.eye(K),) + (np.zeros((K, K)),) * (nlags - 1))
  
          def _make_g(i):
              # p. 111 Lutkepohl
              G = 0.
              for m in range(i):
                  # be a bit cute to go faster
                  idx = i - 1 - m
                  if idx in self._g_memo:
                      apow = self._g_memo[idx]
                  else:
                      apow = la.matrix_power(self._A.T, idx)
                      # apow = np.dot(J, apow)
                      apow = apow[:K]
                      self._g_memo[idx] = apow
  
                  # take first K rows
                  piece = np.kron(apow, self.irfs[m])
                  G = G + piece
  
              return G
  
          return [_make_g(i) for i in range(1, self.periods + 1)]
  
      def _orth_cov(self):
          # Lutkepohl 3.7.8
  
          Ik = np.eye(self.neqs)
          PIk = np.kron(self.P.T, Ik)
          H = self.H
  
          covs = self._empty_covm(self.periods + 1)
          for i in range(self.periods + 1):
              if i == 0:
                  apiece = 0
              else:
                  Ci = np.dot(PIk, self.G[i-1])
                  apiece = chain_dot(Ci, self.cov_a, Ci.T)
  
              Cibar = np.dot(np.kron(Ik, self.irfs[i]), H)
              bpiece = chain_dot(Cibar, self.cov_sig, Cibar.T) / self.T
  
              # Lutkepohl typo, cov_sig correct
              covs[i] = apiece + bpiece
  
          return covs
  
      def cum_effect_cov(self, orth=False):
          """
          Compute asymptotic standard errors for cumulative impulse response
          coefficients
  
          Parameters
          ----------
          orth : boolean
  
          Notes
          -----
          eq. 3.7.7 (non-orth), 3.7.10 (orth)
  
          Returns
          -------
  
          """
          Ik = np.eye(self.neqs)
          PIk = np.kron(self.P.T, Ik)
  
          F = 0.
          covs = self._empty_covm(self.periods + 1)
          for i in range(self.periods + 1):
              if i > 0:
                  F = F + self.G[i - 1]
  
              if orth:
                  if i == 0:
                      apiece = 0
                  else:
                      Bn = np.dot(PIk, F)
                      apiece = chain_dot(Bn, self.cov_a, Bn.T)
  
                  Bnbar = np.dot(np.kron(Ik, self.cum_effects[i]), self.H)
                  bpiece = chain_dot(Bnbar, self.cov_sig, Bnbar.T) / self.T
  
                  covs[i] = apiece + bpiece
              else:
                  if i == 0:
                      covs[i] = np.zeros((self.neqs**2, self.neqs**2))
                      continue
  
                  covs[i] = chain_dot(F, self.cov_a, F.T)
  
          return covs
  
      def cum_errband_mc(self, orth=False, repl=1000,
                            signif=0.05, seed=None, burn=100):
          """
          IRF Monte Carlo integrated error bands of cumulative effect
          """
          model = self.model
          periods = self.periods
          return model.irf_errband_mc(orth=orth, repl=repl,
                                      T=periods, signif=signif, seed=seed, burn=burn, cum=True)
  
      def lr_effect_cov(self, orth=False):
          """
          Returns
          -------
  
          """
          lre = self.lr_effects
          Finfty = np.kron(np.tile(lre.T, self.lags), lre)
          Ik = np.eye(self.neqs)
  
          if orth:
              Binf = np.dot(np.kron(self.P.T, np.eye(self.neqs)), Finfty)
              Binfbar = np.dot(np.kron(Ik, lre), self.H)
  
              return (chain_dot(Binf, self.cov_a, Binf.T) +
                      chain_dot(Binfbar, self.cov_sig, Binfbar.T))
          else:
              return chain_dot(Finfty, self.cov_a, Finfty.T)
  
      def stderr(self, orth=False):
          return np.array([tsa.unvec(np.sqrt(np.diag(c)))
                           for c in self.cov(orth=orth)])
  
      def cum_effect_stderr(self, orth=False):
          return np.array([tsa.unvec(np.sqrt(np.diag(c)))
                           for c in self.cum_effect_cov(orth=orth)])
  
      def lr_effect_stderr(self, orth=False):
          cov = self.lr_effect_cov(orth=orth)
          return tsa.unvec(np.sqrt(np.diag(cov)))
  
      def _empty_covm(self, periods):
          return np.zeros((periods, self.neqs ** 2, self.neqs ** 2),
                          dtype=float)
  
      @cache_readonly
      def H(self):
          k = self.neqs
          Lk = tsa.elimination_matrix(k)
          Kkk = tsa.commutation_matrix(k, k)
          Ik = np.eye(k)
  
          # B = chain_dot(Lk, np.eye(k**2) + commutation_matrix(k, k),
          #               np.kron(self.P, np.eye(k)), Lk.T)
  
          # return np.dot(Lk.T, L.inv(B))
  
          B = chain_dot(Lk,
                        np.dot(np.kron(Ik, self.P), Kkk) + np.kron(self.P, Ik),
                        Lk.T)
  
          return np.dot(Lk.T, L.inv(B))
  
      def fevd_table(self):
          pass
  
  
  
apis: numpy.sqrt numpy.eye numpy.tile numpy.linalg.matrix_power numpy.zeros numpy.dot numpy.diag numpy.copy scipy.linalg.inv numpy.linalg.cholesky numpy.argmax numpy.ravel numpy.size numpy.shape numpy.sort numpy.kron numpy.cov