Asymptotic Covariance of Factor Loadings and Uniquenesses in Classical Factor Analysis

1. Introduction

Factor-analytic representations play a central role in econometric modeling by providing a flexible framework for dimension reduction and structured covariance modeling in settings characterized by many correlated variables. In macroeconomics and finance, factor structures underpin dynamic factor models, structural vector autoregressions (SVARs), and variance-covariance decompositions of large systems, where a small number of latent shocks drive comovement among a high-dimensional set of observables. Beyond likelihood-based approaches, a broad class of non-maximum-likelihood (non-ML) factor extraction methods—such as principal factor, least-squares, alpha, and image factor analysis—remain widely used owing to their computational simplicity, robustness, and minimal distributional assumptions.

Let $Y={(X_1,\dots ,X_p)}^{\top }$ denote a p-dimensional random vector with variance–covariance (or correlation) matrix $\mathrm{\Sigma }$. Classical factor analysis postulates the decomposition

$$\mathrm{\Sigma }=\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }+\mathrm{\Psi },$$

(1)

where $\mathrm{\Lambda }$ is a $p\times k$ loading matrix with $k<p$, and $\mathrm{\Psi }$ is a diagonal matrix of idiosyncratic variances (uniquenesses). In econometric applications, such a structure often emerges as a reduced-form representation of systems driven by clustered or common shocks. Prominent examples include structural VARs with latent structural factors, cross-sectional dependence in macroeconomic panels, and error-component models for residual covariance matrices. Consequently, estimation of $\mathrm{\Lambda }$ and $\mathrm{\Psi }$ is not merely of descriptive interest but constitutes a key input for structural inference, including impulse response analysis, forecast-error variance decompositions, and related functionals.

While point estimation of factor loadings is relatively straightforward across a wide range of extraction methods, formal statistical inference remains considerably more limited. In particular, closed-form standard error formulas for factor loadings and uniquenesses are readily available almost exclusively within the framework of maximum-likelihood (ML) factor analysis, where Fisher-information-based results follow directly from likelihood theory. By contrast, for non-ML extraction methods—despite their widespread use in applied econometrics and standard software packages—explicit asymptotic covariance formulas are largely unavailable or are typically obtained only indirectly through numerical differentiation or resampling procedures.

This gap has important implications for econometric practice. Non-ML factor estimators are often preferred in environments characterized by non-normality, heavy tails, or data heterogeneity, where likelihood-based assumptions are implausible. Observed variables may also be binary, categorical, ordinal, bounded, or integer-valued. Moreover, in time-series and panel settings, the sampling variability of the sample covariance or correlation matrix may itself exhibit serial dependence or heteroskedasticity, rendering distribution-free or sandwich-based inference particularly appealing. However, in the absence of closed-form expressions that link the asymptotic covariance of factor estimators to that of the sample covariance matrix, transparent inference becomes difficult to implement and the resulting estimates are challenging to interpret.

The purpose of this paper is to address this limitation by developing explicit, closed-form expressions for the asymptotic covariance matrices of factor loading and uniqueness estimators obtained from several widely used non-ML extraction methods. Specifically, we consider least-squares factor analysis, principal factor analysis, iterative principal component factor analysis, alpha factor analysis, and image factor analysis. Our approach expresses the asymptotic covariance of these estimators directly in terms of the asymptotic covariance of the sample covariance or correlation matrix, $\mathrm{acov}\left(\widehat{\mathrm{\Sigma }}\right)$. As a result, the proposed formulas remain valid under both normal and non-normal sampling schemes and can be readily combined with existing asymptotic results for $\widehat{\mathrm{\Sigma }}$, including heteroskedasticity- and autocorrelation-robust estimators.

Methodologically, our analysis is based on multivariate implicit differentiation of the estimating equations that characterize each factor solution. By viewing the estimated loading matrix as a differentiable function of the sample covariance matrix, we employ the delta method to derive the asymptotic distribution of the estimators. This unified framework yields analytical expressions for the Jacobian matrices that map perturbations in $\mathrm{\Sigma }$ to changes in factor loadings and uniquenesses, resulting in tractable covariance formulas. Importantly, these results also provide the key inputs required for delta-method-based inference on rotated factor solutions, building directly on the classical contributions of Archer and Jennrich (1973), Jennrich (1973), and related studies.

From an econometric perspective, the proposed results offer key several advantages. First, they enable transparent computation of standard errors without reliance on numerical differentiation, simulation-based resampling, or likelihood-specific assumptions. Second, they apply naturally to settings in which factor analysis serves as an intermediate step—such as SVARs with latent structural shocks—thereby allowing uncertainty in factor estimates to propagate coherently into downstream quantities of interest. Third, because the formulas depend only on $\mathrm{acov}\left(\widehat{\mathrm{\Sigma }}\right)$, they accommodate a broad class of data-generating processes, including non-Gaussian and high-frequency environments.

To illustrate the practical implementation of the proposed approach, we present an empirical application based on least-squares factor analysis and conduct a Monte Carlo simulation to assess the finite-sample performance of the derived standard errors for uniqueness estimators. The simulation results indicate that the empirical standard errors converge to their theoretical asymptotic counterparts, providing strong support for the analytical derivations.

The remainder of the paper is organized as follows. Section 2 reviews exploratory factor analysis. Section 3 derives the asymptotic covariance matrix for iterative principal component and principal factor analysis. Section 4 addresses least-squares factor analysis, while Section 5 examines alpha factor analysis and Section 6 studies image factor analysis. Section 7 derives closed-form expressions for the standard error of uniqueness estimators. Section 8 presents an empirical illustration and simulation study, and Section 9 concludes with a discussion of implications for econometric modeling and inference.

2. Structural Modeling by Exploratory Factor Analysis

For a vector of correlated random variables $Y={(X_1,X_2,\dots ,X_p)}^{\top }$, common (or exploratory) factor analysis posits the existence of k mutually independent common factors that underlie the observed variables Y. Let these latent factors be jointly represented by the random vector F. The relationship between the observed variables and the latent factors is specified as

$$Y=c+\mathrm{\Lambda }F+ϵ,$$

(2)

where c is a constant vector and $ϵ$ is a p-dimensional vector of mutually independent idiosyncratic disturbances. Assuming that the common factors have unit variance, Eq. (2) implies that the covariance structure of Y is given by Eq. (1), where the diagonal elements of $\mathrm{\Psi }$ correspond to the variances of the idiosyncratic disturbances, commonly referred to as uniquenesses.

Eq. (1) implies that the off-diagonal elements of the covariance matrix of Y are entirely determined by the factor loadings. In contrast, the diagonal elements reflect contributions from both common and idiosyncratic components: the common factors capture variation shared across the variables in Y, whereas the idiosyncratic disturbances represent variable-specific noise. Let ${\lambda }_{it}$ denote the element in the i-th row and t-th column of $\mathrm{\Lambda }$. The i-th diagonal element of $\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }$,

