!**********************************************************************************
!
!	SUBROUTINE: AccCheck3
!
!	PURPOSE: Computes DenHaan-Marcet statistic for model.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    8/15/2007
!       REVISED: 9/4/2007
!**********************************************************************************

subroutine AccCheck3
  
  ! inputs:
  !
  ! output:
  ! 
  ! remarks:
  !
  
  use Params     !Contains global variables
  use StaticFOC  !Contains equations that determine labor supply, consumption

  implicit none
  
  !Variable Declaration
  integer :: i,j
  integer,parameter :: num = (N+1) * maxsamp * 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(N) :: Err	!approximation errors
  real(prec),dimension(N) :: MPS 
  real(prec),dimension(2*N*N+3*N+1) :: Inst !Instruments to check for orthogonality
  integer,dimension(N) :: sumi  

  integer,parameter :: matsize = N*(2*N*N+3*N+1) 
  real(prec),dimension(matsize) :: colB
  real(prec),dimension(matsize) :: B_T
  real(prec),dimension(matsize,matsize) :: A_T
  real(prec),dimension(matsize,matsize) :: invA_T
  real(prec),dimension(matsize) :: AB_T
  real(prec) :: dm
  integer :: P_05, P_50, P_95
  
  real(prec),dimension(matsize) :: WK

  !Calculates optimal policy funcs (either cons or cap) given state 
  real(prec),external :: ChebValue
  
  !Variables used in nle 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 nle 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

  !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

  !Start timer
  call CPU_TIME (time_begin)
  
  !Initialize output statistics
  P_05 = 0
  P_50 = 0
  P_95 = 0

  !Construct an index for use in constructing instruments
  sumi(1) = 0
  do i = 2,N
     sumi(i) = sumi(i-1) + i-1
  end do
  
  !Initialize matrices which hold simulation output (global variables)
  SimCapital = 0.0
  SimOutput = 0.0
  SimConsumption = 0.0
  SimTech = 0.0
  SimLabor = 0.0

  !Initialize 'ones' vector
  ones = 1.0

  !Initialize constant instrument
  Inst(1) = 1.0

  !Draw a bunch of random numbers
  call rnset(1)   !If 0, different random numbers drawn each time program is run.
  call DRNNOA(num,nr)	!IMSL random number generator (normal distribution)
  
  !Subdivide random numbers into sample paths
  do sn = 1, maxsamp
     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, maxsamp

     !Initialize inputs to DenHaan-Marcet Statistic
     A_T = 0.0
     B_T = 0.0
     
     !Initial Shock (t=1)
     do i = 1,N
