1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
|
C
C Summary:
C
C SYMSTB
C Fortran program to compute a fast numerical approximation to the Symmetric
C Stable distribution and density functions.
C (Hu McCulloch, mcculloch.2@osu.edu)
C
C Submission:
C SYMSTB (VS FORTRAN)
C Numerical Approximation to Symmetric Stable Distribution and Density.
C Posted by: J. Huston McCulloch
C Economics Dept.
C Ohio State Univ.
C 1945 N. High St.
C Columbus, OH 43210
C (614) 292-0382
C mcculloch.2@osu.edu
C The code may be freely used for non-commercial purposes and freely copied
C and distributed. The computation is described in
C J. Huston McCulloch, "Numerical Approximation of the Symmetric Stable
C Distributions and Densities," Ohio State Univ. Economics Dept., Oct. 1994,
C which is available from the author on request. Any research use
C of this code should cite this working paper. Contact the author for
C assistance.
C
C Accuracy:
C The expected relative density precision is 1.0e-6 for alpha in the
C range [.84, 2.00]. The programs have considerably reduced precision for
C alpha < .84, although no error message is given for lower alpha
C values. The absolute precision of the density is 6.6e-5
C for alpha in the range [..92, 2.00], while that of the distribution is
C 2.2e-5 in the same range. See paper for details.
C
C Speed:
C Changing alpha induces set-up calculations, so submit as many x values
C as you can before changing alpha. Only the commands from 20 down are
C are executed if alpha is unchanged from the previous call. The GAUSS
C version of this program takes as little as 33 microseconds to compute
C the density (only) on a P5/100.
C
C Compatibility:
C DERFC is a built-in VS FORTRAN function that computes the complemented
C error function, which is related to the cumulative normal distribution.
C If your FORTRAN does not have this, but instead has a complemented
C cumulative normal distribution function named (eg) NCDFC, command #1
C may be replaced by
C 1 GFUN(X) = NCDFC(X/SQRT(2))
C
C A GAUSS version of this routine is archived at
C gopher.american.edu/academic.depts/cas/econ/software/gauss
C
C FORTRAN code follows:
C
C
C ******************************************************************************
C
SUBROUTINE SYMSTB(XX,YY,ZZ,NN,ALPHA)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XX(NN),YY(NN),ZZ(NN)
DIMENSION EI(3),U(3),S(4,19),R(19),Q(0:5),P(0:5,0:19)
DIMENSION PD(0:4,0:19),ALF2I(4)
DIMENSION ZNOT(19),ZN4(19),ZN5(19),COMBO(0:5),ZJI(19,0:5)
DATA COMBO/1.0D0,5.0D0,10.0D0,10.0D0,5.0D0,1.0D0/
C
DATA ((S(I, J), I = 1, 4), J = 1, 7) /
1 1.85141 90959 D2, -4.67693 32663 D2,
1 4.84247 20302 D2, -1.76391 53404 D2,
1 -3.02365 52164 D2, 7.63519 31975 D2,
1 -7.85603 42101 D2, 2.84263 13374 D2,
1 4.40789 23600 D2, -1.11811 38121 D3,
1 1.15483 11335 D3, -4.19696 66223 D2,
1 -5.24481 42165 D2, 1.32244 87717 D3,
1 -1.35556 48053 D3, 4.88340 79950 D2,
1 5.35304 35018 D2, -1.33745 70340 D3,
1 1.36601 40118 D3, -4.92860 99583 D2,
1 -4.89889 57866 D2, 1.20914 18165 D3,
1 -1.22858 72257 D3, 4.40631 74114 D2,
1 3.29055 28742 D2, -7.32117 67697 D2,
1 6.81836 41829 D2, -2.28242 91084 D2 /
DATA ((S(I, J), I = 1, 4), J = 8,14) /
1 -2.14954 02244 D2, 3.96949 06604 D2,
1 -3.36957 10692 D2, 1.09058 55709 D2,
1 2.11125 81866 D2, -2.79211 07017 D2,
1 1.17179 66020 D2, 3.43946 64342 D0,
1 -2.64867 98043 D2, 1.19990 93707 D2,
1 2.10448 41328 D2, -1.51108 81541 D2,
1 9.41057 84123 D2, -1.72219 88478 D3,
1 1.40875 44698 D3, -4.24725 11892 D2,
1 -2.19904 75933 D3, 4.26377 20422 D3,
1 -3.47239 81786 D3, 1.01743 73627 D3,
1 3.10474 90290 D3, -5.42042 10990 D3,
1 4.22210 52925 D3, -1.23459 71177 D3,
1 -5.14082 60668 D3, 1.10902 64364 D4,
1 -1.02703 37246 D4, 3.42434 49595 D3 /
DATA ((S(I, J), I = 1, 4), J = 15, 19) /
1 1.12151 57876 D4, -2.42435 29825 D4,
1 2.15360 57267 D4, -6.84909 96103 D3,
1 -1.81206 31586 D4, 3.14301 32257 D4,
1 -2.41642 85641 D4, 6.91268 62826 D3,
1 1.73884 13126 D4, -2.21083 97686 D4,
1 1.33979 99271 D4, -3.12466 11987 D3,
1 -7.24357 75303 D3, 4.35453 99418 D3,
1 2.36161 55949 D2, -7.65716 53073 D2,
1 -8.73767 25439 D3, 1.55108 52129 D4,
1 -1.37897 64138 D4, 4.63874 17712 D3 /
C
PI=3.141592653589793D0
CA=DGAMMA(ALPHA)*DSIN(PI*ALPHA/2.0D0)/PI
SQPI=DSQRT(PI)
A2=DSQRT(2.0D0)-1.0D0
CPXP0=1.0D0/PI
GPXP0=1.0D0/(4.0D0*A2*SQPI)
CPXPP0=CPXP0*2.0D0
GPXPP0=GPXP0*1.5D0
CPPP=CPXPP0*3.0D0-2.0D0/PI
GPPP=GPXPP0*2.5D0-1.0D0/(32.0D0*SQPI*A2**3)
DO J=1,19
ZNOT(J)=J*0.05D0
ZN4(J)=(1-ZNOT(J))**4
ZN5(J)=(1-ZNOT(J))*ZN4(J)
DO I=0,5
ZJI(J,I)=COMBO(I)*(-ZNOT(J))**(5-I)
ENDDO
ENDDO
DO I=1,3
EI(I)=I
U(I)=1
ENDDO
Q(0)=0.0D0
A=2.0D0**(1.0D0/ALPHA)-1.0D0
ALA=ALPHA*A
ALF2=2.0D0-ALPHA
ALF1=ALPHA-1.0D0
PIALF=PI*ALPHA
SP0=DGAMMA(1.0D0/ALPHA)/PIALF
SPPP0=-DGAMMA(3.0D0/ALPHA)/PIALF
XP0=1.0D0/(ALA)
XPP0=XP0*(1.0D0+ALPHA)/ALPHA
XPPP0=XPP0*(1.0D0+2.0D0*ALPHA)/ALPHA
SPZP1=(A**ALPHA)*DGAMMA(ALPHA)*DSIN(PIALF/2.0D0)/PI
RP0=-SP0*XP0+ALF2*CPXP0+ALF1*GPXP0
RPP0=-SP0*XPP0+ALF2*CPXPP0+ALF1*GPXPP0
RPPP0=-SP0*XPPP0-SPPP0*XP0**3+ALF2*CPPP+ALF1*GPPP
RP1=-SPZP1+ALF2/PI
DO I=1,4
ALF2I(I)=ALF2**I-1.0D0
ENDDO
DO J=1,19
R(J)=ALF2*(ALF2I(1)*S(1,J)+ALF2I(2)*S(2,J)
& +ALF2I(3)*S(3,J)+ALF2I(4)*S(4,J))
ENDDO
Q(1)=RP0
Q(2)=RPP0/2.0D0
Q(3)=RPPP0/6.0D0
B=-(U(1)+U(2)+U(3))*Q(1)
DO IP=1,19
B=B-R(IP)*ZN5(IP)
ENDDO
C=RP1-(EI(1)+EI(2)+EI(3))*Q(1)
DO IP=1,19
C=C-5.0D0*R(IP)*ZN4(IP)
ENDDO
Q(4)=5.0D0*B-C
Q(5)=B-Q(4)
DO I=0,5
P(I,0)=Q(I)
DO J=1,19
PRD=0.0D0
DO IP=1,J
PRD=PRD+R(IP)
ENDDO
P(I,J)=Q(I)+PRD*ZJI(1,I)
ENDDO
ENDDO
DO I=1,5
DO J=0,19
PD(I-1,J)=I*P(I,J)
ENDDO
ENDDO
C LOOP OVER ALL DATAPOINTS:
DO II=1,NN
X=XX(II)
XA1=1.0D0+A*DABS(X)
XA1A=XA1**(-ALPHA)
Z=1.0D0-XA1A
ZP=ALA*XA1A/XA1
X1=((1.0D0-Z)**(-1.0D0)-1.0D0)
X2=((1.0D0-Z)**(-0.5D0)-1.0D0)/A2
X1P=1.0D0/((1.0D0+X1)**(-2.0D0))
X2P=1.0D0/(2.0D0*A2*(1.0D0+A2*X2)**(-3.0D0))
J=20*Z
IF (J.GT.19) J=19
RZ=P(0,J)
DO IP=4,0,-1
RZ=RZ*Z+P(0,J)
ENDDO
RPZ=PD(0,J)
DO IP=3,0,-1
RPZ=RPZ*Z+PD(0,J)
ENDDO
C
C *** CUMULATED PROBABILITY FUNCTION:
CFUN=0.5D0-DATAN(X1)/PI
GFUN=0.5D0*DERFC2(X2/2.0D0)
PROBFUN=ALF2*CFUN+ALF1*GFUN+RZ
IF(X.LT.0.0D0) PROBFUN=1.0D0-PROBFUN
YY(II)=1.0D0-PROBFUN
IF (YY(II).LT.10*2.2D-5) THEN
YY(II)=CA*DABS(X)**(-ALPHA)
ENDIF
C
C *** PROBABILITY DENSITY FUNCTION:
CDEN=1.0D0/(PI*(1.0D0+X1*X1))
GDEN=DEXP(-x2*X2/4.0D0)/(2.0D0*SQPI)
PROBDEN=(ALF2*CDEN*X1P+ALF1*GDEN*X2P-RPZ)*ZP
ZZ(II)=PROBDEN
IF (ZZ(II).LT.10*6.6D-5) THEN
ZZ(II)=ALPHA*CA*DABS(X)**(-ALPHA-1.0D0)
ENDIF
C
C END OF X(NN) LOOP:
ENDDO
RETURN
END
C
C ******************************************************************************
C
DOUBLE PRECISION FUNCTION DERFC2(X)
INTEGER JINT
DOUBLE PRECISION X, RESULT
JINT = 1
CALL CALERF(X,RESULT,JINT)
DERFC2 = RESULT
RETURN
END
C
SUBROUTINE CALERF(ARG,RESULT,JINT)
INTEGER I,JINT
DOUBLE PRECISION
1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI,
2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL,
3 Y,YSQ,ZERO
DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5)
DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/,
1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/,
2 SIXTEN/16.0D0/
DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/,
1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/
DATA A/3.16112374387056560D00,1.13864154151050156D02,
1 3.77485237685302021D02,3.20937758913846947D03,
2 1.85777706184603153D-1/
DATA B/2.36012909523441209D01,2.44024637934444173D02,
1 1.28261652607737228D03,2.84423683343917062D03/
DATA C/5.64188496988670089D-1,8.88314979438837594D0,
1 6.61191906371416295D01,2.98635138197400131D02,
2 8.81952221241769090D02,1.71204761263407058D03,
3 2.05107837782607147D03,1.23033935479799725D03,
4 2.15311535474403846D-8/
DATA D/1.57449261107098347D01,1.17693950891312499D02,
1 5.37181101862009858D02,1.62138957456669019D03,
2 3.29079923573345963D03,4.36261909014324716D03,
3 3.43936767414372164D03,1.23033935480374942D03/
DATA P/3.05326634961232344D-1,3.60344899949804439D-1,
1 1.25781726111229246D-1,1.60837851487422766D-2,
2 6.58749161529837803D-4,1.63153871373020978D-2/
DATA Q/2.56852019228982242D00,1.87295284992346047D00,
1 5.27905102951428412D-1,6.05183413124413191D-2,
2 2.33520497626869185D-3/
X = ARG
Y = DABS(X)
IF (Y .LE. THRESH) THEN
YSQ = ZERO
IF (Y .GT. XSMALL) YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO 20 I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
20 CONTINUE
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = DEXP(YSQ) * RESULT
GO TO 800
ELSE IF (Y .LE. FOUR) THEN
XNUM = C(9)*Y
XDEN = Y
DO 120 I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
120 CONTINUE
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT .NE. 2) THEN
YSQ = DINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = DEXP(-YSQ*YSQ) * DEXP(-DEL) * RESULT
END IF
ELSE
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 300
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO 240 I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
240 CONTINUE
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT .NE. 2) THEN
YSQ = DINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = DEXP(-YSQ*YSQ) * DEXP(-DEL) * RESULT
END IF
END IF
300 IF (JINT .EQ. 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X .LT. ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X .LT. ZERO) RESULT = TWO - RESULT
ELSE
IF (X .LT. ZERO) THEN
IF (X .LT. XNEG) THEN
RESULT = XINF
ELSE
YSQ = DINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = DEXP(YSQ*YSQ) * DEXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
800 RETURN
END
C
C*******************************************************************************
C
DOUBLE PRECISION FUNCTION DGAMMA(X)
INTEGER I,N
LOGICAL PARITY
DOUBLE PRECISION
1 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE,
2 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
DIMENSION C(7),P(8),Q(8)
DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
1 SQRTPI/0.9189385332046727417803297D0/,
2 PI/3.1415926535897932384626434D0/
DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
1 XINF/1.79D308/
DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
DATA C/-1.910444077728D-03,8.4171387781295D-04,
1 -5.952379913043012D-04,7.93650793500350248D-04,
2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
3 5.7083835261D-03/
CONV(I) = DBLE(I)
PARITY = .FALSE.
FACT = ONE
N = 0
Y = X
IF (Y .LE. ZERO) THEN
Y = -X
Y1 = DINT(Y)
RES = Y - Y1
IF (RES .NE. ZERO) THEN
IF (Y1 .NE. DINT(Y1*HALF)*TWO) PARITY = .TRUE.
FACT = -PI / DSIN(PI*RES)
Y = Y + ONE
ELSE
RES = XINF
GO TO 900
END IF
END IF
IF (Y .LT. EPS) THEN
IF (Y .GE. XMININ) THEN
RES = ONE / Y
ELSE
RES = XINF
GO TO 900
END IF
ELSE IF (Y .LT. TWELVE) THEN
Y1 = Y
IF (Y .LT. ONE) THEN
Z = Y
Y = Y + ONE
ELSE
N = INT(Y) - 1
Y = Y - CONV(N)
Z = Y - ONE
END IF
XNUM = ZERO
XDEN = ONE
DO 260 I = 1, 8
XNUM = (XNUM + P(I)) * Z
XDEN = XDEN * Z + Q(I)
260 CONTINUE
RES = XNUM / XDEN + ONE
IF (Y1 .LT. Y) THEN
RES = RES / Y1
ELSE IF (Y1 .GT. Y) THEN
DO 290 I = 1, N
RES = RES * Y
Y = Y + ONE
290 CONTINUE
END IF
ELSE
IF (Y .LE. XBIG) THEN
YSQ = Y * Y
SUM = C(7)
DO 350 I = 1, 6
SUM = SUM / YSQ + C(I)
350 CONTINUE
SUM = SUM/Y - Y + SQRTPI
SUM = SUM + (Y-HALF)*LOG(Y)
RES = DEXP(SUM)
ELSE
RES = XINF
GO TO 900
END IF
END IF
IF (PARITY) RES = -RES
IF (FACT .NE. ONE) RES = FACT / RES
900 DGAMMA = RES
RETURN
END
C
C ******************************************************************************
C
|
