Computing First Moment of TransformedDistribution slow. Is there a way to speed things up? If not, can you help me understand why its slow?

Question

I'm a novice Mathematica program, toying with the idea of solving for the first four moments of a linear combination of a multivariate Gaussian mixture distribution. The first two moments I have solved by hand, but the 3rd and 4th moments I believe to be better suited for something like Mathematics given there are hundreds of terms to collect.

The problem is as follows:

Let $s_t \in\{1,2\}$ be a finite-state Markov chain whose transition matrix is $\Pi$ whose rows sum to unity. The $(i,j)^{th}$ element of $\Pi$, $\pi_{i,j} := Pr(s_t = j | s_{t-1} = i)$.

Let the dynamics of $x_t \in \mathbb{R}^2$ by parameterized as a Vector Auto-regression subject to switches:

$$x_t = \mu_{s_t} + \Phi_{s_t}x_{t-1} + \Sigma_{s_t} \varepsilon_{t}$$

such that $\varepsilon_t \sim N(\vec{0}_{2 \times 1}, \mathbb{I}_2)$.

Let $\delta \in \mathbb{R}^2$. Then the goal is the compute the $n$-step ahead conditional expectation of the sum

$$\xi_{t,n} = \delta(x_t + x_{t+1} + \cdots + x_{t+n-1}),$$

and the corresponding coefficients.

$$m_n^{(1)}(s,x) := \mathbb{E}[\xi_{t,n} | s_t = s, x_t = x] = A_n(s) + B_n(s)^\prime x$$

The solution is the recurrence relationship $A_1(s) = 0$, $ B_1(s) = \delta $, $$ A_{n+1}(s) = \sum_{j=1}^2 \pi_{s,j} \big[ A_n(j) + \mu_j^\prime B_n(j)\big]$$

and

$$ B_{n+1}(s) = \delta + \sum_{j=1}^2 \pi_{s,j} \Phi_j^\prime B_n(j) $$

Here is the code for $n=2$, and it is taking multiple hours to run. I'm wondering if there is a way to speed it up, and if the higher-order moments are feasible? Any other tips on my code, and how to eventually collect coefficients A and B systematically would also be helpful.

Π = {{Subscript[π, 1, 1], 
    1 - Subscript[π, 1, 1]}, {1 - Subscript[π, 2, 2], 
    Subscript[π, 2, 2]}};
Subscript[μ, 1] = {Subscript[μ, 1, 1], Subscript[μ, 1, 2]};
Subscript[μ, 2] = {Subscript[μ, 2, 1], Subscript[μ, 2, 2]};
Subscript[Φ, 
   1] = {{Subscript[Φ, 1, 1, 1], 
    Subscript[Φ, 1, 1, 2]}, {Subscript[Φ, 1, 
     2, 1], Subscript[Φ, 1, 2, 2]}};
Subscript[Φ, 
   2] = {{Subscript[Φ, 2, 1, 1], 
    Subscript[Φ, 2, 1, 2]}, {Subscript[Φ, 2, 
     2, 1], Subscript[Φ, 2, 2, 2]}};
Subscript[Ω, 
   1] = {{Subscript[Ω, 1, 1, 1], 
    Subscript[Ω, 1, 1, 2]}, {Subscript[Ω, 
     1, 1, 2], Subscript[Ω, 1, 2, 2]}};
Subscript[Ω, 
   2] = {{Subscript[Ω, 2, 1, 1], 
    Subscript[Ω, 2, 1, 2]}, {Subscript[Ω, 
     2, 1, 2], Subscript[Ω, 2, 2, 2]}};
δ = {Subscript[d, 1], Subscript[d, 2]};
x0 = {x1, x2};
s0 = 1;

ComputeMean[s_Integer, x_] := 
 Subscript[μ, s] + Subscript[Φ, s] . x
ComputeCovariance[s_Integer] := Subscript[Ω, s]
mean\[LetterSpace]1\[LetterSpace]1 = 
 ComputeMean[1, x0] + x0
mean\[LetterSpace]1\[LetterSpace]2 = 
  ComputeMean[2, x0] + x0;
cov\[LetterSpace]1\[LetterSpace]1 = ComputeCovariance[1];
cov\[LetterSpace]1\[LetterSpace]2 = ComputeCovariance[2];
d\[LetterSpace]1\[LetterSpace]1 = 
 MultinormalDistribution[mean\[LetterSpace]1\[LetterSpace]1, 
  cov\[LetterSpace]1\[LetterSpace]1]
d\[LetterSpace]1\[LetterSpace]2 = 
 MultinormalDistribution[mean\[LetterSpace]1\[LetterSpace]2, 
  cov\[LetterSpace]1\[LetterSpace]2]
mix\[LetterSpace]1 = 
 MixtureDistribution[{Π[[1, 1]], Π[[1, 
    2]]}, {d\[LetterSpace]1\[LetterSpace]1, 
   d\[LetterSpace]1\[LetterSpace]2}]
xform\[LetterSpace]d\[LetterSpace]1 = 
 TransformedDistribution[δ . {u1, u2}, {u1, u2} \[Distributed] 
   mix\[LetterSpace]1]
m1\[LetterSpace]1\[LetterSpace]1 = 
 Moment[xform\[LetterSpace]d\[LetterSpace]1, 1]