Expectation of a composite Markov-Gamma distribution

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]

• things like PDF[P[n], k] // PiecewiseExpand or Mean[P[n]] // FullSimplify? – chris Aug 21 '14 at 17:09
• I did not understand the link with R? – chris Aug 21 '14 at 17:23
• R is a gamma distribution which determines the outcome of the series if it is in state '2'. Otherwise the process outputs zero in state '1'. – ben18785 Aug 21 '14 at 17:46

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]]);

• I am impressed you seem to have understood what the OP wanted based on limited evidence :-) – chris Aug 22 '14 at 7:44
• Thanks to Karsten for answering my question - that helps me understand how to approach this problem. – ben18785 Aug 22 '14 at 9:23
• @Karsten. Can I ask how I can simulate draws from the transformed distribution? I have tried RandomFunction (substituting in the required parameters with numbers), but this does not seem to work. – ben18785 Aug 22 '14 at 10:03
• @user2003951 your difficulties seem to be due to the circumstance that for Mathematica there is a bigger difference between distributions and random processes than for you. RandomFunction is for processes, not for distributions. For distributions you should use RandomVariate. I'll make an edit to my answer, in order to give you an example. – Karsten 7. Aug 22 '14 at 11:09
• @Karsten - that's great thank you. Yes, the latter option is much faster. Can I ask though, how do I actually plot the data from the simulation? I have tried ListPlot[sim] but that doesn't seem to work. Sorry! – ben18785 Aug 22 '14 at 13:13

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


• Hi, Thanks for your answer, although it is not quite what I am looking for. I can derive the properties of a simple Markov chain using Mathematica. What I can't do is form the composite process consisting of zeros if the process is in state 1, and a gamma distributed random variable in state 2. I don't know how to make this composite process using Mathematica, then symbollically derive its long run mean and variance. Thanks! – ben18785 Aug 21 '14 at 17:50
• ok I cannot help you there because I know too little about Markov processes. There are ways to compose probabilities in mathematica: may be someone else will help. – chris Aug 21 '14 at 18:41