$$h_i=\sum _{t=1}^k{\lambda }_{it}^2,$$

(3)

represents the proportion of the variance of $X_i$ explained by the common factors and is referred to as the communality of variable $X_i$.

The matrix $\mathrm{\Sigma }$ need not represent the variance-covariance matrix of a single group of observable data series; rather, it may characterize the covariance structure of residuals arising from a system of equations. To illustrate this point, consider the following structural vector autoregressive (SVAR) model:

$$Y_t=c+\sum _{s=1}^{\zeta }{A}_s{Y}_{t-s}+\mathrm{\Lambda }{F}_t+{ϵ}_t.$$

(4)

Here, $Y_t$ is a p-dimensional vector of observed time series at time t; $\zeta$ denotes the lag length; $F_t$ represents k common structural shocks; and ${ϵ}_t$ captures p idiosyncratic shocks. The composite disturbance term $\mathrm{\Lambda }{F}_t+{ϵ}_t$ therefore exhibits a two-component error structure. However, the common factors F in Eq. (2) may possess a nonzero mean. For panel data, Eq. (4) can be extended to

$$Y_{it}=c_i+\sum _{s=1}^{\zeta }{A}_s{Y}_{i,t-s}+{\mathrm{\Lambda }}_i{F}_t+{ϵ}_{it},$$

where ${\mathrm{\Lambda }}_i$ is a cross-section-specific loading matrix and $c_i$ is a cross-section-specific intercept vector .

In much of the SVAR literature, idiosyncratic shocks are often ignored, effectively imposing the restriction $k=p$. Allowing for idiosyncratic shocks, however, substantially reduces both the number of structural shocks and the number of identifying conditions required. From a computational perspective, incorporating an error-component structure—such as a common-factor model—into SVAR systems is straightforward. For example, the parameters c and $A_s$ in Eq. (4) can be estimated by least-squares under the assumption that the variance-covariance matrix of $\mathrm{\Lambda }{F}_t+{ϵ}_t$ is time-invariant. As a byproduct, this procedure also delivers an estimate of the variance-covariance matrix, denoted $\widehat{\mathrm{\Sigma }}$, together with its sampling covariance, $\mathrm{acov}\left(\widehat{\mathrm{\Sigma }}\right)$. In particular, when $\mathrm{\Lambda }{F}_t+{ϵ}_t$ follows a multivariate normal distribution, $\widehat{\mathrm{\Sigma }}$ has a scaled Wishart distribution. Unfortunately, the uncertainty associated with $\widehat{\mathrm{\Sigma }}$ has not been exploited in defining structural shocks nor incorporated into subsequent inference, such as impulse-response analysis.

This form of dimensionality reduction is particularly valuable in environments characterized by highly interconnected data. In macroeconomic modeling, structural shocks are typically assumed to be independent. Yet, although an economy can be examined from multiple perspectives—thereby generating a large number of interrelated observable time series—its dynamics are generally driven by only a small set of underlying shocks, such as demand, supply, geopolitical, and technological disturbances.

A wide range of factor-extraction methods has been proposed in the literature (Harman, 1976, Lawley & Maxwell, 1971, SAS, 2015), including principal-factor, maximum-likelihood, and least-squares approaches. These methods seek to recover the factor model in Eq. (1) from a given sample covariance or correlation matrix $\widehat{\mathrm{\Sigma }}$. From a practical perspective, it is important not only to obtain point estimates of factor loadings and uniquenesses, but also to assess the precision of these estimates and to identify substantively meaningful—often rotated—factor structures. Such assessments require knowledge of the asymptotic covariance matrices of the estimators, as well as an understanding of how these covariances are affected by the chosen rotation criterion.

The asymptotic covariance of the estimators depends on both the factor-extraction method and the covariance matrix of $\widehat{\mathrm{\Sigma }}$. Let n denote the sample size, and suppose that the estimator of the loading matrix can be expressed as $\mathrm{\Lambda }=f\left(\mathrm{\Sigma }\right)$ for a given extraction method f. Assume that f is differentiable at $\mathrm{\Sigma }$, with differential $df$, which defines a linear mapping from the space of $p\times p$ symmetric matrices to the space of $p\times k$ matrices. By the delta method (Rao, 2009), the difference between $n^{1/2}d\mathrm{\Lambda }$ and $df\left(n^{1/2}d\mathrm{\Sigma }\right)$ converges to zero in probability as $n\to \infty$. In this setting, $d\mathrm{\Lambda }\approx \widehat{\mathrm{\Lambda }}-\mathrm{\Lambda }$ and $d\mathrm{\Sigma }\approx \widehat{\mathrm{\Sigma }}-\mathrm{\Sigma }$ represent deviations from the unobservable population parameters $\mathrm{\Sigma }$ and $\mathrm{\Lambda }$, respectively.

Let $\mathrm{\Sigma }=\left[{\sigma }_{ij}\right]$ and define $d{\sigma }_{ij}={\widehat{\sigma }}_{ij}-{\sigma }_{ij}$, where ${\widehat{\sigma }}_{ij}$ denotes the sample covariance (or correlation) between $X_i$ and $X_j$. The asymptotic covariance between elements of the estimated loading matrix $\widehat{\mathrm{\Lambda }}$ can then be expressed as

$$\mathrm{acov}({\widehat{\lambda }}_{ir},{\widehat{\lambda }}_{js})=\sum _{m,n,x,y=1}^p{\displaystyle \frac{\partial {\lambda }_{ir}}{\partial {\sigma }_{mx}}}\mathrm{acov}({\widehat{\sigma }}_{mx},{\widehat{\sigma }}_{ny}){\displaystyle \frac{\partial {\lambda }_{js}}{\partial {\sigma }_{ny}}}.$$

Similarly, let $\psi ={({\psi }_1,\dots ,{\psi }_p)}^{\top }$ denote the vector of diagonal elements of the uniqueness matrix $\mathrm{\Psi }$. The asymptotic covariance of the uniqueness estimators is given by

$$\mathrm{acov}({\widehat{\psi }}_i,{\widehat{\psi }}_j)=\sum _{m,n,x,y=1}^p{\displaystyle \frac{\partial {\psi }_i}{\partial {\sigma }_{mx}}}\mathrm{acov}({\widehat{\sigma }}_{mx},{\widehat{\sigma }}_{ny}){\displaystyle \frac{\partial {\psi }_j}{\partial {\sigma }_{ny}}}.$$

