ps9.f95

MODULE PS9PARA
    IMPLICIT NONE
    REAL(KIND=8), PARAMETER:: PH=5.0, PL=0.0 ! PRICES
    REAL(KIND=8), PARAMETER:: ML=0.0, MH=15.0 ! MASS OF FIRMS
    REAL(KIND=8), PARAMETER:: THETA=0.64, BETA=0.8, GAMMA=0.5772156649
    REAL(KIND=8)::ENT_MASS=100.
    REAL(KIND=8), PARAMETER:: ENT_COST=5., FIX_COST=10.
    REAL(KIND=8),PARAMETER:: A=1./200.
    REAL(KIND=8), DIMENSION(5),PARAMETER:: ENT_DIST=(/0.37,0.4631,0.1102,0.0504,0.0063/) !V(S)
    REAL(KIND=8), DIMENSION(5),PARAMETER:: STATES=(/3.98E-4,3.58,6.82,12.18,18.79/) !5 STATES
    REAL(KIND=8), DIMENSION(5,5),PARAMETER:: TRANS=RESHAPE((/0.6598,0.1997,0.2,0.2,0.2,0.26,0.7201,0.2,0.2,&
    0.2,0.0416,0.042,0.5555,0.2502,0.25,0.0331,&
    0.0326,0.0344,0.3397,0.34,0.0055,0.0056,0.0101,0.0101,0.01/), (/5,5/))
    INTEGER, PARAMETER:: NZ=SIZE(STATES)
    REAL(KIND=8), PARAMETER:: CRIT=1E-2
    REAL(KIND=8), DIMENSION(1):: MUL_TMP ! A VARIABLE TO STORE RESULT OF MATMUL
END MODULE

MODULE PS9RES
    USE PS9PARA
    IMPLICIT NONE
    REAL(KIND=8):: P_HIGH=PH, P_LOW=PL, M_LOW=ML, M_HIGH=MH
    REAL(KIND=8):: P, M
    REAL(KIND=8):: ERROR_P, ERROR_LAB
    REAL(KIND=8), DIMENSION(NZ):: LAB_DEM
    REAL(KIND=8):: LAB_SUP, M_EXIT
    REAL(KIND=8):: AGG_LAB_DEM, INC_LAB_DEM, ENT_LAB_DEM, LAB_ENT_RATIO
    REAL(KIND=8), DIMENSION(NZ):: ENT_VF=0., ENT_VF_NEW, MOTIONS, MOTIONS_NEW
    INTEGER(KIND=8), DIMENSION(NZ):: ENT_PF
END MODULE

MODULE PS9TVRES
    USE PS9PARA
    IMPLICIT NONE
    REAL(KIND=8):: ALPHA
    REAL(KIND=8), DIMENSION(NZ, 2):: SIGMA1, SIGMA2, SIGMA
    REAL(KIND=8), DIMENSION(NZ):: PF_TV1, PF_TV2, VF_TV, VF_TV_NEW

END MODULE

PROGRAM PS9
    USE PS9PARA
    USE PS9RES
    USE PS9TVRES
    IMPLICIT NONE
    INTEGER:: IDX

    CALL STANDARD()
    PRINT*, "P, M_INC, M_ENT, M_EXIT, AGGLAB, INCLAB, ENTLAB, %ENTLAB"
    PRINT*, P, SUM(MOTIONS), M, M_EXIT, AGG_LAB_DEM, INC_LAB_DEM, ENT_LAB_DEM, ENT_LAB_DEM/AGG_LAB_DEM
    PRINT*, "================================================"
    PRINT*, "TV ALPHA=1"
    ALPHA=1.
    P_HIGH=PH
    P_LOW=PL
    M_LOW=ML
    M_HIGH=MH
    CALL TV()
    PRINT*, "P, M_INC, M_ENT, M_EXIT, AGGLAB, INCLAB, ENTLAB, %ENTLAB"
    PRINT*, P, SUM(MOTIONS), M, M_EXIT, AGG_LAB_DEM, INC_LAB_DEM, ENT_LAB_DEM, ENT_LAB_DEM/AGG_LAB_DEM
    PRINT*, "================================================"
    PRINT*, "TV ALPHA=2"
    ALPHA=2.
    P_HIGH=PH
    P_LOW=PL
    M_LOW=ML
    M_HIGH=MH
    CALL TV()
    PRINT*, "P, M_INC, M_ENT, M_EXIT, AGGLAB, INCLAB, ENTLAB, %ENTLAB"
    PRINT*, P, SUM(MOTIONS), M, M_EXIT, AGG_LAB_DEM, INC_LAB_DEM, ENT_LAB_DEM, ENT_LAB_DEM/AGG_LAB_DEM

    OPEN(UNIT=23, FILE='DECISION', STATUS='REPLACE')
        DO IDX=1, NZ
            WRITE(UNIT=23, FMT=*) ENT_PF(IDX), SIGMA1(IDX, 1), SIGMA2(IDX, 1)
        ENDDO
    CLOSE(UNIT=23)
    OPEN(UNIT=24, FILE='PRODUCTIVITY', STATUS='REPLACE')
        DO IDX=1, NZ
            WRITE(UNIT=24, FMT=*) STATES(IDX)
        ENDDO
    CLOSE(UNIT=24)
