Macroeconomic factor model (aka explicit factor model)

In this model, the factors are observed economic/financial time series. The macroeconomic factor model can be estimated through Least-Squares (LS) regression: \[\underset{\boldsymbol{\gamma}_{i}}{\mathsf{minimize}}\quad\Vert\mathbf{x}_{i}-\tilde{\mathbf{F}}\boldsymbol{\gamma}_{i}\Vert^{2}\] where \(\tilde{\mathbf{F}}=\left[\begin{array}{cc} \mathbf{1}_{T} & \mathbf{F}\end{array}\right]\) and \(\boldsymbol{\gamma}_{i}=\left[\begin{array}{c} \alpha_{i}\\ \boldsymbol{\beta}_{i} \end{array}\right]\). The closed-form solution is: \(\hat{\boldsymbol{\gamma}}_{i}=\left(\tilde{\mathbf{F}}^{T}\tilde{\mathbf{F}}\right)^{-1}\tilde{\mathbf{F}}^{T}\mathbf{x}_{i}\). Then simply use the factor model decomposition to get the residual \(\boldsymbol{\epsilon}_{i}=\mathbf{x}_{i}-\tilde{\mathbf{F}}\hat{\boldsymbol{\gamma}}_{i}\).

BARRA Industry factor model (specific case of fundamental factor models)

Normally, fundamental factor model use observable asset specific characteristics (fundamentals) like industry classification, market capitalization, style classification (value, growth), etc., to determine the common risk factors \(\mathbf{F}\). In this function, we only consider one of the cases: BARRA Industry factor model, which assumes that there are \(K\) factors corresponding to \(K\) mutually exclusive industries (aka, sectors). Apart from that, the loadings \(\beta_{i,k}\) are directly defined as \[\beta_{i,k}=\begin{cases} 1 & \textrm{if stock } i \textrm{ is in industry } k\\ 0 & \textrm{otherwise.} \end{cases}\] Using compact combination \(\mathbf{B}=\left[\begin{array}{ccc} \boldsymbol{\beta}_{1} & \cdots & \boldsymbol{\beta}_{N}\end{array}\right]^{T}\), the industry factor model is (note that \(\boldsymbol{\alpha} = \mathbf{0}\)): \[\mathbf{x}_{t} = \mathbf{B} \mathbf{f}_{t} + \boldsymbol{\epsilon}_{t}, \; t=1,\dots,T\] where \(\mathbf{x}_{t} = [x_{1,t},\dots,x_{N,t}]^T\) and \(\mathbf{f}_{t} = [f_{1,t},\dots,f_{K,t}]^T\). Here the LS regression can also be applied to recover the factors (instead of the loadings as before) as \[\underset{\mathbf{f}_{t}}{\mathsf{minimize}}\quad\frac{1}{T}\sum_{t=1}^{T}\Vert\mathbf{x}_{t}-\mathbf{B}\mathbf{f}_{t}\Vert_{2}^{2}\] The solution is \(\hat{\mathbf{f}}_{t}=(\mathbf{B}^{T}\mathbf{B})^{-1}\mathbf{B}^{T}\mathbf{x}_{t}, \; t=1,\dots,T\) and the residual can be simply calculated as \([\hat{\epsilon}_{1,t},\dots,\hat{\epsilon}_{N,t}]^{T}=\mathbf{x}_{t}-\mathbf{B}\hat{\mathbf{f}}_{t}\).

Statistical factor model via PCA (aka implicit factor model)

