Defining a function

[John_Val]

hello,
I set up the following function, but when I try to print a value it comes up as NA

function JumpM n lambda k F X sigma delta gamma T r
type integer n X
type real JumpM lambda K F sigma delta gamma T r
local real bn d1n d2n

com JumpM=0
do i=0,n
com bn = -lambda*K + i*gamma/T
com d1n = ( log(F/X) + bn*T + 0.5*(sigma^2*T + i*delta^2) ) / (sigma^2*T + i*delta^2)^0.5
com d2n = d1n-( sigma^2*T + i*delta^2 )^0.5

com JumpM = JumpM + exp(-r*T)*( (exp(-lambda*T)*( lambda*T)^i)/%factorial(i) )*( F*exp(bn*T)*%CDF(d1n) - X*%CDF(d2n) )

end do
end Function

com a = JumpM(1000,5,0.05,34,20,0.3,0.02,0.04,1,0.02)

dis a
NA

[TomDoan]

The %factorial(i) in this will overflow when i is 1000

(exp(-lambda*T)*( lambda*T)^i)/%factorial(i)

If you look at the calculation in the JumpGARCH procedure in the Chan Maheu JBES 2002 procedure:

Code: Select all

function jumpgarch u h lambda deltasq theta
type real jumpgarch u h lambda deltasq theta
*
local   integer k
local   real    jsd wt jp
*
compute wt=0.0
compute jumpgarch=0.0
do k=0,kmax
   compute jsd=sqrt(h+k*deltasq)
   compute jp =exp(k*log(lambda)-%lngamma(k+1))
   compute wt =wt+jp
   compute jumpgarch=jumpgarch+jp*%density((u-k*theta)/jsd)/jsd
end do k
compute jumpgarch=log(jumpgarch/wt)
end jumpgarch

you'll notice how jp (which is the Poisson density) is computed by exping a log calculation. That's done to avoid overflows from the very large values that are cancelling out in the end.

[John_Val]

If I understand correctly,

exp(K*log(lambda)-%lngamma(K+1)) = (exp(-lambda)*lambda^k)/%factorial(k)

I tryed a few calculations for values of lambda and K, but they don't equal.

[TomDoan]

Sorry for the confusion. In the Poisson density, the exp(-lambda) is the normalizing constant, but it's the normalizing constant only if you sum to infinity. In order to use a standard Poisson density, you have to sum to a really large number to make the tail probability close to machine zero. That can be very expensive. What the JumpGARCH function does is to normalize by the actual sum of the lambda^j/j! over a more restricted length; in effect, using a Poisson conditional on an upper bound. That's what the calculations involving wt do. This is more stable numerically than using the truncated sum combined with the infinite normalizer. If you're going with the long sum route you would want

exp(-lambda+K*log(lambda)-%lngamma(K+1))

which will be the same as:

(exp(-lambda)*lambda^k)/%factorial(k)

[John_Val]

I am now trying to do a little simulation using this function, but it is not evaluating within a do loop.
Any calculations outside the loop work fine.

com draws=100
dec vec stock(draws+1) call(draws+1)
com dt=1.0/365
com stock(1)=20
com X=20
com mu=0.01
com sigma=0.9
com r=0.03

do i=1,draws
com stock(i+1) = stock(i)*exp((mu-0.5*sigma^2)*dt + sqrt(dt)*sigma*%ran(1))
com Term = (float(draws-i)/365)

com Call(i) = jumpM(300,5,0.05,stock(i),X,sigma,0.02,0.04,Term,r)

end do
## MAT13. Store into Out-of-Range Matrix or Series Element
The Error Occurred At Location 0043 of loop/block
Line 3 of loop/block

dis jumpM(300,5,0.05,price(1),X,sigma,0.02,0.04,1,r)
6.19998

At other times it would produce NA for all values

[TomDoan]

The variable "i" in your function isn't declared local so it ends up conflicting with the i that you're using in the loop. Add

local integer i

to your function definition and it should be fine.