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

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;
cov\[LetterSpace]1\[LetterSpace]2 = ComputeCovariance;
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]

• I'm not quite following exactly what you want but changing the last line in your code to m1\[LetterSpace]1\[LetterSpace]1 = Mean[xform\[LetterSpace]d\[LetterSpace]1] gets you the first moment in 0.01 seconds. The Mean, Variance, Skewness, and Kurtosis functions can be used to obtain the central moments pretty quickly (0.01, 0.06, 0.3, and 1.5 seconds, respectively). You could then get the raw moments from those quantities. I would also suggest using indexed variables (\[Mu] = {\[Mu][1, 1], \[Mu][1, 2]};, for example) rather than superscripts or subscripts.
– JimB
Dec 23, 2022 at 5:56