p:=0.7;q:=0.4;
P=(({{(1-p)+p q, p(1-q), 0, 0}, {(1-p)q, (1-p)(1-q) + p q, p (1-q), 0}, {0, (1-p)q, (1-p)(1-q) + p q, p (1-q)}, {0, 0, q, 1-q}}));

nextState[s_]:=Module[{},
sNew=s;
If[sNew< 3, sNew=sNew+RandomInteger[BernoulliDistribution[p]]];
If[sNew>0,sNew=sNew-RandomInteger[BernoulliDistribution[q]]];
Return[sNew];
];

traj=NestList[nextState,0,10];

ListPlot[traj,DataRange->{0,10},PlotStyle->PointSize[Large],Ticks->{Automatic,{0,1,2,3}}]

((Table[Count[traj,i],{i,0,3}])/1000000.)