In the following five sections, we derive explicit expressions for the differentials $d\mathrm{\Lambda }$ and $d\psi$ as linear functions of $d\mathrm{\Sigma }$ for several widely used factor-extraction methods. Given the covariance matrix of $\widehat{\mathrm{\Sigma }}$, these expressions allow practitioners to compute standard errors for the estimated factor loadings $\widehat{\mathrm{\Lambda }}$ and uniquenesses $\widehat{\psi }$. For data drawn from a multivariate normal population, the covariance of $\widehat{\mathrm{\Sigma }}$ is well established (Girshick, 1939). For nonnormal data, asymptotic results for sample covariances or correlations are also available in the literature (de Leeuw, 1983, Hsu, 1949). Moreover, delta-method-based inference for rotated factor loadings (Archer & Jennrich, 1973, Jennrich, 1973, Yung & Hayashi, 2001) relies on these differentials through the chain rule of differentiation. While earlier studies have examined the asymptotic properties of principal component analysis (Anderson, 1963, Girshick, 1939) and derived standard-error formulas for maximum-likelihood factor analysis (Jennrich & Thayer, 1973, Lawley, 1967), our focus is on least-squares, iterative principal component, principal, alpha, and image factor solutions as implemented in SAS (2015).

Despite the availability of maximum-likelihood (ML) factor analysis, non–ML extraction methods remain widely used and broadly implemented in statistical and econometric software (Lin et al., 2023, Raykov et al., 2026). In many empirical applications, factor analysis is applied to discrete or qualitative data, which often exhibit high-frequency dynamics, high-dimensionality, or heavy-tailed distributions. Although recent studies have emphasized robustness, high-dimensional asymptotics, and software-based inference for rotated solutions, explicit closed-form asymptotic covariance formulas for non-ML factor extraction methods remain largely unavailable. This methodological gap substantially limits formal statistical inference and rigorous model comparison.

3. Iterative Principal Component and Principal Factor Analysis

We begin by introducing notation. For a $p\times k$ matrix $\mathrm{\Lambda }$, let $\overrightarrow{\mathrm{\Lambda }}$ denote its vectorization obtained by stacking the columns of $\mathrm{\Lambda }$, that is,

$$\overrightarrow{\mathrm{\Lambda }}={({\lambda }_{11},\dots ,{\lambda }_{p1},{\lambda }_{12},\dots ,{\lambda }_{p2},\dots ,{\lambda }_{1k},\dots ,{\lambda }_{pk})}^{\top }.$$

When no confusion arises, we continue to use double subscripts to refer to individual elements of $\overrightarrow{\mathrm{\Lambda }}$. In particular, the $(i,r)$-th element of $\mathrm{\Lambda }$, denoted by ${\lambda }_{ir}$, corresponds to the $\left(\right(r-1)p+i)$-th element of $\overrightarrow{\mathrm{\Lambda }}$. For a square matrix M, let $\mathrm{diag}\left(M\right)$ denote the vector formed by the diagonal elements of M, and let $\mathrm{Diag}\left(M\right)$ denote the diagonal matrix whose diagonal is given by $\mathrm{diag}\left(M\right)$. We also use ${\delta }_{ij}$ to denote the Kronecker delta, which equals 1 if $i=j$ and 0 otherwise.

Both iterative principal component factor analysis and principal factor analysis construct the loading matrix $\mathrm{\Lambda }$ from the principal eigenvectors of $\mathrm{\Sigma }-\mathrm{\Psi }$:

$$\mathrm{\Sigma }-\mathrm{\Psi }=\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }=\sum _{r=1}^k{\lambda }_r{\lambda }_r^{\top },$$

(5)

where ${\lambda }_r$ denotes the r-th column of $\mathrm{\Lambda }=[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k]$. These columns are ordered according to the eigenvalues of $\mathrm{\Sigma }-\mathrm{\Psi }$, with ${\lambda }_r$ corresponding to the r-th largest eigenvalue. Without loss of generality, we assume that these eigenvalues are distinct, which is typically the case for a sample variance-covariance matrix. Consequently, the inner product between any two distinct columns of $\mathrm{\Lambda }$ is zero. In particular, Eq. (5) implies

$${\sigma }_{ii}-{\psi }_i=\sum _{t=1}^k{\lambda }_{it}^2,\forall i=1,...,p.$$

(6)

For the sample variance-covariance matrix $\widehat{\mathrm{\Sigma }}$, Eq. (5) may, in general, include additional items of the form ${\widehat{\lambda }}_r{\widehat{\lambda }}_r^{\top }$ for $r>k$. These items are omitted in accordance with the model specification in Eq. (1).

Using Eq. (5), we obtain

$$\left(\mathrm{\Sigma }-\mathrm{\Psi }\right){\lambda }_r=\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }{\lambda }_r=\sum _{s=1}^k{\lambda }_s{\lambda }_s^{\top }{\lambda }_r=\sum _{s=1}^k{\lambda }_s\left({\lambda }_s^{\top }{\lambda }_r\right)=\left({\lambda }_r^{\top }{\lambda }_r\right){\lambda }_r.$$

(7)

This result holds for any $1\le r\le k$. Now fix indices $1\le i\le p$ and $1\le r\le k$. Applying Eq. (6) to the i-th components of both sides of Eq. (7) yields

$$\sum _{j\ne i}{\sigma }_{ij}{\lambda }_{jr}+{\lambda }_{ir}\sum _{t=1}^k{\lambda }_{it}^2={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}^2.$$

Taking differentials on both sides, we obtain

$$resizeboxwidth418.45818pt\sum _{j\ne i}\left({\sigma }_{ij}d{\lambda }_{jr}+{\lambda }_{jr}d{\sigma }_{ij}\right)+\left(\sum _{t=1}^k{\lambda }_{it}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}=\left(\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}d{\lambda }_{jr}.$$

After algebraic simplification, this expression reduces to

