Robust Factor Analysis Parameter Estimation
Rui Zhou Junyan Liu Sandeep Kumar Daniel P. Palomar


Robust Factor Analysis Parameter EstimationThis work was supported by the Hong Kong RGC 16208917 research grant.
    Peter M. Oppeneer
    December 30, 2023
================================================================================================================




This paper considers the problem of robustly estimating the parameters
of a heavy-tailed multivariate distribution when the covariance matrix
is known to have the structure of a low-rank matrix plus a diagonal
matrix as considered in factor analysis (FA). By assuming the observed
data to follow the multivariate Student's t distribution, we can
robustly estimate the parameters via maximum likelihood estimation
(MLE). However, the MLE of parameters becomes an intractable problem
when the multivariate Student's t distribution and the FA structure
are both introduced. In this paper, we propose an algorithm based
on the generalized expectation maximization (GEM) method to obtain
estimators. The robustness of our proposed method is further enhanced
to cope with missing values. Finally, we show the performance of our
proposed algorithm using both synthetic data and real financial data.







§ INTRODUCTION




Factor analysis (FA) is of great significance in various fields like
finance, statistics, and cognitive ratio <cit.>.
A basic FA model can be written as 𝐱=μ+𝐁𝐟+ε,
where 𝐱∈ℝ^p is the observed vector, μ∈ℝ^p
is a constant vector, 𝐁∈ℝ^p× r(r≪ p)
is the factors loading matrix, 𝐟∈ℝ^r is a
vector of low-dimensional common factors, and ε∈ℝ^p
is a vector of uncorrelated noise. For example, in a financial market,
𝐱 can be the return of stocks, and 𝐟 can be
some macroeconomic factors like growth rate of the GDP, inflation
rate, unemployment rate, etc. <cit.>.  FA model typically
assumes that 𝐟 and ε are uncorrelated
and both zero-mean, and the covariance matrix of 𝐟 is
an r× r identity matrix, denoted by 𝐈_r. Following
this, the covariance matrix of 𝐱 can be expressed as Σ=𝐁𝐁^T+Ψ,
where Ψ is a p× p diagonal matrix containing
the variance of noise at its diagonal. Note that, with the FA structure,
the number of parameters of the covariance matrix has been greatly
reduced from p(p+1)/2 to p(r+1). Therefore,
the estimation of Σ could be improved due to the
FA structure.





 §.§ Related Works






  §.§.§ Learn From Σ:


A large amount of literature has focused on estimating the covariance
matrix with FA structure. One choice is decomposing the estimated
Σ into the sum of a low-rank matrix and a diagonal
matrix, e.g., constrained minimum trace factor analysis (MTFA) <cit.>.
But a main drawback is imposing the hard equality constraint Σ=𝐁𝐁^T+Ψ
, which does not allow any differences between Σ
and 𝐁𝐁^T+Ψ. Its application is
restricted as the exact Σ is usually not observed.
Another choice is to approximate the target matrix by a FA structured
one. A naive method is to obtain 𝐁̂ firstly via
principal component analysis (PCA) and then Ψ̂
by taking directly the residual's sample variance. A joint estimation
method over 𝐁 and Ψ is usually chosen to
minimize ‖Σ-(𝐁𝐁^T+Ψ)‖_F^2,
where ‖·‖_F denotes the Frobenius norm of a matrix.
This problem can be solved by applying PCA iteratively
<cit.>.







  §.§.§ Learn From Data:


Different from works mentioned above, the MLE for FA directly learns
the parameters from raw data. It assumes that the data are generated
from a certain statistical model, typically the multivariate Gaussian
distribution, and then the parameters are estimated by maximizing
the likelihood function. However, a disadvantage of the estimators
under the Gaussian assumption is sensitiveness to outliers <cit.>.
A popular way to obtain a more robust estimation result is to consider
some robust heavy-tailed distribution, such as multivariate Student's
t or multivariate Skew t <cit.>
instead of Gaussian. The two aforementioned methods both assume that
𝐟 and ε follow the same heavy
tail distribution sharing the same degrees of freedom. As 𝐟
and ε are not observed, such an assumption
is very restrictive and difficult to verify in practice.





 §.§ Contributions




 This paper considers more general and practically verifiable assumptions