The statistical factor model via Principal Component Analysis (PCA) assumes that \(\mathbf{f}_{t}\) is an affine transformation of \(\mathbf{x}_{t}\), i.e., \(\mathbf{f}_{t}=\mathbf{d}+\mathbf{C}^{T}\mathbf{x}_{t}\), where \(\mathbf{d}\in\mathbb{R}^{K}\) and \(\mathbf{C}\in\mathbb{R}^{N\times K}\) are parameters to be estimated. Therefore, we can formulate the problem as \[ \underset{\boldsymbol{\alpha},\mathbf{B},\mathbf{C},\mathbf{d}}{\textrm{minimize}}\quad\frac{1}{T}\sum_{t=1}^{T}\Vert\mathbf{x}_{t}-\boldsymbol{\alpha}-\mathbf{B}\left(\mathbf{C}^{T}\mathbf{x}_{t}+\mathbf{d}\right)\Vert_{2}^{2} \] which is proven equivalent to PCA and thus solved by \[ \hat{\boldsymbol{\alpha}}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{x}_{t},\quad\hat{\mathbf{B}}=\hat{\mathbf{C}}=\hat{\boldsymbol{\Gamma}}^{(K)},\quad\hat{\mathbf{d}}=-\hat{\mathbf{C}}^{T}\hat{\boldsymbol{\alpha}} \] where \(\hat{\boldsymbol{\Gamma}}^{(K)} \in \mathbb{R}^{N\times K}\) with \(k\)-th column being the \(k\)-th largest eigenvector of sample covariance matrix \(\hat{\mathbf{S}}\). The desired structure on the residual covariance matrix is then enforced.

A further refinement of the above PCA method is based on an iterative approach where each iteration performs an improved PCA (the method above corresponds to just one iteration) [see [2] for details]:

Algorithm 1
1. Calculate sample covariance matrix \(\hat{\mathbf{S}}\) and eigen-decomposition (EVD): \(\hat{\boldsymbol{\Gamma}}_{1} \hat{\boldsymbol{\Lambda}}_{1} \hat{\boldsymbol{\Gamma}}^{T}_{1}\)
2. Set index \(s=1\)
3. Compute \(\hat{\mathbf{B}}_{(s)} = \hat{\boldsymbol{\Gamma}}_{(s)}^{(K)} \hat{\boldsymbol{\Lambda}}^{(K)\frac{1}{2}}_{(s)}\), \(\hat{\boldsymbol{\Psi}}_{(s)} = \textrm{struct}(\hat{\boldsymbol{\Sigma}} - \hat{\mathbf{B}}_{(s)} \hat{\mathbf{B}}^{T}_{(s)})\), \(\hat{\boldsymbol{\Sigma}}_{(s)} = \hat{\mathbf{B}}_{(s)} \hat{\mathbf{B}}^{T}_{(s)} + \hat{\boldsymbol{\Psi}}_{(s)}\)
4. Update EVD: \(\hat{\boldsymbol{\Sigma}} - \hat{\boldsymbol{\Psi}}_{(s)} = \hat{\boldsymbol{\Gamma}}_{(s+1)} \hat{\boldsymbol{\Lambda}}_{(s+1)} \hat{\boldsymbol{\Gamma}}^{T}_{(s+1)}\) and \(s \gets s+1\)
5. Repeat Steps 3-4 until convergence.
6. Return \((\hat{\mathbf{B}}_{(s)}, \hat{\boldsymbol{\Psi}}_{(s)}, \hat{\boldsymbol{\Sigma}}_{(s)})\)

where \(\textrm{struct}()\) is to impose certain structure on \(\hat{\boldsymbol{\Psi}}_{(s)}\), one typical option is diagonal. After the algorithm is done, we can calculate \(\hat{\boldsymbol{\alpha}} = \frac{1}{T} \sum_{t=1}^{T} \mathbf{x}_{t}\) and build the statistical factor model using the algorithm output: \[ \hat{\mathbf{B}} = \hat{\boldsymbol{\Gamma}}^{(K)} \hat{\boldsymbol{\Lambda}}^{(K)^{\frac{1}{2}}}, \quad \hat{\mathbf{f}}_{t} = \hat{\boldsymbol{\Lambda}}^{(K)^{-\frac{1}{2}}} \hat{\boldsymbol{\Gamma}}^{(K)^{T}} (\mathbf{x}_{t} - \hat{\boldsymbol{\alpha}}), \quad \hat{\boldsymbol{\epsilon}}_{t} = \mathbf{x}_{t} - \hat{\boldsymbol{\alpha}} - \hat{\mathbf{B}} \hat{\mathbf{f}}_{t}\] (Note that the Algorithm 1 is equivalent to previous method when only one iteration is performed.)

Through the factor decomposition obtained by the function factorModel()