$$\sum _{j\ne i}\left({\sigma }_{ij}-2{\lambda }_{ir}{\lambda }_{jr}\right)d{\lambda }_{jr}+\left(\sum _{t=1}^k{\lambda }_{it}^2-\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+2{\lambda }_{ir}\sum _{t\ne r}{\lambda }_{it}d{\lambda }_{it}=-\sum _{j\ne i}{\lambda }_{jr}d{\sigma }_{ij}.$$

Collecting all linear equations derived above, the system can be written compactly in matrix form as

$$Ad\overrightarrow{\mathrm{\Lambda }}=Bd\overrightarrow{\mathrm{\Sigma }},$$

where the coefficient of $d{\lambda }_{jt}$ in the $(i,r)$-th equation is given by

$$A_{ir,jt}=\left\{\begin{array}{cc}{\displaystyle \sum _{s=1}^k}{\lambda }_{is}^2-{\displaystyle \sum _{s=1}^p}{\lambda }_{sr}^2,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ 2{\lambda }_{ir}{\lambda }_{it},\hfill & \mathrm{if}j=i\mathrm{but}t\ne r,\hfill \\ {\sigma }_{ij}-2{\lambda }_{ir}{\lambda }_{jr},\hfill & \mathrm{if}j\ne i\mathrm{but}t=r,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

for all $1\le i,j\le p$ and $1\le r,t\le k$. Similarly, the coefficient of $d{\sigma }_{xy}$ in the same equation is

$$B_{ir,xy}=-{\delta }_{ix}(1-{\delta }_{xy}){\lambda }_{yr},\forall 1\le i,x,y\le p,1\le r\le k.$$

Finally, we obtain

$$d\overrightarrow{\mathrm{\Lambda }}=A^{-1}Bd\overrightarrow{\mathrm{\Sigma }}$$

and, by the delta method, the asymptotic covariance matrix of $\overrightarrow{\widehat{\mathrm{\Lambda }}}$ is given by

$$\mathrm{acov}\left(\overrightarrow{\widehat{\mathrm{\Lambda }}}\right)=\left(A^{-1}B\right)\mathrm{acov}\left(\overrightarrow{\widehat{\mathrm{\Sigma }}}\right){\left(A^{-1}B\right)}^{\top }.$$

4. Least-Squares Factor Analysis

In least-squares factor analysis (Harman, 1976), the objective is to minimize the sum of squared off-diagonal elements of the residual covariance matrix $\widehat{\mathrm{\Sigma }}-\widehat{\mathrm{\Lambda }}{\widehat{\mathrm{\Lambda }}}^{\top }$. Define the objective function

$$\mathrm{g}\left(\mathrm{\Lambda }\right)=\sum _{1\le y<x\le p}{\left({\sigma }_{xy}-\sum _{t=1}^k{\lambda }_{xt}{\lambda }_{yt}\right)}^2,$$

that is, the sum of squares of the strictly lower-triangular elements of $\mathrm{\Sigma }-\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }$.

The first-order optimality conditions yield $p\times k$ constraints on the loading matrix$\mathrm{\Lambda }$. Specifically, we require

$${\displaystyle \frac{\partial \mathrm{g}\left(\mathrm{\Lambda }\right)}{\partial {\lambda }_{ir}}}=0$$

for all $1\le i\le p$ and $1\le r\le k$. Straightforward differentiation gives

$$\begin{array}{ccc}\hfill -{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial \mathrm{g}}{\partial {\lambda }_{ir}}}\left(\mathrm{\Lambda }\right)& =& {\displaystyle \sum _{1\le y<x\le p}}\left({\sigma }_{xy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{yt}\right){\displaystyle \sum _{t=1}^k}({\delta }_{ix}{\delta }_{rt}{\lambda }_{yt}+{\delta }_{iy}{\delta }_{rt}{\lambda }_{xt})\hfill \\ & =& {\displaystyle \sum _{1\le y<x\le p}}({\delta }_{ix}{\lambda }_{yr}+{\delta }_{iy}{\lambda }_{xr})\left({\sigma }_{xy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{yt}\right)\hfill \\ & =& {\displaystyle \sum _{1\le y<i}}{\lambda }_{yr}\left({\sigma }_{iy}-{\displaystyle \sum _{t=1}^k}{\lambda }_{it}{\lambda }_{yt}\right)+{\displaystyle \sum _{i<x\le p}}{\lambda }_{xr}\left({\sigma }_{xi}-{\displaystyle \sum _{t=1}^k}{\lambda }_{xt}{\lambda }_{it}\right)\hfill \\ & =& {\displaystyle \sum _{j\ne i}}{\lambda }_{jr}\left({\sigma }_{ij}-{\displaystyle \sum _{t=1}^k}{\lambda }_{it}{\lambda }_{jt}\right).\hfill \end{array}$$

For any fixed $1\le i\le p$ and $1\le r\le k$, the first-order condition reduces to

$$\sum _{j\ne i}{\lambda }_{jr}{\sigma }_{ij}=\sum _{j\ne i}{\lambda }_{jr}\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}.$$

(8)

The differential form of Eq. (8) is given by

$$\sum _{j\ne i}\left({\lambda }_{jr}d{\sigma }_{ij}+{\sigma }_{ij}d{\lambda }_{jr}\right)=\sum _{j\ne i}\left(\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}\right)d{\lambda }_{jr}+\sum _{j\ne i}{\lambda }_{jr}\sum _{t=1}^k\left({\lambda }_{jt}d{\lambda }_{it}+{\lambda }_{it}d{\lambda }_{jt}\right).$$

Grouping terms and simplifying yield

$$resizeboxwidth418.45818pt\sum _{j\ne i}{\lambda }_{jr}d{\sigma }_{ij}=\sum _{j\ne i}\left(\sum _{t=1}^k{\lambda }_{it}{\lambda }_{jt}-{\sigma }_{ij}+{\lambda }_{jr}{\lambda }_{ir}\right)d{\lambda }_{jr}+\sum _{t=1}^k\left(\sum _{j\ne i}{\lambda }_{jr}{\lambda }_{jt}\right)d{\lambda }_{it}+\sum _{j\ne i}\sum _{t\ne r}{\lambda }_{jr}{\lambda }_{it}d{\lambda }_{jt},$$

which completes the differential characterization of the first-order condition.

We rewrite the preceding system of equations in matrix form as

$$-Bd\overrightarrow{\mathrm{\Sigma }}=Dd\overrightarrow{\mathrm{\Lambda }},$$

where the coefficient of $d{\lambda }_{jt}$ in the $(i,r)$-th equation is defined as

$$D_{ir,jt}=\left\{\begin{array}{cc}{\displaystyle \sum _{z\ne i}}{\lambda }_{zr}{\lambda }_{zt},\hfill & \mathrm{if}j=i,\hfill \\ {\displaystyle \sum _{z=1}^k}{\lambda }_{iz}{\lambda }_{jz}-{\sigma }_{ij}+{\lambda }_{ir}{\lambda }_{jr},\hfill & \mathrm{if}j\ne i\mathrm{and}t=r,\hfill \\ {\lambda }_{jr}{\lambda }_{it},\hfill & \mathrm{if}j\ne i\mathrm{and}t\ne r.\hfill \end{array}\right.$$

These definitions apply for all $1\le i,j\le p$ and $1\le r,t\le k$. Finally, solving for $d\overrightarrow{\mathrm{\Lambda }}$ yields

$$d\overrightarrow{\mathrm{\Lambda }}=-D^{-1}Bd\overrightarrow{\mathrm{\Sigma }}.$$

5. Alpha Factor Analysis

Let H be the diagonal matrix whose i-th diagonal element is $\sqrt{h_i}$. Clearly, $h_i>0$; otherwise, Eq. (3) would imply that the i-th row of the loading matrix $\mathrm{\Lambda }$ is a zero vector, rendering the corresponding factor redundant. Moreover, $H^{-1}(\mathrm{\Sigma }-\mathrm{\Psi })H^{-1}$ is the correlation matrix for $\mathrm{\Lambda }F$ and Eq (3) implies that

$$d{h}_i=2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}.$$

