Page 471 - A First Course In Stochastic Models

466 APPENDICES

∗
where f is defined by
p
n
b + iλ j + iπ(p/m − 1)
∗ ∗
f = α j Re f , p = 0, 1, . . . , 2m.
p
j=1
bℓ
The approximation of (2e / )g ℓ to f (ℓ ) is extraordinarily accurate. Rather
than calculating from (F.3) the constants g ℓ for ℓ = 0, 1, . . . , M − 1 by direct
summation, it is much better to use the discrete Fast Fourier Transform method to
calculate the constants g ℓ for ℓ = 0, 1, . . . , 2m − 1. More important than speeding
up the calculations, the discrete FFT method has the advantage of its numerical
stability. To see how to apply the discrete FFT method to (F.3), define % g k by
1 ∗ ∗

2 (f + f 2m ), k = 0,
0
% g k =
f , k = 1, . . . , 2m − 1.
∗
k
Then, we can rewrite the expression (F.3) for g ℓ as
& 2m−1 '
1
g ℓ ≈ Re % g k e 2πiℓk/2m (F.4)
2m
k=0
for ℓ = 0, 1, . . . , 2m − 1. The discrete FFT method can be applied to this repre-
sentation. Applying the inverse discrete FFT method to the vector (% g 0 , . . . ,% g 2m−1 )
yields the sought vector (g 0 , . . . , g 2m−1 ). Here is a summary of the algorithm:
Input: M, , b, n and m.
Output: f (k ) for k = 0, 1, . . . , M − 1.
Step 1: Calculate for p = 0, 1, . . . , 2m and 1 ≤ j ≤ n,

b + iλ j + iπ(p/m − 1)
∗
f jp = Re f ∗ .

n 1
α
∗
Next calculate f p ∗ = j=1 j f jp for p = 0, 1, . . . , 2m. Let % g 0 = 2 (f + f 2m )
∗
∗
0
∗
and % g k = f for k = 1, . . . , 2m − 1.
k
Step 2: Apply the inverse discrete FFT method to the vector (% g 0 , . . . ,% g 2m−1 ) in
order to obtain the desired vector (g 0 , . . . , g 2m−1 ).
bℓ
Step 3: Let f (ℓ ) = (2e / )g ℓ for 0 ≤ ℓ ≤ M − 1.
bℓ
In step 3 of the algorithm g ℓ is multiplied by e . In order to avoid numerical
instability, it is important to choose b not too large. Assuming that the ratio m/M
is large enough, say 4, numerical experiments indicate that b = 22/m gives results
that are almost of machine accuracy 2E − 16 (in general, it is best to choose
b somewhat larger than − ln(ξ)/(2m) where ξ is the machine precision). If f
is sufficiently smooth, it usually suffices to take n = 8, otherwise n = 16 is