!        if (randnr(sn,1,i) * sigma < zmin) then
!           shock(sn,1,i) = zmin
!        elseif (randnr(sn,1,i) * sigma > zmax) then
!           shock(sn,1,i) = zmax
!        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
!           elseif (rho * shock(sn, tc-1, i) + randnr(sn,tc,i) * sigma > zmax) then 
!              shock(sn,tc,i) = zmax
!           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

     ! Set initial capital stock
     SimCapital(sn,1,1:N) = kss 

     ! Compute sequence of capital stocks, consumption, and net output
     do tc = 2,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)
        !Compute labor supply and consumption of country 1
        SimConsumption(sn,tc-1,1) = ChebValue(Ttemp,CoeffsCons1(:,:,1), &
             & CoeffsCons2(:,:,:,:,1),CoeffsCons3(:,:,1),YCons0(1))
        guessA = lss
        call DNEQNF(AConsLabOne, RelerrA, countA, itmaxA, guessA, finalA, fnormA)
        SimLabor(sn,tc-1,1) = finalA(1)
        SimOutput(sn,tc-1,1) = F(SimCapital(sn,tc-1,1),SimLabor(sn,tc-1,1),shock(sn,tc-1,1)) 
        
        do i=1,N
           SimCapital(sn,tc,i) = ChebValue(Ttemp,CoeffsCap1(:,:,i), & 
                & CoeffsCap2(:,:,:,:,i),CoeffsCap3(:,:,i),YCap0(i))
           if (SimCapital(sn,tc,i) < kmin(i)) then
              print*, 'Capital below kmin in AccCheck3 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(i)) then
              print*, 'Capital above kmax in AccCheck3 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
        end do

        !Compute consumption, labor supply and output of countries 2 - N.
        do FOCi = 2,N
           guessB(1) = log(SimConsumption(sn,tc-1,1))
           guessB(2) = SimLabor(sn,tc-1,1)
           call DNEQNF(AConsLabOther, RelerrB, countB, itmaxB, guessB, finalB, fnormB)
           SimConsumption(sn,tc-1,FOCi) = exp(finalB(1))
           SimLabor(sn,tc-1,FOCi) = finalB(2)
           SimOutput(sn,tc-1,FOCi) = F(SimCapital(sn,tc-1,FOCi),SimLabor(sn,tc-1,FOCi),shock(sn,tc-1,FOCi))
        end do
     enddo
     
     do tc = 201,samplen-1
        !Calculate Error in Realization of Intertemporal Euler Equation
        do i = 1,N
           MPS(i) = betah * Pareto(i) * Uc(SimConsumption(sn,tc,i),SimLabor(sn,tc,i),i) * & 
                & (1.0 + Fk(SimCapital(sn,tc,i),SimLabor(sn,tc,i),shock(sn,tc,i)) - & 
                & AdjCost_k(SimCapital(sn,tc,i),SimCapital(sn,tc+1,i))) / &
                & (1.0 + AdjCost_ktom(SimCapital(sn,tc-1,i),SimCapital(sn,tc,i)))
           Err(i) = Pareto(i)*Uc(SimConsumption(sn,tc-1,i),SimLabor(sn,tc-1,i),i) - MPS(i)
        end do
        
        !Construct Instruments to be used in DenHaan-Marcet Statistic
        do i = 1,N
           Inst(1+i) = log(SimCapital(sn,tc-1,i))
           Inst(N+1+i) = shock(sn,tc-1,i)
           Inst(2*N+1+i) = Inst(1+i) * Inst(1+i)
           Inst(3*N+1+i) = Inst(N+1+i) * Inst(N+1+i)
        end do
        
        do i = 1,N
           do j = 1,N
              Inst(4*N+1+(i-1)*N+j) = Inst(1+i) * Inst(N+1+j)
           end do
           do j = i+1,N
              Inst(4*N+1+N*N+2*((i-1)*N-sumi(i))+2*(j-i-1)+1) = Inst(1+i) * Inst(1+j)
              Inst(4*N+1+N*N+2*((i-1)*N-sumi(i))+2*(j-i-1)+2) = Inst(N+1+i) * Inst(N+1+j)
           end do
        end do

        !Construct Inputs to DM-Statistic
        do i = 1,N
           do j = 1,2*N*N+3*N+1
              colB((i-1)*(2*N*N+3*N+1) + j) = Err(i)*Inst(j)
              B_T((i-1)*(2*N*N+3*N+1) + j) = colB((i-1)*(2*N*N+3*N+1)+j) + B_T((i-1)*(2*N*N+3*N+1)+j)
           end do
        end do

        do i = 1,matsize
           do j = 1,matsize
              A_T(i,j) = colB(i)*colB(j) + A_T(i,j)
           end do
        end do
     end do
     
     B_T = B_T / real(samplen-201,prec) 
     A_T = A_T / real(samplen-201,prec) 
     
     !Compute inverse of A_T
     CALL DL2NDS (matsize, A_T, matsize, invA_T, matsize, WK)

     !Construct Test Statistic
     AB_T = 0.0
     dm = 0.0
     do i = 1,matsize
        do j = 1,matsize
           AB_T(i) = invA_T(i,j) * B_T(j) + AB_T(i)
        enddo
     enddo
     do i = 1,matsize
        dm = B_T(i)*AB_T(i) + dm
     enddo
     dm = real(samplen-201,prec)*dm
     
     !Compare Test Statistic to chi-sqaured distribution, dof N*(2*N*N+3*N+1)
     if (N == 2) then
        if (dm < 18.49 ) then
           P_05 = P_05 + 1
           P_50 = P_50 + 1
        elseif (dm > 18.49  .and. dm < 29.34 ) then
           P_50 = P_50 + 1
        elseif (dm > 43.77) then
           P_95 = P_95 + 1
        endif
     elseif (N == 4) then
       if (dm < 149.97 ) then
           P_05 = P_05 + 1
           P_50 = P_50 + 1
        elseif (dm > 149.97  .and. dm < 179.33 ) then
           P_50 = P_50 + 1
        elseif (dm > 212.30) then
           P_95 = P_95 + 1
        endif
     elseif (N == 6) then
        if (dm < 492.81 ) then
           P_05 = P_05 + 1
           P_50 = P_50 + 1
        elseif (dm > 492.81  .and. dm < 545.33 ) then
           P_50 = P_50 + 1
        elseif (dm > 601.47) then
           P_95 = P_95 + 1
        endif
     endif
  end do
  
  !Prepare Output Matrix
  Acc3Out(1,1) = real(P_05,prec)/real(maxsamp,prec)
  Acc3Out(1,2) = real(P_50,prec)/real(maxsamp,prec)
  Acc3Out(1,3) = real(P_95,prec)/real(maxsamp,prec) 

  !Calculate Running Time
  call CPU_TIME (time_end)
  Acc3_runtime = time_end - time_begin
  
end subroutine AccCheck3
