!**********************************************************************************
!
!	SUBROUTINE: AccCheck2
!
!	PURPOSE: Simulates model and reports max and mean error along simulated
!                path.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 4/2009
!
!**********************************************************************************

subroutine AccCheck2(check1,check2)
  
  ! 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:
  ! 
  ! remarks:
  !   Simulates economy and checks Euler Equation errors along simulated path.
  !
  
  use Params
  
  implicit none
  
  !Variable Declaration
  integer,intent(in) :: check1,check2
  integer :: i,j	
  integer,parameter :: num = (N+1) * 1 * samplen		
  real(prec),dimension(N) :: ones
  real(prec),dimension(num) :: nr
  real(prec),dimension(d) :: statevec
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  !real(prec),dimension(samplen-1,N) :: sample1con, sample1out, sample1inv, sample1cap
  integer :: flag, capflag1, capflag2
  
  real(prec),dimension(N) :: Err1,MaxErr1,MeanErr1	!intertemporal euler equation errors
  real(prec),dimension(N) :: Err2,MaxErr2,MeanErr2	!consumption smoothing (across countries) errors 
  real(prec) :: Err3,MaxErr3,MeanErr3			!error in world budget constraint
  real(prec),dimension(N) :: mincons, mincap, maxcap 	!min and max values of sim cons and capital
  real(prec),dimension(N) :: minz, maxz 	!min and max values of shocks
  real(prec),dimension(N) :: MPS 
  real(prec) :: IntResult
  real(prec) :: income,expend,incomeA,expendA
  
  !Function used to calculate optimal policy functions (either consumption or capital) given state 
  real(prec),external :: ChebValue
  
  !Start timer
  call CPU_TIME (time_begin)
  
  !Initialize Errors
  Err1 = 0
  MaxErr1 = 0
  MeanErr1 = 0
  Err2 = 0 
  MaxErr2 = 0
  MeanErr2 = 0
  Err3 = 0 
  MaxErr3 = 0
  MeanErr3=0
  
  !Initialize flags that record if state variables go outside of state space
  flag = 0
  capflag1 = 0
  capflag2 = 0
  
  !Initialize minimal and maximal consumption/capital vectors
  mincons = 2
  mincap = 2
  maxcap = 0
  minz = 2
  maxz = -2
  
  !Initialize matrices which hold simulation output
  SimCapital = 0.0
  SimOutput = 0.0
  SimConsumption = 0.0
  SimInvestment = 0.0
  Simtech = 0.0

  !Initialize 'ones' vector
  ones = 1.0

  !Draw a bunch of random numbers
  call rnset(1)
  call DRNNOA(num,nr)	!IMSL random number generator (normal distribution)
  
  !Subdivide random numbers into sample paths
  do sn = 1,1
     do i = 1, samplen
        randnr(sn,i,1:N) = nr((sn-1)*samplen*(N+1) + (i-1)*(N+1) + 2 : (sn-1)*samplen*(N+1) + i*(N+1)) + nr((sn-1)*samplen*(N+1) + (i-1)*(N+1) + 1)*ones
     end do
  enddo
  
  !Simulate the economy
  do sn = 1,1
     
     !Initial Shock (t=1)
     do i = 1,N
!        if (randnr(sn,1,i) * sigma < zmin) then
!           shock(sn,1,i) = zmin
!           flag = flag + 1
!        elseif (randnr(sn,1,i) * sigma > zmax) then
!           shock(sn,1,i) = zmax
!           flag = flag + 1
!        else
!           shock(sn,1,i)= randnr(sn,1,i) * sigma
!        endif
        shock(sn,1,i)= randnr(sn,1,i) * sigma
     end do
     
     !Subsequent shocks (t=2,...,samplen)
     do tc = 2, samplen
        do i = 1,N
