c Copyright 1996 (c) by Marne C. Cario and Barry L. Nelson C Marne Cario C Original code written October, 1994, for use with IMSL C Code updated 1/96 to omit IMSL routines c Subsequent updates: 1/21/96, 2/17/96 c c Current Release: 1.1 5/30/96 program generate integer n,p,nobs,maxobs,p1,distrib,elimit,nempir, &count,ndim,lag,nn,ifault,nullty parameter(n=21,maxobs=20000,elimit=20,nn=231) real zinit(n),y(maxobs),yavg,zavg, &ycorr(n),zcorr(n),yrho(n),c,svary,svarz,dparam(5), &sigma,z(maxobs),empcdf(elimit),prob(elimit), &epsilon(n),sum,zvari real*8 zlt(nn),t(nn),zrho(n),alpha(n) character*10 observ,stats integer*4 seed external chol,yw c open(15,file='temp.out',status='unknown') open(8,file='input.gen',status='old') read(8,1)observ read(8,1)stats 1 format(a10) open(9,file=observ,status='unknown') open(10,file=stats,status='unknown') call input(n,maxobs,elimit,p,zrho,distrib,nobs,dparam, &yrho,nempir,empcdf,prob) seed=748932582 c Calculate the AR parameters, alpha if(p.ne.1)then Call YW(p,zrho,alpha) else alpha(1)=zrho(1) endif c Calculate sigma, the square root of the error variance sum=0.0 do i=1,p sum=sum+sngl(alpha(i)*zrho(i)) enddo if(distrib.ne.1)then sigma=sqrt(1.0-sum) else sigma=sqrt((1.0-sum)*dparam(2)) endif c Construct the implied COVARIANCE matrix (lower triangular portion) c This is just the correlation matrix if the ARTA distribution c is not normal. If the ARTA distribution is normal, then this c is the covariance matrix. count=0 ndim=(p+1)*(p+2)/2 if(distrib.ne.1)then zvari=1.0 else zvari=dparam(2) endif do i=1,p+1 do j=1,i count=count+1 if(j.eq.i)then zlt(count)=dble(zvari) else lag=i-j zlt(count)=zrho(lag)*dble(zvari) endif enddo enddo C GET THE CHOLESKY DECOMPOSITION OF THE IMPLIED COVARIANCE C MATRIX p1=p+1 call chol(zlt,p1,ndim,t,nullty,ifault) C GENERATE A (p+1)x1 VECTOR OF I.I.D. NORMALS call geneps(p1,epsilon,seed) C GENERATE a (p+1)x1 VECTOR OF STD. NORMALS WITH THE IMPLIED C COVARIANCE MATRIX and MEAN ZERO do i=1,p1 zinit(i)=0.0 k=i*(i-1)/2 do m=1,i zinit(i)=zinit(i)+sngl(t(k+m))*epsilon(m) enddo enddo C GENERATE THE OBSERVATIONS if(distrib.ne.1)then call geny(zinit,y,nobs,alpha,p,p1,distrib,sigma,seed, &dparam,z,nempir,empcdf,prob) call avg(y,nobs,yavg) call avg(z,nobs,zavg) call sampvarb(y,nobs,yavg,svary) call sampvarb(z,nobs,zavg,svarz) do i=1,p call corr(y,i,nobs,yavg,svary,c) ycorr(i)=c call corr(z,i,nobs,zavg,svarz,c) zcorr(i)=c enddo call outy(y,nobs) call output(yavg,zavg,svary,svarz,ycorr,zcorr,p) else call genz(zinit,nobs,alpha,p,p1,sigma,seed,z) c ADD THE MEAN OF THE AR PROCESS TO EACH Z do i=1,nobs z(i)=z(i)+dparam(1) enddo call avg(z,nobs,zavg) call sampvarb(z,nobs,zavg,svarz) do i=1,p call corr(z,i,nobs,zavg,svarz,c) zcorr(i)=c enddo call outy(z,nobs) write(10,29) 29 format('***** SUMMARY STATISTICS *****') write(10,30)zavg,svarz 30 format(/,'Average of AR Process: ',f15.5, &/,'Sample Variance of AR Process: ',f15.5) write(10,35) 35 format('Lag',5x,'Correlation of AR Process') do i=1,p write(10,40)i,zcorr(i) 40 format(i2,10x,f10.8) enddo endif close(9) close(10) c close(15) stop end subroutine geneps(p1,epsilon,seed) integer p1,ifault integer*4 seed real epsilon(p1) real*8 u,vnorm2 external vnorm2 Do i=1,p1 u=unif(seed) epsilon(i)=sngl(vnorm2(u,ifault)) Enddo Return End subroutine input(n,maxobs,elimit,p,zrho,distrib,nobs, &dparam,yrho,nempir,empcdf,prob) integer p,n,nobs,distrib,i,nu,nempir,elimit real yrho(n),dparam(5),mean,variance, &empcdf(elimit),prob(elimit),probab real*8 zrho(n) read(8,10)distrib,p,nobs 10 format(i2,/,i2,/,i10) c distrib Distribution c ******* ************ c 1 Normal c 2 t c 3 Uniform c 4 Exponential c 5 Gamma c 6 Weibull c 7 Lognormal c 8 Johnson Unbounded c 9 Johnson Bounded c 10 Empirical c 11 Discrete with Finite Support if((distrib.lt.1).or.(distrib.gt.11))then write(10,15) 15 format('Invalid input for distribution number.' &,/,'Distribution must an integer from 1 to 11.' &,/,'1 = Normal ',/,'2 = t',/,'3 = Continuous Uniform' &,/,'4 = Exponential',/,'5 = Gamma',/,'6 = Weibull' &,/,'7 = Lognormal',/,'8 = Johnson Unbounded' &,/,'9 = Johnson Bounded',/,'10 = Empirical' &,/,'Discrete with Finite Support') close(8) close(9) close(10) stop endif i=n-1 if(p.gt.i)then write(10,17)i 17 format('Maximum number of lags is ',i2) close(8) close(9) close(10) stop endif if(nobs.gt.maxobs)then write(10,18)maxobs 18 format('Maximum number of observations is ',i5) close(8) close(9) close(10) stop endif do i=1,p read(8,20)zrho(i) 20 format(f10.8) enddo if(distrib.ne.1)then do i=1,p read(8,30)yrho(i) 30 format(f10.8) enddo endif nempir=2 mean=0.0 variance=0.0 goto(110,120,130,140,150,160,170,180,190,200,210),distrib c NORMAL DISTRIBUTION c dparam(1) = MEAN c dparam(2) = VARIANCE 110 read(8,540)dparam(1),dparam(2) close(8) if(dparam(2).le.0.0)then write(10,115) 115 format('Variance must be > 0.') close(9) close(10) stop endif goto 1000 c STUDENT'S t DISTRIBUTION c dparam(1) = DEGREES OF FREEDOM 120 read(8,550)nu close(8) dparam(1)=dble(nu) if(dparam(1).lt.3)then write(10,125) 125 format('Degrees of freedom must be > 2.') close(9) close(10) stop endif goto 1000 c CONTINUOUS UNIFORM DISTRIBUTION c dparam(1) = LEFT ENDPOINT c dparam(2) = RIGHT ENDPOINT 130 read(8,540)dparam(1),dparam(2) close(8) goto 1000 c EXPONENTIAL DISTRIBUTION c dparam(1) = RATE 140 read(8,560)dparam(1) close(8) if(dparam(1).le.0.0)then write(10,145) 145 format('Rate must be > 0.') close(9) close(10) stop endif goto 1000 c GAMMA DISTRIBUTION c dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 150 read(8,540)dparam(1),dparam(2) close(8) if((dparam(1).le.0.0).or.(dparam(2).le.0.0))then write(10,155) 155 format('Scale and shape parameters must be > 0.') close(9) close(10) stop endif goto 1000 c WEIBULL DISTRIBUTION C dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 160 read(8,540)dparam(1),dparam(2) close(8) if((dparam(1).le.0.0).or.(dparam(2).le.0.0))then write(10,165) 165 format('Scale and shape parameters must be > 0.') close(9) close(10) stop endif goto 1000 c LOGNORMAL DISTRIBUTION c dparam(1)=scale parameter mu c dparam(2)=square of the shape parameter sigma 170 read(8,540)mean,variance close(8) if(variance.le.0.0)then write(10,175) 175 format('Variance must be > 0.') close(9) close(10) stop endif y=2.0*alog(mean) x=alog(variance+mean*mean) dparam(1)=y-x/2.0 dparam(2)=x-y goto 1000 c JOHNSON UNBOUNDED FAMILY c F(y) = Phi{dparam(1)+dparam(2)*f[(y-dparam(3))/dparam(4)]} c Where f(x)=inverse_sinh(x) c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter 180 do i=1,4 read(8,560)dparam(i) enddo close(8) goto 1000 c JOHNSON BOUNDED FAMILY c F(y) = Phi{dparam(1)+dparam(2)*f[(y-dparam(3))/dparam(4)]} c Where f(x)=log[x/(1-x)] c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter 190 do i=1,4 read(8,560)dparam(i) enddo close(8) goto 1000 c EMPIRICAL CDF c dparam(1) = NUMBER OF POINTS IN CDF 200 read(8,550)nempir close(8) dparam(1)=nempir if(nempir.gt.elimit)then write(10,205)elimit 205 format('Maximum number of points in empirical cdf is ',i4) close(9) close(10) stop endif call emprcal(nempir,empcdf) goto 1000 c DISCRETE DISTRIBUTION WITH FINITE SUPPORT c dparam(1) = NUMBER OF POINTS IN CDF 210 read(8,550)nempir dparam(1)=nempir if(nempir.gt.elimit)then write(10,215)elimit 215 format('Maximum number of points in pmf is ',i4) close(9) close(10) stop endif do i=1,nempir read (8,*) empcdf(i),prob(i) enddo close(8) call bubsort3(empcdf,prob,nempir) C CONVERT PMF TO CDF probab=0.0 do i=1,nempir probab=probab+prob(i) prob(i)=probab enddo if(dabs(probab-dble(1.0)).gt.0.0001)then write(10,220) 220 format(/,'Input probabilities do not add up to one.') close(9) close(10) stop endif goto 1000 500 format(i1) 501 format(i5) 502 format(i2) 510 format(i2) 520 format(i2) 530 format(f15.13,/,f7.5) 540 format(f10.5,/,f10.5) 550 format(i4) 560 format(f10.5) 570 format(i4,/,i4) 580 format(i4,/,f7.5) 590 format(f7.5) 600 format(i10,/,f7.5) 610 format(i4,/,i5,/,i5) 1000 call echoin(p,yrho,zrho,dparam,empcdf,prob, &nempir,distrib,mean,variance) return end subroutine echoin(p,yrho,zrho,dparam,empcdf,prob, &nempir,distrib,mean,variance) integer p,nempir,distrib real yrho(p),dparam(5),empcdf(nempir),prob(nempir), &mean,variance real*8 zrho(p) write(10,5) 5 format('********** USER INPUT **********',/) goto(10,20,30,40,50,60,70,80,90,100,110),distrib 10 write(10,15) 15 format('ARTA Distribution: Normal') write(10,16)dparam(1),dparam(2) 16 format('Mean: ',f10.5,/,'Variance: ',f10.5) goto 165 20 write(10,25) 25 format('ARTA Distribution: t') nu=dparam(1) write(10,26)nu 26 format('Degrees of freedom: ',i4) goto 150 30 write(10,35) 35 format('ARTA Distribution: Continuous Uniform') write(10,36)dparam(1),dparam(2) 36 format('Left Endpoint: ',f10.5,/,'Right Endpoint: ',f10.5) goto 150 40 write(10,45) 45 format('ARTA Distribution: Exponential') write(10,46)dparam(1) 46 format('Rate: ',f10.5) goto 150 50 write(10,55) 55 format('ARTA Distribution: Gamma') write(10,56)dparam(1),dparam(2) 56 format('Scale Parameter: ',f10.5,/,'Shape Parameter: ', &f10.5) goto 150 60 write(10,65) 65 format('ARTA Distribution: Weibull') write(10,66)dparam(1),dparam(2) 66 format('Scale Parameter: ',f10.5,/,'Shape Parameter: ', &f10.5) goto 150 70 write(10,75) 75 format('ARTA Distribution: Lognormal') write(10,76)mean,variance 76 format('Mean: ',f10.5,/,'Variance: ', &f10.5) goto 150 80 write(10,85) 85 format('ARTA Distribution: Johnson Unbounded') write(10,86)dparam(1),dparam(2),dparam(3),dparam(4) 86 format('Shape Parameter 1: ',f10.5,/,'Shape Parameter 2: ', &f10.5,/,'Location Parameter: ',f10.5,/,'Scale Parameter: ', &f10.5) goto 150 90 write(10,95) 95 format('ARTA Distribution: Johnson Bounded') write(10,96)dparam(1),dparam(2),dparam(3),dparam(4) 96 format('Shape Parameter 1: ',f10.5,/,'Shape Parameter 2: ', &f10.5,/,'Location Parameter: ',f10.5,/,'Scale Parameter: ', &f10.5) goto 150 100 write(10,105) 105 format('ARTA Distribution: Empirical') nempir=dparam(1) write(10,106)nempir 106 format('Number of points: ',i4) goto 150 110 write(10,115) 115 format('ARTA Distribution: Discrete with finite support') nempir=dparam(1) write(10,116)nempir 116 format('Number of points: ',i4) write(10,117) 117 format(/,3x,'Value',7x,'Cumulative Probability', &/,3x,'*****',7x,'*********************') do i=1,nempir write(10,118)empcdf(i),prob(i) enddo 118 format(f10.5,10x,f7.5) goto 150 150 write(10,155) 155 format(/,'Lag',5x,'ARTA Corr.',5x,'AR Corr.') do i=1,p write(10,160)i,yrho(i),zrho(i) 160 format(i2,7x,f6.4,9x,f6.4) enddo goto 169 165 write(10,166) 166 format(/'Lag',5x,'AR Corr.') do i=1,p write(10,167)i,zrho(i) 167 format(i2,7x,f6.4) enddo 169 write(10,170) 170 format(' **********************************',/) return end subroutine genz(zinit,nobs,alpha,p,p1,sigma,seed,z) integer*4 seed integer p,p1,nobs,ifault,i real zinit(p1),z(nobs),sigma,e real*8 u,vnorm2,alpha(p) z(1)=zinit(p1) do i=2,p z(i)=0 do j=1,p z(i)=z(i)+sngl(alpha(j))*zinit(p+2-j) enddo u=dble(unif(seed)) e=sngl(vnorm2(u,ifault)) z(i)=z(i)+sigma*e call redoz(zinit,p1,z(i)) enddo do i=p+1,nobs z(i)=0 do j=1,p z(i)=z(i)+sngl(alpha(j))*z(i-j) enddo u=dble(unif(seed)) e=sngl(vnorm2(u,ifault)) z(i)=z(i)+sigma*e enddo return end subroutine geny(zinit,y,nobs,alpha,p,p1,distrib,sigma,seed, &dparam,z,nempir,empcdf,prob) integer*4 seed integer p,p1,nobs,distrib,ifault,nempir,i real zinit(p1),y(nobs),z(nobs),sigma,e, &dparam(5),empcdf(nempir),prob(nempir),yy real*8 z1,u1,q,pdf,u,vnorm2,alpha(p),w(4) w(1)=-1.0 z(1)=zinit(p1) z1=dble(z(1)) call fnorm2(z1,u1,q,pdf) call invert(u1,distrib,dparam,empcdf,prob,nempir,yy,w) y(1)=yy do i=2,p z(i)=0 do j=1,p z(i)=z(i)+sngl(alpha(j))*zinit(p+2-j) enddo u=dble(unif(seed)) e=sngl(vnorm2(u,ifault)) z(i)=z(i)+sigma*e z1=dble(z(i)) call fnorm2(z1,u1,q,pdf) call invert(u1,distrib,dparam,empcdf,prob,nempir,yy,w) y(i)=yy call redoz(zinit,p1,z(i)) enddo do i=p+1,nobs z(i)=0 do j=1,p z(i)=z(i)+sngl(alpha(j))*z(i-j) enddo u=dble(unif(seed)) e=sngl(vnorm2(u,ifault)) z(i)=z(i)+sigma*e z1=dble(z(i)) call fnorm2(z1,u1,q,pdf) call invert(u1,distrib,dparam,empcdf,prob,nempir,yy,w) y(i)=yy enddo return end subroutine invert(u1,distrib,dparam,empcdf,prob,nempir,y,w) integer distrib,nempir,ifault real u1s,y,dparam(5),uinv,expinv,weibinv,lninv, &juinv,jbinv,empinv,discinv,g,prob(nempir),empcdf(nempir) real*8 u1,gaminv,alph,w(4),vstud,df goto(10,20,30,40,50,60,70,80,90,100,110),distrib C NORMAL DISTRIBUTION (Special case - should never get here) 10 print*,'Error - specialized routine for Normal distribution.' stop C STUDENT'S t DISTRIBUTION 20 df=dble(dparam(1)) y=sngl(vstud(u1,df,w,ifault)) c write(15,*)u1,y,df return C UNIFORM DISTRIBUTION 30 u1s=sngl(u1) y=uinv(dparam,u1s) return C EXPONENTIAL DISTRIBUTION 40 u1s=sngl(u1) y=expinv(dparam,u1s) return C GAMMA DISTRIBUTION 50 alph=dble(dparam(2)) g=sngl(gaminv(alph,u1,ifault)) y=dparam(1)*g return C WEIBULL DISTRIBUTION 60 u1s=sngl(u1) y=weibinv(dparam,u1s) return C LOGNORMAL 70 u1s=sngl(u1) y=lninv(dparam,u1s) return C JOHNSON UNBOUNDED FAMILY 80 u1s=sngl(u1) y=juinv(dparam,u1s) return C JOHNSON BOUNDED FAMILY 90 u1s=sngl(u1) y=jbinv(dparam,u1s) return C EMPIRICAL DISTRIBUTION 100 u1s=sngl(u1) y=empinv(u1s,empcdf,nempir) return C DISCRETE WITH FINITE SUPPORT 110 u1s=sngl(u1) y=discinv(u1s,empcdf,prob,nempir) return end subroutine avg(y,nobs,yavg) integer nobs real y(nobs),sum,yavg sum=0 do i=1,nobs sum=sum+y(i) enddo yavg=sum/nobs return end subroutine corr(y,i,nobs,yavg,czero,c) integer nobs,i,n real y(nobs),c,yavg,czero c write(9,*)'Got there 4aaa' n=nobs-i c=0 c write(9,*)i,nobs,yavg,czero do j=1,n c write(9,*)'Got there 4aa' c=c+(y(j)-yavg)*(y(j+i)-yavg) c write(9,*)'Got there 4ab' enddo c write(9,*)'Got there 4ac' c=c/(nobs*czero) c write(9,*)'Got there 4ad' return end subroutine sampvarb(y,nobs,yavg,svar) integer nobs real y(nobs),yavg,svar svar=0 do j=1,nobs svar=svar+(y(j)-yavg)*(y(j)-yavg) enddo svar=svar/nobs return end subroutine redoz(zinit,p1,znew) integer p1,p real zinit(p1) p=p1-1 do i=2,p zinit(i)=zinit(i+1) enddo zinit(p1)=znew return end subroutine outy(y,nobs) integer nobs real y(nobs) do i=1,nobs write(9,*)y(i) enddo return end subroutine output(yavg,zavg,svary,svarz,ycorr,zcorr,p) integer p real yavg,zavg,svary,svarz,ycorr(p),zcorr(p) write(10,5) 5 format('***** SUMMARY STATISTICS *****') write(10,10)yavg,svary 10 format(/,'Average of Output (Y) Process: ',f10.5, &/,'Sample Variance of Output (Y) Process: ',f10.5) write(10,15) 15 format('Lag',5x,' Correlation of Output (Y) Process') do i=1,p write(10,20)i,ycorr(i) 20 format(i2,10x,f10.8) enddo write(10,30)zavg,svarz 30 format(//,'Average of Input (Z) Process: ',f15.5, &/,'Sample Variance of Input (Z) Process: ',f15.5) write(10,35) 35 format('Lag',5x,'Correlation of Input (Z) Process') do i=1,p write(10,40)i,zcorr(i) 40 format(i2,10x,f10.8) enddo return end real function expinv(dparam,u1) c c function to generate negative exponential variate c via inverse cdf c real mean,u1,dparam(5) if(u1.gt.0.999999)then u1=dble(0.999999) endif mean=1.0/dparam(1) expinv = (-alog(1.0-u1))*mean return end real function weibinv(dparam,u) c function to generate weibull variate via with c shape parameter alpha and scale parameter beta c via the inverse cdf real dparam(5), alpha, beta, u, w beta=dparam(1) alpha=dparam(2) w = (-alog(1. - u))**(1./alpha) weibinv = w*beta return end c real function lninv(dparam, u) c Computes the inverse of the lognormal distribution. c Marne C. Cario, 11/9/95 integer ifault real dparam(5),u,norm,mu,sigma2 real*8 vnorm2,uu mu=dparam(1) sigma2=dparam(2) uu=dble(u) norm=sngl(vnorm2(uu,ifault)) norm=mu+sqrt(sigma2)*norm lninv=exp(norm) return end real function juinv(dparam,u) c Function to invert a Johnson Unbounded family distribution. integer ifault real dparam(5),u,s,term real*8 vnorm2,uu c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter uu=dble(u) term=(sngl(vnorm2(uu,ifault))-dparam(1))/dparam(2) s=sinh(term) juinv=dparam(3)+dparam(4)*s return end real function jbinv(dparam,u) c Function to invert a Johnson Bounded family distribution. integer ifault real dparam(5),u,term,num,denom real*8 uu,vnorm2 c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter uu=dble(u) term=(sngl(vnorm2(uu,ifault))-dparam(1))/dparam(2) num=exp(term) denom=1.0+num jbinv=dparam(3)+dparam(4)*num/denom return end FUNCTION UNIF(IX) c c portable random number generator implementing the c recursion: ix = 16807*ix mod (2**(31) - 1) c using only 32 bits, including sign c c some compilers require the declaration: c integer*4 ix, k1 c c input: c ix = integer greater than 0 and less than 2147483647 c c output: c ix = new pseudorandom value, c unif = a uniform fraction between 0 and 1 c c***code from Bratley, Fox and Schrage (1983), page 202 c modified by BLN to generate antithetic variate c when ix is negative c c***seed values 131,072 apart: 748932582, 1985072130, c 1631331038, 67377721, 366304404 c INTEGER*4 IX, K1, AIX AIX = ABS(IX) K1 = AIX/127773 AIX = 16807*(AIX - K1*127773) - K1*2836 IF (AIX .LT. 0) AIX = AIX + 2147483647 IF (IX .GE. 0) THEN UNIF = AIX*4.656612875E-10 IX = AIX ELSE UNIF = 1. - AIX*4.656612875E-10 IX = -AIX ENDIF C RETURN END REAL FUNCTION EMPINV(u,empcdf,nempir) C Inverts an empirical cdf using the binary-search method in C Bratley, Fox, and Schrage (Springer-Verlag, 1987), pp. 147-148. c Marne C. Cario c ARTA Subroutines c Updates: 11/20/95 (Tested) c Subroutines called: None c Calling program: inner c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c u double precision value between 0.0 and 1.0 at which we want to c invert the empirical cdf c empcdf double-prec. vector vector containing the points of the empirical c cdf in ascending order c nempir integer number of points in the empirical cdf c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c empinv double precision value of the inverted empirical cdf at u c ***************************************************************************** integer nempir,l,r,k,ix real u,empcdf(nempir),fx l=0 r=nempir k=r-1 dowhile(l.lt.k) ix=(l+r)/2.0 fx=real(ix)/nempir if(u.gt.fx)then l=ix else r=ix endif k=r-1 enddo empinv=empcdf(r) RETURN END SUBROUTINE BUBSORT3(data,attrib,n) c Does bubble sort on array of length n and carries a real c attribute along with the sorted array. Sorts from low to high. integer j,n real data(n),attrib(n),temp1,temp2 if(n.gt.1)then do i=1,n do j=1,n-i if(data(j).gt.data(j+1))then temp1=data(j) temp2=attrib(j) data(j)=data(j+1) attrib(j)=attrib(j+1) data(j+1)=temp1 attrib(j+1)=temp2 endif enddo enddo endif return end subroutine bubsort(data,n) c c subroutine to do bubble sort on array of length n c sorts from low to high c integer i,j,n real data(n), temp c do i = 1, n do j = 1, n - i if (data(j) .gt. data(j+1)) then temp = data(j) data(j) = data(j+1) data(j+1) = temp endif enddo enddo return end SUBROUTINE EMPRCAL(nempir,empcdf) C Reads in nempir points of an empirical c.d.f. from the file empir.dat C and sorts them into ascending order c Marne C. Cario c ARTA Subroutines c Updates: 8/95, 10/10/95, 11/20/95 (Tested) c Subroutines called: bubsort c Calling program: input c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c nempir integer number of points in the empirical cdf c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c empcdf double-prec. vector points of the empirical cdf in ascending order c ***************************************************************************** integer nempir real empcdf(nempir) external bubsort open(11,file='empir.dat',status='old') do i=1,nempir read(11,*)empcdf(i) enddo call bubsort(empcdf,nempir) c do i=1,10 c write(8,*)empcdf(i) c enddo close(11) RETURN END real function uinv(dparam,u) real dparam(5),u uinv=dparam(1)+u*(dparam(2)-dparam(1)) return end real*8 FUNCTION vnorm2(P, IFAULT) c Name changed by MCC 1/4/96; originally named PPND16 C C ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 C C Produces the normal deviate Z corresponding to a given lower C tail area of P; Z is accurate to about 1 part in 10**16. C C The hash sums below are the sums of the mantissas of the C coefficients. They are included for use in checking C transcription. C real*8 ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, & CONST2, A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, & B4, B5, B6, B7, & C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, & D6, D7, E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, & F4, F5, F6, F7, P, Q, R PARAMETER (ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0, & SPLIT1 = 0.425D0, SPLIT2 = 5.D0, & CONST1 = 0.180625D0, CONST2 = 1.6D0) C C Coefficients for P close to 0.5 C PARAMETER (A0 = 3.38713 28727 96366 6080D0, * A1 = 1.33141 66789 17843 7745D+2, * A2 = 1.97159 09503 06551 4427D+3, * A3 = 1.37316 93765 50946 1125D+4, * A4 = 4.59219 53931 54987 1457D+4, * A5 = 6.72657 70927 00870 0853D+4, * A6 = 3.34305 75583 58812 8105D+4, * A7 = 2.50908 09287 30122 6727D+3, * B1 = 4.23133 30701 60091 1252D+1, * B2 = 6.87187 00749 20579 0830D+2, * B3 = 5.39419 60214 24751 1077D+3, * B4 = 2.12137 94301 58659 5867D+4, * B5 = 3.93078 95800 09271 0610D+4, * B6 = 2.87290 85735 72194 2674D+4, * B7 = 5.22649 52788 52854 5610D+3) C HASH SUM AB 55.88319 28806 14901 4439 C C Coefficients for P not close to 0, 0.5 or 1. C PARAMETER (C0 = 1.42343 71107 49683 57734D0, * C1 = 4.63033 78461 56545 29590D0, * C2 = 5.76949 72214 60691 40550D0, * C3 = 3.64784 83247 63204 60504D0, * C4 = 1.27045 82524 52368 38258D0, * C5 = 2.41780 72517 74506 11770D-1, * C6 = 2.27238 44989 26918 45833D-2, * C7 = 7.74545 01427 83414 07640D-4, * D1 = 2.05319 16266 37758 82187D0, * D2 = 1.67638 48301 83803 84940D0, * D3 = 6.89767 33498 51000 04550D-1, * D4 = 1.48103 97642 74800 74590D-1, * D5 = 1.51986 66563 61645 71966D-2, * D6 = 5.47593 80849 95344 94600D-4, * D7 = 1.05075 00716 44416 84324D-9) C HASH SUM CD 49.33206 50330 16102 89036 C C Coefficients for P near 0 or 1. C PARAMETER (E0 = 6.65790 46435 01103 77720D0, * E1 = 5.46378 49111 64114 36990D0, * E2 = 1.78482 65399 17291 33580D0, * E3 = 2.96560 57182 85048 91230D-1, * E4 = 2.65321 89526 57612 30930D-2, * E5 = 1.24266 09473 88078 43860D-3, * E6 = 2.71155 55687 43487 57815D-5, * E7 = 2.01033 43992 92288 13265D-7, * F1 = 5.99832 20655 58879 37690D-1, * F2 = 1.36929 88092 27358 05310D-1, * F3 = 1.48753 61290 85061 48525D-2, * F4 = 7.86869 13114 56132 59100D-4, * F5 = 1.84631 83175 10054 68180D-5, * F6 = 1.42151 17583 16445 88870D-7, * F7 = 2.04426 31033 89939 78564D-15) C HASH SUM EF 47.52583 31754 92896 71629 C IFAULT = 0 Q = P - HALF IF (ABS(Q) .LE. SPLIT1) THEN R = CONST1 - Q * Q vnorm2 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) * * R + A2) * R + A1) * R + A0) / * (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) * * R + B2) * R + B1) * R + ONE) RETURN ELSE IF (Q .LT. ZERO) THEN R = P ELSE R = ONE - P END IF IF (R .LE. ZERO) THEN IFAULT = 1 vnorm2 = ZERO RETURN END IF R = SQRT(-LOG(R)) IF (R .LE. SPLIT2) THEN R = R - CONST2 vnorm2 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3) * * R + C2) * R + C1) * R + C0) / * (((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) * * R + D2) * R + D1) * R + ONE) ELSE R = R - SPLIT2 vnorm2 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3) * * R + E2) * R + E1) * R + E0) / * (((((((F7 * R + F6) * R + F5) * R + F4) * R + F3) * * R + F2) * R + F1) * R + ONE) END IF IF (Q .LT. ZERO) vnorm2 = - vnorm2 RETURN END IF END SUBROUTINE fnorm2(Z, P, Q, PDF) c Subroutine renamed by MCC on 1/4/96; original name was NORMP C C Normal distribution probabilities accurate to 1.e-15. C Z = no. of standard deviations from the mean. C P, Q = probabilities to the left & right of Z. P + Q = 1. C PDF = the probability density. C C Based upon algorithm 5666 for the error function, from: C Hart, J.F. et al, 'Computer Approximations', Wiley 1968 C C Programmer: Alan Miller C C Latest revision - 30 March 1986 C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DATA P0, P1, P2, P3, P4, P5, P6/220.20 68679 12376 1D0, * 221.21 35961 69931 1D0, 112.07 92914 97870 9D0, * 33.912 86607 83830 0D0, 6.3739 62203 53165 0D0, * .70038 30644 43688 1D0, .35262 49659 98910 9D-01/, * Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7/440.41 37358 24752 2D0, * 793.82 65125 19948 4D0, 637.33 36333 78831 1D0, * 296.56 42487 79673 7D0, 86.780 73220 29460 8D0, * 16.064 17757 92069 5D0, 1.7556 67163 18264 2D0, * .88388 34764 83184 4D-1/, * CUTOFF/7.071D0/, ROOT2PI/2.5066 28274 63100 1D0/ C ZABS = ABS(Z) C C |Z| > 37. C IF (ZABS .GT. 37.D0) THEN PDF = 0.D0 IF (Z .GT. 0.D0) THEN P = 1.D0 Q = 0.D0 ELSE P = 0.D0 Q = 1.D0 END IF RETURN END IF C C |Z| <= 37. C EXPNTL = EXP(-0.5D0*ZABS**2) PDF = EXPNTL/ROOT2PI C C |Z| < CUTOFF = 10/sqrt(2). C IF (ZABS .LT. CUTOFF) THEN P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS + * P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS + * Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS + * Q0) C C |Z| >= CUTOFF. C ELSE P = PDF/(ZABS + 1.D0/(ZABS + 2.D0/(ZABS + 3.D0/(ZABS + 4.D0/ * (ZABS + 0.65D0))))) END IF C IF (Z .LT. 0.D0) THEN Q = 1.D0 - P ELSE Q = P P = 1.D0 - Q END IF RETURN END c This file contains two algorithms for the logarithm of the gamma function. c Algorithm AS 245 is the faster (but longer) and gives an accuracy of about c 10-12 significant decimal digits except for small regions around X = 1 and c X = 2, where the function goes to zero. c The second algorithm is not part of the AS algorithms. It is slower but c gives 14 or more significant decimal digits accuracy, except around X = 1 c and X = 2. The Lanczos series from which this algorithm is derived is c interesting in that it is a convergent series approximation for the gamma c function, whereas the familiar series due to De Moivre (and usually wrongly c called Stirling's approximation) is only an asymptotic approximation, as c is the true and preferable approximation due to Stirling. c c c DOUBLE PRECISION FUNCTION ALNGAM(XVALUE, IFAULT) C C ALGORITHM AS245 APPL. STATIST. (1989) VOL. 38, NO. 2 C C Calculation of the logarithm of the gamma function C INTEGER IFAULT DOUBLE PRECISION ALR2PI, FOUR, HALF, ONE, ONEP5, R1(9), R2(9), + R3(9), R4(5), TWELVE, X, X1, X2, XLGE, XLGST, XVALUE, + Y, ZERO C C Coefficients of rational functions C DATA R1/-2.66685 51149 5D0, -2.44387 53423 7D1, + -2.19698 95892 8D1, 1.11667 54126 2D1, + 3.13060 54762 3D0, 6.07771 38777 1D-1, + 1.19400 90572 1D1, 3.14690 11574 9D1, + 1.52346 87407 0D1/ DATA R2/-7.83359 29944 9D1, -1.42046 29668 8D2, + 1.37519 41641 6D2, 7.86994 92415 4D1, + 4.16438 92222 8D0, 4.70668 76606 0D1, + 3.13399 21589 4D2, 2.63505 07472 1D2, + 4.33400 02251 4D1/ DATA R3/-2.12159 57232 3D5, 2.30661 51061 6D5, + 2.74647 64470 5D4, -4.02621 11997 5D4, + -2.29660 72978 0D3, -1.16328 49500 4D5, + -1.46025 93751 1D5, -2.42357 40962 9D4, + -5.70691 00932 4D2/ DATA R4/ 2.79195 31791 8525D-1, 4.91731 76105 05968D-1, + 6.92910 59929 1889D-2, 3.35034 38150 22304D0, + 6.01245 92597 64103D0/ C C Fixed constants C DATA ALR2PI/9.18938 53320 4673D-1/, FOUR/4.D0/, HALF/0.5D0/, + ONE/1.D0/, ONEP5/1.5D0/, TWELVE/12.D0/, ZERO/0.D0/ C C Machine-dependant constants. C A table of values is given at the top of page 399 of the paper. C These values are for the IEEE double-precision format for which C B = 2, t = 53 and U = 1023 in the notation of the paper. C DATA XLGE/5.10D6/, XLGST/1.D+30/ C X = XVALUE ALNGAM = ZERO C C Test for valid function argument C IFAULT = 2 IF (X .GE. XLGST) RETURN IFAULT = 1 IF (X .LE. ZERO) RETURN IFAULT = 0 C C Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined C IF (X .LT. ONEP5) THEN IF (X .LT. HALF) THEN ALNGAM = -LOG(X) Y = X + ONE C C Test whether X < machine epsilon C IF (Y .EQ. ONE) RETURN ELSE ALNGAM = ZERO Y = X X = (X - HALF) - HALF END IF ALNGAM = ALNGAM + X * ((((R1(5)*Y + R1(4))*Y + R1(3))*Y + + R1(2))*Y + R1(1)) / ((((Y + R1(9))*Y + R1(8))*Y + + R1(7))*Y + R1(6)) RETURN END IF C C Calculation for 1.5 <= X < 4.0 C IF (X .LT. FOUR) THEN Y = (X - ONE) - ONE ALNGAM = Y * ((((R2(5)*X + R2(4))*X + R2(3))*X + R2(2))*X + + R2(1)) / ((((X + R2(9))*X + R2(8))*X + R2(7))*X + + R2(6)) RETURN END IF C C Calculation for 4.0 <= X < 12.0 C IF (X .LT. TWELVE) THEN ALNGAM = ((((R3(5)*X + R3(4))*X + R3(3))*X + R3(2))*X + R3(1)) / + ((((X + R3(9))*X + R3(8))*X + R3(7))*X + R3(6)) RETURN END IF C C Calculation for X >= 12.0 C Y = LOG(X) ALNGAM = X * (Y - ONE) - HALF * Y + ALR2PI IF (X .GT. XLGE) RETURN X1 = ONE / X X2 = X1 * X1 ALNGAM = ALNGAM + X1 * ((R4(3)*X2 + R4(2))*X2 + R4(1)) / + ((X2 + R4(5))*X2 + R4(4)) RETURN END c c c c REAL FUNCTION DISCINV(u,points,prob,nempir) C Inverts a discrete distribution with finite support given C the pmf points and their cdf values using the binary-search C method in Bratley, Fox, and Schrage (Springer-Verlag, 1987), C pp. 147-148. c Marne C. Cario c ARTA Subroutines c Updates: 1/21/96 integer nempir,l,r,k,ix real u,points(nempir),prob(nempir) l=0 r=nempir k=r-1 dowhile(l.lt.k) ix=(l+r)/2.0 if(u.gt.prob(ix))then l=ix else r=ix endif k=r-1 enddo discinv=points(r) RETURN END real*8 function gaminv(alpha,phi,ifault) c modified by MCC 1/4/96 to work with ARTA code c***corrected typing error found by Tao 1/17/89 c c evaluate the percentage points of the gamma(alpha, 1) c probability distribution function. c c adapted from function vchi2 in Bratley, Fox and Schrage (1983) c c inputs: c phi = probability to be inverted c phi should lie in (.000002, .999998) c alpha = shape parameter of gamma distribution > 0 c c outputs: c gaminv = inverse of gamma(alpha) cdf of phi c ifault = 4 if input error, else 0 c c auxiliary routines: c gamain - BLN code c algama - IMSL algorithm for evaluating log of the gamma function c (MCC version uses alngam.for - code obtained from statlib) c mdnris - IMSL algorithm for evaluating the inverse cdf of c the standard normal distribution c (MCC version uses vnorm2.for - code obtained from statlib) c c c desired accuracy and ln(2.): c integer ifault real*8 phi,alpha,x,v,g,xx,c,ch,q,p1,p2, &t,b,a,gamain,s1,s2,s3,s4,s5,s6,vnorm2,e,alngam,aa data e, aa/0.5e-5, 0.6931471805e0/ c MODIFICATION BY MCC if(phi.ge.0.999998)then phi=dble(0.9999979) endif C END OF MODIFICATION ifault = 0 c c***convert to appropriate chi-square parameters for vchi2 v = 2.*alpha c g = alngam(alpha,ifault) c c***check for unusual input if (v .lt. 0.) goto 9000 if (phi .le. .000002 .or. phi .ge. .999998) goto 9000 xx = alpha c = xx - 1.0 c c***starting approximation for small chi-squared if (v .ge. -1.24*dlog(phi)) goto 1 ch = (phi*xx*dexp(g + xx*aa))**(1./xx) if (ch - e) 6, 4, 4 c c***starting approximation for v less than or equal to .32 1 if (v .gt. 0.32) goto 3 ch = 0.4 a = dlog(1.0 - phi) 2 q = ch p1 = 1.0 + ch*(4.67 + ch) p2 = ch*(6.73 + ch*(6.66 + ch)) t = -0.5 + (4.67 +2.*ch)/p1 - $ (6.73 + ch*(13.32 + 3.0*ch))/p2 ch = ch - (1.0 - dexp(a + g + .5*ch + c*aa)*p2/p1)/t if (dabs(q/ch - 1.) - 0.01) 4, 4, 2 c c***get the corresponding normal deviate c 3 call mdnris(phi,x,ier) c MCC change to use different function 3 x=vnorm2(phi,ifault) c c***starting approximation using Wilson and Hilferty estimate p1 = 0.222222/v ch = v*(x*dsqrt(p1) + 1.0 - p1)**3. c c***starting approximation for p tending to 1 if (ch .gt. 2.2*v + 6.0) $ ch = -2.0*(dlog(1.0 - phi) - c*dlog(.5*ch) +g) c c***call to algorithm as32 and calculation of seven term c Taylor series 4 q = ch p1 = .5*ch p2 = phi - gamain(p1,xx,g,ifail) if (ifail .eq. 0) goto 5 ifault = 4 return 5 t = p2*dexp(xx*aa + g + p1 - c*dlog(ch)) b = t/ch a = .5*t - b*c s1 = (210.0+a*(140.0+a*(105.0+a*(84.0+a*(70.0+60.0*a)))))/420.0 s2 = (420.0+a*(735.0+a*(966.0+a*(1141.0+1278.0*a))))/2520.0 s3 = (210.0 + a*(462.0 + a*(707.0 + 932.0*a)))/2520.0 s4 = (252.0+a*(672.0+1182.0*a)+c*(294.0+a*(889.0+1740.0*a)))/5040. s5 = (84.0 + 264.0*a + c*(175.0 + 606.0*a))/2520.0 s6 = (120.0 + c*(346.0 + 127.0*c))/5040.0 ch = ch+t*(1.0+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6)))))) if (dabs(q/ch - 1.) .gt. e) goto 4 c c***return answer; divide by 2 to convert from chi-squared to gamma 6 gaminv = ch/2. return 9000 ifault = 4 c c***despite fault, try to give a sensible response gaminv = 0. if (phi .gt. .999998) gaminv = (v + 4.*dsqrt(2.*v))/2. return end c c real*8 function gamain(x, p, g, ifail) c c algorithm as32 J.R. Statist. Soc. C. (1970), 19, 3. c c computes incomplete gamma ratio for positive values of c arguments x and p. g must be supplied and is equal to c ln(gamma(p)) c c ifail = 1 if p <= 0, else 2 if x < 0, else 0 c uses series expansion if p > x or x <= 1, otherwise a c continued fraction approximation c integer ifail real*8 pn(6),x,p,g,gin,factor,term,rn,a,b,an,dif, &acu,oflo c c***define accuracy and initialize data acu/1.e-8/, oflo/1.e30/ gin = 0.0 ifail = 0 c c***test for admissibility of arguments if (p .le. 0.) ifail = 1 if (x .lt. 0.) ifail = 2 if (ifail .gt. 0 .or. x .eq. 0.) goto 50 factor = dexp(p*dlog(x) - x - g) if (x .gt. 1.0 .and. x .ge. p) goto 30 c c***calculation by series expansion gin = 1. term = 1. rn = p 20 rn = rn + 1. term = term*x/rn gin = gin + term if (term .gt. acu) goto 20 gin = gin*factor/p goto 50 c c***calculation by continued fraction 30 a = 1.0 - p b = a + x + 1.0 term = 0. pn(1) = 1.0 pn(2) = x pn(3) = x + 1.0 pn(4) = x*b gin = pn(3)/pn(4) 32 a = a + 1.0 b = b + 2.0 term = term + 1.0 an = a*term do 33 i = 1, 2 33 pn(i+4) = b*pn(i+2) - an*pn(i) if (pn(6) .eq. 0.0) goto 35 rn = pn(5)/pn(6) dif = dabs(gin - rn) if (dif .gt. acu) goto 34 if (dif .le. acu*rn) goto 42 34 gin = rn 35 do 36 i = 1, 4 36 pn(i) = pn(i+2) if (dabs(pn(5)) .lt. oflo) goto 32 do 41 i = 1, 4 41 pn(i) = pn(i)/oflo goto 32 42 gin = 1. - factor*gin c 50 gamain = gin return end Subroutine yw(p,r,alpha) c Subroutine to solve the Yule-Walker equations c Marne C. Cario c ARTA Subroutines c Updates: 11/15/95, 11/16/95, 11/17/95 (Tested) c Subroutines called: GaussJord c Calling program: main c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c r double-prec. vector px1 vector of AR-process autocorrelations c p integer number of autocorrelation lags c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c alpha double-prec. vector px1 vector of AR parameters c ***************************************************************************** integer p real*8 r(p),alpha(p),h(6,6),temp1(5) external gaussjord do i=1,p temp1(i)=r(i) enddo c CONSTRUCT THE PSI MATRIX do i=2,p do j=1,i-1 k=i-j h(i,j)=r(k) enddo enddo do i=1,p h(i,i)=1.0 enddo do i=1,p-1 do j=i+1,p h(i,j)=h(j,i) enddo enddo call GaussJord(H,temp1,p,alpha) RETURN END subroutine GaussJord(A,b,dim,soln) c Performs Gauss-Jordan elimination to arrive at the solution c of the system of linear equations Ax=b, where A is a dim x dim matrix, and c b is a dim x 1 vector. The solution is returned in the dim x 1 vector c soln. c Marne C. Cario c ARTA Subroutines c Updates: 11/15/95, 11/17/95 (Tested) c Subroutines called: None c Calling program: main integer dim,i,j,k real*8 A(6,6),b(6),soln(dim),factor c ROW-REDUCE THE iTH ROW do i=1,dim C ROW-REDUCE ALL OTHER ROWS j NOT EQUAL TO i AS WELL AS THE RIGHT-HAND SIDE C ON ROW i (PROVIDED A(j,i) IS NOT ALREADY ZERO) do j=1,dim if((j.ne.i).and.(A(j,i).ne.0.0))then factor=A(i,i)/A(j,i) do k=1,dim A(j,k)=factor*A(j,k)-A(i,k) enddo b(j)=factor*b(j)-b(i) endif enddo enddo do i=1,dim soln(i)=b(i)/A(i,i) enddo return end C This file contains AS6 and the enhanced version ASR44. See AS7 also. C C SUBROUTINE CHOL (A,N,NN,U,NULLTY,IFAULT) C C Algorithm AS6, Applied Statistics, vol.17, (1968) C C Given a symmetric matrix order n as lower triangle in a( ) C calculates an lower triangle, u( ), such that u*u' = a. C a must be positive semi-definite. eta is set to multiplying C factor determining effective zero for pivot. C C arguments:- C a() = input, a +ve definite matrix stored in lower-triangula C form. C n = input, the order of a C nn = input, the size of the a and u arrays n*(n+1)/2 C u() = output, a lower triangular matrix such that u*u' = a. C a & u may occupy the same locations. C nullty = output, the rank deficiency of a. C ifault = output, error indicator C = 1 if n < 1 C = 2 if a is not +ve semi-definite C = 3 if nn < n*(n+1)/2 C = 0 otherwise C C*********************************************************************** C integer n,nn,ifault,j,k,ii,icol,l,kk,irow,m,nullty DOUBLE PRECISION A(NN),U(NN),ETA,ETA2,X,W,ZERO C C The value of eta will depend on the word-length of the C computer being used. See introductory text. C DATA ETA,ZERO/1.D-9,0.0D0/ C IFAULT=1 IF (N.LE.0) RETURN IFAULT=3 IF (NN.LT.N*(N+1)/2) RETURN IFAULT=2 NULLTY=0 J=1 K=0 ETA2=ETA*ETA II=0 C C Factorize column by column, icol = column no. C DO 80 ICOL=1,N II=II+ICOL X=ETA2*A(II) L=0 KK=0 C C IROW = row number within column ICOL C DO 40 IROW=1,ICOL KK=KK+IROW K=K+1 W=A(K) M=J DO 10 I=1,IROW L=L+1 IF (I.EQ.IROW) GO TO 20 W=W-U(L)*U(M) M=M+1 10 CONTINUE 20 IF (IROW.EQ.ICOL) GO TO 50 IF (U(L).EQ.ZERO) GO TO 30 U(K)=W/U(L) GO TO 40 30 IF (W*W.GT.DABS(X*A(KK))) RETURN U(K)=ZERO 40 CONTINUE 50 IF (ABS(W).LE.DABS(ETA*A(K))) GO TO 60 IF (W.LT.ZERO) RETURN U(K)=DSQRT(W) GO TO 70 60 U(K)=ZERO NULLTY=NULLTY+1 70 J=J+ICOL 80 CONTINUE IFAULT=0 END c This file contains AS7 and an enhanced alternative - ASR44. See also AS6. c c subroutine syminv(a, n, nn, c, w, nullty, ifault) c c Algorithm AS7, Applied Statistics, vol.17, 1968, p.198. c c Forms in c( ) as lower triangle, a generalised inverse c of the positive semi-definite symmetric matrix a( ) c order n, stored as lower triangle. c c arguments:- c a() = input, the symmetric matrix to be inverted, stored in c lower triangular form c n = input, order of the matrix c nn = input, the size of the a and c arrays n*(n+1)/2 c c() = output, the inverse of a (a generalized inverse if c is c singular), also stored in lower triangular. c c and a may occupy the same locations. c w() = workspace, dimension at least n. c nullty = output, the rank deficiency of a. c ifault = output, error indicator c = 1 if n < 1 c = 2 if a is not +ve semi-definite c = 3 if nn < n*(n+1)/2 c = 0 otherwise c c*************************************************************************** c double precision a(nn), c(nn), w(n), x, zero, one c data zero, one /0.0d0, 1.0d0/ c c cholesky factorization of a, result in c c call chol(a, n, nn, c, nullty, ifault) if(ifault.ne.0) return c c invert c & form the product (cinv)'*cinv, where cinv is the inverse c of c, row by row starting with the last row. c irow = the row number, ndiag = location of last element in the row. c irow=n ndiag=nn 10 l=ndiag if (c(ndiag) .eq. zero) goto 60 do 20 i=irow,n w(i)=c(l) l=l+i 20 continue icol=n jcol=nn mdiag=nn 30 l=jcol x=zero if(icol.eq.irow) x=one/w(irow) k=n 40 if(k.eq.irow) go to 50 x=x-w(k)*c(l) k=k-1 l=l-1 if(l.gt.mdiag) l=l-k+1 go to 40 50 c(l)=x/w(irow) if(icol.eq.irow) go to 80 mdiag=mdiag-icol icol=icol-1 jcol=jcol-1 go to 30 60 do 70 j=irow,n c(l)=zero l=l+j 70 continue 80 ndiag=ndiag-irow irow=irow-1 if(irow.ne.0) go to 10 return end real*8 FUNCTION vstud(phi, df, w, ifault) C GENERATE THE INVERSE CDF OF THE STUDENT T C Taken from Bratley, Fox, and Schrage (1987, Springer-Verlag), C pp. 340-341. c Typed in by Marne C. Cario, 11/9/95 C INPUTS: C phi = Probability, 0 <= phi <= 1 c df = degrees of freedom c w = work vector of length 4 which depends only on df. c tinv will recalculate, e.g., on first call, if c w(1) < 0 on input. c OUTPUTS: c vstud = F inverse of u, i.e., value such that Pr{X <= tinv} = u. c w = work vector of length 4. w retains information between c calls when df does not change. c ifault = 7 if input error, else 0. c Accuracy in decimal digits of about: c 5 if df >= 8 c machine precision if df = 1 or 2. c 3 otherwise c Reference: G.W. Hill, Communications of the ACM, Vol. 13, No. 10, c Oct. 1970. integer ifault real*8 w(4),phi,phin,phix2,p2tail,y,t,x,y2,c,vnorm2,df, &dwarf,xpiovr2,arg,sinarg,vtemp c Machine epsilon: c BFS VALUE c Data dwarf/ .1e-30/ c MCC VALUE 1/4/96 Data dwarf/ .1e-15/ c Pi over 2 Data xpiovr2/1.5707963268e0/ c write(15,*)phi,df,w(1),w(2),w(3),w(4),ifault c Date 30 Sept 1986 c Check for input error ifault=0 if((phi.lt.0.0).or.(phi.gt.1.0))goto 9000 if(df.le.0)goto 9000 c If it is a Cauchy use a special method if(df.eq.1)goto 300 c There is also a special, exact method if df=2 if(df.eq.2)goto 400 c Check to see if we're not too far out in the tails c Modified by MCC 1/4/96 - t=10.e30 is giving abnormal results c in the ARTA code c ORIGINAL CODE c if((phi.gt.dwarf).and.((1.0-phi).gt.dwarf))goto 205 c t=10.e30 c goto 290 c MODIFIED CODE if((phi.gt.dwarf).and.((1.0-phi).gt.dwarf))then goto 205 elseif(phi.le.dwarf)then phi=0.11e-15 goto 205 else phi=1.0-0.11e-15 goto 205 endif c END OF MODIFICATIONS c General Case c First, convert the cumulative probability to the two-tailed c equivalent used by this method. 205 phix2=2.0*phi p2tail=dmin1(phix2, 2.0-phix2) c Begin the approximation c If w(1) > 0 then skip the computations that depend solely on df. c The values for these variables were stored in w during a previous c call. if(w(1).gt.0.0)goto 250 w(2)=1.0/(df-0.5) w(1)=48.0/(w(2)*w(2)) w(3)=((20700.0*w(2)/w(1)-98.)*w(2)-16.)*w(2) + 96.36 w(4)=((94.5/(w(1)+w(3))-3.)/w(1) + 1.)*dsqrt(w(2)*xpiovr2)*df 250 x = w(4)*p2tail c=w(3) y=x**(2./df) if(y.le.(.05 + w(2)))goto 270 c Asymptotic inverse expansion about normal phin=p2tail*0.5 x=vnorm2(phin, ifault) y=x*x if(df.lt.5) c=c+0.3*(df-4.5)*(x+.6) c=(((.05*w(4)*x-5.)*x-7.)*x-2.)*x+w(1)+c y=(((((.4*y+6.3)*y+36.)*y+94.5)/c-y-3.)/w(1)+1.)*x y=w(2)*y*y if(y.le.0.002)goto 260 y=dexp(y)-1. goto 280 260 y=.5*y*y + y goto 280 270 y2=(df+6.)/(df+y) - 0.089*w(4) - .822 y2=y2*(df+2.)*3. y2=(1./y2+.5/(df+4.))*y -1. y2=y2*(df+1.)/(df+2.) y=y2+1./y 280 t=dsqrt(df*y) 290 if(phi.lt.0.5) t=-t vstud=t return c Inverse of standard Cauchy distribution cdf. c Relative accuracy = machine precision, e.g., 5 digits 300 arg=3.1415926*(1.-phi) sinarg=dmax1(dwarf,dsin(arg)) vstud=dcos(arg)/sinarg return c Special case of df=2. c Relative accuracy = Machine precision, e.g., 5 digits 400 if(phi.gt..5)goto 440 t=2.*phi goto 450 440 t=2.*(1.-phi) 450 if(t.le.dwarf) t=dwarf vtemp=dsqrt((2./(t*(2.-t)))-2.) if(phi.le..5)vtemp=-vtemp vstud=vtemp return 9000 ifault=7 return end