!**********************************************************************************
!
!	SUBROUTINE: Simul
!
!	PURPOSE: Simulates model once policy functions have been found.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 7/23/2007
!
!**********************************************************************************

subroutine Simul(check)
  
  ! inputs:
  !   check -- denotes whether integration will be done with monomial formula (check = 1)
  !            or monte-carlo routine (check = 2).
  !
  ! output:
  ! 
  ! remarks:
  !   Simulates economy and checks Euler Equation errors along simulated path.
  !
  
  use Params
  
  implicit none
  
  !Variable Declaration
  integer,intent(in) :: check
  integer :: i,j	
  integer,parameter :: num = N * maxsamp * samplen		
  real(prec),dimension(num) :: nr
  real(prec),dimension(d) :: statevec
  real(prec),dimension(2**(q-d)+1,d) :: Ttemp
  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) :: Err4,MaxErr4,MeanErr4	!labor supply condition error 
  real(prec),dimension(N) :: mincons, mincap, maxcap 	!min and max values of sim cons and capital
  real(prec),dimension(N) :: minlab, maxlab 	        !min and max values of sim lab
  real(prec),dimension(N) :: minz, maxz                 !min and max values of tech shock
  real(prec),dimension(N) :: MPS
  real(prec) :: IntResult
  real(prec) :: income,expend
  real(prec),dimension(N) :: MRS, MPL

  !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
  Err4 = 0 
  MaxErr4 = 0
  MeanErr4 = 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
  minlab = 2
  maxlab = 0
  minz = 2
  maxz = -2

  !Initialize matrices which hold simulation output
  SimCapital = 0.0
  SimOutput = 0.0
  SimConsumption = 0.0
  SimInvestment = 0.0
  SimLabor = 0.0
  Simtech = 0.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,maxsamp
     do i = 1, samplen
        randnr(sn,i,1:N) = nr((sn-1)*samplen*N + (i-1)*N + 1 : (sn-1)*samplen*N + i*N)
     end do
  enddo
  
  !Simulate the economy
  do sn=1,maxsamp
     
     !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
     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
        end do
     end do
     
     !Compute sequences of capital, consumption, investment, output, labor
     ! Set initial capital stock
     SimCapital(sn,1,1:N) = kss 
     
     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)
        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 simul 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 simul 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))
           SimLabor(sn,tc-1,i) = ChebValue(Ttemp,CoeffsLab1(:,:,i), & 
                & CoeffsLab2(:,:,:,:,i),CoeffsLab3(:,:,i),YLab0(i))
           SimOutput(sn,tc-1,i) = F(SimCapital(sn,tc-1,i),SimLabor(sn,tc-1,i),shock(sn,tc-1,i)) 
        enddo
        
        !Update minimal and maximal capital and labor 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)
        minlab = min(SimLabor(sn,tc-1,:), minlab)
        maxlab = max(SimLabor(sn,tc-1,:), maxlab)
        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)		
        MaxErr4 = max(Err4,MaxErr4)		
        
        !Construct vectors to hold capital stocks (tomorrow) and technological 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 (check == 1) then
              call MonIntegrate (i,IntResult,3)	!Integration by Monomial formula
           else
              call MCIntegrate (i,IntResult)	!Monte Carlo 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(i)/Uc(SimConsumption(sn,tc-1,i),SimLabor(sn,tc-1,i)))
        end do
        
        !Calculate Error in consumption-smoothing first-order conditions
        do i = 1,N
           Err2(i) = abs(Pareto(i)*Uc(SimConsumption(sn,tc-1,i),SimLabor(sn,tc-1,i))/ &
                & (Pareto(1)*Uc(SimConsumption(sn,tc-1,1),SimLabor(sn,tc-1,1))) - 1.0)
        end do
        
        !Calculate Error in the world budget constraint
        income = 0
        expend = 0
        
        do i = 1,N
           income = SimOutput(sn,tc-1,i) - AdjCost(SimCapital(sn,tc-1,i),SimCapital(sn,tc,i)) + income 
           expend = SimConsumption(sn,tc-1,i) + SimCapital(sn,tc,i) & 
                &  - SimCapital(sn,tc-1,i) + expend 
        enddo
        
        Err3 = abs (income/expend - 1.0)
        
        !Calculate Error in labor supply condition
        do i = 1,N
           MRS(i) = UlUc(SimConsumption(sn,tc-1,i), SimLabor(sn,tc-1,i))
           MPL(i) = Fl(SimCapital(sn,tc-1,i), SimLabor(sn,tc-1,i), shock(sn,tc-1,i))
           Err4(i) = abs(MRS(i)/MPL(i) + 1.0)
        end do

        MeanErr1 = Err1 + MeanErr1
        MeanErr2 = Err2 + MeanErr2
        MeanErr3 = Err3 + MeanErr3
        MeanErr4 = Err4 + MeanErr4
        
        if (check == 2) then	!Print to provide estimate of amount of time required by monte-carlo integrator.
           print*, tc
        endif
     end do
     
     MaxErr1 = max(Err1,MaxErr1)
     MaxErr2 = max(Err2,MaxErr2)
     MaxErr3 = max(Err3,MaxErr3)
     MaxErr4 = max(Err4,MaxErr4)
     MeanErr1 = MeanErr1/real(samplen-1, prec)
     MeanErr2 = MeanErr2/real(samplen-1, prec)
     MeanErr3 = MeanErr3/real(samplen-1, prec)
     MeanErr4 = MeanErr4/real(samplen-1, prec)
     
     !Prepare output matrix
     do i = 1,N 
        SimulOut(i,1) = MaxErr1(i)
        SimulOut(i,2) = MaxErr2(i)
        SimulOut(i,3) = MaxErr3
        SimulOut(i,4) = MaxErr4(i)
        SimulOut(i,5) = MeanErr1(i)
        SimulOut(i,6) = MeanErr2(i)
        SimulOut(i,7) = MeanErr3
        SimulOut(i,8) = MeanErr4(i)
        SimulOut(i,9) = mincap(i)
        SimulOut(i,10) = maxcap(i)
        SimulOut(i,11) = mincons(i)
        SimulOut(i,12) = minlab(i)
        SimulOut(i,13) = maxlab(i)
        SimulOut(i,14) = minz(i)
        SimulOut(i,15) = maxz(i)
     end do

!     do tc = 1,samplen-1
        ! Compute sequence of investment for each of N countries
!        do i = 1,N
!           SimInvestment(sn,tc,i) = SimCapital(sn,tc+1,i) - (1.0 - delta) * SimCapital(sn,tc,i) 
!        end do
!     end do
     
  end do
  
  !Calculate Running Time
  call CPU_TIME (time_end)
  simul_runtime = time_end - time_begin
  
end subroutine Simul