!           if (rho * shock(sn, tc-1, i) + randnr(sn,tc,i) * sigma < zmin) then
!              shock(sn,tc,i) = zmin
!              flag = flag + 1
!           elseif (rho * shock(sn, tc-1, i) + randnr(sn,tc,i) * sigma > zmax) then 
!              shock(sn,tc,i) = zmax
!              flag = flag + 1
!           else
!              shock(sn,tc,i) = rho * shock(sn, tc-1, i) + randnr(sn,tc,i) * sigma 
!           endif
           shock(sn,tc,i) = rho * shock(sn, tc-1, i) + randnr(sn,tc,i) * sigma 
        end do
     end do
     
     !Compute sequences of capital, consumption, investment, output
     ! Set initial capital stock
     SimCapital(sn,1,1:N) = kss 
     
     ! Compute sequence of capital stocks, consumption, and net output and check errors in Euler Equations
     do tc = 2,201
        
        do i = 1,N
           statevec(i) = SimCapital(sn,tc-1,i)
           statevec(N+i) = shock(sn,tc-1,i)
        end do
        
        call ChebPoly(d, statevec, Ttemp)
        do i=1,N
           SimCapital(sn,tc,i) = ChebValue(Ttemp,CoeffsCap1(:,:,i), & 
                & CoeffsCap2(:,:,:,:,i),CoeffsCap3(:,:,i),YCap0(i))
!           if (SimCapital(sn,tc,i) < kmin) then
!              capflag1 = capflag1 + 1
!              print*, 'Capital below kmin in AccCheck2 for country', i, ' in period', tc
!              print*, 'State', statevec
!              print*, 'Capital', SimCapital(sn,:,i)
!              print*, 'Shock history', shock(sn,:,1) - shock(sn,:,2)
!              pause
!           elseif (SimCapital(sn,tc,i) > kmax) then
!              capflag2 = capflag2 + 1
!              print*, 'Capital above kmax in AccCheck2 for country', i, ' in period', tc
!              print*, 'State', statevec
!              print*, 'Capital', SimCapital(sn,:,i)
!              print*, 'Shock history', shock(sn,:,i) - shock(sn,:,2)
!              pause
!           endif
           SimConsumption(sn,tc-1,i) = ChebValue(Ttemp,CoeffsCons1(:,:,i), & 
                & CoeffsCons2(:,:,:,:,i),CoeffsCons3(:,:,i),YCons0(i))
           SimOutput(sn,tc-1,i) = F(SimCapital(sn,tc-1,i),shock(sn,tc-1,i)) 
        enddo
     enddo

     do tc = 202,samplen
        
        do i = 1,N
           statevec(i) = SimCapital(sn,tc-1,i)
           statevec(N+i) = shock(sn,tc-1,i)
        end do
        
        call ChebPoly(d, statevec, Ttemp)
        do i=1,N
           SimCapital(sn,tc,i) = ChebValue(Ttemp,CoeffsCap1(:,:,i), & 
                & CoeffsCap2(:,:,:,:,i),CoeffsCap3(:,:,i),YCap0(i))