(9)

To obtain the spectral decomposition of this correlation matrix, let ${\theta }_r$ denote its r-th largest eigenvalue, and let ${\beta }_r$ be the corresponding eigenvector with length ${\theta }_r^{1/2}$. Thus,

$${\theta }_r={\beta }_r^{\top }{\beta }_r,$$

and

$$H^{-1}(\mathrm{\Sigma }-\mathrm{\Psi })H^{-1}{\beta }_r={\theta }_r{\beta }_r,\forall 1\le r\le k.$$

In alpha factor analysis (Kaiser & Caffrey, 1965), the loading matrix $\mathrm{\Lambda }=[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k]$ is defined as

$$\mathrm{\Lambda }=H[{\beta }_1,{\beta }_2,\dots ,{\beta }_k].$$

Additionally, let $\Gamma =[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_k]$, where ${\gamma }_r=H^{-1}{\beta }_r$. Then, for each $r=1,\dots ,k$, the following relationships—corresponding to Eqs. (10)-(12)—hold:

$${\lambda }_r=H{\beta }_r=H^2{\gamma }_r,$$

(10)

$${\theta }_r={\beta }_r^{\top }{\beta }_r={\gamma }_r^{\top }{H}^2{\gamma }_r={\gamma }_r^{\top }{\lambda }_r,$$

(11)

$$(\mathrm{\Sigma }-\mathrm{\Psi }){\gamma }_r=(\mathrm{\Sigma }-\mathrm{\Psi })H^{-1}{\beta }_r=H{H}^{-1}(\mathrm{\Sigma }-\mathrm{\Psi })H^{-1}{\beta }_r={\theta }_{r}H{\beta }_r={\gamma }_r^{\top }{\lambda }_r{\lambda }_r.$$

(12)

In the remainder of this section, we first express $d\overrightarrow{\Gamma }$ in terms of $d\overrightarrow{\mathrm{\Lambda }}$ using Eq. (10). We then derive a matrix equation relating $d\overrightarrow{\mathrm{\Sigma }}$, $d\overrightarrow{\Gamma }$, and $d\overrightarrow{\mathrm{\Lambda }}$ from Eq. (12).

From Eq. (10), for any $1\le i\le p$ and $1\le r\le k$, we have

$${\lambda }_{ir}=h_i{\gamma }_{ir},$$

whose differential is given by

$$d{\lambda }_{ir}=h_{i}d{\gamma }_{ir}+2{\gamma }_{ir}\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}.$$

Collecting these expressions yields the following matrix representation:

$$d\overrightarrow{\Gamma }=Ed\overrightarrow{\mathrm{\Lambda }},$$

(13)

where the elements of the matrix E are defined by

$$E_{ir,jt}={\displaystyle \frac{\partial {\gamma }_{ir}}{\partial {\lambda }_{jt}}}={\delta }_{ij}({\delta }_{rt}-2{\gamma }_{ir}{\lambda }_{it})/h_i,$$

for all $1\le i,j\le p$ and $1\le r,t\le k$.

For any fixed $1\le i\le p$ and $1\le r\le k$, the i-th component of Eq. (12) can be written as

$$\sum _{j\ne i}{\sigma }_{ij}{\gamma }_{jr}+h_i{\gamma }_{ir}={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}{\gamma }_{jr}.$$

(14)

Taking differentials of Eq. (14) yields

$$resizeboxwidth418.45818pt\sum _{j\ne i}\left({\sigma }_{ij}d{\gamma }_{jr}+{\gamma }_{jr}d{\sigma }_{ij}\right)+h_{i}d{\gamma }_{ir}+{\gamma }_{ir}d{h}_i={\lambda }_{ir}\sum _{j=1}^p\left({\lambda }_{jr}d{\gamma }_{jr}+{\gamma }_{jr}d{\lambda }_{jr}\right)+\left(\sum _{j=1}^p{\lambda }_{jr}{\gamma }_{jr}\right)d{\lambda }_{ir}.$$

Substituting Eq. (9) into the above expression and rearranging terms, we obtain

$$resizeboxwidth418.45818pt\sum _{j\ne i}{\gamma }_{jr}d{\sigma }_{ij}=\sum _{j\ne i}({\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij})d{\gamma }_{jr}+({\lambda }_{ir}^2-h_i)d{\gamma }_{ir}+\left(\sum _{j\ne i}{\lambda }_{jr}{\gamma }_{jr}\right)d{\lambda }_{ir}-2{\gamma }_{ir}\sum _{t\ne r}{\lambda }_{it}d{\lambda }_{it}+{\lambda }_{ir}\sum _{j\ne i}{\gamma }_{jr}d{\lambda }_{jr}.$$

We express all the differential forms above in matrix notation as

$$Fd\overrightarrow{\mathrm{\Sigma }}=Gd\overrightarrow{\Gamma }+Jd\overrightarrow{\mathrm{\Lambda }}.$$

(15)

The coefficient matrices F, G, and J are defined as follows.

• Matrix F: The coefficient of $d{\sigma }_{xy}$ in the $(i,r)$-th equation of Eq. (15) is

  $$F_{ir,xy}={\delta }_{ix}(1-{\delta }_{iy}){\gamma }_{yr}.$$
• Matrix G: The coefficient of $d{\gamma }_{jt}$ is given by

  $$G_{ir,jt}=\left\{\begin{array}{cc}{\lambda }_{ir}^2-h_i,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ {\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$
• Matrix J: The coefficient of $d{\lambda }_{jt}$ is

  $$J_{ir,jt}=\left\{\begin{array}{cc}{\displaystyle \sum _{s\ne i}}{\lambda }_{sr}{\gamma }_{sr},\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ -2{\gamma }_{ir}{\lambda }_{it},\hfill & \mathrm{if}j=i\mathrm{but}t\ne r,\hfill \\ {\lambda }_{ir}{\gamma }_{jr},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Combining Eqs. (13) and (15), we obtain

$$Fd\overrightarrow{\mathrm{\Sigma }}=Gd\overrightarrow{\Gamma }+Jd\overrightarrow{\mathrm{\Lambda }}=G\left(Ed\overrightarrow{\mathrm{\Lambda }}\right)+Jd\overrightarrow{\mathrm{\Lambda }}=\left(GE+J\right)\overrightarrow{\mathrm{\Lambda }}.$$

Consequently,

$$d\overrightarrow{\mathrm{\Lambda }}={(GE+J)}^{-1}Fd\overrightarrow{\mathrm{\Sigma }}.$$

6. Image Factor Analysis

Let $\mathsf{\Delta }={\left(\mathrm{Diag}\left({\mathrm{\Sigma }}^{-1}\right)\right)}^{-1}$. In image factor analysis (Jöreskog, 1969), the covariance matrix $\mathrm{\Sigma }$ is modeled as

$$\mathrm{\Sigma }=\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top }+\tau \mathsf{\Delta },$$

(16)

where $\tau >0$ is an unknown scalar parameter. Within this framework, one may write $\mathrm{\Psi }=\tau \mathsf{\Delta }$. Jöreskog (1969) proposed a maximum-likelihood procedure for estimating this model, under which the covariance matrix of the estimators follows directly from standard maximum-likelihood theory.

For the principal-factor procedure (SAS, 2015), the loading matrix $\mathrm{\Lambda }$ is obtained from the k largest eigenvalues and their corresponding eigenvectors of $\mathrm{\Sigma }-\tau \mathsf{\Delta }$. The parameter $\tau$ is simultaneously estimated as the coefficient in the regression of $\mathrm{diag}(\mathrm{\Sigma }-\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top })$ on $\mathrm{diag}\left(\mathsf{\Delta }\right)$, namely,

$$\tau ={\displaystyle \frac{\langle \mathrm{diag}\left(\mathsf{\Delta }\right),\mathrm{diag}(\mathrm{\Sigma }-\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top })\rangle }{\langle \mathrm{diag}\left(\mathsf{\Delta }\right),\mathrm{diag}\left(\mathsf{\Delta }\right)\rangle }}$$

