! TLAYOFF.f90 ! Fortran code for: ! Temporary Layoffs and Split Population Models ! K.G. Mavromaras and C.D.Orme program layoff ! LAST MODIFIED ON 19th November 2002 ! F90 Program to obtain MLE for SPLIT or SINGLE POPULATION MODEL ! with application to Temporary Layoffs ! ! Assuming: (a) Weibull hazard for temporary layoffs, or ! Almon flexible hazard specification ! (b) Probit spec for probability of eventual recall ! ! Log-lik contributions: ! l(i)=(1-d(i))*(lng(i;theta)+lnF(i;gamma)) ! + d(i)*ln(1-G(i;theta)*F(i;gamma)) ! where: ! d(i)=1 if no recall observed ! =0 if recall observed ! F(i;gamma) = Pr(R=1;z(i)) ! G(i;theta)=Pr(T<=t;x(i),theta|R=1) ! g(i;theta)=dG(i;theta)/dt ! t=t(i) is observed duration, possibly censored ! X (duration) can have more regressors than Z (recall) implicit double precision (a-h,o-z) implicit integer(i-n) double precision par(90),varb(90,90),loglik,err(5000) double precision dur(5000),x(5000,45),z(5000,45) double precision b(45),v(10),c(45),d(5000),se(90) double precision w(0:100,10) double precision start(90),flex(100) double precision seflex(100) integer mlopt(10),pos(60) character*8 cov(60),shift(100) write(*,*)' I exist' write(*,*)' FLAG 1' !Flags indicate progress of program open(unit=34,file='SPLIT.res',status='unknown')!output file open(unit=44,file='split.info',status='unknown')!Instructions/input open(unit=36,file='check.out',status='unknown')!check data manipulations open(unit=48,file='START.VAL',status='unknown')!read in start values write(*,*)' FLAG 2' ! READ IN DATA from LAYOFF.DAT ! n=numb obs; k1=numb x's; k2=numb z's ! k3=numb of ALMON coefficients (order of polynomial) ! MAX defines the MAX number of intervals on the time axis ! with I[0] to I[MAX]: I[j]={s;c(j)= #Z's if(k2.gt.0)then do 6 j=1,k1 ! fill in any k1-k2 coefficients which should not appaer with zero read(48,'(33x,f9.4)')c(j)! coeff for Z (recall) 6 continue end if ! Almon coefficients: also pre-specified (perhaps all zero?) if(k3.gt.0)then do 2 j=1,k3 read(48,'(10x,f13.4)')v(j) 2 continue end if a=1.0d0 ! set alpha=1 in Weibull hazard spec, since using flexible basline hazard rewind 48 ! construct parameter vector: ! note that here nop will include alpha, whether or not it is fixed at 1 nop=k1+k2+k3+1 ! dim of full paramater vector, theta do 3 j=1,k1 start(j)=b(j) 3 continue start(nop)=a if(k2.gt.0)then do 4 j=1,k2 jj=k1 4 start(jj+j)=c(j) end if if(k3.gt.0)then do 5 j=1,k3 jj=k1+k2 5 start(jj+j)=v(j) end if ! call MAXLIK and estimate model kreps=1 iscr=1 ! mlopt(.) gives various instructions to ML routine MAXLIK mlopt(1)=kreps mlopt(2)=iscr mlopt(3)=ictrl mlopt(4)=iopg mlopt(5)=MAX mlopt(6)=iwid mlopt(7)=ihaz write(*,*)' CALL THE ONE AND ONLY MAXLIK!' call maxlik(start,k1,k2,k3,n,dur,d,x,z,mlopt,& loglik,par,varb,err,het,cov) ! loglik contains maximised log-likelihood ! par contains MLE: dim = NOP ! varb contains est var matrix of MLE's ! dl contains contributions to scores ! d2l contains contributions to hessian ! err contains gen errors : unit exponential in the absence of censoring ! d=0 if censored, d=1 if complete do 100 j=1,nop-ihaz! for generality, this allows alpha N.E. 1 100 se(j)=sqrt(varb(j,j)) if(ihaz.eq.1)se(nop)=0.0d0 if(max.gt.0)then !mmmmmmmmmmmmmmmmmmmmmmmmmmmmm call almon(w,MAX,k3)! create ALMON matrix, W(max,k3): flex=W*v, use parameter v in ML routine kk=k1+k2 do 60 j=1,MAX flex(j)=0.0d0 do 60 k=1,k3 j1=kk+k flex(j)=flex(j)+w(j,k)*par(j1)! construct hazard shifters from estimated v 60 continue do 61 j=1,max seflex(j)=0.0d0! construct StErrs for hazars shifters do 62 k=1,k3 do 62 m=1,k3 seflex(j)=seflex(j)+w(j,k)*varb(kk+k,kk+m)*w(j,m) 62 continue write(*,'(f8.4)')seflex(j) seflex(j)=sqrt(seflex(j)) 61 continue end if !mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm ! Test coefficient restrictions on ALG and ALH ! WH1 N(0,1) test of ALG=ALH coefficients in HAZARD WH1=(par(32)-par(33))/& sqrt(varb(32,32)+varb(33,33)-2.0d0*varb(32,33)) ! only do the following for SPLIT model WR1=(par(65)-par(66))/& sqrt(varb(65,65)+varb(66,66)-2.0d0*varb(65,66)) !******************************************************* ! Now write results to output file ! Unit=34 !******************************************************* write(34,1066) 1066 format(/' MAXIMUM LIKELIHOOD ESTIMATION RESULTS: MALES'/& ' SPLIT POPULATION MODEL') ! ' SINGLE POPULATION MODEL') write(34,1666)n 1666 format(/' N=',i4) if(max.gt.0)then max1=max+1 write(34,2066)max1,iwid,k3 2066 format(/,'NUMBER OF INTERVALS,INCLUDING OPEN END =',I4,/& 'WITH WIDTH OF',I4,' DAYS.'/& 'ORDER OF ALMON POLYNOMIAL=',I3) end if if(ihaz.eq.1)then write(34,*)'*alpha FIXED at ONE*' end if write(34,1067)loglik 1067 format(/,'LOG-LIK=',f15.4,/) write(34,1068)het 1068 format('Standard Normal test for overdispersion=',f10.4/) write(34,1069)WH1,WR1,& par(32),par(65),par(33),par(66) 1069 format('Wald test statistics [standard normal]'/& 'ALG=ALH in HAZARD:',f10.4,2x,/& 'ALG=ALH in RECALL:',f10.4/& 'Unconstrained estimates are:'/& 'Hazard ALG :',f10.4,2x,'Recall ALG :',f10.4/& 'Hazard ALH :',f10.4,2x,'Recall ALH :',f10.4/) write(34,1166) 1166 format(/,16x,'DURATION',10X,'RECALL PROBABILITY'/& 14x,'MLE',5x,'Sterr',10x,'MLE',5x,'Sterr') if(k2.gt.0)then! k2>0 in Split Pop Model do 1266 j=1,k1 jm=k1+j pjm=par(jm) sejm=se(jm) if(j.gt.k2)then pjm=0.0d0 sejm=0.0d0 end if 1266 write(34,1366)cov(j),par(j),se(j),pjm,sejm else do 1966 j=1,k1 1966 write(34,1366)cov(j),par(j),se(j) end if 1366 format(a8,2x,2f9.4,5X,2f9.4) if(max.gt.0)then do 1466 j=1,max 1466 write(34,1566)shift(j),j,flex(j),seflex(j) end if 1566 format(a8,i2,2x,2f9.4) write(34,1766)par(nop),se(nop) 1766 format(' alpha ',2x,2f9.4) write(34,2166) 2166 format(/,'ALMON POLYNOMIAL ESTIMATES') if(max.gt.0)then do 2466 j=1,k3 j1=k1+k2+j 2466 write(34,2566)j,par(j1),se(j1) 2566 format(6x,i2,2x,2f13.4) end if write(34,*)' ===========================================' write(34,'(/)') !**** end of writing results ********** !****************************************************** STOP end PROGRAM LAYOFF ! end of main unit !************************************************************** ! SUBROUTINES FOLLOW !mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm subroutine getdata(dur,d,n,pos,cov,x,k,z,m) ! getdata(dur,d,n,x,k,z,m,wt) implicit double precision (a-h,o-z) implicit integer(i,n) double precision dur(5000),reg(5000,60),z(5000,45) double precision x(5000,45),ave(60),std(60) double precision d(5000),zt(0:3),tau(0:4) integer pos(40) character*8 cov(60),vbl(80) open(unit=40,file='TLAYOFF.DAT',status='old')! original data file open(unit=42,file='test.prog',status='unknown')! check data manipulations open(unit=38,file='descstats.dat',status='unknown')! descriptive statistics ! names of covariates (8 length max) vbl(1)='const ' vbl(2)='national' vbl(3)='edprof2 ' vbl(4)='edhigh2 ' vbl(5)='age17-24' vbl(6)='age25-34' vbl(7)='age35-44' vbl(8)='shortill' vbl(9)='longill ' vbl(10)='sperr ' vbl(11)='public ' vbl(12)='chemical' vbl(13)='metal ' vbl(14)='electric' vbl(15)='food ' vbl(16)='textiles' vbl(17)='trade ' vbl(18)='construc' vbl(19)='agricul ' vbl(20)='S-H ' vbl(21)='Hamburg ' vbl(22)='Nieders ' vbl(23)='Bremen ' vbl(24)='Nrhein-W' vbl(25)='Hessen ' vbl(26)='R-P ' vbl(27)='Baden-W ' vbl(28)='Saarland' vbl(29)='Berlin ' vbl(30)='manual ' vbl(31)='skilled ' vbl(32)='seq ' vbl(33)='train ' vbl(34)='algplot ' vbl(35)='married ' vbl(36)='children' vbl(37)='wwgr ' vbl(38)='pay2d ' vbl(39)='pay3d ' vbl(40)='gender ' vbl(41)='ben0 ' vbl(42)='ben1 ' vbl(43)='ben2 ' vbl(44)='ben3 ' vbl(45)='time0 ' vbl(46)='time1 ' vbl(47)='time2 ' vbl(48)='time3 ' ! specially constructed benefit indicator vbl(50)='ALG ONLY' vbl(51)='ALH ' vbl(52)='ALG&ALH ' vbl(60)=' RECALL ' vbl(80)='duration' vbl(79)='new hire' vbl(78)=' recall ' vbl(77)=' died ' do 11 j=1,k+4 ave(j)=0.0d0 std(j)=0.0d0 11 continue do 100 i=1,n! read in from TLAYOFF.DAT write(*,*)i reg(i,1)=1.0d0 read(40,1)(reg(i,j),j=2,31),died,hn,re,& reg(i,32),reg(i,33),& dur(i),(reg(i,j),j=34,40),& (zt(j),j=0,3),(tau(j),j=1,4),& t1,t2,t3,t4,ftot2 ! zt(j)=ben(j) and tau(j)=time(j) 1 format(5(f8.0,6f10.0,/)& f8.0,5f10.0,f14.0,/f8.0,6f10.0,/f8.0,6f10.0) d(i)=1.0d0-re alg=0.0d0! CONSTRUCT ALG VARIABLE if(t1.eq.1.0d0.and.t2.eq.2.0d0)alg=1.0d0 if(t1.eq.2.0d0)alg=1.0d0 if(t1.eq.0.0d0.and.t2.eq.0.0d0.and.& t3.eq.0.0d0.and.t4.eq.0.0d0)alg=1.0d0 reg(i,50)=alg alh=0.0d0! CONSTRUCT ALH VARIABLE if(t1.eq.1.0d0.and.t2.eq.3.0d0)alh=1.0d0 if(t1.eq.3.0d0)alh=1.0d0 if(t1.eq.1.0d0.and.t2.eq.2.0d0.and.& t3.eq.3.0d0)alh=1.0d0 if(t1.eq.2.0d0.and.t2.eq.3.0d0)alh=1.0d0 reg(i,51)=alh ! create x(50)=1 only if alg=1 and alh=0 if(alh.eq.1.0d0)reg(i,50)=0.0d0! ALG=1 if ALG only!! reg(i,45)=tau(1) reg(i,46)=tau(2) reg(i,47)=tau(3) reg(i,48)=tau(4) rec=dble(i)! data/record indicator do 8 j=1,k ip=pos(j)! get covariates that are going to be sued as regressors cov(j)=vbl(ip) ave(j)=ave(j)+reg(i,ip) std(j)=std(j)+(reg(i,ip)**2) ! X(.) contains the duration part regressors ! Z(.) the recall part regressors ! the first k columns of X must equal Z ! X can have more variables than Z 8 x(i,j)=reg(i,ip) ! note that ALG & ALH are the last regressors in that order ave(k+1)=ave(k+1)+died ave(k+2)=ave(k+2)+re ave(k+3)=ave(k+3)+hn ave(k+4)=ave(k+4)+dur(i) std(k+1)=std(k+1)+(died**2) std(k+2)=std(k+2)+(re**2) std(k+3)=std(k+3)+(hn**2) std(k+4)=std(k+4)+(dur(i)**2) do 7 j=k+1,45 7 x(i,j)=0.0d0! sets unused space in X equal to zero do 6 j=1,m z(i,j)=x(i,j)! Z=X 6 continue 100 continue ! finish off descriptiove stats do 12 j=1,k ave(j)=ave(j)/dble(n) std(j)=sqrt((std(j)/dble(n))-(ave(j)**2)) 12 continue ave(k+1)=ave(k+1)/dble(n) ave(k+2)=ave(k+2)/dble(n) ave(k+3)=ave(k+3)/dble(n) ave(k+4)=ave(k+4)/dble(n) std(k+1)=std(k+1)/dble(n) std(k+2)=std(k+2)/dble(n) std(k+3)=std(k+3)/dble(n) std(k+4)=std(k+4)/dble(n) std(k+1)=sqrt(std(k+1)-(ave(k+1)**2)) std(k+2)=sqrt(std(k+2)-(ave(k+2)**2)) std(k+3)=sqrt(std(k+3)-(ave(k+3)**2)) std(k+4)=sqrt(std(k+4)-(ave(k+4)**2)) ! write desc stats to unit 44 write(38,*)' MEANS and STANDARD DEVIATIONS' write(38,*)' Sample size = 5000' do 15 j=1,k ip=pos(j) write(38,'(a8,2x,2f10.4)')vbl(ip),ave(j),std(j) 15 continue write(38,'(a8,2x,2f10.4)')vbl(77),& ave(k+1),std(k+1) write(38,'(a8,2x,2f10.4)')vbl(78),& ave(k+2),std(k+2) write(38,'(a8,2x,2f10.4)')vbl(79),& ave(k+3),std(k+3) write(38,'(a8,2x,2f10.4)')vbl(80),& ave(k+4),std(k+4) ! !!!!! d=1-cens end getdata !mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm ! MAXLIK subroutine subroutine maxlik(start,k1,k2,k3,n,dur,d,x,z,mlopt,& loglik,par,varb,err,het,cov) implicit double precision (a-h,o-z) implicit integer(i,n) double precision par(90),varb(90,90),loglik,err(5000) double precision grad(90),hess(90,90),dir(90),tag(90) double precision dur(5000),x(5000,45),z(5000,45) double precision b(45),c(45),d(5000) double precision w(0:100,10),one(5000),start(90) double precision ba(90),trat(90),v(10) double precision den(10),score(5000,90) double precision tau(0:101),xint(5000,0:101) integer mlopt(10),mn(5000) character*8 cov(60) external s15abf open(unit=42,file='hess.out',status='unknown')! test file open(unit=48,file='current.VAL',status='unknown')! outputs current iterates ! aide memoire: ! kreps=1 ! iscr=1 ! mlopt(1)=kreps ! mlopt(2)=iscr ! mlopt(3)=ictrl ! mlopt(4)=iopg ! mlopt(5)=MAX ! mlopt(6)=iwid ! mlopt(7)=ihaz ictrl=mlopt(3) iscr=mlopt(2) kreps=mlopt(1) iopg=mlopt(4) MAX=mlopt(5) iwid=mlopt(6) ihaz=mlopt(7) const=1.0d0/sqrt(2.0d0*3.1415926d0)! =1/sqrt(2*pi) write(*,*)' FLAG MAX-LIK 1' nop1=k1+k2+k3+1 do 13 j=1,nop1 write(*,'(I4,2x,8f10.4)')j,start(j) 13 par(j)=start(j) call parvec(par,b,a,v,c,k1,k2,k3)! get full parameter vector, theta nop=nop1-ihaz ! construct ALMON matrix if(max.gt.0)then call almon(w,MAX,k3) end if ! construct the global intervals on the time axis tau(0)=0.0d0 if (max.gt.0)then do 1 j=1,MAX 1 tau(j)=tau(j-1)+dble(iwid) end if stem=0.0d0 write(*,*)' FLAG MAX-LIK 2' ! construct local/individual grids given obs. duration=t ! m is such that t falls in interval I[m], where intervals run ! from I[0],I[1],...,I[MAX+1]; and we calculate tau(m+1)=t, for each i call grid(dur,n,tau,xint,MAX,mn)! get time interval grid ! start iterations oldll=0.0d0 L=0 icount=0 ifin=0 100 write(*,*)' FLAG MAX-LIK 3' ! generate required quantities for individual i do 52 j=1,nop grad(j)=0.0d0 do 51 jj=1,nop 51 hess(j,jj)=0.0d0 52 continue write(*,*)' FLAG MAX-LIK 4' !=============================================== ! REMINDER: Log-lik contributions: ! l(i)=(1-d(i))*(lng(i;theta)+lnF(i;gamma)) ! + d(i)*ln(1-G(i;theta)*F(i;gamma)) ! where: ! d(i)=1 if no recall observed ! =0 if recall observed !=============================================== do 200 i=1,n t=dur(i) r=d(i) ir=int(r) ! get local.individual grids from the xint matrix m=mn(i) do 4 j=0,m+1 4 tau(j)=xint(i,j) ! set up contributions to 1st and 2nd derivatives ! integrated hazard and other components xb=0.0d0 do 2 j=1,k1 xj=x(i,j) bj=b(j) 2 xb=xb+xj*bj exb=exp(xb) if(k2.gt.0)then zc=0.0d0 do 3 j=1,k2 zj=z(i,j) cj=c(j) 3 zc=zc+zj*cj end if ! construct e=-ln(SURVIVOR)=INT(HAZARD) and various derivatives ! to use in GRADIENT subroutine; ! duration bits first if(max.gt.0)then !mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm e=0.0d0 dea=0.0d0 do 10 js=1,k3 10 den(js)=0.0d0 do 5 j=0,m wvj=0.0d0 do 6 js=1,k3 wjs=w(j,js) vs=v(js) wvj=wvj+wjs*vs 6 continue qj=exp(wvj) sj=qj*((tau(j+1)**a)-(tau(j)**a)) do 7 js=1,k3 wjs=w(j,js) den(js)=den(js)+exb*sj*wjs 7 continue tb=(tau(j+1)**a)*log(tau(j+1)) if(j.eq.0)then ta=0.0d0 else ta=(tau(j)**a)*log(tau(j)) end if pj=tb-ta dea=dea+qj*pj e=e+sj 5 continue e=exb*e dea=exb*dea else e=exb*(t**a) dea=e*log(t) ! end if max.gt.0 !mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm end if G=1.0d0-exp(-e) ! now recall probbailities if(k2.gt.0)then ifail=0 F=s15abf(zc,ifail) dF0=const*exp(-0.5d0*(zc**2)) P=1.0d0-(G*F) else F=1.0d0 dF0=0.0d0 P=exp(-e) end if ! construct gradient contributions do 20 j=1,nop if(j.le.k1)then dej=e*x(i,j) dhj=x(i,j) end if if(j.eq.k1+k2+k3+1)then dej=dea dhj=(1.0d0/a)+log(t) end if if(j.gt.k1+k2.and.j.le.k1+k2+k3)then jj=j-(k1+k2) dej=den(jj) dhj=w(m,jj) end if if(j.gt.k1.and.j.le.k1+k2)then jj=j-k1 dlF=(dF0/F)*z(i,jj) dF=dF0*z(i,jj) end if ! score constributions call gradi(dlj,j,k1,k2,k3,dej,dhj,F,dlF,dF,G,P,r) grad(j)=grad(j)+dlj !***************************************** if(ifin.eq.1)then score(i,j)=dlj end if !***************************************** ! Now construct OPG covariance matrix approx. to hessian do 21 j1=j,nop if(j1.le.k1)then dej=e*x(i,j1) dhj=x(i,j1) end if if(j1.eq.k1+k2+k3+1)then dej=dea dhj=(1.0d0/a)+log(t) end if if(j1.gt.k1+k2.and.j.le.k1+k2+k3)then jj=j1-(k1+k2) dej=den(jj) dhj=w(m,jj) end if if(j1.gt.k1.and.j.le.k1+k2)then jj=j1-k1 dlF=(dF0/F)*z(i,jj) dF=dF0*z(i,jj) end if ! write(*,*)' call gradi' call gradi(dlj1,j1,k1,k2,k3,dej,dhj,F,dlF,dF,G,P,r) hess(j,j1)=hess(j,j1)-dlj*dlj1 ! OPG form hess(j1,j)=hess(j,j1) 21 continue 20 continue if(ifin.eq.1)then one(i)=1.0d0 if(r.eq.0.0d0)then score(i,nop+1)=((1.0d0-e)**2)-e end if if(r.eq.1.0d0)then score(i,nop+1)=e*(e-1.0d0)*F*(1.0d0-G)/P end if end if 200 continue write(*,*)' FLAG MAX-LIK 5' ! test calculations write(42,'(10f7.1)')(par(j),j=1,nop) write(42,'(8f10.1)')(grad(j),j=1,nop) do 5000 j=1,nop write(42,'(//8f10.1)')(hess(j,jj),jj=1,nop) 5000 continue ! get log-likelihood @ current parameter values CALL LOGL(loglik,dur,d,x,z,w,par,k1,k2,k3,n,xint,mn) write(*,*)' FLAG MAX-LIK 6' 7095 IF(ISCR.EQ.1)then WRITE(6,7889) end if 7889 FORMAT(' CALLING NEWT RAPHSON FOR UPDATES') CALL NR(DIR,NOP,grad,hess,varb,ESS) ! NR converts HESSIAN to -HESSIAN in construction of DIR ! so that update is par(t+1)=par(t)+step(t)*dir(t) rewind 42 IF(ifin.eq.0)then !**** IFIN=0 ****************** write(*,*)' FLAG MAX-LIK 7' !===================================================== ! temporary storage of current estaimtes in START.VAL do 111 j=1,k1 if(k2.gt.0)then j2=k1+j write(48,'(a8,2x,f9.4,14x,f9.4)')& cov(j),par(j),par(j2) else write(48,'(a8,2x,f9.4,)')cov(j),par(j) end if 111 continue ! Almon coefficients if(k3.gt.0)then do 211 j=1,k3 j3=k1+k2+j write(48,'(10x,F13.4)')par(j3) 211 continue end if rewind 48 !===================================================== IF(ISCR.EQ.1)THEN! write iteration results to screen WRITE(*,7992)LOGLIK WRITE(*,7999) do 7001 j=1,nop 7001 WRITE(*,'(i8,3(2X,F10.4))')j,par(J),grad(j),DIR(j) write(*,7997)L write(*,7800)OLDLL,LOGLIK write(*,7993) WRITE(*,'(2X,f7.3,7x,F15.9,//)')stem,ESS 7800 FORMAT(' OLD LOG-LIK=',1X,F20.9/& ' NEW LOG-LIK=',1X,F20.9) 7997 FORMAT(' ITER=',I3) 7992 FORMAT(' LOG-LIKELIHOOD=',1X,F20.9) 7999 FORMAT(11x,' ESTIMATE',7x,'GRAD',7x,'DIR') 7995 FORMAT(' FAIL TO CONVERGE AFTER 100 ITERATIONS') 7993 FORMAT(' STEP LENGTH & NORMED GRADIENT') END IF ! Now search along a step length grid fashion oldll=loglik write(*,*)' FLAG MAX-LIK 8' if(ictrl.eq.1)then write(*,*)' Input desired step length' read(*,*)stem end if if(ictrl.eq.0)then STEP=0.0d0 STEM=0.1d0 DO 29 J=1,20 STEP=STEP+0.1d0 DO 28 I=1,NOP 28 TAG(I)=par(I)+STEP*DIR(I) if(ihaz.eq.1)tag(nop1)=par(nop1) CALL LOGL(TEMP,dur,d,x,z,w,tag,k1,k2,k3,n,xint,mn) IF(TEMP.GT.logLIK)THEN ! check for improvement logLIK=TEMP STEM=STEP END IF 29 CONTINUE end if write(*,*)' FLAG MAX-LIK 9' ! ending search of step length DO 30 I=1,NOP 30 par(I)=par(I)+STEM*DIR(I) call parvec(par,b,a,v,c,k1,k2,k3) IF(ESS.GT.1.0d-6)THEN L=L+1 if(l.gt.100)then WRITE(6,7995) NafF=1 else GO TO 100 end if END IF ! iteration count mlopt(4)=L GO TO 4000 !*** in calling newt-Raphs sub-rout HEss has changed to (-HESS)^-1 *** 4000 write(*,*)' Convergence obtained' if(icount.ne.L)then write(*,*)' Updating var-cov matrix' ifin=1 icount=L go to 100 end if !********** END OF IFIN=0 ****************** end if ! call OLS routine for test statistic (neg heterogeneity) not=nop+1 call OPG(score,one,ba,trat,NOP,Not,N,ess) het=trat(not) write(*,*)' neg-het: coeff & tratio' write(*,'(2(2x,f10.5))')ba(not),het 4001 write(*,'(//)') end maxlik !------------------------------------------------------------ ! Construct the full parameter vector, theta subroutine parvec(par,b,a,v,c,k1,k2,k3) implicit double precision (a-h,o-z) implicit integer (i-n) double precision par(90),b(45),c(45),v(10) do 11 j=1,k1 11 b(j)=par(j) a=par(k1+k2+k3+1) if(k3.gt.0)then do 12 j=1,k3 jj=k1+k2 12 v(j)=par(jj+j) end if if(k2.gt.0)then do 13 j=1,k2 jj=k1 13 c(j)=par(jj+j) end if end parvec !------------------------------------------------------------- ! Almon Matrix subroutine almon(w,m,k) implicit double precision (a-h,o-z) implicit integer(i-n) double precision w(0:100,10) do 1 j=0,m xj=dble(j) do 1 i=1,k p=0.000001d0 1 w(j,i)=p*(xj**i) end almon !------------------------------------------------------------- ! divide up the time interval, with obs. duration being the upper bound ! on the last interval for individual i subroutine grid(dur,n,tau,x,MAX,mn) implicit double precision (a-h,o-z) implicit integer(i-n) double precision dur(5000),x(5000,0:101) double precision tau(0:101) integer mn(5000) do 2 i=1,n t=dur(i) ! m defines then interval in which t lies m=-1 do 1 j=0,max ! write(*,'(2x,i4,f8.2)')m,tau(j) if(t.gt.tau(j))then x(i,j)=tau(j) m=m+1 end if 1 continue mn(i)=m x(i,m+1)=t 2 continue end grid !------------------------------------------------------------- ! Contributions to score subroutine gradi(dl,j,k1,k2,k3,dej,dhj,F,dlF,dF,G,P,r) implicit double precision (a-h,o-z) implicit integer (i-n) double precision de(4) ! Work out derivative for r=0.0 then r=1.0 IF(r.eq.0.0d0)then !RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR ! ! r=0,if recall to same employer observed ! r=1, if recall not observed (censored by sample frame or ! by exit to another employer) if(j.gt.k1.and.j.le.k1+k2)then ! derivatives wrt gamma r=0 ! dl=dlF else ! derivatives wrt beta,alpha, or eta r=0 dl=dhj-dej end if !RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR else ! if r=1, that is ... S=1.0d0-G if(j.gt.k1.and.j.le.k1+k2)then ! derivatives wrt gamma r=0 ! dl=-G*dF/P else ! derivatives wrt beta,alpha, or eta r=0 dl=-F*S*dej/P end if !RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR end if end gradi !------------------------------------------------------------- ! Construct log-likelihood at current parameter values subroutine LOGL(fn,dur,d,x,z,w,par,k1,k2,k3,n,xint,mn) implicit double precision (a-h,o-z) implicit integer (i-n) double precision dur(5000),x(5000,45),z(5000,45) double precision wv(0:100) double precision b(45),c(45),d(5000),par(90),v(10) double precision w(0:100,10),xint(5000,0:101),tau(0:101) integer mn(5000) external s15abf call parvec(par,b,a,v,c,k1,k2,k3) ! write(*,*)' log-lik routine' ! write(*,'(8f10.4)')(par(j),j=1,k1+k2+k3+1) fn=0.0d0 do 20 i=1,n t=dur(i) r=d(i) ir=int(r) ! write(*,'(2f10.4)')t,r m=mn(i) do 4 j=0,m+1 4 tau(j)=xint(i,j) ! set up contributions xb=0.0d0 do 2 j=1,k1 xj=x(i,j) bj=b(j) 2 xb=xb+xj*bj exb=exp(xb) if(k2.gt.0)then zc=0.0d0 do 3 j=1,k2 zj=z(i,j) cj=c(j) 3 zc=zc+zj*cj end if ! construct e=-ln(SURVIVOR)=INT(HAZARD) if(m.gt.0)then e=0.0d0 do 5 j=0,m wv(j)=0.0d0 do 6 js=1,k3 wjs=w(j,js) vs=v(js) wv(j)=wv(j)+wjs*vs 6 continue qj=exp(wv(j)) sj=qj*((tau(j+1)**a)-(tau(j)**a)) e=e+sj 5 continue e=exb*e haz=log(a)+a*log(t)+xb+wv(m) else e=exb*(t**a) haz=log(a)+a*log(t)+xb end if G=1.0d0-exp(-e) ! now recall probbailities if(k2.gt.0)then ifail=0 F=s15abf(zc,ifail) Fcont=log(F) P=1.0d0-(G*F) else Fcont=0.0d0 P=exp(-e) end if if(r.eq.0.0d0)then fn=fn+haz-e+Fcont else fn=fn+log(P) end if 20 continue end logl !------------------------------------------------------------------ ! generalised errors from PH spec of model subroutine generri(e,G,i,t,x,b,a,k) implicit double precision (a-h,o-z) implicit integer (i-n) double precision x(5000,45),b(45) xb=0.0d0 do 1 j=1,k 1 xb=xb+x(i,j)*b(j) e=exp(xb)*(t**a) G=1.0d0-exp(-e) end generri !----------------------------------------------------------------- ! initialise matrix for inversion subroutine inita(a) double precision a(100,100) DO 804 i=1,100 DO 805 j=1,100 805 A(i,j)=0.0d0 804 A(i,i)=1.0d0 return end inita !----------------------------------------------------------------- ! Newton Raphson iterative scheme (possibly BHHH) SUBROUTINE NR(BA,NOX,grad,HESS,VARB,ESS) double precision HESS(90,90),GRAD(90),BA(90),VARB(90,90) double precision AI(100,100),CI(100,100) double precision ess,tss,step INTEGER NOX,n,i,j,f call inita(ai) DO 800 J=1,NOX DO 801 F=1,NOX AI(J,F)=-HESS(J,F) 801 CONTINUE 800 CONTINUE CALL INV(AI,CI,NOX) ESS=0.0d0 DO 870 J=1,NOX BA(J)=0.0d0 DO 871 F=1,NOX VARB(J,F)=CI(J,F) ESS=ESS+GRAD(J)*CI(J,F)*GRAD(F) 871 BA(J)=BA(J)+CI(J,F)*GRAD(F) 870 CONTINUE END nr !------------------------------------------------------------------------ ! BASIC OLS routine used to generate NEG HET test statistic SUBROUTINE OPG(X,Y,b,trat,NOP,NOX,N,ess) implicit double precision (a-h,o-z) double precision X(5000,90),Y(5000),trat(90) double precision AI(100,100),CI(100,100),B(90),res(5000) double precision XTX(90,90),XTY(90),VAR(90,90) double precision ess,tss INTEGER NOX,n,i,j,f call inita(ai) DO 800 J=1,NOX DO 801 F=1,NOX 801 XTX(J,F)=0.0d0 800 XTY(J)=0.0d0 tss=0.0d0 DO 850 I=1,N tss=tss+(y(i)**2) DO 851 J=1,NOX dj=x(i,j)*y(i) DO 852 F=1,NOX 852 XTX(J,F)=XTX(J,F)+X(I,J)*X(I,F) 851 XTY(J)=XTY(J)+dj 850 CONTINUE DO 860 J=1,NOX DO 861 F=1,NOX 861 AI(J,F)=XTX(J,F) B(J)=0.0d0 860 CONTINUE CALL INV(AI,CI,NOX)! invert matrix ESS=0.0d0 DO 870 J=1,NOX DO 871 F=1,NOX VAR(J,F)=CI(J,F) ESS=ESS+XTY(J)*CI(J,F)*XTY(F) 871 B(J)=B(J)+CI(J,F)*XTY(F) trat(j)=b(j)/sqrt(ci(j,j)) 870 CONTINUE rss=tss-ess do 900 i=1,n res(i)=y(i) do 900 j=1,nox 900 res(i)=res(i)-x(i,j)*b(j) return END opg !------------------------------------------------------------------------ ! Matrix inversion (method 1) SUBROUTINE INV(A,C,K) double precision X(100,100),Y(100,100),Z(100,100) double precision W(100),XX(100,100),YY(100,100) double precision A(100,100),C(100,100) INTEGER I,J,K,M,SING DO 2110 I=1,K Y(I,I)=1.0d0 DO 2120 J=1,K 2120 X(I,J)=A(I,J) DO 2115 J=I+1,K Y(I,J)=0.0d0 2115 Y(J,I)=0.0d0 2110 CONTINUE IFAIL=0 CALL F04AEf(X,100,Y,100,K,K,Z,100,W,XX,100,YY,100,IFAIL) IF(IFAIL.EQ.1)THEN SING=IFAIL RETURN END IF DO 2130 I=1,K DO 2140 J=1,K 2140 C(I,J)=Z(I,J) 2130 CONTINUE END inv !................................................................... ! Matrix Inversion (Method 2) subroutine invert(a,n) implicit none double precision a(100,100),amax,eps integer ipiv(65),n,n5,i,j,k,icol do 1 i=1,n 1 ipiv(i)=0 do 14 i=1,n amax=0.0d0 do 2 j=1,n if(ipiv(j))3,4,2 3 continue stop 4 if(dabs(a(j,j))-amax)6,6,7 7 icol=j amax=dabs(a(j,j)) 6 continue 2 continue ipiv(icol)=1 n5=n+5 eps=0.01d0**n5 if(amax-eps)8,8,10 8 write(32,9) 9 format(1h ,' singular matrix ') go to 200 ! stop 10 continue amax=a(icol,icol) a(icol,icol)=1.0d0 do 11 k=1,n 11 a(icol,k)=a(icol,k)/amax do 14 j=1,n if(j-icol)12,14,12 12 amax=a(j,icol) a(j,icol)=0.0d0 do 13 k=1,n 13 a(j,k)=a(j,k)-a(icol,k)*amax 14 continue 200 continue return end invert !++++++++++++++++++++++++++++++++++++++++++++++++++++ ! THAT'S ALL FOLKS! !++++++++++++++++++++++++++++++++++++++++++++++++++++