The first approach is based on the factor model decomposition. As mentioned above, the factor model can be expressed as: \[\mathbf{x}_{t} = \boldsymbol{\alpha} + \mathbf{B} \mathbf{f}_{t} + \boldsymbol{\epsilon}_{t}, \; t = 1, \dots, T\] Assuming \(\{\mathbf{f}_{t}\}\) and \(\{\boldsymbol{\epsilon}_{t}\}\) are uncorrelated, the covariance matrix \(\boldsymbol{\Sigma}\) can be written as \[\boldsymbol{\Sigma} = \mathbf{B} \boldsymbol{\Sigma}_{\mathbf{f}} \mathbf{B}^{T} + \boldsymbol{\Psi}\] where \(\boldsymbol{\Sigma}_{\mathbf{f}} = \mathsf{Cov}[\mathbf{x}_{t}]\) and \(\boldsymbol{\Psi} = \mathsf{Cov}[\boldsymbol{\epsilon}_{t}]\). Therefore, we can simply use result from function factorModel() to estimate covariance matrix \(\boldsymbol{\Sigma}\) as: \[\hat{\boldsymbol{\Sigma}} = \hat{\mathbf{B}} \hat{\boldsymbol{\Sigma}}_{\mathbf{f}} \hat{\mathbf{B}}^{T} + \hat{\boldsymbol{\Psi}}\] where \(\hat{\boldsymbol{\Sigma}}_{\mathbf{f}}\) and \(\hat{\boldsymbol{\Psi}}\) are the sample covariance matrix of \(\{\mathbf{\mathbf{f}}_{t}\}\) and \(\{\boldsymbol{\epsilon}_{t}\}\). Besides, the \(\boldsymbol{\Psi}\) is expected to follow a special structure, i.e., \[\hat{\boldsymbol{\Sigma}} = \hat{\mathbf{B}} \hat{\boldsymbol{\Sigma}}_{\mathbf{f}} \hat{\mathbf{B}}^{T} + \textrm{struct}\{ \hat{\boldsymbol{\Psi}} \}.\] In the statistical factor model by PCA of function factorModel(), the estimate \(\hat{\boldsymbol{\Sigma}}\) is actually available when building the model. Therefore the algorithm output \(\hat{\boldsymbol{\Sigma}}_{(s)}\) is directly extracted as the covariance matrix estimation.

Through maximum likelihood (ML) estimation under the Gaussian assumption

Another popular statistical factor model covariance matrix estimation is based on maximum likelihood (ML) estimation. Under the Gaussian assumption on the returns, it can be formulated as \[ \begin{aligned}\underset{\mathbf{B},\boldsymbol{\Psi}}{\mathsf{minimize}}\,\,\, & \,\,\,-\log\lvert\boldsymbol{\Sigma}^{-1}\rvert+\textrm{Tr}\left(\boldsymbol{\Sigma}^{-1}\hat{\mathbf{S}}\right)\\ \mathsf{subject\,\,to} & \,\,\,\boldsymbol{\Sigma}=\mathbf{B}\mathbf{B}^{T}+\boldsymbol{\Psi}\\ & \,\,\,\mathbf{B}\in\mathbb{R}^{N\times K}\\ & \,\,\,\boldsymbol{\Psi}=\textrm{diag}\left(\psi_{1},\dots,\psi_{p}\right)\ge\epsilon\mathbf{I} \end{aligned} \] To solve this problem, we implement a very efficient algorithm from [3], which is much faster than the one used in stats::factanal(). We also give the complete procedure of the algorithm for reference, see Algorithm 2.

