!   Last change:  SKP  26 Nov 2001   10:11 pm
!   This program finds the optimal replacement ratio for T=1


!-----------------

module calibration

REAL :: y, h, beta, pie, piu
! Note: sigma and rho are defined in the function uf(a,b)
! This is the no moral hazard case.
parameter (y=1.0_8, h=0.45_8, beta=0.995_8, pie=0.0, piu=0.0)


INTEGER :: step, bt
parameter (step=300, bt=1)

REAL :: er_crit1, er_crit2, er_crit3, bc_err

parameter (er_crit1=0.01_8, er_crit11=0.000000001_8, er_crit2=0.000000000000001_8, er_crit3=0.0000001_8, bc_err=0.000000001_8)

end module

!----------------

module global
save
  INTEGER, PARAMETER ::tstep=300, sont=1,npop=20
  real :: m(0:tstep+1), v(0:tstep,0:(sont+1),0:1), tv(0:tstep,0:(sont+1),0:1)
  real :: i1(0:tstep,0:sont), i2(0:tstep,0:sont), i3(0:tstep,0:sont)
  INTEGER :: indextv(0:tstep,0:sont), indexti1(0:tstep,0:sont),sigma
  INTEGER :: indexti2(0:tstep,0:sont), indexti3(0:tstep,0:sont), ne(0:tstep,0:sont)
  INTEGER :: optm(0:tstep,0:sont,1:3), optn(0:tstep,0:sont,1:3)
  real :: x(0:1,0:1), tau, temptau, theta(0:sont), pi(0:sont)
  real :: lamda(0:tstep,0:sont,3),evolve(1:npop,0:sont+5)
  real :: erate(0:sont), unerate(0:sont), bc , avut, avas, stdc, emprate, tau_l, tau_h
  REAL:: maxj
end module

!----------------

real function uf(a,b)
implicit none
REAL sigma, rho
parameter (sigma=0.67, rho=2.5)
REAL, INTENT(IN) :: a,b

uf=  (((a**(1-sigma))*(b**sigma))**(1-rho)-1)/(1-rho)
end

!----------------


subroutine first_guess
USE global
INTEGER:: ibu,n1


evolve=0
do n1=1,20
    evolve(n1,0) = 0.05*(n1-1);
    evolve(n1,1) = 0.05*(n1-1);
end do

do ibu=1,10
   PRINT*,'evolve=',evolve(ibu,:)
end do

end subroutine

!----------------

subroutine initialize
USE global
USE calibration
implicit none

INTEGER :: i,j,l

maxj=0

  x(1,1)=0.9681_8
  x(1,0)=0.0319_8
  x(0,0)=0.5_8
  x(0,1)=0.5_8

  m(0)=0.0_8
  m(step+1)=100.0_8
  do i=1, (step)
   m(i)=m(i-1)+(8.0/REAL(step))
  end do

  !initialize value function

  do i=0,1
    do j=0,(sont+1)
      do l=0,tstep
        v(l,j,i)=-(500)+l
      end do
    end do
  end do

 ! initialize prob. dbn.
  lamda=(1.0)/( REAL( (step+1)*(bt+1)*3))

  pi(0)=pie
  pi(1:bt)=piu
end subroutine


!----------------



! this subroutine makes only one iteration of the value function
! and finds decision rules

subroutine valuemap
USE calibration
USE global
implicit none
INTEGER :: i, j, t
real :: tu, tu1, ti, ti1, uf