(17)

where $\langle ·,·\rangle$ denotes the standard inner product.

Let ${\sigma }^{xy}$ denote the element at the x-th row and y-th column of ${\mathrm{\Sigma }}^{-1}$. Since $\mathrm{\Sigma }{\mathrm{\Sigma }}^{-1}$ is the $p\times p$ identity matrix, differentiating this identity using the product rule yields

$$d{\mathrm{\Sigma }}^{-1}=-{\mathrm{\Sigma }}^{-1}\left(d\mathrm{\Sigma }\right){\mathrm{\Sigma }}^{-1}.$$

In particular, for each $i=1,...,p$,

$$d{\sigma }^{ii}=-\sum _{x,y=1}^p{\sigma }^{ix}\left(d{\sigma }_{xy}\right){\sigma }^{yi}=-\sum _{x,y=1}^p{\sigma }^{ix}{\sigma }^{iy}d{\sigma }_{xy}.$$

Consequently,

$$d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=-{\displaystyle \frac{1}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }^{ii}=\sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy}.$$

(18)

Eq. (17) implies that

$$\tau \langle \mathrm{diag}\left(\mathsf{\Delta }\right),\mathrm{diag}\left(\mathsf{\Delta }\right)\rangle =\langle \mathrm{diag}\left(\mathsf{\Delta }\right),\mathrm{diag}(\mathrm{\Sigma }-\mathrm{\Lambda }{\mathrm{\Lambda }}^{\top })\rangle ,$$

and therefore

$$\tau \sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2=\sum _{i=1}^p{\displaystyle \frac{1}{{\sigma }^{ii}}}({\sigma }_{ii}-h_i).$$

Differentiating both sides yields

$$\left[\sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2\right]d\tau +\sum _{i=1}^p{\displaystyle \frac{2\tau }{{\sigma }^{ii}}}d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=\sum _{i=1}^p\left[{\displaystyle \frac{1}{{\sigma }^{ii}}}(d{\sigma }_{ii}-d{h}_i)+({\sigma }_{ii}-h_i)d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)\right].$$

Substituting Eq. (18) into the above expression, we obtain

$$resizeboxwidth418.45818pt\left[\sum _{i=1}^p{\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)}^2\right]d\tau =\sum _{x,y=1}^p\left[\sum _{i=1}^p\left({\sigma }_{ii}-h_i-{\displaystyle \frac{2\tau }{{\sigma }^{ii}}}\right){\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}\right]d{\sigma }_{xy}+\sum _{i=1}^p{\displaystyle \frac{1}{{\sigma }^{ii}}}\left(d{\sigma }_{ii}-2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it}\right).$$

These expressions can be written compactly as

$$d\tau ={\mu }^{\top }d\overrightarrow{\mathrm{\Sigma }}+{\eta }^{\top }d\overrightarrow{\mathrm{\Lambda }},$$

(19)

where $\mu$ and $\eta$ are vectors of dimensions $p^2$ and $pk$, respectively. The coefficient of $d{\sigma }_{xy}$ is given by

$${\mu }_{xy}={\displaystyle \frac{\frac{{\delta }_{xy}}{{\sigma }^{xx}}+{\displaystyle \sum _{i=1}^p}\left({\sigma }_{ii}-h_i-\frac{2\tau }{{\sigma }^{ii}}\right)\frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}{{\displaystyle \sum _{i=1}^p}{\left(\frac{1}{{\sigma }^{ii}}\right)}^2}}$$

and the coefficient of $d{\lambda }_{ir}$ is

$${\eta }_{ir}=-{\displaystyle \frac{2{\lambda }_{ir}}{{\sigma }^{ii}{\displaystyle \sum _{z=1}^p}{\left(\frac{1}{{\sigma }^{zz}}\right)}^2}}.$$

Next, let ${\theta }_r$ denote the r-th largest eigenvalue of $\mathrm{\Sigma }-\tau \mathsf{\Delta }$, and let ${\lambda }_r$ be the associated eigenvector with length ${\theta }_r^{1/2}$. Then,

$${\lambda }_r^{\top }{\lambda }_r={\theta }_r,$$

and

$$\left(\mathrm{\Sigma }-\tau \mathsf{\Delta }\right){\lambda }_r={\theta }_r{\lambda }_r={\lambda }_r^{\top }{\lambda }_r{\lambda }_r.$$

For any fixed $1\le i\le p$, the i-th component of this equation is given by

$$\sum _{j=1}^p{\sigma }_{ij}{\lambda }_{jr}-{\displaystyle \frac{\tau {\lambda }_{ir}}{{\sigma }^{ii}}}={\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}^2.$$

Differentiating both sides with respect to the underlying parameters yields

$$\sum _{j=1}^p\left({\sigma }_{ij}d{\lambda }_{jr}+{\lambda }_{jr}d{\sigma }_{ij}\right)-{\displaystyle \frac{\tau }{{\sigma }^{ii}}}d{\lambda }_{ir}-{\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}}d\tau -\tau {\lambda }_{ir}d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)=2{\lambda }_{ir}\sum _{j=1}^p{\lambda }_{jr}d{\lambda }_{jr}+\left(\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}.$$

Applying Eq. (18) and collecting terms involving $d\mathrm{\Sigma },d\Gamma$, and $d\tau$, we obtain

