docs: def irf(self, ts_length=100, shock=None):
          """
          Create Impulse Response Functions
  
          Parameters
          ----------
  
          ts_length : scalar(int)
              Number of periods to calculate IRF
  
          Shock : array_like(float)
              Vector of shocks to calculate IRF to. Default is first element of w
  
          """
  
          if type(shock) != np.ndarray:
              # Default is to select first element of w
              shock = np.vstack((np.ones((1, 1)), np.zeros((self.nw - 1, 1))))
  
          self.c_irf = np.empty((ts_length, self.nc))
          self.s_irf = np.empty((ts_length, self.nb))
          self.i_irf = np.empty((ts_length, self.ni))
          self.k_irf = np.empty((ts_length, self.nk))
          self.h_irf = np.empty((ts_length, self.nh))
          self.g_irf = np.empty((ts_length, self.ng))
          self.d_irf = np.empty((ts_length, self.nd))
          self.b_irf = np.empty((ts_length, self.nb))
  
          for i in range(ts_length):
              self.c_irf[i, :] = self.Sc.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.s_irf[i, :] = self.Ss.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.i_irf[i, :] = self.Si.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.k_irf[i, :] = self.Sk.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.h_irf[i, :] = self.Sh.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.g_irf[i, :] = self.Sg.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.d_irf[i, :] = self.Sd.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
              self.b_irf[i, :] = self.Sb.dot(
                  np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
  
          return