on the FA model: we only assume that the observation 𝐱
follows a multivariate Student's t distribution and the FA structure
with no additional restrictions on 𝐟 and ε.
For this more general model, we propose an efficient algorithm to
estimate the parameters based on the generalized expectation maximization
(GEM) <cit.> method.   In addition, we use the
PX-EM <cit.> method to accelerate the GEM. Our proposed
algorithm can be easily extended to other situations, e.g., when observed
data contains missing values <cit.> or when it follows
the multivariate Skew t distribution.  With synthetic data, our
proposed algorithm shows great estimation accuracy and robustness
to outliers and missing data, which is very meaningful in practical
applications. We also consider real market data in the numerical results,
where the global minimum variance portfolio is designed using our
estimator and compared with those using other estimators.







§ GAUSSIAN FA PROBLEM




Given the sample covariance matrix of the observed data as 𝐒,
the MLE problem for FA under the Gaussian distribution assumption
(GFA) is given as below:





    Σ,𝐁,Ψ𝗆𝖺𝗑𝗂𝗆𝗂𝗓𝖾       log|Σ^-1|-Tr(Σ^-1𝐒)
    𝗌𝗎𝖻𝗃𝖾𝖼𝗍  𝗍𝗈     Σ=𝐁𝐁^T+Ψ
         Ψ=Diag(ψ_1,…,ψ_p)≻0.

The solution to problem (<ref>) would
be Σ^⋆=𝐒 if the structure constraints
were ignored, but becomes intractable when the FA structure is imposed.
Here we introduce two algorithms for solving problem (<ref>).





  §.§.§ Alternating Algorithm:


Problem (<ref>) can be solved by an
alternating optimization approach, which can be performed by alternately
updating 𝐁 and Ψ (note that this actually
corresponds to an alternating optimization over a factor model decomposition
of Σ^-1 instead of an alternating optimization
 over 𝐁 and Ψ, cf. <cit.>).
For fixed Ψ, the optimal update for 𝐁 is
given next. 




 (<cit.>)
Given a feasible Ψ, the optimal 𝐁^⋆
maximizing problem (<ref>) is 𝐁^⋆=Ψ^1/2𝐔𝐃^1/2
where 𝐔Λ𝐔^T is the eigenvalue
decomposition (EVD) of Ψ^-1/2𝐒Ψ^-1/2
and 𝐃=Diag(d_1,…,d_r,0,…,0)
with d_i=max(λ_i-1,0).




For fixed 𝐁, the update for Ψ
is set as Diag(𝐒-𝐁𝐁^T).






  §.§.§ MM Algorithm:


Recently, a majorization-minimization (MM) <cit.>
based method has been proposed in <cit.> to
obtain the optimal Σ. Plugging the optimal 𝐁^⋆
from Lemma <ref> in (<ref>),
we can achieve a concentrated version of (<ref>).


 


 (<cit.>) Denoting
Φ=Ψ^-1, the problem (<ref>)
is equivalent to minimizing f(ϕ), where
f(ϕ)=f_1(ϕ)-f_2(ϕ)
with f_1(ϕ)=∑_i=1^p(-logϕ_i+S_iiϕ_i),
f_2(ϕ)=-∑_i=1^r(log(max{ 1,λ_i^*})-max{ 1,λ_i^*} +1)
and {λ_i^*} _i=1^r are the top r
eigenvalues of 𝐒^*=Φ^1/2𝐒Φ^1/2.
Besides, f_1(ϕ) and f_2(ϕ)
are both convex in ϕ.



By linearizing the f_2(ϕ)
using its sub-gradient, we can majorize f(ϕ)
by f̅(ϕ)=∑_i=1^p(-logϕ_i+S_iiϕ_i-▿_i^(k)ϕ_i),
where ▿^(k) is a subgradient
of f_2(ϕ) at the kth iteration.
The ▿^(k) can be calculated
as ▿_i^(k)=(Φ^-1/2𝐔^*𝐃_1𝐔^*TΦ^1/2𝐒)_ii
where 𝐔^*Λ^*𝐔^*T is the
EVD of Φ^1/2𝐒^(k)Φ^1/2
and 𝐃_1=Diag(d_1,…,d_r,0,…0)
with d_i=max{ 0,1-1/λ_i^*} <cit.>.
Then the update of Φ can be easily obtained as ϕ_i^(k+1)=(S_ii^(k)-▿_i^(k))^-1
for i=1,…,p. By iteratively taking the above procedure, we
can get a converged sequence of ϕ^(k).
Finally we can set Ψ=Φ^-1 and compute
𝐁 via Lemma <ref>. 

