!Estimation was conducted with the IMSL routine UMINF.  
!The following subroutine contains the likelihood function that was maximized by UMINF.
!The "Use firm_data" command refers to a Fortran module that contained the proprietary data.


Subroutine lik_fcn(n, x, f)
    Use firm_data
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: n
    REAL, INTENT(IN) :: x(n)
    REAL, INTENT(OUT) :: f
 
  !Local Variable Declarations
    INTEGER :: obs,i,g,j
    INTEGER, DIMENSION(7) :: good
    REAL, DIMENSION(85391,3,0:127) :: lambdal, lambdan, lambda, psia, psip, mup, eqc, omega, muomega, &
                                  & sigmap, fc, delta, dep, l, lprod, prob, weight
    REAL, DIMENSION(18595) :: lik
   
 
   lik=0
   F = 0
  
   Do obs = 1,85391
   Do i = 1,3
   Do g = 0, 127

  ! Compute lambdas for claims prob
	lambdal(obs,i,g) = x(i)
	lambdan(obs,i,g)=0

  	Do j=4,9
	 lambdal(obs,i,g)=lambdal(obs,i,g)+x(j)*RC(j,obs)
	 lambdan(obs,i,g)=lambdan(obs,i,g)+x(6+j)*RC(j,obs)
 	End Do

	lambda(obs,i,g) = lambdal(obs,i,g)+lambdan(obs,i,g)

  ! Parse g-step to determine good service by claim
	Good(7) = Floor(g/2.0**6)
	Good(6) = Floor((g-(2.0**6)*good(7))/2.0**5)
	Good(5) = Floor((g-(2.0**6)*good(7)-(2.0**5)*good(6))/2.0**4)
	Good(4) = Floor((g-(2.0**6)*good(7)-(2.0**5)*good(6)-(2.0**4)*good(5))/2.0**3)
	Good(3) = Floor((g-(2.0**6)*good(7)-(2.0**5)*good(6)-(2.0**4)*good(5)-(2.0**3)*good(4))/2.0**2)
	Good(2) = Floor((g-(2.0**6)*good(7)-(2.0**5)*good(6)-(2.0**4)*good(5)-(2.0**3)*good(4)-(2.0**2)*good(3))/2.0)
	Good(1) = Floor((g-(2.0**6)*good(7)-(2.0**5)*good(6)-(2.0**4)*good(5)-(2.0**3)*good(4)-(2.0**2)*good(3)-(2.0**1)*good(2)))

  ! Compute claims service expectation, given learning from good service draws by claim
	psia(obs,i,g) = x(16)*(age(obs)-16) + x(17)*((age(obs)-16)**2)
	psip(obs,i,g) = x(18) + psia(obs,i,g)
	mup(obs,i,g) = (x(18)*x(19) + psia(obs,i,g)*x(20))/psip(obs,i,g)

	eqc(obs,i,g) = ((mup(obs,i,g)*psip(obs,i,g))/(psip(obs,i,g)+tot_clms(obs)))
	Do j=1,tot_clms(obs)
	 eqc(obs,i,g) = eqc(obs,i,g) + Good(j)/(psip(obs,i,g)+tot_clms(obs))
	End Do

   ! Compute Departure Probability
	omega(obs,i,g) = x(21) + x(22)*srcht(obs)
	muomega(obs,i,g) = x(19) + x(23)*(age(obs)-16) + x(24)*((age(obs)-16)**2)
  
  	sigmap(obs,i,g) = x(25)*pmarket(obs)
  
 	fc(obs,i,g) = x(26) + x(27)*mcar(obs)+x(28)*fp(obs)+x(29)*lambdal(obs,i,g)+x(30)*tennc(obs)
	do j=1,15
	 fc(obs,i,g) = fc(obs,i,g) + x(30+j)*tenid(j,obs)
	End Do

	delta(obs,i,g) = x(46)*lambda(obs,i,g)*(eqc(obs,i,g) - muomega(obs,i,g)) + fc(obs,i,g)
  	dep(obs,i,g) = 1 - ((ANORDF(((pmarket(obs) - pc(obs)) + delta(obs,i,g))/sigmap(obs,i,g)))**omega(obs,i,g))

  !Compute likelihood contribution
 	l(obs,i,g) = (dep(obs,i,g)**drop(obs)) * ((1-dep(obs,i,g))**fin(obs)) * (lambdal(obs,i,g)**loss(obs)) * &
		     &((1-lambdal(obs,i,g))**noloss(obs)) * (lambdan(obs,i,g)**nloss(obs)) * ((1-lambdan(obs,i,g))**nonloss(obs))


  !Compute likelihood function
  	IF (obs==1) THEN
	    lprod(obs,i,g) = l(obs,i,g)
	ELSE IF (id(obs) /= id(obs-1)) THEN
	    lprod(obs,i,g) = l(obs,i,g)
  	ELSE
    	    lprod(obs,i,g) = lprod(obs-1,i,g)*l(obs,i,g)
 	END IF

	IF (obs < 85391) THEN
          IF (id(obs) /= id(obs+1)) THEN
	    IF (i==1 .OR. i==2) THEN
       		prob(obs,i,g) = x(46+i)
	    ELSE
       		prob(obs,i,g) = 1 - x(47) - x(48)
    	    END IF
    	    weight(obs,i,g)=prob(obs,i,g)*x(20)**good(1)*(1-x(20))**(1-good(1))*x(20)**good(2)*(1-x(20))**(1-good(2))*&
					 &x(20)**good(3)*(1-x(20))**(1-good(3))*x(20)**good(4)*(1-x(20))**(1-good(4))*&
					 &x(20)**good(5)*(1-x(20))**(1-good(5))*x(20)**good(6)*(1-x(20))**(1-good(6))*&
					 &x(20)**good(7)*(1-x(20))**(1-good(7))
	    lik(id(obs)) = lik(id(obs))+weight(obs,i,g)*lprod(obs,i,g)
            if (i==3 .and. g==127) F = F + lik(id(obs))
          END IF
 
        ELSE IF (obs == 85391) THEN
          IF (i==1 .OR. i==2) THEN
             prob(obs,i,g) = x(46+i)
          ELSE
              prob(obs,i,g) = 1 - x(47) - x(48)
          END IF
          weight(obs,i,g)=prob(obs,i,g)*x(20)**good(1)*(1-x(20))**(1-good(1))*x(20)**good(2)*(1-x(20))**(1-good(2))*&
					 &x(20)**good(3)*(1-x(20))**(1-good(3))*x(20)**good(4)*(1-x(20))**(1-good(4))*&
					 &x(20)**good(5)*(1-x(20))**(1-good(5))*x(20)**good(6)*(1-x(20))**(1-good(6))*&
					 &x(20)**good(7)*(1-x(20))**(1-good(7))
          lik(id(obs)) = lik(id(obs))+weight(obs,i,g)*lprod(obs,i,g)
          if (i==3 .and. g==127) F = F + lik(id(obs))
        END IF

  END DO
  END DO
  END DO

END SUBROUTINE lik_fcn



