Great, thanks you very much! I did not know that creations of random number was very expensive.

Hi Alessio,

The main problem is that you create q everytime (in two loops). If you move q outside the inner Do loop it goes 100x faster:

ClearAll[RecurrenceDistributionRandomGraph,y]
RecurrenceDistributionRandomGraph[J_,H_,\[Beta]_,y_,Pf_,P0_,n_,m_]:=
    Block[{x0,q,qs,x,z},
        x0=RandomVariate[P0,n];
        R={x0};
        Do[
            x=ConstantArray[0.0,n];
            qs=Round[RandomVariate[Pf,n]];
            Do[
                (*q=Round[RandomVariate[Pf]];*)
                q=qs[[j]];
                z=RandomSample[R[[i]],q];
                x[[j]]=N[y[z,J,H,\[Beta]]]
                ,
                {j,1,n}
            ];
            R=Append[R,x];
            ,
            {i,1,m}
        ]
    ]
P0=PoissonDistribution[5];
Pf=ProbabilityDistribution[q*(E^-3 3^q/q!)/3,{q,0,+\[Infinity]}];
y[z_,J_,H_,\[Beta]_]:=(E^(\[Beta] (J+H)) Product[i,{i,z}]+E^(-\[Beta] (J+H)))/(E^(\[Beta] (-J+H)) Product[i,{i,z}]+E^(\[Beta] (J-H)))

executing:

RecurrenceDistributionRandomGraph[1., 1., 0.2, y, Pf, P0, 100, 20] // AbsoluteTiming

runs in 1.12 seconds on my laptop (down from 108 seconds). Note that I also used ConstantArray rather than Table, and I put 0.0 in there so the data-type does not change during evaluation (from integer to doubles).