It should be noted that the two algorithms are the same and can be
verified by matrix algebra.





§ PROBLEM STATEMENT




The p-dimensional multivariate Student's t distribution, denoted
as t_p(μ,Σ,ν),
has the probability density function (pdf)



    f(𝐱|θ)=Γ(ν+p/2)/Γ(ν/2)ν^p/2π^p/2|Σ|^1/2[1+1/ν(𝐱-μ)^TΣ^-1(𝐱-μ)]^-ν+p/2

where θ=(μ,Σ,ν),
ν is the degrees of freedom, Σ is the scale
p× p positive definite matrix, μ is the
p-dimensional mean vector, and Γ(a)=∫_0^∞t^(a-1)exp(-t)dt
is the gamma function. Note that the covariance matrix of 𝐱
is ν/ν-2Σ, and it is not defined for
ν≤2. Interestingly, the above multivariate Student's t distribution
can be represented in a hierarchical structure as 𝐱|τi.i.d∼𝒩_p(μ,1/τΣ)
with τi.i.d∼Gamma(ν/2,ν/2),
where 𝒩_p(μ,Σ)
is the multivariate Gaussian distribution with mean vector μ
and covariance matrix Σ. Gamma(a,b)
means gamma distribution with shape a and rate b, whose pdf
is f(τ)=b^aτ^(a-1)exp(-bτ)/Γ(a).

We consider that the observed p-dimensional data 𝐱_t,t=1,…,T
follows the independent and identical distributed (i.i.d.) multivariate
Student's t distribution, i.e., 𝐱_t∼t_p(μ,Σ,ν).
Besides, we assume that 𝐱_t follows the FA model, which
means Σ=𝐁𝐁^T+Ψ. Note
that we omit a scaling factor in order to simplify the notation (recall
that here Σ refers to the scale matrix).  A natural
approach is to obtain the parameter estimation through MLE method,
i.e., maximizing L(θ|𝐗)=∑_tlog f(𝐱_t|θ)
w.r.t. θ, where 𝐗∈ℝ^T× p
with 𝐱_t along the t-th row.





§ PROBLEM SOLUTION




It is very difficult to directly solve the above MLE problem as the
objective function and constraints are both non-convex. The expectation
maximization (EM) algorithm is a powerful iterative method to handle
such problem <cit.>. By incorporating the latent
data 𝐙, EM can be employed to convert the maximization
for L(𝐗|θ) to the maximization
for a sequence of simpler and solvable problems. In each iteration,
it requires Q(θ|θ^(k)),
which is the expected log-likelihood function of L(𝐗|θ)
with respect to the current conditional distribution of 𝐙
given the 𝐗 and the current estimate of the parameter
θ^(k). Then it finds θ^(k+1)
by maximizing Q(θ|θ^(k)).
However, the computational cost of solving the subproblem might still
be rather heavy and make the whole EM algorithm impractical. That
is when the GEM algorithm can help. The GEM is an iterative method
based on the EM philosophy but requiring an improvement at each iteration
instead of a full maximization as in EM.





 §.§ The RFA Algorithm




In this section, we propose a robust factor analysis (RFA) algorithm
to solve the above problem. By incorporating the latent variables
τ from the multivariate Student's t hierarchical structure,
the log-likelihood function of the complete data (𝐗,τ)
is given in (<ref>). Note that τ
corresponds to 𝐙 in our application.



    L(θ|𝐱,τ)     =T/2log|Σ^-1|-1/2Tr(Σ^-1∑_t=1^Tτ_t(𝐱_t-μ)(𝐱_t-μ)^T)
         +Tν/2logν/2+ν/2∑_t=1^T(logτ_t-τ_t)-TlogΓ(ν/2)+const.







  §.§.§ Expectation Step:


