Inverse Mills Ratio in Matlab

Question

I am trying to program in Matlab a conditional expectation of the form:

E[x|A<=x<=B] where X~N(u,s^2) (sorry, apparently the math editing here isn't what I am used to)

In Matlab, I have written up the following code:

Y=u+s*(normpdf(A,u,s)-normpdf(B,u,s))/(normcdf(B,u,s)-normcdf(A,u,s))

the problem is that it breaks down at higher values of A and B. For example, let u=0, s=1, A=10 and B=11. Simple logic says the answer should be between 10 and 11, but Matlab gives me back Inf because the denominator essentially becomes 0 while the numerator is 10^-23.

Any suggestions to make the formula give real numbers for all inputs?

Answer 1

One way is to do the numerical integration yourself:

x = linspace(A,B,1000);
trapz(x,x.*normpdf(x,u,s)) / trapz(x,normpdf(x,u,s))

With your example values, this gives 10.0981, and it is pretty fast