Mixture of normals, dependent on past I have the following probability model:
$(X_k|\text{PastHistory}_{k-1}, \theta_0,\theta_1,\theta_2) \sim (\pi\cdot N(\theta_1+\theta_0\cdot X_{k-1},1)+(1-\pi)\cdot N(\theta_2+\theta_0\cdot X_{k-1},1))$
and $X_k=0.5X_{k-1}+\epsilon_k$ with $\epsilon_k \sim N(0,1)$
I want to build a function that gives a sample of size n, from the mixture.
I've tried something like this in Mathematica:
  Sample[n_] := Module[{lst},
    lst = {RandomVariate[NormalDistribution[0, 1/(1 - 0.5^2)]]};
    element = 1;
    While[element <= n,
      AppendTo[lst, 
        RandomVariate[
          TransformedDistribution[
           p*X1 + (1 - p) X2, {
           X1 \[Distributed] 
           NormalDistribution[ θ1 + θ0*lst[[element]], 1], 
           X2 \[Distributed] 
           NormalDistribution[ θ2 + θ0*lst[[element]], 1]
        }]]];
     element = element + 1;
    ];
   lst
   ];

Is this correct?
Many thanks for any help. ;)
 A: We can write
$$X_k = \theta_0 X_{k-1} + \theta_1 U_{k-1} + \theta_2 (1 - U_{k-1}) + Z_{k-1}$$
where $U_i$ are iid Bernoulli$(\pi)$ variables and $Z_i$ are iid Normal$(0,1)$ variables.  This form of expression allows us to capitalize on vectorized methods of generating large lists of random variates at once, thereby obtaining $(U_i)$ and $(Z_i)$ quickly and from them the vector $Y = \theta_1 U + \theta_2 (1-U) + Z.$  All that remains is to compute the $X_k$ iteratively, which in Mathematica-speak is a "fold" of the bivariate function
$$f(X, Y) = \theta_0 X + Y.$$
That leads directly to code like the following, whose arguments are the length of the simulated process $n$, the vector of parameters $(\theta_0, \theta_1, \theta_2),$ the Bernoulli parameter $\pi$ (here written $p$), and the starting value $x_0$:
simulate[n_, {θ0_, θ1_, θ2_}, p_: 0.5, x0_: 0] := 
 With[
  {u = RandomVariate[BernoulliDistribution[p], n],
   z = RandomVariate[NormalDistribution[0, 1], n]}, 
  FoldList[θ0 #1 + #2 &, x0, θ1 u + θ2 (1 - u) + z]
  ]

To test it, consider the process with no autocorrelation and widely separated means in the mixture:
n = 10^5;
AbsoluteTiming[sim = simulate[n, {0, -3, 3}, 1/3];]


{0.100006, Null}

That's a million values per second--not too bad.  Let's look at some of these values; I plot every hundredth one in order to keep the plot visually recognizable:
ListPlot[sim[[1 ;; n ;; 100]], DataRange -> {1, n}]


Histogram[sim]


As intended, one-third ($p$) of the values are near $\theta_1=-3$ and the remaining two-thirds are near $\theta_2=3.$
