!***********************************************************************************
!
!    FUNCTION: ExpectFOC
!
!    PURPOSE: Specifies integrand on the rhs of the euler-equation.
!              Will be integrated by MonIntegrate subroutine, when called in 
!              EulerEq subroutine. 
!
!    AUTHOR:  Ben Malin 
!
!    DATE:  6/2004
!    REVISED: 7/23/2007
!
!***********************************************************************************

function ExpectFOC(dim1, techshock, nc)

  ! inputs:
  !    dim1 -- number of countries: N
  !    techshock -- holds the technological innovation for each of N countries and aggregate shock
  !    nc -- index denotes which country's euler equation is being evaluated
  !
  ! output:
  ! 
  ! remarks:
  !
  
  use Params
  implicit none
  
  !Variable Declaration
  integer, intent(in) :: dim1
  real(prec), dimension(dim1+1), intent(in) :: techshock
  integer, intent(in) :: nc
  real(prec) :: ExpectFOC
  real(prec) :: constom, captomtom, labtom
  real(prec), dimension(N) :: newstate, techstate
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  integer :: i,j
  real(prec) :: MPKtom, MPCtom, density
  
  !Returns policy choice (either consumption or capital) using Chebyshev polynomial for given state
  real(prec), external :: ChebValue  
  
  !Construct z shock
  do i = 1,dim1
     newstate(i) = rho * z(stind(N+i),i) + sigma * sqrt(2.0) * (techshock(i) + techshock(dim1+1))  !Constucts tech shock on [zmin, zmax] grid
!     if (newstate(i) < zmin) then
!        newstate(i) = zmin
!        print*, 'Tech shock below zmin in ExpectFOC', stind
!        pause
!     else if (newstate(i) > zmax) then
!        newstate(i) = zmax
!        print*, 'Tech shock above zmax in ExpectFOC', stind
!        pause
!     endif
  enddo
  
  do i = 1, dim1
     techstate(i) = (newstate(i) - zmin)*2.0/(zmax-zmin) - 1.0	!Converts tech shock to [-1,1] grid
  enddo
  
  !Construct Chebyshev polynomials from capital stocks
  Ttemp(:,1:N) = unitTtemp(:,:)
  
  !Construct Chebyshev polynomials from technological shocks
  do i = 1,dim1
     Ttemp(1,i+N) = 1.0
     Ttemp(2,i+N) = techstate(i)
     do j = 3,5			!Since q = d+2, 2**(q-d)+1 = 5
        Ttemp(j,i+N) = 2.0 * techstate(i) * Ttemp(j-1,i+N) - Ttemp(j-2,i+N)
     end do
  enddo
  
  !Compute tomorrow's consumption, labor and capital stock in t+2 using Chebyshev polynomial
  constom = ChebValue(Ttemp,CoeffsCons1(:,:,nc),CoeffsCons2(:,:,:,:,nc),CoeffsCons3(:,:,nc),OldyCons0(nc))
  captomtom = ChebValue(Ttemp,CoeffsCap1(:,:,nc),CoeffsCap2(:,:,:,:,nc),CoeffsCap3(:,:,nc),OldyCap0(nc))
  labtom = (Lfoc * exp(newstate(nc)) * ktemp(nc)**alpha * constom**(-xgamma))**(Lelast)

  !Compute MPKtom
  MPKtom = 1.0 + Fk(ktemp(nc),labtom,newstate(nc)) - AdjCost_k(ktemp(nc),captomtom)
  
  !Compute MPCtom
  MPCtom = Pareto(nc) * Uc(constom)
  
  !Compute Density -- Note (7.5.8) in Judd
  density = pie ** (-(dim1+1) * 0.5)
  
  !Specify Integrand
  ExpectFOC = MPKtom * MPCtom * density
  
end function ExpectFOC
  
