I found a problem in using the integral function to compute a partial expectation for a lognormal distribution. Here is the code:
mu = 2.5652;
sig = 0.0055;
int = @(x) x.*lognpdf(x,mu,sig);
xx = (0.00001:0.1:20)';
f1 = zeros(length(xx),1);
f2 = zeros(length(xx),1);
for i = 1:length(xx)
f1(i) = integral(int,xx(i),Inf);
f2(i) = exp(mu+0.5*sig^2)*logncdf((mu+sig^2-log(xx(i)))/sig);
end
plot(xx,f1)
title('Using the integral function');
plot(xx,f2);
title('Using the formula for partial expectation');
Note that for that particular parametrization (with low variance), the partial expetation computed using the integral function drops abruptly to zero for some values of the low bound of the integral. This should no occur. In fact, when computing the partial expetation using the algebraic version (available for the lognomal) that problem does not ocurr. Is there a bug in the integral function?