Algorithm 2
1. Require \(\boldsymbol{\Psi}^{\left(0\right)}\) and set \(\boldsymbol{\Phi}^{\left(0\right)}=\left(\boldsymbol{\Psi}^{\left(0\right)}\right)^{-1}\), \(t=0\)
2. Repeat the following
3. Compute \(\bigtriangledown_{i}=\left(\left(\boldsymbol{\Phi}^{\left(t\right)}\right)^{-\frac{1}{2}}\mathbf{U}\mathbf{D}_{1}\mathbf{U}^{T}\left(\boldsymbol{\Phi}^{\left(t\right)}\right)^{\frac{1}{2}}\mathbf{S}\right)_{ii}\) for \(i=1,\dots,p\), where \(\mathbf{U}\mathrm{diag}\left(\lambda_{1}^{*},\dots,\lambda_{p}^{*}\right)\mathbf{U}^{T}=\left(\boldsymbol{\Phi}^{\left(t\right)}\right)^{\frac{1}{2}}\mathbf{S}\left(\boldsymbol{\Phi}^{\left(t\right)}\right)^{\frac{1}{2}}\) and \(\mathbf{D}_{1}=\mathrm{diag}\left(\delta_{1},\dots,\delta_{p}\right)\) with \(\delta_{i}=\begin{cases} \max\left\{ 0,1-\frac{1}{\lambda_{i}^{*}}\right\} & 1\le i\le r\\ 0 & \textrm{otherwise} \end{cases}\)
4. Update \(\Phi_{ii}^{\left(t+1\right)}=\min\left\{ \frac{1}{S_{ii}-\bigtriangledown_{i}},\frac{1}{\epsilon}\right\}\) for \(i=1,\dots,p\) and \(t\leftarrow t+1\)
5. Until convergence.
6. Recovery \(\boldsymbol{\Psi}^{\star}=\left(\boldsymbol{\Phi}^{\left(t\right)}\right)^{-1}\)
7. Compute \(\mathbf{B}^{\star}=\left(\boldsymbol{\Psi}^{\star}\right)^{\frac{1}{2}}\left[\mathbf{z}_{1},\mathbf{z}_{2},\dots,\mathbf{z}_{r}\right]\) where \(\mathbf{z}_{i}\) are the largest \(r\) eigenvectors of \(\left(\boldsymbol{\Psi}^{\star}\right)^{-\frac{1}{2}}\mathbf{S}\left(\boldsymbol{\Psi}^{\star}\right)^{-\frac{1}{2}}\) rescaled by corresponding largest \(r\) eigenvalues \(\lambda_{i}^{*},i\le r\), with \(\lVert\mathbf{z}_{i}\rVert^{2}=\max\left\{ 1,\lambda_{i}^{*}\right\} -1\)

Group Lasso

Basically, the group lasso problem is \[ \begin{aligned}\underset{\boldsymbol{\beta}}{\textrm{minimize}} & \quad\frac{1}{2}\Bigg\lVert\mathbf{x}-\sum_{l}^{m}\mathbf{F}^{\left(l\right)}\boldsymbol{\beta}^{\left(l\right)}\Bigg\rVert+\lambda\sum_{l}^{m}\sqrt{p_{l}}\big\lVert\boldsymbol{\beta}^{\left(l\right)}\big\rVert_{2}\end{aligned} \] where \(\mathbf{F}^{\left(l\right)}\) is the submatrix of \(\mathbf{F}\) with columns corresponding to the factors in group \(l\), \(\boldsymbol{\beta}^{\left(l\right)}\) the coefficient vector of that group and \(p_{l}\) is the length of \(\boldsymbol{\beta}^{\left(l\right)}\). We can reformulate our macroeconomic factor model estimation problem by using group lasso \[ \begin{aligned}\underset{\textrm{vec}\left(\mathbf{B}^{T}\right)}{\textrm{minimize}} & \quad\frac{1}{2}\Bigg\lVert\textrm{vec}\left(\mathbf{X}\right)-\left(\mathbf{I}\otimes\mathbf{F}\right)\textrm{vec}\left(\mathbf{B}^{T}\right)\Bigg\rVert+\lambda\sum_{l}^{m}\sqrt{p_{l}}\big\lVert\textrm{vec}\left(\mathbf{B}^{T}\right)^{\left(l\right)}\big\rVert_{2}\end{aligned} \] As our expectation, that \(\mathbf{B}\) should be column-wise sparse, we cam simply pass the information that \(i\)-th and \(\left(i+K\right)\)-th factors are of same group. We found a package SGL which fits a linear regression model of lasso and group lasso regression, i.e., \[ \begin{aligned}\underset{\textrm{vec}\left(\mathbf{B}^{T}\right)}{\textrm{minimize}} & \quad\frac{1}{2}\Bigg\lVert\textrm{vec}\left(\mathbf{X}\right)-\left(\mathbf{I}\otimes\mathbf{F}\right)\textrm{vec}\left(\mathbf{B}^{T}\right)\Bigg\rVert+\left(1-\alpha\right)\lambda\sum_{l}^{m}\sqrt{p_{l}}\big\lVert\textrm{vec}\left(\mathbf{B}^{T}\right)^{\left(l\right)}\big\rVert_{2}+\alpha\lambda\big\lVert\textrm{vec}\left(\mathbf{B}^{T}\right)^{\left(l\right)}\big\rVert_{1}\end{aligned} \] where \(\alpha\) is the turning parameter for a convex combination of the lasso and group lasso penalties. In our case, we realize what we want by the following R codes

