SHARE
TWEET

Untitled

a guest Jun 27th, 2019 64 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. r"""
  2.    Markov switching vector autoregressive model
  3.  
  4.    Parameters
  5.    ----------
  6.    endo: pandas DataFrame
  7.        The endogenous variables.
  8.    order: integer
  9.        The order of MSVAR. Default to auto-selection according to information criteria.
  10.    regime: integer
  11.        The number of regimes. Default is 2.
  12.    switch_omega: boolean
  13.        Whether or not there is regime-specific heteroskedasticity, i.e.
  14.        whether or not the error term has a switching variance. Default is
  15.        False.
  16.    """
  17.  
  18. def __init__(self, endo, regime=2, order=None, switch_omega=False):
  19.     self.endo = endo; self.regime = regime; self.order = order
  20.     self.vars = len(self.endo.T)
  21.     self._base_model = VAR(self.endo)
  22.     lags_to_select = np.array(list(self._base_model.select_order().selected_orders.values()))
  23.     mod = mode(lags_to_select); min = lags_to_select.min()
  24.     if self.order is None:
  25.         self.order = mod[0][0] if mod[0][0] < 12 * (len(self.endo)/100)**0.25 and mod[1][0] > 1 else min
  26.     self.switch_omega = switch_omega
  27.  
  28.     # Sanity checks
  29.     if not isinstance(self.endo, pd.core.frame.DataFrame):
  30.         raise TypeError('Use pandas DataFrame with rows as observations and columns as variables.')
  31.     if self.regime < 2:
  32.         raise ValueError('Must have at least 2 regimes.')
  33.     if len(self.endo) <= self.order:
  34.         raise ValueError('Must have more observations than the order of the MSVAR.')
  35.     if not isinstance(self.switch_omega, bool):
  36.         raise ValueError('Switching options must be boolean.')
  37.  
  38. def _initial_values(self):
  39.     # Initial regime transition probability
  40.     P = np.random.rand(self.regime, self.regime)
  41.     P_cols = np.array([sum(P[:, j]) for j in range(len(P))])
  42.     for j, k in zip(range(len(P)), P_cols):
  43.         for i in range(len(P)):
  44.             P[i, j] = P[i, j] / k
  45.     P[-1] = 1 - sum(P[:-1])
  46.     '''you could check the initP for if it's good or bad, maybe could find something.'''
  47.     self.initP = P
  48.     # Initial unconditional regime probability
  49.     rho = np.ones(self.regime) / self.regime
  50.     rho[-1] = 1 - sum(rho[:-1])
  51.     # Initial parameter values
  52.     omega = np.tile(self.endo.cov().values.flatten(), self.regime).reshape(self.regime, self.vars, self.vars)
  53.     res = self._base_model.fit(self.order)
  54.     params = np.tile(res.params.T.values.flatten(), self.regime)
  55.         .reshape(self.regime, self.vars, self.order * self.vars + 1)
  56.     self.__P = P; self.__omega = omega; self.__params = params; self.__rho = rho
  57.     return self.__P, self.__omega, self.__params, self.__rho
  58.  
  59. def _filtered_prob(self):
  60.     filtered = []; filtered1 = []; likelihood = []
  61.     self.__boundary = (1e-06, 1-1e-06)
  62.     lagged = lagmat(self.endo, self.order, use_pandas=True)
  63.     lagged.insert(0, 'intercept', 1)
  64.     self.__lagged = lagged
  65.     for i in range(len(self.endo)):
  66.         if i < self.order:
  67.             # TODO use unconditional
  68.             kstt = self.__rho; kst1t = self.__rho
  69.             filtered.append(kstt); filtered1.append(kst1t); # likelihood.append(sum(dens))
  70.         else:
  71.             dens = np.array([density.pdf(self.endo.iloc[i, :], self.__params[j].dot
  72.             (np.array(lagged.iloc[i, :].T)), self.__omega[j]) for j in range(self.regime)])
  73.             kstt = (kst1t * dens) / sum(kst1t * dens)
  74.             for i in range(len(kstt)):
  75.                 if kstt[i] < min(self.__boundary):
  76.                     kstt[i] = min(self.__boundary)
  77.                 elif kstt[i] > max(self.__boundary):
  78.                     kstt[i] = max(self.__boundary)
  79.             kst1t = self.__P.dot(kstt)
  80.             for i in range(len(kst1t)):
  81.                 if kst1t[i] < min(self.__boundary):
  82.                     kst1t[i] = min(self.__boundary)
  83.                 elif kst1t[i] > max(self.__boundary):
  84.                     kst1t[i] = max(self.__boundary)
  85.             filtered.append(kstt); filtered1.append(kst1t); likelihood.append(np.sum(kst1t * dens));
  86.     self.__filtered = np.array(filtered); self.__filtered1 = np.array(filtered1)
  87.     self.__log_likelihood = np.sum(np.array([np.log(i) for i in likelihood]))
  88.     print('The log-likelihood is {}.'.format(self.__log_likelihood))
  89.     assert len(self.__filtered[self.__filtered >= 1]) is 0
  90.     assert len(self.__filtered[self.__filtered <= 0]) is 0
  91.     assert len(self.__filtered1[self.__filtered1 >= 1]) is 0
  92.     assert len(self.__filtered1[self.__filtered1 <= 0]) is 0
  93.     return self.__filtered, self.__filtered1, self.__log_likelihood, self.__lagged
  94.  
  95. def _smoothed_prob(self):
  96.     smoothed = []
  97.     smoothed.append(self.__filtered[::-1][0])
  98.     for i, j in zip(self.__filtered[::-1][1:], self.__filtered1[::-1][1:]):
  99.         kstT = i * (self.__P.T.dot(smoothed[-1] / j))
  100.         for i in range(len(kstT)):
  101.             if kstT[i] < min(self.__boundary):
  102.                 kstT[i] = min(self.__boundary)
  103.             elif kstT[i] > max(self.__boundary):
  104.                 kstT[i] = max(self.__boundary)
  105.         smoothed.append(kstT)
  106.     self.__smoothed = np.array(smoothed[::-1])
  107.     assert len(self.__smoothed[self.__smoothed >= 1]) is 0
  108.     assert len(self.__smoothed[self.__smoothed <= 0]) is 0
  109.     return self.__smoothed
  110.  
  111. def _em(self):
  112.     # Estimate regime transition probability
  113.     new_P = self.__P.copy()
  114.     for i in range(self.regime):
  115.         for j in range(self.regime):
  116.             new_P[i, j] = np.sum(self.__filtered[self.order:][:, i] * self.__P[i, j] *
  117.                                  self.__smoothed[self.order:][:, j] / self.__filtered1[self.order:][:, j]) /
  118.                           np.sum(self.__smoothed[self.order:][:, i])
  119.             if new_P[i, j] > max(self.__boundary):
  120.                 new_P[i, j] = max(self.__boundary)
  121.             elif new_P[i, j] < min(self.__boundary):
  122.                 new_P[i, j] = min(self.__boundary)
  123.     '''I think this is where the problem at. If I remove the add up to 1 constraint, it yields different results
  124.     at least in the first several recursion. This implies somehow, the 1 - distorts the converging path.
  125.     I don't know if this formula should be theoretically in (0, 1) which can be proven
  126.     analytically with every element of it is in (0, 1). '''
  127.     new_P[-1] = 1 - sum(new_P[:-1])
  128.     print(new_P)
  129.     # assert len(new_P[new_P >= 1]) is 0
  130.     # assert len(new_P[new_P <= 0]) is 0
  131.     # Estimate unconditional regime probability
  132.     new_rho = self.__rho.copy()
  133.     for i in range(self.regime):
  134.         '''According to Hamilton, we should update it with the self.orderth smooth prob (which
  135.         in python is self.order-1 cuz python counts from 0). If the plausible P and params I get is somehow resonable,
  136.         then this one is completely rubbish with value like near 0 for one state and near 1 for another.'''
  137.         new_rho[i] = self.__smoothed[0, i]
  138.         assert len(new_rho[new_rho >= 1]) is 0
  139.         assert len(new_rho[new_rho <= 0]) is 0
  140.         for i in range(len(new_rho)):
  141.             if new_rho[i] < min(self.__boundary):
  142.                 new_rho[i] = min(self.__boundary)
  143.             elif new_rho[i] > max(self.__boundary):
  144.                 new_rho[i] = max(self.__boundary)
  145.     print(new_rho)
  146.     # Estimate parameters
  147.     '''These are estimated using the weighted OLS which maximizes the prob w.r.t smoothed prob (the expectation of EM
  148.     ). Hamilton gives a scalar valued example in his book. I just changed it into the vector version. His 1990 paper gives
  149.     vector version examples but with no autoregressive dynamics. Anyway, it's just weighted OLS estimator.'''
  150.     new_params = self.__params.copy()
  151.     for i in range(self.regime):
  152.         new_params[i] = np.array(
  153.             sum([(np.mat(self.endo.iloc[j, :].values).T * self.__smoothed[j, i]**0.5) *
  154.                     (np.mat(self.__lagged.iloc[j, :].values) * self.__smoothed[j, i]**0.5)
  155.                 for j in range(len(self.endo))]) * np.linalg.inv(
  156.             sum([(np.mat(self.__lagged.iloc[j, :].values).T * self.__smoothed[j, i]**0.5) *
  157.                  (np.mat(self.__lagged.iloc[j, :].values) * self.__smoothed[j, i]**0.5)
  158.                 for j in range(len(self.endo))]))
  159.                                      )
  160.     residuals = np.array([(self.endo.T.values - new_params[i].dot(self.__lagged.T.values)).T
  161.                           for i in range(self.regime)])
  162.     new_omega = self.__omega.copy()
  163.     if self.switch_omega:
  164.         for i in range(self.regime):
  165.             new_omega[i] = sum([np.outer(residuals[i][j], residuals[i][j]) * self.__smoothed[j, i]
  166.                                for j in range(len(self.endo))]) /
  167.                            sum(self.__smoothed[:, i])
  168.     else:
  169.         for i in range(self.regime):
  170.             new_omega[i] = sum([np.outer(residuals[i][j], residuals[i][j]) for j in range(len(self.endo))]) /
  171.                            len(self.endo) - self.order
  172.     self.__new_P = new_P; self.__new_rho = new_rho; self.__new_params = new_params; self.__new_omega = new_omega
  173.     self.__residuals = residuals
  174.     # assert sum(new_P[:, 0]) == 1
  175.     # assert sum(new_P[:, 1]) == 1
  176.     return self.__new_P, self.__new_rho, self.__new_params, self.__new_omega, self.__residuals
  177.  
  178. def fit(self):
  179.     self._initial_values()
  180.     self._filtered_prob()
  181.     self._smoothed_prob()
  182.     self._em()
  183.     convergence_criteria = 1e-06
  184.     P_conditions = np.array([abs(self.__new_P - self.__P) > convergence_criteria]).flatten().all()
  185.     rho_conditions = np.array([abs(self.__new_rho - self.__rho) > convergence_criteria]).flatten().all()
  186.     params_conditions = np.array([abs(self.__new_params - self.__params) > convergence_criteria]).flatten().all()
  187.     omega_conditions = np.array([abs(self.__new_omega - self.__omega) > convergence_criteria]).flatten().all()
  188.     conditions = P_conditions or rho_conditions or params_conditions or omega_conditions
  189.     while conditions:
  190.         '''When the conditions are printed, it turns out the P never satisfies the criteria. This is to say even
  191.         the local maximum is achieved, the P do not converge. I think that's something
  192.         related to the 1- constraint. By the way, you could get the mle result if you set the convergence criteria to
  193.         self.__log_likelihood < 'some specific value depends on your data', many of the times you could end the iteration.'''
  194.         print(P_conditions, rho_conditions, params_conditions, omega_conditions)
  195.         self.__P = self.__new_P; self.__rho = self.__new_rho
  196.         self.__params = self.__new_params; self.__omega = self.__new_omega
  197.         self._filtered_prob()
  198.         self._smoothed_prob()
  199.         self._em()
  200.     else:
  201.         self.P = self.__new_P; self.rho = self.__new_rho
  202.         self.params = self.__new_params; self.omega = self.__new_omega; self.residuals = self.__residuals
  203.         self.filtered, non1, self.log_likelihood, non2 = self._filtered_prob()
  204.         self.smoothed = self._smoothed_prob()
  205.     return MSVARResult
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top