!**********************************************************************************
!
!	SUBROUTINE:  DecRule2
!
!	PURPOSE: Calculates decision rules for comparison project.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 12/10/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
  
  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
  
  !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)
     OwnLab1(i) = 1
     OthCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OthConsN(i) = ChebValue(Ttemp,CoeffsCons1(:,:,N),CoeffsCons2(:,:,:,:,N),CoeffsCons3(:,:,N),YCons0(N))
     OthInvN(i) = OthCapN(i) - (1.0 - delta)*DecPoint(N)
     OthLabN(i) = 1
  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)
     OthLab1(i) = 1
     OwnCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OwnConsN(i) = ChebValue(Ttemp,CoeffsCons1(:,:,N),CoeffsCons2(:,:,:,:,N),CoeffsCons3(:,:,N),YCons0(N))
     OwnInvN(i) = OwnCapN(i) - (1.0 - delta)*DecPoint(N)
     OwnLabN(i) = 1
  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)
     OwnSLab1(i) = 1
     OthSCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OthSConsN(i) = ChebValue(Ttemp,CoeffsCons1(:,:,N),CoeffsCons2(:,:,:,:,N),CoeffsCons3(:,:,N),YCons0(N))
     OthSInvN(i) = OthSCapN(i) - (1.0 - delta)*DecPoint(N)
     OthSLabN(i) = 1
  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)
     OthSLab1(i) = 1
     OwnSCapN(i) = ChebValue(Ttemp,CoeffsCap1(:,:,N),CoeffsCap2(:,:,:,:,N),CoeffsCap3(:,:,N),YCap0(N))
     OwnSConsN(i) = ChebValue(Ttemp,CoeffsCons1(:,:,N),CoeffsCons2(:,:,:,:,N),CoeffsCons3(:,:,N),YCons0(N))
     OwnSInvN(i) = OwnSCapN(i) - (1.0 - delta)*DecPoint(N)
     OwnSLabN(i) = 1
  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

  OthOut(:,1) = OthCons1
  OthOut(:,2) = OwnConsN
  OthOut(:,3) = OthLab1
  OthOut(:,4) = OwnLabN
  OthOut(:,5) = OthCap1
  OthOut(:,6) = OwnCapN
  OthOut(:,7) = OthInv1
  OthOut(:,8) = OwnInvN

  OwnSOut(:,1) = OwnSCons1
  OwnSOut(:,2) = OthSConsN
  OwnSOut(:,3) = OwnSLab1
  OwnSOut(:,4) = OthSLabN
  OwnSOut(:,5) = OwnSCap1
  OwnSOut(:,6) = OthSCapN
  OwnSOut(:,7) = OwnSInv1
  OwnSOut(:,8) = OthSInvN

  OthSOut(:,1) = OthSCons1
  OthSOut(:,2) = OwnSConsN
  OthSOut(:,3) = OthSLab1
  OthSOut(:,4) = OwnSLabN
  OthSOut(:,5) = OthSCap1
  OthSOut(:,6) = OwnSCapN
  OthSOut(:,7) = OthSInv1
  OthSOut(:,8) = OwnSInvN
end subroutine DecRule2
