Page 1 of 1
Defining a function
Posted: Tue Jun 01, 2010 4:13 pm
by 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
Re: Defining a function
Posted: Wed Jun 02, 2010 3:21 am
by 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.
Re: Defining a function
Posted: Wed Jun 02, 2010 1:24 pm
by 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.
Re: Defining a function
Posted: Wed Jun 02, 2010 1:51 pm
by 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)
Re: Defining a function
Posted: Wed Jun 02, 2010 3:20 pm
by 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
Re: Defining a function
Posted: Wed Jun 02, 2010 3:38 pm
by 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.