both of your integrals can be done analytically with appropriate assumptions, note the way I formulated things both integrals are in non-delayed defintions.

c[eloc_] = 1/2*eloc + 1;
d[y_, x_] = x*y^2;
f[e_ /; e > 0, y_] = Assuming[{0 < x < 1, e > 0},
   Simplify[
    Integrate[Simplify[PDF[BetaDistribution[2, e], x]*d[y, x] + c[e]],
     {x, 0, 1}]]];
efun[y_?NumericQ] := e /. NMinimize[{f[e, y], e >= 0}, e][[2]]
(* this is the value of e that minimizes f *)
U[z_ /; z > 0, y_, p_] = p - Assuming[{0 < x < 1, z > 0},
   Integrate[
    PDF[BetaDistribution[2, z], x]*d[y, x] + c[z], {x, 0, 1}]]

ContourPlot[ U[efun[y], y, p] == 0, {y, 0, 2}, {p, 1, 3}]

I'd suggest you review this and make sure its right then work on the NMaximize step.

both of your integrals can be done analytically with appropriate assumptions, note the way I formulated things both integrals are in non-delayed defintions.

c[eloc_] = 1/2*eloc + 1;
d[y_, x_] = x*y^2;
f[e_ /; e > 0, y_] = Assuming[{0 < x < 1, e > 0},
   Simplify[
    Integrate[Simplify[PDF[BetaDistribution[2, e], x]*d[y, x] + c[e]],
     {x, 0, 1}]]];
efun[y_?NumericQ] := e /. NMinimize[{f[e, y], e >= 0}, e][[2]]
(* this is the value of e that minimizes f *)
U[z_ /; z > 0, y_, p_] = p - Assuming[{0 < x < 1, z > 0},
   Integrate[
    PDF[BetaDistribution[2, z], x]*d[y, x] + c[z], {x, 0, 1}]]

ContourPlot[ U[efun[y], y, p] == 0, {y, 0, 2}, {p, 1, 3}]

I'd suggest you review this and make sure its right..

FindMaximum[{y - p, U[efun[y], y, p] == 0, y > 0,   p > 0}, {{y, 1}, {p, 1}}]

{-0.75, {y -> 0.5, p -> 1.25}}