The expectation step of the GEM algorithm is to find the conditional
expectation of L(θ|𝐗,τ)
over τ given the observed 𝐗 and the
current estimation of θ, i.e., θ^(k).
Since the conditional expectation of τ_t and logτ_t
for t=1,…,T can be directly calculated as



    e_1,t^(k)    =E(τ_t|𝐱_t,θ^(k))=ν^(k)+p/ν^(k)+d(𝐱_t,μ^(k),Σ^(k))
    
    e_2,t^(k)    =E(logτ_t|𝐱_t,θ^(k))=ψ(ν^(k)+p/2)-logν^(k)+d(𝐱_t,μ^(k),Σ^(k))/2

where d(𝐱,μ,Σ)=(𝐱-μ)^T(Σ)^-1(𝐱-μ)
is the Mahalanobis distance between 𝐱_t and μ
<cit.>, then the expectation of the complete data log-likelihood
(<ref>) is 



    Q(θ|θ^(k))     =T/2log|Σ^-1|-1/2Tr(Σ^-1(∑_t=1^Te_1,t^(k)(𝐱_t-μ)(𝐱_t-μ)^T))
         +Tν/2logν/2+ν/2∑_t=1^T(e_2,t^(k)-e_1,t^(k))-TlogΓ(ν/2)+const.







  §.§.§ Maximization Step:


Here we devide the parameters update into two parts.





  
Update of μ and ν:

 It is easy to see that μ and ν are actually
decoupled in Q(θ|θ^(k))
and thus can be easily obtained by setting their derivative to zero.
Then the update scheme for (μ,ν) is



    μ^(k+1)=∑_t=1^Te_1,t^(k)𝐱_t/∑_t=1^Te_1,t^(k)





    ν^(k+1)=ν>0argmax{Tν/2logν/2+ν/2∑_t=1^T(e_2,t^(k)-e_1,t^(k))-TlogΓ(ν/2)} .

According to Proposition 1 in <cit.>, ν^(k+1)
always exists and can be found by bisection search. Interestingly,
the update for μ and ν are independent from
each other and irrelevant to Σ. 





  
Update of Σ:

Fixing μ and ν, maximizing Q(θ|θ^(k))
w.r.t. Σ is reduced to a GFA problem with 𝐒^(k)=1/T∑_t=1^Te_1,t^(k)(𝐱_t-μ^(k+1))(𝐱_t-μ^(k+1))^T,
which can be solved by any of the two iterative methods described
in Section <ref>, which require several
iterations until convergence. Considering that solving Σ
exactly would require several iterations and could be time-consuming,
we can instead only run the algorithms for one round, which would
correspond to implementing the GEM instead of EM.





 §.§ An Acceleration Scheme: PX-EM 




A drawback of the EM algorithm is the slow convergence. The parameter
expanded EM (PX-EM) <cit.> was proposed as an efficient
method to accelerate the classical EM method and can be applied here.
A well-known application of PX-EM on multivariate Student's t case
is to assume that we have τi.i.d∼αGamma(ν/2,ν/2)
in the Student's t hierarchical structure, where α is the
expanded parameter. 





  §.§.§ PX-E


step: The PX-E step needs only a few modifications on the original
expectation step. Specifically, the Mahalanobis distance should be
calculated as d(𝐱_t,μ_*,Σ_*/α),
where (x)_* is the corresponding notation of x
in the parameter expanded statistical model.





  §.§.§ PX-M


step: The update schemes for parameters is similar to those in EM
with only few changes. The update of μ and ν
keeps unchanged but the update of Σ depends on (𝐒^(k))_*=α^(k)𝐒^(k).
The update of the new parameter α is α^(k+1)=α^(k)/T∑_t=1^Te_1,t^(k).
After the algorithm achieves convergence, the real parameters should
be recovered as μ=μ_*, Σ=Σ_*/α
and ν=ν_*.





 §.§ Robust Enhancement to Missing Data 




