Expectation of a composite Markov-Gamma distribution

Question

In a model I have a discrete two-state first order Markov process, defined by a (2x2) transition matrix with two free parameters. If the first state occurs then the process outputs zero for that period. If the second state occurs then the process outputs a Gamma-distributed random variable characterised by another two parameters.

I can work out the theoretical mean and variance of the overall stationary process using pen and paper, but I am unsure how to check my results (and also I want to practice my Mathematica) using symbolic calculations. Put another way, I would like to know how to calculate the mean and variance to which the composite Markov-Gamma process will converge. My code is below - all I have currently is a way to calculate the mean of the Markov process, P, and its stationary mean (Mean[S]). However, I would like to be able to derive the mean and variance of R inserted into P for times when the process is in state 2, and zeros for when P is in state 1. Many thanks! Ben

M = {{1 - pdw, pdw}, {pwd, 1 - pwd}}
P  = DiscreteMarkovProcess[{1, 0}, M]
S = StationaryDistribution[P]
R = GammaDistribution[a,b]

Answer 1

With

$Assumptions = {0 <= pdw <= 1, 0 <= pwd <= 1};
trM = {{1 - pdw, pdw}, {pwd, 1 - pwd}};
dmP = DiscreteMarkovProcess[{1, 0}, trM];
gD = GammaDistribution[a, b];

the mean and variance of the described process for a time slice at time t are

mean = Expectation[(m[t] - 1)*g, {m \[Distributed] dmP, g \[Distributed] gD}] 
variance = Expectation[((m[t] - 1)*g - mean)^2, {m \[Distributed] dmP, g \[Distributed] gD}] // Simplify

-((a b pdw (-1 + (1 - pdw - pwd)^t))/(pdw + pwd))
-((a b^2 pdw (-1 + (1 - pdw - pwd)^t) (pdw + a pdw (1 - pdw - pwd)^t + pwd + a pwd))/(pdw + pwd)^2)

The stationary mean and variance are

stationaryMean = Limit[mean, t -> ∞]
stationaryVariance = Limit[variance, t -> ∞]

(a b pdw)/(pdw + pwd)
(a b^2 pdw (pdw + pwd + a pwd))/(pdw + pwd)^2

---

An alternative approach is to define a TransformedDistribution

trD = TransformedDistribution[(m - 1)* g, {m \[Distributed] SliceDistribution[dmP, t], 
       g \[Distributed] gD}]

and then calculate the stationary mean and variance with

Mean@trD ~Limit~ (t -> ∞)
Variance@trD ~Limit~ (t -> ∞)

---

---

Performing simulations:
One can use the TransformedDistribution to perform simulations:

rvSim = RandomVariate[trD /. {pdw -> 0.2, pwd -> 0.7, a -> 1, b -> 2, t -> #}] & /@ Range[1000];

Nevertheless, the calculation time is much lower, if one first uses RandomFunction with the DiscreteMarkovProcess and than does the state dependent transformation:

rf = RandomFunction[dmP /. {pdw -> 0.2, pwd -> 0.7}, {1, 10000}];

f[x_Integer] := 0 /; x == 1
f[x_Integer] := RandomVariate[GammaDistribution[1, 2]] /; x == 2

sim = f /@ (rf["Path"][[All, 2]]);

Answer 2

Following closely the documentation for DiscreteMarkovProcess (I am no expert though…)

M = {{1 - w, w}, {d, 1 - d}} P = DiscreteMarkovProcess[{1, 0}, M]

Mean[P[n]] // FullSimplify

Variance[P[n]] // FullSimplify

Graph[P]

data = RandomFunction[P /. w -> 1/2 /. d -> 1/4, {0, 10}]

ListPlot[data, Filling -> Axis, Ticks -> {Automatic, {1, 2, 3}}]

MarkovProcessProperties[P /. w -> 1/2 /. d -> 1/4]

Now for the stationary process

S = StationaryDistribution[P] // FullSimplify;

PDF[S, n]

Mean[S] // FullSimplify

Variance[S] // FullSimplify