END PROGRAM

SUBROUTINE STANDARD()
    USE PS9PARA
    USE PS9RES
    IMPLICIT NONE

    ERROR_P=100.
    P = 0.8
    DO WHILE (ERROR_P>CRIT) ! ENTRANT CONDITION
        CALL ENTRY_CLEAR()
    ENDDO

    ERROR_LAB=100.
    M=1. ! AS STATED IN THE PROBLEM SET, INITIALIZE
    DO WHILE (ERROR_LAB>CRIT)
        CALL LAB_CLEAR(.TRUE.)
    ENDDO
END SUBROUTINE

SUBROUTINE ENTRY_CLEAR()
    USE PS9PARA
    USE PS9RES
    IMPLICIT NONE
    INTEGER:: IDX, NEXTI
    INTEGER:: DEBUG=0

    REAL(KIND=8):: NSTAR ! OPT EMPLOY
    REAL(KIND=8):: ERROR_VFI
    REAL(KIND=8):: VF_STAY, VF_EXIT
    REAL(KIND=8):: ENT_COND

    ERROR_VFI=100.
    ENT_VF = 0.
    DO WHILE(ERROR_VFI>CRIT)
        VF_STAY = 0.
        VF_EXIT=0.
        DO IDX=1, NZ
            NSTAR = (THETA*P*STATES(IDX))**(1/(1-THETA))
            LAB_DEM(IDX)=NSTAR
            VF_EXIT = P*STATES(IDX)*(NSTAR**THETA)-NSTAR-P*FIX_COST
            VF_STAY = VF_EXIT
            DO NEXTI=1,NZ
                 VF_STAY = VF_STAY+ BETA*TRANS(IDX, NEXTI)*ENT_VF(NEXTI)
            ENDDO
            ENT_VF_NEW(IDX)=MAX(VF_EXIT, VF_STAY)
            MUL_TMP = MAXLOC((/VF_STAY, VF_EXIT/))
            ENT_PF(IDX) = MUL_TMP(1) -1.
        ENDDO
        ERROR_VFI= MAXVAL(ABS(ENT_VF-ENT_VF_NEW))
        ENT_VF = ENT_VF_NEW
    ENDDO

    ! CALCUALTE ENTRY CONDITIONS
    ENT_MASS = 0.
    DO IDX=1, NZ
        ENT_MASS= ENT_MASS + ENT_DIST(IDX)*ENT_VF(IDX)
    ENDDO
    ENT_COND = ENT_MASS/(P*ENT_COST)-1

    IF (ENT_COND>0.) THEN
        P_HIGH = P
    ELSE
        P_LOW = P
    ENDIF
    P = (P_LOW+P_HIGH)*0.5
    ERROR_P= ABS(ENT_COND)
END SUBROUTINE