set.seed(123)
library(SGL)

# implement the function with package SGL
linreg_row_sparse <- function(Factor, X, lambda = 0.01, alpha = 0.85) {
  N <- ncol(X)
  K <- ncol(Factor)
  index <- rep(1:K, N)
  data <- list(x = diag(N) %x% Factor, y = as.vector(X))
  beta <- SGL(data, index, type = "linear", lambdas = lambda / N, alpha = alpha,
              thresh = 1e-5, standardize = FALSE)$beta
  B <- t(matrix(beta, K, N, byrow = FALSE))
  return(B)
}

Then, we generate data with some factors only influencing one among \(N\) stocks, i.e., these factors are low-impact.

n_lowimp <- 4
F_ <- mvrnorm(n = T, mu = rep(0, n_lowimp+K), Sigma = diag(n_lowimp+K))

B_ <- cbind(B_real, matrix(0, N, n_lowimp))
for (i in 1:n_lowimp) {
  B_[i, K+i] <- 0.5 
}
print(B_)
#>             [,1]       [,2]       [,3] [,4] [,5] [,6] [,7]
#>  [1,] -0.7152422  1.3011760  0.0000000  0.5  0.0  0.0  0.0
#>  [2,] -0.7526890  0.7567748  1.2181086  0.0  0.5  0.0  0.0
#>  [3,] -0.9385387 -1.7267304 -1.3387743  0.0  0.0  0.5  0.0
#>  [4,] -1.0525133 -0.6015067  0.6608203  0.0  0.0  0.0  0.5
#>  [5,]  0.0000000  0.0000000 -0.5229124  0.0  0.0  0.0  0.0
#>  [6,]  0.0000000  0.7035239  0.6837455  0.0  0.0  0.0  0.0
#>  [7,] -2.0142105  0.0000000  0.0000000  0.0  0.0  0.0  0.0
#>  [8,]  0.0000000 -1.2586486  0.6329607  0.0  0.0  0.0  0.0
#>  [9,]  1.2366750  1.6844357  1.3355176  0.0  0.0  0.0  0.0
#> [10,]  2.0375740  0.9113913  0.0000000  0.0  0.0  0.0  0.0
X_ <- F_ %*% t(B_) + 0.01 * matrix(rnorm(T*N), T, N)

We then compare the differences between element-wise sparse and row-wise sparse regression.

B_elesparse <- linreg_ele_sparse(F_, X_, penalty = "lasso", lambda = 0.3)$B
B_rowsparse <- linreg_row_sparse(F_, X_, lambda = 0.3, alpha = 0.2)

