When estimating VAR(1) by OLS, I get transpose of AR matrix. Why? Model: $y_t = \Phi y_{t-1} + (\mathbb{I_n}-\Phi)\mu + \Sigma \varepsilon_t$
s.t. $y_t$ is a $(n \times 1)$ vector, and $\varepsilon_t \sim N(\vec{0},\mathbb{I}_n)$. $\mathbb{I}$ is the identity matrix. $\Phi$ and $\Sigma$ are $(n \times n)$ matrices.

Consider the following MATLAB code.

Sigma = eye(2); %sqrtm(Covariance matrix)
Phi{1} = [0.9,0.12;0,0.9]; %Autoregressive matrix
mu = [0;0]; %Constant vector
nTs = 1000; %Length of simulation

%Specify Model
mdl = varm(2,1); 
mdl.AR = Phi;
mdl.Constant = mu;
mdl.Covariance = Sigma^2;

%Simulate Model
simTs = simulate(mdl,nTs); %Simulate VAR(1) (nTs x 2)

%Estimate Model
X = simTs(1:end-1,:);
Y = simTs(2:end,:); 
hatPhi = (X'*X)\(X'*Y)

&gt;&gt; hatPhi =
       0.8939  -0.0057
       0.1721   0.8732


Notice the output is the transpose of the AR matrix supplied. Why is this? I feel dumb, lol.
I see that $n=2$ and $\mathbf \mu = \mathbb0$. 
At single (2-dimensional) observation level, we usually write (ignoring the disturbance vector)

$$\mathbf y_t = \Phi \mathbf y_{t,-1}$$

where dimensions are for your case $(2 \times 1) = (2 \times 2) \times (2 \times 1)$.

But notice that the standard regression formulation has the coefficient vector to be "to the right" of the regressor matrix... so in order to be able to apply the standard OLS estimation formulas, we transpose and obtain

$$\mathbf y'_t = \mathbf y_{t,-1}'\Phi'$$

and at full sample level we get

$$\mathbf Y' = \mathbf Y_{-1}'\Phi'$$

where with $T+1$ 2-dimensional observations available, the dimensions are ($T \times 2$) = $(T \times 2) \times  (2 \times 2)$

So, when applying the standard OLS formulas as your code does, the estimator estimates 
$\Phi'$, not $\Phi$.