Due to measurement problems or transmission/storage errors, the observed
data 𝐱_t might contain some missing values, which has
been well studied, cf. <cit.>. It turns out that the missing
values can be regarded as latent data like τ. The
new Q(θ|θ^(k))
at expectation step has the same expression w.r.t. θ
as in (<ref>) <cit.>. Therefore the
maximization step can also be achieved by first updating μ
and ν, and then Σ with FA structure imposed.







§ NUMERICAL EXPERIMENTS






 §.§ Synthetic Data




We generate synthetic data following the multivariate Student's t
distribution. The basic dimension setting is p=100 and r=5.
The true distribution parameters are chosen as follows: ν_true
is set to be 7, the elements of μ_true
are drawn i.i.d. from 𝒩(0,1), 𝐁_true
comes from BARRA industry model <cit.> with each
sector has size 20, diagonal elements of Ψ_true
are generated independently from an exponential distribution with
mean 10. The initial μ^(0), ν^(0)
are set as the sample mean and 10, and Ψ^(0),
𝐁^(0) are given by the naive PCA method.
The convergence condition is set as | L^(k+1)-L^(k)|≤10^-6| L^(k)|.





  §.§.§ Convergence Illustration:


In Fig. <ref>, we compare the convergence
of our proposed methods. The log-likelihood function of the observed
data increases monotonically with the iterations and can finally converge.
The PX-EM can significantly accelerate the convergence in this case.
Then, in Fig. <ref>, we show the
parameter estimation convergence. The parameter estimation normalized
errors are defined as NE(μ)=‖μ̂-μ_true‖_2/‖μ_true‖_2,
NE(Σ)=‖Σ̂-Σ_true‖_F/‖Σ_true‖_F
and NE(ν)=| s(ν̂)-s(ν_true)|/| s(ν_true)|
for μ, Σ and ν, where s(ν)=ν/ν-2.
In Figure <ref>, we can find all
the errors are decreasing and finally converge.





  §.§.§ Robustness Illustration:


We compare our proposed RFA methods in covariance matrix estimation
with sample covariance matrix (SCM), Student's t (Stu-t) estimation
(without FA structure) <cit.>, GFA, and iterative PCA (Iter-PCA)
estimation. The results shown below are averaged over 100 different
realizations of 𝐗 following Gaussian distribution. In
Fig. <ref>, we change the sample number
n but fixing p=100. All methods show better performance when
n goes large while our proposed RFA method always gives the best
result. In Fig. <ref>, we randomly
pick some rows of 𝐗 and element-wisely add outliers drawn
from 𝒩(0,50). It is significant that our proposed
RFA method can still hold a good estimation while the results from
non-robust estimation methods are totally destroyed. The results owe
to the robustness of Student's t assumption in resisting the outliers.
In Fig. <ref>, we randomly pick
some rows of 𝐗 and randomly set 10% values be missing
for each row. Our proposed robust FA algorithm will be directly fed
with incomplete data while for other methods we need to manually remove
the rows containing the missing values. It is impressive that our
proposed robust FA method can keep the lowest and almost unchanged
performance.











 §.§ Real Data




In this Section, we show the performance of our proposed algorithm
based on the backtest with real financial data. We randomly choose
50 stocks for 10 times from Standard & Poor's 500 list and
2 years (2×252) continuous historical daily prices data between
01 Dec 2008 and 01 Dec 2018.  Then for each resampling dataset, we
perform the rolling window backtest with lookback window length set
to 100 days and test window be 5 days. The rebalance is assumed
to be done everyday without transaction cost. To fairly compare the
estimation performance, we are particularly interested in the global
minimum variance portfolio (GMVP): minimize 𝐰^TΣ𝐰
with constraint 1^T𝐰=1, where Σ
is the covariance matrix obtained from various methods. We respectively
set r=2 and r=4 in Fig. <ref>
and <ref>. It turns
out that our proposed RFA method can achieve smaller out-of-sample
risk with less uncertainty.








§ CONCLUSION




In this paper, we have proposed the RFA algorithm to obtain the MLE
of the multivariate Student's t distribution with FA structure
imposed. The algorithm was based on the EM framework and had great
estimation accuracy and robustness to outliers and missing values.
The backtest over real financial data has shown advantages and practical
usefulness of our proposed RFA algorithm.

splncs04


