We consider n i.i.d. observations from the structural equation model. We denote by $X := (X_1,...,X_p)$ the $n \times p$ data matrix with n i.i.d. rows, each of them being $N(0, \Sigma)$-distributed, where $\Sigma$ is a nonsingular covariance matrix. The relations between the variables in a row can be represented as
$$X=XB + E,$$
where $B:= (\beta_{k,j})$ is a $p \times p$ matrix with $\beta_{j,j} =0$ for all j, and where E as an $n \times p$ matrix of noise vectors $E:=(\epsilon_1,...,\epsilon_p)$, with $\epsilon_j$ independent of $X_k$ whenever $\beta_{k,j} \neq 0$. Furthermore, E has n i.i.d rows which are $N(0, \Omega)$-distributed, with $\Omega := \text{diag} (|\omega_1|^2,...,|\omega_p|^2 ) $ a $p \times p$ diagonal matrix. The model implies that
$$\Sigma = [(I-B)^{-1}]^T \Omega [(I-B)^{-1}]$$ How to get this implication?
It's from a paper, but it doesn't have any steps that how to get the implication.
The paper is this one: https://projecteuclid.org/journals/annals-of-statistics/volume-41/issue-2/ell_0-penalized-maximum-likelihood-for-sparse-directed-acyclic-graphs/10.1214/13-AOS1085.full
You have $X = E(I-B)^{-1}$. It actually suffices to look at the first row of both sides. We have $$\begin{bmatrix} X_{1,1}, \ldots, X_{1, p} \end{bmatrix} = \begin{bmatrix} E_{1,1}, \ldots, E_{1, p} \end{bmatrix} (I-B)^{-1}.$$
We know the covariance matrix of $\begin{bmatrix} X_{1,1}, \ldots, X_{1, p} \end{bmatrix}$ and of $\begin{bmatrix} E_{1,1}, \ldots, E_{1, p}\end{bmatrix}$ are $\Sigma$ and $\Omega$ respectively.
Using the formula for the covariance matrix of a linear transformation of a random vector, we have $$\Sigma = ((I-B)^{-1})^\top \Omega (I-B)^{-1}.$$