Code:
* Example generated by -dataex-. To install: ssc install dataex
clear
input byte(y1 y2) float(dum1 dum2 y1F y2F)
0 1 0 1 .         .
3 0 1 0 .         .
0 2 0 1 .         .
1 0 0 1 .         .
2 2 1 0 .         .
3 1 0 1 .         .
2 0 1 0 .         .
4 1 1 0 .         .
5 3 1 0 .         .
2 0 1 0 .         .
1 2 0 1 .         .
8 1 1 0 .         .
1 1 0 1 .         .
5 1 1 0 .         .
3 3 0 1 .         .
4 0 1 0 .         .
. . 1 1 5 1.2698413
end

I estimate
gsem (y1 y2 <-dum1 dum2), poisson noconstant
then I obtain the prediction for the last observation (Note in older Satata versions I think following predict command gives an error)
predict y1F y2F if _n==_N
keeping in mind I need to transform coefficients to counts, now I want to know
1) probability of y1F -y2F for the last observation being greater than 0
I thought of using nlcom exp([y1]dum1 + [y1]dum2) - exp([y2]dum1 + [y2]dum2)
which produces:
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_nl_1 | 3.730159 1.963384 1.90 0.057 -.1180025 7.57832

so 3.730159 is the difference in predicted counts for obs _N but I am not sure if this is the way to go. I also looked into predictnl but with two poisson equations could not figure it out.
2) probability of y1F -y2F for the last observation being less than -1
3) probability of y1F -y2F for the last observation being less than +1

Any help would be appreciated.