c c--------------------------------------------------------------------- c c 9/26/1997 "m-estim.f" M. Palumbo c c--------------------------------------------------------------------- c c c The program estimates means and covariances for an artificial c dataset to be made available to persons interested in implementing c our Imputation-by-Simulation Algorithm, as described in the c Advances in Econometrics volume of 1998. c c c This fortran code reads monte carlo generated data from "monte.dat" c and takes estimation parameters from "monte.parms" c c c The fortran program "m-simul.f" takes this program's output and c uses simulation procedures to impute missing values. c c c--------------------------------------------------------------------- c c main program c implicit real*8(a-h,o-z) parameter(nfacs=1000,nrow=3504,nidx=6,nxt=50,ntype=4,nparis=14, $ nvars=122,norder=7,ncut=norder+1) real anormc,anormd,arg,aback dimension ndev(nvars),nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xmean(nvars), $ xcut(nvars,ncut),xsum(nvars), $ xcov(nvars,nvars),xcovt(nvars,nvars),xcovn(nvars,nvars), $ xnum(nvars),xfreq(nvars,norder),xreal(nfacs,nvars) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut common/stand/ystand(nvars,203),xstand(nvars,203) c external anormc,anormd,fn,bivn,corr,datget,golden, $ func2,func3,func4,func5,func6,standard,undo c data pi/3.1415927/ data tol/0.000001/ c c--------------------------------------------------------------------- c c set parameter "imonte" = 0 for run using real Jamaican data; c = 1 for first monte carlo simulation; c = 13 for second monte carlo simulation. c c cmonte imonte = 0 cmonte imonte = 13 c imonte = 1 c c--------------------------------------------------------------------- c if (imonte .eq. 0) open(unit=10,file="m-estim.parms") if (imonte .eq. 13) open(unit=10,file="m-estim.parms") if (imonte .eq. 1) open(unit=10,file="monte.parms") if (imonte .eq. 0) open(unit=11,file="jamreal.dat") if (imonte .eq. 1) open(unit=11,file="monte.dat") if (imonte .eq. 13) open(unit=11,file="jamtest.dat") c c--------------------------------------------------------------------- c c read Monte Carlo data from "monte.dat". c call datget c c--------------------------------------------------------------------- c c program needs fresh array, "nxkind." c rewind(10) read(10,*) iestflg,iacrj,itype,nfac,nvar read(10,*) (nxkind(i),i=1,nvar) read(10,*) (nxcut(i),i=1,nvar) c c--------------------------------------------------------------------- c c open some output or input files. c c open different files depending on "itype" of facility. c if(iestflg.eq.1)then if (itype .eq. 1) then open(unit=8,file='xcov.datfirst') open(unit=12,file="xmean.dat") open(unit=13,file="xcut.dat") open(unit=14,file="xcov.dat") open(unit=15,file="m-estim.out") endif endif c c--------------------------------------------------------------------- c c if "iestflg = 0", then mean/covariances already computed. c if (iestflg.eq.0) goto 650 c c--------------------------------------------------------------------- c c estimate 1st and 2nd moements of variables by type only. c do 600 i = itype,itype c it = i ip = 99 c c------------------------------------------------------------------------- c c initialize moment arrays. c do 240 k = 1,nvars xmean(k) = 0. ndev(k) = 0 do 220 l = 1,norder 220 xfreq(k,l) = 0. do 230 l = 1,ncut 230 xcut(k,l) = 0. do 240 l = 1,nvars 240 xcov(k,l) = 0. c do 400 k = 1,nvar xsum(k) = 0. xnum(k) = 0. xmean(k) = 0. c nxki = nxkind(k) c do 300 l = 1,nfac c if (nid(l,1) .eq. it) then c nonm = nonmis(l,k) if (nonm .eq. 1) then c c count number of facilities with nonmissing data. c xnum(k) = xnum(k) + 1. c c add up values of binary variables. c if (nxki .eq. 1) then if (nxbnry(l,k) .eq. 1) then xsum(k) = xsum(k) + 1. endif endif c c count frequencies for ordered discrete variables. c if (nxki .eq. 2) then nord = nxordd(l,k) xfreq(k,nord) = xfreq(k,nord) + 1. endif c c add up values of continuous variables. c if (nxki .eq. 3) then xval = xreal(l,k) xsum(k) = xsum(k) + xval endif endif endif c 300 continue c c compute the means for each variable. c if (xnum(k) .gt. 0.5) then if (nxki .eq. 1) then arg = xsum(k)/xnum(k) if (arg .lt. 0.0000001) arg = 0.0000001 if (arg .gt. 0.9999999) arg = 0.9999999 call mdnris(arg,aback,ier) xmean(k) = aback endif if (nxki .eq. 2) then xsum(k) = xfreq(k,1) arg = xsum(k)/xnum(k) if (arg .lt. 0.0000001) arg = 0.0000001 if (arg .gt. 0.9999999) arg = 0.9999999 call mdnris(arg,aback,ier) xmean(k) = -aback endif if (nxki .eq. 3) then xmt = xsum(k)/xnum(k) xmean(k) = xmt endif endif c c compute cut-off points for binary and ordered discrete variables. c if (xnum(k) .gt. 0.5) then if (nxki .eq. 1) then xcut(k,1) = -20. xcut(k,2) = 0. xcut(k,3) = 20. endif if (nxki .eq. 2) then xcut(k,1) = -20. xcut(k,2) = 0. do 310 l = 3,nxcut(k)-1 arg = xcut(k,l-1) - xmean(k) xmt1 = anormc(arg) xmt2 = xfreq(k,l-1)/xnum(k) arg = xmt1 + xmt2 if (arg .lt. 0.0000001) arg = 0.0000001 if (arg .gt. 0.9999999) arg = 0.9999999 call mdnris(arg,aback,ier) xcut(k,l) = xmean(k) + aback 310 continue xcut(k,nxcut(k)) = 20. do 312 l = nxcut(k)-1,3, -1 xcutl = xcut(k,l-1) xcutu = xcut(k,l) xcutd = xcutu - xcutl if (xcutd .lt. 0.0000001) then xcut(k,l) = 20.0 endif if (xcutd .gt. 0.0000001) goto 400 312 continue endif endif c 400 continue c c--------------------------------------------------------------------- c c write some output. c 325 format(//,4x,'output from m-estim.f',//,'Variable Means') 335 format(/,4x,'Facility Type: ',i3,/) 345 format(6x,' Var ',8x,' Kind ',8x,' Mean ',8x,' N ') 355 format(7x,i3,9x,2x,i2,2x,8x,f8.4,8x,f5.0) c write(15,325) write(15,335) i write(15,345) write(15,355) (k,nxkind(k),xmean(k),xnum(k),k=1,nvar) c write(12,625) (xmean(k),k=1,nvar) write(13,625) ((xcut(k,l),l=1,ncut),k=1,nvar) c c--------------------------------------------------------------------- c c compute variance/covariance matrix for each type/parish. c c if testing means only, set "istop1 = 1". c istop1 = 0 if (istop1 .eq. 1) goto 600 c do 500 k = 1,nvar do 500 l = k,1, -1 ii = k ij = l c nxki = nxkind(k) nxkj = nxkind(l) c c if testing variances only, set "istop2 = 1". c istop2 = 0 if (istop2 .eq. 1) then if (k .ne. l) goto 500 endif c c count numbers of observations for covariances. c count = 0. do 420 m = 1,nfac if (nid(m,1).eq.it) then inmi = nonmis(m,k) if (inmi .eq. 1) then inmj = nonmis(m,l) if (inmj .eq. 1) count = count + 1. endif endif 420 continue c c case 0: assign unit (or zero) variances to discrete random variables. c if (k .eq. l) then c c discrete rv's can get zero variances if they don't vary. c if ((nxki .eq. 1) .or. (nxki .eq. 2)) then if (count .gt. 0.) then nbig = -100 nlit = 100 do 430 m = 1,nfac if (nid(m,1).eq.it) then if (nonmis(m,k) .eq. 0) go to 430 if (nxki .eq. 1) nval = nxbnry(m,k) if (nxki .eq. 2) nval = nxordd(m,k) if (nval .ge. nbig) nbig = nval if (nval .le. nlit) nlit = nval endif 430 continue ndev(k) = nbig - nlit if (ndev(k) .lt. 0.01) then print *,'non-varying variables follow: ' print *,'var. is: ',k,' nval is: ',nval xcov(k,k) = 1.0 do 431 m = 1,nxcut(k) if (m .le. nval) xcut(k,m) = -20.0 if (m .gt. nval) xcut(k,m) = 20.0 431 continue endif if (ndev(k) .ge. 0.01) xcov(k,k) = 1.0 endif endif c endif c c case 1: both variables approximately continuous. c if ((nxki .eq. 3) .and. (nxkj .eq. 3)) then xcount = 0. do 440 m = 1,nfac if (nid(m,1).eq.it) then xk = xreal(m,k) - xmean(k) xnk = dfloat(nonmis(m,k)) xl = xreal(m,l) - xmean(l) xnl = dfloat(nonmis(m,l)) xcount = xcount + xnk*xnl xcov(k,l) = xcov(k,l) + xk*xl*xnk*xnl endif 440 continue if (xcount .gt. 1.5) then xcov(k,l) = xcov(k,l)/(xcount - 1.) endif if (xcount .lt. 1.5) xcov(k,l) = 0. endif c c following complicated cases only involve covariances. c if ((k .ne. l) .and. (count .gt. 1.5)) then c c covariances should equal zero if either var. doesn't vary. c if ((ndev(k) .eq. 0) .or. (ndev(l) .eq. 0)) $ goto 500 c xsigi = xcov(k,k)**(0.5) xsigj = xcov(l,l)**(0.5) if (xsigi .lt. 0.0001) go to 500 if (xsigj .lt. 0.0001) go to 500 c ax = -1. bx = 1. c c case 2: both variables binary. c if ((nxki .eq. 1) .and. (nxkj .eq. 1)) then call mnbrak(ax,bx,cx,fa,fb,fc,func2) xmin = ax if (ax .lt. 10.) then fmin = golden(ax,bx,cx,func2,tol,xmin) endif r = corr(xmin) xcov(k,l) = r endif c c case 3: one continuous variable, one binary variable. c if (((nxki .eq. 3) .and. (nxkj .eq. 1)) .or. $ ((nxki .eq. 1) .and. (nxkj .eq. 3))) then call mnbrak(ax,bx,cx,fa,fb,fc,func3) xmin = ax if (ax .lt. 10.) then fmin = golden(ax,bx,cx,func3,tol,xmin) endif r = corr(xmin) xsig = dsqrt(xcov(k,k)) xcov(k,l) = xsig*r endif c c case 4: both variables ordered discrete. c if ((nxki .eq. 2) .and. (nxkj .eq. 2)) then call mnbrak(ax,bx,cx,fa,fb,fc,func4) xmin = ax if (ax .lt. 10.) then fmin = golden(ax,bx,cx,func4,tol,xmin) endif r = corr(xmin) xcov(k,l) = r endif c c case 5: one continuous variable, one ordered discrete. c if (((nxki .eq. 3) .and. (nxkj .eq. 2)) .or. $ ((nxki .eq. 2) .and. (nxkj .eq. 3))) then call mnbrak(ax,bx,cx,fa,fb,fc,func5) xmin = ax if (ax .lt. 10.) then fmin = golden(ax,bx,cx,func5,tol,xmin) endif r = corr(xmin) xsig = dsqrt(xcov(k,k)) xcov(k,l) = xsig*r endif c c case 6: one binary variable, one ordered discrete. c if (((nxki .eq. 1) .and. (nxkj .eq. 2)) .or. $ ((nxki .eq. 2) .and. (nxkj .eq. 1))) then call mnbrak(ax,bx,cx,fa,fb,fc,func6) xmin = ax if (ax .lt. 10.) then fmin = golden(ax,bx,cx,func6,tol,xmin) endif r = corr(xmin) xcov(k,l) = r endif c endif c 500 continue c c--------------------------------------------------------------------- c c fill in the "other" side of the covariance matrix. c do 510 k = 1,nvar do 510 l = k, 1, -1 510 xcov(l,k) = xcov(k,l) c c--------------------------------------------------------------------- c c impose positive definiteness on estimated covariance matrix. c c open(unit=8,file="xcov.dat2first") write(8,625) ((xcov(k,l),l=1,nvar),k=1,nvar) do 511 k = 1,nvar do 511 l = 1,nvar 511 xcovt(k,l) = xcov(k,l) if (itype .eq. 1) call eigchk1(itype,xcovt,xcovn) 513 continue do 512 k = 1,nvar do 512 l = 1,nvar 512 xcov(k,l) = xcovn(k,l) c c--------------------------------------------------------------------- c c write estimated variances to output file. c 535 format(/,4x,'Facility Type: ',i3,/) 545 format(6x,' Var ',8x,' Kind ',8x,'Variance',8x,' N ') 555 format(7x,i3,9x,2x,i2,2x,6x,f10.4,8x,f5.0) write(15,535) i write(15,545) write(15,555) (k,nxkind(k),xcov(k,k),xnum(k),k=1,nvar) c write(14,625) ((xcov(k,l),l=1,nvar),k=1,nvar) c c--------------------------------------------------------------------- c 600 continue c c--------------------------------------------------------------------- c 625 format(/,5(2x,g13.6)) c c--------------------------------------------------------------------- c c run steve's program to simulate missing observations. c close(8) close(12) close(13) close(14) close(15) c 650 call simul c charry 650 continue c c--------------------------------------------------------------------- c stop end c c====================================================================== c function anormc(arg) c c---------------------------------------------------------------------- c c this function evaluates the standard normal distribution c function c implicit real(a-h,o-z) c c standardize the argument of the distribution function c xt=arg*.707107 c c evaluate using the error function c if(xt.le.-20.)anormc=0. if((xt.lt.0.).and.(xt.gt.-20.))anormc=1.-(.5*erfc(xt)) if((xt.ge.0.).and.(xt.lt.20.))anormc=.5*erfc(-xt) if(xt.ge.20.)anormc=1. return end c c---------------------------------------------------------------------- c function anormd(arg) c c---------------------------------------------------------------------- c c this function evaluates the standard normal density function c implicit real(a-h,o-z) data sqpi2/2.506628275/ x=-(arg*arg)/2. if(x.lt.-50.)x=-50. anormd=exp(x)/sqpi2 return end c c---------------------------------------------------------------------- c c John Antel mailed this file, "BIVNRM.F". c c it is not actually the distribution function but rather c the often used l function as in johnson and kotz c function returns prob(x>a,y>b) where x,y standard normal c rvs with corr rho c c the first subroutine fn calcs the standard normal density and c distribution function. c c checked out close to owen's 1959 tables bureau of standards c subroutine fn(x,d,dd) implicit real*8 (a-h,o-z) d=0 if(x.gt.5.0d0) x=5.0d0 if(x.le.-5.0d0) x=-5.0d0 arg1=x/dsqrt(2.0d0) dd =(derf(arg1)+1.0d0)/2.0d0 d=dexp(-.5d0*x*x)/dsqrt(6.28318530717959d0) return end c c double precision function bivn(ah,ak,r) implicit real*8(a-h,o-z) twopi=6.283185307179587d0 if(ah.le.-5.0d0) ah=-5.0d0 if(ah.ge.5.0d0) ah=5.0d0 if(ak.le.-5.0d0) ak=-5.0d0 if(ak.ge.5.0d0) ak=5.0d0 if(r.gt.0.996d0) r=.996d0 if(r.lt.-0.996d0) r=-.996d0 b=0.d0 x1=-ah x2=-ak call fn(x1,d1,gh) call fn(x2,d1,gk) gh=gh/2.d0 gk=gk/2.d0 if (r) 10, 30, 10 10 rr=1.d0-r**2 if(rr) 20,40,100 20 print 9999,r 9999 format(12h bivnor r is, d26.16) stop 30 b=4.d0*gh*gk go to 350 40 if(r) 50,70,70 50 if(ah+ak) 60,350,350 60 b=2.d0*(gh+gk)-1.d0 go to 350 70 if(ah-ak) 80,90,90 80 b=2.d0*gk go to 350 90 b=2.d0*gh go to 350 100 sqr=dsqrt(rr) 110 con=twopi*1.d-15/2 go to 140 140 if(ah) 170,150,170 150 if(ak) 190,160,190 160 b=datan(r/sqr)/twopi+.25d0 go to 350 170 b=gh if(ah*ak) 180,200,190 180 b=b-.5d0 190 b=b+gk if(ah) 200,340,200 200 wh=-ah wk=(ak/ah-r)/sqr gw=2.d0*gh is=-1 210 sgn=-1.d0 t=0.d0 if(wk)220,230,220 220 if(dabs(wk)-1.d0) 270,230,240 230 t=wk*gw*(1.d0-gw)/2.d0 go to 310 240 sgn=-sgn wh=wh*wk xz=wh call fn(xz,d1,g2) wk=1.d0/wk if(wk)250,260,260 250 b=b+.5d0 260 b=b-(gw+g2)/2.d0+gw*g2 270 h2=wh*wh a2=wk*wk h4=h2/2.d0 if(h4.gt. 180.) h4=180. ex=dexp(-h4) w2=h4*ex ap=1.d0 s2=ap-ex sp=ap s1=0.d0 sn=s1 conex=dabs(con/wk) go to 290 280 sn=sp sp=sp+1.d0 s2=s2-w2 w2=w2*h4/sp ap=-ap*a2 290 cn=ap*s2/(sn+sp) s1=s1+cn if (dabs(cn)-conex) 300,300,280 300 t=(datan(wk)-wk*s1)/twopi 310 b=b+sgn*t 320 if (is) 330,350,350 330 if (ak) 340,350,340 340 wh=-ak wk=(ah/ak-r)/sqr gw=2.0d0*gk is=1 go to 210 350 if (b) 360,360,370 360 bivn=1.d-16 370 if(b-1.d0) 390,390,380 380 bivn=1.d0-1.d-16 390 bivn=b if(b.lt.1.d-16)bivn=1.d-16 401 return end c c---------------------------------------------------------------------- c c this function insures that all correlation coefficient guesses c fall between -1 and 1. c real*8 function corr(x) implicit real*8(a-h,o-z) c s = dexp(x) t = s/(1. + s) corr = (2.*t) - 1. return end c c---------------------------------------------------------------------- c subroutine datget implicit real*8(a-h,o-z) parameter(nfacs=1000,nrow=3504,nidx=6,nxt=50,ntype=4,nparis=14, $ nvars=122,norder=7,ncut=norder+1) character*8 arange(7) dimension nread(nidx),nxidx(nidx), $ nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars),nmt(nfacs) dimension xt(nvars),xreal(nfacs,nvars), $ xmean(nvars),xcov(nvars,nvars),xcut(nvars,ncut), $ xwrite(nvars),xtemp(nfacs),ytemp(nfacs) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut common/stand/ystand(nvars,203),xstand(nvars,203) c data arange/' 0% ',' 1-20 ',' 21-40 ',' 41-60 ', $ ' 61-80 ',' 81-99 ',' 100% '/ data nread/ 28, 30, 24, 40, 13, 9/ data nxidx/ 27, 27, 14, 37, 10, 7/ c c--------------------------------------------------------------------- c c read the array "nxkind" to describe each variable as: c "binary" (nxkind = 1), c "ordered discrete" (nxkind = 2), c "continuous" (nxkind = 3). c read(10,*) iestflg,iacrj,itype,nfac,nvar read(10,*) (nxkind(i),i=1,nvar) c c read an array, nxcut, which contains the number of c cut-off points for the discrete variables: c "binary" (nxcut = 3) c "ordered discrete" c (no. physicians, atwork: nxcut = 6 + 1) c (time since problems: nxcut = 4 + 1) c (no. vehicles: nxcut = 4 + 1) c "continuous" (nxcut = 0) c read(10,*) (nxcut(i),i=1,nvar) c c c--------------------------------------------------------------------- c 1 format(3(1x,i2),9(1x,f4.1)) c--------------------------------------------------------------------- c c initialize data arrays (nid, nonmis, xreal). c do 10 i = 1,nfacs do 8 j = 1,3 8 nid(i,j) = 0 do 10 j = 1,nvars nonmis(i,j) = 0 nxbnry(i,j) = -13 nxordd(i,j) = -13 10 xreal(i,j) = -13. c c--------------------------------------------------------------------- c c Steve has set up a monte carlo simulation to test "jamsimul.f". c run programs using simulated data ONLY if "imonte = 1" . c A new, second monte carlo test runs ONLY if "imonte = 13". c cmonte imonte = 0 cmonte imonte = 13 c imonte = 1 c if (imonte .eq. 1) goto 11 if (imonte .eq. 13) goto 11 if (imonte .eq. 0) goto 101 c 11 continue do 100 iii = 1,nfac c c initialize temporary array for writing raw data, "xwrite". c do 12 j = 1,nvars 12 xwrite(j) = -99. c c initialize temporary simulated data array, "xt". c do 13 j = 1,nxt 13 xt(j) = -13. c c read simulated data for this observation. c read(11,*) (xt(j),j=1,nvar) c c assign all simulated observations facility types, c parishes and numbers. c nid(iii,1) = 1 nid(iii,2) = 1 nid(iii,3) = iii c c code nonmissing value array, "nonmis", for simulated data. c also, code binary, ord. disc., real-valued data arrays. c do 468 j=1,nvar if (xt(j) .ge. -8.99) then nonmis(iii,j) = 1 if (nxkind(j) .eq. 1) nxbnry(iii,j) = xt(j) if (nxkind(j) .eq. 2) nxordd(iii,j) = xt(j) if (nxkind(j) .eq. 3) xreal(iii,j) = xt(j) endif c if (xt(j) .lt. -8.99) then nonmis(iii,j) = 0 if (nxkind(j) .eq. 1) nxbnry(iii,j) = -9 if (nxkind(j) .eq. 2) nxordd(iii,j) = -9 if (nxkind(j) .eq. 3) xreal(iii,j) = -9. endif 468 continue c c finished preparing simulated data for this observation. c 100 continue c c--------------------------------------------------------------------- c c finished preparing simulated data for all observations. c goto 201 c c--------------------------------------------------------------------- c c c All code for reading real Jamaican Data has been deleted. c 101 continue c c--------------------------------------------------------------------- c 201 continue c c--------------------------------------------------------------------- c c call subroutine to standardize continuously distributed vars. c if "istan = 1". c istan = 0 if (istan .eq. 1) then do 520 i = itype,itype do 510 j = 1,nvar if (nxkind(j) .eq. 3) then do 507 k = 1,nfac c c want to standardize by type of facility. c nmt(k) = -9 xtemp(k) = -9. if (nid(k,1) .eq. i) then nmt(k) = nonmis(k,j) xtemp(k) = xreal(k,j) endif 507 continue call standard(nmt,xtemp,j,ytemp,nfac) do 508 k = 1,nfac if (nid(k,1) .eq. i) then xreal(k,j) = ytemp(k) endif 508 continue endif 510 continue 520 continue endif c c------------------------------------------------------------------------- c c if iundo = 1, test subroutine undo. c iundo = 0 if (iundo .eq. 1) then do 620 i = 1,nvar if (nxkind(i) .eq. 3) then do 610 j = 1,nfac if (nonmis(j,i) .eq. 1) then xu = xreal(j,i) call undo(xu,i,xut) xreal(j,i) = xut endif 610 continue endif 620 continue endif c c--------------------------------------------------------------------- c return end c c--------------------------------------------------------------------- c c this function is used to compute the covariance between two c binary variables. c real*8 function func2(rt) implicit real*8(a-h,o-z) parameter(nfacs=1000,nvars=122,ntype=4,nparis=14, $ norder=7,ncut=norder+1) dimension nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xreal(nfacs,nvars),xmean(nvars), $ xcov(nvars,nvars), $ xcut(nvars,ncut) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut c data pi/3.1415927/ c r = corr(rt) c r = (2./pi)*(datan(rt)) c func2 = 0. do 200 k = 1,nfac c c if this facility is of the appropriate type and parish, then ... c c if ((nid(k,1) .eq. it) .and. (nid(k,2) .eq. ip)) then c if (nid(k,1) .eq. it) then c inmi = nonmis(k,ii) inmj = nonmis(k,ij) c c if both variables are observed for this facility, then ... c if ((inmi .eq. 1) .and. (inmj .eq. 1)) then c xmi = xmean(ii) xmj = xmean(ij) nxi = nxbnry(k,ii) nxj = nxbnry(k,ij) xi = dfloat(nxi) xj = dfloat(nxj) xmi = ((2.*xi) - 1.)*xmi xmj = ((2.*xj) - 1.)*xmj r = ((2.*xi) - 1.)*((2.*xj) - 1.)*r c c recall, function bivn returns Prob(x>a,y>b|r). c xmi = -xmi xmj = -xmj b = bivn(xmi,xmj,r) if (b .lt. 0.00000001) b = 0.00000001 bb = dlog(b) func2 = func2 + bb c endif c endif c 200 continue c return end c c---------------------------------------------------------------------- c c this function is used to compute the covariance between a c continuous variable and a binary variable. c real*8 function func3(rt) implicit real*8(a-h,o-z) parameter(nfacs=1000,nvars=122,ntype=4,nparis=14, $ norder=7,ncut=norder+1) real anormc,anormd,arg1,arg2 dimension nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xreal(nfacs,nvars),xmean(nvars), $ xcov(nvars,nvars), $ xcut(nvars,ncut) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut c data pi/3.1415927/ c r = corr(rt) c r = (2./pi)*(datan(rt)) c func3 = 0. do 200 k = 1,nfac c c if this facility is of the appropriate type and parish, then ... c c if ((nid(k,1) .eq. it) .and. (nid(k,2) .eq. ip)) then c if (nid(k,1) .eq. it) then c inmi = nonmis(k,ii) inmj = nonmis(k,ij) c c if both variables are observed for this facility, then ... c if ((inmi .eq. 1) .and. (inmj .eq. 1)) then c c if xi is the continuous variable and xj is binary. c if (nxki .eq. 3) then xmi = xmean(ii) xmj = xmean(ij) xi = xreal(k,ii) nxj = nxbnry(k,ij) xj = dfloat(nxj) xsigi = xcov(ii,ii)**(0.5) if (xsigi .lt. 0.0001) go to 200 xa1 = (xi - xmi)/(xsigi) xa2 = xmj rtt = r if (xj .gt. 0.5) then xa2 = -xmj rtt = -r endif arg1 = (xa2 - (rtt*xa1))/(dsqrt(1 - (rtt*rtt))) t1 = anormc(arg1) arg2 = xa1 t2 = anormd(arg2) endif c c if xi is the binary variable and xj is continuous. c if (nxkj .eq. 3) then xmi = xmean(ii) xmj = xmean(ij) nxi = nxbnry(k,ii) xi = dfloat(nxi) xj = xreal(k,ij) xsigj = xcov(ij,ij)**(0.5) if (xsigj .lt. 0.0001) go to 200 xa1 = (xj - xmj)/(xsigj) xa2 = xmi rtt = r if (xi .gt. 0.5) then xa2 = -xmi rtt = -r endif arg1 = (xa2 - (rtt*xa1))/(dsqrt(1 - (rtt*rtt))) t1 = anormc(arg1) arg2 = xa1 t2 = anormd(arg2) endif c b = t1*t2 if (b .lt. 0.00000001) b = 0.00000001 bb = dlog(b) func3 = func3 + bb c endif c endif c 200 continue c return end c c---------------------------------------------------------------------- c c this function is used to compute the covariance between two c ordered discrete variables. c real*8 function func4(rt) implicit real*8(a-h,o-z) parameter(nfacs=1000,nvars=122,ntype=4,nparis=14, $ norder=7,ncut=norder+1) dimension nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xreal(nfacs,nvars),xmean(nvars), $ xcov(nvars,nvars), $ xcut(nvars,ncut) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut c data pi/3.1415927/ c r = corr(rt) c r = (2./pi)*(datan(rt)) c func4 = 0. do 200 k = 1,nfac c c if this facility is of the appropriate type and parish, then ... c c if ((nid(k,1) .eq. it) .and. (nid(k,2) .eq. ip)) then c if (nid(k,1) .eq. it) then c inmi = nonmis(k,ii) inmj = nonmis(k,ij) c c if both variables are observed for this facility, then ... c if ((inmi .eq. 1) .and. (inmj .eq. 1)) then iti = nxordd(k,ii) itj = nxordd(k,ij) xciu = xcut(ii,iti+1) xcil = xcut(ii,iti) xcju = xcut(ij,itj+1) xcjl = xcut(ij,itj) c c recall, function bivn returns Prob(x>a,y>b|r). c xciu = -xciu xcil = -xcil xcju = -xcju xcjl = -xcjl b1 = bivn(xciu,xcju,r) b2 = bivn(xcil,xcju,r) b3 = bivn(xciu,xcjl,r) b4 = bivn(xcil,xcjl,r) b = b1 - b2 - b3 + b4 if (b .lt. 0.00000001) b = 0.00000001 bb = dlog(b) func4 = func4 + bb endif c endif c 200 continue c return end c c---------------------------------------------------------------------- c c this function is used to compute the covariance between a c continuous variable and an ordered discrete variable. c real*8 function func5(rt) implicit real*8(a-h,o-z) parameter(nfacs=1000,nvars=122,ntype=4,nparis=14, $ norder=7,ncut=norder+1) real anormc,anormd,arg dimension nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xreal(nfacs,nvars),xmean(nvars), $ xcov(nvars,nvars), $ xcut(nvars,ncut) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut c data pi/3.1415927/ c r = corr(rt) c r = (2./pi)*(datan(rt)) c func5 = 0. do 200 k = 1,nfac c c if this facility is of the appropriate type and parish, then ... c c if ((nid(k,1) .eq. it) .and. (nid(k,2) .eq. ip)) then c if (nid(k,1) .eq. it) then c inmi = nonmis(k,ii) inmj = nonmis(k,ij) c c if both variables are observed for this facility, then ... c if ((inmi .eq. 1) .and. (inmj .eq. 1)) then c c if xi is ordered discrete and xj continuous. c if (nxkj .eq. 3) then xj = xreal(k,ij) xmj = xmean(ij) xsigj = dsqrt(xcov(ij,ij)) if (xsigj .lt. 0.00001) go to 200 xa1 = (xj - xmj)/(xsigj) iti = nxordd(k,ii) xcut1 = xcut(ii,iti+1) xa2 = xcut1 arg = (xa2 - (r*xa1))/(dsqrt(1. - (r*r))) t1 = anormc(arg) arg = xa1 t2 = anormd(arg) b1 = t1*t2 xcut2 = xcut(ii,iti) xa2 = xcut2 arg = (xa2 - (r*xa1))/(dsqrt(1. - (r*r))) t1 = anormc(arg) b2 = t1*t2 endif c c if xi is continuous and xj ordered discrete. c if (nxki .eq. 3) then xi = xreal(k,ii) xmi = xmean(ii) xsigi = dsqrt(xcov(ii,ii)) if (xsigi .lt. 0.00001) go to 200 xa1 = (xi - xmi)/(xsigi) itj = nxordd(k,ij) xcut1 = xcut(ij,itj+1) xa2 = xcut1 arg = (xa2 - (r*xa1))/(dsqrt(1. - (r*r))) t1 = anormc(arg) arg = xa1 t2 = anormd(arg) b1 = t1*t2 xcut2 = xcut(ij,itj) xa2 = xcut2 arg = (xa2 - (r*xa1))/(dsqrt(1. - (r*r))) t1 = anormc(arg) b2 = t1*t2 endif c b = b1 - b2 if (b .lt. 0.00000001) b = 0.00000001 bb = dlog(b) func5 = func5 + bb c endif c endif c 200 continue c return end c c---------------------------------------------------------------------- c c this function is used to compute the covariance between a c binary variable and an ordered discrete variable. c real*8 function func6(rt) implicit real*8(a-h,o-z) parameter(nfacs=1000,nvars=122,ntype=4,nparis=14, $ norder=7,ncut=norder+1) dimension nid(nfacs,3),nonmis(nfacs,nvars), $ nxbnry(nfacs,nvars),nxordd(nfacs,nvars), $ nxkind(nvars),nxcut(nvars) dimension xreal(nfacs,nvars),xmean(nvars), $ xcov(nvars,nvars), $ xcut(nvars,ncut) c common/ishare/it,ip,ii,ij,nxki,nxkj common/nshare/nid,nonmis,nxbnry,nxordd,nxkind,nxcut, $ iestflg,iacrj,itype,nvar,nfac common/xshare/xreal,xmean,xcov,xcut c data pi/3.1415927/ c r = corr(rt) c r = (2./pi)*(datan(rt)) c func6 = 0. do 200 k = 1,nfac c c if this facility is of the appropriate type and parish, then ... c c if ((nid(k,1) .eq. it) .and. (nid(k,2) .eq. ip)) then c if (nid(k,1) .eq. it) then c inmi = nonmis(k,ii) inmj = nonmis(k,ij) c c if both variables are observed for this facility, then ... c if ((inmi .eq. 1) .and. (inmj .eq. 1)) then c c if xi is binary and xj is ordered discrete. c if (nxki .eq. 1) then iti = nxbnry(k,ii) + 1 itj = nxordd(k,ij) endif c c if xi is ordered discrete and xj is binary. c if (nxki .eq. 2) then iti = nxordd(k,ii) itj = nxbnry(k,ij) + 1 endif c xciu = xcut(ii,iti+1) xcil = xcut(ii,iti) xcju = xcut(ij,itj+1) xcjl = xcut(ij,itj) c c recall, function bivn returns Prob(x>a,y>b|r). c xciu = -xciu xcil = -xcil xcju = -xcju xcjl = -xcjl b1 = bivn(xciu,xcju,r) b2 = bivn(xcil,xcju,r) b3 = bivn(xciu,xcjl,r) b4 = bivn(xcil,xcjl,r) b = b1 - b2 - b3 + b4 if (b .lt. 0.00000001) b = 0.00000001 bb = dlog(b) func6 = func6 + bb c endif c endif c 200 continue c return end c c---------------------------------------------------------------------- c function golden(ax,bx,cx,f,tol,xmin) implicit real*8(a-h,o-z) c c This code is from Numerical Recipes, p. 282. c given a function f, and given a bracketing triplet of c abscissas ax, bx, cx (such that bx is between ax and c cx, and f(bx) is less than both f(ax) and f(cx), this c routine performs a golden section search for the mi - c nimum, isolating it to a fractional precision of about c tol. the abscissa of the minimum is returned as xmin, c and the minimum function value is returned as golden, c the returned fucntion value. c c golden ratios. c parameter (r=.61803399,c=1.-r) c c at any given time we will keep track of four points, x0, c x1, x2, x3. c x0=ax x3=cx c c make x0 to x1 the small segment c if (abs(cx-bx).gt.abs(bx-ax)) then x1=bx c c and fill in the new point to be tried. c x2=bx+c*(cx-bx) else x2=bx x1=bx-c*(bx-ax) endif c c the initial function evaluations. note that we never need to c evaluate the function at the original endpoints. c f1=-f(x1) f2=-f(x2) c c do-while loop: we keep returning here. c 1 if (abs(x3-x0).gt.tol*(abs(x1)+abs(x2))) then c c one possible outcome, c if (f2.lt.f1) then c c its housekeeping, c x0=x1 x1=x2 x2=r*x1+c*x3 f1=f2 c c and a new function evaluation. c f2=-f(x2) c c the other outcome, c else x3=x2 x2=x1 x1=r*x2+c*x0 f2=f1 c c and its new function evaluation c f1=-f(x1) endif c c back to see if we are done. c goto 1 endif c c we are done. output the best of the two current values. c if (f1.lt.f2) then golden=f1 xmin=x1 else golden=f2 xmin=x2 endif c return end c c---------------------------------------------------------------------- c subroutine mnbrak (ax, bx, cx, fa, fb, fc, func) implicit real*8(a-h,o-z) c c this routine comes from Numerical Recipes, p. 281. c given a function func, and given distinct initial points c ax and bx, this routine searches in the downhill direc - c tion (defined by the function as evaluated at the initi- c al points returns new points ax, bx, cx which bracket a c minimum of the function. Also returned are the function c values at the three points, fa, fb, and fc. c parameter (gold=1.618034, glimit=100., tiny=1.e-20) c c the first parameter is the default ratio by which succe- c sive intervals are magnified; the second is the maximum c magnification allowed for a parabolic fit-step. c fa=-func(ax) fb=-func(bx) c c switch roles of a and b so that we can go downhill in the di- c rection from a to b. c if (fb.gt.fa) then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif c c first guess for c. c cx=bx+gold*(bx-ax) fc=-func(cx) c c "do while": keep returning here until we bracket. c 1 if (fb.ge.fc) then c c compute u by parabolic extrapolation from a,b,c. tiny is used c to prevent any possible division by zero. c r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/ $ (2.*sign(max(abs(q-r),tiny),q-r)) c c we won't go farther than this. now to test various possibi- c lities: c ulim=bx+glimit*(cx-bx) c c parabolic u is between b and c: try it. c if ((bx-u)*(u-cx).gt.0.) then fu=-func(u) c c got a minimum between b and c. c if (fu.lt.fc) then ax=bx fa=fb bx=u fb=fu return c c got a minimum between a and u. c else if (fu.gt.fb) then cx=u fc=fu return endif c c parabolic fit was no use. use default magnification. c u=cx+gold*(cx-bx) fu=-func(u) c c parabolic fit is between c and its allowed limit. c else if ((cx-u)*(u-ulim).gt.0.) then fu=-func(u) if (fu.lt.fc) then bx=cx cx=u u=cx+gold*(cx-bx) fb=fc fc=fu fu=-func(u) endif c c limit parabolic u to maximum allowed value. c else if ((u-ulim)*(ulim-cx).ge.0.) then u=ulim fu=-func(u) c c reject parabolic u, use default magnification. c else u=cx+gold*(cx-bx) fu=-func(u) endif c c eliminate oldest point and continue. c ax=bx bx=cx cx=u fa=fb fb=fc fc=fu go to 1 endif c return end c c---------------------------------------------------------------------- c c imsl routine name - mdnris c c----------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - november 1, 1984 c c purpose - inverse standard normal (gaussian) c probability distribution function c c usage - call mdnris (p,y,ier) c c arguments p - input value in the exclusive range (0.0,1.0) c y - output value of the inverse normal (0,1) c probability distribution function c ier - error parameter (output) c terminal error c ier = 129 indicates p lies outside the legal c range. plus or minus machine infinity is c given as the result (sign is the sign of c the function value of the nearest legal c argument). c c precision/hardware - single/all c c reqd. imsl routines - merfi c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1982 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c----------------------------------------------------------------------- c subroutine mdnris (p,y,ier) c specifications for arguments real p,y integer ier c specifications for local variables real eps,g0,g1,g2,g3,h0,h1,h2,a,w,wi,sn,sd real sigma,sqrt2,x,xinf data xinf/.12650140831e+323/ data sqrt2/1.4142135623731/ data eps/.710543e-14/ data g0/.18511591e-3/,g1/-.20281520e-2/ data g2/-.14983844/,g3/.10786386e-1/ data h0/.99529751e-1/,h1/.52117329/ data h2/-.68883009e-1/ c first executable statement ier = 0 if (p .gt. 0.0 .and. p .lt. 1.0) go to 5 ier = 129 sigma = sign(1.0,p) y = sigma * xinf go to 9000 5 if(p.le.eps) go to 10 x = 1.0 -(p + p) call merfi (x,y,ier) y = -sqrt2 * y go to 9005 c p too small, compute y directly 10 a = p+p w = (-log(a+(a-a*a)))**.5 c use a rational function in 1./w wi = 1./w sn = ((g3*wi+g2)*wi+g1)*wi sd = ((wi+h2)*wi+h1)*wi+h0 y = w + w*(g0+sn/sd) y = -y*sqrt2 go to 9005 9000 continue 9005 return end c c---------------------------------------------------------------------- c c imsl routine name - merfi c c---------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - january 1, 1978 c c purpose - inverse error function c c usage - call merfi (p,y,ier) c c arguments p - input value in the exclusive range (-1.0,1.0) c y - output value of the inverse error function c ier - error parameter (output) c terminal error c ier = 129 indicates p lies outside the lega c range. plus or minus machine infinity is c given as the result (sign is the sign of c the function value of the nearest legal c argument). c c precision/hardware - single/all c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1978 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c---------------------------------------------------------------------- c subroutine merfi (p,y,ier) implicit real(a-h,o-z) c specifications for arguments real p,y integer ier c specifications for local variables real a(65),h1,h2,h3,h4,xinf,x,sigma,z,w,x3,x4,x5, * x6,b integer n,ipp,l,lb2 data a(1),a(2),a(3),a(4),a(5),a(6),a(7),a(8),a(9), * a(10),a(11),a(12),a(13),a(14),a(15),a(16), * a(17),a(18),a(19),a(20),a(21),a(22),a(23) * /.99288537662,.12046751614, * .16078199342e-01,.26867044372e-02, * .49963473024e-03,.98898218599e-04, * .20391812764e-04,.4327271618e-05, * .938081413e-06,.206734721e-06, * .46159699e-07,.10416680e-07, * .2371501e-08,.543928e-09, * .125549e-09,.29138e-10, * .6795e-11,.1591e-11, * .374e-12,.88e-13, * .21e-13,.5e-14, * .1e-14/ data a(24),a(25),a(26),a(27),a(28),a(29),a(30), * a(31),a(32),a(33),a(34),a(35),a(36),a(37), * a(38),a(39),a(40) * /.91215880342e00,-.16266281868e-01, * .43355647295e-03,.21443857007e-03, * .2625751076e-05,-.302109105e-05, * -.12406062e-07,.62406609e-07, * -.540125e-09,-.142328e-08, * .34384e-10,.33585e-10, * -.1458e-11,-.81e-12, * .53e-13,.2e-13, * -.2e-14/ data a(41),a(42),a(43),a(44),a(45),a(46),a(47), * a(48),a(49),a(50),a(51),a(52),a(53),a(54), * a(55),a(56),a(57),a(58),a(59),a(60),a(61), * a(62),a(63),a(64),a(65) * /.95667970902,-.023107004309, * -.43742360975e-02,-.57650342265e-03, * -.10961022307e-04,.25108547025e-04, * .10562336068e-04,.275441233e-05, * .432484498e-06,-.20530337e-07, * -.43891537e-07,-.1768401e-07, * -.3991289e-08,-.186932e-09, * .272923e-09,.132817e-09, * .31834e-10,.1670e-11, * -.2036e-11,-.965e-12, * -.22e-12,-.1e-13, * .13e-13,.6e-14, * .1e-14/ data h1,h2,h3,h4/-1.5488130424, * 2.5654901231,-.55945763133, * 2.2879157163/ data xinf/.12650140831e+323/ c first executable statement x = p ier = 0 sigma = sign(1.,x) if (.not.(x.gt.-1..and.x.lt.1.)) go to 35 z = abs(x) if(z.gt..8) go to 20 w = z*z/.32-1. n = 22 ipp = 1 l = 1 5 lb2 = 1 x3 = 1. x4 = w x6 = a(ipp) 10 x6 = x6 + a(ipp+lb2) * x4 x5 = x4 * w * 2.-x3 x3 = x4 x4 = x5 lb2 = lb2 + 1 if (lb2 .le. n) go to 10 go to (15,30),l 15 y = z * x6 * sigma go to 9005 20 b = (-log(1.-z*z))**.5 if (z .gt. .9975) go to 25 w = h1*b+h2 ipp = 24 l = 2 n = 16 go to 5 25 w = h3 * b + h4 ipp = 41 n = 24 l = 2 go to 5 30 y = b * x6 * sigma go to 9005 35 y = sigma*xinf ier = 129 9000 continue 9005 return end c c---------------------------------------------------------------------- c subroutine standard(nmis,x,istand,y,nobs) c c this subroutine translates x into y so that y is normal c (with mean zero, unit variance) c and stores the translation matrix in standmat. c parameter(nobss=1000,nxs=122) implicit real*8(a-h,o-z) real frat,yrat dimension ixc(nobss),nmis(nobss),x(nobss),y(nobss), & xmesh(203),yco(203), & ymesh(203) common/stand/ystand(nxs,203),xstand(nxs,203) xmin=1000000. xmax=-xmin nobsa=nobs do 2 i=1,nobs ixc(i)=0 c if(x(i).ge.0.)then if (nmis(i) .eq. 1) then if(xmin.gt.x(i))xmin=x(i) if(xmax.lt.x(i))xmax=x(i) endif c if(x(i).lt.0.)nobsa=nobsa-1 if (nmis(i) .eq. 0) nobsa=nobsa-1 2 continue xmesh(2)=xmin del=(xmax-xmin)/200 do 3 i=3,202 3 xmesh(i)=xmesh(i-1)+del xmesh(1)=xmesh(2)-(del*50.) xmesh(203)=xmesh(202)+(del*50.) do 6 i=1,203 6 yco(i)=0. do 4 i=1,nobs if((x(i).gt.(xmin-0.001)).and.(x(i).le.xmin))then yco(2)=yco(2)+1. ixc(i)=2 endif if(x(i).gt.xmin)then do 5 j=2,202 if((x(i).gt.xmesh(j)).and.(x(i).le.xmesh(j+1)))then yco(j)=yco(j)+1. ixc(i)=j goto 4 endif 5 continue endif 4 continue ymesh(1)=-20. ymesh(203)=20. ysum=0. nobsb=nobsa+1 do 7 i=2,202 if(yco(i).eq.0.)ymesh(i)=ymesh(i-1) if(yco(i).gt.0.)then ysum=ysum+yco(i) yrat=ysum/nobsb if (yrat .lt. 0.0000001) yrat = 0.0000001 if (yrat .gt. 0.9999999) yrat = 0.9999999 call mdnris(yrat,frat,ier) ymesh(i)=frat+10. endif 7 continue do 9 i=1,203 xstand(istand,i)=xmesh(i) 9 ystand(istand,i)=ymesh(i) do 8 i=1,nobs c if(x(i).ge.0.)y(i)=ymesh(ixc(i)) c if(x(i).lt.0.)y(i)=x(i) if (nmis(i) .eq. 1) y(i) = ymesh(ixc(i)) if (nmis(i) .eq. 0) y(i) = x(i) 8 continue return end c c---------------------------------------------------------------------- c subroutine undo(y,istand,x) c c this subroutine undoes standard. c parameter(nobss=1000,nxs=122) implicit real*8(a-h,o-z) common/stand/ystand(nxs,203),xstand(nxs,203) if(y.lt.0.)x=y if(y.ge.0.)then do 3 i=1,202 if((y.gt.ystand(istand,i)).and.(y.le. & ystand(istand,(i+1))))then ybot=ystand(istand,i) ytop=ystand(istand,(i+1)) xbot=xstand(istand,i) xtop=xstand(istand,(i+1)) if(xtop.eq.xbot)then x=xtop return endif if(xtop.gt.xbot)then ydif=ytop-ybot xdif=xtop-xbot yfrac=(y-ybot)/ydif x=xbot+(yfrac*xdif) cx x=xbot+(0.5*xdif) endif endif 3 continue endif return end c c--------------------------------------------------------------------- c c this subroutine checks eigenvalues of the covariance matrix. c c if estimating monte carlo data #1, "nrowe" = 15 (eigchk1). c if estimating monte carlo data #2, "nrowe" = 122 (eigchk3). c if estimating moments for parish #1, "nrowe" = 89 (eigchk1); c if estimating moments for parish #2, "nrowe" = 103 (eigchk2). c c subroutine eigchk1(itype,covt,xcovn) implicit real*8(a-h,o-z) parameter(nrowes=15,nrowe2s=(nrowes+1)*nrowes/2) dimension cov(nrowes,nrowes),covc(nrowe2s), & covn(nrowes,nrowes), & covt(122,122),eigval(nrowes),eigvec(nrowes,nrowes), & ichk(nrowes),variance(nrowes),work(nrowes),xcovn(122,122) data alammin/.2/ c if (itype .eq. 1) print *, 'got to subroutine eigchk1' if (itype .eq. 2) then print *,'got to wrong subroutine -- eigchk1' stop endif c nrowe=0 do 10 i=1,15 if(covt(i,i).le.0.)then print *, i, covt(i,i) endif c if(covt(i,i).gt.0.)then nrowe=nrowe+1 ichk(nrowe)=i endif if(covt(i,i).le.0.)then do 12 j=1,15 if((covt(i,j).ne.0.).or.(covt(j,i).ne.0.))then write(6,103)i,i,j,covt(i,i),covt(i,j),covt(j,i) 103 format(1x,'problem with var ',i3,' and cov ',i3,2x, & i3,/,3x,3(1x,g15.8)) stop endif 12 continue endif 10 continue c c shrink "covt" to "cov" by deleting variables with zero variance. c do 11 i=1,nrowe ichkt=ichk(i) do 11 j=1,nrowe jchkt=ichk(j) 11 cov(i,j)=covt(ichkt,jchkt) c c keep variances and copy "cov" to "covc" for subroutine call. ij=0 do 2 i=1,nrowe variance(i)=cov(i,i) charry: "cy" next 6 lines for real runs of program. if(abs(variance(i)).lt.alammin)then write(6,101)i,variance(i),alammin 101 format(1x,'var for ',i3,' is ',g15.8, & ' which is less than ',f8.5) stop endif do 2 j=1,i ij=ij+1 2 covc(ij)=cov(i,j) c c c subroutine eigchks returns eigenvalues and eigenvectors of covc. c call eigrs(covc,nrowe,1,eigval,eigvec,15,work,ier) c c charry: "cy" next 2 lines for real runs of program. write(6,100)(eigval(i),i=1,nrowe) 100 format(1x,'eigenvalues',/,25(5(1x,g15.8),/)) c c check that product of eigenvalues and eigenvector equals c the original covariance matrix. c do 5 i=1,nrowe do 5 j=1,nrowe covn(i,j)=0. do 6 k=1,nrowe 6 covn(i,j)=covn(i,j)+(eigvec(i,k)*eigvec(j,k)*eigval(k)) if(abs(covn(i,j)-cov(i,j)).gt..01)then write(6,105)i,j,covn(i,j),cov(i,j) 105 format(1x,'problem with ',i3,2x,i3,': ',2(1x,g15.8)) endif 5 continue c c replace very small and negative eigenvalues with positive value. c (this ensures positive definiteness of final covariance matrix.) c do 3 i=1,nrowe if(eigval(i).lt.alammin)eigval(i)=alammin if(eigval(i).gt.alammin)goto 4 3 continue 4 continue c c positive definite covariance matrix results from multiplying c new eigenvalues (all positive) with original eigenvectors. c do 15 i=1,nrowe do 15 j=1,nrowe covn(i,j)=0. do 16 k=1,nrowe 16 covn(i,j)=covn(i,j)+(eigvec(i,k)*eigvec(j,k)*eigval(k)) 15 continue c c divide all rows and columns by standard deviation to return to c proper normalization for variances and correlations (unit c variance for binary and ordered discrete random variables). c do 7 i=1,nrowe rat=covn(i,i)/variance(i) rat=rat**.5 do 8 j=1,nrowe covn(i,j)=covn(i,j)/rat 8 covn(j,i)=covn(j,i)/rat 7 continue c charry: "cy" next 10 lines for real runs of program. c ij=0 do 9 i=1,nrowe do 9 j=1,i ij=ij+1 9 covc(ij)=covn(i,j) call eigrs(covc,nrowe,0,eigval,eigvec,15,work,ier) write(6,100)(eigval(i),i=1,nrowe) cy do 25 i=1,103 cy do 26 j=1,i cy 26 covt(i,j)=covn(i,j)-cov(i,j) cy 25 write(6,110)i,(covt(i,j),j=1,i) cy 110 format(1x,i3,5(2x,f8.4),17(/,4x,5(2x,f8.4))) c do 30 i = 1,15 do 30 j = 1,15 if (i.eq.j) xcovn(i,j) = 0.1000 if (i.ne.j) xcovn(i,j) = 0.0000 30 continue cmike 30 xcovn(i,j) = 0.0000 do 32 i = 1,nrowe ichkt = ichk(i) do 32 j = 1,nrowe jchkt = ichk(j) 32 xcovn(ichkt,jchkt) = covn(i,j) c return end c c---------------------------------------------------------------------- c c c imsl routine name - ehobks c c---------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - june 1, 1982 c c purpose - nucleus called only by imsl routine eigrs c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c reqd. imsl routines - none required c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1982 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c---------------------------------------------------------------------- c subroutine ehobks (a,n,m1,m2,z,iz) implicit real*8(a-h,o-z) c dimension a(1),z(iz,1) c first executable statement if (n .eq. 1) go to 30 do 25 i=2,n l = i-1 ia = (i*l)/2 h = a(ia+i) if (h.eq.0.) go to 25 c derives eigenvectors m1 to m2 of c the original matrix from eigenvector c m1 to m2 of the symmetric c tridiagonal matrix do 20 j = m1,m2 s = 0.0 do 10 k = 1,l s = s+a(ia+k)*z(k,j) 10 continue s = s/h do 15 k=1,l z(k,j) = z(k,j)-s*a(ia+k) 15 continue 20 continue 25 continue 30 return end c c c imsl routine name - ehouss c c---------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - november 1, 1984 c c purpose - nucleus called only by imsl routine eigrs c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c reqd. imsl routines - none required c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1982 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c---------------------------------------------------------------------- c subroutine ehouss (a,n,d,e,e2) implicit real*8(a-h,o-z) c dimension a(1),d(n),e(n),e2(n) real*8 a,d,e,e2,zero,h,scale,f,g,hh data zero/0.0/ c first executable statement np1 = n+1 nn = (n*np1)/2-1 nbeg = nn+1-n do 70 ii = 1,n i = np1-ii l = i-1 h = zero scale = zero if (l .lt. 1) go to 10 c scale row (algol tol then not needed nk = nn do 5 k = 1,l scale = scale+abs(a(nk)) nk = nk-1 5 continue if (scale .ne. zero) go to 15 10 e(i) = zero e2(i) = zero go to 65 15 nk = nn do 20 k = 1,l a(nk) = a(nk)/scale h = h+a(nk)*a(nk) nk = nk-1 20 continue e2(i) = scale*scale*h f = a(nn) g = -sign(sqrt(h),f) e(i) = scale*g h = h-f*g a(nn) = f-g if (l .eq. 1) go to 55 f = zero jk1 = 1 do 40 j = 1,l g = zero ik = nbeg+1 jk = jk1 c form element of a*u do 25 k = 1,j g = g+a(jk)*a(ik) jk = jk+1 ik = ik+1 25 continue jp1 = j+1 if (l .lt. jp1) go to 35 jk = jk+j-1 do 30 k = jp1,l g = g+a(jk)*a(ik) jk = jk+k ik = ik+1 30 continue c form element of p 35 e(j) = g/h f = f+e(j)*a(nbeg+j) jk1 = jk1+j 40 continue hh = f/(h+h) c form reduced a jk = 1 do 50 j = 1,l f = a(nbeg+j) g = e(j)-hh*f e(j) = g do 45 k = 1,j a(jk) = a(jk)-f*e(k)-g*a(nbeg+k) jk = jk+1 45 continue 50 continue 55 do 60 k = 1,l a(nbeg+k) = scale*a(nbeg+k) 60 continue 65 d(i) = a(nbeg+i) a(nbeg+i) = h*scale*scale nbeg = nbeg-i+1 nn = nn-i 70 continue return end c c c imsl routine name - eigrs c c---------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - june 1, 1980 c c purpose - eigenvalues and (optionally) eigenvectors of c a real symmetric matrix c c usage - call eigrs (a,n,jobn,d,z,iz,wk,ier) c c arguments a - input real symmetric matrix of order n, c whose eigenvalues and eigenvectors c are to be computed. input a is c destroyed if ijob is equal to 0 or 1. c n - input order of the matrix a. c jobn - input option parameter. if jobn.ge.10 c a is assumed to be in full storage mode c (in this case, a must be dimensioned exactly c n by n in the calling program). c if jobn.lt.10 then a is assumed to be in c symmetric storage mode. define c ijob=mod(jobn,10). then when c ijob = 0, compute eigenvalues only c ijob = 1, compute eigenvalues and eigen- c vectors. c ijob = 2, compute eigenvalues, eigenvectors c and performance index. c ijob = 3, compute performance index only. c if the performance index is computed, it is c returned in wk(1). the routines have c performed (well, satisfactorily, poorly) if c wk(1) is (less than 1, between 1 and 100, c greater than 100). c d - output vector of length n, c containing the eigenvalues of a. c z - output n by n matrix containing c the eigenvectors of a. c the eigenvector in column j of z corres- c ponds to the eigenvalue d(j). c if ijob = 0, z is not used. c iz - input row dimension of matrix z exactly as c specified in the dimension statement in the c calling program. c wk - work area, the length of wk depends c on the value of ijob, when c ijob = 0, the length of wk is at least n. c ijob = 1, the length of wk is at least n. c ijob = 2, the length of wk is at least c n(n+1)/2+n. c ijob = 3, the length of wk is at least 1. c ier - error parameter (output) c terminal error c ier = 128+j, indicates that eqrt2s failed c to converge on eigenvalue j. eigenvalues c and eigenvectors 1,...,j-1 have been c computed correctly, but the eigenvalues c are unordered. the performance index c is set to 1000.0 c warning error (with fix) c in the following, ijob = mod(jobn,10). c ier = 66, indicates ijob is less than 0 or c ijob is greater than 3. ijob set to 1. c ier = 67, indicates ijob is not equal to c zero, and iz is less than the order of c matrix a. ijob is set to zero. c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c reqd. imsl routines - ehobks,ehouss,eqrt2s c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1980 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c---------------------------------------------------------------------- c subroutine eigrs (a,n,jobn,d,z,iz,wk,ier) implicit real*8(a-h,o-z) c specifications for arguments integer n,jobn,iz,ier real*8 a(1),d(1),wk(1),z(iz,1) c specifications for local variables integer ijob,ir,jr,ij,ji,np1 integer jer,na,nd,iiz,ibeg,il,kk,lk,i,j,k,l real*8 anorm,asum,pi,sumz,sumr,an,s,ten,rdelp,zero, 1 one,thous data rdelp/.710543e-14/ data zero,one/0.0,1.0/,ten/10.0/,thous/1000.0/ c initialize error parameters c first executable statement ier = 0 jer = 0 if (jobn.lt.10) go to 15 c convert to symmetric storage mode k = 1 ji = n-1 ij = 1 do 10 j=1,n do 5 i=1,j a(k) = a(ij) ij = ij+1 k = k+1 5 continue ij = ij + ji ji = ji - 1 10 continue 15 ijob = mod(jobn,10) if (ijob.ge.0.and.ijob.le.3) go to 20 c warning error - ijob is not in the c range ier = 66 ijob = 1 go to 25 20 if (ijob.eq.0) go to 35 25 if (iz.ge.n) go to 30 c warning error - iz is less than n c eigenvectors can not be computed, c ijob set to zero ier = 67 ijob = 0 30 if (ijob.eq.3) go to 75 35 na = (n*(n+1))/2 if (ijob.ne.2) go to 45 do 40 i=1,na wk(i) = a(i) 40 continue c save input a if ijob = 2 45 nd = 1 if (ijob.eq.2) nd = na+1 c reduce a to symmetric tridiagonal c form call ehouss (a,n,d,wk(nd),wk(nd)) iiz = 1 if (ijob.eq.0) go to 60 iiz = iz c set z to the identity matrix do 55 i=1,n do 50 j=1,n z(i,j) = zero 50 continue z(i,i) = one 55 continue c compute eigenvalues and eigenvectors 60 call eqrt2s (d,wk(nd),n,z,iiz,jer) if (ijob.eq.0) go to 9000 if (jer.gt.128) go to 65 c back transform eigenvectors call ehobks (a,n,1,n,z,iz) 65 if (ijob.le.1) go to 9000 c move input matrix back to a do 70 i=1,na a(i) = wk(i) 70 continue wk(1) = thous if (jer.ne.0) go to 9000 c compute 1 - norm of a 75 anorm = zero ibeg = 1 do 85 i=1,n asum = zero il = ibeg kk = 1 do 80 l=1,n asum = asum+abs(a(il)) if (l.ge.i) kk = l il = il+kk 80 continue if(anorm.lt.asum)anorm=asum c anorm = amax1(anorm,asum) ibeg = ibeg+i 85 continue if (anorm.eq.zero) anorm = one c compute performance index pi = zero do 100 i=1,n ibeg = 1 s = zero sumz = zero do 95 l=1,n lk = ibeg kk = 1 sumz = sumz+abs(z(l,i)) sumr = -d(i)*z(l,i) do 90 k=1,n sumr = sumr+a(lk)*z(k,i) if (k.ge.l) kk = k lk = lk+kk 90 continue s = s+abs(sumr) ibeg = ibeg+l 95 continue if (sumz.eq.zero) go to 100 ssumz=s/sumz if(pi.lt.ssumz)pi=ssumz c pi = amax1(pi,s/sumz) 100 continue an = n pi = pi/(anorm*ten*an*rdelp) wk(1) = pi if (jobn.lt.10) go to 9000 c convert back to full storage mode np1 = n+1 ij = (n-1)*np1 + 2 k = (n*(np1))/2 do 110 jr=1,n j = np1-jr do 105 ir=1,j ij = ij-1 a(ij) = a(k) k = k-1 105 continue ij = ij-jr 110 continue ji = 0 k = n-1 do 120 i=1,n ij = i-n do 115 j=1,i ij = ij+n ji = ji+1 a(ij) = a(ji) 115 continue ji = ji + k k = k-1 120 continue 9000 continue if (jer.eq.0) go to 9005 ier = jer 9005 return end c c c imsl routine name - eqrt2s c c---------------------------------------------------------------------- c c computer - cdcft5/single c c latest revision - november 1, 1984 c c purpose - eigenvalues and (optionally) eigenvectors of c a symmetric tridiagonal matrix using the c ql method. c c usage - call eqrt2s (d,e,n,z,iz,ier) c c arguments d - on input, the vector d of length n contains c the diagonal elements of the symmetric c tridiagonal matrix t. c on output, d contains the eigenvalues of c t in ascending order. c e - on input, the vector e of length n contains c the sub-diagonal elements of t in position c 2,...,n. on output, e is destroyed. c n - order of tridiagonal matrix t.(input) c z - on input, z contains the identity matrix of c order n. c on output, z contains the eigenvectors c of t. the eigenvector in column j of z c corresponds to the eigenvalue d(j). c iz - input row dimension of matrix z exactly as c specified in the dimension statement in the c calling program. if iz is less than n, the c eigenvectors are not computed. in this case c z is not used. c ier - error parameter c terminal error c ier = 128+j, indicates that eqrt2s failed c to converge on eigenvalue j. eigenvalues c and eigenvectors 1,...,j-1 have been c computed correctly, but the eigenvalues c are unordered. c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1978 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c---------------------------------------------------------------------- c subroutine eqrt2s (d,e,n,z,iz,ier) implicit real*8(a-h,o-z) c dimension d(1),e(1),z(iz,1) data rdelp/.710543e-14/ data zero,one/0.0,1.0/ c move the last n-1 elements c of e into the first n-1 locations c first executable statement ier = 0 if (n .eq. 1) go to 9005 do 5 i=2,n e(i-1) = e(i) 5 continue e(n) = zero b = zero f = zero do 60 l=1,n j = 0 h = rdelp*(abs(d(l))+abs(e(l))) if (b.lt.h) b = h c look for small sub-diagonal element do 10 m=l,n k=m if (abs(e(k)) .le. b) go to 15 10 continue 15 m = k if (m.eq.l) go to 55 20 if (j .eq. 30) go to 85 j = j+1 l1 = l+1 g = d(l) p = (d(l1)-g)/(e(l)+e(l)) r = abs(p) if (rdelp*abs(p) .lt. 1.0) r = sqrt(p*p+one) d(l) = e(l)/(p+sign(r,p)) h = g-d(l) do 25 i = l1,n d(i) = d(i)-h 25 continue f = f+h c ql transformation p = d(m) c = one s = zero mm1 = m-1 mm1pl = mm1+l if (l.gt.mm1) go to 50 do 45 ii=l,mm1 i = mm1pl-ii g = c*e(i) h = c*p if (abs(p).lt.abs(e(i))) go to 30 c = e(i)/p r = sqrt(c*c+one) e(i+1) = s*p*r s = c/r c = one/r go to 35 30 c = p/e(i) r = sqrt(c*c+one) e(i+1) = s*e(i)*r s = one/r c = c*s 35 p = c*d(i)-s*g d(i+1) = h+s*(c*g+s*d(i)) if (iz .lt. n) go to 45 c form vector do 40 k=1,n h = z(k,i+1) z(k,i+1) = s*z(k,i)+c*h z(k,i) = c*z(k,i)-s*h 40 continue 45 continue 50 e(l) = s*p d(l) = c*p if ( abs(e(l)) .gt.b) go to 20 55 d(l) = d(l) + f 60 continue c order eigenvalues and eigenvectors do 80 i=1,n k = i p = d(i) ip1 = i+1 if (ip1.gt.n) go to 70 do 65 j=ip1,n if (d(j) .ge. p) go to 65 k = j p = d(j) 65 continue 70 if (k.eq.i) go to 80 d(k) = d(i) d(i) = p if (iz .lt. n) go to 80 do 75 j = 1,n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 75 continue 80 continue go to 9005 85 ier = 128+l 9000 continue 9005 return end c c--------------------------------------------------------------------- c