!**********************************************************************************
!
!	SUBROUTINE:  DecRule
!
!	PURPOSE: After policy function has converged, calculates decision rules.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 8/3/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 countries 1 & N.
  !

  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,N) :: DecOwnCap, DecOwnCons, DecOthCap, DecOthCons
  real(prec),dimension(nz,N) :: DecOwnSCap, DecOwnSCons, DecOthSCap, DecOthSCons
  real(prec),dimension(nk,N) :: DecOwnLab, DecOthLab
  real(prec),dimension(nz,N) :: DecOwnSLab, DecOthSLab
  
  !Returns policy choice (either consumption or capital) given state
  real(prec), external :: ChebValue
  
  !Construct Decision Rule over 100 evenly spaced points, varying country-1 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(1))/(kmax(1)-kmin(1)) - 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,1) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOwnCons(i,1) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOwnLab(i,1) = (Lfoc(1) * exp(DecPoint(N+1)) * DecPoint(1)**alpha &
          & * DecOwnCons(i,1)**(-xgamma(1)))**(Lelast(1))
     DecOthCap(i,N) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     DecOthCons(i,N) = DecOwnCons(i,1)**(gamma(N)/gamma(1)) * (Pareto(N)/Pareto(1))**gamma(N)
     DecOthLab(i,N) = (Lfoc(N) * exp(DecPoint(N+N)) * DecPoint(N)**alpha &
          & * DecOthCons(i,N)**(-xgamma(N)))**(Lelast(N))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying country-N capital
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do
  
  do i = 1,nk	
     
     DecPoint(N) = knorm(i,N)	
     
     !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(N))/(kmax(N)-kmin(N)) - 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,1) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOthCons(i,1) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOthLab(i,1) = (Lfoc(1) * exp(DecPoint(N+1)) * DecPoint(1)**alpha &
          & * DecOthCons(i,1)**(-xgamma(1)))**(Lelast(1))
     DecOwnCap(i,N) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     DecOwnCons(i,N) = DecOthCons(i,1)**(gamma(N)/gamma(1)) * (Pareto(N)/Pareto(1))**gamma(N)
     DecOwnLab(i,N) = (Lfoc(N) * exp(DecPoint(N+N)) * DecPoint(N)**alpha &
          & * DecOwnCons(i,N)**(-xgamma(N)))**(Lelast(N))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying country-1 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,1) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOwnSCons(i,1) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOwnSLab(i,1) = (Lfoc(1) * exp(DecPoint(N+1)) * DecPoint(1)**alpha &
          & * DecOwnSCons(i,1)**(-xgamma(1)))**(Lelast(1))
     DecOthSCap(i,N) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     DecOthSCons(i,N) = DecOwnSCons(i,1)**(gamma(N)/gamma(1)) * (Pareto(N)/Pareto(1))**gamma(N)
     DecOthSLab(i,N) = (Lfoc(N) * exp(DecPoint(N+N)) * DecPoint(N)**alpha &
          & * DecOthSCons(i,N)**(-xgamma(N)))**(Lelast(N))
  enddo
  
  !Construct Decision Rule over 100 evenly spaced points, varying country-N tech shock
  do j = 1,N
     DecPoint(j) = kss
     DecPoint(N+j) = 0
  end do
  
  do i = 1,nz	
     
     DecPoint(N+N) = znorm(i,N)	
     
     !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,1) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     DecOthSCons(i,1) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     DecOthSLab(i,1) = (Lfoc(1) * exp(DecPoint(N+1)) * DecPoint(1)**alpha &
          & * DecOthSCons(i,1)**(-xgamma(1)))**(Lelast(1))
     DecOwnSCap(i,N) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     DecOwnSCons(i,N) = DecOthSCons(i,1)**(gamma(N)/gamma(1)) * (Pareto(N)/Pareto(1))**gamma(N)
     DecOwnSLab(i,N) = (Lfoc(N) * exp(DecPoint(N+N)) * DecPoint(N)**alpha &
          & * DecOwnSCons(i,N)**(-xgamma(N)))**(Lelast(N))
  enddo
  
  !Save decision rules to be output to file
  VarOut(:,1) = DecOwnCap(:,1)
  VarOut(:,2) = DecOthCap(:,1)
  VarOut(:,3) = DecOwnSCap(:,1)
  VarOut(:,4) = DecOthSCap(:,1)
  VarOut(:,5) = DecOwnCons(:,1)
  VarOut(:,6) = DecOthCons(:,1)
  VarOut(:,7) = DecOwnSCons(:,1)
  VarOut(:,8) = DecOthSCons(:,1)
  VarOut(:,9) = DecOwnLab(:,1)
  VarOut(:,10) = DecOthLab(:,1)
  VarOut(:,11) = DecOwnSLab(:,1)
  VarOut(:,12) = DecOthSLab(:,1)
  VarOut(:,13) = knorm(:,1) 
  VarOut(:,14) = znorm(:,1) !Assumes countries have same z-grid and that nz = nk
  VarOut(:,15) = DecOwnCap(:,N)
  VarOut(:,16) = DecOthCap(:,N)
  VarOut(:,17) = DecOwnSCap(:,N)
  VarOut(:,18) = DecOthSCap(:,N)
  VarOut(:,19) = DecOwnCons(:,N)
  VarOut(:,20) = DecOthCons(:,N)
  VarOut(:,21) = DecOwnSCons(:,N)
  VarOut(:,22) = DecOthSCons(:,N)
  VarOut(:,23) = DecOwnLab(:,N)
  VarOut(:,24) = DecOthLab(:,N)
  VarOut(:,25) = DecOwnSLab(:,N)
  VarOut(:,26) = DecOthSLab(:,N)
  VarOut(:,27) = knorm(:,N) 

end subroutine DecRule