!           if (SimCapital(sn,tc,i) < kmin) then
!              capflag1 = capflag1 + 1
!              print*, 'Capital below kmin in AccCheck2 for country', i, ' in period', tc
!              print*, 'State', statevec
!              print*, 'Capital', SimCapital(sn,:,i)
!              print*, 'Shock history', shock(sn,:,1) - shock(sn,:,2)
!              pause
!           elseif (SimCapital(sn,tc,i) > kmax) then
!              capflag2 = capflag2 + 1
!              print*, 'Capital above kmax in AccCheck2 for country', i, ' in period', tc
!              print*, 'State', statevec
!              print*, 'Capital', SimCapital(sn,:,i)
!              print*, 'Shock history', shock(sn,:,i) - shock(sn,:,2)
!              pause
!           endif
           SimConsumption(sn,tc-1,i) = ChebValue(Ttemp,CoeffsCons1(:,:,i), & 
                & CoeffsCons2(:,:,:,:,i),CoeffsCons3(:,:,i),YCons0(i))
           SimOutput(sn,tc-1,i) = F(SimCapital(sn,tc-1,i),shock(sn,tc-1,i)) 
        enddo
        
        !Update minimal and maximal capital and minimal consumption vectors
        mincap = min(SimCapital(sn,tc-1,:), mincap)
        maxcap = max(SimCapital(sn,tc-1,:), maxcap)
        mincons = min(SimConsumption(sn,tc-1,:), mincons)
        minz = min(shock(sn,tc-1,:), minz)
        maxz = max(shock(sn,tc-1,:), maxz)
        
        !Check errors in Euler Equation along the simulated path
        MaxErr1 = max(Err1,MaxErr1)
        MaxErr2 = max(Err2,MaxErr2)
        MaxErr3 = max(Err3,MaxErr3)		
        
        !Construct vectors to hold cap stocks (tom) and tech shocks (today) of countries 1-N
        do i = 1,N
           ktemp(i) = SimCapital(sn,tc,i)
           Simtech(i) = shock(sn,tc-1,i)
        enddo
        
        call ChebPoly(N, ktemp, unitTtemp(:,:))
        
        !Compute MPS -- i.e., RHS of Intertemporal Euler Equation
        do i = 1,N
           if (check1 == 1) then
              call MonIntegrate (i,IntResult,3)	!Integration by Monomial formula
           elseif (check1 == 2) then
              call MCIntegrate (i,IntResult)	!Monte Carlo Integration
           else
              call GHIntegrate (i,IntResult,3)	!Gauss-Hermite Integration
           endif
           
           MPS(i) = betah*IntResult/(1 + AdjCost_ktom(SimCapital(sn,tc-1,i),SimCapital(sn,tc,i)))
        end do
        
        !Calculate Error in Intertemporal Euler Equation
        do i = 1,N
           Err1(i) = abs(1.0 - MPS(i)/(Pareto(1)*Uc(SimConsumption(sn,tc-1,1))))
        end do
        
        !Calculate Error in consumption-smoothing first-order conditions
        do i = 1,N
           Err2(i) = abs(Pareto(1)*Uc(SimConsumption(sn,tc-1,1))/ & 
                & (Pareto(i)*Uc(SimConsumption(sn,tc-1,i))) - 1.0)
        end do
     
        !Calculate Error in the world budget constraint
        income = 0
        expend = 0
        incomeA = 0
        expendA = 0

        do i = 1,N
           income = SimOutput(sn,tc-1,i) + SimCapital(sn,tc-1,i) + income 
           expend = SimConsumption(sn,tc-1,i) + SimCapital(sn,tc,i) + & 
                & AdjCost(SimCapital(sn,tc-1,i),SimCapital(sn,tc,i)) + expend 
           incomeA = SimOutput(sn,tc-1,i) - AdjCost(SimCapital(sn,tc-1,i),SimCapital(sn,tc,i)) + incomeA
           expendA = SimConsumption(sn,tc-1,i) + SimCapital(sn,tc,i) - & 
                & SimCapital(sn,tc-1,i) + expendA 
        enddo
         
        if (check2 == 1) then
           Err3 = abs(incomeA/expendA - 1.0)
        else
           Err3 = abs(income/expend - 1.0)
        end if
                
        MeanErr1 = Err1 + MeanErr1
        MeanErr2 = Err2 + MeanErr2
        MeanErr3 = Err3 + MeanErr3
        
        if (check1 == 2) then	!Provide estimate of amount of time required by mc integrator.
           print*, tc
        endif
     end do
     
     MaxErr1 = max(Err1,MaxErr1)
     MaxErr2 = max(Err2,MaxErr2)
     MaxErr3 = max(Err3,MaxErr3)
     MeanErr1 = MeanErr1/real(samplen-201, prec)
     MeanErr2 = MeanErr2/real(samplen-201, prec)
     MeanErr3 = MeanErr3/real(samplen-201, prec)
     
     !Prepare output matrix
     do i = 1,N 
        Acc2Out(i,1) = MaxErr1(i)
        Acc2Out(i,2) = MaxErr2(i)
        Acc2Out(i,3) = MaxErr3
        Acc2Out(i,4) = MeanErr1(i)
        Acc2Out(i,5) = MeanErr2(i)
        Acc2Out(i,6) = MeanErr3
        Acc2Out(i,7) = mincap(i)
        Acc2Out(i,8) = maxcap(i)
        Acc2Out(i,9) = mincons(i)
        Acc2Out(i,10) = minz(i)
        Acc2Out(i,11) = maxz(i)
     end do

  end do
  
!  !Print number of times tech shocks go outside of bounds
!  print*,"Flag = ", flag

  !Calculate Running Time
  call CPU_TIME (time_end)
  Acc2_runtime = time_end - time_begin
  
end subroutine AccCheck2