$$resizeboxwidth418.45818pt\sum _{j=1}^p{\lambda }_{jr}d{\sigma }_{ij}-\tau {\lambda }_{ir}\sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy}=\sum _{j\ne i}(2{\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij})d{\lambda }_{jr}+\left(2{\lambda }_{ir}^2-{\sigma }_{ii}+{\displaystyle \frac{\tau }{{\sigma }^{ii}}}+\sum _{j=1}^p{\lambda }_{jr}^2\right)d{\lambda }_{ir}+{\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}}d\tau .$$

The matrix representation of the above system of equations can be written as

$$Ld\overrightarrow{\mathrm{\Sigma }}=Pd\overrightarrow{\mathrm{\Lambda }}+\pi d\tau$$

(20)

where L, P, and $\pi$ are appropriately defined coefficient matrices. In particular, the coefficient of $d{\sigma }_{xy}$ in the $(i,r)$-th equation is given by

$$L_{ir,xy}={\delta }_{ix}{\lambda }_{yr}-{\displaystyle \frac{\tau {\lambda }_{ir}{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}.$$

The coefficient of $d{\lambda }_{jt}$ is

$$P_{ir,jt}=\left\{\begin{array}{cc}2{\lambda }_{ir}^2-{\sigma }_{ii}+{\displaystyle \frac{\tau }{{\sigma }^{ii}}}+{\displaystyle \sum _{z=1}^p}{\lambda }_{zr}^2,\hfill & \mathrm{if}j=i\mathrm{and}t=r,\hfill \\ 2{\lambda }_{ir}{\lambda }_{jr}-{\sigma }_{ij},\hfill & \mathrm{if}t=r\mathrm{but}j\ne i,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

and the coefficient of $d\tau$ is

$${\pi }_{ir}={\displaystyle \frac{{\lambda }_{ir}}{{\sigma }^{ii}}},$$

for any $1\le i,j,x,y\le p$ and $1\le r,t\le k$. Note that $\pi$ is a vector of dimension $pk$.

Combining Eqs. (19) and (20) yields

$$d\overrightarrow{\mathrm{\Lambda }}={\left(P+\pi {\eta }^{\top }\right)}^{-1}\left(L-\pi {\mu }^{\top }\right)d\overrightarrow{\mathrm{\Sigma }}.$$

Substituting this expression back into Eq. (19) gives the differential of $\tau$ in terms of $d\overrightarrow{\mathrm{\Sigma }}$ :

$$d\tau =\left[{\mu }^{\top }+{\eta }^{\top }{\left(P+\pi {\eta }^{\top }\right)}^{-1}\left(L-\pi {\mu }^{\top }\right)\right]d\overrightarrow{\mathrm{\Sigma }}.$$

(21)

7. Standard Error of Uniqueness Estimator $\widehat{\psi }$

Assume that $d\overrightarrow{\mathrm{\Lambda }}=Md\overrightarrow{\mathrm{\Sigma }}$ for some $pk\times p^2$ matrix M. To derive the standard errors of the uniqueness estimator $\widehat{\psi }$, we begin with the generic relationship from Eq. (6): ${\psi }_i={\sigma }_{ii}-\sum _{t=1}^k{\lambda }_{it}^2.$ Differentiating both sides yields

$$d{\psi }_i=d{\sigma }_{ii}-2\sum _{t=1}^k{\lambda }_{it}d{\lambda }_{it},i=1,2,\dots ,p.$$

(22)

Eq. (22) can be written more compactly in matrix form as

$$d\psi =Td\overrightarrow{\mathrm{\Sigma }}+Zd\overrightarrow{\mathrm{\Lambda }},$$

(23)

where the coefficient of $d{\sigma }_{xy}$ in the i-th equation of Eq. (23) is given by

$$T_{i,xy}={\delta }_{ix}{\delta }_{iy}$$

and the coefficient of $d{\lambda }_{jt}$ in the same equation is

$$Z_{i,jt}=-2{\delta }_{ij}{\lambda }_{jt},$$

for all $1\le i,j,x,y\le p$ and $1\le t\le k$. Consequently, we obtain

$$d\psi =(T+ZM)d\overrightarrow{\mathrm{\Sigma }}.$$

When $\mathrm{\Sigma }$ is a correlation matrix, $d{\sigma }_{ii}=0$, and the expression simplifies further to $d\psi =ZMd\overrightarrow{\mathrm{\Sigma }}$.

In image factor analysis, the relationship ${\widehat{\psi }}_i={\widehat{\sigma }}_{ii}-{\displaystyle \sum _{t=1}^k}{\widehat{\lambda }}_{it}^2$ may no longer hold. Instead, we rewrite Eq. (21) as $d\tau ={\varphi }^{\top }d\overrightarrow{\mathrm{\Sigma }}$ for some vector $\varphi$. Using the identity ${\psi }_i={\displaystyle \frac{\tau }{{\sigma }^{ii}}}$ and applying Eq. (18), we obtain

$$d{\psi }_i={\displaystyle \frac{d\tau }{{\sigma }^{ii}}}+\tau d\left({\displaystyle \frac{1}{{\sigma }^{ii}}}\right)={\displaystyle \frac{1}{{\sigma }^{ii}}}\sum _{x,y=1}^p{\varphi }_{xy}d{\sigma }_{xy}+\tau \sum _{x,y=1}^p{\displaystyle \frac{{\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}}d{\sigma }_{xy},$$

(24)

for all $1\le i\le p$. In matrix form, this expression can be written compactly as

$$d\psi =Qd\overrightarrow{\mathrm{\Sigma }},$$

where the coefficient of $d{\sigma }_{xy}$ in the i-th equation is given by

$$Q_{i,xy}={\displaystyle \frac{\partial {\psi }_i}{\partial {\sigma }_{xy}}}={\displaystyle \frac{{\sigma }^{ii}{\varphi }_{xy}+\tau {\sigma }^{ix}{\sigma }^{iy}}{{\left({\sigma }^{ii}\right)}^2}},\forall 1\le i,x,y\le p.$$

8. An Empirical and Simulation Study

In this section, we apply the least-squares factor analysis method to a sample correlation matrix consisting of nine variables measured on 211 subjects. The data were originally reported in page 43 of Lawley and Maxwell (1971), and the sample is assumed to be drawn from an approximately multivariate normal population. For simplicity, we set the number of common factors to $k=2$. This choice is supported by both Horn (1965)’s parallel analysis and Velicer (1976)’s minimum average partial test . Estimates of the factor loading matrix $\widehat{\mathrm{\Lambda }}$ and the uniqueness vector $\widehat{\psi }$ are obtained using the iterative procedure described by Harman (1976).

To compute the asymptotic covariances of these estimators, we first apply Girshick (1939)’s formula (3.23) to derive the covariances of the sample correlation coefficients. We then use the results developed in Section 4 and Section 7 to compute the asymptotic covariance matrix of $\widehat{\mathrm{\Lambda }}$, as well as the standard errors of $\widehat{\psi }$. To illustrate the effect of factor rotation, we orthogonally rotate $\widehat{\mathrm{\Lambda }}$ using the varimax criterion and apply the method of Jennrich and Thayer (1973) to obtain the standard errors of the rotated factor loadings. Table 1 reports the parameter estimates for the two-factor least-squares model; standard errors are reported in parentheses.

To verify the theoretical derivation of the asymptotic covariance, $\mathrm{acov}\left(\widehat{\psi }\right)$, presented in Section 7, we conduct a simulation study. The sample correlation matrix reported above is treated as the population correlation matrix $\mathrm{\Sigma }$. Using this matrix, we generate a large number m of random correlation matrices from a Wishart distribution with parameter $(\mathrm{\Sigma },211)$. For each simulated correlation matrix, the uniqueness parameters are estimated using the least-squares method with $k=2$.

Based on the resulting collection of simulated estimates of $\psi$, we compute empirical standard errors for the uniqueness estimates. For $m=50,100,500,1000,2000$, the corresponding empirical standard errors are reported in Table 2. The final column of the table presents the theoretical asymptotic standard errors computed using the formulas developed in Section 7. The results demonstrate clear convergence of the empirical standard errors toward their theoretical counterparts, although the convergence occurs at a relatively slow rate.

9. Conclusion with Discussions

Recent years have witnessed renewed interest in formal statistical inference for factor-analytic models. This resurgence is driven by several practical considerations, including the prevalence of non-normal data, heightened concerns about robustness, and the growing complexity of covariance structures encountered in empirical applications. Broadly speaking, the contemporary literature can be organized into three main strands. The first strand emphasizes likelihood-based approaches to factor analysis and covariance structure modeling, building on classical frameworks and their subsequent extensions (Jennrich & Thayer, 1973, Lawley, 1967, Raykov et al., 2026). A second strand focuses on high-dimensional and dynamic factor models, in which both cross-sectional and temporal dimensions expand, giving rise to new challenges for estimation and inference (Bai, 2003, Fan et al., 2021). The third strand concentrates on software-oriented implementations of standard error estimation for rotated factor solutions, typically relying on numerical differentiation or sandwich-type approximations (Archer & Jennrich, 1973, Yung & Hayashi, 2001).

In contrast, explicit asymptotic covariance formulas for widely used non-maximum-likelihood extraction methods—such as least-squares, principal, alpha, and image factor analysis—remain relatively scarce, despite their continued prevalence in empirical research and standard statistical software (Harman, 1976, SAS, 2015). This paper helps bridge this gap by deriving closed-form, analytically tractable asymptotic covariance expressions for several classical no-ML factor solutions. By casting these results within a unified implicit-differentiation framework, the proposed approach accommodates non-normal data and can be readily integrated with existing delta-method results for factor rotation (Archer & Jennrich, 1973, Jennrich, 1973, Yung & Hayashi, 2001). From a practical standpoint, the framework enables transparent computation of standard errors without reliance on numerical differentiation or resampling methods. From a methodological perspective, it complements recent advances in robust and high-dimensional factor modeling by providing a rigorous inferential foundation for classical extraction methods that remain central to applied factor analysis.

In the preceding sections, the derived asymptotic covariances of $\widehat{\mathrm{\Lambda }}$ and $\widehat{\psi }$ are conditional on the assumed number of factors. Consequently, anomalous standard errors may reflect model misspecification or systematic errors—such as an incorrect choice of factor dimensionality—rather than genuine sampling variability. In particular, such distortions may arise from either over-factoring (when k is chosen too large) or under-factoring (when k is chosen too small).

For expositional simplicity, we do not explicitly impose the symmetry condition ${\sigma }_{ij}={\sigma }_{ji}$. In practical implementations, however, exploiting this property is important not only for reducing memory requirements but also for improving computational efficiency, especially when working with large-scale datasets.

For the SVAR model in Eq. (4), target or Procrustean rotation is unnecessary when the structural shocks $F_t$ are properly aligned with the endogenous variables $Y_t$ in the loading matrix. Such alignment naturally imposes identifying restrictions—such as zero or sign restrictions—on the elements of $\mathrm{\Lambda }$, which may require corresponding adjustments in its estimation. Nevertheless, the formulas developed in this paper remain valid, as they treat the estimated $\mathrm{\Lambda }$ as given inputs. Moreover, these formulas not only assess the precision of coefficient estimates associated with structural shocks, but also facilitate the evaluation of uncertainty in impulse response functions, forecast-error variance decompositions, historical decompositions, and the estimated factors $F_t$.

From an econometric perspective, the results developed in this paper are particularly relevant for empirical settings in which factor extraction serves as an intermediate step rather than an end in itself. In applied macroeconomics and finance, factor-analytic structures commonly arise in dynamic factor models, large Bayesian or classical VARs, factor-augmented regressions, and error-component representations of high-dimensional systems. Within these frameworks, factor loadings and idiosyncratic components are typically treated as nuisance parameters, and the uncertainty associated with their estimation is rarely propagated into subsequent stages of inference.

The explicit asymptotic covariance formulas derived in this paper enable researchers to systematically account for factor-estimation uncertainty in econometric inference. In particular, they facilitate the construction of delta-method-based standard errors for impulse response functions, forecast-error variance decompositions, counterfactual simulations, and other nonlinear functionals that depend on estimated factors. This consideration is especially important in structural VAR applications with latent common shocks, where factor estimates play a central role in both the identification and interpretation of structural disturbances.

Moreover, because the proposed expressions take the covariance of the sample covariance or correlation as inputs, they remain valid under heteroskedasticity, serial dependence, and non-Gaussian sampling schemes—features that are pervasive in macroeconomic and financial data. As a result, the framework naturally aligns with heteroskedasticity- and autocorrelation-robust covariance estimators widely used in econometrics, without relying on likelihood-based assumptions. For example, when the data are multivariate binary, the asymptotic covariance matrix, $\mathrm{acov}\left(\widehat{\mathrm{\Lambda }}\right)$, remains well defined. More generally, for categorical data—common in applied factor-analysis settings—the computation of $\mathrm{acov}\left(\widehat{\mathrm{\Lambda }}\right)$ is unlikely to pose practical difficulties.

Finally, the availability of closed-form expressions provides a valuable diagnostic tool for model specification in econometric practice. Abnormally large estimated standard errors may indicate misspecification of factor dimensionality, weak factor structure, or inadequate separation between common and idiosyncratic components. Moreover, a large proportion of statistically insignificant coefficients in $\widehat{\mathrm{\Lambda }}$ may simply reflect over-factoring. From this perspective, asymptotic covariance analysis not only supports formal statistical inference but also offers practical guidance for model selection and validation in factor-based econometric models.