!**********************************************************************************
!
!	SUBROUTINE:  DecRule2
!
!	PURPOSE: Calculates decision rules for comparison project.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 12/03/09
!**********************************************************************************

subroutine DecRule2
  
  ! 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.
  !      The decision rules are computed from the perspective of country 1 and country N.
  !

  use Params
  use StaticFOC
    
  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
  
  !Country1
  real(prec),dimension(nk) :: OwnCap1, OwnCons1, OwnInv1, OwnLab1
  real(prec),dimension(nk) :: OthCap1, OthCons1, OthInv1, OthLab1
  real(prec),dimension(nz) :: OwnSCap1, OwnSCons1, OwnSInv1, OwnSLab1
  real(prec),dimension(nz) :: OthSCap1, OthSCons1, OthSInv1, OthSLab1

  !Country N
  real(prec),dimension(nk) :: OwnCapN, OwnConsN, OwnInvN, OwnLabN
  real(prec),dimension(nk) :: OthCapN, OthConsN, OthInvN, OthLabN
  real(prec),dimension(nz) :: OwnSCapN, OwnSConsN, OwnSInvN, OwnSLabN
  real(prec),dimension(nz) :: OthSCapN, OthSConsN, OthSInvN, OthSLabN
   
  !Returns policy choice (either consumption or capital) given state
  real(prec), external :: ChebValue
  
  !Variables used in nonlinear equation solver to solve for country 1's labor supply
  integer, parameter :: itmaxA = 300
  integer, parameter :: countA = 1
  real(prec) :: RelerrA, fnormA
  real(prec), dimension(countA) :: guessA, finalA
  
  !Variables used in nonlinear equation solver to solve for cons and labor of countries 2-N
  integer, parameter :: itmaxB = 300
  integer, parameter :: countB = 2
  real(prec) :: RelerrB, fnormB
  real(prec), dimension(countB) :: guessB, finalB
  
  !Initialize inputs to nle solver, DNEQNF
  RelerrA = 0.0001
  fnormA = 0
  guessA = 1.0
  finalA = 1.0
  
  !Initialize inputs to nle solver, DNEQNF
  RelerrB = 0.0001
  fnormB = 0
  guessB(1) = A_tfp
  guessB(2) = 1.0
  finalB = 1.0

  !Construct Decision Rule over nk 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) = knorm2(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
        
        !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
     
     OwnCap1(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     OwnCons1(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     OwnInv1(i) = OwnCap1(i) - (1.0 - delta)*DecPoint(1)

     guessA = lss
     DecCons1 = OwnCons1(i)
     DecCap1 = DecPoint(1)
     DecShock1 = DecPoint(N+1)
     call DNEQNF(DConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
     OwnLab1(i) = finalA(1)

     OthCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OthInvN(i) = OthCapN(i) - (1.0 - delta)*DecPoint(N)

     DecLab1 = OwnLab1(i)
     DecCapN = DecPoint(N)
     DecShockN = DecPoint(N+N)
     guessB(1) = log(OwnCons1(i))
     guessB(2) = OwnLab1(i)
     call DNEQNF(DConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
     OthConsN(i) = exp(finalB(1))
     OthLabN(i) = finalB(2)
  enddo
  
  !Construct Decision Rule over nk 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) = knorm2(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
        
        !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
     
     OthCap1(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     OthCons1(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     OthInv1(i) = OthCap1(i) - (1.0 - delta)*DecPoint(1)

     guessA = lss
     DecCons1 = OthCons1(i)
     DecCap1 = DecPoint(1)
     DecShock1 = DecPoint(N+1)
     call DNEQNF(DConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
     OthLab1(i) = finalA(1)

     OwnCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OwnInvN(i) = OwnCapN(i) - (1.0 - delta)*DecPoint(N)

     DecLab1 = OthLab1(i)
     DecCapN = DecPoint(N)
     DecShockN = DecPoint(N+N)
     guessB(1) = log(OthCons1(i))
     guessB(2) = OthLab1(i)
     call DNEQNF(DConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
     OwnConsN(i) = exp(finalB(1))
     OwnLabN(i) = finalB(2)
  enddo
  
  !Construct Decision Rule over nz 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) = znorm2(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
        
        !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
     
     OwnSCap1(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     OwnSCons1(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     OwnSInv1(i) = OwnSCap1(i) - (1.0 - delta)*DecPoint(1)

     guessA = lss
     DecCons1 = OwnSCons1(i)
     DecCap1 = DecPoint(1)
     DecShock1 = DecPoint(N+1)
     call DNEQNF(DConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
     OwnSLab1(i) = finalA(1)

     OthSCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OthSInvN(i) = OthSCapN(i) - (1.0 - delta)*DecPoint(N)

     DecLab1 = OwnSLab1(i)
     DecCapN = DecPoint(N)
     DecShockN = DecPoint(N+N)
     guessB(1) = log(OwnSCons1(i))
     guessB(2) = OwnSLab1(i)
     call DNEQNF(DConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
     OthSConsN(i) = exp(finalB(1))
     OthSLabN(i) = finalB(2)
  enddo
  
  !Construct Decision Rule over nz 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) = znorm2(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
        
        !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
     
     OthSCap1(i) = ChebValue(Ttemp,CoeffsCap1(:,:,1),CoeffsCap2(:,:,:,:,1),CoeffsCap3(:,:,1),YCap0(1))
     OthSCons1(i) = ChebValue(Ttemp,CoeffsCons1(:,:,1),CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     OthSInv1(i) = OthSCap1(i) - (1.0 - delta)*DecPoint(1)

     guessA = lss
     DecCons1 = OthSCons1(i)
     DecCap1 = DecPoint(1)
     DecShock1 = DecPoint(N+1)
     call DNEQNF(DConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
     OthSLab1(i) = finalA(1)

     OwnSCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OwnSInvN(i) = OwnSCapN(i) - (1.0 - delta)*DecPoint(N)

     DecLab1 = OthSLab1(i)
     DecCapN = DecPoint(N)
     DecShockN = DecPoint(N+N)
     guessB(1) = log(OthSCons1(i))
     guessB(2) = OthSLab1(i)
     call DNEQNF(DConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
     OwnSConsN(i) = exp(finalB(1))
     OwnSLabN(i) = finalB(2)
  enddo
  
  !Save decision rules to be output to file
  OwnOut(:,1) = OwnCons1
  OwnOut(:,2) = OthConsN
  OwnOut(:,3) = OwnLab1
  OwnOut(:,4) = OthLabN
  OwnOut(:,5) = OwnCap1
  OwnOut(:,6) = OthCapN
  OwnOut(:,7) = OwnInv1
  OwnOut(:,8) = OthInvN
  OwnOut(:,9) = knorm2(:,1)
  OwnOut(:,10) = kss
  OwnOut(:,11) = exp(0.0)
  OwnOut(:,12) = exp(0.0)

  OthOut(:,1) = OthCons1
  OthOut(:,2) = OwnConsN
  OthOut(:,3) = OthLab1
  OthOut(:,4) = OwnLabN
  OthOut(:,5) = OthCap1
  OthOut(:,6) = OwnCapN
  OthOut(:,7) = OthInv1
  OthOut(:,8) = OwnInvN
  OthOut(:,9) = kss
  OthOut(:,10) = knorm2(:,1)
  OthOut(:,11) = exp(0.0)
  OthOut(:,12) = exp(0.0)

  OwnSOut(:,1) = OwnSCons1
  OwnSOut(:,2) = OthSConsN
  OwnSOut(:,3) = OwnSLab1
  OwnSOut(:,4) = OthSLabN
  OwnSOut(:,5) = OwnSCap1
  OwnSOut(:,6) = OthSCapN
  OwnSOut(:,7) = OwnSInv1
  OwnSOut(:,8) = OthSInvN
  OwnSOut(:,9) = kss
  OwnSOut(:,10) = kss
  OwnSOut(:,11) = exp(znorm2(:,1))
  OwnSOut(:,12) = exp(0.0)

  OthSOut(:,1) = OthSCons1
  OthSOut(:,2) = OwnSConsN
  OthSOut(:,3) = OthSLab1
  OthSOut(:,4) = OwnSLabN
  OthSOut(:,5) = OthSCap1
  OthSOut(:,6) = OwnSCapN
  OthSOut(:,7) = OthSInv1
  OthSOut(:,8) = OwnSInvN
  OthSOut(:,9) = kss
  OthSOut(:,10) = kss
  OthSOut(:,11) = exp(0.0)
  OthSOut(:,12) = exp(znorm2(:,1))

end subroutine DecRule2
