I will note before starting that this question is related to my own work, in which i am aiming to extend and replicate the methods used in Attanasio et. al. (2015), "Human Capital Development and Parental Investment in India"**.
I am carrying out a large confirmatory factor analysis where the relationship between observed measurements and latent variables can be written as follows:
$$ y_{i,j} = \lambda_{i,j} \theta_j + \varepsilon_{i,j} $$
where $y_{i,j}$ is measurement $i$ on latent variable $j$, $\lambda_{i,j}$ is a factor loading, and $\varepsilon_{i,j}$ a measurement error. My aim is to estimate the joint distribution of the latent factors - the $\theta_j$. For my application i cannot assume normality of this distribution, so instead i assume their joint distribution to be a mixture of two normals. I also assume that the latent variables are independent of measurement errors and that the latter are normally distributed. First, writing the measurement system in matrix form, we have
$$ \textbf{y} = \Lambda \boldsymbol{\theta} + \boldsymbol{\varepsilon} $$
where $\textbf{y}$ is a vector of all measurements on each latent factor, $\boldsymbol{\Lambda}$ is a matrix containing all the factor loadings, $\boldsymbol{\theta}$ a vector containing all latent variables, and $\boldsymbol{\varepsilon}$ is a vector containing all the measurement errors. Using basic rules about the distribution of sums of independent random variables and the convolution operator, the theoretical joint distribution of the of the measurements, $p(\textbf{y})$ is given by:
$$p(\textbf{y}) = \tau \underbrace{\int g(\textbf{y} -\boldsymbol{\Lambda }\boldsymbol{\theta})f_A(\boldsymbol{\Lambda}\boldsymbol{\theta})d\boldsymbol{\theta}}_{p_A({\cdot})} + (1-\tau) \underbrace{ \int g(\textbf{y} -\boldsymbol{\Lambda }\boldsymbol{\theta})f_B(\boldsymbol{\Lambda}\boldsymbol{\theta})d\boldsymbol{\theta}}_{p_B({\cdot})}$$
with $f_A(\cdot)$ and $f_B(\cdot)$ being multivariate normal probability density functions representing the two mixture components of the joint distribution of the latent factors and $g(\cdot)$ is the joint density of the measurement errors. With the assumption made $p(\textbf{y})$ is also a mixture of two normals.
I then use a ML approximation of the joint distribution of measurements (assuming it is in fact a mixture of two normals) and obtain estimates of the moments of $p_A(\cdot)$ and $p_B({\cdot})$, denoted $\boldsymbol{\mu}^y_A$, $\boldsymbol{\mu}^y_B$, $\boldsymbol{\Sigma}^A_y$, $\boldsymbol{\Sigma}^B_y$, and $\tau$ where A,and B represent the first a second mixture components..
Once these are estimated, i would like to set them equal to the implied structure of the distribution given in the above equation. As such i am left with the following equalities:
$$ \tau \boldsymbol{\mu}_A^{\boldsymbol{\theta}} + (1-\tau)\boldsymbol{\mu}_b^{\boldsymbol{\theta}} = 0 \nonumber \\[5pt] \boldsymbol{\Lambda}\boldsymbol{\mu}^{\boldsymbol{\theta}}_A = \hat{\boldsymbol{\mu}}_A^{\tilde{\textbf{y}}} \nonumber \\[5pt] \label{momentconditions} \boldsymbol{\Lambda}\boldsymbol{\mu}^{\boldsymbol{\theta}}_B = \hat{\boldsymbol{\mu}}_B^{\tilde{\textbf{y}}} \\[5pt] \boldsymbol{\Lambda} \boldsymbol{\Theta}_A \boldsymbol{\Lambda}' + \boldsymbol{\Psi} = \hat{\boldsymbol{\Sigma}}_A \nonumber \\[5pt] \boldsymbol{\Lambda} \boldsymbol{\Theta}_B \boldsymbol{\Lambda}' + \boldsymbol{\Psi} = \hat{\boldsymbol{\Sigma}}_B \nonumber $$
Where $\boldsymbol{\Lambda}$ is still the matrix of factor loadings, and $\boldsymbol{\mu}_j^{\boldsymbol{\theta}}$ $\boldsymbol{\Theta}_j$ for $j \in \{ A,B \}$ are the vector of means and the covariance matrix from the corresponding component of the mixed joint distribution of the latent factors. All of the quantities on the right hand side of the above equations and $tau$ have been estimated.
The matrices $\boldsymbol{\Lambda}$, $\boldsymbol{\Theta}_j$ and $\boldsymbol{\Psi}$ have a structure implied by the measurement equation, and so my question is: can their parameters be estimated via minimum distance or any other suitable method in Matlab or STATA? If anyone knows of existing functions that can be used to do so, or how to go about this manually, then any comments or advice would be much appreciated.
Thanks in advance.
** URL to the original paper: http://egcenter.economics.yale.edu/sites/default/files/files/cdp1052.pdf
This question is linked to a previous question of mine that inquired about the first stage of this procedure - deriving and estimating the mixture of normal distribution for the observed measurements. Link:
How to calculate the sum of a "function" of a Gaussian mixture and a Gaussian variable?