PROGRAM PORTFOLIO


! This fortran code includes:
!    a- Numerical solution to the life cycle problem of consumption and investment in the risky asset market
!       * for given values of: the discount rate, expectation on asset returns, deterministic and stochastic parameters of labor income 
!       * over a double dimensional grid of: the CRRA and fixed participation costs 
!       * over a three-dimensional grid for the following variables: permanent income, ratio of financial wealth to permanent income and transitory income 
!    b- Evaluation of each individual contribution to the pseudo-likelihood and the pseudo-likelihood of different subsamples 


! Autor: Graciela Sanroman
! This version: August, 15 2012



USE IMSL_LIBRARIES


!11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111!
! 1- DEFINITION OF VARIABLES

IMPLICIT NONE


! 1.1 PARAMETERS
CHARACTER*34 :: RUTA = '/home/gsanroman/decon/portfolio15/'

! DECLARE THE SPECIFICATION FOR RETURNS TO BE USED

CHARACTER*26 :: RETURNS_DATA = './INPUTS/RETORNOS_ESP1.TXT'
! 1. SPECIFICATION 1 IS 3 POINTS GAUSS-HERMITE FOR MEAN AND VARIANCE 1991-2011
! 2. SPECIFICATION 2 IS 3 POINTS GAUSS-HERMITE FOR Pelizzon & Weber
! 3. SPECIFICATION 3 IS 5 POINTS GAUSS-HERMITE FOR MEAN AND VARIANCE 1991-2011
! 4. SPECIFICATION 4 IS 5 POINTS FOR MEAN, VARIANCE, SKEWNESS AND KURTOSIS 1991-2011

INTEGER, PARAMETER :: LDR=3  !LDR=5  Number of grid point to the gross return of the risky asset
REAL*8, PARAMETER :: Beta0=0.94d0   !0.78d0  ;!0.8700000000D0; !

!Note: in adittion one cell for the gross return of risk-free asset with probab. 1 is included 


! DECLARE THE EDUCATIONAL GROUP 
!   - 2 elementary school
!   - 3 basic secondary
!   - 4 high school
!   - 5 college
CHARACTER*25 :: IND_DATAA = './INPUTS/DATA_INC_A_5.TXT'
CHARACTER*25 :: IND_DATAB = './INPUTS/DATA_INC_B_5.TXT'
INTEGER, PARAMETER :: E_GROUP = 4    !E_GROUP NUMBER IS EDUCATIONAL GROUP MINUS 1
INTEGER, PARAMETER :: SIZE_E_GROUP = 510
!Samples sizes
! Edu=2  784                  /L.384
! Edu=3  2890                 /L.1770
! Edu=4  2193                 /L.1273
! Edu=5  510                  /L.284


INTEGER, PARAMETER :: G=40, LCCRRA=20
! G is the number of evaluation points for the participation costs / LCCRRA is the number of evaluation points for the CRRA

REAL*8, PARAMETER :: vini_gg=0.0001d0, vini_crra=1.3d0
!vini_gg is the first value of the grid for participation costs, vini_crra is the first value of the grid for CRRA

REAL*8, PARAMETER :: stepgg=0.001d0, stepcrra=0.050d0
! stepgg is the step of the grid for participation costs, stepcrra is the step of the grid for CRRA

INTEGER, PARAMETER :: typec1=1, typec2=4  
! typec is the type of cost: 1. permanent income 2. transitory income 3. income 4. monetary


! 1.2 VARIABLES
Real*8 :: gamma0  !CRRA
 

INTEGER :: IDATAB(9,SIZE_E_GROUP) 
!variables are 
!    id time 
!    age 
!    SOUTH: Regional DUMMY =1 IF SOUTH OR ISLANDS
!    Doccup: OCCUPATION DUMMY =1 if BLUE-COLLAR =0 IF WHITE-COLLAR
!    Young
!    NoHouse 
!    Dw =1: above median wit 
!    Dr1 =1: above median r1

INTEGER, PARAMETER :: t0=1, T=35
INTEGER, PARAMETER :: eretiro=21
! t is the minimun age 0 corresponds to 24-25 years
! T is the maximum age 34 corresponds to 91-92 years

INTEGER, PARAMETER :: J=100
! J is the number of evaluation point in the dimension of Financial Wealth 
INTEGER, PARAMETER :: GHP=7, GHT=5
INTEGER, PARAMETER :: ISI=GHP*(GHT)
 ! GHP: 5  GAUSS_HERMITE POINTS FOR PERMANENT SHOCKS + 2 CELLS FOR PERCENTILES 1 AND 9 OF THE OBSERVED DISTRIBUTION
 ! GHT: 3 GAUSS_HERMITE POINTS FOR TRANSITORY SHOCKS +  1 CELL FOR THE DISASTER EVENT + 1 CELL FOR THE PERCENTILE 99 OF THE OBSERVED DISTRIBUTION


INTEGER, PARAMETER :: FIRSTJJ = 3
! FIRSTJJ IS THE NUMBER OF INITIAL POINTS OF THE W GRID THAT ARE COMPUTED SIMPLY ADDING 0.001 TO THE PREVIOUS POINT (RECALL GRID(1)=0)


INTEGER, PARAMETER :: M=1, N=2, N2=1, MCON=1,  MAXITERA=100
INTEGER, PARAMETER :: LDC=MCON
INTEGER, PARAMETER :: IBTYPE=0
INTEGER ::  ISTAT, NOUT, IPARA, IRTYPE(MCON)
INTEGER :: IPARAM(7)

! SPECIFICATIONS FOR LOCAL VARIABLES TO BE USED BY THE OPTIMIZATION ALGORITHM



INTEGER :: AGE


INTEGER :: pp, pp2, vv, vv2, ii, iip, iit, ii_prime, tt, jj, gg, ww, rar, hh, llll, cccc, ccrra
! variables to be used as indexes for loops or whithin loops

INTEGER :: orden(ISI), OBS1(13), OBS2(15), INDg(LCCRRA), INDcrra
INTEGER :: young, Miip, Miit

REAL*8, PARAMETER :: MinMaxWgrid=10.000000d0
 !MinMaxWgrid is the minimun of the maximum for the W GRID (is expressed as the ratio FWealth/Eincome)


REAL*8, PARAMETER :: PRECISION=1000.00D0
REAL*8, PARAMETER :: min_feli_arg=0.0001d0

!PARTICIPATION COST
REAL*8 :: G_(2,2), g_coef(G,2), llli(4,g+1,LCCRRA,SIZE_E_GROUP), max_pseudo(LCCRRA), g_pseudo(LCCRRA), g2_pseudo(LCCRRA)


!DATA
REAL*8 :: IDATAA(5,SIZE_E_GROUP)  
!variables are 
!    1- coh1: ratio consumption(t)+financial asset at the end of t to labor income(t)
!    2- coh2: ratio labor income(t) + financial asset at the end of t to labor income(t)
!    3- ratio r1: etai/mean(Y by age2) where etai is the average income of i in the panel
!    4- ratio r2: yi(t)/etai
!    5- alpha(t)
!    6- alpha(t-1) (-99 if missing)
   
REAL*8 ::  wit,Alphai, LagAlphai, Z(15), ratio1, ratio2
REAL*8 :: pseudo1_ll(13,G,LCCRRA),  outprod1_ll(13,13,G,LCCRRA), hessian1_ll(13,13,G,LCCRRA)
REAL*8 :: pseudo2_ll(15,G,LCCRRA),  outprod2_ll(15,15,G,LCCRRA), hessian2_ll(15,15,G,LCCRRA)
REAL*8 :: var1(13), var2(15)
Real*8 :: lli, lli0, scorei, DEv1_v0, lant
  

