When estimating VAR(1) by OLS, I get transpose of AR matrix. Why?
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$.