print(B_)
#>             [,1]       [,2]       [,3] [,4] [,5] [,6] [,7]
#>  [1,] -0.7152422  1.3011760  0.0000000  0.5  0.0  0.0  0.0
#>  [2,] -0.7526890  0.7567748  1.2181086  0.0  0.5  0.0  0.0
#>  [3,] -0.9385387 -1.7267304 -1.3387743  0.0  0.0  0.5  0.0
#>  [4,] -1.0525133 -0.6015067  0.6608203  0.0  0.0  0.0  0.5
#>  [5,]  0.0000000  0.0000000 -0.5229124  0.0  0.0  0.0  0.0
#>  [6,]  0.0000000  0.7035239  0.6837455  0.0  0.0  0.0  0.0
#>  [7,] -2.0142105  0.0000000  0.0000000  0.0  0.0  0.0  0.0
#>  [8,]  0.0000000 -1.2586486  0.6329607  0.0  0.0  0.0  0.0
#>  [9,]  1.2366750  1.6844357  1.3355176  0.0  0.0  0.0  0.0
#> [10,]  2.0375740  0.9113913  0.0000000  0.0  0.0  0.0  0.0
print(B_elesparse)
#>             [,1]       [,2]       [,3]      [,4]      [,5]      [,6]       [,7]
#>  [1,] -0.4672228  1.0303872  0.0000000 0.2074709 0.0000000 0.0000000 0.00000000
#>  [2,] -0.4967683  0.3758835  0.8502721 0.0000000 0.1161491 0.0000000 0.00000000
#>  [3,] -0.5907305 -1.3630490 -1.0550118 0.0000000 0.0000000 0.1904491 0.00000000
#>  [4,] -0.7142322 -0.2500819  0.2563827 0.0000000 0.0000000 0.0000000 0.08978659
#>  [5,]  0.0000000  0.0000000 -0.2029512 0.0000000 0.0000000 0.0000000 0.00000000
#>  [6,]  0.0000000  0.3654966  0.3474925 0.0000000 0.0000000 0.0000000 0.00000000
#>  [7,] -1.7212507  0.0000000  0.0000000 0.0000000 0.0000000 0.0000000 0.00000000
#>  [8,]  0.0000000 -0.9518036  0.3266415 0.0000000 0.0000000 0.0000000 0.00000000
#>  [9,]  0.9110410  1.3026263  1.0140192 0.0000000 0.0000000 0.0000000 0.00000000
#> [10,]  1.7007091  0.5442624  0.0000000 0.0000000 0.0000000 0.0000000 0.00000000
print(B_rowsparse)
#>             [,1]       [,2]        [,3] [,4] [,5] [,6] [,7]
#>  [1,] -0.5672868  0.9521714  0.00164297    0    0    0    0
#>  [2,] -0.5698816  0.4169104  0.70405147    0    0    0    0
#>  [3,] -0.6291273 -1.1877945 -0.84138534    0    0    0    0
#>  [4,] -0.7616124 -0.3750200  0.30986008    0    0    0    0
#>  [5,]  0.0000000  0.0000000 -0.29007523    0    0    0    0
#>  [6,]  0.0000000  0.4614279  0.38507470    0    0    0    0
#>  [7,] -1.5377723  0.0000000  0.00000000    0    0    0    0
#>  [8,]  0.0000000 -0.8912436  0.37153231    0    0    0    0
#>  [9,]  0.8944246  1.1460135  0.80219627    0    0    0    0
#> [10,]  1.5247966  0.5781963  0.00000000    0    0    0    0

Obviously, we can obtain the row-sparse \(\mathbf{B}\) using sparse-group lasso by properly choosing penalty coefficient \(\lambda\) and \(\alpha\).

Subset selection

In machine learning, there exists a classical method called subset selection. We found a function best.r.sq() from a recently released R package mvabund, which implements a forward selection in a multivariate linear model.

library(mvabund)
best.r.sq( X_~F_ )
#> $xs
#> [1] 1 2 3
#> 
#> $r2Step
#>   F_[, 1]  +F_[, 2]  +F_[, 3] 
#> 0.3379401 0.6656976 0.9675958 
#> 
#> $r2Matrix
#>             Step 1    Step 2    Step 3
#> F_[, 1] 0.33794014        NA        NA
#> F_[, 2] 0.31615795 0.6656976        NA
#> F_[, 3] 0.29150681 0.6309482 0.9675958
#> F_[, 4] 0.02937063 0.3650031 0.6868668
#> F_[, 5] 0.01793328 0.3543163 0.6734623
#> F_[, 6] 0.02076671 0.3561918 0.6799328
#> F_[, 7] 0.01915951 0.3556588 0.6843262

Then, we can perform the trivial linear factor model regression with chosen factors.