REAL*8 :: IRATIOS(2,SIZE_E_GROUP)  ! 1: RINC ratio INCOME(AGE)/EINCOME(AGE) 2: RAF (AF(AGE)/EINCOME(AGE)

REAL*8 :: BETAV(3,4)
! BETA IS A 4*3 matrix contains:
!            - rows: 1-coefficient of age of (log) labor income process 
!                    2-coefficient of age^2/10 of (log) labor income process
!            - columns; educational groups

REAL*8 :: MINMAXI(3,4) !MIN AND MAX PERMANENT INCOME GRID, MAX TRANSITORY INCOME GRID

REAL*8 :: EINCOME(T), MIT(T)
!EINCOME(T) vector of expected value of labor income at periods t0 to T 
! MIT(t) INCOME deterministic growth rate at each period


REAL*8 :: GRIDW(J,T)
!matrix of evaluation points for the financial wealth at each age-period


REAL*8 :: VALUEFUNCTION(J,GHP,GHT,T), EVALUEFUNCTION(J,T,3)
!VALUEFUNCTION at each evaluation point of financial wealth and permanent and transitory income at each age-period

REAL*8 :: VALUEFUNCTION1(J,GHP,GHT,T)
!VALUEFUNCTION if the agent invests/pay the cost in the risky asset market at the period at each evaluation point of the financial wealth and permanent and transitory income at each age-period

REAL*8 :: VALUEFUNCTION0(J,GHP,GHT,T)
!Matrix with VALUEFUNCTION if the agent does not invest in the risky asset market at each period at each evaluation points for the financial wealth and permanent and transitory income at each age-period


REAL*8 :: ALPHAFUNCTION(J,GHP,GHT,T)
!OPTIMAL RISKY ASSET SHARES at each evaluation point of financial wealth and permanent and transitory income at each age-period

REAL*8 :: ALPHAFUNCTION0(J,GHP,GHT,T)
REAL*8 :: ALPHAFUNCTION1(J,GHP,GHT,T)

REAL*8 :: SAVINGSFUNCTION(J,GHP,GHT,T),  ESAVINGSFUNCTION(J,T)
!OPTIMAL CONSUMPTION at each evaluation point of financial wealth and permanent and transitory income at each age-period
!EXPECTED VALUE OF  at each evaluation point of financial wealth and and permanent and transitory income at each age-period
REAL*8 :: SAVINGSFUNCTION0(J,GHP,GHT,T)
REAL*8 :: SAVINGSFUNCTION1(J,GHP,GHT,T)


REAL*8 :: SIGMA(2,4), GAUSSHERMITEP(GHP-2,2), GAUSSHERMITET(GHT-2,2)
REAL*8 :: PSHOCKS(GHP,2,2), TSHOCKS(GHT,2,2)
! SIGMA IS A 2X4 matrix 
!        - each column educational group
!        - row 1: the variance of the permanent shocks PHI
!        - row 2:  the variance of the transitory shocks V
! GAUSSHERMITEt/p is a GHx2 matrix with
!               - first columm: values at evaluation points
!               - second column: weights
!               - rows: for each point of gauss-hermite
! PSHOCKS  is the GHP,2,2 ARRAY contained values and probabilities for the PERMANENT INCOME SHOCK-
!                 GHP,:,1  the first and final point correspond to the 1 and 99 percentiles of the observed distribution, with probabilities equal 0
!                   the second to the ghp-1 points correspond to the gauss-hermite quadrature 
 !                GHP=1,:,2= 1 IS FOR EVALUATION AFTER AGE 60 ( TT=eretiro) TAKES ONLY VALUE 1 WITH PROB. 1

! TSHOCKS  is the (GHT,2,2 contained values and probabilities for the TRANSITORY INCOME SHOCK-
!         GHT,:,1 1- value 0.05: probability 0.005
!                 2 TO GHT-1  given by Gauss-Hermite algorithm times 0.995
!                 3 GHT max value of TSHOCKS as the 99th percentile of observed distribution in E_group 
!                GHT=1,:,2= 1 IS FOR EVALUATION AFTER AGE 60 ( TT=eretiro) TAKES ONLY VALUE 0.80 WITH PROB. 1 

REAL*8 :: R(LDR+1), PrR(LDR+1), RETURNS(LDR+1,2)
!R(1) GROSS RETURN ON THE RISKYFREE ASSET, 
! R(2), R(3) GROSS RETURN OF THE RISKY ASSET ON THE BAD STATE AND GOOD STATE RESPECTIVELY 
!PrBad= PROBABILITY OF THE BAD STATE 


Real*8 :: PDV(T), insPDV(T), MaxW(T), UMGCPRIME, vPRIME, dvPRIME, nWP, EV_PRIME, W_PRIME, WPRIME_INC, CONS_PRIME(ISI,3)
Real*8 :: I_VPRIME, factor_psi_mit, I_Ve 
! PDV(tt): Present Discounted Value (at rate R(1)) of the future expected stream of labor income
! insPDV(tt): equivalent installment of PDV(tt)
! I_VPRIME INTERPOLATED VALUE FUNCTION
! CONS_PRIME(ISI,3), correspond to consumption at ii_prime=ISI and returns to risky asset in the bad stat rar=2  and good state rar=3 (CONS_PRIME(ISI,1) is not used)

REAL*8 :: BL(MCON), C(MCON,N), RNORM, X(N), XGUESS(N), XLB(N), XUB(N), XX(N), FV !Variables to be used in the optimization algorithm, 
!   X(1) / optimal savings
!   X(2) / optimal risky asset shares

REAL*8, PARAMETER :: XSCALE(N)=1000D0, FSCALE=1.00D0
REAL*8 :: RPARAM(7)

REAL*8 :: IX2  ! = 1 IF X(2) GT 0

REAL*8 :: umgCtt_nosavings, EumgCtt_prime_nosavings, LAMBDA, uctt
    ! umgCtt_nosavings: Marginal utility of consumption at t when no savings
    ! EumgCtt_prime_nosavings Expected value of Marginal utility of consumption at t+1 when no savings (obtained by the numerical derivative of E(valuefunction)
    
REAL*8 :: temp1, temp2, temp3, temp4, temp5, temp6, acum1, acum2, acum3, aux1, aux2, aux3, aux4, W, Y, COHI, INCOMET, v(2), uct, u_c(3), err(2), DifOPT0, DifOPT1, DifOPT2, ch1, ch2,ch3, ch4, auxgrid1, auxgrid2
    
REAL*8 :: v1, xx1, xx2
! Temporary and auxiliar variables


REAL*8 :: time_begin0, time_end0


!2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222!
! DECLARATION OF INPUT FILES

OPEN (UNIT=1, FILE = RUTA//'INPUTS/BETA.TXT')
OPEN (UNIT=2, FILE=  RUTA//'INPUTS/SHOCKSRENTA.TXT')
OPEN (UNIT=31, FILE= RUTA//'INPUTS/GAUSS_HERMITEP.TXT')
OPEN (UNIT=32, FILE= RUTA//'INPUTS/GAUSS_HERMITET.TXT')
OPEN (UNIT=33, FILE= RUTA//'INPUTS/MINMAXY.TXT')
OPEN (UNIT=4, FILE= IND_DATAA )
OPEN (UNIT=5, FILE= IND_DATAB )
OPEN (UNIT=44, FILE= RETURNS_DATA )

! DECLARATION OF OUTPUT FILES
OPEN (UNIT=104, FILE= RUTA//'OUTPUTS/EDU5/INCDATA.TXT')

OPEN (UNIT=105, FILE= RUTA//'OUTPUTS/EDU5/INCSHOCKS.TXT')
OPEN (UNIT=101, FILE= RUTA//'OUTPUTS/EDU5/WGRID.TXT')
OPEN (UNIT=1001, FILE= RUTA//'OUTPUTS/EDU5/VALUEFUNCTION.TXT')
OPEN (UNIT=1002, FILE= RUTA//'OUTPUTS/EDU5/VALUEFUNCTION0.TXT')
OPEN (UNIT=1003, FILE= RUTA//'OUTPUTS/EDU5/VALUEFUNCTION1.TXT')


OPEN (UNIT=1004, FILE= RUTA//'OUTPUTS/EDU5/ALPHAFUNCTION.TXT')
OPEN (UNIT=1005, FILE= RUTA//'OUTPUTS/EDU5/ALPHAFUNCTION0.TXT')
OPEN (UNIT=1006, FILE= RUTA//'OUTPUTS/EDU5/ALPHAFUNCTION1.TXT')


OPEN (UNIT=1007, FILE= RUTA//'OUTPUTS/EDU5/SAVINGSFUNCTION.TXT')
OPEN (UNIT=1008, FILE= RUTA//'OUTPUTS/EDU5/SAVINGSFUNCTION0.TXT')
OPEN (UNIT=1009, FILE= RUTA//'OUTPUTS/EDU5/SAVINGSFUNCTION1.TXT')

OPEN (UNIT=1010, FILE= RUTA//'OUTPUTS/EDU5/EVALUEFUNCTION.TXT')

OPEN (UNIT=508, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCESTYPE0.TXT')
OPEN (UNIT=509, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCESTYPEA.TXT')
OPEN (UNIT=510, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCESTYPEB.TXT')
OPEN (UNIT=511, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCESTYPE1.TXT')
OPEN (UNIT=512, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCESTYPE2.TXT')
OPEN (UNIT=513, FILE= RUTA//'OUTPUTS/EDU5/INCONSISTENCES.TXT')

OPEN (UNIT=13, FILE= RUTA//'OUTPUTS/EDU5/controles.TXT')
OPEN (UNIT=1313, FILE= RUTA//'OUTPUTS/EDU5/residuos.TXT')

OPEN (UNIT=231313, FILE= RUTA//'OUTPUTS2/EDU5/OBSERVATIONS.TXT')
OPEN (UNIT=331313, FILE= RUTA//'OUTPUTS2/EDU5/GROUPANDTYPE.TXT')

OPEN (UNIT=1313131, FILE= RUTA//'OUTPUTS2/EDU5/pseudo1_loglik1.TXT')
OPEN (UNIT=1313132, FILE= RUTA//'OUTPUTS2/EDU5/pseudo1_loglik2.TXT')
OPEN (UNIT=1313133, FILE= RUTA//'OUTPUTS2/EDU5/pseudo1_loglik3.TXT')
OPEN (UNIT=1313134, FILE= RUTA//'OUTPUTS2/EDU5/pseudo1_loglik4.TXT')

OPEN (UNIT=2313131, FILE= RUTA//'OUTPUTS2/EDU5/pseudo2_loglik1.TXT')
OPEN (UNIT=2313132, FILE= RUTA//'OUTPUTS2/EDU5/pseudo2_loglik2.TXT')
OPEN (UNIT=2313133, FILE= RUTA//'OUTPUTS2/EDU5/pseudo2_loglik3.TXT')
OPEN (UNIT=2313134, FILE= RUTA//'OUTPUTS2/EDU5/pseudo2_loglik4.TXT')

OPEN (UNIT=4313131, FILE= RUTA//'OUTPUTS2/EDU5/pseudo_loglik1_IND_1.TXT')
OPEN (UNIT=4313132, FILE= RUTA//'OUTPUTS2/EDU5/pseudo_loglik2_IND_2.TXT')
OPEN (UNIT=4313133, FILE= RUTA//'OUTPUTS2/EDU5/pseudo_loglik3_IND_3.TXT')
OPEN (UNIT=4313134, FILE= RUTA//'OUTPUTS2/EDU5/pseudo_loglik4_IND_4.TXT')

!333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333!
! DECLARATION OF FORMATS
101 FORMAT ( <40> f12.5)
102 FORMAT ( <1> i6, <4> f12.5)
103 FORMAT ( <1> i6, <J+1> f15.11)
104 FORMAT ( <2> i6, <J+1> f22.5)
1014 FORMAT ( <2> i6, <2> f12.5, <40> f25.4)
105 FORMAT ( <2> A12, <J+1> f12.5)
106 FORMAT ( <2> i12, <ISI+1> f25.11)
107 FORMAT ( <1> a12, <J+1> i12)
108 FORMAT ( <3> A12, <J+1> f15.11)
109 FORMAT ( <20> a15)
110 FORMAT ( <4> i12, <J+1> f15.5)
111 FORMAT ( <3> i12, <8> f12.5)
1331 FORMAT ( <5> i12, <8> f15.11)
1332 FORMAT ( <11> i25, <30> f25.11)
1313 FORMAT ( <20> i12)
1314 FORMAT ( <10> a25, <G+1> f25.11 )

CALL CPU_TIME(time_begin0) 

!444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444!
! READ INPUT FILES, VALUES OF OTHER PARAMETERES


READ (1,*) BETAV
CLOSE (1, STATUS="KEEP")
!print *, "BETA", BETAV(4,3)
READ (2,*) SIGMA
!print *, "SIGMA", SIGMA(1,E_GROUP), SIGMA(2,E_GROUP)
CLOSE (2, STATUS="KEEP")
READ (31,*) GAUSSHERMITEP
!print *, "GAUSSP", GAUSSHERMITEP
CLOSE (3, STATUS="KEEP")

READ (32,*) GAUSSHERMITET
!print *, "GAUSST", GAUSSHERMITET
CLOSE (32, STATUS="KEEP")

READ (33,*) MINMAXI
CLOSE (33, STATUS="KEEP")
!print *, "MINMAXI", MINMAXI

READ (5,*) IDATAB
CLOSE (5, STATUS="KEEP")

READ (4,*) IDATAA
CLOSE (4, STATUS="KEEP")


READ (44,*) RETURNS
CLOSE (44, STATUS="KEEP")

!FINANCIAL ASSETS RETURNS


! R(1) GROSS RETURN ON THE RISKYFREE ASSET,
! R(2) a R(Ldr+1) and PrR(2) a PrR(Ldr+1) GROSS RETURN and respectives probabilities OF THE RISKY ASSET 
! ON THE BAD STATE AND
do rar=1,LDR+1
R(rar)=RETURNS(rar,1)
PrR(rar)=RETURNS(rar,2)
enddo


!66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666
!6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\6.1\\

!First line of outputs VALUEFUNCTION, ALPHAFUNCTION AND CONSUMPTION FUNCTION
do ww = 1001,1009
WRITE (ww, 108) "tt", "JJ", "GRIDW"  
enddo

WRITE (1010, 109) "tt", "JJ", "EVALUEFUNCTION", "EVALUEFUNCTION0", "EVALUEFUNCTION1"
WRITE (508, 109) "jj", "iip", "iit", "t", "DifOPT", "VF(j,ip,it,t)", "VF(j,ip,it-1,t)", "SF(j,ip,it,t)", "AF(j,ip,it,t)" 
WRITE (509, 109) "jj", "iip", "iit", "t", "DifOPT", "VF(j,ip,it,t)", "VF(j,ip-1,it,t)", "SF(j,ip,it,t)", "AF(j,ip,it,t)" 
WRITE (510, 109) "jj", "iip", "iit", "t", "DifOPT", "VF(j,ip,it,t)", "VF(j-1,ip,it,t)", "SF(j,it,ip,t)", "AF(j,ip,it,t)" 
WRITE (511, 109) "jj", "iip", "iit", "t", "DifOPT", "VF(j,ip,it,t)", "VF0(j,ip,it,t)", "SF(j,ip,it,t)", "SF0(j,ip,it,t)" 
WRITE (512, 109) "jj", "iip", "iit", "t", "DifOPT", "VF(j,ip,it,t)", "VF1(j,ip,it,t)", "SF(j,ip,it,t)", "SF1(j,ip,it,t)", "AF(j,ip,it,t)", "AF1(j,ip,it,t)"  
WRITE (513, 109) "jj", "iip", "iit", "t", "DifOPT", "AF(j,ip,it,t)", "VF0(j,ip,it,t)" , "VF(j,ip,it,t)", "VF1(j,ip,it,t)", "SF(j,ip,it,t)", "SF1(j,ip,it,t)" 
!
!WRITE (1313, 109) "t", "jj", "iip", "iit", "itera", "x1", "x2", "d1", "d2"


!
!do hh=2344,2349
!print 111,  idatab(:,hh), idataa(:,hh)
!pause
!enddo


!PREFERENCE PARAMETERS
DO CCRRA=1,LCCRRA
    print *, "lccrra", CCRRA
    !CRRA

    Gamma0=vini_crra+(stepcrra*(CCRRA*1.D0-1.D0))
    Lambda=DEXP( (1.000D0/GAMMA0)*DLOG(Beta0*R(1))-DLOG(R(1)) )



!55555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555!
! COMPUTE ALL THE VARIABLES RELATED WITH THE NON-FINANCIAL INCOME PROCESS

!5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//
! VALUES OF SHOCKS AND CORRESPONDING PROBABILITIES

!5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//5.1//
! VALUES OF SHOCKS AND CORRESPONDING PROBABILITIES


! points of the Gauss-Hermite quadrature

             PSHOCKS(1,1,1)=MINMAXI(1,E_GROUP)
             PSHOCKS(1,2,1)=0.00D0
             write (105, *),  "pShocks", PSHOCKS(1,:,1)
             DO pp = 2,GHP-1
             pp2=pp-1
      	                  PSHOCKS(pp,1,1)=dexp(dsqrt(2d0*(SIGMA(1,E_GROUP)))*GAUSSHERMITEP(pp2,1))/dexp(0.5d0*(SIGMA(1,E_GROUP)))
      	                  PSHOCKS(pp,2,1)=GAUSSHERMITEP(pp2,2)
                        write (105, *),  "pShocks", PSHOCKS(pp,:,1)
             ENDDO !PP

             PSHOCKS(GHP,1,1)=MINMAXI(2,E_GROUP)
             PSHOCKS(GHP,2,1)=0.00D0
             write (105, *),  "pShocks", PSHOCKS(GHP,:,1)
             
    
         
             
               
                            
                 ! 0.05 event with probability 0.005
                  TSHOCKS(1,1,1)=0.05d0
                  TSHOCKS(1,2,1)=0.005d0   
                 
                           write (105, *),  "TShocks",1, TSHOCKS(1,:,1)   
                          
             DO vv = 2,GHT-1
                          vv2=vv-1
      	                  TSHOCKS(vv,1,1)=dexp(dsqrt(2d0*(SIGMA(2,E_GROUP)))*GAUSSHERMITET(vv2,1))/dexp(0.5d0*(SIGMA(2,E_GROUP)))
      	                  TSHOCKS(vv,2,1)=GAUSSHERMITET(vv2,2)*(1.000000000D0-TSHOCKS(1,2,1))
      	                 write (105, *),  "TShocks", VV, TSHOCKS(vv,:,1)  
            ENDDO    
            
                  TSHOCKS(GHT,1,1)=MINMAXI(3,E_GROUP)
                  TSHOCKS(GHT,2,1)=0.00d0  
                    write (105, *),  "TShocks", VV, TSHOCKS(GHT,:,1)   
    !FOR AGE +eretiro    
      PSHOCKS(1,:,2)=1.0d0
      TSHOCKS(1,1,2)=0.80d0; TSHOCKS(1,2,2)=1.0d0 

!5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//5.2//
! Expected value of labor income, PDV of the expected labor income,  equivalent installment of PDV



!Expected value for each period labor income at the initial period
  EINCOME(t0)=1.0000D0 
  MIT(t0)=1.000D0
 !z(it)=Beta1*age(it)+Beta2*age2(it)+Beta3+Beta4*t
DO tt = t0+1, eretiro
  EINCOME(tt)=EINCOME(tt-1)* dexp(BETAV(1 ,E_GROUP)*(tt*1.0d0))* dexp(BETAV(2,E_GROUP)*((tt*1.0d0)**2.0d0))
  MIT(tt)= dexp(BETAV(1 ,E_GROUP)*(tt*1.0d0))* dexp(BETAV(2,E_GROUP)*((tt*1.0d0)**2.0d0))
ENDDO

temp1=0.000D0

do tt=1,5
temp1=temp1+EINCOME(eretiro+1-tt)/5.0d0
enddo
  
  EINCOME(eretiro+1:T)=temp1*0.8d0
  MIT(eretiro+1:T)=1.0000D0
  EINCOME(:)=EINCOME(:) !*dexp(BETAV(4,E_GROUP))

!
!PRINT 101, EINCOME
!pause

!!!RATIOS 1 Y 2 INCOME AND FINANCIAL WEALTH DIVIDED EXPECTED INCOME FOR CORRESPONDING AGES

!do ii=1,SIZE_E_GROUP
!   pp=IDATAA(3,ii)
!   print *, "pp", pp
!   IRATIOS(1,ii)=dexp(IDATAB(1,ii)*1.0d0)/EINCOME(pp)
!   IRATIOS(2,ii)=IDATAB(2,ii)/EINCOME(pp)
!   WRITE (104, 111), IDATAA(:,II), IDATAB(:,II),IRATIOS(:,II)
!enddo

! PDV(tt): Present Discounted Value (at rate R(1)) of the future expected stream of non-financial income
! insPDV(tt): equivalent installment of PDV(tt)
PDV(T)=EINCOME(T)
insPDV(T)=EINCOME(T)
tt=T-1
    DO WHILE (  tt .ge. 1)
     PDV(tt)=EINCOME(tt)+PDV(tt+1)/R(1)
      temp1=R(1)**(-(T-tt+1))
        temp2=(T-tt+1)*1.0d0
        insPDV(tt)=PDV(tt)*(LAMBDA-1.000000000D0)/(LAMBDA**temp2-1.000000000D0)  
!             print *, PDV(tt), insPDV(tt)
             tt=tt-1
            
    ENDDO

do tt=1,T
  WRITE (104, 102), tt, EINCOME(tt), PDV(tt), insPDV(tt)
ENDDO

!5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//5.3//
! FINANCIAL WEALTH GRIDS


! Max Value for the Financial Wealth at each period
MaxW(1)=MinMaxWgrid
 !print 102, 1, EINCOME(1), PDV(1), insPDV(1), MaxW(1)

DO tt = 2,T
temp1=MaxW(tt-1)*EINCOME(tt-1)*R(1) + EINCOME(tt)-insPDV(tt)/1.5d0
temp2=temp1/Eincome(tt)
MaxW(tt)=max(MinMaxWgrid,temp2)
!  print 102, tt, EINCOME(tt), insPDV(tt), MaxW(tt)
ENDDO
!maxw(:)=1.000000000D0
!pause
! FINANCIAL WEALTH GRIDS for each period
DO tt = 1,T
  GRIDW(1,tt)=0.00000000000d0
  
  do jj=2,FIRSTJJ
    GRIDW(jj,tt)=GRIDW(jj-1,tt)+0.001000d0
  ENDDO
 
  
  temp1=J*1.00000000D0-FIRSTJJ*1.00000000D0
  DO jj = FIRSTJJ+1 , J
     temp2=(jj*1.00000000D0-FIRSTJJ*1.00000000D0)
     GRIDW(jj,tt)=GRIDW(FIRSTJJ,tt)+temp2*temp2/(temp1*temp1)*(MaxW(tt)-GRIDW(FIRSTJJ,tt))
  ENDDO !jj
  
  
  WRITE (101, 103) tt, GRIDW(:,tt)
END DO !tt


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777!
! SOLUTION TO THE ECONOMIC PROBLEM

!7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//
! PARTICIPATION COST
 do cccc=typec1,typec2  !COST TYPE  

   do gg = 1, G  !VALUE OF G_COEF
       print *, "GG", GG
                    g_coef(gg,1)=vini_gg+(gg*1.d0-1.d0)*stepgg
                     g_coef(gg,2)=0.00D0


                    !PARTICIPATION COSTS
                    !   Type 1     G_(1,2)=0d0 !g_coef(gg,1)*EINCOME(T-1)
                    !        G_(1,1)=g_coef(1,1)
                    !        G_(2,2)=0D0
                    !        G_(2,1)=0D0
       

        do llll=1,1 !VALUE OF smoothing
           !lant=0.001
           lant=0.001d0+0.0005d0*(llll*1.d0-1.d0) !SMOOTHING PARAMETER
          
        
        
                print *, "CCCC", CCCC
                  if ( gg .eq. 1) then; WRITE (331313, 1313), E_GROUP, cccc; else; endif
        
                   

  
call erset(0,0,0) 




!7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//7.0//
! PROBLEM AT PERIOD T

! Expected value of utility from consumption at each point of the FINANCIAL GRID
             jj=J
             DO WHILE (jj .ge.1)
                    
                    
                    !Value function at each point, jj, iip, iit
                    
                    
                    !Expected Value Function at each jj
                         temp2= GRIDW(jj,T)+TSHOCKS(1,1,2)
		                 EVALUEFUNCTION(jj,T,1)=feli(temp2,gamma0)
!                   WRITE (1010, 106), T, jj, GRIDW(jj,T), EVALUEFUNCTION(jj,T,1)
              jj=jj-1
            ENDDO !jj
            
            
 
!7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//7.1//          
!SOLUTION TO PERIOD tt 


                  
	      tt=T-1    
             DO WHILE ( tt .ge. 1 )
                   print *, "TT", TT
              if ( tt .gt. eretiro ) then; young=2; Miip=1; Miit=1; else; young=1; Miip=GHP; Miit=GHT;; endif                 
	          jj=J
              DO WHILE ( jj .ge. 1 )
             
               W=GRIDW(jj,tt)
	             
                  DO iip=1,Miip
                     factor_psi_mit=dexp((1.000d0-gamma0)*dlog(PSHOCKS(iip,1,young)*MIT(tt)))
                          
                          DO iit=1,Miit
                             !Participation costs
                             if (cccc .eq. 1) then; G_(1,1)=g_coef(gg,1); else; endif
                             if (cccc .eq. 2) then; G_(1,1)=g_coef(gg,1)*TSHOCKS(iit,1,young)/PSHOCKS(iip,1,young); else; endif
                             if (cccc .eq. 3) then; G_(1,1)=g_coef(gg,1)*TSHOCKS(iit,1,young); else; endif
                             if (cccc .eq. 4) then; G_(1,1)=g_coef(gg,1)/PSHOCKS(iip,1,young); else; endif
                             !Normalized cash-on-hand
                             COHi=((W/(PSHOCKS(iip,1,young)*MIT(tt))+TSHOCKS(iit,1,young)))                   
                              
                                !7.1.a OBVIOUS CORNER SOLUTION NO SAVINGS, ALPHA=0
						        !Mg utility Consumption at t, no savings at t
						           temp1=COHi-0.001d0
                                   temp2=COHi
                                   umgCtt_nosavings=(feli(temp2,Gamma0)-feli(temp1,Gamma0))/(0.001d0)
						       
						       !DERIVATIVE OF THE EXPECTED (at t) VALUE FUNCTION (tt+1) no savings at t  
						                         
                                                 temp1=EVALUEFUNCTION(1,tt+1,1)
                                                 
                                                 CALL INTERPOLAV(tt,0.001d0,vPRIME)
                                                 temp2=vPRIME
                                                 EumgCtt_prime_nosavings=Beta0*(temp2-temp1)/0.001d0
                            
		       	                      if (( umgCtt_nosavings .gt. EumgCtt_prime_nosavings) ) then   ! OBVIOUS CORNER SOLUTION NO SAVINGS
							   	        X(1)=0d0
							   	        X(2)=0d0
							   	        v(1)=feli(COHi,gamma0)+Beta0*EVALUEFUNCTION(1,tt+1,1)
							   	        
							   	        SAVINGSFUNCTION0(jj,iip,iit,tt)=x(1)
							   	        SAVINGSFUNCTION(jj,iip,iit,tt)=x(1)
							   	        SAVINGSFUNCTION1(jj,iip,iit,tt)=x(1)
							   	        ALPHAFUNCTION1(jj,iip,iit,tt)=x(2)
							   	      
							   	        VALUEFUNCTION0(jj,iip,iit,tt)=factor_psi_mit*v(1)
							   	        VALUEFUNCTION(jj,iip,iit,tt)=factor_psi_mit*v(1)
							   	           temp1=max(COHi-G_(1,1), 0.0d0)
							   	        v(1)=feli(temp1,gamma0)+Beta0*EVALUEFUNCTION(1,tt+1,1)
			                            VALUEFUNCTION1(jj,iip,iit,tt)=factor_psi_mit*v(1) 
			                            
			                          else; !NON OBVIOUS CORNER SOLUTION NO SAVINGS  
			                  
			                        !7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\7.1b\\
                                    !PARTB. SOLUTION AT CORNER ALPHA=0 (NO RISKY ASSET) TO PERIOD tt 
			                                 
			                                          !DEFINITION OF CONSTRAINTS 
			                                           CALL BOUNDSTYPE3  ! BORROWING CONSTRAINED 0<= x(1) (AND X(1)<=0.95 NON RESTRICTIVE)
			                                           
			                                           !Set initial values
                                                       IF ( JJ .EQ. J) THEN
                                                        Xguess(1)=0.50d0
!                                                       print *, "X(1)", X(1)
!                                                       pause
                                                       ELSE
                                                        Xguess(1)=SAVINGSFUNCTION0(jj+1,iip,iit,tt)
                                                       ENDIF
                                                           
                                                       !CALL OPTIMIZATION SUBROUTINE
!                                                           print *, "7.1 Vo"
		                                                      CALL DBCONF (PROBLEM_tt_b, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
                                                        
                                                            
				                                               call VALUEF(tt,cohi,G_(1,1),X,I_Ve)
!                                                            print *, cohi, x(1), I_Ve, -Fv/precision, iparam(3), iparam(4)
!                                                            pause;
                                                            !computed from consumption and its utility
                                                            v(1)=- factor_psi_mit*Fv/precision    
                                                       
                                                       
                                                            SAVINGSFUNCTION0(jj,iip,iit,tt)=x(1)
							   	        	   	            VALUEFUNCTION0(jj,iip,iit,tt)=v(1)
                                                    
                                                        
                                                        !COMPUTING Value function AT PERIOD tt
                                                               
                                                        
							   	           
				   	                   !7.1c  PART c- OBVIOUS CORNER SOLUTION ALPHA=0
				   	                   
				   	                   !EXPECTED GROSS RETURN OF THE RISKY ASSET
				   	                   temp1=0.000d0
				   	                   DO rar=2,LDR+1    
                                       temp1=temp1+R(rar)*PrR(rar)
                                       ENDDO !RR
                                       
                                       !EXPECTED GROSS RETURN IF INVEST ALL CASH ON HAND IN THE RISKY/RISKLESS ASSET (IMPLIES NO CONSUMPTION AT tt)
                                       temp3=(COHi-G_(1,1))*temp1
                                       temp2=COHi*R(1)
                                      
                                       if ( temp3 .lt. temp2 ) then 
                                              
                                              SAVINGSFUNCTION(jj,iip,iit,tt)=SAVINGSFUNCTION0(jj,iip,iit,tt)
				   	        	   	          VALUEFUNCTION(jj,iip,iit,tt)=VALUEFUNCTION0(jj,iip,iit,tt)
				   	        	   	          ALPHAFUNCTION(jj,iip,iit,tt)=0.000D0
	                                     
                                       else; 
                                                   
                                                     !7.1d\ PART D. SOLUTION TO PERIOD T-1 IF NON OBVIOUS CORNER SOLUTIONS
                                    
                                                     !DEFINITION OF CONSTRAINTS 
		                                                
                !			                         !Set initial values
                                                                  !X(1)
                                                                  XGUESS(1)=SAVINGSFUNCTION0(jj,iip,iit,tt)
                                                                  !X(2)
                                                                  IF ( JJ .EQ. J) THEN
                                                                  XGUESS(2)=0.600D0
                                                                  ELSE
                                                                  XGUESS(2)=ALPHAFUNCTION(jj+1,iip,iit,tt)
                                                                  ENDIF
                                                            
                                                            CALL BOUNDSTYPE0   ! The agent can decide whether to invest in R.A s.t. borrowing constraints                          
			                                                             ! 0<= x(1)<=0.95
			                                                             ! 0<= x(2)<=1.000100D0                   
                                                            !CALL OPTIMIZATION SUBROUTINE
!                                                                     print *, "7.1d 1"  
                                                                     temp1=0.000001d0
!                                                                      CALL DBCPOL (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, temp1, MAXITERA,X, FV)
              	                                                      CALL DBCONF (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
!		                                                             print *, "7.1d 2", iparam(3) 
                                                                       v(1)=- factor_psi_mit*FV/PRECISION
                                                                        call VALUEF(tt,cohi,G_(1,1),X,I_Ve)
!                                                                      print *, "X", x(1), x(2), -fv/PRECISION, iparam(3), iparam(4)
!		                                                              pause
                                                                    
                                                                      ALPHAFUNCTION(jj,iip,iit,tt)=X(2)
                                                                      SAVINGSFUNCTION(jj,iip,iit,tt)=X(1)
				   	        	   	                                  VALUEFUNCTION(jj,iip,iit,tt)=V(1)
!				   	        	   	                                print *, "7.1d  3"               
	                                                                   !new check corner X2=0
	                                                                   IPARA=0
	                                                                   DO WHILE ( (IPARA .Le. 12) .AND. ( VALUEFUNCTION(jj,iip,iit,tt) .lt. VALUEFUNCTION0(jj,iip,iit,tt) ) )	                                                                  
!	                                                                              
!	                                                                               print *, ipara, TT, JJ, IIP, IIT
!	                                                                              print *, savingsfunction0(jj,iip,iit,tt), savingsfunction(jj,iip,iit,tt)
!	                                                                              print *, alphafunction(jj,iip,iit,tt)     
!	                                                                              print *, valuefunction0(jj,iip,iit,tt), valuefunction(jj,iip,iit,tt)
!                                                                                  pause
                                                                                 
                                                                                 ! INITIAL VALUES AT CORNER SOLUTION ALPHA=0
                                                                                  !X(1)
                                                                                  XGUESS(1)=MAX(SAVINGSFUNCTION0(jj,iip,iit,tt)+(IPARA*1.0D0)*0.01D0,0.0D0)
                                                                                  !X(2)
                                                                                  XGUESS(2)=0.00D0+(IPARA*1.0D0)*0.04D0
                                                                                     !CALL OPTIMIZATION SUBROUTINE
                                                                                      CALL BOUNDSTYPE0
                                                                                      CALL DBCONF (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
                		                                                                 
                		                                                               temp4=- factor_psi_mit*FV/PRECISION
                		                                                             
                		                                                            
                		                                                               XX(1)=X(1)
                		                                                               XX(2)=X(2)
                                                                                   
                                                                                  ! INITIAL VALUES AT CORNER SOLUTION ALPHA=1 
                                                                                    !X(1)
                                                                                     XGUESS(1)=MAX(SAVINGSFUNCTION(jj,iip,iit,tt)-(IPARA*1.0D0)*0.01D0,0.0D0)
                                                                                     !X(2)
                                                                                     XGUESS(2)=1.000D0-(IPARA*1.0D0)*0.04D0
                                                                                    !CALL OPTIMIZATION SUBROUTINE
                                                                                      CALL BOUNDSTYPE0
                                                                                      CALL DBCONF (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
                		                                                              
                		                                                               temp5=- factor_psi_mit*FV/PRECISION
!                                                                                  print *, savingfunction0(jj,iip,iit,tt), xx(1), x(1), xx(2), x(2)
!                                                                                  print *, VALUEFUNCTION0(jj,iip,iit,tt), temp4, temp5
!                                                                                  pause
                                                                                  If ( temp4 .ge. temp5) then
                                                                                      ALPHAFUNCTION(jj,iip,iit,tt)=XX(2)
                                                                                      SAVINGSFUNCTION(jj,iip,iit,tt)=XX(1)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=temp4
				   	        	   	                                              else
				   	        	   	                                                  ALPHAFUNCTION(jj,iip,iit,tt)=X(2)
                                                                                      SAVINGSFUNCTION(jj,iip,iit,tt)=X(1)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=temp5
				   	        	   	                                              endif
				   	        	   	                                          
				   	        	   	                                              if ( (IPARA .eq. 12) .AND. ( VALUEFUNCTION(jj,iip,iit,tt) .lt. VALUEFUNCTION0(jj,iip,iit,tt) ) ) then
				   	        	   	                                                  SAVINGSFUNCTION(jj,iip,iit,tt)=SAVINGSFUNCTION0(jj,iip,iit,tt)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=VALUEFUNCTION0(jj,iip,iit,tt)
				   	        	   	                                                  ALPHAFUNCTION(jj,iip,iit,tt)=0.000D0
				   	        	   	                                              else; endif
				   	        	   	                                    
				   	        	   	                                        ipara=ipara+1
				   	        	   	                                    ENDDO
!				   	        	   	                                    print *, "7.1d post-check "
!                                                                      
!                                                                      !check corners X2=0 X2=1, given x1 gt 0.01
!	                                                                   if ( (x(1) .gt. 0.01d0) .and. ( x(2) .lt. 0.01d0 ) .or. ( x(2) .gt. 0.999d0 ) ) then
!	                                                                           !X2=0
!	                                                                           XX(1)=x(1)
!                                                                               XX(2)= 0.000d0
!                                                                               call VALUEF(tt,cohi,G_(1,1),XX,I_Ve)
!                                                                               v(1)= factor_psi_mit*I_Ve
!                                                                               if ( VALUEFUNCTION(jj,iip,iit,tt) .lt. v(1) ) then; 
!                                                                                       PRINT *, "X2=0",  x(1), VALUEFUNCTION(jj,iip,iit,tt), v(1)
!!!                                                                                       pause
!                                                                                     ALPHAFUNCTION(jj,iip,iit,tt)=0.00d0
!	                                                                                 VALUEFUNCTION(jj,iip,iit,tt)=v(1) 
!	                                                                                 ISTAT=1 
!	                                                                                 else; endif
!    	                                                                       
!	                                                                           !X2=1
!	                                                                           XX(1)=x(1)
!                                                                               XX(2)= 1.000d0
!                                                                               call VALUEF(tt,cohi,G_(1,1),XX,I_Ve)
!                                                                               v(1)= factor_psi_mit*I_Ve      
!	                                                                           if ( VALUEFUNCTION(jj,iip,iit,tt) .lt. v(1) ) then;   
!	                                                                                 PRINT *, "X2=1",  x(1), x(2), VALUEFUNCTION(jj,iip,iit,tt), v(1)
!!                                                                                     pause
!	                                                                                 ALPHAFUNCTION(jj,iip,iit,tt)=1.00d0
!	                                                                                 VALUEFUNCTION(jj,iip,iit,tt)=v(1) 
!	                                                                                 ISTAT=1
!	                                                                           else; endif      
!	                                                                     else; endif; !check corners X2
                                                             
!                                                           
!                                                     write (1313, 1331) tt, jj, iip, iit, IPARA, x(1), x(2), err(1), err(2) 
!                                         
                                                
!				                                        !           ATYPE=1
!                                                print *, "AUX4", factor_psi_mit, uctt, ev_prime
!                                                print *, "v(1)", v(1), X(1), x(2)
!                                                pause
            
      endif  ! NON OBVIOUS CORNER SOLUTIONS NO SAVINGS OR ALPHA=0
                  
                                      !7.1e// SOLUTION TO PERIOD tt , "
	                                  !"COMPULSORY INVESTMENST IN THE RISKY ASSETS
	                                  !THE AGENT HAVE TO PAY PARTICIPATION COSTS
                        	          ! PART A: IF OPTIMAL ALPHA >0 THEN 
                                      if ( ALPHAFUNCTION(jj,iip,iit,tt) .gt. 0.00100D0 ) then 
                                              !If unrestricted solution is at x(2)>0, restricted value and savings function are given by unrestricted solution 
                                                  SAVINGSFUNCTION1(jj,iip,iit,tt)=SAVINGSFUNCTION(jj,iip,iit,tt)
                                                  ALPHAFUNCTION1(jj,iip,iit,tt)=ALPHAFUNCTION(jj,iip,iit,tt)
	                                              VALUEFUNCTION1(jj,iip,iit,tt)=VALUEFUNCTION(jj,iip,iit,tt) 
!	                                              print *, "VF1", VALUEFUNCTION1(jj,iiP,IIT,TT)
                                    
                                       else
                                                   
                                                
!                                                          
                                                               !CALL THE SUBROUTINE TO FIND THE OPTIMAL DECISION IN TERMS OF X(1) AND X(2)
                                                               !FOR THOSE AGENTS FOR WHOM IT IS NON OPTIMAL TO INVEST IN RISKY ASSET IN ORDER TO OBTAIN THE CORRESPONDING VALUE FUNCTION
                                                               !NOTICE THAT THE AGENT HAVE TO PAY THE PARTICIPATION COST, OPTIMAL SAVINGS AND INVESTING IN R.A ARE OBTAINED (NOTICE THAT ALPHA CAN BE 0, AND ALSO SAVINGS), AND 
                                                               !DEFINITION OF CONSTRAINTS 
		                                                          CALL BOUNDSTYPE2   ! The agent can decide whether to invest in R.A s.t. borrowing constraints                          
			                                                                         ! 0<= x(1)<=0.95
			                                                                       ! 0.001<= x(2)<=1.00000D0        

                                                 !Set initial values
                                                                  !X(1)
                                                                  XGUESS(1)=SAVINGSFUNCTION0(jj,iip,iit,tt)
                                                                  !X(2)
                                                                  IF ( JJ .EQ. J) THEN
                                                                  XGUESS(2)=0.600D0
                                                                  ELSE
                                                                  XGUESS(2)=ALPHAFUNCTION(jj+1,iip,iit,tt)
                                                                  ENDIF
                                                              
                                                            !CALL OPTIMIZATION SUBROUTINE
                                                                    
              	                                                      CALL DBCONF (PROBLEM_tt_c, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
		                                                         
                                                                       v(1)=- factor_psi_mit*FV/PRECISION
                                                                        call VALUEF(tt,cohi,G_(1,1),X,I_Ve)
!                                                                      print *, "XMAS", x(1), x(2), -fv/PRECISION, iparam(3), iparam(4)
!		                                                              pause
                                                                      SAVINGSFUNCTION1(jj,iip,iit,tt)=X(1)
                                                                      ALPHAFUNCTION1(jj,iip,iit,tt)=X(2)
	                                                                  VALUEFUNCTION1(jj,iip,iit,tt)=v(1) 
!	                                                                     print *, "7.1e"
                                             ENDIF
                                             
                                       
      
    endif  ! NON OBVIOUS CORNER SOLUTION NO SAVINGS        

             
                                             !new check  X2>0
	                                                                   IPARA=0
	                                                                   DO WHILE ( (IPARA .Le. 12) .AND. ( VALUEFUNCTION(jj,iip,iit,tt) .lt. VALUEFUNCTION1(jj,iip,iit,tt) ) )	                                                                  
!	                                                                              
!	                                                                               print *, ipara, TT, JJ, IIP, IIT
!	                                                                              print *, savingsfunction0(jj,iip,iit,tt), savingsfunction(jj,iip,iit,tt)
!	                                                                              print *, alphafunction(jj,iip,iit,tt)     
!	                                                                              print *, valuefunction0(jj,iip,iit,tt), valuefunction(jj,iip,iit,tt)
!                                                                                  pause
                                                                                 
                                                                                 ! INITIAL VALUES AT CORNER SOLUTION ALPHA=0
                                                                                  !X(1)
                                                                                  XGUESS(1)=MAX(SAVINGSFUNCTION1(jj,iip,iit,tt)+(IPARA*1.0D0)*0.01D0,0.0D0)
                                                                                  !X(2)
                                                                                  XGUESS(2)=ALPHAFUNCTION1(jj,iip,iit,tt)+(IPARA*1.0D0)*0.04D0
                                                                                     !CALL OPTIMIZATION SUBROUTINE
                                                                                      CALL BOUNDSTYPE0
                                                                                      CALL DBCONF (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
                		                                                               temp4=- factor_psi_mit*FV/PRECISION
                		                                                             
                		                                                            
                		                                                               XX(1)=X(1)
                		                                                               XX(2)=X(2)
                                                                                   
                                                                                  ! INITIAL VALUES AT CORNER SOLUTION ALPHA=1 
                                                                                    !X(1)
                                                                                     XGUESS(1)=MAX(SAVINGSFUNCTION(jj,iip,iit,tt)-(IPARA*1.0D0)*0.01D0,0.0D0)
                                                                                     !X(2)
                                                                                     XGUESS(2)=ALPHAFUNCTION(jj,iip,iit,tt)+(IPARA*1.0D0)*0.04D0
                                                                                    !CALL OPTIMIZATION SUBROUTINE
                                                                                      CALL BOUNDSTYPE0
                                                                                      CALL DBCONF (PROBLEM_tt, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FV)
!                		                                                               print *, "2", fv 
                		                                                               temp5=- factor_psi_mit*FV/PRECISION
!                                                                                  print *, savingfunction0(jj,iip,iit,tt), xx(1), x(1), xx(2), x(2)
!                                                                                  print *, VALUEFUNCTION0(jj,iip,iit,tt), temp4, temp5
!                                                                                  pause
                                                                                  If ( temp4 .ge. temp5) then
                                                                                      ALPHAFUNCTION(jj,iip,iit,tt)=XX(2)
                                                                                      SAVINGSFUNCTION(jj,iip,iit,tt)=XX(1)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=temp4
				   	        	   	                                              else
				   	        	   	                                                  ALPHAFUNCTION(jj,iip,iit,tt)=X(2)
                                                                                      SAVINGSFUNCTION(jj,iip,iit,tt)=X(1)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=temp5
				   	        	   	                                              endif
				   	        	   	                                            
				   	        	   	                                              if ( (IPARA .eq. 12) .AND. ( VALUEFUNCTION(jj,iip,iit,tt) .lt. VALUEFUNCTION1(jj,iip,iit,tt) ) ) then
				   	        	   	                                                  SAVINGSFUNCTION(jj,iip,iit,tt)=SAVINGSFUNCTION1(jj,iip,iit,tt)
				   	        	   	                                                  VALUEFUNCTION(jj,iip,iit,tt)=VALUEFUNCTION1(jj,iip,iit,tt)
				   	        	   	                                                  ALPHAFUNCTION(jj,iip,iit,tt)=ALPHAFUNCTION1(jj,iip,iit,tt)
				   	        	   	                                              else; endif
				   	        	   	                                      ipara=ipara+1
				   	        	   	                                    
				   	        	   	                                    ENDDO
!                                                               print *, "7.1e fin"
!         print *, "coh", cohi                                              
!         print *,   SAVINGSFUNCTION0(jj,iip,iit,tt), 0.0d0, VALUEFUNCTION0(jj,iip,iit,tt)   
!         print *,   SAVINGSFUNCTION(jj,iip,iit,tt), alphaFUNCTION(jj,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt)
!         print *,   SAVINGSFUNCTION1(jj,iip,iit,tt), ALPHAFUNCTION1(jj,iip,iit,tt), VALUEFUNCTION1(jj,iip,iit,tt)                                     
!         pause                                        
             
                 ENDDO !IIT
  
               
             ENDDO !IIP
!   
!                        WRITE (1001, 106), tt, jj, GRIDW(jj,tt),  VALUEFUNCTION(jj,1,:,tt), VALUEFUNCTION(jj,2,:,tt), VALUEFUNCTION(jj,3,:,tt), VALUEFUNCTION(jj,4,:,tt), VALUEFUNCTION(jj,5,:,tt),VALUEFUNCTION(jj,6,:,tt), VALUEFUNCTION(jj,7,:,tt)
!                  
!                         WRITE (1004, 106), tt, jj, GRIDW(jj,tt), ALPHAFUNCTION(jj,1,:,tt), ALPHAFUNCTION(jj,2,:,tt), ALPHAFUNCTION(jj,3,:,tt), ALPHAFUNCTION(jj,4,:,tt), ALPHAFUNCTION(jj,5,:,tt),ALPHAFUNCTION(jj,6,:,tt), ALPHAFUNCTION(jj,7,:,tt)
!                         WRITE (1007, 106), tt, jj, GRIDW(jj,tt), SAVINGSFUNCTION(jj,1,:,tt), SAVINGSFUNCTION(jj,2,:,tt), SAVINGSFUNCTION(jj,3,:,tt), SAVINGSFUNCTION(jj,4,:,tt), SAVINGSFUNCTION(jj,5,:,tt),SAVINGSFUNCTION(jj,6,:,tt), SAVINGSFUNCTION(jj,7,:,tt)
!                
!                       
!                         WRITE (1003, 106), tt, jj, GRIDW(jj,tt),  VALUEFUNCTION1(jj,1,:,tt), VALUEFUNCTION(jj,2,:,tt), VALUEFUNCTION1(jj,3,:,tt), VALUEFUNCTION1(jj,4,:,tt), VALUEFUNCTION1(jj,5,:,tt),VALUEFUNCTION1(jj,6,:,tt), VALUEFUNCTION1(jj,7,:,tt)
!                         WRITE (1006, 106), tt, jj, GRIDW(jj,tt),  ALPHAFUNCTION1(jj,1,:,tt), ALPHAFUNCTION1(jj,2,:,tt), ALPHAFUNCTION1(jj,3,:,tt), ALPHAFUNCTION1(jj,4,:,tt), ALPHAFUNCTION1(jj,5,:,tt),ALPHAFUNCTION1(jj,6,:,tt), ALPHAFUNCTION1(jj,7,:,tt)
!                         WRITE (1009, 106), tt, jj, GRIDW(jj,tt),  SAVINGSFUNCTION1(jj,1,:,tt), SAVINGSFUNCTION1(jj,2,:,tt), SAVINGSFUNCTION1(jj,3,:,tt), SAVINGSFUNCTION1(jj,4,:,tt), SAVINGSFUNCTION1(jj,5,:,tt),SAVINGSFUNCTION1(jj,6,:,tt), SAVINGSFUNCTION1(jj,7,:,tt)
!                          
!                         WRITE (1002, 106), tt, jj, GRIDW(jj,tt),  VALUEFUNCTION0(jj,1,:,tt), VALUEFUNCTION0(jj,2,:,tt), VALUEFUNCTION0(jj,3,:,tt), VALUEFUNCTION0(jj,4,:,tt), VALUEFUNCTION0(jj,5,:,tt),VALUEFUNCTION0(jj,6,:,tt), VALUEFUNCTION0(jj,7,:,tt)
!                         WRITE (1005, 106), tt, jj, GRIDW(jj,tt),  ALPHAFUNCTION0(jj,1,:,tt), ALPHAFUNCTION0(jj,2,:,tt), ALPHAFUNCTION0(jj,3,:,tt), ALPHAFUNCTION0(jj,4,:,tt), ALPHAFUNCTION0(jj,5,:,tt),ALPHAFUNCTION0(jj,6,:,tt), ALPHAFUNCTION0(jj,7,:,tt)
!                         WRITE (1008, 106), tt, jj, GRIDW(jj,tt),  SAVINGSFUNCTION0(jj,1,:,tt), SAVINGSFUNCTION0(jj,2,:,tt), SAVINGSFUNCTION0(jj,3,:,tt), SAVINGSFUNCTION0(jj,4,:,tt), SAVINGSFUNCTION0(jj,5,:,tt),SAVINGSFUNCTION0(jj,6,:,tt), SAVINGSFUNCTION0(jj,7,:,tt)
!!                      
                         
                  EVALUEFUNCTION(jj,tt,:)=0.00d0
                  
                  if ( young .eq. 1) then;
                          DO iip = 2, Miip-1
                               DO iit=1,Miit-1
                                EVALUEFUNCTION(jj,tt,1)=EVALUEFUNCTION(jj,tt,1)+VALUEFUNCTION(jj,iip,iit,tt)*PSHOCKS(iip,2,young)*TSHOCKS(iit,2,young)
                                EVALUEFUNCTION(jj,tt,2)=EVALUEFUNCTION(jj,tt,2)+VALUEFUNCTION0(jj,iip,iit,tt)*PSHOCKS(iip,2,young)*TSHOCKS(iit,2,young)
                                EVALUEFUNCTION(jj,tt,3)=EVALUEFUNCTION(jj,tt,3)+VALUEFUNCTION1(jj,iip,iit,tt)*PSHOCKS(iip,2,young)*TSHOCKS(iit,2,young)
                               ENDDO !iit
                          ENDDO !iip
                 else;
                                EVALUEFUNCTION(jj,tt,1)=VALUEFUNCTION(jj,1,1,tt)
                                EVALUEFUNCTION(jj,tt,2)=VALUEFUNCTION0(jj,1,1,tt)
                                EVALUEFUNCTION(jj,tt,3)=VALUEFUNCTION1(jj,1,1,tt)
                              
                endif; 
                 
               WRITE (1010, 106), tt, jj, GRIDW(jj,tt), EVALUEFUNCTION(jj,tt,:)
            
    jj=jj-1
  enddo !jj
     
    
!CHECKING RESULTS
if ( young .eq. 1) then
                        ! CHECK1: VALUE FUNCTION IF INCREASING OVER II FOR A GIVEN JJ. INCONSISTENCE TYPE A
                            DO jj=1,J
                            DO iip=2,GHP
                            DO iit=1,GHT
                                  IF (VALUEFUNCTION(jj,iip,iit,tt) .Lt. VALUEFUNCTION(jj,iip-1,iit,tt)) then
                                  WRITE (509,110), jj, iip, iit, tt, VALUEFUNCTION(jj,iip,iit,tt)/VALUEFUNCTION(jj,iip-1,iit,tt), VALUEFUNCTION(jj,iip,iit,tt), VALUEFUNCTION(jj,iip-1,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), ALPHAFUNCTION(jj,iip,iit,tt)
                                  ELSE; ENDIF
                            ENDDO; ENDDO; ENDDO

                        ! CHECK2: VALUE FUNCTION IF INCREASING OVER JJ FOR A GIVEN II INCONSISTENCE TYPE B
                            DO iip=1,GHP
                            DO iit=1,GHT
                            DO jj=2,J
                            IF (VALUEFUNCTION(jj,iip,iit,tt) .Lt. 1.000001D0*VALUEFUNCTION(jj-1,iip,iit,tt)) then
                            WRITE (510,110), jj, iip, iit, tt, VALUEFUNCTION(jj,iip,iit,tt)/VALUEFUNCTION(jj-1,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt), VALUEFUNCTION(jj-1,iip,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), ALPHAFUNCTION(jj,iip,iit,tt)
                            ELSE; ENDIF
                            ENDDO; ENDDO; ENDDO

                            
                            
                        ! CHECK3: VALUE FUNCTION OF UNRESTRICTED PROBLEM >= VALUE FUNCTION IMPOSING ALPHA=0 INCONSISTENCE TYPE 1
                            DO iip=1,GHP
                            DO iit=1,GHT
                            DO jj=1,J
                            IF (VALUEFUNCTION(jj,iip,iit,tt) .Lt. 1.00001D0*VALUEFUNCTION0(jj,iip,iit,tt)) then
                            WRITE (511,110), jj,iip,iit,tt, VALUEFUNCTION(jj,iip,iit,tt)/VALUEFUNCTION0(jj,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt), VALUEFUNCTION0(jj,iip,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), SAVINGSFUNCTION0(jj,iip,iit,tt)
                            ELSE; ENDIF
                            ENDDO; ENDDO; ENDDO
                          
                         !CHECK4: VALUE FUNCTION OF UNRESTRICTED PROBLEM >= VALUE FUNCTION FORCING THE AGENT TO PAY PARTICIPATION COST INCONSISTENCE TYPE 2
                            DO iip=1,GHP
                            DO iit=1,GHT
                            DO jj=1,J
                            IF (VALUEFUNCTION(jj,iip,iit,tt) .Lt. 1.00001D0*VALUEFUNCTION1(jj,iip,iit,tt)) then
                            WRITE (512,110), jj,iip,iit, tt, VALUEFUNCTION(jj,iip,iit,tt)/VALUEFUNCTION1(jj,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt), VALUEFUNCTION1(jj,iip,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), SAVINGSFUNCTION1(jj,iip,iit,tt), ALPHAFUNCTION(jj,iip,iit,tt), ALPHAFUNCTION1(jj,iip,iit,tt)
                            !PAUSE
                            ELSE; ENDIF
                            ENDDO; ENDDO; ENDDO;
                          
                          ! CHECK5: 
                          !   IF ALPHA=0 VF0>VF1
                           !  IF ALPHA>0 VF1>VFO INCONSISTENCE TYPE 3 
                            DO iip=1,GHP
                            DO iit=1,GHT
                            DO jj=1,J
                            IF ( (ALPHAFUNCTION(jj,iip,iit,tt) .gt. 0.0d0 ) .and. ( VALUEFUNCTION1(jj,iip,iit,tt) .Lt. 0.99999d0*VALUEFUNCTION0(jj,iip,iit,tt) ) ) then
                            WRITE (513,110), jj,iip,iit, tt, VALUEFUNCTION1(jj,iip,iit,tt)/VALUEFUNCTION0(jj,iip,iit,tt), ALPHAFUNCTION(jj,iip,iit,tt), VALUEFUNCTION0(jj,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt),VALUEFUNCTION1(jj,iip,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), SAVINGSFUNCTION1(jj,iip,iit,tt)  
                            ELSE; ENDIF
                            IF ( (ALPHAFUNCTION(jj,iip,iit,tt) .eq. 0.0000000d0) .and. (VALUEFUNCTION0(jj,iip,iit,tt) .Lt. 0.99999d0*VALUEFUNCTION1(jj,iip,iit,tt)) )then
                            WRITE (513,110), jj,iip,iit, tt, VALUEFUNCTION0(jj,iip,iit,tt)/VALUEFUNCTION1(jj,iip,iit,tt), ALPHAFUNCTION(jj,iip,iit,tt), VALUEFUNCTION0(jj,iip,iit,tt), VALUEFUNCTION(jj,iip,iit,tt),VALUEFUNCTION1(jj,iip,iit,tt), SAVINGSFUNCTION(jj,iip,iit,tt), SAVINGSFUNCTION1(jj,iip,iit,tt)   
                            ELSE; ENDIF
                            ENDDO; ENDDO; ENDDO;
  else; endif; !CHECKS for youngs
    
 JJ=0
 II=0  
  tt=tt-1
  ENDDO !tt
  
 !9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.9.
 ! Calcule pseudo log-likelihood 

  pseudo2_ll(:,gg,ccrra)=0.d0

  pseudo1_ll(:,gg,ccrra)=0.d0

      
      
!NUMBER OF OBSERVATION WHOLE SAMPLE AND THOSE WITHIN KNOWN Lagged alpha
 OBS1(:)=0
 OBS2(:)=0
do hh=1,SIZE_E_GROUP
      Z(1)=1.d0
      Z(2)=IDATAB(4,hh)*1.d0  !SOUTH=1
      Z(3)=1.d0-Z(2) !SOUTH=0
      Z(4)=IDATAB(5,hh)*1.d0  !Blue-collar=1
      Z(5)=1.d0-Z(4) !Blue-collar=1
      Z(6)=IDATAB(6,hh)*1.d0  !Young=1
      Z(7)=1.d0-Z(6)   !Young=0
      Z(8)=IDATAB(7,hh)*1.d0  !NO house owner=1
      Z(9)=1.d0-Z(8)    !NO house owner=0
      Z(10)=1-IDATAB(8,hh)*1.d0  !Below median wit
      Z(11)=1.d0-Z(10)   !Above median wit
      Z(12)=1-IDATAB(9,hh)*1.d0  !Below median r1it
      Z(13)=1.d0-Z(12)   !Above median r1it
      
                               do vv=1,13
                              OBS1(vv)=OBS1(vv)+Z(vv)
                              enddo
                              !Those with L.alpha>0 / L.alpha==0  
                              
                              if ( (IDATAA(5,hh) .gt. 0d0 ) .OR. ( IDATAA(5,hh) .eq. 0d0 ) ) then; 
                              if (IDATAA(5,hh) .gt. 0d0 ) then; Z(14)=1.d0; Z(15)=0.d0; else; Z(14)=0.d0; Z(15)=1.d0; endif;    
                                do vv=1,15
                                OBS2(vv)=OBS2(vv)+Z(vv)
                                enddo
                              else; endif 

enddo 

!PSEUDO LIKELIHOOD
  do hh=1,SIZE_E_GROUP
                           
      wit=IDATAA(1,hh)
      ratio1=IDATAA(2,hh)
      ratio2=IDATAA(3,hh)
      AGE=IDATAB(3,hh)
      

      
      Alphai=IDATAA(4,hh)
      LagAlphai=IDATAA(5,hh)
      Z(1)=1.d0
      Z(2)=IDATAB(4,hh)*1.d0  !SOUTH=1
      Z(3)=1.d0-Z(2) !SOUTH=0
      Z(4)=IDATAB(5,hh)*1.d0  !Blue-collar=1
      Z(5)=1.d0-Z(4) !Blue-collar=1
      Z(6)=IDATAB(6,hh)*1.d0  !Young=1
      Z(7)=1.d0-Z(6)   !Young=0
      Z(8)=IDATAB(7,hh)*1.d0  !NO house owner=1
      Z(9)=1.d0-Z(8)    !NO house owner=0
      Z(10)=1-IDATAB(8,hh)*1.d0  !Below median wit
      Z(11)=1.d0-Z(10)   !Above median wit
      Z(12)=1-IDATAB(9,hh)*1.d0  !Below median r1it
      Z(13)=1.d0-Z(12)   !Above median r1it
      
       if ( (IDATAA(5,hh) .gt. 0d0 ) .OR. ( IDATAA(5,hh) .eq. 0d0 ) ) then; 
                              if (IDATAA(5,hh) .gt. 0d0 ) then; Z(14)=1.d0; Z(15)=0.d0; else; Z(14)=0.d0; Z(15)=1.d0; endif;    
       else; endif 
    
!      print *, age
!      print *, OBS_wit,ratio1,ratio2
        call  INTERPOLALL(AGE,wit,ratio1,ratio2,lli0)
     
      DEv1_v0=lli0/lant*dexp((1.d0-gamma0)*dlog(ratio1))
      if (Alphai .gt. 0d0 ) then; lli=dexp(Dev1_v0)/(1.d0+dexp(Dev1_v0)); else; lli=1.0d0/(1.d0+dexp(Dev1_v0)); endif
      
!       print *, lli,  DEv1_v0, IDATAA(5,hh) 
      
      DO vv=1,13 
     
        pseudo1_ll(vv,gg,ccrra)=pseudo1_ll(vv,gg,ccrra)+dlog(lli)*Z(vv)
          
      ENDDO !vv 
                 
      !Those with L.alpha>0 / L.alpha==0 
      if ( (IDATAA(5,hh) .gt. 0d0 ) .OR. ( IDATAA(5,hh) .eq. 0d0 ) ) then;   
          DO vv=1,15
             pseudo2_ll(vv,gg,ccrra)=pseudo2_ll(vv,gg,ccrra)+dlog(lli)*Z(vv)
          ENDDO !vv 
      else; endif            


! WRITE FILES WITH INDIVIDUAL CONTRIBUTION TO THE PLL, FOR DIF. SPEC. FOR PARTICIPATION COST 
            if ( cccc .eq. 1) then;    write (4313131, 1332), ccrra, gg, IDATAB(:,hh), IDATAA(4:5,hh), gamma0,  g_coef(gg,1),  lli; else; endif
            if ( cccc .eq. 2) then;    write (4313132, 1332), ccrra, gg, IDATAB(:,hh), IDATAA(4:5,hh), gamma0,  g_coef(gg,1), lli; else; endif
            if ( cccc .eq. 3) then;    write (4313133, 1332), ccrra, gg, IDATAB(:,hh), IDATAA(4:5,hh), gamma0,  g_coef(gg,1), lli; else; endif
            if ( cccc .eq. 4) then;    write (4313134, 1332), ccrra, gg, IDATAB(:,hh), IDATAA(4:5,hh), gamma0,  g_coef(gg,1), lli; else; endif

enddo  !hh
 
! WRITE FILES WITH PLL FOR DIF. SPEC. FOR PARTICIPATION COST 

if ( cccc .eq. 1) then; write (1313131, 104),ccrra, gg, gamma0, g_coef(gg,1), pseudo1_ll(:,gg,ccrra); else; endif  
if ( cccc .eq. 2) then; write (1313132, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo1_ll(:,gg,ccrra); else; endif
if ( cccc .eq. 3) then; write (1313133, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo1_ll(:,gg,ccrra);  else; endif
if ( cccc .eq. 4) then; write (1313134, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo1_ll(:,gg,ccrra); else; endif

if ( cccc .eq. 1) then; write (2313131, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo2_ll(:,gg,ccrra); else; endif  
if ( cccc .eq. 2) then; write (2313132, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo2_ll(:,gg,ccrra); else; endif
if ( cccc .eq. 3) then; write (2313133, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo2_ll(:,gg,ccrra);  else; endif
if ( cccc .eq. 4) then; write (2313134, 104), ccrra, gg, gamma0, g_coef(gg,1), pseudo2_ll(:,gg,ccrra); else; endif



      
   enddo !llll  
  enddo !GG
 enddo !cccc
enddo ! CCRRA  


CONTAINS
!10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//10.1//
SUBROUTINE PROBLEM_tt (N, X, F)
  ! SPECIFICATIONS FOR ARGUMENTS
      INTEGER N
        INTEGER cc, iii
		   REAL*8 X(*), F(*), Z(3), D(2), ZF(3), Z2(3)
	 	     REAL*8 W, cmax, COH, vc(3), VFmg, x2, nWPrime, uc(3), v2, vPRIME, EVPRIME(LDR+1,2), D_vPRIME, dEVPRIME(2), dev
	 	      REAL*8 CONSUM_tt, umg_Ctt, Eumg_Cttprime, Eumg_alpha, Consum_T(ghp,GHT,LDR+1), u_ctt, Eu_Cttprime
	 	      REAL*8 AUX1_tt, AUX2_tt, AUX3_tt, AUX4_tt, AUX1_T_1, AUX2_T_1, AUX3_T_1, AUX4_T_1, Eumg_CT,  Ev_alpha1, Ev_alpha0, tempf1(2), tempx1(2), tempx2(2) 
                
102 FORMAT ( <1> i6, <4> f12.5)

105 FORMAT ( <2> i6, <15> f15.10)
     
                 COH=COHi   !Cash on hand at tt
                    tempx1(1)=tempx1(2)
                    tempx2(1)=tempx2(2)
                       if ( ( X(1) .ne. X(1) ) .or. (x(2).ne. x(2) ) )then; print *, "problem_tt", tt,jj,iip,iit; x(1)=max(0.1d0,min(tempx1(1)-0.01d0,0.6d0)); x(2)=min(tempx2(1)+0.01d0,0.50d0); endif
        	                
!			                       !UTILITY OF CONSUMPTION AT tt
!			                       
						           AUX1_tt=COH*(1.00000000D0-X(1))
                                   u_Ctt=feli(AUX1_tt,Gamma0)
!						         
						          !PERIOD tt+1  EXPECTED VALUE FUNCTION              
						             if ( X(2) .gt. 0.000000000d0) then; x2=1.00d0; else; x2=0.0000d0; endif;
                                  			 Eu_Cttprime=0.000d0
                                  	         AUX2_tt=COH*X(1)-G_(1,1)*x2
                                  	         
                                  	         if ( AUX2_tt .lt. 0.0d0) then
                                  	         AUX1_tt=COH*(1.00000000D0-X(1))+AUX2_tt
                                             u_Ctt=feli(AUX1_tt,Gamma0)
                                  	         AUX2_tt=0.d0
                                  	         else; endif
                                  	         
                                  			  DO rar=2,LDR+1
						                         !WPRIME=dlog(dexp(COH*X(1)*R(1))*dexp(COH*X(1)*X(2)*R(cc))/dexp(COH*X(1)*X(2)*R(1))/(dexp(G_(1,1)*R(1))*dexp(G_(1,1)*X(2)*R(cc))/dexp(G_(1,1)*X(2)*R(1))))
						                         AUX3_tt=R(1)+X(2)*(R(rar)-R(1)); 
						                         nWPRIME=AUX3_tt*AUX2_tt;
						                         CALL INTERPOLAV(tt,nWPRIME,vPRIME)
                                                 EVPRIME(rar,2)=vPRIME
                                                 Eu_Cttprime= Eu_Cttprime+EVPRIME(rar,2)*PrR(rar)
!                                                 print *, "EV", EVPRIME(rar,2),EVPRIME(rar,1), (EVPRIME(rar,2)-EVPRIME(rar,1))/(coh*0.002d0)
!                                                 pause
                                             ENDDO  
                                  			 
                                  			
!!        	                      
        	                      
!        	                    
!        	                   if ( ( ( tt .eq. 27) .and. (jj .eq. 36) ) ) then
!                               print *, tt, iip, iit, coh, x(1), x(2)
!                               print *, F(1),  u_Ctt, Eu_Cttprime
!                               print *, AUX1_tt, AUX2_tt
!                               pause
!                               else; endif
                                  
        	                     !value function
        	                      F(1)=-PRECISION*(u_Ctt+Beta0*Eu_Cttprime)
!!        	                      if (F(1) .lt. 0.01d0) then; F(1)=1/sqrt(epsilon(1.d0)); else; endif
        	                      tempx1(2)=x(1)
        	                      tempx2(2)=x(2) 
!        	                   if ( ( ( tt .eq. 28) .and. (jj .eq. 77) ) ) then
!                               print *, tt, iip, iit, coh, x(1), x(2)
!!                               pause
!                               else; endif
                               
                                                                     
END SUBROUTINE PROBLEM_tt

!10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\10.2\\

!SUBROUTINE TO FIND X(1) WHEN X(2) EQUALS 0 at T-1
SUBROUTINE PROBLEM_tt_b (N2, X, F)
 ! SPECIFICATIONS FOR ARGUMENTS
      INTEGER M, N2
        INTEGER cc, iii
		   REAL*8 X(*), F(*), Z(3), D(2), ZF(3)
	 	     REAL*8 W, cmax, var1, var2, COH, uc(3), vc(2), WPRIME, nWPRIME,v2, Eumg_Cttprime, umg_Ctt, Eu_Cttprime, u_Ctt
	 	      REAL*8 CONSUM_tt_b,   EVPRIME(1,2), D_vPRIME, AUX1_tt, AUX2_tt , AUX3_tt , tempx1(2)
                
102 FORMAT ( <1> i6, <4> f12.5)
     
                 COH=COHi
                 		   !UTILITY OF CONSUMPTION AT tt
!			                       
						           AUX1_tt=COH*(1.00000000D0-X(1))
                                   u_Ctt=feli(AUX1_tt,Gamma0)
!						         
						          !PERIOD tt+1  EXPECTED VALUE FUNCTION              
						                         nWPRIME=COH*X(1)*R(1)
						               
                                                 CALL INTERPOLAV(tt,nWPRIME,vPRIME)
                                                 Eu_Cttprime= VPRIME
!                                                
!!        	                      !value function
        	                     
        	                       if ( ( X(1) .ne. X(1) ) .or. (x(2).ne. x(2) ) )then; print *, "problem_ttc", tt,jj,iip,iit; else; endif
                                  F(1)=-PRECISION*(u_Ctt+Beta0*Eu_Cttprime)
                                     
END SUBROUTINE PROBLEM_tt_b 



!10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\10.3\\

!SUBROUTINE TO FIND X(1) WHEN X(2) EQUALS 1 at T-1
SUBROUTINE PROBLEM_tt_c (N, X, F)
  !SPECIFICATIONS FOR ARGUMENTS
      INTEGER M, N
        INTEGER cc, iii
		   REAL*8 X(*), F(*), Z(3), D(2), ZF(3), Z2(3)
	 	     REAL*8 W, cmax, COH, vc(3),  VFmg, x2, WPrime, nWPrime, uc(3), v2, D_vPRIME, EVPRIME(LDR+1,2)
	 	      REAL*8 CONSUM_T_1, CONSUM_T(ISI,3), u_Ctt, Eu_Cttprime, Eumg_alpha
	 	      REAL*8 AUX1_tt, AUX2_tt, AUX3_tt, AUX4_tt  
                
102 FORMAT ( <1> i6, <4> f12.5)

105 FORMAT ( <2> i6, <15> f15.10)
     
                 COH=COHi   !Cash on hand at tt
                 !UTILITY OF CONSUMPTION AT tt
!			                       
						           AUX1_tt=COH*(1.00000000D0-X(1))
                                   u_Ctt=feli(AUX1_tt,Gamma0)
!						         
						          !PERIOD tt+1  EXPECTED VALUE FUNCTION              
						                      Eu_Cttprime=0.000d0
						                      AUX2_tt=COH*X(1)-G_(1,1)
						                     
						                     if ( AUX2_tt .lt. 0.0d0) then
                                  	         AUX1_tt=max(0.0001d0,COH*(1.00000000D0-X(1))+AUX2_tt)
                                             u_Ctt=feli(AUX1_tt,Gamma0)
                                  	         AUX2_tt=0.d0
                                  	         else; endif
						                      
                                  			  DO rar=2,LDR+1
						                         !WPRIME=dlog(dexp(COH*X(1)*R(1))*dexp(COH*X(1)*X(2)*R(cc))/dexp(COH*X(1)*X(2)*R(1))/(dexp(G_(1,1)*R(1))*dexp(G_(1,1)*X(2)*R(cc))/dexp(G_(1,1)*X(2)*R(1))))
						                         AUX1_tt=R(1)+X(2)*(R(rar)-R(1)); 
						                         nWPRIME=AUX1_tt*AUX2_tt;
						               
                                                 CALL INTERPOLAV(tt,nWPRIME,vPRIME)
                                                 EVPRIME(rar,2)=vPRIME
                                                 Eu_Cttprime= Eu_Cttprime+EVPRIME(rar,2)*PrR(rar)
!                                                 print *, "EV", EVPRIME(rar,2),EVPRIME(rar,1), (EVPRIME(rar,2)-EVPRIME(rar,1))/(coh*0.002d0)
!                                                 pause
                                             ENDDO  
                                  			 
                                  			
!!        	                      
        	                      !D(1)- value function
        	                   F(1)=-PRECISION*(u_Ctt+Beta0*Eu_Cttprime) 
        	                      if ( ( X(1) .ne. X(1) ) .or. (x(2).ne. x(2) ) )then; print *, "problem_ttc", x(1), x(2), tt,jj,iip,iit; else; endif
        	                      
!        	                   if ( ( f(1) .lt. -1.d0/epsilon(1.d0)) .or. (( tt .eq. 26) .and. (jj .eq. 19)) ) then
!                               print *, tt, iip, iit, coh, x(1), x(2), -f(1)/precision
!                               pause
!                               else; endif
                               
!                                if ( iparam(3) .gt. 100) then
!                                print *, tt, jj, iip, iit, x(1), x(2), F(1), iparam(3)
!                                pause
!                                else; endif
						                                       
END SUBROUTINE PROBLEM_tt_c 

SUBROUTINE VALUEF(ttt,coh,gv,Xs,Ve)
INTEGER :: ttt
REAL*8 :: COH, gv, Xs(N), x2, EV_P, W_P, I_VP, uctt0, Ve

!COMPUTING Value function AT PERIOD tt
                                                            if ( Xs(2) .gt. 0.0000100D0) then; x2=1.000000000D0; else; x2=0d0; endif
        				                                    if ( Xs(2) .gt. 1.0000000D0) then; Xs(2)= 1.0000000d0; else; endif
				                                            !FELICITY FUNCTION AT T-1
				                                             uCtt0=feli(COH*(1.d0-Xs(1)), gamma0)
                    				                        
            !				                                !EXPECTED VALUE FUNCTION(W_PRIME)
            !				                                      !computed by interpolating the expected value function
				                                                   EV_P=0.000D0
                                                                    DO rar=2,LDR+1
						                                             W_P=(COH*Xs(1)-gv*x2)*(R(1)+ x2*Xs(2)*(R(rar)-R(1)))
                                                                     CALL INTERPOLAV(tt,W_P,I_VP)
                                                                     EV_P=EV_P+I_VP*PrR(rar)
            !                                                         print *, cohi, x(1), x(2) 
            !                                                         print *, w_prime,I_Vprime, EV_prime
                                                                    ENDDO  
                                                                    Ve=uctt0+Beta0*EV_P !        				                     
!                                                               print *, "coh", coh, xs(1),xs(2)
!                                                               print *, uctt0, EV_p, ve 
				                                            !VALUE FUNCTION AT T-1,jj,ii 
END SUBROUTINE

!12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1//12.1
!Interpolation of expected value function
SUBROUTINE INTERPOLAV(ttt,nWPRIME0,VPRIME0)
  INTEGER :: ttt, INDJ, iiii 
   REAL*8 :: nWPRIME0, UW, LW
   REAL*8 :: VPRIME0
     
                    
                     if ( nWPRIME0 .lt. 0.0000000100D0) then
                     
                                      !  Expected value of Value Function at W(t+1)=0 (no savings at t)
                                      vPRIME0=EVALUEFUNCTION(1,ttt+1,1)
                     else
                          !FINDING THE CORRESPONDING CELL FOR w(tt+1) AT THE GRIDW OF tt+1
                          INDJ=1
                          DO WHILE ( nWPRIME0 .gT. GRIDW(INDJ,ttt+1) .and. (INDJ .lt. J) )
                             INDJ=INDJ+1
                          ENDDO
!                         print *, "indj", indj
                    	                  !  Expected value of Value Function at W(t+1) if 0<W(t+1)<GRIDW(J,tt)
	                                      IF ( ( INDJ .GT. 1)  .AND. (INDJ .Le. J)  )THEN
                                                            LW=GRIDW(INDJ-1,ttt+1)
                                                            UW=GRIDW(INDJ,ttt+1)
                                                            temp1=EVALUEFUNCTION(INDJ-1,ttt+1,1)
                                                            temp2=EVALUEFUNCTION(INDJ,ttt+1,1)
                                                            vPRIME0=temp1+(nWPRIME0-LW)/(UW-LW)*(temp2-temp1)                                                                                        
                                            ENDIF
                                        
                                             !  Expected value of Value Function at W(t+1) if W(t+1)>GRIDW(J,tt) (extrapolation)
                                                               ! linear extrapolation corrected by the factor exp(-(W(t+1)-UW)/((W(t+1)+UW)/2))
                                             IF ( nWPRIME0 .gT. GRIDW(J,ttt+1) )THEN
                                                                LW=GRIDW(J-1,ttt+1)
                                                                UW=GRIDW(J,ttt+1)
                                                                temp1=EVALUEFUNCTION(J-1,ttt+1,1)
                                                                temp2=EVALUEFUNCTION(J,ttt+1,1)
                                                                vPRIME0=temp2+(nWPRIME0-UW)/(UW-LW)*(temp2-temp1)*dexp(-4*(nWPRIME0-UW)/((nWPRIME0+UW)))
                                                                
!                                                                if ( tt .eq. T-2) THEN
!                                                                print *, INDJ, nWPRIME0, UW, vPRIME0, temp6
!                                                                pause
!                                                                ELSE; ENDIF
                                              endif
                                         
                                                    
           endif 
! print *, nwprime0, indj, vprime0      
ENDSUBROUTINE INTERPOLAV



SUBROUTINE INTERPOLALL(ttt,nWPRIME0,ratio10, ratio20,LL0)
  INTEGER :: ttt, INDJ, INDJJ,INDP, INDT, iiii, jjjj 
   REAL*8 :: nWPRIME0,  ratio10, ratio20, Ev1, Ev0
   REAL*8 :: UW, LW, UPI, LPI, UTI, LTI
   REAL*8 :: AUXL_0(2), AUXL_1(2), AUXL_2(2), AUXL_3, AUXL_4, V0jp(2), V1jp(2), v0j(2), v1j(2), v0,v1 
   REAL*8 :: LL0
     
      INDP=1
      DO WHILE ( ratio10 .gT. PSHOCKS(INDP,1,young) .and. (INDP .lt. GHP) )
      INDP=INDP+1
      ENDDO
      if ( ratio10 .gT. PSHOCKS(GHP,1,young)) then; INDP=GHP+1; else; endif
      INDT=1
      DO WHILE ( ratio20 .gT. TSHOCKS(INDT,1,young) .and. (INDT .lt. GHT) )
      INDT=INDT+1
      ENDDO
      
      if ( ratio20 .gT. TSHOCKS(GHT,1,young)) then; INDT=GHT+1; else; endif


    ! Compute lli for W=0, iip, iit
                     if ( nWPRIME0 .lt. 0.0000000100D0) then
                                                             
                                                       ! ratio1 < min(grid P INCOME)
                                                         IF ( ( INDP .eq. 1) ) then   
                                                                                  
                                                                                  ! ratio2 < min(grid T INCOME)
                                                  
                                                                                  IF ( INDT .EQ. 1) THEN; 
                                                                                      LTI=TSHOCKS(1,1,young)
                                                                                      UTI=TSHOCKS(2,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,1,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,1,2,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,1,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,1,2,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,2,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,2,2,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,2,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,2,2,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                 
                                                                                  ELSE; ENDIF !INDT .EQ. 1
                                                                                  
                                                                                  
                                                       
                                                                               ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                                                       
                                                                                  IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                                                    LTI=TSHOCKS(INDT-1,1,young)
                                                                                     UTI=TSHOCKS(INDT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,1,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,1,INDT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,1,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,1,INDT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,2,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,2,INDT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,2,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,2,INDT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                  ELSE; ENDIF ! 1<INDT<=GHT  
                                                                                  
                                                                                  
                                                                               !  ratio2 > max(grid T INCOME)
                                                                                  IF ( INDT .EQ. GHT+1) THEN; 
                                                                                      LTI=TSHOCKS(GHT-1,1,young)
                                                                                      UTI=TSHOCKS(GHT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,1,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,1,GHT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,1,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,1,GHT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,2,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,2,GHT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,2,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,2,GHT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                  
                                                                                  ELSE; ENDIF !INDT .EQ. GHT+1
                                                              
                                                                 LPI=PSHOCKS(INDP-1,1,young)
                                                                 UPI=PSHOCKS(INDP,1,young)
                                                                                   V0j(1)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                                                   V1j(1)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                                                LL0=(V1j(1)-V0j(1))/abs(min(V1j(1),V0j(1)))
                                                                !LL0=(V1j(1)-V0j(1))
                                                         else  
                                                       endif ! INDP=1
                                                                      
                                                       
                                                       
                                                       ! min(grid P INCOME) < ratio1 < max(grid P INCOME) 
                                                       
                                                         IF ( ( INDP .gt. 1) .and. ( INDP .le. GHP) ) then   
                                                                                  
                                                                                  ! ratio2 < min(grid T INCOME)
                                                       
                                                                                  IF ( INDT .EQ. 1) THEN; 
                                                                                      LTI=TSHOCKS(1,1,young)
                                                                                      UTI=TSHOCKS(2,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP-1,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP-1,2,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP-1,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP-1,2,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP,2,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP,2,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                  
                                                                                  ELSE; ENDIF !INDT .EQ. 1
                                                                                  
                                                                                  
                                                       
                                                                               ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                                                       
                                                                                  IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                                                    LTI=TSHOCKS(INDT-1,1,young)
                                                                                     UTI=TSHOCKS(INDT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP-1,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP-1,INDT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP-1,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP-1,INDT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP,INDT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP,INDT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                  ELSE; ENDIF ! 1<INDT<=GHT  
                                                                                  
                                                                                  
                                                                               !  ratio2 > max(grid T INCOME)
                                                                                  IF ( INDT .EQ. GHT+1) THEN; 
                                                                                      LTI=TSHOCKS(GHT-1,1,young)
                                                                                      UTI=TSHOCKS(GHT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP-1,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP-1,GHT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP-1,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP-1,GHT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,INDP,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,INDP,GHT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,INDP,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,INDP,GHT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                  
                                                                                  ELSE; ENDIF !INDT .EQ. GHT+1
                                                              
                                                                 LPI=PSHOCKS(INDP-1,1,young)
                                                                 UPI=PSHOCKS(INDP,1,young)
                                                                                   V0j(1)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                                                   V1j(1)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                                                !LL0=V1j(1)-V0j(1)
                                                                LL0=(V1j(1)-V0j(1))/abs(min(V1j(1),V0j(1)))
                                                               
                                                         else  
                                                       endif ! 1<INDP<=GHP
                                                 
                                                  IF ( ( INDP .EQ. GHP+1) ) then   
                                                                                  
                                                                                  ! ratio2 < min(grid T INCOME)
                                                       
                                                                                  IF ( INDT .EQ. 1) THEN; 
                                                                                      LTI=TSHOCKS(1,1,young)
                                                                                      UTI=TSHOCKS(2,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP-1,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP-1,2,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP-1,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP-1,2,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP,1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP,2,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP,1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP,2,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                  
                                                                                  ELSE; ENDIF !INDT .EQ. 1
                                                                                  
                                                                                  
                                                       
                                                                               ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                                                       
                                                                                  IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                                                    LTI=TSHOCKS(INDT-1,1,young)
                                                                                     UTI=TSHOCKS(INDT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP-1,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP-1,INDT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP-1,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP-1,INDT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP,INDT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP,INDT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP,INDT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP,INDT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                  ELSE; ENDIF ! 1<INDT<=GHT  
                                                                                  
                                                                                  
                                                                               !  ratio2 > max(grid T INCOME)
                                                                                  IF ( INDT .EQ. GHT+1) THEN; 
                                                                                      LTI=TSHOCKS(GHT-1,1,young)
                                                                                      UTI=TSHOCKS(GHT,1,young)
                                                                                       
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP-1,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP-1,GHT,ttt)
                                                                                        V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP-1,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP-1,GHT,ttt)
                                                                                        V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                        
                                                                                        AUXL_0(1)=VALUEFUNCTION0(1,GHP,GHT-1,ttt)
                                                                                        AUXL_0(2)=VALUEFUNCTION0(1,GHP,GHT,ttt)
                                                                                        V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                                                        
                                                                                        AUXL_1(1)=VALUEFUNCTION1(1,GHP,GHT-1,ttt)
                                                                                        AUXL_1(2)=VALUEFUNCTION1(1,GHP,GHT,ttt)
                                                                                        V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                                                  
                                                                                  ELSE; ENDIF !INDT .EQ. GHT+1
                                                              
                                                                 LPI=PSHOCKS(GHP-1,1,young)
                                                                 UPI=PSHOCKS(GHP,1,young)
                                                                                   V0j(1)=V0jp(1)+(ratio10-UPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                                                   V1j(1)=V1jp(1)+(ratio10-UPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                                                !LL0=V1j(1)-V0j(1)
                                                                LL0=(V1j(1)-V0j(1))/abs(min(V1j(1),V0j(1)))
                                                         else  
                                                       endif ! INDP==GHP+1
                                                 
         else ! w greater than 0
                          !FINDING THE CORRESPONDING CELL FOR w(tt+1) AT THE GRIDW OF tt+1
          INDJ=1
          DO WHILE ( nWPRIME0 .gT. GRIDW(INDJ,ttt) .and. (INDJ .lt. J) )
             INDJ=INDJ+1
          ENDDO
          
         
!                         print *, "indj", indj
                	                  !  Expected value of Value Function at W(t+1) if 0<W(t+1)<GRIDW(J,tt)
	 IF ( ( INDJ .GT. 1)  .AND. (INDJ .Le. J)  )THEN
                    
                    DO jjjj=1,2
                       if (jjjj .eq. 1) then; indjj=INDJ-1; else; indjj=INDJ; endif                                                
                           ! ratio1 < min(grid P INCOME)
                             IF ( ( INDP .eq. 1) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                           
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(1,1,young)
                                     UPI=PSHOCKS(2,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                  else  
                           endif ! INDP=1
                                          
                           
                           
                           ! min(grid P INCOME) < ratio1 < max(grid P INCOME) 
                           
                             IF ( ( INDP .gt. 1) .and. ( INDP .le. GHP) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(INDP-1,1,young)
                                     UPI=PSHOCKS(INDP,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                     else  
                           endif ! 1<INDP<=GHP
                     
                      IF ( ( INDP .EQ. GHP+1) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(GHP-1,1,young)
                                     UPI=PSHOCKS(GHP,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-UPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-UPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                   else  
                            endif ! INDP==GHP+1
                    
!                    print *, "v0jp", v0jp(1), v0jp(2)
!                    print *, "v1jp", v1jp(1), v1jp(2)
!                    print *, "v0j", v0j(1), v0j(2)
!                    print *, "v1j", v1j(1), v1j(2)
!                            
                    enddo  !jjjj=1,2
                    LW=GRIDW(INDJ-1,ttt)
                    UW=GRIDW(INDJ,ttt)      
                    V0=V0j(1)+(nWprime0-LW)/(UW-LW)*(V0j(2)-V0j(1))
                    V1=V1j(1)+(nWprime0-LW)/(UW-LW)*(V1j(2)-V1j(1))
                    
!                    print *, "v", v0, v1
!                    pause
                    !LL0=V1-V0 
                    LL0=(V1-V0)/abs(min(V1,V0))
                    !LL0=-(V1-V0)/V0                                                                                          
      ENDIF  ! 1< indj<= J
                                        
                                             !  Expected value of Value Function at W(t+1) if W(t+1)>GRIDW(J,tt) (extrapolation)
                                                               ! linear extrapolation corrected by the factor exp(-(W(t+1)-UW)/((W(t+1)+UW)/2))
                                             IF ( nWPRIME0 .gT. GRIDW(J,ttt) )THEN
                                                                 DO jjjj=1,2
                                              if (jjjj .eq. 1) then; indjj=J-1; else; indjj=J; endif                                                
                           ! ratio1 < min(grid P INCOME)
                             IF ( ( INDP .eq. 1) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                           
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,2,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,2,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,2,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,2,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(1,1,young)
                                     UPI=PSHOCKS(2,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                  else  
                           endif ! INDP=1
                                          
                           
                           
                           ! min(grid P INCOME) < ratio1 < max(grid P INCOME) 
                           
                             IF ( ( INDP .gt. 1) .and. ( INDP .le. GHP) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP-1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP-1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP-1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP-1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,INDP,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,INDP,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,INDP,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,INDP,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(INDP-1,1,young)
                                     UPI=PSHOCKS(INDP,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-LPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-LPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                     else  
                           endif ! 1<INDP<=GHP
                     
                      IF ( ( INDP .EQ. GHP+1) ) then   
                                                      
                                                      ! ratio2 < min(grid T INCOME)
                           
                                                      IF ( INDT .EQ. 1) THEN; 
                                                          LTI=TSHOCKS(1,1,young)
                                                          UTI=TSHOCKS(2,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,2,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,2,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,2,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,2,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. 1
                                                      
                                                      
                           
                                                   ! min(grid T INCOME) < ratio2 < max(grid T INCOME) 
                           
                                                      IF ( ( INDT .gt. 1) .and. ( INDT .le. GHT) ) then  
                                                        LTI=TSHOCKS(INDT-1,1,young)
                                                         UTI=TSHOCKS(INDT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,INDT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,INDT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,INDT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,INDT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,INDT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,INDT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-LTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                      ELSE; ENDIF ! 1<INDT<=GHT  
                                                      
                                                      
                                                   !  ratio2 > max(grid T INCOME)
                                                      IF ( INDT .EQ. GHT+1) THEN; 
                                                          LTI=TSHOCKS(GHT-1,1,young)
                                                          UTI=TSHOCKS(GHT,1,young)
                                                           
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP-1,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP-1,GHT,ttt)
                                                            V0jp(1)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP-1,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP-1,GHT,ttt)
                                                            V1jp(1)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                            
                                                            AUXL_0(1)=VALUEFUNCTION0(indjj,GHP,GHT-1,ttt)
                                                            AUXL_0(2)=VALUEFUNCTION0(indjj,GHP,GHT,ttt)
                                                            V0jp(2)=AUXL_0(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_0(2)-AUXL_0(1))
                                                            
                                                            AUXL_1(1)=VALUEFUNCTION1(indjj,GHP,GHT-1,ttt)
                                                            AUXL_1(2)=VALUEFUNCTION1(indjj,GHP,GHT,ttt)
                                                            V1jp(2)=AUXL_1(1)+(ratio20-UTI)/(UTI-LTI)*(AUXL_1(2)-AUXL_1(1))
                                                      
                                                      ELSE; ENDIF !INDT .EQ. GHT+1
                                  
                                     LPI=PSHOCKS(GHP-1,1,young)
                                     UPI=PSHOCKS(GHP,1,young)
                                                       V0j(jjjj)=V0jp(1)+(ratio10-UPI)/(UPI-LPI)*(V0jp(2)-V0jp(1))
                                                       V1j(jjjj)=V1jp(1)+(ratio10-UPI)/(UPI-LPI)*(V1jp(2)-V1jp(1))
                                   else  
                            endif ! INDP==GHP+1
                    
!                    print *, "v0jp", v0jp(1), v0jp(2)
!                    print *, "v1jp", v1jp(1), v1jp(2)
!                    print *, "v0j", v0j(1), v0j(2)
!                    print *, "v1j", v1j(1), v1j(2)
!                            
                    enddo  !jjjj=1,2
                    LW=GRIDW(J-1,ttt)
                    UW=GRIDW(J,ttt)      
                    V0=V0j(1)+(nWprime0-UW)/(UW-LW)*(V0j(2)-V0j(1))
                    V1=V1j(1)+(nWprime0-UW)/(UW-LW)*(V1j(2)-V1j(1))
                    
!                    print *, "v", v0, v1
!                    pause
                    !LL0=V1-V0 
                   LL0=(V1-V0)/abs(min(V1,V0))                       
!                                                              
                                              endif
                                         
                                                    
           endif 
! print *, nwprime0, indj, vprime0      
ENDSUBROUTINE INTERPOLALL

!12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//12.2//
!Calculus of numerical derivative of the Value Function at period tt+1
SUBROUTINE INTERPOLA_DV(ttt,nWPRIME0,dVPRIME0)
  INTEGER :: ttt, INDJ, iiii, hhh, hhh2 
   REAL*8 :: nWPRIME0, UW, LW
   REAL*8 :: dVPRIME0, dv0, dv1, LU(4), DV(3), XTX(2,2), XTY(2,1)
   LOGICAL :: MATILL
   105 FORMAT ( <2> A12, <J+1> f15.11)
  
                   
                    if ( nWPRIME0 .lt. 0.0000000100D0) then
                   
                                                            !Derivative at W(t+1)=0 (no savings at t)
                                                            LW=GRIDW(1,ttt+1)
                                                            UW=GRIDW(2,ttt+1)
                                                            temp1=EVALUEFUNCTION(1,ttt+1,1)
                                                            temp2=EVALUEFUNCTION(2,ttt+1,1)
                                                            dvPRIME0=(temp2-temp1)/(UW-LW)
                    else
        
                          !FINDING THE CORRESPONDING CELL FOR w(tt+1) AT THE GRIDW OF tt+1
                          INDJ=1
                          DO WHILE ( ( nWPRIME0 .gT. GRIDW(INDJ,ttt+1) ) .and. ( INDJ .lt. J) )
                          INDJ=INDJ+1
                          ENDDO
                                	                       !Derivative at W(t+1) if 0<W(t+1)<GRIDW(J,tt) if JJJ < FIRSTJJ
	                                                       IF ( (INDJ .gt. 1) .and. ( INDJ .le. 4)  )THEN
	                                                        LW=GRIDW(INDJ,ttt+1)
                                                            UW=GRIDW(INDJ+1,ttt+1)
                                                            temp1=EVALUEFUNCTION(INDJ,ttt+1,1)
                                                            temp2=EVALUEFUNCTION(INDJ+1,ttt+1,1)
                                                            dvPRIME0=(temp2-temp1)/(UW-LW)
!!                                                            	LU(1)=(GRIDW(1,ttt+1)+GRIDW(2,ttt+1))/2.d0
!!                                                            	LU(2)=(GRIDW(3,ttt+1)+GRIDW(2,ttt+1))/2.d0
!!																dv(1)=(EVALUEFUNCTION(2,ttt+1)-EVALUEFUNCTION(1,ttt+1))/(GRIDW(2,ttt+1)-GRIDW(1,ttt+1))
!!																dv(2)=(EVALUEFUNCTION(3,ttt+1)-EVALUEFUNCTION(2,ttt+1))/(GRIDW(3,ttt+1)-GRIDW(2,ttt+1))
!!																dvPrime0=dv(1)+(dv(2)-dv(1))*(nwprime0-lu(1))/(lu(2)-lu(1))
														   ENDIF
                                         
                                         
                                         
                                                           !Derivative at W(t+1) if 0<W(t+1)<GRIDW(J,tt) if JJJ < FIRSTJJ
	                                                       IF ( ( INDJ .GT. 4)  .AND. (INDJ .LE. J-1)  )THEN
                                                            	LU(1)=(GRIDW(INDJ-2,ttt+1)+GRIDW(INDJ-1,ttt+1))/2.d0
                                                            	LU(2)=(GRIDW(INDJ-1,ttt+1)+GRIDW(INDJ,ttt+1))/2.d0
																LU(3)=(GRIDW(INDJ,ttt+1)+GRIDW(INDJ+1,ttt+1))/2.d0
																dv(1)=(EVALUEFUNCTION(INDJ-1,ttt+1,1)-EVALUEFUNCTION(INDJ-2,ttt+1,1))/(GRIDW(INDJ-1,ttt+1)-GRIDW(INDJ-2,ttt+1))
																dv(2)=(EVALUEFUNCTION(INDJ,ttt+1,1)-EVALUEFUNCTION(INDJ-1,ttt+1,1))/(GRIDW(INDJ,ttt+1)-GRIDW(INDJ-1,ttt+1))
																dv(3)=(EVALUEFUNCTION(INDJ+1,ttt+1,1)-EVALUEFUNCTION(INDJ,ttt+1,1))/(GRIDW(INDJ+1,ttt+1)-GRIDW(INDJ,ttt+1))
																
																if (nWPRIME0 .GT. LU(2)) then
																dvPrime0=dv(2)+(dv(3)-dv(2))*(nwprime0-lu(2))/(lu(3)-lu(2))
																else
																dvPrime0=dv(1)+(dv(2)-dv(1))*(nwprime0-lu(1))/(lu(2)-lu(1))
																endif
														   ENDIF
                                                        
                                                           !Derivative at W(t+1) if W(t+1)>=GRIDW(J,tt) (extrapolation)
                                                           !linear extrapolation corrected by the factor exp(-(W(t+1)-UW)/((W(t+1)+UW)/2))
                                                             IF ( ( INDJ .eq. J)  )THEN
                                                            	LU(1)=GRIDW(INDJ-2,ttt+1)
                                                            	LU(2)=GRIDW(INDJ-1,ttt+1)
																LU(3)=GRIDW(INDJ,ttt+1)
																dv(1)=(EVALUEFUNCTION(INDJ-1,ttt+1,1)-EVALUEFUNCTION(INDJ-2,ttt+1,1))/(LU(2)-LU(1))
																dv(2)=(EVALUEFUNCTION(INDJ,ttt+1,1)-EVALUEFUNCTION(INDJ-1,ttt+1,1))/(LU(3)-LU(2))
																if ( nWPRIME0 .lt.  GRIDW(J,ttt+1)) then 
																dvPrime0=dv(1)+(dv(2)-dv(1))*(nwprime0-LU(2))/(LU(3)-LU(2))
																else
																dvPrime0=dv(2)+(dv(2)-dv(1))*(nwprime0-LU(3))/(LU(3)-LU(2))*dexp((nWPRIME0-lu(3))/((nWPRIME0+lu(3))/2.D0))                
                                                                endif
                                                              ENDIF !( INDJ .ge. J) 
                                         
!                      print *, INDJ                                     
					endif !( nWPRIME0 .lt. 0.0000000100D0)

ENDSUBROUTINE INTERPOLA_DV




!13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//13//
!DEFINE CONSTRAINTS FOR THE DIFERENT ALGORITHMS

SUBROUTINE BOUNDSTYPE0 !Define constraints for the UNRESTRICTED optimization problem
      
               !Set lower bounds on variables
               
                XLB(1) = 0.00000000D0
                XLB(2) = 0.00D0
              !Set upper bounds on variables
                XUB(1) = 0.95d0
                XUB(2) = 1.00000D0
                      ISTAT=0
			       
					  IPARAM(1)=0
					
ENDSUBROUTINE BOUNDSTYPE0

SUBROUTINE BOUNDSTYPE2 !Define constraints for the UNRESTRICTED optimization problem
        
               !Set lower bounds on variables
                XLB(1) = 0.00000000D0
                XLB(2) = 0.000D0
              !Set upper bounds on variables
                XUB(1) = 0.95d0
                XUB(2) = 1.00000D0
                      ISTAT=0
			          
					  IPARAM(1)=0
					
ENDSUBROUTINE BOUNDSTYPE2

SUBROUTINE BOUNDSTYPE3
!CORNER SOLUTION AT ALPHA=0
              
               !Set lower bounds on variables
                XLB(1) = 0.00000000D0
                XLB(2) = 0.00D0
              !Set upper bounds on variables
                XUB(1) = 0.95d0
                XUB(2) = 0.00000D0
                      ISTAT=0
			         
					  IPARAM(1)=0
					 
ENDSUBROUTINE BOUNDSTYPE3	

DOUBLE PRECISION function feli(a,r1)
        
        REAL*8  :: a, r1, f2, b
        b=0.001d0
                   
                if ( ( r1 .gt. 1d0 ) .and. (a .gt. b ) )then
         feli=(1d0/(1d0-r1))*dexp((1d0-r1)*dlog(a) )
        else
        if ( ( r1 .gt. 1d0 ) .and. (a .le. b ) )then
         feli=(1d0/(1d0-r1))*dexp((1d0-r1)*dlog(b))-100d0*(b-a)
        else
        endif
        endif
        RETURN



!        ENTRY FUMG(a,r1)
!        if ( ( r1 .gt. 1d0 ) .and. (a .gt. b ) )then
!        FUMG=dexp((-r1)*dlog(a) )
!        else
!    if ( ( r1 .gt. 1d0 ) .and. (a .le. b) )then
!        FUMG=dexp((-r1)*dlog(b))+100000d0*(b-a)
!        endif 
!        endif
!        !print *, "fumg", r1, a, fumg, feli, f2
!        !pause 
!        RETURN

ENDFUNCTION feli


END PROGRAM PORTFOLIO
