!**********************************************************************************
!
!	SUBROUTINE:  DecRule
!
!	PURPOSE: After policy function has converged, calculates decision rules.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 8/6/2007
!**********************************************************************************

subroutine DecRule
  
  ! inputs:
  !
  ! output:
  ! 
  ! remarks:
  !     Calculates decision rules in 4 directions of interest while holding other 3
  !      constant at their steady-state value.  4 directions are own and other 
  !      country capital stock and own and other country tech shock.
  !      All decision rules are computed from the perspective of country 1.
  !

  use Params
  
  implicit none
  
  !Variable Declaration
  integer :: h,i,j
  real(prec),dimension(d) :: DecPoint
  real(prec) :: NewDecPoint
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  
  real(prec),dimension(nk) :: DecOwnCap, DecOwnCons, DecOthCap, DecOthCons
  real(prec),dimension(nz) :: DecOwnSCap, DecOwnSCons, DecOthSCap, DecOthSCons
  real(prec),dimension(nk) :: DecOwnLab, DecOthLab
  real(prec),dimension(nz) :: DecOwnSLab, DecOthSLab
  
  !Returns policy choice (either consumption or capital) given state
  real(prec), external :: ChebValue
  
  !Construct Decision Rule over 100 evenly spaced points, varying own-country capital
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do
  
  do i = 1,nk	
     DecPoint(1) = knorm(i,1)	
     
     !Constructs policy functions based on Coeffs (both Cons and Cap) 
     if (i == 1) then
        call ChebPoly(d, DecPoint, Ttemp)
     elseif (i > 1) then
        !Convert DecPoint(1) to [-1,1] grid
        NewDecPoint = 2.0 * (DecPoint(1) - kmin)/(kmax-kmin) - 1.0
        
        if (NewDecPoint < -1) then
           print*, 'error: capital below kmin in DecRule'
           pause
        elseif (NewDecPoint > 1) then
           print*, 'error: capital above kmax in DecRule'
           pause
        end if
        
        !Construct Chebyshev Polynomials at newstate
        Ttemp(1,1) = 1
        Ttemp(2,1) = NewDecPoint
        do h = 3,5	!Since q = d+2, 2**(q-d)+1 = 5
           Ttemp(h,1) = 2.0 * NewDecPoint * Ttemp(h-1,1) - Ttemp(h-2,1)
        end do
     endif
     
     DecOwnCap(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOwnCons(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOwnLab(i) = ChebValue(Ttemp,CoeffsLab1(:,:,1),CoeffsLab2(:,:,:,:,1),CoeffsLab3(:,:,1),YLab0(1))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying other-country capital
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do
  
  do i = 1,nk	
     
     DecPoint(N) = knorm(i,1)	
     
     !Constructs policy functions based on Coeffs (both Cons and Cap) 
     if (i == 1) then
        call ChebPoly(d, DecPoint, Ttemp)
     elseif (i > 1) then
        !Convert DecPoint(N) to [-1,1] grid
        NewDecPoint = 2.0 * (DecPoint(N) - kmin)/(kmax-kmin) - 1.0
        
        if (NewDecPoint < -1) then
           print*, 'error: capital below kmin in DecRule'
           pause
        elseif (NewDecPoint > 1) then
           print*, 'error: capital above kmax in DecRule'
           pause
        end if
        
        !Construct Chebyshev Polynomials at newstate
        Ttemp(1,N) = 1
        Ttemp(2,N) = NewDecPoint
        do h = 3,5			!Since q = d+2, 2**(q-d)+1 = 5
           Ttemp(h,N) = 2.0 * NewDecPoint * Ttemp(h-1,N) - Ttemp(h-2,N)
        end do
     endif
     
     DecOthCap(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOthCons(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOthLab(i) = ChebValue(Ttemp,CoeffsLab1(:,:,1),CoeffsLab2(:,:,:,:,1),CoeffsLab3(:,:,1),YLab0(1))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying own-country tech shock
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do

  do i = 1,nz	
     
     DecPoint(N+1) = znorm(i,1)	
     
     !Constructs policy functions based on Coeffs (both Cons and Cap) 
     if (i == 1) then
        call ChebPoly(d, DecPoint, Ttemp)
     elseif (i > 1) then
        !Convert DecPoint(N+1) to [-1,1] grid
        NewDecPoint = 2.0 * (DecPoint(N+1) - zmin)/(zmax-zmin) - 1.0
        
        if (NewDecPoint < -1) then
           print*, 'error: tech shock below zmin in DecRule'
           pause
        elseif (NewDecPoint > 1) then
           print*, 'error: tech shock above zmax in DecRule'
           pause
        end if
        
        !Construct Chebyshev Polynomials at newstate
        Ttemp(1,N+1) = 1
        Ttemp(2,N+1) = NewDecPoint
        do h = 3,5			!Since q = d+2, 2**(q-d)+1 = 5
           Ttemp(h,N+1) = 2.0 * NewDecPoint * Ttemp(h-1,N+1) - Ttemp(h-2,N+1)
        end do
     endif
     
     DecOwnSCap(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOwnSCons(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOwnSLab(i) = ChebValue(Ttemp,CoeffsLab1(:,:,1),CoeffsLab2(:,:,:,:,1),CoeffsLab3(:,:,1),YLab0(1))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying other-country tech shock
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do
  
  do i = 1,nz	
     
     DecPoint(N+N) = znorm(i,1)	
     
     !Constructs policy functions based on Coeffs (both Cons and Cap) 
     if (i == 1) then
        call ChebPoly(d, DecPoint, Ttemp)
     elseif (i > 1) then
        !Convert DecPoint(N+N) to [-1,1] grid
        NewDecPoint = 2.0 * (DecPoint(N+N) - zmin)/(zmax-zmin) - 1.0
        
        if (NewDecPoint < -1) then
           print*, 'error: tech shock below zmin in DecRule'
           pause
        elseif (NewDecPoint > 1) then
           print*, 'error: tech shock above zmax in DecRule'
           pause
        end if
        
        !Construct Chebyshev Polynomials at newstate
        Ttemp(1,N+N) = 1
        Ttemp(2,N+N) = NewDecPoint
        do h = 3,5			!Since q = d+2, 2**(q-d)+1 = 5
           Ttemp(h,N+N) = 2.0 * NewDecPoint * Ttemp(h-1,N+N) - Ttemp(h-2,N+N)
        end do
     endif
     
     DecOthSCap(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOthSCons(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOthSLab(i) = ChebValue(Ttemp,CoeffsLab1(:,:,1),CoeffsLab2(:,:,:,:,1),CoeffsLab3(:,:,1),YLab0(1))
  enddo
  
  !Save decision rules to be output to file
  VarOut(:,1) = DecOwnCap
  VarOut(:,2) = DecOthCap
  VarOut(:,3) = DecOwnSCap
  VarOut(:,4) = DecOthSCap
  VarOut(:,5) = DecOwnCons
  VarOut(:,6) = DecOthCons
  VarOut(:,7) = DecOwnSCons
  VarOut(:,8) = DecOthSCons
  VarOut(:,9) = DecOwnLab
  VarOut(:,10) = DecOthLab
  VarOut(:,11) = DecOwnSLab
  VarOut(:,12) = DecOthSLab
  VarOut(:,13) = knorm(:,1)  !Assumes all countries have same grid
  VarOut(:,14) = znorm(:,1) !Assumes above and that nz = nk

end subroutine DecRule
