how to make an infinite sum

Hi guys,

I have an infintite sum, I write the following code, but it doesn't compute the values.

function y=Phiplus(p,k,delta,alpha)
h=0;    
    for n=1:1000,
    f=2^p*exp(-alpha/2)*(alpha/2)^n*(gamma(k/2)-gamma(p+delta/2+n)*gammainc(p+delta/2+n,k/2))/(factorial(n)*gamma(delta/2+n));
    h=h+f;
    end
y=h;    
end

Pleas help! I use this function to calculate the price of a call option:

Call(i,j,m,n)=exp(-q(j)*T(i))*S*Phiplus(0,k(i,j,n)^2/tau(i,j), delta(m), x^2/tau(i,j))-exp(-(r(j)+b(j))*T(i))*K(n)*(x^2/tau(i,j))^(1/(2*abs(beta)))* Phiplus(-0.5/abs(beta),k(i,j,n)^2/tau(i,j),delta(m),x^2/tau(i,j));

If you need any additional info I'm ready to provide it.

[dpb]

Instead of using a fixed upper limit on the loop, consider using a conditional loop, stopping when the additional term no longer adds more than some preset tolerance (say 0.1% or so?) to the value.

factorial(1000) (and long before then) is going to be inf so at that point the term is going to be identically zero and since there are negative exponentials, they're going to zero as well and so I presume you're getting NaN as a result owing to adding 0/0 to the result.

The answer is to only compute the minimum number of terms you need and to probably eliminate the factorial function entirely but instead simply divide each time by the next N in the sequence.

I'd think you probably will need only a few terms at most to get the result but that will depend upon the magnitude of the various coefficients. If both numerator and denominator are small from the git-go, you may need to recast entirely, but try the obvious first.

[Daniel]

Great! Many thanks!

Answer 1

[Roger Stafford]

I'll expand on dpb's good advice. With regard to the factor contained in the terms you are adding which I will call F(n)

 F(n) = (alpha/2)^n/factorial(n)/gamma(delta/2+n)

you can obtain the next value using the iteration

 F(n+1) = F(n)*alpha/2/(n+1)/(delta/2+n)

If you carry out the increase of n by 1 at each step using this technique, you will not encounter the difficulties of a very large (alpha/2)^n divided by a very large factorial(n) and a large gamma(delta/2+n) which will all too soon give you "inf" values.

This way you can hopefully go sufficiently far out in your series so that they become small enough to legitimately halt your iteration without prematurely encountering inf values from the computer.

As for the factor involving the incomplete gamma function perhaps a similar strategy can be devised. Incidentally in case there is any confusion here, I will remind you that matlab's incomplete gamma function is both normalized and gives its arguments in the reverse order from the way it is often defined in mathematics.

[Daniel]

Thanks a lot! The idea with the F(n) is quite good, but I don't know how to simplify the gamma part, so that it would be just f(n+1)=f(n)*x for the whole function... Any ideas? In my equation gamma inverse is not normalized and the limits of the integral are zero and x.

[Roger Stafford]

My question is, in "gammainc(p+delta/2+n,k/2)" which argument is to be the parameter and which the limit of integration? Also how is it that both quantities appear as parameters in (complete) gamma functions elsewhere in your series? That seems a little curious.

[Daniel]

Yes, you are right, it should be the other way around. gammainc(x,a)=gammainc(k/2,p+delta/2+n)

[Roger Stafford]

If k/2 is the upper limit of integration in the incomplete gamma function, it strikes as curious that in your expression you have also called on the complete gamma function using k/2 as its parameter. Is that what you intended?

[Daniel]

The idea is that the gammainv function in Matlab is normalized, and I need an not normalized gamma inv