mavats.factormodel

 1from typing import Tuple, Union
 2
 3import numpy as np
 4from numpy import linalg as LA
 5from scipy.linalg import eigh
 6
 7
 8def estimate_factor_model(
 9    X: np.ndarray, h0: int, k1: Union[int, None] = None, k2: Union[int, None] = None
10) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
11    r"""
12    Estimates the high-dimensional matrix factor model in Wang, Liu, Chen 2019 (https://doi.org/10.1016/j.jeconom.2018.09.013)
13    where $X_t = R F_t C^T + E_t = Q_1 Z_t Q_2^T + E_t$
14
15    Parameters
16    ----------
17    X : (T, p1, p2) ndarray
18        The observed data.
19    h0 : int
20        Lag parameter.
21    k1 : int, optional
22        Rank of the front loading matrix. If None, the rank is estimated.
23    k2 : int, optional
24        Rank of the back loading matrix. If None, the rank is estimated.
25
26    Returns
27    -------
28    Z : (T, k1, k2) ndarray
29        The estimated transformed factor matrix.
30    S : (T, p1, p2) ndarray
31        The estimated dynamic signal $S_t = Q_1 Z_t Q_2^T$.
32    Q1 : (p1, k1) ndarray
33        The estimated front loading matrix.
34    Q2 : (p2, k2) ndarray
35        The estimated back loading matrix.
36
37    """
38    M1 = _compute_M(X, h0)
39    Q1 = _compute_Q(M1, k1)
40    M2 = _compute_M(X.transpose(0, 2, 1), h0)
41    Q2 = _compute_Q(M2, k2)
42    Z = Q1.T @ X @ Q2
43    S = Q1 @ Z @ Q2.T
44    return Z, S, Q1, Q2
45
46
47def _compute_omega_hat(X: np.ndarray, h: int) -> np.ndarray:
48    T, p1, p2 = X.shape
49    omega_hat = np.zeros((p2, p2, p1, p1))
50    X_t = X[: T - h]
51    X_th = X[h:T]
52    for i in range(p2):
53        for j in range(p2):
54            omega_hat[i, j] = np.einsum("tk,tl->kl", X_t[:, :, i], X_th[:, :, j])
55    return omega_hat / (T - h)
56
57
58def _compute_M(X: np.ndarray, h0: int) -> np.ndarray:
59    T, p1, p2 = X.shape
60    M = np.zeros((p1, p1))
61    for h in range(1, h0 + 1):
62        omega_hat = _compute_omega_hat(X, h)
63        for i in range(p2):
64            for j in range(p2):
65                M += omega_hat[i, j] @ omega_hat[i, j].T
66    return M
67
68
69def _compute_Q(M: np.ndarray, k: Union[int, None]) -> np.ndarray:
70    if k is None:
71        w, v = eigh(M)
72        k = _estimate_k(w)
73        w, v = w[-k:], v[:, -k:]
74    else:
75        w, v = eigh(M, subset_by_index=(M.shape[0] - k, M.shape[0] - 1))
76    v = _ensure_positive_eigenvecs(v)
77    return v
78
79
80def _estimate_k(w: np.ndarray) -> int:
81    w = np.flipud(w)[: len(w) // 2 + 1]
82    ratios = w[1:] / w[:-1]
83    return int(np.argmin(ratios)) + 1
84
85
86def _ensure_positive_eigenvecs(v: np.ndarray) -> np.ndarray:
87    v_summed = np.sum(v, axis=0)
88    lt_zero = v_summed < 0
89    v[:, lt_zero] *= -1
90    return v