sup = (x > 0 && y > 0 && 1 > θ > 0)

gfxy = ((1 + θ x) (1 + θ y) - θ) 
        Exp[-(x + y + θ x y)]

gfx = 
      Integrate[gfxy, {y, 0, Infinity}, 
        Assumptions -> (sup)]

gfy = 
      Integrate[gfxy, {x, 0, Infinity}, 
        Assumptions -> (sup)]

gfxdy = (gfxy/gfy);
gfydx = (gfxy/gfx);

"E(X|Y=y)"
gexdy = 
      Integrate[x*gfxdy, {x, 0, Infinity}, 
        Assumptions -> (sup)]
"E(Y|X=x)"
geydx = 
      Integrate[y*gfydx, {y, 0, Infinity}, 
        Assumptions -> (sup)]

gexxy = 
      Integrate[(x^2) (gfxy/gfy), {x, 0, Infinity}, 
        Assumptions -> (sup)]
gexy = 
      Integrate[x  (gfxy/gfy), {x, 0, Infinity}, 
        Assumptions -> (sup)]
"Var(X|Y=y)"
gvar_xy = Simplify[gexxy - (gexy^2)]

geyyx = 
      Integrate[(y^2) (gfxy/gfx), {y, 0, Infinity}, 
        Assumptions -> (sup)]
geyx = 
      Integrate[y  (gfxy/gfx), {y, 0, Infinity}, 
        Assumptions -> (sup)]
"Var(Y|X=x)"
gvar_yx = Simplify[geyyx - (geyx^2)]