c Copyright 1996 (c) by Marne C. Cario and Barry L. Nelson c Marne C. Cario c ARTA Search Program c Version: 8/30/95, 11/7/95, 11/16/95, 11/20/95, 12/3/95, 1/3/96, c 1/7/96, 2/9/96, 2/14/96 c Current Release: 1.1 6/4/96 PROGRAM MAIN c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c empcdf double-prec. vector vector of points in the empirical distribution c (if applicable) c nempir integer number of points in the empirical distribution c (if applicable) c p integer number of autocorrelation lags to match c rho double-prec. vector px1 vector of autocorrelations to match c errrel double-prec. vector px1 vector of desired relative errors for the c elements of rho c distrib integer code corresponding to the marginal distribution c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c errcode integer error code c 0 = no errors c 1 = distribution not available c 2 = maximum number of autocorrelation lags c exceeded c 3 = distribution parameter infeasible c 4 = correlation out of bounds c 5 = correlation matrix not positive definite c 6 = AR process non-stationary c 7 = roundoff error detected c 8 = maximum number of subintervals exceeded c 9 = integral not converging c 10 = subinterval too small c c r double-prec. vector px1 vector of the AR-process autocorrelations c that result in rho to within errrel c estrho double-prec. vector px1 vector of the calculated values of the ARTA- c process autocorrelations for the values given c in r c esterr double-prec. vector px1 vector of the estimated relative errors of c estrho c alpha double-prec. vector px1 vector of the AR parameters c ***************************************************************************** c *********** VARIABLES ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c laglimit integer maximum number of lags that the c program can match c elimit integer maximum number of points for an empirical c distribution c isym integer 0 if the ARTA distribution is not symmetric c 1 if the ARTA distribution is symmetric c minrho double precision smallest ARTA autocorrelation c maxrho double precision largest ARTA autocorrelation c minfeas double precision minimum feasible bivariate correlation c for distrib c dparam double-prec. vector vector of ARTA-distribution parameters c mean double precision mean of the ARTA distribution c variance double precision variance of the ARTA distribution c rhoincr double-prec. vector px1 vector of the ARTA-process a.c.'s to match c in increasing order c ordrho integer vector px1 vector of the lags corresponding to c the members of rhoincr c remrho integer vector vector of the positions in rhoincr c corresponding to ARTA-process a.c.'s that c remain to be "matched" c iterno integer number of iterations so far c lagsrem integer number of ARTA-process a.c.'s that have not c been achieved within the specified relative c error c dimgrid integer number of r-values for which to carry out c the integration (at the current iteration) c rgrid double-prec. vector dimgridx1 vector of the r-values for which c to carry out the integration (at the current c iteration) c result double-prec. vector vector of the rho values corresponding to c each member in rgrid c a double-prec. vector vector of left brackets for the points in rgrid c b double-prec. vector vector of right brackets for the pts in rgrid c rhoa double-prec. vector vector of ARTA-process a.c.'s corresponding c to the r-values in a c rhob double-prec. vector vector of ARTA-process a.c.'s corresponding c to the r-values in b c gspace1 integer 1.0/gspace1 is the initial spacing between c grid points c gspace2 integer 1.0/gspace2 is the spacing between grid points c after the first iteration c case integer 1 if all ARTA-process a.c.'s are positive c 2 if all ARTA-process a.c.'s are negative c 3 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is not c symmetric c 4 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is c symmetric c maxgrid integer maximum number of grid points c tolout double precision absolute error tolerance for the outer integral c tolin double precision absolute error tolerance for the inner integral c idone integer 0 if at least one ARTA-process a.c. has not c been matched c 1 if all ARTA-process a.c.'s have been matched c istat integer 0 if AR process is stationary c 1 if AR process is not stationary c ***************************************************************************** implicit real*8 (x) external input,echoin,chkpd,subchkpd,chkbnd,emprcal,correl, &order,output,bubsort,initgrid,outer,inner,gkpoints,putheap, &remheap,vstud,empinv,bvnpdf,uniform,normal, &yw,station,procerr,calcerr,search,update,gaussjord integer laglimit,elimit,p,distrib,errcode,isym,nempir,dimgrid, &idone,case,gspace1,gspace2,lagsrem,ordrho(5), &remrho(5),iterno,istat,igrid(5) real*8 rho(5),errrel(5),r(5),estrho(5),esterr(5),alpha(5), &empcdf(300),mean,variance,rgrid(50),maxrho,minrho,minfeas, &tolout,tolin,result(50),rhoincr(5),a(50),b(50),rhoa(50), &rhob(50),dparam(5),z(300),prob(300) character*10 fname parameter(laglimit=5, &elimit=300, &tolout=0.0001, &tolin=0.00005, &gspace1=5, &gspace2=5, &maxsubs=50) common/gkwts/xgwts(15),xkwts(15) common/gkpnts/xgk(15) C GAUSS AND KRONROD WEIGHT VALUES c xgwts(1)=0.0 c xgwts(2)=0.129484966168869693270611432679082 c xgwts(3)=0.0 c xgwts(4)=0.279705391489276667901467771423780 c xgwts(5)=0.0 c xgwts(6)=0.381830050505118944950369775488975 c xgwts(7)=0.0 c xgwts(8)=0.417959183673469387755102040816327 c xgwts(9)=0.0 c xgwts(10)=0.381830050505118944950369775488975 c xgwts(11)=0.0 c xgwts(12)=0.279705391489276667901467771423780 c xgwts(13)=0.0 c xgwts(14)=0.129484966168869693270611432679082 c xgwts(15)=0.0 c xkwts(1)=0.022935322010529224963732008058970 c xkwts(2)=0.063092092629978553290700663189204 c xkwts(3)=0.104790010322250183839876322541518 c xkwts(4)=0.140653259715525918745189590510238 c xkwts(5)=0.169004726639267902826583426598550 c xkwts(6)=0.190350578064785409913256402421014 c xkwts(7)=0.204432940075298892414161999234649 c xkwts(8)=0.209482141084727828012999174891714 c xkwts(9)=0.204432940075298892414161999234649 c xkwts(10)=0.190350578064785409913256402421014 c xkwts(11)=0.169004726639267902826583426598550 c xkwts(12)=0.140653259715525918745189590510238 c xkwts(13)=0.104790010322250183839876322541518 c xkwts(14)=0.063092092629978553290700663189204 c xkwts(15)=0.022935322010529224963732008058970 C GAUSS-KRONROD ABSCISSAE c xgk(1)=0.991455371120812639206854697526329 c xgk(2)=0.949107912342758524526189684047851 c xgk(3)=0.864864423359769072789712788640926 c xgk(4)=0.741531185599394439863864773280788 c xgk(5)=0.586087235467691130294144838258730 c xgk(6)=0.405845151377397166906606412076961 c xgk(7)=0.207784955007898467600689403773245 c xgk(8)=0.000000000000000000000000000000000 c xgk(9)=-0.207784955007898467600689403773245 c xgk(10)=-0.405845151377397166906606412076961 c xgk(11)=-0.586087235467691130294144838258730 c xgk(12)=-0.741531185599394439863864773280788 c xgk(13)=-0.864864423359769072789712788640926 c xgk(14)=-0.949107912342758524526189684047851 c xgk(15)=-0.991455371120812639206854697526329 idone=0 iterno=0 errcode=0 Open(7,file='input.dat',status='old') read(7,5)fname 5 format(a10) open(8,file=fname,status='unknown') Call Input(laglimit,elimit,distrib,p,rho,errrel,isym,maxrho, &minrho,empcdf,nempir,mean,variance,minfeas,rhoincr,remrho, &ordrho,dparam,z,prob) lagsrem=p Call InitGrid(isym,minrho,maxrho,gspace1, &dimgrid,rgrid,case) Dowhile(idone.eq.0) iterno=iterno+1 if(distrib.lt.10)then Call Outer(rgrid,dimgrid,distrib,tolin,tolout,dparam, & mean,variance,p,maxsubs,minfeas,result) else call outerd(empcdf,z,nempir,rgrid,dimgrid,mean, &variance,result) endif c write(8,*)iterno,dimgrid c do i=1,dimgrid c write(8,*)rgrid(i),result(i) c enddo Call Update(rgrid,dimgrid,gspace1,gspace2,case,minfeas,result, &rhoincr,ordrho,errrel,lagsrem,p,remrho,iterno,igrid, &estrho,esterr,r,idone,a,b,rhoa,rhob) Enddo if(p.ne.1)then Call YW(p,r,alpha) Call Station(alpha,p,istat) else alpha(1)=r(1) istat=0 endif if(istat.eq.1)then call output(6,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) else call output(errcode,laglimit,maxsubs,elimit,p,r,estrho, &esterr,alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif close(8) STOP END Block Data WEIGHTS implicit real*8 (x) common/gkwts/xgwts(15),xkwts(15) C GAUSS AND KRONROD WEIGHT VALUES data xgwts(1),xgwts(2),xgwts(3),xgwts(4),xgwts(5), & xgwts(6),xgwts(7),xgwts(8),xgwts(9),xgwts(10), & xgwts(11),xgwts(12),xgwts(13),xgwts(14),xgwts(15), & xkwts(1),xkwts(2),xkwts(3),xkwts(4),xkwts(5), & xkwts(6),xkwts(7),xkwts(8),xkwts(9),xkwts(10), & xkwts(11),xkwts(12),xkwts(13),xkwts(14),xkwts(15) / & 0.0, & 0.129484966168869693270611432679082, & 0.0, & 0.279705391489276667901467771423780, & 0.0, & 0.381830050505118944950369775488975, & 0.0, & 0.417959183673469387755102040816327, & 0.0, & 0.381830050505118944950369775488975, & 0.0, & 0.279705391489276667901467771423780, & 0.0, & 0.129484966168869693270611432679082, & 0.0, & 0.022935322010529224963732008058970, & 0.063092092629978553290700663189204, & 0.104790010322250183839876322541518, & 0.140653259715525918745189590510238, & 0.169004726639267902826583426598550, & 0.190350578064785409913256402421014, & 0.204432940075298892414161999234649, & 0.209482141084727828012999174891714, & 0.204432940075298892414161999234649, & 0.190350578064785409913256402421014, & 0.169004726639267902826583426598550, & 0.140653259715525918745189590510238, & 0.104790010322250183839876322541518, & 0.063092092629978553290700663189204, & 0.022935322010529224963732008058970 / end block data ABSC C GAUSS-KRONROD ABSCISSAE implicit real*8 (x) common/gkpnts/xgk(15) data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), &xgk(9),xgk(10),xgk(11),xgk(12),xgk(13),xgk(14),xgk(15) / & 0.991455371120812639206854697526329, & 0.949107912342758524526189684047851, & 0.864864423359769072789712788640926, & 0.741531185599394439863864773280788, & 0.586087235467691130294144838258730, & 0.405845151377397166906606412076961, & 0.207784955007898467600689403773245, & 0.000000000000000000000000000000000, & -0.207784955007898467600689403773245, & -0.405845151377397166906606412076961, & -0.586087235467691130294144838258730, & -0.741531185599394439863864773280788, & -0.864864423359769072789712788640926, & -0.949107912342758524526189684047851, & -0.991455371120812639206854697526329 / end SUBROUTINE INPUT(laglimit,elimit,distrib,p,rho,errrel,isym, &maxrho,minrho,empcdf,nempir,mean,variance,minfeas,rhoincr, &remrho,ordrho,dparam,z,prob) c Reads user input from the file "input.dat". Checks input for errors. C Marne C. Cario C ARTA Subroutines C Updates: 8/95,9/16/95,10/10/95,11/8/95,11/15/95,11/16/95,11/20/95(Tested) c 1/3/96, 1/7/96, 1/15/96 c Release 1.1 6/3/96 c Subroutines called: output, chkpd, chkbnd, gamfxn, emprcal, echoin, correl, c order C Called by: main program c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c laglimit integer maximum number of lags that the c program can match c elimit integer maximum number of points for an empirical c distribution c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c distrib integer code corresponding to the marginal distribution c p integer number of autocorrelation lags to match c rho double-prec. vector px1 vector of autocorrelations to match c errrel double-prec. vector px1 vector of desired relative errors for the c elements of rho c isym integer 0 if the ARTA distribution is not symmetric c 1 if the ARTA distribution is symmetric c maxrho double precision largest ARTA autocorrelation c minrho double precision smallest ARTA autocorrelation c empcdf double-prec. vector vector of points in the empirical distribution c (if applicable) c nempir integer number of points in the empirical distribution c (if applicable) c mean double precision mean of the ARTA distribution c variance double precision variance of the ARTA distribution c minfeas double precision minimum feasible bivariate correlation c for distrib c rhoincr double-prec. vector px1 vector of the ARTA-process a.c.'s to match c in increasing order c remrho integer vector vector of positions in rhoincr corresponding c to ARTA-process a.c.'s that remain to be c matched c ordrho integer vector px1 vector of the lags corresponding to the c members of rhoincr c dparam double-prec. vector vector of ARTA-distribution parameters c ***************************************************************************** c *********** VARIABLE ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c ipd integer 0 if correlation matrix is pos. definite c 1 if correlation matrix is not pos. definite c ifeas integer 0 if all correlations are within the bounds c 1 if one or more correlations are out of bounds c icor integer 1 if the user wants empirical time-series c data read in and a.c.'s calculated c 0 if not c ***************************************************************************** integer laglimit,distrib,p,isym,elimit,nempir,ipd,ifeas,icor, &remrho(5),ordrho(5),ifault real*8 rho(5),errrel(5),maxrho,minrho,dparam(5), &empcdf(300),mean,variance,rhoincr(5),minfeas,prob(300), &z(300),probab,vnorm2,pr,x,y,alngam,mu,sigma2,sumsq, &r(5),estrho(5),esterr(5),alpha(5) external output,chkbnd,chkpd,emprcal,correl,order nempir=2 read(7,500)icor if(icor.eq.1)then read(7,501)nobs read(7,502)nlags call correl(nobs,nlags) close(8) close(7) stop endif read(7,510)distrib 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 close(7) call output(1,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif read(7,520)p if(p.gt.laglimit)then close(7) call output(2,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif do i=1,p read(7,530)rho(i),errrel(i) enddo call chkpd(rho,p,ipd) if(ipd.eq.1)then close(7) call output(5,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif call order(rho,p,rhoincr,ordrho,remrho) minrho=rhoincr(1) maxrho=rhoincr(p) goto(10,20,30,40,50,60,70,80,90,100,110),distrib c NORMAL DISTRIBUTION c dparam(1) = MEAN c dparam(2) = VARIANCE 10 read(7,540)dparam(1),dparam(2) close(7) if(dparam(2).le.0.0)then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif mean=dparam(1) variance=dparam(2) call chkbnd(1,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif c NORMAL MARGINALS ARE MODELED DIRECTLY AS AN AR(p) PROCESS. c CALL SUBROUTINE TO CALCULATE AND OUTPUT THE PARAMETERS OF THE AR(p) PROCESS call normal(p,rho,dparam,errrel) stop c STUDENT'S t DISTRIBUTION c dparam(1) = DEGREES OF FREEDOM 20 read(7,550)nu close(7) dparam(1)=nu if(dparam(1).lt.3)then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif mean=0.0 variance=dparam(1)/(dparam(1)-2.0) call chkbnd(2,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif isym=1 goto 1000 c CONTINUOUS UNIFORM DISTRIBUTION c dparam(1) = LEFT ENDPOINT c dparam(2) = RIGHT ENDPOINT 30 read(7,540)dparam(1),dparam(2) close(7) mean=(dparam(1)+dparam(2))/2.0 variance=(dparam(2)-dparam(1))*(dparam(2)-dparam(1))/12.0 call chkbnd(3,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif c CALL A SUBROUTINE THAT SOLVES FOR THE AR_PROCESS A.C.'s USING c THE CLOSED-FORM SOLUTION IN LI AND HAMMOND (1975) call uniform(p,rho,dparam,errrel) stop c EXPONENTIAL DISTRIBUTION c dparam(1) = RATE 40 read(7,560)dparam(1) close(7) if(dparam(1).le.0.0)then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif mean=1.0/dparam(1) variance=1.0/(dparam(1)*dparam(1)) call chkbnd(4,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif goto 1000 c GAMMA DISTRIBUTION c dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 50 read(7,540)dparam(1),dparam(2) close(7) if((dparam(1).le.0.0).or.(dparam(2).le.0.0))then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif mean=dparam(1)*dparam(2) variance=dparam(2)*dparam(1)*dparam(1) call chkbnd(5,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif goto 1000 c WEIBULL DISTRIBUTION C dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 60 read(7,540)dparam(1),dparam(2) close(7) if((dparam(1).le.0.0).or.(dparam(2).le.0.0))then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif x=1.0+1.0/dparam(2) x=alngam(x,ifault) mean=dparam(1)*dexp(x) x=1.0+2.0/dparam(2) x=alngam(x,ifault) variance=dparam(1)*dparam(1)*dexp(x)-mean*mean call chkbnd(6,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif goto 1000 c LOGNORMAL DISTRIBUTION c dparam(1)=scale parameter mu c dparam(2)=square of the shape parameter sigma 70 read(7,540)mean,variance close(7) y=2*dlog(mean) x=dlog(variance+mean*mean) dparam(1)=y-x/2.0 dparam(2)=x-y if(dparam(2).le.0.0)then call output(3,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif call chkbnd(7,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif 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 80 do i=1,4 read(7,560)dparam(i) enddo close(7) sigma2=dble(1.0)/(dparam(2)*dparam(2)) mu=dble(-1.0)*dparam(1)/dparam(2) mean=dparam(3)+dparam(4)*dexp(sigma2/dble(2.0))*dsinh(mu) variance=dparam(4)*dparam(4)*(dexp(sigma2)-dble(1.0))* &(dble(1.0)+dcosh(dble(2.0)*mu)*dexp(sigma2))/dble(2.0) call chkbnd(8,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif 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 90 do i=1,4 read(7,560)dparam(i) enddo close(7) c STILL NEED MEAN AND VARIANCE FOR THIS DISTRIBUTION print*, ' Johnson bounded distribution not yet avaliable' stop call chkbnd(9,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif c mean= c variance= goto 1000 c EMPIRICAL CDF c dparam(1) = NUMBER OF POINTS IN CDF 100 read(7,550)nempir close(7) dparam(1)=nempir if(nempir.gt.elimit)then call output(9,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif call emprcal(nempir,empcdf) mean=0.0 sumsq=0.0 x=1.0/nempir do i=1,nempir mean=mean+empcdf(i)*x sumsq=sumsq+empcdf(i)*empcdf(i)*x enddo variance=sumsq-mean*mean call chkbnd(10,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif pr=1.0/nempir probab=0.0 do i=1,nempir-1 prob(i)=pr probab=probab+prob(i) z(i)=vnorm2(probab,ifault) enddo prob(nempir)=pr goto 1000 c DISCRETE DISTRIBUTION WITH FINITE SUPPORT c dparam(1) = NUMBER OF POINTS IN CDF 110 read(7,550)nempir dparam(1)=nempir do i=1,nempir read (7,*) empcdf(i),prob(i) if(prob(i).le.0.0)then write(8,112) 112 format('Input probability must be between 0.0 and 1.0') stop endif enddo close(7) if(nempir.gt.elimit)then call output(9,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif call bubsort3(empcdf,prob,nempir) mean=0.0 sumsq=0.0 probab=0.0 do i=1,nempir probab=probab+prob(i) mean=mean+empcdf(i)*prob(i) sumsq=sumsq+empcdf(i)*empcdf(i)*prob(i) enddo if(dabs(probab-dble(1.0)).gt.0.0001)then write(8,113) 113 format(/,'ERROR: Probabilities do not add up to 1.') close(8) stop endif variance=sumsq-mean*mean call chkbnd(11,minrho,maxrho,dparam,empcdf,nempir,prob,mean, &variance,minfeas,ifeas) if(ifeas.eq.1)then call output(4,laglimit,maxsubs,elimit,p,r,estrho,esterr, &alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif probab=0.0 do i=1,nempir-1 probab=probab+prob(i) z(i)=vnorm2(probab,ifault) enddo 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(distrib,dparam,p,rho,errrel,empcdf,prob, &nempir,mean,variance,minfeas) 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*8 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 echoin(distrib,dparam,p,rho, &errrel,empcdf,prob,nempir,mean,variance,minfeas) C Echoes user input to output file "ARTA.out" c Marne C. Cario c ARTA subroutines c Updates: 10/10/95, 11/8/95, 11/20/95 (Tested), 1/3/96, 1/7/96, c 2/9/96 c Subroutines called: None c Calling program: input integer distrib,p,nempir real*8 dparam(5),rho(p),errrel(p),prob(nempir), &empcdf(nempir),mean,variance,minfeas write(8,5) 5 format('********** USER INPUT **********',/) goto(10,20,30,40,50,60,70,80,90,100,110),distrib 10 write(8,15) 15 format('ARTA Distribution: Normal') write(8,16)dparam(1),dparam(2) 16 format('Mean: ',f10.5,/,'Variance: ',f10.5,/) goto 150 20 write(8,25) 25 format('ARTA Distribution: t') nu=dparam(1) write(8,26)nu 26 format('Degrees of freedom: ',i4,/) goto 150 30 write(8,35) 35 format('ARTA Distribution: Continuous Uniform') write(8,36)dparam(1),dparam(2) 36 format('Left Endpoint: ',f10.5,/,'Right Endpoint: ',f10.5,/) goto 150 40 write(8,45) 45 format('ARTA Distribution: Exponential') write(8,46)dparam(1) 46 format('Rate: ',f10.5,/) goto 150 50 write(8,55) 55 format('ARTA Distribution: Gamma') write(8,56)dparam(1),dparam(2) 56 format('Scale Parameter: ',f10.5,/,'Shape Parameter: ', &f10.5,/) goto 150 60 write(8,65) 65 format('ARTA Distribution: Weibull') write(8,66)dparam(1),dparam(2) 66 format('Scale Parameter: ',f10.5,/,'Shape Parameter: ', &f10.5,/) goto 150 70 write(8,75) 75 format('ARTA Distribution: Lognormal') write(8,76)mean,variance 76 format('Mean: ',f10.5,/,'Variance: ', &f10.5,/) goto 150 80 write(8,85) 85 format('ARTA Distribution: Johnson Unbounded') write(8,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(8,95) 95 format(' ARTA Distribution: Johnson Bounded') write(8,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(8,105) 105 format('ARTA Distribution: Empirical') nempir=dparam(1) write(8,106)nempir 106 format('Number of points: ',i4,/) goto 150 110 write(8,115) 115 format('ARTA Distribution: Discrete with finite support') nempir=dparam(1) write(8,116)nempir 116 format('Number of points: ',i4,/) write(8,117) 117 format(3x,'Value',7x,'Probability',/,3x,'*****',7x,'***********', &/) do i=1,nempir write(8,118)empcdf(i),prob(i) enddo 118 format(f10.5,5x,f7.5) goto 150 150 write(8,155)p 155 format('Number of ARTA autocorrelation lags to match: ',i2,/) write(8,160) 160 format('Lag',3x,'Desired Autocorrelation', &3x,'Desired Relative Error', &/,' ***',3x,'***********************',3x,'**********************') do i=1,p write(8,170)i,rho(i),errrel(i) 170 format(1x,i2,7x,f15.13,10x,f7.5) enddo write(8,175)minfeas 175 format(/,'Minimum feasible bivariate correlation is: ',f8.5,/) write(8,180) 180 format(/,' **********************************') return end SUBROUTINE CHKPD(rho,p,ipd) c Checks whether the correlation matrix defined by rho is positive c definite. Based upon Lemma 3.3.11 and the superdiagonalization algorithm on c p. 95 and pp. 96-98, respectively, of "Nonlinear Programming Theory and c Algorithms" by Bazaraa, Sherali, and Shetty (John Wiley & Sons, 1983). C Note: The algorithm may only be used for symmetric matrices. c Marne C. Cario c ARTA Subroutines c Updates: 11/9/95, 11/16/95, 11/17/95 (Tested) c Subroutines called: subchkpd c Calling program: input c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c rho double-prec. vector px1 vector of autocorrelations to match c p integer number of autocorrelation lags to match c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c ipd integer 0 if correlation matrix is pos. definite c 1 if correlation matrix is not pos. definite c ***************************************************************************** integer p,ipd,i,ifin,n,m real*8 h(6,6),rho(p) external subchkpd n=p+1 c CONSTRUCT THE CORRELATION MATRIX H IMPLIED BY RHO do i=2,n do j=1,i-1 k=i-j h(i,j)=rho(k) enddo enddo do i=1,n h(i,i)=1.0 enddo do i=1,n-1 do j=i+1,n h(i,j)=h(j,i) enddo enddo C CODE TO WRITE OUT H TO FILE 7 FOR DEBUGGING PURPOSES c do i=1,n c do j=1,n c write(7,*)i,j,h(i,j) c enddo c enddo C IMPLEMENT RECURSIVE ROUTINE TO DETERMINE IF H IS POSITIVE DEFINITE m=0 ifin=0 dowhile(ifin.eq.0) m=m+1 call subchkpd(h,n,m,ifin,ipd) enddo RETURN END SUBROUTINE SUBCHKPD(h,n,m,ifin,ipd) C Checks the current submatrix of H to determine if it is NONpositive c definite. If it is not determined to be NONpositive definite, then c the submatrix is further row-reduced and the routine may be called c again for this submatrix. When the submatrix is of dimension 2x2, then c a determination can be made quickly. For completeness, the test for c a 1x1 matrix is included, but shouldn't be needed for the ARTA program. c Marne C. Cario c ARTA Subroutines c Updates: 11/9/95, 11/17/95 (Tested) c Subroutines called: None c Calling program: chkpd integer n,m,ifin,ipd,dim,rowno,j,jj real*8 h(6,6),factor C CODE TO WRITE OUT H TO FILE 7 FOR DEBUGGING PURPOSES c do i=1,n c do j=1,n c write(7,*)i,j,h(i,j) c enddo c enddo c Calculate dimension of the submatrix dim=n-m+1 c If the submatrix is 1x1 or 2x2, then test for p.d. and we are done. if(dim.eq.1)then c write(7,*)'Dim = 1' ifin=1 if(h(n,n).gt.0.0)then ipd=0 else ipd=1 endif goto 100 elseif(dim.eq.2)then c write(7,*)'Dim = 2' diff=h(n-1,n-1)*h(n,n) - h(n-1,n)*h(n,n-1) c write(7,*)diff ifin=1 if((h(n-1,n-1).gt.0.0).and.(h(n,n).gt.0.0).and.(diff.gt.0.0))then ipd=0 else ipd=1 endif goto 100 endif c If the submatrix is 3x3 or larger, then we perform a test on the diagonal c elements. If any diagonal element is <= 0, then H is NOT positive definite c and we are done. Otherwise, row-reduce the submatrix to get a new submatrix c which is one less in dimension than the previous submatrix. do j=m,n if(h(j,j).le.0.0)then c write(7,*)'diag. element le 0' ifin=1 ipd=1 goto 100 endif enddo c write(7,*)'row reducing' do j=1,n-m rowno=m+j factor=h(rowno,m)/h(m,m) do jj=m,n h(rowno,jj)=h(rowno,jj)-factor*h(m,jj) enddo enddo ifin=0 100 return end SUBROUTINE CHKBND(distrib,minrho,maxrho,dparam, &empcdf,nempir,prob,mean,variance,minfeas,ifeas) c Checks whether all of the elements of rho are within the minimum and c maximum feasible bivariate correlations for distrib c Marne C. Cario c ARTA subroutines c Updates: 8/31/95, 11/8/95, 11/16/95, 11/20/95, 1/3/96, 1/7/96 c 2/8/96, 2/9/96 (Tested) c Subroutines called: None c Called by: Input c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c distrib integer code corresponding to the marginal distribution c minrho double precision smallest ARTA autocorrelation c maxrho double precision largest ARTA autocorrelation c dparam double-prec. vector vector of parameters for the ARTA distribution c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c minfeas double precision minimum feasible bivariate correlation c for distrib c ifeas integer 0 if all correlations are within the bounds c 1 if one or more correlations are out of bounds c ***************************************************************************** integer ifeas,distrib,nempir real*8 dparam(5),minfeas,minrho,maxrho,expmin,pi, &empcdf(nempir),prob(nempir),mean,variance pi=dacos(dble(-1.0)) expmin=1.0-pi*pi/6.0 ifeas=0 if(maxrho.gt.1.0)then ifeas=1 goto 300 endif c write(8,5)mean,variance c 5 format('mean=',f10.3,/,'variance=',f10.3) goto (10,10,10,20,30,30,30,30,30,40,50),distrib C Minimum bivariate correlation is -1.0 for this distribution c (Normal, t, and uniform) 10 minfeas=-1.0 goto 210 C Exponential 20 minfeas=expmin goto 210 C Continuous distribution c Calculate int_{0}^{1} Finv(u)Finv(1-u) du 30 call rhomin(distrib,dparam,minfeas) goto 200 c Empirical cdf - Use CHR results 40 call erhomin(empcdf,nempir,minfeas) c write(8,45)mean,variance,minfeas c 45 format('mean=',f10.3,/,'variance=',f10.3,/,'minfeas=',f10.3) goto 200 c Discrete with finite support - Use CHR results 50 call drhomin(empcdf,prob,nempir,minfeas) 200 minfeas=(minfeas-mean*mean)/variance 210 if(minrho.lt.minfeas)then ifeas=1 endif 300 RETURN END SUBROUTINE ORDER(rho,p,rhoincr,ordrho,remrho) c Orders the ARTA-process autocorrelations and puts them in c increasing order in rhoincr and puts their corresponding lags c in ordrho. Also, initializes the remrho vector. c Marne C. Cario c ARTA Subroutines c Updates: 11/8/95, 11/17/95 (Tested) c Subroutines Called: bubsort2 c Calling program: input integer p,ordrho(p),remrho(p) real*8 rho(p),rhoincr(p) do i=1,p remrho(i)=i ordrho(i)=i rhoincr(i)=rho(i) enddo call bubsort2(rhoincr,ordrho,p) C WRITE STATEMENTS FOR DEBUGGING c write(7,*)' rhoincr ordrho remrho' c do i=1,p c write(7,*)rhoincr(i),ordrho(i),remrho(i) c enddo return end SUBROUTINE BUBSORT2(data,attrib,n) c Does bubble sort on array of length n and carries an integer c attribute along with the sorted array. Sorts from low to high. integer j,n,attrib(n),temp2 real*8 data(n),temp1 C WRITE STATEMENTS FOR DEBUGGING c write(7,*)' rhoincr ordrho ' c do i=1,n c write(7,*)'got there' c write(7,*)data(i),attrib(i) c enddo 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 OUTPUT(errcode,laglimit,maxsubs,elimit,p,r, &estrho,esterr,alpha,minfeas,distrib,rho,dparam,mean, &variance,nempir,empcdf,prob) c Writes program output to the file "ARTA.out" c Marne C. Cario c ARTA Subroutines c Updates: 11/8/95, 11/16/95, 11/20/95 (Tested), 1/3/96, 1/7/96, c 2/9/96, 2/19/96 c Subroutines called: none c Calling programs: input,outer,inner,main integer errcode,laglimit,maxsubs,elimit,p,distrib, &nempir real*8 r(p),estrho(p),esterr(p),alpha(p),minfeas, &rho(p),dparam(5),mean,variance,prob(nempir),empcdf(nempir) write(8,5) 5 format(//,' ********** ARTA PROGRAM OUTPUT **********',/) if(errcode.eq.0)then goto 200 else goto(10,20,30,40,50,60,70,80,90),errcode endif 10 write(8,15) 15 format(' ERROR: Distribution not available. Please refer to &the',/,9x,'file readme.txt for the available distributions and', &/,9x,'their corresponding code numbers') close(8) STOP 20 write(8,25)laglimit 25 format(' ERROR: Maximum number of autocorrelation lags is ', &i2,'.',/,9x,'Please refer to the file readme.txt for ', &/,9x,'instructions on changing the maximum number of lags.') close(8) STOP 30 write(8,35) 35 format(' ERROR: Parameter infeasible for the chosen ARTA &distribution.', &/,9x,'Please refer to the file readme.txt for constraints on ', &/,9x,'distribution parameters.') close(8) STOP 40 write(8,45)minfeas 45 format(' ERROR: Autocorrelation out of bounds for one or more &lags.',/,' Minimum feasible autocorrelation is ',f9.6) close(8) STOP 50 write(8,55) 55 format(' ERROR: The correlation matrix implied by the ARTA', &/,9x,'autocorrelation structure is not positive definite.') close(8) STOP 60 write(8,65) 65 format(' ERROR: The AR output process is not stationary.') goto 200 STOP 70 write(8,75)maxsubs 75 format(' ERROR: Maximum number ',i3,' of subintervals exceeded', &/,9x,'for the outer integral.') close(8) STOP 80 write(8,85)maxsubs 85 format(' ERROR: Maximum number ',i3,' of subintervals exceeded', &/,9x,'for the inner integral.') close(8) STOP 90 write(8,95)elimit 95 format(' ERROR: Limit on number of points in cdf', &/,9x,'exceeded.') close(8) STOP 200 write(8,205) 205 format(/,' Lag',10x,'ARTA Correl.',10x,'AR Correl.', &10x,'Rel. Error',/) do j=1,p write(8,210) j,estrho(j),r(j),esterr(j) 210 format(2x,i2,9x,f8.6,14x,f9.6,14x,f8.6) enddo write(8,215) 215 format(//,' i ',10x,' Alpha(i)') do j=1,p write(8,220) j,alpha(j) 220 format(i2,10x,f10.4) enddo C WRITE INPUT FILE FOR ARTAGEN open(15,file='input.gen',status='unknown') write(15,500) 500 format('obs.ag'/,'stats.ag') write(15,510)distrib 510 format(i2) write(15,520)p 520 format(i2) ii=1000 write(15,530)ii 530 format(i4) do i=1,p write(15,540)r(i) 540 format(f9.5) enddo do i=1,p write(15,540)rho(i) enddo goto(610,620,630,640,650,660,670,680,690,700,710),distrib c NORMAL DISTRIBUTION c NOTE: Should never get here. c dparam(1) = MEAN c dparam(2) = VARIANCE 610 write(15,615)dparam(1),dparam(2) 615 format(f10.5,/,f10.5) goto 1000 c STUDENT'S t DISTRIBUTION c dparam(1) = DEGREES OF FREEDOM 620 ii=int(dparam(1)) write(15,625)ii 625 format(i4) goto 1000 c CONTINUOUS UNIFORM DISTRIBUTION c NOTE: Should never get here c dparam(1) = LEFT ENDPOINT c dparam(2) = RIGHT ENDPOINT 630 write(15,635)dparam(1),dparam(2) 635 format(f10.5,/,f10.5) goto 1000 c EXPONENTIAL DISTRIBUTION c dparam(1) = RATE 640 write(15,645)dparam(1) 645 format(f10.5) goto 1000 c GAMMA DISTRIBUTION c dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 650 write(15,655)dparam(1),dparam(2) 655 format(f10.5,/,f10.5) goto 1000 c WEIBULL DISTRIBUTION C dparam(1) = SCALE PARAMETER c dparam(2) = SHAPE PARAMETER 660 write(15,665)dparam(1),dparam(2) 665 format(f10.5,/,f10.5) goto 1000 c LOGNORMAL DISTRIBUTION c dparam(1)=scale parameter mu c dparam(2)=square of the shape parameter sigma 670 write(15,675)mean,variance 675 format(f10.5,/,f10.5) 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 680 do i=1,4 write(15,685)dparam(i) 685 format(f10.5) enddo 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 690 do i=1,4 write(15,695)dparam(i) 695 format(f10.5) enddo goto 1000 c EMPIRICAL CDF c dparam(1) = NUMBER OF POINTS IN CDF 700 write(15,705)nempir 705 format(i4) goto 1000 c DISCRETE DISTRIBUTION WITH FINITE SUPPORT c dparam(1) = NUMBER OF POINTS IN CDF 710 write(15,715)nempir 715 format(i4) do i=1,nempir write(15,717)empcdf(i),prob(i) 717 format(f10.5,5x,f6.5) enddo 1000 write(15,1010) 1010 format(//) close(15) RETURN END SUBROUTINE InitGrid(isym,minrho,maxrho,gspace1, &dimgrid,rgrid,case) c Calculates the grid points (r-values) for which the double c integral will be calculated initially c Marne C. Cario c ARTA subroutine c Updates: 8/31/95, 9/21/95, 10/10/95, 11/8/95, 11/16/95, 11/20/95 (Tested) c Subroutines called: None c Called by: Main program c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c isym integer 0 if distrib is not symmetric, 1 if symmetric c minrho double precision smallest ARTA autocorrelation c maxrho double precision largest ARTA autocorrelation c gspace1 integer 1.0/gspace1 is the initial spacing between c grid points c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c dimgrid integer number of r-values for which the integration c will be carried out initially c rgrid double-prec. vector dimgridx1 vector of the r-values for which the c integration will be carried out initially c case integer 1 if all ARTA-process a.c.'s are positive c 2 if all ARTA-process a.c.'s are negative c 3 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is not c symmetric c 4 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is c symmetric c ***************************************************************************** integer dimgrid,isym,gspace1,n,case real*8 maxrho,minrho,rgrid(50) n=gspace1-1 if(minrho.gt.0.0)then case=1 elseif(maxrho.lt.0.0)then case=2 elseif((minrho.lt.0.0).and.(maxrho.gt.0.0).and.(isym.ne.1))then case=3 else case=4 endif c If the ARTA autocorrelations are all negative, then only consider negative c r-values. Note: This will rarely be feasible. c If the ARTA distribution is symmetric or if the ARTA autocorrelations are all c positive, then only consider positive r-values c Otherwise, consider both positive and negative r-values. if(case.eq.2)then do i=1,n rgrid(i)=-real((gspace1-i))/gspace1 enddo dimgrid=n elseif((case.eq.1).or.(case.eq.4))then do i=1,n rgrid(i)=real(i)/gspace1 enddo dimgrid=n else do i=1,n rgrid(i)=-real((gspace1-i))/gspace1 rgrid(i+n)=real(i)/gspace1 enddo dimgrid=2*n endif RETURN END SUBROUTINE GKPOINTS(l,r,gkpts) c Given an interval (l,r) on the real line, calculate the points c at which the integrand evaluation will take place, based upon c the technique of a 15-point Gauss-Kronrod quadrature rule. c Marne C. Cario c ARTA subroutine c Updates: 8/31/95, 11/8/95, 11/16/95 c Subroutines called: None c Calling programs: outer, inner c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c l double precision left endpoint of the interval over which the c quadrature calculation will take place c r double precision right endpoint of the interval over which the c quadrature calculation will take place c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c gkpoints double-prec. vector 15x1 vector of the points at which the c integrand evaluation will take place c ***************************************************************************** implicit real*8 (x) real*8 l,r,gkpts(15),center,halfwid common/gkpnts/xgk(15) center=(l+r)/2.0 halfwid=(r-l)/2.0 do i=1,15 gkpts(i)=center+halfwid*xgk(i) enddo RETURN END SUBROUTINE PUTHEAP(maxsubs,error,interval,eheap,iheap,last) C Subroutine to put an error and its interval number in a c partially-ordered binary heap (top element of heap has maximum error) c Marne C. Cario c ARTA Subroutines c Updates: 9/6/95, 10/10/95, 11/8/95, 11/16/95, 11/20/95 (Tested), c 2/21/96 c Subroutines Called: None c Called by: outer, inner c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c error double precision value of the error that is to be placed on the c heap c interval integer subinterval number corresponding to error c eheap double-prec. vector heap array of maximum errors for each c subinterval that has not yet been bisected c iheap integer array array of subinterval numbers corresponding to c the members of eheap c last integer pointer to the last element of eheap c ***************************************************************************** c *********** OUTPUT ********************************************************* c eheap, iheap, and last get updated and output c ***************************************************************************** c *********** VARIABLES *************************************************** c VARIABLE TYPE DESCRIPTION c ******** **** *********** c newnode integer node in which the error value that is to be c placed n the heap currently resides c parent integer parent node of newnode c iflag integer = 0 if the error is not yet in the correct node c = 1 if the error is in the correct node c ***************************************************************************** integer maxsubs,last,interval,iflag,itemp,newnode,parent real*8 error,eheap(maxsubs),tempe integer iheap(maxsubs) C INCREMENT last last=last+1 C TEST TO SEE IF THE HEAP IS FULL c if(last.gt.500)then c print*,'ERROR - Exceeded dimension of heap array.' c stop c endif C IF THE HEAP IS NOT FULL, THEN PROCEED TO PUT THE NEW ERROR IN THE C HEAP; PLACE IT AT THE BOTTOM OF THE HEAP INITIALLY eheap(last)=error iheap(last)=interval newnode=last C SWAP THE ERROR WITH THE PARENT NODE UNTIL THE PARENT NODE HAS A C LARGER ERROR OR UNTIL THE ERROR IS AT THE TOP OF THE HEAP (PARENT=0) iflag=0 dowhile(iflag.ne.1) parent=newnode/2.0 if(parent.gt.0)then if(eheap(newnode).gt.eheap(parent))then tempe=eheap(parent) itemp=iheap(parent) eheap(parent)=eheap(newnode) iheap(parent)=iheap(newnode) eheap(newnode)=tempe iheap(newnode)=itemp newnode=parent else iflag=1 endif else iflag=1 endif enddo return end SUBROUTINE REMHEAP(maxsubs,eheap,iheap,last,ierr) C Subroutine to remove the maximum error from the top of the heap and c rearrange the heap so that it is still a partially-ordered, packed c heap c Marne C. Cario c ARTA Subroutines c Updates: 9/6/95, 11/8/95, 11/16/95, 11/20/95 (Tested), c 2/21/96 c Subroutines Called: None c Called by: outer, inner c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c eheap double-prec. vector heap array of maximum errors for each c subinterval that has not yet been bisected c iheap integer array array of subinterval numbers corresponding to c the members of eheap c last integer pointer to the last element of eheap c ***************************************************************************** c *********** OUTPUT ********************************************************* c ierr integer interval number corresponding to the largest c error c NOTE: eheap,iheap, and last get updated and output c ***************************************************************************** c *********** VARIABLES *************************************************** c VARIABLE TYPE DESCRIPTION c ******** **** *********** c node integer node in which the error value that was last on c the heap at the time the subroutine was called c currently resides c node1 integer first "child" of node c node2 integer second "child" of node c iflag integer = 0 if the error is not yet in the correct node c = 1 if the error is in the correct node c ***************************************************************************** integer maxsubs,last,iflag,itemp,node,node1,node2,ierr real*8 eheap(maxsubs),tempe integer iheap(maxsubs) C TEST TO MAKE SURE THAT THE HEAP IS NOT EMPTY if(last.eq.0)then print*,'ERROR: Trying to remove from an empty heap' stop endif C SET THE VALUE OF THE MAXIMUM-ERROR SUBINTERVAL NUMBER TO THE INTERVAL C NUMBER CORRESPONDING TO THE ELEMENT AT THE TOP OF THE HEAP ierr=iheap(1) c write(15,1000)ierr c 1000 format('ierr in remheap = ',i2) C IF THERE WAS MORE THAN ONE ELEMENT ON THE HEAP, THEN REARRANGE THE C HEAP TO ACCOUNT FOR "REMOVAL" OF THE TOP ELEMENT if(last.gt.1)then C PUT THE ELEMENT THAT WAS AT THE BOTTOM OF THE HEAP AT THE TOP eheap(1)=eheap(last) iheap(1)=iheap(last) last=last-1 node=1 C WHILE THE ERROR OF AT LEAST ONE OF ITS CHILDREN IS LARGER, SWAP C PLACES WITH THE "LARGEST-ERROR" CHILD iflag=0 node1=2.0*node node2=node1+1 dowhile(iflag.ne.1) if(eheap(node1).ge.eheap(node2))then if(eheap(node1).gt.eheap(node))then tempe=eheap(node) itemp=iheap(node) eheap(node)=eheap(node1) iheap(node)=iheap(node1) eheap(node1)=tempe iheap(node1)=itemp node=node1 node1=2.0*node node2=node1+1 if(node1.gt.last)then iflag=1 endif else iflag=1 endif else if(eheap(node2).gt.eheap(node))then tempe=eheap(node) itemp=iheap(node) eheap(node)=eheap(node2) iheap(node)=iheap(node2) eheap(node2)=tempe iheap(node2)=itemp node=node2 node1=2.0*node node2=node1+1 if(node1.gt.last)then iflag=1 endif else iflag=1 endif endif enddo else last=last-1 endif return end real*8 function bvnpdf(z1,z2,rho) c Calculates the value of the standard bivariate normal pdf with c correlation rho at the point (z1,z2) c Marne C. Cario, 11/9/95, 11/16/95 real*8 z1,z2,rho,pi,diff,term1,term2 pi=dacos(dble(-1.0)) diff=dble(1.0)-rho*rho term1=dble(1.0)/(dble(2.0)*pi*dsqrt(diff)) term2=dexp((-1.0/(2.0*diff))*(z1*z1-2.0*rho*z1*z2+z2*z2)) bvnpdf=term1*term2 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 STATION(alpha,p,istat) c Determines whether the AR process defined by alpha is stationary c Marne C. Cario c ARTA subroutine c Updates: 11/94, 8/31/95, 10/10/95, 11/17/95 (Tested) c Subroutines called: None c Calling program: Main c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c alpha double-prec. vector px1 vector of the AR parameters c p integer number of autocorrelation lags to match c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c istat integer 0 if AR process is stationary c 1 if AR process is not stationary c ***************************************************************************** integer istat,p real*8 alpha(p),a(6),b(6) do i=2,p+1 a(i)=alpha(i-1) enddo a(1)=-1.0 do k=p,1,-1 do j=1,k b(j)=a(1)*a(j)-a(k+1)*a(k+2-j) enddo if(b(1).le.0)then istat=1 goto 10 endif do j=1,k a(j)=b(j) enddo enddo istat=0 10 RETURN END SUBROUTINE UPDATE(rgrid,dimgrid,gspace1,gspace2,case,minfeas, &result,rhoincr,ordrho,errrel,lagsrem,p,remrho,iterno, &igrid,estrho,esterr,r,idone,a,b,rhoa,rhob) c Determines which ARTA-process autocorrelations have been achieved c to within the desired relative error. c For ARTA-process autocorrelations that have not been achieved, calculates c new grid points. c c Performs record-keeping functions c Marne C. Cario c ARTA Subroutines c Updates: 10/4/95, 10/5/95, 10/10/95, 11/8/95, 11/16/95, 1/3/96, 1/11/96, c 1/16/96 c Subroutines called: search, calcerr, procerr c Calling program: Main program c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c gspace1 integer number of grid points between bracketing c AR-process a.c.'s during the first iteration c gspace2 integer number of grid points between bracketing c AR-process a.c.'s during subsequent iterations c result double-prec. vector vector of the rho values corresponding to c each member of rgrid c p integer number of a.c. lags to match c rhoincr double-prec. vector px1 vector of ARTA-process a.c.'s to match c in increasing order c ordrho integer vector px1 vector of the lags corresponding to c the members of rhoincr c remrho integer vector vector of the positions in rhoincr c corresponding to the ARTA-process a.c.'s that c remain to be matched c errrel double-prec. vector px1 vector of desired relative errors for the c elements of rho c case integer 1 if all ARTA-process a.c.'s are positive c 2 if all ARTA-process a.c.'s are negative c 3 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is not c symmetric c 4 if ARTA-process a.c.'s are mixed in sign c and the ARTA-process distribution is c symmetric c minfeas double precision minimum feasible bivariate correlation c for distrib c iterno integer number of iterations so far c ***************************************************************************** c *********** VARIABLES THAT ARE UPDATED ******************************** c VARIABLE TYPE DESCRIPTION c ******** **** *********** c rgrid double-prec. vector vector of current grid points - gets updated c within subroutine c dimgrid integer number of elements of rgrid - gets updated c lagsrem integer number of ARTA-process autocorrelations that c have not been achieved within the specified c relative error - gets updated c r double-prec. vector px1 vector of the AR-process a.c.'s c that result in ARTA-process a.c.'s that are c closest to the desired value of those c considered so far c estrho double-prec. vector vector of the calculated values of the ARTA- c process a.c.'s corresponding to the c members of r c esterr double-prec. vector vector of the relative errors of estrho c a double-prec. vector vector of left brackets for the AR-process c grid points to be used during the next c iteration c b double-prec. vector vector of right brackets for the AR-process c grid points to be used during the next c iteration c rhoa double-prec. vector vector of ARTA-process a.c.'s corresponding c to a c rhob double-prec. vector vector of ARTA-process a.c.'s corresponding c to b c idone integer 0 if at least one ARTA-process a.c. has not been c matched c 1 if all ARTA-process a.c.'s have been matched c ***************************************************************************** c *********** VARIABLES ****************************************************** c VARIABLE TYPE DESCRIPTION c ******** **** *********** c pos integer position in ordrho of the current ARTA-process c a.c. lag c lftbrack double precision lower AR-process a.c. bracket for the current c lag c rtbrack double precision upper AR-process a.c. bracket for the current c lag c lvalue double precision value of the ARTA-process a.c. at lftbrack c rvalue double precison value of the ARTA-process a.c. at rtbrack c lftpos integer position in rgrid of lftbrack c rtpos integer position in rgrid of rtbrack c closerho double precision of rvalue and lvalue, the one that is closest c to the desired ARTA-process a.c. (for the c current lag) c closer double precison of rtbrack and lftbrack, the one corresponding c to closerho c error double precision relative error of closerho c newgrid double-prec. vector AR-process grid points to be used during the c next iteration (at the end of the subroutine, c rgrid gets set equal to newgrid) c dimnew integer dimension of newgrid c newlr integer current count of the number of ARTA-process c lags that remain to be matched (lagsrem gets c set equal to this at the end of the subroutine) c (counts down from lagsrem) c countrem integer current count of the number of ARTA-process c lags that remain to be matched (counts up from c 0) c llimit integer lowest value of the position of the AR-process c a.c. for the current lag c rlimit integer highest value of the position of the AR-process c a.c. for the current lag c llimr double precision lower bound on the AR-process a.c. for the c current lag c rlimr double precision upper bound on the AR-process a.c. for the c current lag c llimrho double precision lower bound on the ARTA-process a.c for the c current lag c rlimrho double precision upper bound on the ARTA-process a.c. for the c current lag c yrho double precision Current ARTA-process a.c. C isign integer -1 if the ARTA-process a.c. is negative and c the ARTA distribution is symmetric c pcount integer 0 if the ARTA-process a.c.'s are mixed in sign c and the first positive a.c. has not yet been c processed (during the first iteration only) c 1 if the ARTA-process a.c.'s are mixed in sign c and the first positive a.c. has been processed c (during the first iteration only) c ***************************************************************************** integer p,dimgrid,gspace1,gspace2,ordrho(p),lagsrem,remrho(p), &idone,newlr,countrem,dimnew,case,llimit,rlimit,lftpos,rtpos,pos, &pcount,isign,remrho2(5),igrid(p) real*8 rgrid(dimgrid),result(dimgrid),rhoincr(p),errrel(p), &estrho(p),esterr(p),r(p),minfeas,llimr,rlimr,llimrho,rlimrho, &yrho,abserror,closer,closerho,error,erel,newgrid(50),a(p), &b(p),rhoa(p),rhob(p),lftbrack,rtbrack,lvalue,rvalue external search,calcerr,procerr newlr=lagsrem countrem=0 dimnew=0 do i=1,lagsrem remrho2(i)=remrho(i) enddo if(iterno.eq.1)then if(case.eq.1)then llimit=0 rlimit=gspace1 llimr=0.0 rlimr=1.0 llimrho=0.0 rlimrho=1.0 elseif(case.eq.2)then llimit=0 rlimit=gspace1 llimr=-1.0 rlimr=0.0 llimrho=minfeas rlimrho=0.0 endif if(case.lt.3)then lftpos=0 rtpos=gspace1 do i=1,lagsrem yrho=rhoincr(remrho(i)) pos=remrho(i) rtpos=gspace1 call search(yrho,llimit,rlimit,llimr,rlimr,llimrho,rlimrho, & lftpos,rtpos,case,isign,rgrid,dimgrid,result, & lftbrack,rtbrack,lvalue,rvalue) call calcerr(lvalue,rvalue,yrho,lftbrack,rtbrack,closer, & closerho,error) r(ordrho(pos))=closer estrho(ordrho(pos))=closerho esterr(ordrho(pos))=error abserror=abs(error) erel=errrel(ordrho(pos)) call procerr(abserror,erel,remrho2,lftbrack,rtbrack, & lvalue,rvalue,gspace2,case,isign,p,a,b,rhoa, & rhob,newgrid,dimnew,countrem,newlr,igrid) enddo elseif(case.eq.3)then pcount=0 lftpos=0 rtpos=gspace1 do i=1,lagsrem yrho=rhoincr(remrho(i)) pos=remrho(i) if(yrho.lt.0.0)then rtpos=gspace1 llimit=0 rlimit=gspace1 llimr=-1.0 rlimr=0.0 llimrho=minfeas rlimrho=0.0 call search(yrho,llimit,rlimit,llimr,rlimr, & llimrho,rlimrho,lftpos,rtpos,case,isign, & rgrid,dimgrid,result,lftbrack,rtbrack, & lvalue,rvalue) else if(pcount.eq.0)then lftpos=gspace1-1 pcount=1 endif rtpos=2.0*lftpos+1 llimit=gspace1-1 rlimit=2*llimit+1 llimr=0.0 rlimr=1.0 llimrho=0.0 rlimrho=1.0 call search(yrho,llimit,rlimit,llimr,rlimr, & llimrho,rlimrho,lftpos,rtpos,case,isign, & rgrid,dimgrid,result,lftbrack,rtbrack, & lvalue,rvalue) endif call calcerr(lvalue,rvalue,yrho,lftbrack,rtbrack,closer, & closerho,error) r(ordrho(pos))=closer estrho(ordrho(pos))=closerho esterr(ordrho(pos))=error abserror=abs(error) erel=errrel(ordrho(pos)) call procerr(abserror,erel,remrho2,lftbrack,rtbrack, & lvalue,rvalue,gspace2,case,isign,p,a,b,rhoa, & rhob,newgrid,dimnew,countrem,newlr,igrid) enddo elseif(case.eq.4)then pcount=0 lftpos=0 rtpos=gspace1 do i=1,lagsrem yrho=rhoincr(remrho(i)) pos=remrho(i) if(yrho.lt.0.0)then lftpos=0 llimit=0 rlimit=gspace1 llimr=-1.0 rlimr=0.0 llimrho=minfeas rlimrho=0.0 case=4 isign=-1 call search(yrho,llimit,rlimit,llimr,rlimr, & llimrho,rlimrho,lftpos,rtpos,case,isign, & rgrid,dimgrid,result,lftbrack,rtbrack, & lvalue,rvalue) else if(pcount.eq.0)then lftpos=0 pcount=1 endif rtpos=gspace1 llimit=0 rlimit=gspace1 llimr=0.0 rlimr=1.0 llimrho=0.0 rlimrho=1.0 case=4 isign=1 call search(yrho,llimit,rlimit,llimr,rlimr, & llimrho,rlimrho,lftpos,rtpos,case,isign, & rgrid,dimgrid,result,lftbrack,rtbrack, & lvalue,rvalue) endif call calcerr(lvalue,rvalue,yrho,lftbrack,rtbrack,closer, & closerho,error) r(ordrho(pos))=closer estrho(ordrho(pos))=closerho esterr(ordrho(pos))=error abserror=abs(error) erel=errrel(ordrho(pos)) call procerr(abserror,erel,remrho2,lftbrack, & rtbrack,lvalue,rvalue,gspace2,case,isign,p,a,b, & rhoa,rhob,newgrid,dimnew,countrem,newlr,igrid) enddo endif else dimnew=0 do i=1,lagsrem yrho=rhoincr(remrho(i)) pos=remrho(i) llimit=igrid(i) rlimit=igrid(i)+gspace2 if((case.ne.4).or.(yrho.ge.0.0))then llimr=a(i) rlimr=b(i) llimrho=rhoa(i) rlimrho=rhob(i) isign=1 else isign=-1 llimr=-b(i) rlimr=-a(i) llimrho=-rhob(i) rlimrho=-rhoa(i) endif lftpos=llimit rtpos=rlimit call search(yrho,llimit,rlimit,llimr,rlimr,llimrho,rlimrho, &lftpos,rtpos,case,isign,rgrid,dimgrid,result,lftbrack,rtbrack, &lvalue,rvalue) call calcerr(lvalue,rvalue,yrho,lftbrack,rtbrack,closer, & closerho,error) r(ordrho(pos))=closer estrho(ordrho(pos))=closerho esterr(ordrho(pos))=error abserror=abs(error) erel=errrel(ordrho(pos)) call procerr(abserror,erel,remrho2,lftbrack,rtbrack, & lvalue,rvalue,gspace2,case,isign,p,a,b,rhoa,rhob, & newgrid,dimnew,countrem,newlr,igrid) enddo endif lagsrem=newlr if(lagsrem.eq.0)then idone=1 else dimgrid=dimnew do i=1,dimnew rgrid(i)=newgrid(i) enddo do i=1,lagsrem remrho(i)=remrho2(i) enddo endif 10 RETURN END Subroutine procerr(abserror,erel,remrho2,lftbrack,rtbrack, &lvalue,rvalue,gspace2,case,isign,p,a,b,rhoa,rhob,newgrid,dimnew, &countrem,newlr,igrid) c Subroutine to determine whether the user-input relative error c has been met for the current lag. If it has, then update the c vector of remaining autocorrelations. If it has not, then c update the grid of AR-process autocorrelations to try during c the next iteration. c Marne C. Cario c ARTA Subroutines c Updates: 10/4/95, 11/8/95, 1/3/96, 1/11/96, 1/16/96 c Subroutines called: None c Calling program: Update integer p,newlr,gspace2,dimnew,n,last,countrem, &case,isign,remrho2(p),igrid(p) real*8 abserror,erel,a(p),b(p),newgrid(50), &lftbrack,rtbrack, &lvalue,rvalue,rhoa(p),rhob(p),temp if(abserror.gt.erel)then n=gspace2-1 last=dimnew if((case.eq.4).and.(isign.eq.-1))then temp=lftbrack lftbrack=-rtbrack rtbrack=-temp endif if((countrem.eq.0).or.(lftbrack.ne.a(countrem)))then do i=1,n xincr=(rtbrack-lftbrack)/gspace2 newgrid(last+i)=lftbrack+i*xincr enddo dimnew=dimnew+n igrid(countrem+1)=last else igrid(countrem+1)=last-n endif countrem=countrem+1 a(countrem)=lftbrack b(countrem)=rtbrack if((case.eq.4).and.(isign.eq.-1))then rhoa(countrem)=-rvalue rhob(countrem)=-lvalue else rhoa(countrem)=lvalue rhob(countrem)=rvalue endif else newlr=newlr-1 do j=countrem+1,newlr remrho2(j)=remrho2(j+1) enddo endif return end Subroutine calcerr(lvalue,rvalue,yrho,lftbrack,rtbrack,closer, &closerho,error) c Subroutine to calculate the relative error of each bracketing value c of the ARTA-process autocorrelation and return the closest value c and the corresponsing AR-process autocorrelation. c Marne C. Cario c ARTA Subroutines c Updates: 10/4/95, 11/8/95, 1/3/96, 1/16/96 c Subroutines called: None c Calling program: update real*8 lvalue,rvalue,yrho,lftbrack,rtbrack,closer,closerho,error, &eleft,eright,aleft eleft=(lvalue-yrho)/dabs(yrho) aleft=-eleft eright=(rvalue-yrho)/dabs(yrho) if(aleft.lt.eright)then closer=lftbrack closerho=lvalue error=eleft else closer=rtbrack closerho=rvalue error=eright endif return end Subroutine search(yrho,llimit,rlimit,llimr,rlimr, &llimrho,rlimrho,lftpos,rtpos,case,isign,rgrid, &dimgrid,result,lftbrack,rtbrack, &lvalue,rvalue) C Subroutine to search for the bracketing values of the AR-process c autocorrelations for the desired ARTA-process autocorrelation. c Uses a binary search procedure. c Marne C. Cario c ARTA Subroutines c Updates: 10/4/95, 11/8/95 c Subroutines called: none c Calling program: update integer llimit,rlimit,case,isign,dimgrid,lftpos, &rtpos,dummy,mid real*8 yrho,llimr,rlimr,llimrho,rlimrho,rgrid(dimgrid), &result(dimgrid),lftbrack,rtbrack,lvalue,rvalue,value dummy=rtpos-1 dowhile(lftpos.lt.dummy) mid=(lftpos+rtpos)/2.0 if((case.ne.4).or.(isign.ne.-1))then value=result(mid) if(yrho.gt.value)then lftpos=mid else rtpos=mid dummy=rtpos-1 endif else value=-1.0*result(mid) if(yrho.gt.value)then rtpos=mid dummy=rtpos-1 else lftpos=mid endif endif enddo c Account for extreme values for the brackets and for symmetric c distributions. if((case.ne.4).or.(isign.ne.-1))then if(lftpos.eq.llimit)then lvalue=llimrho lftbrack=llimr else lvalue=result(lftpos) lftbrack=rgrid(lftpos) endif if(rtpos.eq.rlimit)then rvalue=rlimrho rtbrack=rlimr else rvalue=result(rtpos) rtbrack=rgrid(rtpos) endif else if(lftpos.eq.llimit)then lvalue=-1.0*result(rtpos) rvalue=rlimrho lftbrack=-1.0*rgrid(rtpos) rtbrack=rlimr elseif(rtpos.eq.rlimit)then lvalue=llimrho rvalue=-1.0*result(lftpos) lftbrack=llimr rtbrack=-1.0*rgrid(lftpos) else lvalue=-1.0*result(rtpos) rvalue=-1.0*result(lftpos) lftbrack=-1.0*rgrid(rtpos) rtbrack=-1.0*rgrid(lftpos) endif endif 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 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*8 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 CORREL(nobs,nlags) C Reads nobs observations of a time series from an input file called c "tsdata.dat" and calculates the sample average, sample variance, c and sample autocorrelations through lag nlags. Also calculates an c approximate estimate of the standard deviation of each autocorrelation c estimate. Write the output to the data file "sumstats.out". C Marne C. Cario c ARTA Subroutines c Updates: 9/12/95, 9/16/95, 11/8/95, 11/16/95, 11/20/95 (Tested) C Subroutines called: None C Called by: Input c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c nobs integer number of time-series observations c nlags integer number of autocorrelation lags to estimate c ***************************************************************************** c *********** VARIABLES ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c y double-prec. vector vector containing the nlags previous c observations read in from the file (circular c array) c crosssum double-prec. vector vector containing the sums of the cross- c product terms needed for calculating the c sample autocorrelations c upsum double-prec. vector vector containing the sum of the (k+1)st c to nobs'th observations (needed for calculating c the kth sample autocorrelation) c lowsum double-prec. vector vector containing the sum of the 1st to (n-k)th c observations (needed for calculating the c kth sample autocorrelation) c ysum double precision sum of the observations c yobs double precision current observation - just read in from file c ysumsq double precision sum of squares of the observations c prev integer location of the previous observation in the c circular array y c sampcorr double-prec. vector vector of sample autocorrelations c sampavg double precision sample average c sampvar double precision sample variance c stderr double-prec. vector vector of approximate standard-error estimates c for the elements of sampcorr; based upon c Box and Jenkins pp. 30-35 (Bartlett's approx.) c iflag integer 1 if nlags > nobs/4 c 0 if nlags <= nobs/4 c ***************************************************************************** integer nobs,nlags,prev,m,iflag,n(25),maxlags real*8 sampcorr(25),sampavg,sampvar,y(25),crosssum(25), &upsum(25),lowsum(25),ysum,yobs,ysumsq,stderr(25), &lower,upper,bottom character*1 nonzero(30) maxlags=20 C OPEN INPUT FILE AND OUTPUT FILE open(10,file='tsdata.dat',status='old') C CHECK THAT THE NUMBER OF SAMPLE AUTOCORRELATION LAGS TO BE CALCULATED C IS APPROPRIATE FOR THE SIZE OF THE DATA SET (GUIDELINE GIVEN IN BOX C AND JENKINS pp. 30-35) if(nlags.ge.nobs)then write(8,5) 5 format('Number of autocorrelation lags exceeds sample size') close(8) close(10) stop elseif(nlags.gt.maxlags)then write(8,7)maxlags 7 format('Maximum number of autocorrelation lags is ',i4) close(10) close(8) stop endif iflag=1 div=nobs/4.0 if(nlags.gt.div)then iflag=0 endif C INITIALIZE SUMS AND VARIABLES ysum=0 ysumsq=0 prev=0 do j=1,nlags y(j)=0.0 crosssum(j)=0.0 upsum(j)=0.0 lowsum(j)=0.0 n(j)=nobs-j enddo C READ IN THE OBSERVATIONS AND COMPUTE THE RELEVANT SUMS do i=1,nobs read(10,*)yobs do j=1,nlags C CALCULATE THE LOCATION (m) IN THE CIRCULAR ARRAY OF THE TERM TO MULTIPLY C WITH THE CURRENT OBSERVATION TO GET THE APPROPRIATE CROSS-PRODUCT FOR C LAG j; UPDATE THE CROSS-PRODUCT SUM m=prev-j+1 if(m.le.0)then m=m+nlags endif crosssum(j)=crosssum(j)+yobs*y(m) C UPDATE THE UPPER AND LOWER SUMS if(i.gt.j)then upsum(j)=upsum(j)+yobs endif if(i.le.n(j))then lowsum(j)=lowsum(j)+yobs endif enddo C UPDATE THE SUM AND SUM OF SQUARES OF THE OBSERVATIONS ysum=ysum+yobs ysumsq=ysumsq+yobs*yobs C UPDATE THE POINTER TO THE PREVIOUS OBSERVATION; CURRENT OBSERVATION C BECOMES PREVIOUS ONE if(prev.lt.nlags)then prev=prev+1 else prev=1 endif y(prev)=yobs enddo close(10) C CALCULATE THE OUTPUT QUANTITIES sampavg=ysum/nobs sampvar=ysumsq/nobs-sampavg*sampavg bottom=sampvar*nobs do j=1,nlags sampcorr(j)=(crosssum(j)-sampavg*upsum(j)-sampavg*lowsum(j) + & (nobs-j)*sampavg*sampavg)/bottom enddo stderr(1)=1.0/nobs do j=2,nlags sum=0.0 do i=1,j-1 sum=sum+sampcorr(i)*sampcorr(i) enddo stderr(j)=sqrt((1.0+2.0*sum)/nobs) enddo C DETERMINE WHETHER THE AUTOCORRELATIONS ARE NON-ZERO ACCORDING TO THE C +/- TWO ESTIMATED STANDARD ERRORS CRITERION do j=1,nlags lower=sampcorr(j)-2.0*stderr(j) upper=sampcorr(j)+2.0*stderr(j) if((lower.le.0.0).and.(0.0.le.upper))then nonzero(j)='N' else nonzero(j)='Y' endif enddo C OUTPUT THE SUMMARY STATISTICS TO AN OUTPUT FILE write(8,10) 10 format('SUMMARY STATISTICS FOR TIME-SERIES DATA',//) write(8,20)nobs,sampavg,sampvar 20 format('Observations',3x,'Sample Average',3x,'Sample Variance', &/,'************',3x,'**************',3x,'***************',/, &4x,i6,5x,f14.5,3x,f15.5) write(8,30) 30 format(//,'Lag',3x,'Sample Autocorrelation',3x,'Est. Std. Error', &3x,'Non-Zero?',/,'***',3x,'***********************', &3x,'***************',3x,'********') do j=1,nlags write(8,40)j,sampcorr(j),stderr(j),nonzero(j) 40 format(i3,12x,f5.3,13x,f10.5,9x,a1) enddo write(8,50) 50 format(//,' NOTE: Autocorrelations are considered non-zero if', &/,' the interval defined by the sample average plus or minus two', &/,' estimated standard errors does not contain zero.') if(iflag.eq.0)then write(8,60) 60 format(/,' NOTE: The number of lags for which sample', &/,' autocorrelations were requested is large relative to the', &/,' number of observations. The higher-lag autocorrelations may &be',/,' based upon a small amount of data.') endif return end 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*8 u,empcdf(nempir),fx l=0 r=nempir k=r-1 c do i=1,10 c write(8,*)empcdf(i) c enddo 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 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*8 empcdf(nempir) external bubsort open(9,file='empir.dat',status='old') do i=1,nempir read(9,*)empcdf(i) enddo call bubsort(empcdf,nempir) c do i=1,10 c write(8,*)empcdf(i) c enddo close(9) RETURN END Subroutine normal(p,rho,dparam,errrel) c Subroutine to calculate the parameters of an AR(p) process c with mean=dparam(1) and variance=dparam(2). c Marne C. Cario c ARTA Subroutines c Updates: 11/15/95, 11/16/95, 11/20/95 (Tested), 1/3/96, 1/10/96, c 2/19/96 c Subroutines called: YW, Station c Calling program: input integer p,istat real*8 rho(p),errrel(p),alpha(5),dparam(5),mean,variance, &sigma2,sum external yw,station mean=dparam(1) variance=dparam(2) write(8,5) 5 format(' ********** USER INPUT **********',//) 10 write(8,15) 15 format(' ARTA Distribution: Normal') write(8,16)dparam(1),dparam(2) 16 format(' Mean: ',f10.5,/,' Variance: ',f10.5,/) 150 write(8,155)p 155 format(' Number of ARTA autocorrelation lags to match: ',i2,/) write(8,160) 160 format(' Lag',3x,'Desired Autocorrelation', &3x,'Desired Relative Error', &/,' ***',3x,'***********************',3x,'**********************') do i=1,p write(8,170)i,rho(i),errrel(i) 170 format(1x,i2,7x,f15.13,10x,f7.5) enddo write(8,180) 180 format(/,' **********************************') c Solve for the AR parameters using the Yule Walker equations. c Check stationarity of the AR process. call YW(p,rho,alpha) call Station(alpha,p,istat) write(8,200) 200 format(/,' NOTE: The normal distribution is modeled as an', &/,8x,'AR(p) process with mean mu and error variance sigma^2.', &/,8x,'Relative errors are essentially zero. The parameters of', &/,8x,'the AR(p) process follow.') c CALCULATE THE VARIANCE OF THE ERROR TERM sum=0.0 do i=1,p sum=sum+alpha(i)*rho(i) enddo sigma2=variance*(1.0-sum) write(8,210)mean,sigma2 210 format(//,' Mu = ',f8.4,5x,'Sigma^2 = ',f8.4) write(8,220) 220 format(//,' i ',10x,' Alpha(i)') do i=1,p write(8,230)i,alpha(i) 230 format(1x,i2,10x,f10.5) enddo if(istat.eq.1)then write(8,240) 240 format(//,' ERROR: The AR output process is not stationary.') stop endif C WRITE INPUT FILE FOR ARTAGEN open(15,file='input.gen',status='unknown') write(15,500) 500 format('obs.ag'/,'stats.ag') ii=1 write(15,510)ii 510 format(i1) write(15,520)p 520 format(i2) ii=1000 write(15,530)ii 530 format(i4) do i=1,p write(15,540)rho(i) 540 format(f9.5) enddo c NORMAL DISTRIBUTION c dparam(1) = MEAN c dparam(2) = VARIANCE 610 write(15,615)dparam(1),dparam(2) 615 format(f10.5,/,f10.5) close(15) RETURN END subroutine uniform(p,rho,dparam,errrel) c Solves for the AR-process a.c.'s using the closed-form solution c given in Li and Hammond, "Generation of Pseudorandom Numbers with c Specified Univariate Distributions and Correlation Coefficients", c IEEE Trans. on Systems, Man, and Cybernetics, September, 1975, c pp. 557-561. c Marne C. Cario c ARTA Subroutines c Updates: 11/15/95, 11/20/95 (Tested), 1/3/96, 1/10/96, c 2/19/96 c Subroutines called: YW, station c Calling program: input integer p real*8 rho(p),r(5),alpha(5),piovr6,dparam(5),errrel(p) external yw,station piovr6=acos(-1.0)/6.0 c Echo the user input to the output file. write(8,5) 5 format('********** USER INPUT **********',/) 30 write(8,35) 35 format('ARTA Distribution: Continuous Uniform') write(8,36)dparam(1),dparam(2) 36 format('Left Endpoint: ',f10.5,/,'Right Endpoint: ',f10.5,/) 150 write(8,155)p 155 format('Number of ARTA autocorrelation lags to match: ',i2,/) write(8,160) 160 format('Lag',3x,'Desired Autocorrelation', &3x,'Desired Relative Error', &/,' ***',3x,'***********************',3x,'**********************') do i=1,p write(8,170)i,rho(i),errrel(i) 170 format(1x,i2,7x,f15.13,10x,f7.5) enddo write(8,180) 180 format(/,'**********************************') do i=1,p r(i)=2.0*sin(piovr6*rho(i)) enddo c Solve for the AR parameters using the Yule Walker equations. c Check stationarity of the AR process. call YW(p,r,alpha) call Station(alpha,p,istat) write(8,200) 200 format(///, &'A closed-form solution is used for the uniform distribution.' &,/, &'Relative errors are essentially zero.') write(8,210) 210 format(//,'Lag',10x,'ARTA Correl.',10x,'AR Correl.') do i=1,p write(8,220)i,rho(i),r(i) 220 format(2x,i2,9x,f8.6,14x,f8.6) enddo write(8,230) 230 format(//,' i ',10x,' Alpha(i)') do i=1,p write(8,240)i,alpha(i) 240 format(1x,i2,10x,f10.5) enddo if(istat.eq.1)then write(8,250) 250 format(//,'ERROR: The AR output process is not stationary.') close(8) stop endif C WRITE INPUT FILE FOR ARTAGEN open(15,file='input.gen',status='unknown') write(15,500) 500 format('obs.ag'/,'stats.ag') ii=3 write(15,510)ii 510 format(i1) write(15,520)p 520 format(i2) ii=1000 write(15,530)ii 530 format(i4) do i=1,p write(15,540)r(i) 540 format(f9.5) enddo do i=1,p write(15,540)rho(i) enddo c CONTINUOUS UNIFORM DISTRIBUTION c dparam(1) = LEFT ENDPOINT c dparam(2) = RIGHT ENDPOINT 630 write(15,635)dparam(1),dparam(2) 635 format(f10.5,/,f10.5) close(15) return end SUBROUTINE OUTER(rgrid,dimgrid,distrib,tolin,tolout, &dparam,mean,variance,p,maxsubs,minfeas,result) c Calculates an approximation to the value of the double integral c for each element of rgrid c Marne C. Cario c ARTA subroutine c Updates: 8/31/95, 9/6/95, 10/10/95, 11/8/95, 11/16/95, 11/28/95, c 2/9/96 c Release 1.1 5/30/96 c SUBROUTINES CALLED: gkpoints, inner, putheap, remheap, output c CALLED BY: Main program c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c rgrid double-prec. vector vector of AR-correlations for which to c approximate the outer integral c dimgrid integer number of AR-correlations in grid c distrib integer code corresponding to the marginal distribution c tolin double precision absolute-error tolerance for the inner integral c tolout double precision absolute-error tolerance for the outer integral c dparam double-prec. vector vector of ARTA-distribution parameters c mean double precision mean of the ARTA distribution c variance double precision variance of the ARTA distribution c p integer number of a.c. lags to match c maxsubs integer maximum number of subintervals allowed c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c result double-prec. vector dimgrid x 1 vector of the integral c approximations for each value in rgrid c ***************************************************************************** c *********** VARIABLES ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c err double-prec. vector dimgrid x 1 vector of the error approximations c corresponding to the members of result c outera double-prec. vector vector of the left endpoints of the integration c intervals considered so far c outerb double-prec. vector vector of the right endpoints of the integration c intervals considered so far c error double-prec. array error estimates for each interval/correlation c combination (rows=interval, cols=correlation) c considered so far c eheap double-prec. vector heap array of maximum errors (over all c correlations) for each subinterval that has c not yet been bisected c iheap integer array array of subinterval numbers corresponding to c the members of eheap c last integer pointer to the last element of eheap c ierr integer interval number corresponding to the maximum c error c outsubs integer number of subintervals considered so far c (for the outer integral) c itol integer 0 if the outer integral has not yet c been calculated to sufficient accuracy for c at least one correlation value c 1 if the accuracy requirement is met for c each correlation vlue c l double-prec. vector 2x1 vector of the left endpoints for the c current bisected subinterval; for the first c iteration there is only 1 left endpoint (-1.0) c rt double-prec. vector 2x1 vector of the right endpoints for the c current bisected subinterval; for the first c iteration there is only 1 rt. endpoint (1.0) c gkpts double-prec. vector 15x1 vector of the Gauss-Kronrod abscissae c for the current subinterval c gsum double-prec. array Gauss sums for each interval/correlation c combination (row=interval, col=correlation) c considered so far c ksum double-prec. array Kronrod sums for each interval/correlation c combination (row=interval, col=correlation) c considered so far c valin double-prec. vector dimgrid x 1 vector of the values of the c inner integral at each correlation value c for the current value of the outer c integration variable c errmax double-precision maximum absolute error over all correlation c values for the current interval c width double-precision width of the current interval c ***************************************************************************** implicit real*8 (x) integer dimgrid,distrib,itol,outsubs,maxsubs,k, &iheap(50),last,ierr,p, &nempir, elimit real*8 rgrid(dimgrid),tolin,tolout,result(dimgrid),err(50), &dparam(5),mean,variance,outera(50),outerb(50), &l(2),rt(2),gkpts(15),error(50,50),valin(50), &gsum(50,50),ksum(50,50),errmax,width,eheap(50),rangeo, &minfeas, &rho(5),empcdf(300),prob(300),r(5),estrho(5),esterr(5),alpha(5) common/gkwts/xgwts(15),xkwts(15) external gkpoints,inner,putheap,remheap,output C INITIALIZE THE NUMBER OF SUBINTERVALS AND THE LAST ELEMENT OF THE HEAP outsubs=0 last=0 C INITIALIZE THE INTEGRAL AND ERROR VALUES TO ZERO do j=1,dimgrid result(j)=0.0 err(j)=0.0 enddo C CALCULATE THE VALUE OF THE OUTER INTEGRALS TO WITHIN THE GIVEN TOLERANCE itol=0 dowhile(itol.eq.0) C IF THIS IS THE FIRST TIME THROUGH, THEN CALCULATE THE INTEGRAL OVER C THE ENTIRE INTERVAL (-10,10); OTHERWISE, BISECT THE INTERVAL WITH THE C MAXIMUM ERROR AND SUBTRACT OFF THE VALUES FOR THE INTERVAL THAT WILL BE C BISECTED FROM THE RESULT AND THE ERROR (FOR EACH CORRELATION VALUE) if(outsubs.eq.0)then l(1)=-10.0 rt(1)=10.0 k=1 else width=(outerb(ierr)-outera(ierr))/2.0 l(1)=outera(ierr) rt(1)=outera(ierr)+width l(2)=outera(ierr)+width rt(2)=outerb(ierr) k=2 do j=1,dimgrid result(j)=result(j)-ksum(ierr,j) err(j)=err(j)-error(ierr,j) enddo endif do m=1,k C INCREMENT THE NUMBER OF INTERVALS CONSIDERED SO FAR AND CALCULATE C THE G-K POINTS FOR THE SUBINTERVAL BEING CONSIDERED outsubs=outsubs+1 if(outsubs.gt.maxsubs)then call output(7,laglimit,maxsubs,elimit,p,r,estrho, &esterr,alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif c write(8,*)outsubs outera(outsubs)=l(m) outerb(outsubs)=rt(m) call gkpoints(l(m),rt(m),gkpts) rangeo=rt(m)-l(m) C INITIALIZE THE QUADRATURE SUMS do j=1,dimgrid gsum(outsubs,j)=0.0 ksum(outsubs,j)=0.0 enddo C FOR EACH G-K POINT, CALCULATE THE VALUE OF THE INNER INTEGRAL C (AT EACH CORRELATION VALUE IN rgrid) do i=1,15 call inner(gkpts(i),rangeo,rgrid,dimgrid,tolin, &distrib,dparam,p,maxsubs,minfeas,valin) C UPDATE THE QUADRATURE SUM FOR THE OUTER INTEGRAL AT EACH CORRELATION C VALUE IN rgrid do j=1,dimgrid gsum(outsubs,j)=gsum(outsubs,j)+xgwts(i)*valin(j) ksum(outsubs,j)=ksum(outsubs,j)+xkwts(i)*valin(j) enddo enddo C UPDATE THE INTEGRAL AND ERROR APPROXIMATIONS FOR EACH CORRELATION VALUE C APPROXIMATE THE QAUDRATURE ERRORS BY THE DIFFERENCE BETWEEN THE C GAUSS AND KRONROD SUMS; KEEP TRACK OF THE MAXIMUM ERROR FOR THE C CURRENT INTERVAL errmax=0.0 do j=1,dimgrid result(j)=result(j)+ksum(outsubs,j) error(outsubs,j)=abs(ksum(outsubs,j)-gsum(outsubs,j)) err(j)=err(j)+error(outsubs,j) if(error(outsubs,j).gt.errmax)then errmax=error(outsubs,j) endif enddo C PUT THE ERROR (ABSOLUTE VALUE) ON THE BINARY HEAP OF C ERROR VALUES TO BE CONSIDERED IN THE FUTURE call putheap(maxsubs,errmax,outsubs,eheap,iheap,last) enddo C CHECK WHETHER ALL THE ERRORS ARE WITHIN THE TOLERANCE C IF THEY ARE, THEN ITOL=1 AND THE DO-WHILE LOOP TERMINATES C IF AT LEAST ONE IS NOT, THEN REMOVE THE LARGEST ERROR FROM THE C HEAP AND CONTINUE. do j=1,dimgrid itol=1 if(err(j).gt.tolout)then itol=0 call remheap(maxsubs,eheap,iheap,last,ierr) goto 10 endif enddo 10 enddo c CHANGE THE RESULT VECTOR FROM EXPECTED VALUES TO CORRELATIONS c write(8,*)mean,variance do j=1,dimgrid result(j)=(result(j)-mean*mean)/variance enddo RETURN 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*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 INNER(z1,rangeo,rgrid,dimgrid,tolin,distrib, &dparam,p,maxsubs,minfeas,valin) c Calculates an approximation to the inner integral at a fixed value c of the outer integration variable, for each AR-correlation value c in a grid c Marne C. Cario c ARTA subroutine c Updates: 8/31/95, 9/6/95, 10/10/95, 11/8/95, 11/15/95, 11/16/95, c 11/28/95, 1/3/96, 1/7/96, 1/10/96, 1/15/96, 1/16/96, c 2/9/96 c Release 1.1 6/4/96 c SUBROUTINES CALLED: gkpoints, putheap, remheap, fnorm2, bvnpdf, output, c vstud, expinv, gaminv, weibinv, juinv, jbinv c CALLED BY: outer c *********** INPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c rgrid double-prec. vector vector of AR-correlations for which to c approximate the outer integral c dimgrid integer number of AR-correlations in grid c tolin double precision absolute-error tolerance for the inner integral c distrib integer code corresponding to the marginal distribution c dparam double-prec. vector vector of ARTA-distribution parameters c (if applicable) c p integer number of a.c. lags to match c maxsubs integer maximum number of subintervals allowed c ***************************************************************************** c *********** OUTPUT ********************************************************* c VARIABLE TYPE DESCRIPTION c ******** **** *********** c valin double-prec. vector dimgrid x 1 vector of the inner integral c approximations for each value in rgrid c ***************************************************************************** c*********** VARIABLES **************************************************** c VARIABLE TYPE DESCRIPTION c ******** **** *********** c err double-prec. vector dimgrid x 1 vector of the error approximations c corresponding to the members of valin c innera double-prec. vector vector of the left endpoints of the integration c intervals considered so far c innerb double-prec. vector vector of the right endpoints of the integration c intervals considered so far c error double-prec. array error estimates for each interval/correlation c combination (rows=interval, cols=correlation) c considered so far c eheap double-prec. vector heap array of maximum errors (over all c correlations) for each subinterval that has c not yet been bisected c iheap integer array array of subinterval numbers corresponding to c the members of eheap c last integer pointer to the last element of eheap c ierr integer interval number corresponding to the maximum c error c insubs integer number of subintervals considered so far c (for the inner integral) c itol integer 0 if the inner integral has not yet c been calculated to sufficient accuracy for c at least one correlation value c 1 if the accuracy requirement is met for c each correlation vlue c l double-prec. vector 2x1 vector of the left endpoints for the c current bisected subinterval; for the first c iteration there is only 1 left endpoint (-1.0) c r double-prec. vector 2x1 vector of the right endpoints for the c current bisected subinterval; for the first c iteration there is only 1 rt. endpoint (1.0) c gkpts double-prec. vector 15x1 vector of the Gauss-Kronrod abscissae c for the current subinterval c gsum double-prec. array Gauss sums for each interval/correlation c combination (row=interval, col=correlation) c considered so far c ksum double-prec. array Kronrod sums for each interval/correlation c combination (row=interval, col=correlation) c considered so far c errmax double-precision maximum absolute error over all correlation c values for the current interval c width double-precision width of the current interval c ***************************************************************************** implicit real*8 (x) integer dimgrid,distrib,ierr,insubs,k, &last,iheap(25),itol,maxsubs,p,ifault,df, &nempir,elimit real*8 rgrid(dimgrid),tolin,valin(dimgrid),err(25), &innera(25),innerb(25),error(25,25),l(2),rt(2),gkpts(15), &gsum(25,25),ksum(25,25),errmax,width, &z1,z2,u1,u2,term1,term2,term3,term,dparam(5), &eheap(25),w(4),bvnpdf,rangeo,rangei,q,pdf,g, &vstud,expinv,weibinv,lninv,juinv,jbinv,gaminv, &minfeas,mean,variance, &rho(5),empcdf(300),prob(300),r(5),estrho(5),esterr(5),alpha(5) common/gkwts/xgwts(15),xkwts(15) external gkpoints,putheap,remheap,fnorm2,bvnpdf,output, &vstud,expinv,weibinv,lninv,juinv,jbinv,gaminv c open(15,file='t.out',status='unknown') C INITIALIZE THE NUMBER OF SUBINTERVALS AND THE LAST ELEMENT OF THE HEAP c write(15,2)distrib c 2 format('distrib=',i2) insubs=0 last=0 C CALCULATE THE QUANTITIES THAT DEPEND UPON THE OUTER VARIABLE ONLY call fnorm2(z1,u1,q,pdf) c write(15,*)dimgrid c write(15,5)z1,u1 c 5 format(' z1= ',f10.5,5x,' u1 = ',f10.5) goto(10,20,30,40,50,60,70,80,90),distrib C NORMAL DISTRIBUTION (Special case - should never get here) 10 print*,'Error - do not need to integrate Normal distribution.' stop C STUDENT'S t DISTRIBUTION 20 df=idint(dparam(1)) term1=vstud(u1,df,w,ifault) c write(15,21)df,term1 c 21 format('df = ',i3,5x,'term1 = ',f10.5) goto 200 C UNIFORM DISTRIBUTION 30 print*,'Error - do not need to integrate Uniform distribution.' stop C EXPONENTIAL DISTRIBUTION 40 term1=expinv(dparam,u1) goto 200 C GAMMA DISTRIBUTION 50 g=gaminv(dparam,u1,ifault) term1=dparam(1)*g goto 200 C WEIBULL DISTRIBUTION 60 term1=weibinv(dparam,u1) goto 200 C LOGNORMAL 70 term1=lninv(dparam,u1) goto 200 C JOHNSON UNBOUNDED FAMILY 80 term1=juinv(dparam,z1) goto 200 C JOHNSON BOUNDED FAMILY 90 term1=jbinv(dparam,z1) goto 200 C INITIALIZE THE INTEGRAL AND ERROR VALUES TO ZERO 200 do j=1,dimgrid valin(j)=0.0 err(j)=0.0 enddo C CALCULATE THE VALUE OF THE INNER INTEGRALS TO WITHIN THE GIVEN TOLERANCE itol=0 dowhile(itol.eq.0) C IF THIS IS THE FIRST TIME THROUGH, THEN CALCULATE THE INTEGRAL OVER C THE ENTIRE INTERVAL (-1,1); OTHERWISE, BISECT THE INTERVAL WITH THE C MAXIMUM ERROR AND SUBTRACT OFF THE VALUES FOR THE INTERVAL THAT WILL BE C BISECTED FROM THE RESULT AND THE ERROR (FOR EACH CORRELATION VALUE) c write(15,201)insubs,ierr c 201 format('insubs=',i3,5x,'ierr=',i3) if(insubs.eq.0)then l(1)=-10.0 rt(1)=10.0 k=1 c write(15,202)l(1),rt(1) c 202 format('l(1) = ',f10.5,5x,'rt(1) = ',f10.5) else width=(innerb(ierr)-innera(ierr))/2.0 l(1)=innera(ierr) rt(1)=innera(ierr)+width l(2)=innera(ierr)+width rt(2)=innerb(ierr) c write(15,203)l(1),rt(1) c 203 format('l(1) = ',f10.5,5x,'rt(1) = ',f10.5) c write(15,202)l(2),rt(2) c 204 format('l(2) = ',f10.5,5x,'rt(2) = ',f10.5) k=2 do j=1,dimgrid valin(j)=valin(j)-ksum(ierr,j) err(j)=err(j)-error(ierr,j) enddo endif do m=1,k C INCREMENT THE NUMBER OF INTERVALS CONSIDERED SO FAR AND CALCULATE C THE G-K POINTS FOR THE SUBINTERVAL BEING CONSIDERED insubs=insubs+1 if(insubs.gt.maxsubs)then call output(8,laglimit,maxsubs,elimit,p,r,estrho, &esterr,alpha,minfeas,distrib,rho,dparam,mean,variance,nempir, &empcdf,prob) endif innera(insubs)=l(m) innerb(insubs)=rt(m) call gkpoints(l(m),rt(m),gkpts) rangei=rt(m)-l(m) c do i=1,15 c write(15,205)i,gkpts(i) c 205 format(' i = ',i2,5x,'gkpts(i) = ',f10.5) c enddo C INITIALIZE THE QUADRATURE SUMS do j=1,dimgrid gsum(insubs,j)=0.0 ksum(insubs,j)=0.0 enddo C FOR EACH G-K POINT, CALCULATE THE VALUE OF THE INNER INTEGRAL C (AT EACH CORRELATION VALUE IN rgrid) do i=1,15 C CALCULATE THE QUANTITIES THAT ARE THE SAME FOR ALL CORRELATIONS z2=gkpts(i) call fnorm2(z2,u2,q,pdf) c write(15,206)z2,u2 c 206 format(' z2= ',f10.5,5x,' u2 = ',f10.5) goto(210,220,230,240,250,260,270,280,290),distrib C NORMAL DISTRIBUTION (Special case - should never get here) 210 print*,'Error - do not need to integrate Normal distribution.' stop C STUDENT'S t DISTRIBUTION 220 df=idint(dparam(1)) term2=vstud(u2,df,w,ifault) c write(15,222)df,term2 c 222 format('df = ',i3,5x,'term2 = ',f10.4) goto 400 C UNIFORM DISTRIBUTION 230 print*,'Error - do not need to integrate Uniform distribution.' stop C EXPONENTIAL DISTRIBUTION 240 term2=expinv(dparam,u2) goto 400 C GAMMA DISTRIBUTION 250 term2=dparam(1)*gaminv(dparam,u2,ifault) goto 400 C WEIBULL DISTRIBUTION 260 term2=weibinv(dparam,u2) goto 400 C LOGNORMAL DISTRIBUTION 270 term2=lninv(dparam,u2) goto 400 C JOHNSON UNBOUNDED FAMILY 280 term2=juinv(dparam,z2) goto 400 C JOHNSON BOUNDED FAMILY 290 term2=jbinv(dparam,z2) goto 400 400 term=(rangeo/2.0)*(rangei/2.0)*term1*term2 C UPDATE THE QUADRATURE SUM FOR THE INNER INTEGRAL AT EACH CORRELATION C VALUE IN rgrid do j=1,dimgrid term3=bvnpdf(z1,z2,rgrid(j)) c write(15,402)term1,term2,term3,term4 c 402 format('terms: ',f10.5,3x,f10.5,3x,f10.5,3x,f10.5) gsum(insubs,j)=gsum(insubs,j)+xgwts(i)*term3*term ksum(insubs,j)=ksum(insubs,j)+xkwts(i)*term3*term c write(15,405)j,gsum(insubs,j),ksum(insubs,j) c 405 format('j=',i2,3x,'gsum(j)=',f10.5,3x,'ksum(j)=',f10.5) enddo enddo C UPDATE THE INTEGRAL AND ERROR APPROXIMATIONS FOR EACH CORRELATION VALUE C APPROXIMATE THE QAUDRATURE ERRORS BY THE DIFFERENCE BETWEEN THE C GAUSS AND KRONROD SUMS; KEEP TRACK OF THE MAXIMUM ERROR FOR THE C CURRENT INTERVAL errmax=0.0 do j=1,dimgrid valin(j)=valin(j)+ksum(insubs,j) error(insubs,j)=abs(ksum(insubs,j)-gsum(insubs,j)) err(j)=err(j)+error(insubs,j) if(error(insubs,j).gt.errmax)then errmax=error(insubs,j) endif enddo C PUT THE ERROR (ABSOLUTE VALUE) ON THE BINARY HEAP OF C ERROR VALUES TO BE CONSIDERED IN THE FUTURE c write(15,898)insubs c 898 format(' insubs = ',i2) call putheap(maxsubs,errmax,insubs,eheap,iheap,last) c write(15,899)errmax c 899 format('errmax = ',f10.5) c write(15,900)last c 900 format('last = ',i2) enddo C CHECK WHETHER ALL THE ERRORS ARE WITHIN THE TOLERANCE C IF THEY ARE, THEN ITOL=1 AND THE DO-WHILE LOOP TERMINATES C IF AT LEAST ONE IS NOT, THEN REMOVE THE LARGEST ERROR FROM THE C HEAP AND CONTINUE. itol=1 do j=1,dimgrid if(abs(err(j)).gt.tolin)then call remheap(maxsubs,eheap,iheap,last,ierr) c write(15,905)last,ierr c 905 format('last = ',i2,5x,'ierr=',i2) itol=0 goto 500 endif enddo 500 enddo c close(15) RETURN END subroutine innerd(zlow,zhigh,z,empcdf,rgrid,dimgrid, &nempir,valin) c Marne C. Cario c ARTA Subroutines c Updates: 2/10/96 c Calculates the inner integral at a fixed value of the outer c integral when the ARTA distribution is discrete. c Called by: outerd c Subroutines/Functions called: bvnarea integer dimgrid,nempir,nn,inf(7),ifault,n real*8 zlow,zhigh,z(nempir),empcdf(nempir),rgrid(dimgrid), &valin(dimgrid),y2,zlow2,zhigh2,a(7),b(7),sig(15),eps,prob, &bound c Initialize the sums. do j=1,dimgrid valin(j)=0.0 enddo inf(1)=2 inf(2)=2 n=2 eps=dble(0.0001) c Calculate the first portion of the inner integral. zlow2=dble(-10.0) zhigh2=z(1) y2=empcdf(1) do j=1,dimgrid a(1)=zhigh a(2)=zhigh2 b(1)=zlow b(2)=zlow2 sig(1)=rgrid(j) call mulnor(a,b,sig,eps,n,inf,prob,bound,ifault) valin(j)=valin(j)+prob*y2 enddo nn=nempir-1 c Calculate the next (nempir-1) portions of the inner integral. do i=2,nn zlow2=z(i-1) zhigh2=z(i) y2=empcdf(i) do j=1,dimgrid a(1)=zhigh a(2)=zhigh2 b(1)=zlow b(2)=zlow2 sig(1)=rgrid(j) call mulnor(a,b,sig,eps,n,inf,prob,bound,ifault) valin(j)=valin(j)+prob*y2 enddo enddo c Calculate the final portion of the inner integral y2=empcdf(nempir) zlow2=z(nn) zhigh2=dble(10.0) do j=1,dimgrid a(1)=zhigh a(2)=zhigh2 b(1)=zlow b(2)=zlow2 sig(1)=rgrid(j) call mulnor(a,b,sig,eps,n,inf,prob,bound,ifault) valin(j)=valin(j)+prob*y2 enddo return end subroutine outerd(empcdf,z,nempir,rgrid,dimgrid, &mean,variance,result) c Marne C. Cario c ARTA Subroutines c Updates: 2/10/96 c Calculates the correlation for the members of rgrid when c the ARTA distribution is a discrete distribution. integer nempir,dimgrid,n real*8 empcdf(nempir),z(nempir),rgrid(dimgrid), &result(dimgrid),zhigh,valin(25),y1,zlow,mean, &variance c Initialize the sums do j=1,dimgrid result(j)=0.0 enddo c Integrate the first portion of the outer integral zlow=dble(-10.0) zhigh=z(1) y1=empcdf(1) call innerd(zlow,zhigh,z,empcdf,rgrid,dimgrid, &nempir,valin) do j=1,dimgrid result(j)=result(j)+valin(j)*y1 enddo c Integrate the next (nempir-2) portions of the outer integral n=nempir-1 do i=2,n y1=empcdf(i) zlow=z(i-1) zhigh=z(i) call innerd(zlow,zhigh,z,empcdf,rgrid,dimgrid, &nempir,valin) do j=1,dimgrid result(j)=result(j)+valin(j)*y1 enddo enddo c Integrate the final portion of the outer integral y1=empcdf(nempir) zlow=z(n) zhigh=dble(10.0) call innerd(zlow,zhigh,z,empcdf,rgrid,dimgrid, &nempir,valin) do j=1,dimgrid result(j)=result(j)+valin(j)*y1 enddo c Change from expected values to correlations do j=1,dimgrid result(j)=(result(j)-mean*mean)/variance enddo return end real*8 function gaminv(dparam,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,dparam(5),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/ alpha=dparam(2) 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 real*8 function expinv(dparam,u1) c c function to generate negative exponential variate c via inverse cdf c real*8 mean,u1,dparam(5) if(u1.gt.0.999999)then u1=0.999999 endif mean=1.0/dparam(1) expinv = (-dlog(1.0-u1))*mean return end real*8 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*8 dparam(5), alpha, beta, u, w, u2 beta=dparam(1) alpha=dparam(2) if(u.gt.0.999999)then u=0.999999 endif u2=1.0-u w = -1.0*dlog(u2) w=w**(1.0/alpha) weibinv = w*beta return end real*8 function lninv(dparam, u) c Computes the inverse of the lognormal distribution. c Marne C. Cario, 11/9/95 real*8 dparam(5),u,norm,mu,sigma2,vnorm2 mu=dparam(1) sigma2=dparam(2) norm=vnorm2(u,ifault) norm=mu+dsqrt(sigma2)*norm lninv=dexp(norm) return end real*8 function juinv(dparam,z) c Function to invert a Johnson Unbounded family distribution. real*8 dparam(5),z,s,term1,term2,term c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter c z = standard normal r.v. c write(8,*)dparam(1),dparam(2),dparam(3),dparam(4) term=(z-dparam(1))/dparam(2) term1=dexp(term) term2=dexp(dble(-1.0)*term) s=dble(0.5)*(term1-term2) c write(8,*)term,s juinv=dparam(3)+dparam(4)*s return end real*8 function jbinv(dparam,z) c Function to invert a Johnson Bounded family distribution. real*8 dparam(5),z,term,num c dparam(1)=first shape parameter c dparam(2)=second shape parameter c dparam(3)=location parameter c dparam(4)=scale parameter term=(z-dparam(1))/dparam(2) num=dble(1.0)+dexp(dble(-1.0)*term) num=dble(1.0)/num jbinv=dparam(3)+dparam(4)*num 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 df real*8 w(4),phi,phin,phix2,p2tail,y,t,x,y2,c,vnorm2,arg,dwarf, &xpiovr2,sinarg 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 Date 30 Sept 1986 c Check for input error c write(15,1)df,phi c 1 format('df in vstud =',i3,'phi=',f10.8) 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=dble(0.000001) goto 205 else phi=dble(0.999999) c write(15,3)phi c 3 format('phi after modific = ',f18.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 write(15,206)p2tail c 206 format('p2tail = ',f18.15) 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 c write(15,*)'got there 1' w(2)=1.0/(dble(df)-0.5) c write(15,*)'got there 2' 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)* &dble(df) 250 x = w(4)*p2tail c=w(3) c write(15,251)x,w(4) c 251 format('x=',f18.15,5x,'w4=',f15.8) c write(15,*)'got there 3' y=x**(2./dble(df)) c write(15,*)'got there 4' 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*(dble(df)-4.5)*(x+.6) c=(((.05*w(4)*x-5.)*x-7.)*x-2.)*x+w(1)+c c write(15,*)'got there 5' 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=(dble(df)+6.)/(dble(df)+y) - 0.089*w(4) - .822 c write(15,*)'got there 6' y2=y2*(dble(df)+2.)*3. y2=(1./y2+.5/(dble(df)+4.))*y -1. c write(15,*)'got there 7' y2=y2*(dble(df)+1.)/(dble(df)+2.) c write(15,*)'got there 8' c write(15,*)y y=y2+1./y c write(15,*)'got there 9' 280 t=dsqrt(dble(df)*y) c write(15,*)'got there 10' 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 subroutine rhomin(distrib,dparam,minfeas) implicit real*8 (x) integer maxsubs parameter(maxsubs=200) integer distrib,ierr,insubs,k,last,iheap(maxsubs),itol, &m,d,ifault,i real*8 tol,minfeas,err,a(maxsubs),b(maxsubs), &error(maxsubs),l(2), &r(2),gkpts(15),gsum(maxsubs),ksum(maxsubs),width,u1,u2, &dparam(5),eheap(maxsubs),range,expinv,weibinv, &lninv,juinv,jbinv,term1,term2,term,gaminv,vnorm2, &z1,z2 common/gkwts/xgwts(15),xkwts(15) external gkpoints,putheap,remheap,gaminv,expinv,weibinv, &lninv,juinv,jbinv,vnorm2 c write(8,1)dparam(1),dparam(2),dparam(3),dparam(4),distrib c 1 format('Subroutine Rhomin',/,f10.5,/,f10.5,/,f10.5,/,f10.5, c &/,i2) d=distrib-4 tol=0.001 insubs=0 last=0 minfeas=0.0 err=0.0 itol=0 dowhile(itol.eq.0) if(insubs.eq.0)then l(1)=dble(0.0) r(1)=dble(1.0) k=1 else width=(b(ierr)-a(ierr))/dble(2.0) l(1)=a(ierr) r(1)=a(ierr)+width l(2)=a(ierr)+width r(2)=b(ierr) k=2 minfeas=minfeas-ksum(ierr) err=err-error(ierr) endif c write(8,7)minfeas c 7 format('minfeas:',f10.5) do m=1,k insubs=insubs+1 if(insubs.gt.maxsubs)then write(8,10) 10 format('Maximum subintervals exceeded & while finding rhomin') close(8) stop endif a(insubs)=l(m) b(insubs)=r(m) call gkpoints(l(m),r(m),gkpts) range=r(m)-l(m) gsum(insubs)=0.0 ksum(insubs)=0.0 do i=1,15 u1=gkpts(i) u2=dble(1.0)-u1 goto(110,120,130,140,150),d C GAMMA DISTRIBUTION 110 term1=dparam(1)*gaminv(dparam,u1,ifault) term2=dparam(1)*gaminv(dparam,u2,ifault) goto 200 C WEIBULL DISTRIBUTION 120 term1=weibinv(dparam,u1) term2=weibinv(dparam,u2) goto 200 C LOGNORMAL DISTRIBUTION 130 term1=lninv(dparam,u1) term2=lninv(dparam,u2) goto 200 C JOHNSON UNBOUNDED FAMILY 140 z1=vnorm2(u1,ifault) z2=vnorm2(u2,ifault) term1=juinv(dparam,z1) term2=juinv(dparam,z2) c write(8,142)u1,z1,term1,u2,z2,term2 c 142 format(f10.5,3x,f10.5,3x,f10.5,/,f10.5,3x, c &f10.5,3x,f10.5) goto 200 C JOHNSON BOUNDED FAMILY 150 term1=jbinv(dparam,z1) term2=jbinv(dparam,z2) 200 term=(range/dble(2.0))*term1*term2 gsum(insubs)=gsum(insubs)+xgwts(i)*term ksum(insubs)=ksum(insubs)+xkwts(i)*term enddo minfeas=minfeas+ksum(insubs) error(insubs)=dabs(ksum(insubs)-gsum(insubs)) err=err+error(insubs) c write(8,205)insubs,ksum(insubs),minfeas,error(insubs) c 205 format('insubs:',i4,/,'ksum:',f10.5,/,'minfeas:', c &f10.5,/,'error: ',f10.5) call putheap(maxsubs,error(insubs),insubs,eheap, &iheap,last) c write(8,*)last,eheap,iheap enddo if(dabs(err).gt.tol)then call remheap(maxsubs,eheap,iheap,last,ierr) c write(8,215)ierr c 215 format('ierr:',i4) c write(8,*)last,eheap,iheap else itol=1 endif enddo return end subroutine erhomin(empcdf,nempir,minfeas) integer nempir,i,k real*8 empcdf(nempir),minfeas,x minfeas=0.0 x=1.0/nempir k=nempir/2.0 do i=1,k minfeas=minfeas+empcdf(i)*empcdf(nempir-i+1)*x*2.0 enddo if(mod(nempir,2).ne.0)then k=(nempir+1)/2.0 minfeas=minfeas+empcdf(k)*empcdf(k)*x endif return end subroutine drhomin(vals,prob,nempir,minfeas) c Calculates the minimum feasible bivariate expected value c for a discrete bivariate random variable with finite support. c (Common marginals) integer nempir,i,j real*8 minfeas,vals(nempir),prob(nempir),a,b, &p,q minfeas=0.0 i=nempir j=1 p=prob(nempir) q=prob(1) dowhile((j.le.nempir).and.(i.ge.1)) a=vals(i) b=vals(j) if(p.gt.q)then minfeas=minfeas+a*b*q p=p-q j=j+1 if(j.le.nempir)then q=prob(j) endif elseif(q.gt.p)then minfeas=minfeas+a*b*p q=q-p i=i-1 if(i.gt.0)then p=prob(i) endif else minfeas=minfeas+a*b*p i=i-1 j=j+1 if(j.le.nempir)then q=prob(j) endif if(i.gt.0)then p=prob(i) endif endif enddo return 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 c c This file contains two versions of the algorithm. The first is that c published in the journal, including several auxiliary routines which c are not in the journal, but converted to double precision. c The second version is one which calls IMSL routines. c subroutine mulnor(a, b, sig, eps, n, inf, prob, bound, ifault) c c Algorithm AS 195 appl. statist. (1984) vol.33, no.1 c c Computes multivariate normal distribution function and computes c the probability that a multivariate normal vector falls in a c rectangle in n-space with error less than eps. c c Auxiliary routines required: BIVNOR (CACM algorithm 462) which c needs DERF, PPND = AS111 (or PPND7 or PPND16 from AS241), and c SYMINV = AS7 which calls CHOL = AS6. BIVNOR and DERF are c included in this file. c double precision * a(*), b(*), sig(*), c(7), d(7), co(25), sd(2), coef(5, 3), * binc(7, 5), bl(7, 5), br(7, 35, 5), r(36, 5), s(5, 27), * xss(6, 6), pr2(6), prep(6), preb(6), cv(5, 15), sinv(28), * condl(5), xm(5), cond(5), beta(5, 5), bh(5), fe(5), sigma(28), * ep(5), del(3, 5), bou4(5), bou5(5), ans(6), fact(5), prod(5), * bint(6), bcs(5), bcn(7), eo(13), do(25) integer inf(7), intvl(5), ind(5), ksa(5), num(5), itype(5),ifault logical simps, ipr double precision eps, bound, prob, cons2, cons3, cons4, epsmin, + epssim, zero, p05, p15, half, one, two, three, six, eight, + ten, twelve, fiften, forty5, ninety, nine45, bou1, bou2, bou3, + cup, det, eplos, epsi, ept, fac, fsa, rho, sigc, t, tem, temb, + temp, wb, wl, wt, wu, x1, x2, xs, y1, y2, z double precision bivnor, f2, f3, f4, f5, f6, vnorm2, phi, x, + y,pp,qq,pdf,pp2 data coef * /0.311111111111111D0, 1.422222222222222D0, 0.533333333333333D0, * 1.422222222222222D0, 0.311111111111111D0, 0.333333333333333D0, * 0.0D0, 1.333333333333333D0, 0.0D0, 0.333333333333333D0, 0.5D0, * 0.0D0, 0.0D0, 0.0D0, 0.5D0/ data bcs /1.0D0, 4.0D0, 6.0D0, 4.0D0, 1.D0/ data bcn /1.0D0, 6.0D0, 15.0D0, 20.0D0, 15.0D0, 6.0D0, 1.D0/ data do /-3.75043972D0,-3.3242574D0, -2.85697D0, -2.36675941D0, * -2.3344142D0, -1.8891759D0, -1.7320508D0, -1.3556262D0, * -1.15440539D0, -1.D0, -0.74196378D0, -0.61670659D0,0.0D0, * 0.61670659D0, 0.74196378D0, 1.0D0, 1.15440539D0,1.3556262D0, * 1.7320508D0, 1.8891759D0, 2.3344142D0, 2.36675941D0, 2.85697D0, * 3.3242574D0,3.75043972D0/ data eo / -2.85697D0, -2.3344142D0, -1.7320508D0, -1.3556262D0, * -1.0D0, -0.74196378D0, 0.0D0, 0.74196378D0, 1.0D0, 1.3556262D0, * 1.7320508D0, 2.3344142D0, 2.85697D0/ data cons2 /6.879833D0/, cons3 /4.517004D0/, cons4 /76.371214D0/, * epsmin /1.0d-8/, epssim /6.0d-5/ data zero, p05, p15, half, one, two, three, six, eight, ten, * twelve, fiften, forty5, ninety, nine45 /0.0D0, 0.05D0, 0.15D0, * 0.5D0, 1.0D0, 2.0D0, 3.0D0, 6.0D0, 8.0D0, 10.0D0, 12.0D0, * 15.0D0, 45.0D0, 90.D0, 945.D0/ c c Functions needed to calculate the first six derivatives of c the normal density c f2(x, y) = abs((x * x - one) / (y * y)) f3(x, y) = abs(-x * (x * x - three) / (y * y * y)) f4(x, y) = abs((three + x * x * (-six + x * x)) / (y ** 4)) f5(x, y) = abs(x * (-fiften + x * x * (ten - x * x))) / (y ** 5) f6(x, y) = abs(-fiften + x * x * (forty5 - x * x * (-fiften + * x * x))) / (y ** 6) c Checking for faulty data c ifault = 0 if (eps .le. epsmin) ifault = 2 if (n .le. 0 .or. n .gt. 7) ifault = 3 if (ifault .ne. 0) return do 1 i = 1, n 1 if (inf(i) .eq. 2 .and. a(i) .lt. b(i)) ifault = 100 + i if (ifault .ne. 0) return bound = eps c Finding z such that p(n(0,1).gt.z) .lt. 0.15*eps/n c co will contain the roots of the first 5 or 7 Hermite c polynomials ept = eps * p15 / float(n) z = -vnorm2(ept, ifault) + epsmin if (ifault .ne. 0) return call fnorm2(z,pp,qq,pdf) cup = qq c c inverting sig and the n-2 lower right hand principal minors c ik = 0 ij = 0 do 3 i = 1, n do 3 j = 1, i ik = ik + 1 if (i .eq. j) goto 2 ij = ij + 1 sigma(ik) = sig(ij) goto 3 2 sigma(ik) = one 3 continue if (n .le. 2) goto 4 call invert(sigma, sinv, cv, n, det, ifault) if (ifault .ne. 0) return simps = .true. if (det .lt. p05 .or. eps .le. epssim) simps = .false. prob = zero det = sinv(1) * sinv(3) - sinv(2) * sinv(2) sd(1) = sqrt(sinv(3) / det) sd(2) = sqrt(sinv(1) / det) rho = -sinv(2) / (sd(1) * sd(2) * det) if (abs(rho) .gt. one) goto 400 4 nm2 = n - 2 nm1 = n - 1 c c checking whether upper and lower endpoints are too big c eplos = zero do 6 l = 1, n c(l) = max(b(l), -z) d(l) = min(a(l), z) if (inf(l) .eq. 0) d(l) = z if (inf(l) .eq. 1) c(l) = -z if (a(l) .gt. z .or. inf(l) .eq. 0) eplos = eplos + cup if (b(l) .lt. -z .or. inf(l) .eq. 1) eplos = eplos + cup if (c(l) .ge. d(l)) return 6 continue if (n .eq. 1) goto 350 fac = one ipr = .false. if (inf(1) .ne. 1 .or. inf(2) .ne. 1) goto 7 ipr = .true. eplos = eplos - two * cup goto 8 7 if (inf(1) .ne. 0 .or. inf(2) .ne. 0) goto 8 fac = -one ipr = .true. d(1) = c(1) d(2) = c(2) eplos = eplos - two * cup 8 if (n .eq. 2) goto 360 ifault = 5 epsi = (eps - eplos) / float(nm2) c Finding regression coefficients (beta,bh) and bounds on the c conditional integrals (binc) c cond(l)=conditional variance of variable n-l+1 given later c variables c do 15 l = 1, nm2 cond(l) = one / sqrt(cv(l, 1)) condl(l) = log(cond(l)) do 10 i = 1, l beta(l, i) = zero do 10 j = 1, l jk = (l - i + 1) * (l - i) / 2 + j if (j .gt. l - i + 1) jk = j * (j - 1) / 2 + l - i + 1 jn = (n - l + j) * (n - l + j - 1) / 2 + n - l beta(l, i) = beta(l, i) + sigma(jn) * cv(l, jk) 10 continue k = n - l - 1 bh(k + 1) = beta(l, l) do 11 i = 1, k bh(i) = zero do 11 j = 1, l jn = (j + n - l) * (j + n - l - 1) / 2 + n - l ijk = 1 + j * (j - 1) / 2 bh(i) = bh(i) + sigma(jn) * cv(l, ijk) 11 continue k = 0 sigc = zero do 12 j = 1, l do 12 i = 1, j k = k + 1 sigc = sigc + bh(i) * bh(j) * sinv(k) 12 continue binc(1, l) = one binc(2, l) = sqrt(sigc) binc(3, l) = two * sigc binc(4, l) = cons2 * sigc * binc(2, l) binc(5, l) = twelve * sigc * sigc if (simps) goto 13 binc(6, l) = sigc * sigc * binc(1, l) * cons3 binc(7, l) = (sigc ** 3) * cons4 13 if (l .lt. nm2) goto 15 do 14 i = 1, nm2 bh(i) = zero do 14 j = 1, nm2 jk = (l - i + 1) * (l - i) / 2 + j if (j .gt. l - i + 1) jk = j * (j - 1) / 2 + l - i + 1 jn = (2 + j) * (1 + j) / 2 + 1 bh(i) = bh(i) + sigma(jn) * cv(l, jk) 14 continue 15 continue l = 1 c c co will contain the roots of the first 5 or 7 Hermite c polynomials c if (simps) goto 50 do 40 i = 1, 25 40 co(i) = do(i) iend = 25 ien = 7 goto 60 50 do 55 i = 1, 13 55 co(i) = eo(i) iend = 13 ien = 5 c c Initialising values. xss contains partial sums used for c calculating conditional means. c 60 do 70 i = 1, nm1 70 xss(i, 1) = zero xm(1) = zero prod(1) = one pr2(1) = one do 80 i = 1, nm2 ni = n - i + 1 pr2(i + 1) = pr2(i) * (d(ni) - c(ni)) 80 continue c c bint(l) is a bound on the error accumulated at levels l and c deeper. c ans(l) is the accumulated integral at level l. c prep(l) contains the integrand at level l. c preb(l) bounds the error accumulated at levels deeper than l. c bint(nm1) = zero 90 intvl(l) = 2 ans(l) = zero bou4(l) = zero bint(l) = zero prep(l) = zero preb(l) = zero k = 1 c c Finding which of the co are in the current interval. c s(l,.) are the endpoints of intervals on which the integrand c at level l and its derivatives are monotone. c num(l) is the number of such intervals. c nl = n - l + 1 s(l, 1) = c(nl) - xm(l) s(l, iend + 2) = d(nl) - xm(l) num(l) = iend + 2 do 91 i = 1, iend njs = i if (s(l, 1) .lt. co(i) * cond(l)) goto 92 num(l) = num(l) - 1 91 continue 92 if (num(l) .eq. 2) goto 99 do 94 i = njs, iend mjs = iend - i + njs if (s(l, iend + 2) .ge. co(mjs) * cond(l)) goto 96 num(l) = num(l) - 1 94 continue 96 if (num(l) .eq. 2) goto 99 do 98 i = njs, mjs inj = i - njs + 2 s(l, inj) = co(i) * cond(l) 98 continue 99 numl = num(l) s(l, numl) = s(l, iend + 2) c c ep(l) is an upper limit on the allowable error at level l. c r(k,l) is the right end-point of the current sub-interval. c fe(l) is the left end-point of the current sub-interval. c ind(l)-1 indicates which point of the Newton-Cotes formula c we are dealing with. c ep(l) = epsi / pr2(l + 1) r(1, l) = s(l, 2) ind(l) = 6 c c Bounding derivatives at left end-point of current interval. c bl(i,l) is a bound on the ith derivative of the normal c density at level l at the left end-point. c fe(l) = s(l, 1) t = fe(l) / cond(l) bl(1, l) = phi(t, condl(l)) bl(2, l) = bl(1, l) * abs(t / cond(l)) bl(3, l) = bl(1, l) * f2(t, cond(l)) bl(4, l) = bl(1, l) * f3(t, cond(l)) bl(5, l) = bl(1, l) * f4(t, cond(l)) if (simps) goto 100 bl(6, l) = bl(1, l) * f5(t, cond(l)) bl(7, l) = bl(1, l) * f6(t, cond(l)) c Bounding derivatives at right end-point of sub-interval. c br(i,l) is a bound on the ith derivative of the normal c density at level l at the right end-point. c 100 t = r(k, l) / cond(l) br(1, k, l) = phi(t, condl(l)) br(2, k, l) = br(1, k, l) * abs(t / cond(l)) br(3, k, l) = br(1, k, l) * f2(t, cond(l)) br(4, k, l) = br(1, k, l) * f3(t, cond(l)) br(5, k, l) = br(1, k, l) * f4(t, cond(l)) if (simps) goto 104 br(6, k, l) = br(1, k, l) * f5(t, cond(l)) br(7, k, l) = br(1, k, l) * f6(t, cond(l)) 104 r(k + 1, l) = (fe(l) + r(k, l)) * half bou5(l) = ep(l) * (r(k, l) - s(l, 1)) del(2, l) = r(k + 1, l) - fe(l) c c Checking the bound for the trapezoidal rule c del(3, l) = two * del(2, l) bou1 = max(br(1, k, l), bl(1, l)) * binc(3, l) + * two * max(br(2, k, l), bl(2, l)) * binc(2, l) + * max(br(3, k, l), bl(3, l)) * binc(1, l) bou3 = bou4(l) + bou1 * (del(3, l) ** 3) * prod(l) / twelve itype(l) = 3 if (bou3 .le. bou5(l)) goto 200 c c Checking the bound for Simpsons rule. c bou1 = zero do 110 ij = 1, 5 jk = 6 - ij bou2 = max(br(ij, k, l), bl(ij, l)) bou1 = bou1 + bou2 * binc(jk, l) * bcs(ij) 110 continue bou3 = bou4(l) + bou1 * (del(2, l) ** 5) * prod(l) / ninety itype(l) = 2 if (bou3 .le. bou5(l)) goto 200 if (simps) goto 130 c c Checking the bound for boules rule, if necessary. c del(1, l) = half * del(2, l) bou1 = zero itype(l) = 1 do 120 ij = 1, 7 jk = 8 - ij bou2 = max(br(ij, k, l), bl(ij, l)) bou1 = bou1 + bou2 * binc(jk, l) * bcn(ij) 120 continue bou3 = bou4(l) + bou1 * (del(1, l) ** 7) * prod(l) * eight / * nine45 if (bou3 .le. bou5(l)) goto 200 c c Sub-dividing further at level l when the bound is too big. c 130 k = k + 1 if (k .gt. 35) return goto 100 200 bint(l) = bint(l) + bou3 - bou4(l) bou4(l) = bou3 ksa(l) = k if (ind(l) .eq. 6) goto 202 if (itype(l) - 2) 205, 206, 210 c c The next 30 lines condition on the value xs and go to level l+1 c 202 ind(l) = 5 xs = fe(l) fact(l) = bl(1, l) 203 xss(nm1, l + 1) = xss(nm1, l) + bh(l) * (xs + xm(l)) do 204 ll = l, nm2 204 xss(ll, l + 1) = xss(ll, l) + beta(ll, l) * (xs + xm(l)) if (l .eq. nm2) goto 300 c c xm is the mean of the next variable given those fixed so far. c xm(l + 1) = xss(l, l + 1) prod(l + 1) = prod(l) * fact(l) l = l + 1 goto 90 205 ind(l) = 4 k = ksa(l) xs = half * (fe(l) + r(k + 1, l)) goto 207 206 ind(l) = 3 k = ksa(l) xs = r(k + 1, l) 207 t = xs / cond(l) fact(l) = phi(t,condl(l)) goto 203 208 ind(l) = 2 k = ksa(l) xs = half * (r(k, l) + r(k + 1, l)) goto 207 210 ind(l) = 1 k = ksa(l) xs = r(k, l) fact(l) = br(1, k, l) goto 203 c c Evaluate conditional bivariate probabilities at deepest level c 300 x1 = fac * (xss(nm1, nm1) - d(1)) / sd(1) x2 = fac * (xss(nm2, nm1) - d(2)) / sd(2) l = nm1 ans(l) = bivnor(x1, x2, rho) if (ipr) goto 310 y1 = (xss(nm1, nm1) - c(1)) / sd(1) y2 = (xss(nm2, nm1) - c(2)) / sd(2) wu = bivnor(y1, y2, rho) wt = bivnor(x1, y2, rho) wb = bivnor(y1, x2, rho) ans(l) = ans(l) + wu - wt - wb 310 if (l .eq. 1) goto 340 l = l - 1 c c Advancing the integration at the current level c indl = ind(l) numl = num(l) ity = itype(l) temp = fact(l) * ans(l + 1) temb = bint(l + 1) fsa = one if (indl .ne. 1) goto 315 tem = temp temp = temp + prep(l) temb = temb + preb(l) prep(l) = tem preb(l) = bint(l + 1) fsa = two 315 ans(l) = ans(l) + coef(indl, ity) * temp * del(ity, l) bint(l) = bint(l) + coef(indl, ity) * temb * del(ity, l) c c Making use of the error which did not accumulate at level l+1 c ep(l) = ep(l) + del(ity, l) * coef(indl, ity) * (fsa * float(nm2 - * l) * epsi / pr2(l + 1) - temb) / (s(l, numl) - s(l, 1)) if (indl .eq. 1) goto 320 igo = indl - (1 + (itype(l) * (itype(l) - 1)) / 2) goto (210, 208, 206, 205), igo c c Un-subdividing at level l. 320 k = ksa(l) do 322 i = 1, ien 322 bl(i, l) = br(i, k, l) ind(l) = 5 fe(l) = r(k, l) if (k .eq. 1) goto 326 k = k - 1 goto 104 326 if (intvl(l) .eq. num(l)) goto 310 intvl(l) = intvl(l) + 1 intl = intvl(l) r(1, l) = s(l, intl) goto 100 c c Completion of integration and bounding. c 340 ifault = 0 prob = ans(1) bound = bint(1) + eplos return c c special cases -- c label 350 - n=1 c label 360 - n=2 c 350 call fnorm2(d(1),pp,qq,pdf) call fnorm2(c(1),pp2,qq,pdf) prob = pp-pp2 return 360 rho = sigma(2) if (abs(rho) .gt. one) goto 400 y1 = -d(1) * fac y2 = -d(2) * fac prob = bivnor(y1, y2, rho) if (ipr) return x1 = -c(1) x2 = -c(2) wl = bivnor(x1, x2, rho) wt = bivnor(x1, y2, rho) wb = bivnor(y1, x2, rho) prob = wl - wt - wb + prob return c c Error return for covariance not positive definite. c 400 ifault = 4 return end c c c double precision function phi(x, y) c c algorithm as 195.1 appl. statist. (1984) vol.33, no.1 c c computes univariate normal density c xlow=log(smallest floating point number) c sq2p=log(sqrt(two*pi)) c double precision arg, half, sq2p, x, xlow, y, zero c data xlow /-87.0D0/, sq2p /0.91893853320467274D0/, zero /0.D0/, * half /0.5D0/ phi = zero arg = -half * x * x - sq2p - y if (arg .gt. xlow) phi = exp(arg) return end c c c DOUBLE PRECISION FUNCTION BIVNOR(AH, AK, R) C C BIVNOR is a controlled precision Fortran function to calculate C the bivariate normal upper right area, viz. the probability for C two normal variates X and Y whose correlation is R, that X > AH C and Y > AK. C The accuracy is specified as the number of decimal digits, IDIG. C C Reference: C Donnelly, T.G. (1973) Algorithm 462, Bivariate normal C distribution, Comm. A.C.M., vol.16, p. 638. C C Corrected for the case AH = 0, AH non-zero by reversing AH & AK C in such cases. C DOUBLE PRECISION AH, AK, R C C Local variables C DOUBLE PRECISION TWOPI, B, XAH, XAK, GH, GK, RR, GAUSS, DERF, H2, + A2, H4, EX, W2, AP, S2, SP, S1, SN, SGN, SQR, CON, WH, WK, GW, + T, G2, CONEX, CN, TWO, ZERO, ONE, FOUR, QUART, HALF, EXPLIM INTEGER IDIG, IS DATA TWO/2.D0/, ZERO/0.D0/, ONE/1.D0/, FOUR/4.D0/, QUART/0.25D0/, + HALF/0.5D0/ GAUSS(T) = (ONE + DERF(T/SQRT(TWO)))/TWO C C GAUSS is a univariate lower normal tail area calculated here from C the central error function, DERF. C It may be replaced by the algorithm in Hill, I.D. and Joyce, S.A. C Algorithm 304, Normal curve integral (S15), Comm. A.C.M. (10), C (June 1967), p.374 or with Applied Statistics algorithm AS66. C Source: Owen, D.B. Ann. Math. Statist., vol.27 (1956), p.1075. C DATA TWOPI/6.2831 85307 17958 7D0/, IDIG/15/, EXPLIM/80.D0/ C B = ZERO IF (AH .EQ. ZERO) THEN XAH = AK XAK = AH ELSE XAH = AH XAK = AK END IF GH = GAUSS(-XAH) / TWO GK = GAUSS(-XAK) / TWO IF (R .EQ. ZERO) THEN B = FOUR * GH * GK GO TO 350 END IF RR = ONE - R*R IF (RR .LT. ZERO) THEN WRITE(*, *)'Error in BIVNOR, R = ', R GO TO 390 END IF IF (RR .GT. ZERO) GO TO 100 C C R^2 = 1.0 C IF (R .GE. ZERO) GO TO 70 IF (XAH + XAK .GE. ZERO) GO TO 350 B = TWO * (GH + GK) - ONE GO TO 350 70 IF (XAH - XAK .LT. ZERO) THEN B = TWO * GK ELSE B = TWO * GH END IF GO TO 350 C C Regular case, R^2 < 1 C 100 SQR = SQRT(RR) IF (IDIG .EQ. 15) THEN CON = TWOPI * 1.D-15 / TWO ELSE CON = TWOPI / TWO / 10**IDIG END IF IF (XAH .NE. ZERO) GO TO 170 IF (XAK .NE. ZERO) GO TO 190 B = ATAN(R/SQR) / TWOPI + QUART GO TO 350 170 B = GH IF (XAH*XAK) 180, 200, 190 180 B = B - HALF 190 B = B + GK 200 WH = -XAH WK = (XAK/XAH - R)/SQR GW = TWO * GH IS = -1 210 SGN = -ONE T = ZERO IF (WK .EQ. ZERO) GO TO 320 IF (ABS(WK) - ONE) 270, 230, 240 230 T = WK * GW * (ONE - GW) / TWO GO TO 310 240 SGN = -SGN WH = WH * WK G2 = GAUSS(WH) WK = ONE / WK IF (WK .LT. ZERO) B = B + HALF B = B - (GW + G2)/TWO + GW*G2 270 H2 = WH * WH A2 = WK * WK H4 = H2 / TWO IF (H4 .LT. EXPLIM) THEN EX = EXP(-H4) ELSE EX = ZERO END IF W2 = H4 * EX AP = ONE S2 = AP - EX SP = AP S1 = ZERO SN = S1 CONEX = ABS(CON / WK) GO TO 290 280 SN = SP SP = SP + ONE S2 = S2 - W2 W2 = W2 * H4 / SP AP = -AP*A2 290 CN = AP * S2 / (SN + SP) S1 = S1 + CN IF (ABS(CN) - CONEX .GT. ZERO) GO TO 280 T = (ATAN(WK) - WK*S1) / TWOPI 310 B = B + SGN*T 320 IF (IS .GE. 0) GO TO 350 IF (XAK .NE. ZERO) THEN WH = -XAK WK = (XAH/XAK - R) / SQR GW = TWO * GK IS = 1 GO TO 210 END IF 350 IF (B .LT. ZERO) B = ZERO IF (B .GT. ONE) B = ONE 390 BIVNOR = B C RETURN END C C C SUBROUTINE INVERT(A, AINV, C, N, DET, IFAULT) C C Invert the NxN symmetric matrix, A, of which only the lower C triangle is stored by rows. The inverse is returned in AINV. C The inverse of the last I rows and columns, for I = 1, 2, ..., C N-2, is returned in C(I,*). DET = the determinant of A. C IFAULT = 4 if A is not positive definite. C DOUBLE PRECISION A(*), AINV(*), C(5,*), DET INTEGER N, IFAULT C C Local variables C DOUBLE PRECISION ONE, WK(15), W(7) INTEGER NULLTY, NN, I, ROW1, WKPOS, ROW, APOS, COL DATA ONE/1.D0/ C DO 50 I = 1, N-2 C C Copy the last I rows and columns of A into WK. C NN = I * (I + 1) / 2 ROW1 = N + 1 - I WKPOS = 1 DO 20 ROW = ROW1, N APOS = ROW * (ROW - 1) / 2 + ROW1 DO 10 COL = ROW1, ROW WK(WKPOS) = A(APOS) WKPOS = WKPOS + 1 APOS = APOS + 1 10 CONTINUE 20 CONTINUE C C Call SYMINV to invert WK in situ. C CALL SYMINV(WK, I, NN, WK, W, NULLTY, IFAULT) IF (IFAULT .NE. 0) GO TO 70 C C Copy the inverse into the I-th row of C. C DO 40 APOS = 1, NN 40 C(I, APOS) = WK(APOS) 50 CONTINUE C C Use AS6 (CHOL) to calculate the determinant. C This is inefficient as AS7 will call AS6 again to repeat the same C calculations. AS6/7 can easily be converted to calculate the C determinant. C NN = N * (N + 1) / 2 CALL CHOL(A, N, NN, AINV, NULLTY, IFAULT) IF (IFAULT .NE. 0) GO TO 70 DET = ONE APOS = 1 DO 60 ROW = 1, N DET = DET * AINV(APOS) APOS = APOS + ROW + 1 60 CONTINUE DET = DET * DET C C Invert the full matrix. C CALL SYMINV(A, N, NN, AINV, W, NULLTY, IFAULT) IF (IFAULT .EQ. 0) RETURN C 70 IFAULT = 4 RETURN END C C C SUBROUTINE CALERF(ARG,RESULT,JINT) C------------------------------------------------------------------ C C THIS PACKET COMPUTES THE ERROR AND COMPLEMENTARY ERROR FUNCTIONS C FOR REAL ARGUMENTS ARG. IT CONTAINS TWO FUNCTION TYPE C SUBPROGRAMS, ERF AND ERFC (OR DERF AND DERFC), AND ONE C SUBROUTINE TYPE SUBPROGRAM, CALERF. THE CALLING STATEMENTS C FOR THE PRIMARY ENTRIES ARE C C Y=ERF(X) (OR Y=DERF(X) ) C AND C Y=ERFC(X) (OR Y=DERFC(X) ). C C THE ROUTINE CALERF IS INTENDED FOR INTERNAL PACKET USE ONLY, C ALL COMPUTATIONS WITHIN THE PACKET BEING CONCENTRATED IN THIS C ROUTINE. THE FUNCTION SUBPROGRAMS INVOKE CALERF WITH THE C STATEMENT C CALL CALERF(ARG,RESULT,JINT) C WHERE THE PARAMETER USAGE IS AS FOLLOWS C C FUNCTION PARAMETERS FOR CALERF C CALL ARG RESULT JINT C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 C ERFC(ARG) ABS(ARG) .LT. XMAX ERFC(ARG) 1 C C THE MAIN COMPUTATION EVALUATES NEAR MINIMAX APPROXIMATIONS C FROM "RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE ERROR FUNCTION" C BY W. J. CODY, MATH. COMP., 1969, PP. 631-638. THIS C TRANSPORTABLE PROGRAM USES RATIONAL FUNCTIONS THAT THEORETICALLY C APPROXIMATE ERF(X) AND ERFC(X) TO AT LEAST 18 SIGNIFICANT C DECIMAL DIGITS. THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC C SYSTEM, THE COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER C SELECTION OF THE MACHINE-DEPENDENT CONSTANTS. C C******************************************************************* C C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS C C XSMALL = ARGUMENT BELOW WHICH ERF(X) MAY BE REPRESENTED C BY 2*X/SQRT(PI) AND ABOVE WHICH X*X WILL C NOT UNDERFLOW. A CONSERVATIVE VALUE IS THE C LARGEST X SUCH THAT 1.0 + X = 1.0 TO MACHINE C PRECISION. C XMAX = LARGEST ARGUMENT ACCEPTABLE TO ERFC; SOLUTION TO C EQUATION: W(X) * (1-0.5/X**2) = XMIN, WHERE C W(X) = EXP(-X*X)/(X*SQRT(PI)), AND XMIN IS THE C SMALLEST POSITIVE MACHINE NUMBER (SEE TABLE BELOW). C C APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: C C XSMALL XMAX XMIN C C IBM 195 (D.P.) 1.39D-17 13.306 5.40D-79 C CDC 7600 (S.P.) 7.11E-15 25.922 3.13E-294 C CRAY-1 (S.P.) 7.11E-15 75.326 4.58E-2467 C UNIVAC 1108 (D.P.) 1.73D-18 26.582 2.78D-309 C VAX 11/780 (S.P.) 5.96E-8 9.269 2.94E-39 C VAX 11/780 (D.P.) 1.39D-17 9.269 2.94D-39 C IBM PC (S.P.) 5.96E-8 9.194 1.18E-38 C IBM PC (D.P.) 1.11D-16 26.543 2.23D-308 C C******************************************************************* C C ERROR RETURNS C C THE PROGRAM RETURNS ERFC = 0 FOR ARG .GT. XMAX. C C C OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) C C ABS, EXP C C OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) C C DABS, DEXP C C C AUTHOR: W. J. CODY C MATHEMATICS AND COMPUTER SCIENCE DIVISION C ARGONNE NATIONAL LABORATORY C ARGONNE, IL 60439 C C LATEST MODIFICATION: JANUARY 8, 1985 C C------------------------------------------------------------------ INTEGER I,JINT CS REAL A,ARG,B,C,D,FOUR,HALF,P,ONE,Q,RESULT,SQRPI, CS 1 TWO,THRESH,X,XMAX,XDEN,XNUM,XSMALL,Y,YSQ,ZERO DOUBLE PRECISION A,ARG,B,C,D,FOUR,HALF,P,ONE,Q,RESULT,SQRPI, 1 TWO,THRESH,X,XMAX,XDEN,XNUM,XSMALL,Y,YSQ,ZERO DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) C------------------------------------------------------------------ C MATHEMATICAL CONSTANTS C------------------------------------------------------------------ CS DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/ CS DATA SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/ DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/ DATA SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/ C------------------------------------------------------------------ C MACHINE-DEPENDENT PARAMETERS C------------------------------------------------------------------ CS DATA XSMALL/4.2E-16/, XMAX/9.269E0/ DATA XSMALL/4.2D-16/, XMAX/9.269D0/ C------------------------------------------------------------------ C COEFFICIENTS FOR APPROXIMATION TO DERF IN FIRST INTERVAL C------------------------------------------------------------------ CS DATA A/3.16112374387056560E00,1.13864154151050156E02, CS 1 3.77485237685302021E02,3.20937758913846947E03, CS 2 1.85777706184603153E-1/ CS DATA B/2.36012909523441209E01,2.44024637934444173E02, CS 1 1.28261652607737228E03,2.84423683343917062E03/ 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/ C------------------------------------------------------------------ C COEFFICIENTS FOR APPROXIMATION TO DERFC IN SECOND INTERVAL C------------------------------------------------------------------ CS DATA C/5.64188496988670089E-1,8.88314979438837594E0, CS 1 6.61191906371416295E01,2.98635138197400131E02, CS 2 8.81952221241769090E02,1.71204761263407058E03, CS 3 2.05107837782607147E03,1.23033935479799725E03, CS 4 2.15311535474403846E-8/ CS DATA D/1.57449261107098347E01,1.17693950891312499E02, CS 1 5.37181101862009858E02,1.62138957456669019E03, CS 2 3.29079923573345963E03,4.36261909014324716E03, CS 3 3.43936767414372164E03,1.23033935480374942E03/ 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/ C------------------------------------------------------------------ C COEFFICIENTS FOR APPROXIMATION TO DERFC IN THIRD INTERVAL C------------------------------------------------------------------ CS DATA P/3.05326634961232344E-1,3.60344899949804439E-1, CS 1 1.25781726111229246E-1,1.60837851487422766E-2, CS 2 6.58749161529837803E-4,1.63153871373020978E-2/ CS DATA Q/2.56852019228982242E00,1.87295284992346047E00, CS 1 5.27905102951428412E-1,6.05183413124413191E-2, CS 2 2.33520497626869185E-3/ 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/ C------------------------------------------------------------------ X = ARG CS Y = ABS(X) Y = DABS(X) IF (Y .GT. FOUR) GO TO 200 IF (Y .GT. THRESH) GO TO 100 C------------------------------------------------------------------ C EVALUATE ERF FOR ABS(X) .LE. 0.46875 C------------------------------------------------------------------ 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 GO TO 800 C------------------------------------------------------------------ C EVALUATE ERFC FOR 0.46875 .LT. ABS(X) .LE. 4.0 C------------------------------------------------------------------ 100 YSQ = Y * Y XNUM = C(9)*Y XDEN = Y DO 120 I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y 120 CONTINUE CS RESULT = EXP(-YSQ) * (XNUM + C(8)) / (XDEN + D(8)) RESULT = DEXP(-YSQ) * (XNUM + C(8)) / (XDEN + D(8)) GO TO 300 C------------------------------------------------------------------ C EVALUATE ERFC FOR ABS(X) .GT. 4.0 C------------------------------------------------------------------ 200 RESULT = ZERO IF (Y .GE. XMAX) GO TO 300 220 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)) CS RESULT = (EXP(-Y*Y) / Y) * (SQRPI - RESULT) RESULT = (DEXP(-Y*Y) / Y) * (SQRPI - RESULT) C------------------------------------------------------------------ C FIX UP FOR NEG. ARG., ERF, ETC. C------------------------------------------------------------------ 300 IF (JINT .EQ. 0) GO TO 350 IF (X .LT. ZERO) RESULT = TWO - RESULT GO TO 800 350 RESULT = (HALF - RESULT) + HALF IF (X .LT. ZERO) RESULT = -RESULT C------------------------------------------------------------------ 800 RETURN C---------- LAST CARD OF CALERF ---------- END CS REAL FUNCTION ERF(X) DOUBLE PRECISION FUNCTION DERF(X) C------------------------------------------------------------------ C C PROGRAM TO COMPUTE THE ERROR FUNCTION C C AUTHOR - W. J. CODY C C DATE - JANUARY 8, 1985 C C------------------------------------------------------------------ INTEGER JINT CS REAL X, RESULT DOUBLE PRECISION X, RESULT C------------------------------------------------------------------ JINT = 0 CALL CALERF(X,RESULT,JINT) CS ERF = RESULT DERF = RESULT RETURN C---------- LAST CARD OF DERF ---------- END CS REAL FUNCTION ERFC(X) DOUBLE PRECISION FUNCTION DERFC(X) C------------------------------------------------------------------ C C PROGRAM TO COMPUTE THE COMPLEMENTARY ERROR FUNCTION C C AUTHOR - W. J. CODY C C DATE - JANUARY 8, 1985 C C------------------------------------------------------------------ INTEGER JINT CS REAL X, RESULT DOUBLE PRECISION X, RESULT C------------------------------------------------------------------ JINT = 1 CALL CALERF(X,RESULT,JINT) CS ERFC = RESULT DERFC = RESULT RETURN C---------- LAST CARD OF DERFC ---------- END c c------------------------------------------------------------------------- c 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 upper triangle, u( ), such that uprime * 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 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.ABS(X*A(KK))) RETURN U(K)=ZERO 40 CONTINUE 50 IF (ABS(W).LE.ABS(ETA*A(K))) GO TO 60 IF (W.LT.ZERO) RETURN U(K)=SQRT(W) GO TO 70 60 U(K)=ZERO NULLTY=NULLTY+1 70 J=J+ICOL 80 CONTINUE IFAULT=0 END C