!**********************************************************************************
!
!	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), or Gauss-Hermite quadrature (check1 = 3).
  !   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 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, Out
  real(prec) :: AggCons
 
  real(prec),dimension(N) :: Err1,MaxErr1 !intertemporal euler equation errors
  real(prec),dimension(N) :: Err2,MaxErr2 !consumption smoothing errors 
  real(prec) :: Err3,MaxErr3	          !error in world budget constraint
  real(prec),dimension(N) :: MPS
  real(prec) :: IntResult
  real(prec) :: income,expend,incomeA,expendA
  real(prec) :: MaxError1, MaxError2
  
  !Calculate optimal policy functions (either consumption or capital) given state 
  real(prec),external :: ChebValue
  
  !Variables used in DNEQNF function - in solving Euler Equations to determine policy function
  integer,parameter :: itmax = 200
  integer,parameter :: count = 1
  real(prec) :: Relerr,fnorm
  real(prec),dimension(count) :: guess, final

  ! 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
  
  !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
  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', rad
!        else if (statevec(j+N) > zmax) then
!           statevec(j+N) = zmax
!           print*, 'Shock above zmax in AccCheck1sub', rad
!        end if
     end do

     call ChebPoly(d, statevec, Ttemp)
!     AggCons = 0.0
     do j = 1,N
        Cap(j) = ChebValue(Ttemp,CoeffsCap1(:,:,j), & 
             & CoeffsCap2(:,:,:,:,j),CoeffsCap3(:,:,j),YCap0(j))
        Cons(j) = ChebValue(Ttemp,CoeffsCons1(:,:,j), & 
             & CoeffsCons2(:,:,:,:,j),CoeffsCons3(:,:,j),YCons0(j))
        Out(j) = F(statevec(j),statevec(j+N)) 
!        AggCons = Out(j) + statevec(j) - Cap(j) - AdjCost(statevec(j),Cap(j)) + AggCons
     enddo
!     do j = 1,N
!        Cons(j) = AggCons/scale*(Pareto(j)/Pareto(1))**gamma
!     enddo
     
     !Construct vectors for capital (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
        elseif (check1 == 2) then
           call MCIntegrate (j,IntResult)	!Monte Carlo Integration
        else 
           call GHIntegrate (j,IntResult,3)	!Gauss-Hermite Quadrature
        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))))
     end do
     
     !Calculate Error in consumption-smoothing first-order conditions
     do j = 1,N
        Err2(j) = abs(Pareto(1)*Uc(Cons(1))/(Pareto(j)*Uc(Cons(j))) - 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
     
     MaxErr1 = max(Err1,MaxErr1)
     MaxErr2 = max(Err2,MaxErr2)
     MaxErr3 = max(Err3,MaxErr3)
  end do

  MaxError1 = maxval(MaxErr1)
  MaxError2 = maxval(MaxErr2)
  MaxError = max(MaxError1, MaxError2, MaxErr3)
print*, MaxError1, MaxError2, MaxErr3
     
end subroutine AccCheck1sub
