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