The expression in question is (we have replaced sigma by s)
g = Integrate[ Exp[-(x - y)^2/(2 y s^2)] f[y] 1/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}] $\int_0^{\infty } \frac{e^{-\frac{(x-y)^2}{2 s^2 y}} f[y]}{\sqrt{2 \pi } \sqrt{s^2 y}} \, dy$
First of all we notice that for y>0
"Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], s -> 0] = DiracDelta[x - y]" Formally
Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], s -> 0, Assumptions -> {y > 0, x > 0, y != x}] (* Out[37]= 0 *) In[39]:= Limit[Exp[-(x - y)^2/(2 y s^2)] 1/Sqrt[2 \[Pi] y s^2], s -> 0, Assumptions -> {y > 0, x > 0, y == x}] (* Out[39]= \[Infinity] *) so that
Limit[g, s -> 0] = f[x] Assuming now that f[y] can be expanded into a power series about y = 0 and consider
g1 = Integrate[ Exp[-(x - y)^2/(2 y s^2)] y^k/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}, Assumptions -> {s > 0, x^2 > 0}] (* Out[16]= (E^(x/s^2) Sqrt[2/\[Pi]] Abs[x]^(1/2 + k) BesselK[1/2 + k, Abs[x]/s^2])/s *) Now the series expansion about s = 0 becomes
sg1 = Series[g1, {s, 0, 2}, Assumptions -> s > 0] // Normal (* Out[27]= E^(x/s^2 - Abs[x]/s^2) (1/2 k (1 + k) s^2 Abs[x]^(-1 + k) + Abs[x]^k) *) Let us distinguish the region x<0 and x>0:
Simplify[sg1, x < 0] (* Out[28]= 1/2 E^((2 x)/s^2) (k (1 + k) s^2 - 2 x) (-x)^(-1 + k) *) Obviously, the exponential factor cannot be expanded about s=0.
Simplify[sg1, x > 0] (* Out[31]= 1/2 x^(-1 + k) (k s^2 + k^2 s^2 + 2 x) *) In the limit
Limit[%, s -> 0] (* Out[32]= x^k *) i.e. we rediscover the effect of the delta function.
Edit (Till Hoffmann): I thought I'd add to this answer to complete it.
Note that the expansion $\mathtt{sg1}=x^k + \sigma^2 \frac{k(1+k)}{2} x^{k-1}$ can be related to the Taylor expansion of $$ f(x) + \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} \left(x f(x)\right), $$ which is the desired series expansion in terms of the original function.
EDIT (Dr. Wolfgang Hintze, 28.02.15)
Let me add the general expansion in the form of Till Hoffmann, which I call s-expansion in the following
$$fs = \sum _{j=0}^{\infty } \frac{s^{2 j} }{j!2^{j}}\frac{\partial ^{2 j}\left(f(x) x^j\right)}{\partial x^{2 j}}$$
It is very interesting that this expression seems to hold not only for the powers of y considered up to now but for more general functions f(y).
In order to verify this hypothesis we shall study examples for other functions than simple powers of the s-expansion fs of our integral.
The s-expansion with a function f is given by
fs[x_, s_, f_, n_] := Sum[s^(2 j) 1/(j! 2^j) D[x^(j) f, {x, 2 j}], {j, 0, n}] Example 1: Simple pole on the negative real y-axis
The integral is numerically
gg1[x_, s_] := NIntegrate[ Exp[-(x - y)^2/(2 y s^2)] 1/(1 + y) 1/Sqrt[2 \[Pi] y s^2], {y,0, \[Infinity]}] The first 3 terms of the s-expansion are
fs1[x_, s_] = Table[fs[x, s, 1/(1 + x), k], {k, 0, 2}] // Simplify (* {1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, ( 3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5} *) Graphs up to order s^4 for two values of s are
With[{s = 1}, Plot[{gg1[x, s], {1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, ( 3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5}}, {x, -1, 2}, PlotRange -> {0, 1.5}, ImageSize -> 400, PlotLabel -> Style["s-expansion of integral\nf(y) = 1/(1+y), s = " <> ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg1[x]"}, Epilog -> {Text[ Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = f(x)\n\ brown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = \ O(\!\(\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]] (* 150228_s-expansion_f1 _s1.jpg *) 
With[{s = 0.3}, Plot[{gg1[x, s], {1/(1 + x), (-s^2 + (1 + x)^2)/(1 + x)^3, ( 3 s^4 - s^2 (1 + x)^2 + (1 + x)^4)/(1 + x)^5}}, {x, -1, 2}, PlotRange -> {0, 1.5}, ImageSize -> 400, PlotLabel -> Style["s-expansion of integral\nf(y) = 1/(1+y), s = " <> ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg1[x]"}, Epilog -> {Text[ Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = \ f(x)\nbrown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = O(\!\(\ \*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]] (* 150228_s-expansion_f1 _s0-3.jpg *) 
Example 2: Simple conjugate poles on the imaginary y-axis
The integral is numerically
gg2[x_, s_] := NIntegrate[ Exp[-(x - y)^2/(2 y s^2)] 1/(1 + y^2) 1/Sqrt[2 \[Pi] y s^2], {y, 0, \[Infinity]}] The first 3 terms of the s-expansion are
fs2[x, s]
(* {1/(1 + x^2), (s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, ( s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5} *) Graphs up to order s^4 for two values of s are
With[{s = 1}, Plot[{gg2[x, s], {1/(1 + x^2), ( s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, ( s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5}}, {x, -1, 2}, PlotRange -> {0, 1.5}, ImageSize -> 400, PlotLabel -> Style["s-expansion of integral\nf(y) = 1/(1+\!\(\*SuperscriptBox[\(y\), \ \(2\)]\)), s = " <> ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg2[x]"}, Epilog -> {Text[ Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = f(x)\n\ brown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = \ O(\!\(\*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]] (* 150228_s-expansion_f2 _s0-9.jpg *) 
With[{s = 0.4}, Plot[{gg2[x, s], {1/(1 + x^2), ( s^2 x (-3 + x^2) + (1 + x^2)^2)/(1 + x^2)^3, ( s^2 x (-3 + x^2) (1 + x^2)^2 + (1 + x^2)^4 - 3 s^4 (1 - 10 x^2 + 5 x^4))/(1 + x^2)^5}}, {x, -1, 2}, PlotRange -> {0, 1.5}, ImageSize -> 400, PlotLabel -> Style["s-expansion of integral\nf(y) = \ 1/(1+\!\(\*SuperscriptBox[\(y\), \(2\)]\)), s = " <> ToString[s] <> "\n", 14], AxesLabel -> {"x", "gg2[x]"}, Epilog -> {Text[ Style["Legend of curves:\nblue = integral, numeric (gg1)\nred = \ f(x)\nbrown = O(\!\(\*SuperscriptBox[\(s\), \(2\)]\))\ngreen = O(\!\(\ \*SuperscriptBox[\(s\), \(4\)]\))", Medium], {1, 1.2}]}]] (* 150228_s-expansion_f2 _s0-4.jpg *) 
We see that the agreement is what is to expected in the sense of an asymptotic expansion.
Conclusion: the formula for $fs$ seems to be valid for a broader class of function.