!****************************************************************************
!
!	MODULE:  Params
!
!	PURPOSE:  Declares and initializes global variables
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:  5/30/2007 (orginally 6/2003)
!       REVISED: 8/23/2007
!****************************************************************************

module Params
  
  ! inputs:
  !
  ! output:
  ! 
  ! remarks:
  
  implicit none
  
  integer,parameter :: prec = selected_real_kind(15,307)
  real(prec),parameter :: pie = 3.1415926535897932
  
  !------------------------------
  !Number of 'state' variables
  !------------------------------
  integer,parameter :: N = 2 	     !Number of capital state variables in model (ie, # of countries)  
  integer,parameter :: d = 4 !2*N!   !Dimension of state space (capital and tech shock)
  integer,parameter :: q = 6 !d+2!   !Used in smolyak interpolation algorithm
  integer,parameter :: Nplus = 3 !N+1!  !# of equations in system solved by nonlinear equation solver

  !----------------------
  !Parameters of model
  !----------------------
  real(prec) :: gamma_min
  real(prec) :: gamma_max
  real(prec),dimension(N) :: gamma	!Risk aversion parameter (country dependent)
  real(prec),dimension(N) :: xgamma     !Intertemporal Elasticity of Substitution
  real(prec) :: eta_min
  real(prec) :: eta_max
  real(prec),dimension(N) :: eta        !Labor supply elasticity parameter (country dependent)
  real(prec),dimension(N) :: xeta       !Frisch elasticity of labor supply
  real(prec),parameter :: betah = 0.99	!Subjective discount factor
  real(prec),parameter :: alpha = 0.36	!Capital Share
  real(prec),parameter :: delta = 0.025	!Depreciation
  real(prec) :: kappa = 0		!Adjustment Cost parameter  ("phi" in organizers' notation)
    ! Always set kappa = 0 for 'hot-start'.  Kappa will then be updated in MAIN.
  
  real(prec) :: A_tfp = (1-betah)/(alpha*betah)	 !Total Factor Productivity
  real(prec),dimension(N) :: b   !Utility parameter for labor
  real(prec),parameter :: sigma = 0.01	 !Parameterizes possible values of tech. shock
  real(prec),parameter :: rho = 0.95     !Persistence of tech. shock
  
  real(prec),dimension(N) :: Pareto    !Pareto weights for each of N countries (set in Grid)

  real(prec),dimension(N) :: Lfoc !Parameter in intratemporal first order condition
  real(prec),dimension(N) :: Lelast !Used in intratemporal first order condition
  
  !--------------------------------------------
  !Weights used in formulas for interpolation 
  !--------------------------------------------
  integer :: coc1	!Used in weighted sum of cheb polynomials in Smolyak interpolation
  integer :: coc2	!coc1 is (d-1) choose 1 and coc2 is (d-1) choose 2
  
  !----------------------------------------
  !Weights and points used in integration
  !----------------------------------------
  !Weights used in integration -- See Judd pg. 275, formula (7.5.11)
  real(prec) :: V_mono
  real(prec) :: A_mono
  real(prec) :: B_mono
  real(prec) :: C_mono
  
  !Points used for integration -- See Judd pg. 275, formula (7.5.11)
  real(prec), dimension(N+1) :: point0_mono		!The Origin
  real(prec), dimension(N+1, N+1) :: points1_mono		!"2nd term" points
  real(prec), dimension(N, N+1, 2, N+1) :: points2_mono	!"3rd term" points
  
  real(prec), dimension(N+1,N+1) :: e_mono	!Used to construct grid of points used in integration
                                        !e(i,:) is defined on bottom of pg. 271
  real(prec) :: r_mono		!Used to construct grid of points used in integration (7.5.11)
  real(prec) :: s_mono		!''
  
  !-------------------------------
  !Points used for interpolation
  !-------------------------------
  real(prec),dimension(2**(q-d)+1,N) :: k	! Grid for chebyshev points on [kmin,kmax] x N
  real(prec),dimension(2**(q-d)+1,N) :: z	! Grid for chebyshev points on [zmin, zmax] x N
  real(prec),dimension(2**(q-d)+1) :: x		! Grid for chebyshev points on [-1,1]
  
  !Smolyak points for d dimensional grid -- used for interpolation
  real(prec),dimension(d,2**(q-d),d) :: points3	!Smolyak pts (excl origin) for 4th order polynomial in dir d
  real(prec),dimension(d-1,d,2**(q-d-1),2**(q-d-1),d) :: points2 !Smolyak pts for 2nd order tensor product
  real(prec),dimension(d) :: point0				! Smolyak point -- the origin
  
  !-----------------------------------------------------------------
  !Variables that hold values of policy functions at Smolyak points
  !-----------------------------------------------------------------
  real(prec),dimension(d,2**(q-d),N) :: yCons3	!Policy func (cons) (for N countries) evaluated at points3
  real(prec),dimension(d,2**(q-d),N) :: yCap3	!Policy function (capital) evaluated at points3
  real(prec),dimension(d,2**(q-d),N) :: yLab3	!Policy function (labor) evaluated at points3
  real(prec),dimension(d-1,d,2**(q-d-1),2**(q-d-1),N) :: yCons2	!Policy function (cons) evaluated at points2
  real(prec),dimension(d-1,d,2**(q-d-1),2**(q-d-1),N) :: yCap2	!Policy function (cap) evaluated at points2
  real(prec),dimension(d-1,d,2**(q-d-1),2**(q-d-1),N) :: yLab2	!Policy func (labor) evaluated at points2
  real(prec),dimension(N) :: yCons0		!Policy function (consumption) evaluated at point0
  real(prec),dimension(N) :: yCap0		!Policy function (capital) evaluated at point0
  real(prec),dimension(N) :: yLab0		!Policy function (labor) evaluated at point0

  real(prec),dimension(N) :: OldYCons0		!Policy function (cons) at point0 from previous iteration
  real(prec),dimension(N) :: OldYCap0		!Policy function (capital) at point0 from previous iteration
  real(prec),dimension(N) :: OldYLab0		!Policy function (labor) at point0 from previous iteration
  
  !--------------------------------------------------------
  !Variables that hold coefficients of policy functions 
  !--------------------------------------------------------
  real(prec),dimension(d,2**(q-d)+1,N) :: CoeffsCons3	!Coeffs for 4th order polynomials for each of d dims
  real(prec),dimension(d,2**(q-d)+1,N) :: CoeffsCap3    !Coeffs for 4th order polynomials for each of d dims
  real(prec),dimension(d,2**(q-d)+1,N) :: CoeffsLab3    !Coeffs for 4th order polynomials for each of d dims
  real(prec),dimension(d,2**(q-d-1)+1,N) :: CoeffsCons1	!Coeffs for 2nd order polynomials for each of d dims
  real(prec),dimension(d,2**(q-d-1)+1,N) :: CoeffsCap1	!Coeffs for 2nd order polynomials for each of d dims
  real(prec),dimension(d,2**(q-d-1)+1,N) :: CoeffsLab1	!Coeffs for 2nd order polynomials for each of d dims
  real(prec),dimension(d-1,d,2**(q-d-1)+1,2**(q-d-1)+1,N) :: CoeffsCons2 !Coeffs for 2nd order tensor prods
  real(prec),dimension(d-1,d,2**(q-d-1)+1,2**(q-d-1)+1,N) :: CoeffsCap2	 !Coeffs for 2nd order tensor prods
  real(prec),dimension(d-1,d,2**(q-d-1)+1,2**(q-d-1)+1,N) :: CoeffsLab2	 !Coeffs for 2nd order tensor prods
  
  !----------------------------------------------------
  !Matrices that hold values of chebyshev polynomials 
  !----------------------------------------------------
  real(prec),dimension(2**(q-d)+1,2**(q-d)+1,d) :: T3		! Chebyshev values for 4th order polynomials 
  real(prec),dimension(2**(q-d-1)+1,2**(q-d-1)+1,d) :: T2	! Chebyshev values for 2nd order polynomials
                                       ! and for 2 dimensional tensor products.
  
  !--------------------------------
  !Variables for grid subroutine
  !--------------------------------
  real(prec) :: kss	! deterministic steady state capital stock
  real(prec),dimension(N) :: kmax	! used to put upper bound on capital grid
  real(prec),dimension(N) :: kmin	! used to put lower bound on capital grid

  real(prec) :: lss	! deterministic steady state labor supply
  
  real(prec) :: zmax	! largest positive shock is e^(zmax)
  real(prec) :: zmin	! largest negative shock is e^(zmin)
  
  integer,parameter :: nk=100	! Number of points in capital grid for each country
  integer,parameter :: nz=100	! Number of points in tech. shock grid for each country
  integer,parameter :: nd=100	! Number of points in error checking grid in any dimension
  
  real(prec),dimension(nk,N) :: knorm	! Grid for k -- used to check convergence
  real(prec),dimension(nz,N) :: znorm	! Grid for z
  
  !--------------------------------
  !General global variables
  !--------------------------------
  integer :: kc,zc	!counters for do loops associated with capital stocks or technology shocks
  integer,parameter :: maxiter = 3000	 !Maximum number of iterations
  real(prec),parameter :: tol = 0.000001 !Error tolerance
  real(prec),parameter :: acc = 1.0   !coefficient used for acceleration methods
  integer,dimension(d) :: stind	 ! holds index of state (capital(N), shock(N))
  integer :: iter       !counter for number of policy function iterations
  integer :: iter1	!counter for number of iterations of Nonlinear Equation Solver
  
  real(prec),dimension(N) :: ktemp		!t+1 capital stocks used in EulerEquation and ExpectFOC
  real(prec),dimension(2**(q-d)+1,N)::unitTtemp	!Chebyshev values used in EulerEquation and ExpectFOC
  
  real(prec) :: aggC   !Passes aggregate consumption between PolicyInitial & InitCons
  
  !-------------------------------
  !Variables for Accuracy Check 1
  !-------------------------------
  real(prec),dimension(100,d) :: accpt  
  real(prec),dimension(100,d) :: accptr
  real(prec) :: rad
  integer :: pti
  integer :: ctryj
  
  !-----------------------------------------------
  !Variables for simulation, Accuracy Check 2 & 3
  !-----------------------------------------------
  integer,parameter :: maxsamp = 200, samplen = 1201
  integer :: sn,tc
  real(prec),dimension(maxsamp,samplen,N) :: randnr
  real(prec),dimension(maxsamp,samplen,N) :: shock
  real(prec),dimension(maxsamp,samplen,N) :: SimCapital,SimConsumption,SimOutput,SimInvestment
  real(prec),dimension(maxsamp,samplen,N) :: SimLabor
  real(prec),dimension(N) :: SimTech	!Holds vector of technological shocks along simulated path
  
  !-------------------------------
  !Variables to compute run times
  !-------------------------------
  real(prec) :: time_begin
  real(prec) :: time_end
  real(prec) :: policy_runtime
  real(prec) :: Acc1_runtime, Acc2_runtime, Acc3_runtime
  
  !--------------------------------
  !Variables to store model output
  !--------------------------------
  real(prec), dimension(nk,27) :: VarOut
  real(prec), dimension(1,7) :: Acc1Out
  real(prec), dimension(N,15) :: Acc2Out
  real(prec), dimension(1,3) :: Acc3Out

contains
  
  !*************************
  !** Primative Functions **
  !*************************
  
  !-------------------------------------------------------
  ! function Ucc
  !
  !	2nd Derivative of (CRRA) utility function with
  !	coefficient of relative risk aversion 1/gamma
  !
  ! Inputs:	c -- real, consumption
  !             i -- integer, country index
  ! Returns:	real -- U"(c) 
  !
  !-------------------------------------------------------
  
  function Ucc(c,i)
    implicit none
    real(prec) :: Ucc
    real(prec), intent(in) :: c
    integer, intent(in) :: i
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (c >= epsilon) then
       Ucc = -xgamma(i) * c**(-xgamma(i) - 1.0)
    else
       print*, 'c is below lower bound in the function Ucc, c = ', c	
       !pause
       Ucc = -xgamma(i) * epsilon**(-xgamma(i) - 1.0)
    end if
    
  end function Ucc
  
  !---------------------------------------------------
  ! function Uc
  !
  !	(CRRA) marginal utility function with
  !	coefficient of relative risk aversion 1/gamma
  !
  ! Inputs:	c -- real, consumption 
  !             i -- integer, country index
  ! Returns:	real, marginal utility 
  !
  !---------------------------------------------------
  
  function Uc(c,i)
    implicit none
    real(prec) :: Uc
    real(prec), intent(in) :: c
    integer, intent(in) :: i
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (c >= epsilon) then
       Uc = c**(-xgamma(i))
    else
       print*, 'c is below lower bound in the function Uc, c = ', c
       !pause
       Uc = epsilon**(-xgamma(i)) + Ucc(epsilon,i) * (c - epsilon)
    end if
    
  end function Uc
  
  !-------------------------------------------------------
  ! function Ull
  !
  !	2nd Derivative of (CRRA) utility function wrt
  !	labor
  !
  ! Inputs:	lab -- real, labor supply 
  !             i -- integer, country index
  ! Returns:	real -- U"(lab) 
  !
  !-------------------------------------------------------
  
  function Ull(lab,i)
    implicit none
    real(prec) :: Ull
    real(prec), intent(in) :: lab
    integer, intent(in) :: i
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (lab >= epsilon) then
       Ull = -b(i) * xeta(i) * lab**(xeta(i) - 1.0)
    else
       print*, 'lab below lower bound in the function Ull, lab = ', lab	
       !pause
       Ull = -b(i) * xeta(i) * epsilon**(xeta(i) - 1.0)
    end if
    
  end function Ull
  
  !--------------------------------------------------------
  ! function UlUc
  !
  !	Marginal Rate of Substitution between consumption 
  !     and labor
  !
  ! Inputs:	c -- real, consumption
  !             lab -- real, labor supply
  !             i -- integer, country index
  ! Returns:	real, marginal rate of substitution
  !
  !--------------------------------------------------------
  
  function UlUc(c,lab,i)
    implicit none
    
    real(prec) :: UlUc	
    real(prec), intent(in) :: c
    real(prec), intent(in) :: lab
    integer, intent(in) :: i
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (c >= epsilon .and. lab >= epsilon) then
       UlUc = b(i) * lab**(xeta(i)) * c**(xgamma(i))
    elseif (c < epsilon .and. lab >= epsilon) then
       print*, "Error - consumption below lower bound in UlUc, c = ", c
       pause
       UlUc = b(i) * lab**(xeta(i))  * &
            & (epsilon**(-xgamma(i)) + Ucc(epsilon,i) * (c - epsilon))**(-1.0)
    elseif (c >= epsilon .and. lab < epsilon) then
       print*, "Error - labor below lower bound in UlUc, l = ", lab
       pause
       UlUc = -1.0 * (-b(i) * epsilon**(xeta(i)) + Ull(epsilon,i) * (lab - epsilon)) * &
            & c**(xgamma(i))
    else
       print*, "Error - consumption below lower bound in UlUc, c = ", c
       print*, "Error - labor below lower bound in UlUc, l = ", lab
       pause
       UlUc = -1.0 * (-b(i) * epsilon**(xeta(i)) + Ull(epsilon,i) * (lab - epsilon)) * &
            & (epsilon**(-xgamma(i)) + Ucc(epsilon,i) * (c - epsilon))**(-1.0)
    end if

  end function UlUc

  !------------------------------------------------------------------------------------
  ! function F
  !
  !	Production Function - Takes country-specific capital stock and 
  !	technology shock and returns output.
  !
  ! Inputs:	capital -- real, capital level
  !             lab -- real, labor supply
  !	        sh -- real, technological shock   
  ! Returns:	real, output
  !
  !------------------------------------------------------------------------------------
  
  function F(capital,lab,sh)
    implicit none
    real(prec) :: F
    real(prec), intent(in) :: capital,lab,sh
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (capital >= epsilon .and. lab >= epsilon) then
       F =  A_tfp * exp(sh) * capital**alpha * lab**(1.0-alpha)
    elseif (capital < epsilon .and. lab >= epsilon) then
       print*, "Error - Capital Stock is negative in F, cap = ", capital 
       !pause
       F = 0
    elseif (capital >= epsilon .and. lab < epsilon) then
       print*, "Error - Labor is negative in F, lab = ", lab 
       !pause
       F = 0
    else
       print*, "Error - Capital Stock is negative in F, cap = ", capital 
       print*, "Error - Labor is negative in F, lab = ", lab 
       !pause
       F = 0
    end if
    
  end function F
  
  !----------------------------------------------------------------------------
  ! function Fk
  !
  !	Marginal Product of Capital - Takes country specific capital stock and 
  !	technology shock and returns mpk.
  !
  ! Inputs:	capital -- real, capital level
  !             lab -- real, labor supply
  !		sh -- real, technological shock   
  ! Returns:	real, mpk
  !
  !----------------------------------------------------------------------------
  
  function Fk(capital,lab,sh)
    implicit none
    real(prec) :: Fk
    real(prec), intent(in) :: capital,lab,sh
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (capital >= epsilon .and. lab >= epsilon) then
       Fk =  alpha * A_tfp * exp(sh) * (capital/lab)**(alpha-1.0)
    elseif (capital < epsilon .and. lab >= epsilon) then
       print*, "Error - Capital Stock is negative in Fk, cap = ", capital  
       !pause
       Fk =  alpha * A_tfp * exp(sh) * epsilon**(alpha-1.0) * lab**(1.0-alpha)
    elseif (capital >= epsilon .and. lab < epsilon) then
       print*, "Error - Labor Supply is negative in Fk, lab = ", lab  
       !pause
       Fk =  alpha * A_tfp * exp(sh) * capital**(alpha-1.0) * epsilon**(1.0-alpha)
    else
       print*, "Error - Capital Stock is negative in Fk, cap = ", capital  
       print*, "Error - Labor Supply is negative in Fk, lab = ", lab  
       !pause
       Fk =  alpha * A_tfp * exp(sh) * epsilon**(alpha-1.0) * epsilon**(1.0-alpha)
    end if
    
  end function Fk
  
  !----------------------------------------------------------------------------
  ! function Fl
  !
  !	Marginal Product of Labor - Takes country specific capital stock and 
  !	technology shock and returns mpl.
  !
  ! Inputs:	capital -- real, capital level
  !             lab -- real, labor supply
  !		sh -- real, technological shock   
  ! Returns:	real, mpl
  !
  !----------------------------------------------------------------------------
  
  function Fl(capital,lab,sh)
    implicit none
    real(prec) :: Fl
    real(prec), intent(in) :: capital,lab,sh
    real(prec) :: epsilon
    
    epsilon = (10.0)**(-6)
    
    if (capital >= epsilon .and. lab >= epsilon) then
       Fl = (1.0-alpha) * A_tfp * exp(sh) * (capital/lab)**(alpha) 
    elseif (capital < epsilon .and. lab >= epsilon) then
       print*, "Error - Capital Stock is negative in Fl, cap = ", capital  
       !pause
       Fl = (1.0-alpha) * A_tfp * exp(sh) * epsilon**(alpha) * lab**(-alpha)
    elseif (capital >= epsilon .and. lab < epsilon) then
       print*, "Error - Labor Supply is negative in Fl, lab = ", lab  
       !pause
       Fl = (1.0-alpha) * A_tfp * exp(sh) * capital**(alpha) * epsilon**(-alpha)
    else
       print*, "Error - Capital Stock is negative in Fl, cap = ", capital  
       print*, "Error - Labor Supply is negative in Fl, lab = ", lab  
       !pause
       Fl = (1.0-alpha) * A_tfp * exp(sh) * epsilon**(alpha) * epsilon**(-alpha)
    end if
    
  end function Fl
  
  !-------------------------------------------------------------------------
  ! function AdjCost
  !
  ! Given country-specific captial stocks for today and tomorrow, calculates
  ! Adjustment costs.
  !
  ! Inputs:	captod -- real, capital level today
  !			captom -- real, capital level tomorrow   
  ! Returns:	real, Adjustment Cost
  !
  !-------------------------------------------------------------------------
  
  function AdjCost(ktod, ktom)
    implicit none
    real(prec) :: AdjCost
    real(prec), intent(in) :: ktod,ktom
    real(prec) :: i,j
    real(prec) :: epsilon
    real(prec) :: captod, captom
    
    epsilon = (10.0)**(-6)
    
    captod = ktod
    captom = ktom
    
    if (captod < epsilon) then
       print*, "Error - Today's Capital Stock is negative in AdjCost"        
       !pause
       captod = epsilon
    elseif (captom < epsilon) then
       print*, "Error - Tomorrow's Capital Stock is negative in AdjCost"
       !pause
       captom = epsilon
    else
       j = captom/captod - 1.0
       Adjcost = 0.5*kappa * j * j * captod
    end if
    
  end function AdjCost
  
  !---------------------------------------------------------------------------
  ! function AdjCost_k
  !
  ! Given country-specific captial stocks for today and tomorrow, calculates
  ! derivative of adjustment cost function wrt capital today.
  !
  ! Inputs:	captod -- real, capital level today
  !			captom -- real, capital level tomorrow   
  ! Returns:	real, Derivative of AdjCost wrt captod
  !
  !---------------------------------------------------------------------------
  
  function AdjCost_k(ktod, ktom)
    implicit none
    real(prec) :: AdjCost_k
    real(prec), intent(in) :: ktod,ktom
    real(prec) :: i,j,j1
    real(prec) :: epsilon
    real(prec) :: captod, captom
    
    epsilon = (10.0)**(-6)
    captod = ktod
    captom = ktom
    
    if (captod < epsilon) then
       print*, "Error - Today's Capital Stock is negative in AdjCost_k"
       !pause
       captod = epsilon
    elseif (captom < epsilon) then
       print*, "Error - Tomorrow's Capital Stock is negative in AdjCost_k"
       !pause
       captom = epsilon
    else
       j = captom/captod - 1.0
       j1 = captom/captod + 1.0
       AdjCost_k = (-0.5)*kappa*j*j1
    end if
    
  end function AdjCost_k
  
  !--------------------------------------------------------------------------
  ! function AdjCost_ktom
  !
  ! Given country-specific captial stocks for today and tomorrow, calculates
  ! derivative of adjustment cost function wrt capital tomorrow.
  !
  ! Inputs:	captod -- real, capital level today
  !			captom -- real, capital level tomorrow   
  ! Returns:	real, Derivative of AdjCost wrt captom
  !
  !--------------------------------------------------------------------------
  
  function AdjCost_ktom(ktod, ktom)
    implicit none
    real(prec) :: AdjCost_ktom
    real(prec), intent(in) :: ktod,ktom
    real(prec) :: i,j
    real(prec) :: epsilon
    real(prec) :: captod, captom
    
    epsilon = (10.0)**(-6)
    captod = ktod
    captom = ktom
    
    if (captod < epsilon) then
       print*, "Error - Today's Capital Stock is negative in AdjCost_ktom"
       !pause
       captod = epsilon
    elseif (captom < epsilon) then
       print*, "Error - Tomorrow's Capital Stock is negative in AdjCost_ktom"
       !pause
       captom = epsilon
    else
       j = captom/captod - 1.0
       AdjCost_ktom = kappa * j
    end if

  end function AdjCost_ktom
  
end module Params
