!**********************************************************************************
!
!	SUBROUTINE: AccCheck1sub
!
!	PURPOSE: Computes conditional error functions at a point in the state
!                space.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    7/24/2007
!       REVISED: 4/2009
!**********************************************************************************

subroutine AccCheck1sub(check1,check2,MaxError)
  
  ! inputs:
  !   check1 -- denotes whether integration will be done with monomial formula (check1 = 1)
  !            or monte-carlo routine (check1 = 2).
  !   check2 -- denotes whether world budget constraint is done with baseline (1) or
  !             alternative (2) specification
  !
  ! output:
  !   MaxError -- Maximum error across all equilibrium conditions at all points.
  !
  ! remarks:
  !
  
  use Params
  use StaticFOC
  use AltVector

  implicit none
  
  !Variable Declaration
  integer,intent(in) :: check1,check2
  real(prec), intent(out) :: MaxError

  integer :: h,i,j	
  real(prec) :: dotprod
  real(prec),dimension(d) :: statevec
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  real(prec),dimension(N) :: Cap, Cons, Lab, Out
 
  real(prec),dimension(N) :: Err1,MaxErr1	!intertemporal euler equation errors
  real(prec),dimension(N) :: Err2,MaxErr2	!consumption smoothing (across countries) errors 
  real(prec) :: Err3,MaxErr3	          !error in world budget constraint
  real(prec),dimension(N) :: Err4,MaxErr4 !labor supply errors 
  real(prec),dimension(N) :: MPS
  real(prec),dimension(N) :: MRS, MPL
  real(prec) :: IntResult
  real(prec) :: income,expend,incomeA,expendA
  real(prec) :: MaxError1, MaxError2, MaxError4

  integer :: temp
  
  !Function used to calculate optimal policy functions (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

  !Variables used in nle solver in case tech shock is outside of grid
  integer,parameter :: itmax = 200
  integer,parameter :: count = 1
  real(prec) :: Relerr,fnorm
  real(prec),dimension(count) :: guess, final

  !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

  ! Initialize inputs to nonlinear equation solver, DNEQNF
  Relerr = 0.0001
  fnorm = 0
  guess = -6
  final = 0.1
  
  !Initialize Errors
  Err1 = 0
  MaxErr1 = 0
  Err2 = 0 
  MaxErr2 = 0
  Err3 = 0 
  MaxErr3 = 0
  Err4 = 0 
  MaxErr4 = 0
  
  !Scale random numbers so that points is rad units away from deterministic steady-state
  do i = 1,100
     dotprod = 0.0
     do j = 1,d
        dotprod = accpt(i,j)**2 + dotprod
     end do
     dotprod = sqrt(dotprod)
     accptr(i,:) = rad / dotprod * accpt(i,:) 
  end do

  !Make sure exogenous state variable is not outside of interval for tech shock
  !If it is, choose a different point on the ellipse.
  do i = 1,100
     do j = 1,N
        statevec(j+N) = log(1.0 + accptr(i,j+N))
!        if (statevec(j+N) < zmin) then
!           accptr(i,j+N) = exp(zmin) - 1.0 + 0.001
!           ctryj = j
!           pti = i
!           call DNEQNF(AltVec, Relerr, count, itmax, guess, final, fnorm)
!           do h = 1,N
!              if (accptr(i,h) < 0) then
!                 accptr(i,h) = accptr(i,h) - exp(final(1))
!              else
!                 accptr(i,h) = accptr(i,h) + exp(final(1))
!             endif
!              if (h .ne. j) then
!                 if (accptr(i,h+N) < 0) then
!                    accptr(i,h+N) = accptr(i,h+N) - exp(final(1))
!                 else
!                    accptr(i,h+N) = accptr(i,h+N) + exp(final(1))
!                 endif
!              endif
!          enddo
!        else if (statevec(j+N) > zmax) then
!           accptr(i,j+N) = exp(zmax) - 1.0 - 0.0001
!           ctryj = j
!           pti = i
!           call DNEQNF(AltVec, Relerr, count, itmax, guess, final, fnorm)
!           do h = 1,N
!              if (accptr(i,h) < 0) then
!                 accptr(i,h) = accptr(i,h) - exp(final(1))
!              else
!                 accptr(i,h) = accptr(i,h) + exp(final(1))
!             endif
!              if (h .ne. j) then
!                 if (accptr(i,h+N) < 0) then
!                    accptr(i,h+N) = accptr(i,h+N) - exp(final(1))
!                 else
!                    accptr(i,h+N) = accptr(i,h+N) + exp(final(1))
!                 endif
!              endif
!           enddo
!       end if
     end do
  end do

  !Compute Conditional Error functions at each point
  sn = 1  !Used for passing global variables to nonlinear equation solver
  tc = 2  !''
  do i = 1,100
     
     do j = 1,N
        statevec(j) = 1.0 + accptr(i,j)
        statevec(j+N) = log(1.0 + accptr(i,j+N))
!        if (statevec(j+N) < zmin) then
!           statevec(j+N) = zmin
!           print*, 'Shock below zmin in AccCheck1sub.'
!        else if (statevec(j+N) > zmax) then
!           statevec(j+N) = zmax
!           print*, 'Shock above zmax in AccCheck1sub.'
!        endif
        SimCapital(sn,tc-1,j) = statevec(j)
        shock(sn,tc-1,j) = statevec(j+N)
     end do

     call ChebPoly(d, statevec, Ttemp)
     !Compute labor supply and consumption of country 1
     Cons(1) = ChebValue(Ttemp,CoeffsCons1(:,:,1), &
          & CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
     guessA = lss
     SimConsumption(sn,tc-1,1) = Cons(1)
     call DNEQNF(AConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
     Lab(1) = finalA(1)
     SimLabor(sn,tc-1,1) = Lab(1)
     Out(1) = F(statevec(1),Lab(1),statevec(1+N)) 
     do j = 1,N
        Cap(j) = ChebValue(Ttemp,CoeffsCap1(:,:,j), & 
             & CoeffsCap2(:,:,:,:,j),CoeffsCap3(:,:,j),YCap0(j))
     enddo
     !Compute consumption, labor supply and output of countries 2 - N.
     do FOCi = 2,N
        guessB(1) = log(Cons(1))
        guessB(2) = Lab(1)
        call DNEQNF(AConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
        Cons(FOCi) = exp(finalB(1))
        Lab(FOCi) = finalB(2)
        Out(FOCi) = F(statevec(FOCi),Lab(FOCi),statevec(FOCi+N))
     end do
     
     !Construct vectors for capital stocks (tom) and tech shocks (today) of countries 1,...,N
     do j = 1,N
        ktemp(j) = Cap(j)
        Simtech(j) = statevec(j+N)
     enddo
     
     call ChebPoly(N, ktemp, unitTtemp(:,:))
     
     !Compute MPS -- i.e., RHS of Intertemporal Euler Equation
     do j = 1,N
        if (check1 == 1) then
           call MonIntegrate (j,IntResult,3)	!Integration by Monomial formula
        else
           call MCIntegrate (j,IntResult)	!Monte Carlo Integration
        endif
        
        MPS(j) = betah*IntResult/(1.0 + AdjCost_ktom(statevec(j),Cap(j)))
     end do
     
     !Calculate Error in Intertemporal Euler Equation
     do j = 1,N
        Err1(j) = abs(1.0 - MPS(j)/(Pareto(1)*Uc(Cons(1),Lab(1),1)))
     end do
     
     !Calculate Error in consumption-smoothing first-order conditions
     do j = 1,N
        Err2(j) = abs(Pareto(j)*Uc(Cons(j),Lab(j),j)/(Pareto(1)*Uc(Cons(1),Lab(1),1)) - 1.0)
     end do
     
     !Calculate Error in the world budget constraint
     income = 0
     expend = 0
     incomeA = 0
     expendA = 0
     
     do j = 1,N
        income = Out(j) + statevec(j) + income 
        expend = Cons(j) + Cap(j) + AdjCost(statevec(j),Cap(j)) + expend 
        incomeA = Out(j) - AdjCost(statevec(j),Cap(j)) + incomeA 
        expendA = Cons(j) + Cap(j) - statevec(j) + expendA 
     enddo
     
     if (check2 == 1) then
        Err3 = abs(incomeA/expendA - 1.0)
     else
        Err3 = abs(income/expend - 1.0)
     end if
     
     !Calculate Error in labor supply condition
     do j = 1,N
        MRS(j) = UlUc(Cons(j), Lab(j),j)
        MPL(j) = Fl(statevec(j), Lab(j), statevec(j+N))
        Err4(j) = abs(MPL(j)/MRS(j) + 1.0)
     end do

     MaxErr1 = max(Err1,MaxErr1)
     MaxErr2 = max(Err2,MaxErr2)
     MaxErr3 = max(Err3,MaxErr3)
     MaxErr4 = max(Err4,MaxErr4)
  end do

  MaxError1 = maxval(MaxErr1)
  MaxError2 = maxval(MaxErr2)
  MaxError4 = maxval(MaxErr4)
  MaxError = max(MaxError1, MaxError2, MaxErr3, MaxError4)

!  print*, MaxErr1, MaxErr2, MaxErr3, MaxErr4    

end subroutine AccCheck1sub