SUBROUTINE LAB_CLEAR(STAND)
    USE PS9PARA
    USE PS9RES
    USE PS9TVRES
    IMPLICIT NONE
    INTEGER:: IDX, NIDX
    REAL(KIND=8):: ERROR_VFI, PI
    LOGICAL, INTENT(IN):: STAND
    REAL(KIND=8), DIMENSION(NZ):: DECISION

    IF (STAND) THEN
        DECISION = ENT_PF
    ELSE
        DECISION = SIGMA(:,1)
    ENDIF

    ERROR_VFI=100.
    MOTIONS=0.1 ! JUST RANDOM GUESS
    DO WHILE(ERROR_VFI> CRIT)
        M_EXIT = 0.
        MOTIONS_NEW = 0.
        DO IDX=1,NZ
            DO NIDX=1,NZ
                MOTIONS_NEW(NIDX) = MOTIONS_NEW(NIDX)+(1-DECISION(IDX))*TRANS(NIDX, IDX)*MOTIONS(NIDX)+&
                (1-DECISION(NIDX))*TRANS(NIDX,IDX)*M*ENT_DIST(NIDX)
                ! MASS OF EXIT IS NOT HOW MUCH STATE 1 (BECAUSE THEY HAVE ALREADY EXIT), BUT HOW MUCH MASS GOING TO STAGE ONE.
                M_EXIT = M_EXIT+ DECISION(NIDX)*TRANS(IDX,NIDX)*(MOTIONS(IDX)+M*ENT_DIST(IDX))
            ENDDO
        ENDDO
        ERROR_VFI = MAXVAL(ABS(MOTIONS_NEW-MOTIONS))
        MOTIONS = MOTIONS_NEW
    ENDDO
    AGG_LAB_DEM=0.
    ENT_LAB_DEM= 0.
    INC_LAB_DEM = 0.
    DO IDX=1,NZ
        ! LABOR DEMAND AGGREGATE
        AGG_LAB_DEM= AGG_LAB_DEM+ MOTIONS(IDX)*LAB_DEM(IDX)+ M*ENT_DIST(IDX)*LAB_DEM(IDX)
        ! LABOR DEMAND ENTRANT
        ENT_LAB_DEM = ENT_LAB_DEM+M*ENT_DIST(IDX)*LAB_DEM(IDX)
        ! LABOR DEMAND INCUMBENT
        INC_LAB_DEM= INC_LAB_DEM+MOTIONS(IDX)*LAB_DEM(IDX)
    ENDDO

    PI=0.
    ! LABOR SUPPLY
    DO IDX=1, NZ
        PI =PI+(P*STATES(IDX)*(LAB_DEM(IDX)**THETA)-LAB_DEM(IDX)-P*FIX_COST)*MOTIONS(IDX)+&
        M*(P*STATES(IDX)*(LAB_DEM(IDX)**THETA)-LAB_DEM(IDX)-P*ENT_COST)*ENT_DIST(IDX)
    ENDDO

    LAB_SUP = (P/A)-PI ! BASED ON FOC

    ERROR_LAB = ABS(LAB_SUP-AGG_LAB_DEM)
    IF (LAB_SUP>AGG_LAB_DEM) THEN
        M_LOW = M
    ELSE
        M_HIGH = M
    ENDIF
    M=(M_HIGH+M_LOW)/2.
END SUBROUTINE

SUBROUTINE TV()
    USE PS9PARA
    USE PS9RES
    USE PS9TVRES
    IMPLICIT NONE
    REAL(KIND=8):: ERROR_TV_VFI, NSTAR, ERROR_TV, ENT_COND
    INTEGER:: IDX, NIDX
    REAL(KIND=8), DIMENSION(NZ):: VF_STAY, VF_EXIT
    REAL(KIND=8), DIMENSION(2):: EXP_VAL

    ERROR_TV = 100.
    P = 0.5
    DO WHILE (ERROR_TV> CRIT)
        VF_TV =1. ! INITIAL GUESS
        ERROR_TV_VFI=100.

        DO WHILE (ERROR_TV_VFI>CRIT) ! ENTRANT CONDITION
            DO IDX=1, NZ
                NSTAR = (THETA*P*STATES(IDX))**(1/(1-THETA))
                LAB_DEM(IDX)=NSTAR
                VF_EXIT(IDX) = P*STATES(IDX)*(NSTAR**THETA)-NSTAR-P*FIX_COST
                VF_STAY(IDX) = VF_EXIT(IDX)
                DO NIDX=1,NZ
                    VF_STAY(IDX) = VF_STAY(IDX) + BETA*TRANS(IDX,NIDX)*VF_TV(NIDX)
                ENDDO
                EXP_VAL = (/EXP(ALPHA*VF_EXIT(IDX)),EXP(ALPHA*VF_STAY(IDX))/)
                VF_TV_NEW(IDX) = (GAMMA/ALPHA)+LOG(SUM(EXP_VAL))/ALPHA
                SIGMA(IDX,:) = EXP_VAL/SUM(EXP_VAL) ! DECISION RULES 1-EXIT 2-STAY
            ENDDO
            ERROR_TV_VFI = MAXVAL(ABS(VF_TV-VF_TV_NEW))
            VF_TV = VF_TV_NEW
        ENDDO

        ! CALCUALTE ENTRY CONDITIONS
        ENT_MASS = 0.
        DO IDX=1, NZ
            ENT_MASS= ENT_MASS + ENT_DIST(IDX)*VF_TV(IDX)
        ENDDO
        ENT_COND = ENT_MASS/(P*ENT_COST)-1

        IF (ENT_COND>0.) THEN
            P_HIGH = P
        ELSE
            P_LOW = P
        ENDIF
        P = (P_LOW+P_HIGH)*0.5
        ERROR_TV= ABS(ENT_COND)
    ENDDO
    ! EXACTLY THE SAME
    ERROR_LAB=100.
    M=1. ! AS STATED IN THE PROBLEM SET, INITIALIZE
    DO WHILE (ERROR_LAB>CRIT)
        CALL LAB_CLEAR(.FALSE.)
    ENDDO


    ! JUST FOR SAVING AND PLOT GRAPH
    IF (ALPHA==1.) THEN
        SIGMA1 = SIGMA
    ELSE
        SIGMA2 = SIGMA
    ENDIF
END SUBROUTINE