! for each t & i, find maximizer and update the value function
  do t=0,bt
      !for a given t, find maximizers for all m(i)
      tu1=tu
      ti1=ti
      ! s=u case
      ! for m(0), find maximizer
      j=0
      if ( ((1-tau)*theta(t)*y) .gt. 0.0_8) then
         tu=uf( ((1-tau)*theta(t)*y),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
         j=j+1
      else
         tu=-10000000.0_8
         j=j+1
      endif

      if ( (m(0)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) then
         tu1=uf( (m(0)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
      end if

      !!!!!!!!
do1:  do while ( (tu1 .ge. tu) .and. ( (m(0)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) )
         tu1=uf( (m(0)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
         if (tu1 .ge. tu) then
            tu=tu1
            j=j+1
         else
             EXIT do1
         endif
      enddo do1

      j=j-1
      indextv(0,t)=j
      tv(0,t,0)=tu


      !for m(i) find maximizer
      do i=1, step
         !!
         tu=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
         j=j+1

         if ( (m(i)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) then
            tu1=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
         end if
         !!!!!
      !   tu1=tu+1
do2:     do while ( (tu1 .ge. tu) .and. ( (m(i)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) )
            tu1=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(0,0)*v(j,t+1,0)+x(0,1)*v(j,t+1,1) )
            if (tu1 .ge. tu) then
               tu=tu1
               j=j+1
            else
                EXIT do2
            endif
         enddo do2
           
         j=j-1
         indextv(i,t)=j
         tv(i,t,0)=tu

      enddo


      !s=e case, if the worker accepts the job offer
      ! for each m(i), find maximizer
      j=0
      do i=0, step
         !!
         ti=uf( (m(i)+(1-tau)*y-m(j)),(1-h))+ beta*( x(1,0)*v(j,0,0)+x(1,1)*v(j,0,1) )
         j=j+1
         if (m(i)+(1-tau)*y-m(j) .GT. 0.0_8 ) then
            ti1=uf( (m(i)+(1-tau)*y-m(j)),(1-h))+ beta*( x(1,0)*v(j,0,0)+x(1,1)*v(j,0,1) )
         end if
         !!!!!!
         !ti1=ti+1

do3:     do while ( (ti1 .ge. ti) .and. ( (m(i)+ (1-tau)*y-m(j)) .gt. 0.0_8) )
            ti1=uf( (m(i)+(1-tau)*y-m(j)),(1-h))+ beta*( x(1,0)*v(j,0,0)+x(1,1)*v(j,0,1) )
            if (ti1 .ge. ti) then
               ti=ti1
               j=j+1
            else
                EXIT do3
            endif
         enddo do3
         j=j-1
         indexti1(i,t)=j
         i1(i,t)=ti
      enddo

      ! s=e case, if the worker rejects job offer, and gets insurance
      ! for m(0) find maximizer
      j=0
      if ( ((1-tau)*theta(t)*y) .gt. 0.0_8) then
         ti=uf( ((1-tau)*theta(t)*y),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         j=j+1
      else
         ti=-10000000.0_8
         j=j+1
      endif

      if ( (m(0)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) then
         ti1=uf( (m(0)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
      end if
      !!!!!!
      !ti1=ti+1
do4:  do while ( (ti1 .ge. ti) .and. ( (m(0)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) )
         ti1=uf( (m(0)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         if (ti1 .ge. ti) then
            ti=ti1
            j=j+1
         else
             EXIT do4
         endif   
      enddo  do4
      j=j-1
      indexti2(0,t)=j
      i2(0,t)=ti

      ! for each m(i) find maximizer
      do i=1, step
         !!
         ti=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         j=j+1
         if ( (m(i)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) then
            ti1=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         end if
         !ti1=ti+1
do5:     do while ( (ti1 .ge. ti) .and. ( (m(i)+ (1-tau)*theta(t)*y-m(j)) .gt. 0.0_8) )
            ti1=uf( (m(i)+(1-tau)*theta(t)*y-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
            if (ti1 .ge. ti) then
               ti=ti1
               j=j+1
            else
                EXIT do5
            endif       
         enddo do5
         j=j-1
         indexti2(i,t)=j
         i2(i,t)=ti
      enddo
         
      ! if the worker rejects job offer and does not get any insurance

      i3(0,t)=-10000000.0_8
      indexti3(0,t)=0
      j=0
      
      do i=1, step

         ti=uf( (m(i)-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         j=j+1
         if ((m(i)-m(j)) .gt. 0.0_8) then
            ti1= uf( (m(i)-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
         end if

do6:     do while ( (ti1 .ge. ti) .and. ( (m(i)-m(j)) .gt. 0.0_8) )
            ti1=uf( (m(i)-m(j)),1.0)+ beta*( x(1,0)*v(j,t+1,0)+x(1,1)*v(j,t+1,1) )
            if (ti1 .ge. ti) then
               ti=ti1
               j=j+1
            else
               EXIT do6
            endif     
         enddo do6
         j=j-1
         indexti3(i,t)=j
         i3(i,t)=ti
      enddo
   END do ! end of loop which changes 't'


end subroutine


!----------------



! This subroutine makes value function iteration until a certain convergence
! criteria is satisfied

subroutine vconverge
USE calibration
USE global
implicit none
INTEGER :: i, t, sayac, maxsayac
REAL ::  error1

sayac=0
maxsayac=2000
optm=1
error1=20


do while ((error1 .gt. er_crit1) .AND. (sayac<=maxsayac))
   sayac=sayac+1
   call valuemap

   do t=0,bt
      do i=0, step
        if (i1(i,t) .gt. ( pi(t)*i2(i,t)+ (1-pi(t))*i3(i,t) )    ) then 
           tv(i,t,1)=i1(i,t)
           ne(i,t)=1
        else
           tv(i,t,1)=pi(t)*i2(i,t)+ (1-pi(t))*i3(i,t)
           ne(i,t)=0
        endif
      enddo

   end do
   tv(:,bt+1,:)=tv(:,bt,:)

   error1=SUM( (v-tv)*(v-tv))
   v=tv

enddo  ! Value fnc iteration is over.

     do t=0,bt
      optn(:,t,1)=t+1
      optn(:,t,2)=(1-ne(:,t))*(t+1)
      optn(:,t,3)=t+1
     END do
     optm(:,:,1)=indextv
     optm(:,:,2)=ne*indexti1+(1-ne)*indexti3
     optm(:,:,3)=indexti2

    ! update for the last period
    optn(:,bt,1)=bt
    optn(:,bt,2)=(1-ne(:,bt))*bt
    optn(:,bt,3)=bt

PRINT*,sayac,'th iteration---vconverge is over.'
end subroutine



! This subroutine makes value function iteration until a certain convergence
! criteria is satisfied

subroutine vconverge2
USE calibration
USE global
implicit none
INTEGER :: i, t , sayac, maxsayac
REAL ::  error11

optm=1
error11=20
sayac=0
maxsayac=10000

do while ((error11 .gt. er_crit11) .AND. (sayac<=maxsayac))

sayac=sayac+1

   call valuemap
   do t=0,bt
      do i=0, step
        if (i1(i,t) .gt. ( pi(t)*i2(i,t)+ (1-pi(t))*i3(i,t) )    ) then 
           tv(i,t,1)=i1(i,t)
           ne(i,t)=1
        else
           tv(i,t,1)=pi(t)*i2(i,t)+ (1-pi(t))*i3(i,t)
           ne(i,t)=0
        endif
      enddo
   end do
   tv(:,bt+1,:)=tv(:,bt,:)

   !!!!!!!!!!!! VALUE FUNCTION CONVERGENCE
   error11=SUM( (v-tv)**2 )
   v=tv

enddo  ! Value fnc iterationis over

     do t=0,bt
      optn(:,t,1)=t+1
      optn(:,t,2)=(1-ne(:,t))*(t+1)
      optn(:,t,3)=t+1
     END do
     optm(:,:,1)=indextv
     optm(:,:,2)=ne*indexti1+(1-ne)*indexti3
     optm(:,:,3)=indexti2

    ! update for the last period
    optn(:,bt,1)=bt
    optn(:,bt,2)=(1-ne(:,bt))*bt
    optn(:,bt,3)=bt

PRINT*, sayac,'th iteration---vconverge is over'

end subroutine


!---------------------


! This subroutine makes convergence of probability distributions

subroutine probdbn
USE global
USE calibration
implicit none
INTEGER :: k,t,i, tm, tn, tn2 ,sayac, maxsayac
REAL :: tlamda(0:step,0:bt,1:3), error2

  sayac=1
  maxsayac=10000
  tlamda=0.0_8
  error2=20
  do while  ((error2 .gt. er_crit2) .AND. (sayac<=maxsayac))
  sayac=sayac+1
      k=1
      do t=0, bt
         do i=0, step
            tm=optm(i,t,k)
            tn=optn(i,t,k)

            ! sprime=u mprime=1
            tlamda(tm,tn,1)=tlamda(tm,tn,1)+lamda(i,t,k)*x(0,0)

            ! sprime=e mprime=0
            tn2=ne(tm,tn)
            tlamda(tm,tn,2)=tlamda(tm,tn,2)+lamda(i,t,k)*x(0,1)*( tn2+(1-tn2)*(1-pi(tn)) )

            ! sprime=e mprime=1
            tlamda(tm,tn,3)=tlamda(tm,tn,3)+lamda(i,t,k)*x(0,1)* ( (1-tn2)* pi(tn) )
         enddo
      enddo
      

      do k=2, 3
        do t=0, bt
           do i=0, step
              tm=optm(i,t,k)
              tn=optn(i,t,k)

              ! sprime=u mprime=1
              tlamda(tm,tn,1)=tlamda(tm,tn,1)+lamda(i,t,k)*x(1,0)

              ! sprime=e mprime=0
              tn2=ne(tm,tn)
              tlamda(tm,tn,2)=tlamda(tm,tn,2)+lamda(i,t,k)*x(1,1)*(tn2+(1-tn2)*(1-pi(tn)))

              ! sprime=e mprime=1
              tlamda(tm,tn,3)=tlamda(tm,tn,3)+lamda(i,t,k)*x(1,1)* ( (1-tn2)* pi(tn) )
           enddo
        enddo     
      enddo     

      error2=SUM((lamda-tlamda)*(lamda-tlamda))
      lamda=tlamda
      tlamda=0.0_8



  enddo
PRINT*,'probdbn iteration=',sayac

end subroutine


!----------------



subroutine avstat
USE global
USE calibration
implicit none
INTEGER :: t,i
REAL :: c(0:step,0:bt,3), meanc, meanc2

  avut=(1-beta)* (SUM( lamda(:,:,1)*v(:,0:bt,0) )+ SUM( lamda(:,:,2)* (ne*i1+(1-ne)*i3) )+ SUM( lamda(:,:,3)*i2 ))

  avas=0
  do t=0,bt
     avas=avas+SUM(m(0:step)*( lamda(:,t,1)+lamda(:,t,2)+lamda(:,t,3) ) )
  end do

  meanc=0
  meanc2=0
  do i=0,step
     do t=0,bt
        c(i,t,1)=m(i)-m(optm(i,t,1))+(1-tau)*theta(t)*y
        c(i,t,2)=m(i)-m(optm(i,t,2))+ne(i,t)*(1-tau)*y
        c(i,t,3)=m(i)-m(optm(i,t,3))+(1-tau)*theta(t)*y
        meanc=meanc+c(i,t,1)*lamda(i,t,1)+c(i,t,2)*lamda(i,t,2)+ c(i,t,3)*lamda(i,t,3)
        meanc2=meanc2+ (c(i,t,1)**2)*lamda(i,t,1)+ (c(i,t,2)**2)*lamda(i,t,2)+  (c(i,t,3)**2)*lamda(i,t,3)
     end do
  end do

  stdc=(meanc2- (meanc**2) )

  do t=0,bt
     erate(t)=sum(lamda(:,t,2)*ne(:,t))
  end do
  emprate=SUM(erate)

  ! This part writes the average consumption to a file

end subroutine



subroutine updatetau
USE global
USE calibration
implicit none
INTEGER :: t

  ! BC computation
  do t=0,bt
     erate(t)=sum(lamda(:,t,2)*ne(:,t))
     unerate(t)=SUM(lamda(:,t,1)+lamda(:,t,3))
  end do
  bc=(1-tau)*y* SUM(theta*unerate)-SUM(erate)*y*tau
  ! BC computed
if (ABS(bc)>bc_err) then

  if (bc>0) then
      tau_l=tau
  else
      tau_h=tau
  end if

  ! update tau
  tau= (y*SUM(theta*unerate))/(SUM(erate)*y+y*SUM(theta*unerate) )
  if ((tau<=tau_l) .or. (tau>=tau_h)) then
     tau=(tau_h+tau_l)/2
  end if
end if
end subroutine


!-------------------

subroutine sort
USE global
INTEGER :: is,idx(1:npop),tempidx
REAL :: tempavut(1:npop),swap,tempevolve(1:npop,0:sont+5)


do is=1,npop
  idx(is)=is
end do

tempavut=evolve(:,sont+2)
swap=1
do WHILE (swap > 0.0)
  swap=0.0
  do is=1,npop-1
    if (tempavut(idx(is))<tempavut(idx(is+1))) then
      swap=1.0
      tempidx=idx(is+1)
      idx(is+1)=idx(is)
      idx(is)=tempidx
    end if
  end do
end do
do is=1,npop
  tempevolve(is,:)=evolve(idx(is),:)
end do
evolve=tempevolve
end subroutine

!------------------------

subroutine evolve_theta2
USE global
INTEGER :: is,iss
REAL :: temptheta,rr


do is=npop/2+1,npop

  do iss=0,sont
    call RANDOM_NUMBER(rr)
    rr=rr-0.5
    temptheta=evolve(is-npop/2,iss)+INT(sigma*rr)*0.05
    rr=MAX(0.0,temptheta)
    evolve(is,iss)=MIN(1.0,rr)
  if (evolve(is,iss)<0.05) then
    evolve(is,iss)=0.0
  end if

  end do
end do
end subroutine

!---------------------

program optimal
USE calibration
USE global
INTEGER :: id,hoho

call first_guess
sigma=7


tau=0.0252
do id=1,npop
  theta=evolve(id,0:sont)
  tau_l=0
  tau_h=1
  call initialize

  PRINT*,id,'th iteration.....'

  PRINT*,'vconverge'
  call vconverge
  PRINT*,'probdbn'
  call probdbn

  tau_l=0
  tau_h=1
  bc=20



  do WHILE (((tau_h-tau_l)>0.00001_8) .AND. (ABS(bc) > bc_err))

   PRINT*,id,'th iteration'
   call updatetau
   PRINT*,'tau=',tau
   PRINT*,'vconverge2'
   call vconverge2
   PRINT*,'probdbn'
   call probdbn
   PRINT*,'bc=',bc
   PRINT*,'difference=',tau_h-tau_l
  END do

  call avstat
  PRINT*,'avut=',avut
  evolve(id,sont+1)=tau
  evolve(id,sont+2)=avut
  evolve(id,sont+3)=emprate
  evolve(id,sont+4)=stdc
  evolve(id,sont+5)=avas

  PRINT*,'The best so far:'
  do hoho=1,20
    PRINT*,hoho,'>',evolve(hoho,:)
  end do

  open(unit=7,file="Temp",status="replace")
  write (unit=7, fmt=*) 'The best so far:'
  do hoho=1,20
        write (unit=7, fmt=*) hoho,'>',evolve(hoho,:)
        write (unit=7, fmt=*) 'avut=',evolve(hoho,bt+2)
        write (unit=7, fmt=*) 'tax=',evolve(hoho,bt+1)
        write (UNIT=7, FMT=*) 'employment rate=', evolve(hoho,bt+3)
        write (UNIT=7, FMT=*) 'variance of consumption=',  evolve(hoho,bt+4)
        write (UNIT=7, FMT=*) 'average asset holdings=',evolve(hoho,bt+5)
  end do
  close(unit=7)

end do

call sort
PRINT*,'The best:'
do hoho=1,20
  PRINT*,hoho,'>',evolve(hoho,:)
end do

  open(unit=7,file="Result",status="replace")
  write (unit=7, fmt=*) 'The sorted replacement ratios:'
  do hoho=1,20
        write (unit=7, fmt=*) hoho,'>',evolve(hoho,:)
        write (unit=7, fmt=*) 'avut=',evolve(hoho,bt+2)
        write (unit=7, fmt=*) 'tax=',evolve(hoho,bt+1)
        write (UNIT=7, FMT=*) 'employment rate=', evolve(hoho,bt+3)
        write (UNIT=7, FMT=*) 'variance of consumption=',  evolve(hoho,bt+4)
        write (UNIT=7, FMT=*) 'average asset holdings=',evolve(hoho,bt+5)
  end do
  close(unit=7)
end program