! 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)<s<=c(j+1)} & c(MAX+1)=infinity
! iwid=interval width (common)
! MAX<=100,k3<=10
      ! From SPLIT.INFO
      read(44,*)n,k1,k2
      read(44,*)MAX,iwid,k3
      read(44,*)ihaz
      read(44,*)iopg
      read(44,*)ictrl
      read(44,*)(pos(j),j=1,k1)

      write(*,*)' get data now ?'
! generate data: "dur, X, Z, d(=recall or not) etc from original data file TLAYOFF.DAT
      call getdata(dur,d,n,pos,cov,x,k1,z,k2)
      
      do 111 j=1,max
 111  shift(j)='shift   '   ! label hazard shifters for results file
 
      write(*,*)' Have got data?'

! set starting parameter values, perhaps from old results in START.VAL
       do 1 j=1,k1
        read(48,'(10x,f9.4,)')b(j)! coeff for X (duration)
        c(j)=0.0d0
 1     continue
       rewind 48
       ! #X's >= #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!
!++++++++++++++++++++++++++++++++++++++++++++++++++++