MODULE mod_transition

  IMPLICIT NONE

  !Other internal variables
  REAL(8) :: svcorigcost, sellcost, sellcost_hh, sellcost_REO, buycost, mortrev, mortcost, maint, maint_hh, maint_REO, cons_rent, cons_hh, bonds_hh, netmortpay, pl, sh, hownstock, nethousebuy, REOsales, adjustcost
  REAL(8), DIMENSION(:), ALLOCATABLE :: Jreo_tomorrow

  !********* Variables for stationary equilibria *********!!!
  REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: mu3own_init, bown_init, mu3own_start, bown_start, bown_today
  INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: locmown_init, locmown_start, locmown_today
  REAL(8), DIMENSION(:,:), ALLOCATABLE :: mu3rent_init, brent_init, mu3rent_start, brent_start, brent_today
  REAL(8), DIMENSION(:), ALLOCATABLE :: htot_reo_init

  !************ Variables for transition path ************!!!
  REAL(8), DIMENSION(:,:,:,:,:), ALLOCATABLE :: qm0_t, default_t, psell_t, Wown_t, Vown_t, bown_sp1_weight_t
  INTEGER, DIMENSION(:,:,:,:,:), ALLOCATABLE :: locmown_t, locbown_sp1_low_t, locbown_sp1_upp_t
  REAL(8), DIMENSION(:,:,:), ALLOCATABLE :: Wrent_t, Vrent_t
  REAL(8), DIMENSION(:,:), ALLOCATABLE :: Jreo_t, series, series_full, prob_ss_steady
  REAL(8), DIMENSION(:), ALLOCATABLE :: w_t, tfp_t, rf_t, pi_t, pinew_t, nom_i_t, rscost_t, ph_t, phnew_t, lumpsum_t, lumpsumnew_t, svcorigcost_t, govspend_t, bonds_hh_t, cons_hh_t, cons_hhnew_t, cons_rent_t, buycost_t, &
       sellcost_t, sellcost_hh_t, sellcost_REO_t, maint_t, maint_hh_t, maint_REO_t, sh_t, pl_t, adjustcost_t, aggdemand_t, aggsupply_t, taylorshock_t, dividends_t, mu_s_init, mu_s_new, price_t, pricenew_t, price_steady, &
       StockMarket_t, haircut_t, haircutnew_t, nominal_debt_t, ChinaDemand_nominal_t, ChinaPayment_t, ChinaShare_t, BondDemand_t
  INTEGER, DIMENSION(:), ALLOCATABLE :: phase_t
  REAL(8) :: hstock_prev, qm0price, max_ltv_steady, rocost_steady, min_xs_steady, curv_ms_steady, min_xb_steady, max_xb_steady, MBS_steady, ChinaDemand_steady, nom_i_steady, qbDebt_Steady
  INTEGER :: T, T_init, T_fwdguide, it, it_actual, in

  !House price loop variables
  REAL(8) :: toler_excessdemand, min_dph, hsupply, hdemand, hactual, supply, demand, diff, excess_demand1, excess_demand0, excessdemand_upp, excessdemand_low, dph, supply_upp, &
       supply_low, demand_upp, demand_low, phweight, ph_low, ph_upp, haircut_new, ChinaPay_Steady
  REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: mu3own_low, mu3own_upp
  REAL(8), DIMENSION(:,:), ALLOCATABLE :: mu3rent_low, mu3rent_upp, htot_reo_temp
  REAL(8), DIMENSION(:), ALLOCATABLE :: htot_reo_new_low, htot_reo_new_upp
  INTEGER :: id, max_iterph, sign_excessdemand1, sign_excessdemand0, iter_ph

  !******************** Panel output ********************!!!
  REAL(8), DIMENSION(:,:), ALLOCATABLE :: cons_nt, earn_nt, nonlaborincome_nt, xsell_nt, psell_nt, sell_nt, default_nt, foreclose_nt, phsim_nt, rfsim_nt, rmsim_nt, pisim_nt, wsim_nt, bonds_nt

  !Other internal variables
  REAL(8), DIMENSION(:,:), ALLOCATABLE :: rand_psell_nt, rand_foreclose_nt, rand_pbuy_nt, rand_hseparate_nt, rand_hdeprec_nt, rand_survive_nt, rand_weight_nt, rand_s_nt, rand_e_nt, rand_all_nt
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: iy_nt, im_nt, ig_nt, is_nt, ie_nt
  INTEGER :: num_id, itsim

  !******* External functions to be called *******
  REAL(8), EXTERNAL :: agglaborsupply, wagedetermine

CONTAINS

  SUBROUTINE sub_solvetransition

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !This program computes perfect foresight "impulse responses" of the economy

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime, isprime

    !************ Variables for transition path ************!!!
    REAL(8) :: toler_ph, toler_pi, toler_price, toler_lumpsum, toler_haircut, err_ph, err_price, err_lumpsum, increment_ph, increment_price, increment_lumpsum, err_haircut, increment_haircut, increment_conshh, toler_conshh, err_conshh
    INTEGER :: max_itertransition, iter_transition

    !Other internal variables
    REAL(8), DIMENSION(5) :: initialize
    INTEGER :: numvars, index

    !********** Variables related to smoothing (((((()))))) !
    INTEGER :: 	iopt, nest, spline_order, lwrk, ier, num_knots, smooth_m
    INTEGER, ALLOCATABLE :: iwrk(:)

    REAL(8), ALLOCATABLE :: weights_t(:), t_array(:), c_array(:), time_t(:), wrk(:), tempprice_t(:), temppricenew_t(:)
    REAL(8) :: smoothing, ssr, realT



    IF (genseries == 1) THEN
       max_itertransition = 1
    ELSE
       max_itertransition = 100
    ENDIF

    IF (algorithm == 4) THEN
       !Set dimensions of panel simulation (length = T_init + T)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_dimensions.txt')
       READ (11, *) T_init
       READ (11, *) num_id
       CLOSE(11)
       T_init = MAX(T_init,1)
    ENDIF

    pi_steady = pi

    !Load transition path
    CALL sub_loadtransition


    DO it = 1, T
       PRINT*,'t =',it,'C =',cons_hh_t(it)
    ENDDO

    iopt = 0
    smooth_m = T + 5
    spline_order = 3
    smoothing = 10.0d0
    nest = smooth_m+spline_order+1 + 2*spline_order+2
    lwrk = (smooth_m*(spline_order+1)+nest*(7+3*spline_order))+100

    ALLOCATE(time_t(T + 5))
    ALLOCATE(tempprice_t(T + 5))
    ALLOCATE(temppricenew_t(T + 5))
    ALLOCATE(t_array(nest))
    ALLOCATE(c_array(nest))
    ALLOCATE(wrk(lwrk))
    ALLOCATE(iwrk(nest))

    ALLOCATE(weights_t(T + 5))


    DO it = 1, T + 5
       time_t(it) = REAL(it,8)
       weights_t(it) = 1.0d0
       !         print*,it,time_t(it)
    ENDDO
    !       pause
    realT=dble(time_t(T+1))
    ! price_t(1)=1.0d0+pi_t(1)
    ! DO it=2,T+1
    !    price_t(it)=price_t(it)*(1.0d0+pi_t(it))
    ! ENDDO



    !********* Compute initial stationary equilibrium *********!!!

    !Allocate variables for initial equilibrium
    ALLOCATE(htot_reo_init(nh))
    ALLOCATE(mu3own_init(ny,nm,ng,ns))
    ALLOCATE(mu3rent_init(ny,ns))
    ALLOCATE(bown_init(ny,nm,ng,ns))
    ALLOCATE(brent_init(ny,ns))
    ALLOCATE(locmown_init(ny,nm,ng,ns))
    ALLOCATE(mu3own_start(ny,nm,ng,ns))
    ALLOCATE(mu3rent_start(ny,ns))
    ALLOCATE(bown_start(ny,nm,ng,ns))
    ALLOCATE(brent_start(ny,ns))
    ALLOCATE(locmown_start(ny,nm,ng,ns))
    ALLOCATE(bown_today(ny,nm,ng,ns))
    ALLOCATE(brent_today(ny,ns))
    ALLOCATE(locmown_today(ny,nm,ng,ns))
    ALLOCATE(mu_s_init(ns))
    ALLOCATE(mu_s_new(ns))

    !Compute equilibrium
    CALL sub_equilibrium_steady
    CALL sub_marketclear_transition
    input_tfp = 1
    PRINT*,'Steady state computed'
!!!Calculate initial values for time series
    IF (genseries == 1) THEN
       numvars = 25
       ALLOCATE(series_full(T + 1,numvars))
       it = 0
       CALL sub_timeseries
    ENDIF
    nominal_debt_t = price_steady*govdebt
    rf_steady = rf
    nom_i_steady  = (1.d0 + rf_steady)*(1.d0 + pi_steady) - 1.d0
    qbDebt_Steady = qb*govdebt
    print*,'rf =',rf_steady,'pi =',pi_steady,'i=',nom_i_steady
    print*,'qb=',qb,'1/1+rf=',1.0d0/(1.0d0+rf),'1+pi/1+i',(1.0d0+pi_steady)/(1.0d0+nom_i_steady)
    print*,'SS trans = ',lumpsum
    print*,'Real trans=',labortax*w*agglaborsupply(w) - (govspend + govdebt*(1.0d0-qb)),', Nom trans=',labortax*w*agglaborsupply(w) - (govspend + (nominal_debt_t(1) - nominal_debt_t(2)/((1.0d0 + nom_i_steady)*(1.0d0 + pi_steady)))/price_steady(1))
    print*, labortax*w*agglaborsupply(w) - (govspend + (nominal_debt_t(1) - nominal_debt_t(2)/((1.0d0 + nom_i_steady)*(1.0d0 + pi_steady)))/price_steady(1))
    print*,govdebt,nominal_debt_t(1)/price_steady(1)
    print*,govdebt*(1.0d0 - qb),(nominal_debt_t(1) - (1.0d0+pi_steady)*nominal_debt_t(1)/(1.0d0 + nom_i_steady))/price_steady(1)
    print*,govdebt*(1.0d0 - qb),(nominal_debt_t(1) - (1.0d0+pi_steady)*nominal_debt_t(2)/(1.0d0 + nom_i_steady))/price_steady(1)
    print*,govdebt*(1.0d0 - qb),(nominal_debt_t(1) - (1.0d0 + pi_steady)*nominal_debt_t(2)/((1.0d0 + nom_i_steady)*(1.0d0 + pi_steady)))/price_steady(1)



!    stop





    ALLOCATE(series(T + 1,17))
    series(1,1)  = ph
    series(1,2)  = (1.d0 + rf)**4.d0 - 1.d0
    series(1,3)  = (1.d0 + pi)**4.d0 - 1.d0
    series(1,4)  = w
    series(1,5)  = lumpsum
    series(1,6)  = aggsupply
    series(1,7)  = aggdemand
    series(1,8)  = cons_hh
    series(1,9)  = pa*cons_rent
    series(1,10) = govspend
    series(1,11) = sh
    series(1,12) = pl
    series(1,13) = buycost
    series(1,14) = sellcost
    series(1,15) = svcorigcost
    series(1,16) = maint
    series(1,17) = adjustcost

    !Store objects from initial equilibrium
    mu_s_init     = mu_s
    mu3own_init	  = mu3own
    mu3rent_init  = mu3rent
    bown_init	  = bown
    brent_init	  = brent
    locmown_init  = locmown
    htot_reo_init = htot_reo
    ! pi_steady     = pi
    rf_steady     = rf
    rf_t(1)       = rf_steady
    MBS_steady    = MBS
    ChinaDemand_steady = qb*ChinaDemand*(govdebt + MBS_steady)
    ChinaDemand_nominal_t = price_steady*ChinaDemand_steady
    ChinaPay_Steady = (1.0d0-qb)*ChinaDemand*(govdebt + MBS_steady)

    IF(use_steady_guess .eq. 1) THEN
       cons_hh_t    = cons_hh
!       haircut_t(1) = haircut / 10.0d0
    ENDIF
    print*,haircut
    DO it = 1, T
       PRINT*,'t =',it,'C =',cons_hh_t(it)
    ENDDO


    ! PRINT*,'Checking mutual fund'
    ! PRINT*,'Receipts: ',mortrev + dividends + govdebt + StockMarketValue*qb*tfp*glabor
    ! PRINT*,'Liabilities: ', psurvive*bonds_hh + ChinaDemand_steady
    ! PRINT*,'Check stock:', StockMarketValue*tfp*glabor, qb*StockMarketValue*tfp*glabor+dividends
    ! PRINT*,qb*mortrev,mortrev,mortcost
    ! PRINT*,qb*mortrev - mortcost
    ! PAUSE
    nom_i_t(1)    = (1.d0 + rf_steady)*(1.d0 + pi_steady) - 1.d0

    IF (deprec_h > 0.d0) THEN
       hstock_prev = yh/deprec_H
    ELSE
       hstock_prev = hstock
    ENDIF
    print*, minval(locmown_init)
    IF(minval(locmown_init) .LT. 1) THEN
      PRINT*,'locmown_init fucked'
      PAUSE
    ENDIF

    !Rescale TFP
    DO it = 1, T
       tfp_t(it) = tfp*tfp_t(it)
    ENDDO
    tfp_t(T + 1) = tfp_t(T)

    IF (algorithm == 4) THEN
       !Initialize household panel
       CALL sub_initialize_panel
    ENDIF

    !********* Compute (actual or initially perceived) final stationary equilibrium *********!!!

!!!UNCOMMENT OUT AND ***EDIT*** THIS SECTION TO ALLOW FOR DIFFERENT INITIAL AND FINAL STEADY STATES
!!!NOTE: sub_equilibrium_steady

    !Set paths for exogenous transition variables
    !		tfp    = tfp_t(T + 1)
    !		rscost = rscost_t(T + 1)
    !		!Load house price guess
    !		open(unit=11,	file='Input'//sep//'Final'//sep//'initialize.txt')
    !		read (11, *) initialize
    !		close(11)
    !		ph = initialize(3)
    !		pi = initialize(4)
    !		rf = initialize(5)
    !		pi = (1.d0 + pi)**0.25d0 - 1.d0
    !		rf = (1.d0 + rf)**0.25d0 - 1.d0
    !		!Compute equilibrium
    !		equil_num = 2
    !		call sub_equilibrium_steady

    !Allocate variables for transition path
    ALLOCATE(Jreo_t(nh,T + 1))
    ALLOCATE(qm0_t(nm,nb,ng,ns,T + 1))
    ALLOCATE(default_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(psell_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(locmown_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(locbown_sp1_low_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(locbown_sp1_upp_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(bown_sp1_weight_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(Wown_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(Wrent_t(ny,ns,T + 1))
    ALLOCATE(Vown_t(ny,nm,ng,ns,T + 1))
    ALLOCATE(Vrent_t(ny,ns,T + 1))
    IF (algorithm == 2 .OR. algorithm == 3 .OR. algorithm == 5) THEN
       ALLOCATE(mu3own_low(ny,nm,ng,ns))
       ALLOCATE(mu3rent_low(ny,ns))
       ALLOCATE(mu3own_upp(ny,nm,ng,ns))
       ALLOCATE(mu3rent_upp(ny,ns))
       ALLOCATE(htot_reo_temp(nh,num_threads_simul))
       ALLOCATE(htot_reo_new_low(nh))
       ALLOCATE(htot_reo_new_upp(nh))
    ENDIF

    !Store objects for transition path
    ph_t(T + 1)						 = ph
    phnew_t(T + 1)					 = ph
    pi_t(T + 1)						 = pi
    pinew_t(T + 1)                   = pi
    w_t(T + 1)						 = w
    aggsupply_t(T + 1)				 = aggsupply
    aggdemand_t(T + 1)				 = aggdemand
    lumpsum_t(T + 1)				 = lumpsum
    lumpsumnew_t(T + 1)				 = lumpsum
    bonds_hh_t(T + 1)				 = bonds_hh
    cons_rent_t(T + 1)				 = cons_rent
    cons_hh_t(T + 1)				 = cons_hh
    maint_t(T + 1)					 = maint
    maint_hh_t(T + 1)				 = maint_hh
    maint_REO_t(T + 1)				 = maint_REO
    sellcost_t(T + 1)				 = sellcost
    sellcost_hh_t(T + 1)			 = sellcost_hh
    sellcost_REO_t(T + 1)			 = sellcost_REO
    buycost_t(T + 1)				 = buycost
    svcorigcost_t(T + 1)			 = svcorigcost
    adjustcost_t(T + 1)				 = 0.d0
    dividends_t(T + 1)				 = dividends
    StockMarket_t(T + 1)                         = dividends/(1.0d0 - qb)
    sh_t(T + 1)						 = sh
    pl_t(T + 1)						 = pl
    qm0_t(:,:,:,:,T + 1)			 = qm0
    Jreo_t(:,T + 1)				     = Jreo
    default_t(:,:,:,:,T + 1)         = default
    psell_t(:,:,:,:,T + 1)           = psell
    locmown_t(:,:,:,:,T + 1)         = locmown
    locbown_sp1_low_t(:,:,:,:,T + 1) = locbown_sp1_low
    locbown_sp1_upp_t(:,:,:,:,T + 1) = locbown_sp1_upp
    bown_sp1_weight_t(:,:,:,:,T + 1) = bown_sp1_weight
    Wown_t(:,:,:,:,T + 1)			 = Wown
    Wrent_t(:,:,T + 1)				 = Wrent
    Vown_t(:,:,:,:,T + 1)			 = Vown
    Vrent_t(:,:,T + 1)				 = Vrent

    !Initialize price level paths
    ! price_steady(1) = 1.d0 + pi_steady
    ! price_t(1)      = 1.d0 + pi_t(1)
    ! DO it = 2, T + 1
    !    price_steady(it) = price_steady(it - 1)*(1.d0 + pi_steady)
    !    price_t(it)      = price_t(it - 1)*(1.d0 + pi_t(it))
    ! ENDDO
    pricenew_t(T + 1) = price_t(T + 1)

    !************ Compute transition path ************!!!

    IF (fwdguide == 1) THEN
       !Determine length of time where Taylor_pi = 0
       T_fwdguide = 0
       DO it = 1, T
          IF (ABS(taylorshock_t(it)) > 1e-6) THEN
             T_fwdguide = it - 1
             EXIT
          ENDIF
       ENDDO
    ENDIF

    !Compute transition path of house prices, transfers, and inflation
    IF (algorithm == 2) THEN
       toler_ph          = 1e-4
       toler_pi          = 1e-6
       toler_price       = 1e-5
       toler_lumpsum     = 1e-1
       toler_haircut     = 2e-4
       toler_conshh      = 1.0d0 !1e-4
       increment_ph 	   = 0.5d0
       IF(shoot_pi .eq. 1) THEN
         increment_price = 0.5d0
       ELSE
         increment_price = 0.2d0
       ENDIF
       increment_lumpsum = 0.5d0
       increment_haircut = 0.0001d0
       increment_conshh  = 0.5d0 !Test to see if that's the issue !0.9d0
    ELSEIF (algorithm == 3 .OR. algorithm == 5) THEN
       toler_ph          = 1e-4
       toler_pi          = 1e-6
       toler_price       = 1e-5
       toler_lumpsum     = 1e-4
       toler_haircut     = 2e-4
       toler_conshh      = 1e-4
       increment_ph 	   = 0.5d0
       increment_price   = 0.5d0
       increment_lumpsum = 0.5d0
       increment_haircut = 0.0001d0
       increment_conshh  = 0.5d0
    ENDIF
    PRINT*,'Starting Transition'
    PRINT*,'Algorithm = ',algorithm
    IF (algorithm == 2 .OR.algorithm == 3 .OR. algorithm == 5) THEN
!!!Solve for converged transition path
       DO iter_transition = 1, max_itertransition
!!!Back solve
          PRINT*,'Backsolve VFN and qm0'
          CALL sub_backsolve
!!!Forward iterate
          PRINT *, ' '
          PRINT *, 'Forward iterate, finding equilibrium prices conditional on future beliefs'
          CALL sub_forwarditerate
          IF(shoot_pi .eq. 0) THEN
             err_price   = MAXVAL(ABS(pricenew_t - price_t))
          ELSE
             err_price   = MAXVAL(ABS(pinew_t - pi_t))
          ENDIF
          err_lumpsum = MAXVAL(ABS(lumpsumnew_t - lumpsum_t))
          err_ph      = MAXVAL(ABS(phnew_t - ph_t))
          err_haircut = MAXVAL(ABS(haircutnew_t - haircut_t))
          err_conshh  = MAXVAL(ABS(cons_hhnew_t - cons_hh_t))
          it = 1
          PRINT*,'t = 1, HC =',haircut_t(it),' HCnew=',haircutnew_t(it)
          IF(shoot_pi .eq. 1) THEN
            DO it = 1, T
               PRINT*,'t =',it,' pi =',pi_t(it),' pinew_t=',pinew_t(it)
            ENDDO
          ELSE
            DO it = 1, T
               PRINT*,'t =',it,' price =',price_t(it),' pricenew_t=',pinew_t(it),'price_s =',price_steady(it)
            ENDDO
          ENDIF

          PRINT *, ' '
          WRITE(*, '(a20, i4, a10, f11.8, a15, f11.8, a10, f11.8)'), '   iter_transition =' ,iter_transition, ', err_price =' ,err_price, ', err_lumpsum =' ,err_lumpsum, ', err_ph =' ,err_ph, ', err_hc =',err_haircut, ', err_cons =',err_conshh
          IF (iter_transition == 1) THEN
             OPEN(unit=11,	file='log.txt', status = 'replace')
             WRITE(11, '(a20, i4, a13, f11.8, a15, f11.8, a10, f11.8)'), '   iter_transition =' ,iter_transition, ', err_price =' ,err_price, ', err_lumpsum =' ,err_lumpsum, ', err_ph =' ,err_ph
             CLOSE(11)
          ELSE
             OPEN(unit=11,	file='log.txt', access = 'append')
             WRITE(11, '(a20, i4, a13, f11.8, a15, f11.8, a10, f11.8)'), '   iter_transition =' ,iter_transition, ', err_price =' ,err_price, ', err_lumpsum =' ,err_lumpsum, ', err_ph =' ,err_ph
             CLOSE(11)
          ENDIF
          !Print output
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pi.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), ((1.d0 + pi_t(it))**4.d0 - 1.d0, it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_i.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), ((1.d0 + nom_i_t(it))**4.d0 - 1.d0, it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_r.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), ((1.d0 + rf_t(it))**4.d0 - 1.d0, it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_price.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (price_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pricenew.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (pricenew_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pricedev.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (price_t(it)/price_steady(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_lumpsum.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (lumpsum_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_ph.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (ph_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pinew.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), ((1.d0 + pinew_t(it))**4.d0 - 1.d0, it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_lumpsumnew.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (lumpsumnew_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_phnew.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (phnew_t(it), it = 1,T)
          CLOSE(11)
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_aggc.txt', status = 'replace')
          WRITE(11, '(1f20.10)'), (cons_hh_t(it), it = 1,T)
          CLOSE(11)
          !Print time series
          DO it = 1, T
             series(it + 1,1)  = ph_t(it)
             series(it + 1,2)  = (1.d0 + rf_t(it + 1))**4.d0 - 1.d0
             series(it + 1,3)  = (1.d0 + pi_t(it))**4.d0 - 1.d0
             series(it + 1,4)  = w_t(it)
             series(it + 1,5)  = lumpsum_t(it)
             series(it + 1,6)  = aggsupply_t(it)
             series(it + 1,7)  = aggdemand_t(it)
             series(it + 1,8)  = cons_hh_t(it)
             series(it + 1,9)  = pa*cons_rent_t(it)
             series(it + 1,10) = govspend_t(it)
             series(it + 1,11) = sh_t(it)
             series(it + 1,12) = pl_t(it)
             series(it + 1,13) = buycost_t(it)
             series(it + 1,14) = sellcost_t(it)
             series(it + 1,15) = svcorigcost_t(it)
             series(it + 1,16) = maint_t(it)
             series(it + 1,17) = adjustcost_t(it)
          ENDDO
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'series.txt', status = 'replace')
          DO it = 1, T
             WRITE (11, '(17f14.8)') (series(it,index), index = 1,17)
          ENDDO
          CLOSE(11)
          IF ((err_price < toler_price .AND. err_lumpsum < toler_lumpsum .AND. err_ph < toler_ph .AND. err_haircut < toler_haircut .AND. err_conshh < toler_conshh) .OR. iter_transition == max_itertransition) THEN
             EXIT
          ENDIF
          !Update paths for the price level, inflation, lump sum transfers, and house prices
          IF (algorithm == 2) THEN
             ph_t = ph_t + increment_ph*(phnew_t - ph_t)
            !  IF (err_ph < 1e-3) THEN
              IF(err_price .gt. 1.0d-3 .OR. err_ph .gt. 1.0d-3) THEN
                 increment_haircut = 0.1d0
              ! ELSEIF(err_price .gt. 5.0d-5) THEN
              !    increment_haircut = 0.0025d0
                increment_price    = 0.2d0
              ELSE
                 increment_haircut = 0.2d0
                 increment_price   = 0.5d0
              ENDIF

                IF(shoot_pi .eq. 0) THEN
                   price_t   = price_t + increment_price*(pricenew_t - price_t)
                ELSE
                   pi_t   = pi_t      + increment_price*(pinew_t - pi_t)
                ENDIF
                lumpsum_t = lumpsum_t + increment_lumpsum*(lumpsumnew_t - lumpsum_t)

                haircut_t = haircut_t + increment_haircut*(haircutnew_t - haircut_t)
                cons_hh_t = cons_hh_t + increment_conshh*(cons_hhnew_t - cons_hh_t)
            !  ENDIF
          ELSEIF (algorithm == 3) THEN
             IF(shoot_pi .eq. 0) THEN
                price_t   = price_t + increment_price*(pricenew_t - price_t)
             ELSE
                pi_t   = pi_t      + increment_price*(pinew_t - pi_t)
             ENDIF
             lumpsum_t = lumpsum_t + increment_lumpsum*(lumpsumnew_t - lumpsum_t)
             ph_t      = ph_t + increment_ph*(phnew_t - ph_t)
             IF(err_price .gt. 1.0d-3) THEN
                increment_haircut = 0.0001d0
             ELSEIF(err_price .gt. 5.0d-5) THEN
                increment_haircut = 0.0025d0
             ELSE
                increment_haircut = 0.05d0
             ENDIF

             haircut_t = haircut_t + increment_haircut*(haircutnew_t - haircut_t)
             cons_hh_t = cons_hh_t + increment_conshh*(cons_hhnew_t - cons_hh_t)
           ELSEIF (algorithm == 5) THEN
              IF(shoot_pi .eq. 0) THEN
                 price_t   = price_t + increment_price*(pricenew_t - price_t)
              ELSE
                 pi_t   = pi_t      + increment_price*(pinew_t - pi_t)
              ENDIF
              lumpsum_t = lumpsum_t + increment_lumpsum*(lumpsumnew_t - lumpsum_t)
              ph_t      = ph_t + increment_ph*(phnew_t - ph_t)
              IF(err_price .gt. 1.0d-3 .OR. err_ph .gt. 1.0d-3) THEN
                 increment_haircut = 0.1d0
              ELSEIF(err_price .gt. 1.0d-4 .OR. err_ph .gt. 1.0d-4) THEN
                 increment_haircut = 0.25d0
              ELSE
                 increment_haircut = 0.5d0
              ENDIF

              haircut_t = haircut_t + increment_haircut*(haircutnew_t - haircut_t)
              cons_hh_t = cons_hh_t + increment_conshh*(cons_hhnew_t - cons_hh_t)
          ENDIF
          price_t(T + 1) = price_steady(T + 1)

          IF(shoot_pi .eq. 0) THEN
             pi_t(1) = price_t(1) - 1.d0
             DO it = 2, T + 1
                pi_t(it) = price_t(it)/price_t(it - 1) - 1.d0
             ENDDO
          ELSE
             price_t(1) = 1.0d0 + pi_t(1)
             DO it = 2, T + 1
                price_t(it) = price_t(it - 1)*(1.0d0 + pi_t(it))
             ENDDO
          ENDIF

          IF(smooth_guess .eq. 1 .AND. shoot_pi .eq. 0) THEN

            tempprice_t(1:T+1) = price_t(1:T+1)
            DO it = 2, 5
              tempprice_t(T + it) = tempprice_t(T + it - 1)*(1.0d0 + pi_steady)
            ENDDO
            realT = dble(T) + 5.0d0
            print*,iopt
            print*,nest
            print*,price_t
            print*,time_t(1)
            print*,time_t(T+1)
            print*,spline_order
            print*,smooth_m
            print*,smoothing
            print*,lwrk
            print*,realT

            smoothing = 50.0d0
            ! CALL curfit(iopt,smooth_m,time_t,price_t,weights_t,time_t(1),time_t(T + 5),spline_order,smoothing,nest,num_knots,t_array,c_array,ssr,wrk,lwrk,iwrk,ier)
            CALL curfit(iopt,smooth_m,time_t,tempprice_t,weights_t,time_t(1),time_t(T + 5),spline_order,smoothing,nest,num_knots,t_array,c_array,ssr,wrk,lwrk,iwrk,ier)

            PRINT*,'ier=',ier

            CALL splev(t_array,num_knots,c_array,spline_order,time_t,temppricenew_t,smooth_m,ier)
            PRINT*,'ier=',ier



            OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_smoothprice.txt', status = 'replace')
            WRITE(11, '(1f20.10)'), (temppricenew_t(it), it = 1,T+1)
            CLOSE(11)
            OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_nonsmoothprice.txt', status = 'replace')
            WRITE(11, '(1f20.10)'), (tempprice_t(it), it = 1,T+1)
            CLOSE(11)

            price_t = price_t + increment_price*(temppricenew_t(1:T) - price_t(1:T))
            pi_t(1) = price_t(1) - 1.d0
            DO it = 2, T + 1
               pi_t(it) = price_t(it)/price_t(it - 1) - 1.d0
            ENDDO

            ! smoothing = 25.0d0
            !
            ! CALL curfit(iopt,smooth_m,time_t,lumpsum_t,weights_t,0.0d0,realT,spline_order,smoothing,nest,num_knots,t_array,c_array,ssr,wrk,lwrk,iwrk,ier)
            !
            ! PRINT*,'ier=',ier
            !
            ! CALL splev(t_array,num_knots,c_array,spline_order,time_t,temppricenew_t,smooth_m,ier)
            ! PRINT*,'ier=',ier
            !
            !
            !
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_smoothlumpsum.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), (temppricenew_t(it), it = 1,T+1)
            ! CLOSE(11)
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_nonsmoothlumpsum.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), (lumpsum_t(it), it = 1,T+1)
            ! CLOSE(11)
            !
            ! lumpsum_t = temppricenew_t
            !
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pinext.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), ((1.d0 + pi_t(it))**4.d0 - 1.d0, it = 1,T)
            ! CLOSE(11)
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pricenext.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), (price_t(it), it = 1,T)
            ! CLOSE(11)
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_lumpsumnext.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), (lumpsum_t(it), it = 1,T)
            ! CLOSE(11)
            ! OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_phnext.txt', status = 'replace')
            ! WRITE(11, '(1f20.10)'), (ph_t(it), it = 1,T)
            ! CLOSE(11)
          ENDIF

       ENDDO


       IF (genseries == 1) THEN
!!!Print files
          OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'series_full.txt', status = 'replace')
          DO it = 1, T + 1
             WRITE (11, '(25f14.8)') (series_full(it,index), index = 1,numvars)
          ENDDO
          CLOSE(11)
       ENDIF
    ELSEIF (algorithm == 4) THEN
!!!Simulate panel
       CALL sub_backsolve
       PRINT *, ' '
       PRINT *, 'Simulate panel'
       bown     = bown_init
       brent    = brent_init
       locmown  = locmown_init
       DO it = 1, T
          itsim = T_init + it
          CALL sub_transition_phase(phase_t(it))
          pi         = pi_t(it)
          ph         = ph_t(it)
          rf	       = rf_t(it + 1)
          qb	       = 1.d0/(1.d0 + rf_t(it + 1))
          qb_annuity = psurvive*qb
          qm         = 1.d0/(1.d0 + pi_t(it + 1))*1.d0/(1.d0 + rf_t(it + 1))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(it + 1))
          qm0	       = qm0_t(:,:,:,:,it)
          w          = w_t(it + 1)
          dividends  = dividends_t(it + 1)
          lumpsum    = lumpsum_t(it + 1)
          Jreo       = Jreo_t(:,it + 1)
          Wown	   = Wown_t(:,:,:,:,it + 1)
          Wrent      = Wrent_t(:,:,it + 1)
          Vown       = Vown_t(:,:,:,:,it + 1)
          Vrent      = Vrent_t(:,:,it + 1)
          CALL sub_backsolve_reo
          CALL sub_backsolve_ylb
          CALL sub_backsolve_valfn
          !Simulate panel from periods t to t + 1
          CALL sub_simulate
          WRITE(*, '(a7, i4)'), ' tsim =', itsim
       ENDDO
       !Print panel output
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_out.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 2, T_init + T
             IF (ig_nt(in,itsim) == 0) THEN
                WRITE(11, '(i7, a2 i5, a2, f10.4, a2, f10.4, a2, f10.4, a2, i5, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4)'), in, ', ', itsim - 1, ', ', yy(iy_nt(in,itsim)), ', ', mm(im_nt(in,itsim)), ', ', 0.d0, ', ', is_nt(in,itsim), ', ', cons_nt(in,itsim), ', ', earn_nt(in,itsim), ', ', nonlaborincome_nt(in,itsim), ', ', xsell_nt(in,itsim), ', ', psell_nt(in,itsim), ', ', sell_nt(in,itsim), ', ', default_nt(in,itsim), ', ', foreclose_nt(in,itsim), ', ', wsim_nt(in,itsim), ', ', phsim_nt(in,itsim), ', ', rfsim_nt(in,itsim), ', ', rmsim_nt(in,itsim), ', ', pisim_nt(in,itsim)
             ELSE
                WRITE(11, '(i7, a2 i5, a2, f10.4, a2, f10.4, a2, f10.4, a2, i5, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4, a2, f10.4)'), in, ', ', itsim - 1, ', ', yy(iy_nt(in,itsim)), ', ', mm(im_nt(in,itsim)), ', ', hh(hfun(ig_nt(in,itsim))), ', ', is_nt(in,itsim), ', ', cons_nt(in,itsim), ', ', earn_nt(in,itsim), ', ', nonlaborincome_nt(in,itsim), ', ', xsell_nt(in,itsim), ', ', psell_nt(in,itsim), ', ', sell_nt(in,itsim), ', ', default_nt(in,itsim), ', ', foreclose_nt(in,itsim), ', ', wsim_nt(in,itsim), ', ', phsim_nt(in,itsim), ', ', rfsim_nt(in,itsim), ', ', rmsim_nt(in,itsim), ', ', pisim_nt(in,itsim)
             ENDIF
          ENDDO
       ENDDO
       CLOSE(11)
    ENDIF

  END SUBROUTINE sub_solvetransition

  SUBROUTINE sub_loadtransition

    USE mod_globvar

    !Set transition length
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_T.txt')
    READ (11, *) T
    CLOSE(11)


!!!Allocate transition path variables
    !Endogenous
    ALLOCATE(Jreo_tomorrow(nh))
    ALLOCATE(aggsupply_t(T + 1))
    ALLOCATE(aggdemand_t(T + 1))
    ALLOCATE(ph_t(T + 1))
    ALLOCATE(phnew_t(T + 1))
    ALLOCATE(pi_t(T + 1))
    ALLOCATE(pinew_t(T + 1))
    ALLOCATE(price_t(T + 1))
    ALLOCATE(pricenew_t(T + 1))
    ALLOCATE(price_steady(T + 1))
    ALLOCATE(w_t(T + 1))
    ALLOCATE(rf_t(T + 1))
    ALLOCATE(nom_i_t(T + 1))
    ALLOCATE(lumpsum_t(T + 1))
    ALLOCATE(lumpsumnew_t(T + 1))
    ALLOCATE(bonds_hh_t(T + 1))
    ALLOCATE(cons_rent_t(T + 1))
    ALLOCATE(cons_hh_t(T + 1))
    ALLOCATE(cons_hhnew_t(T + 1))
    ALLOCATE(maint_t(T + 1))
    ALLOCATE(maint_hh_t(T + 1))
    ALLOCATE(maint_REO_t(T + 1))
    ALLOCATE(sellcost_t(T + 1))
    ALLOCATE(sellcost_hh_t(T + 1))
    ALLOCATE(sellcost_REO_t(T + 1))
    ALLOCATE(buycost_t(T + 1))
    ALLOCATE(svcorigcost_t(T + 1))
    ALLOCATE(adjustcost_t(T + 1))
    ALLOCATE(dividends_t(T + 1))
    ALLOCATE(sh_t(T + 1))
    ALLOCATE(pl_t(T + 1))
    ALLOCATE(haircut_t(T + 1))
    ALLOCATE(haircutnew_t(T + 1))
    ALLOCATE(nominal_debt_t(T + 1))
    ALLOCATE(ChinaPayment_t(T + 1))
    ALLOCATE(ChinaDemand_nominal_t(T + 1))
    ALLOCATE(ChinaShare_t(T + 1))
    ALLOCATE(BondDemand_t(T + 1))
    !Exogenous
    ALLOCATE(phase_t(T + 1))
    ALLOCATE(tfp_t(T + 1))
    ALLOCATE(rscost_t(T + 1))
    ALLOCATE(govspend_t(T + 1))
    ALLOCATE(taylorshock_t(T + 1))
    ALLOCATE(StockMarket_t(T + 1))

    !Set relevant transition path variables
    govspend_t = govspend
    !Load relevant transition path variables
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_phase.txt')
    DO it = 1, T
       READ (11, *) phase_t(it)
    ENDDO
    CLOSE(11)
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_ph.txt')
    DO it = 1, T
       READ (11, *) ph_t(it)
    ENDDO


    price_steady(1) = 1.d0 + pi_steady
    DO it = 2, T + 1
       price_steady(it) = price_steady(it - 1)*(1.d0 + pi_steady)
    ENDDO



    IF(shoot_pi .eq. 0) THEN
       OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pricedev.txt')
       DO it = 1, T
          READ (11, *) price
          price_t(it) = price*price_steady(it)
       ENDDO
       CLOSE(11)
       price_t(T + 1) = price_steady(T + 1)

       pi_t(1) = price_t(1) - 1.0d0
       DO it = 2, T + 1
          pi_t(it) = price_t(it)/price_t(it - 1) - 1.0d0
       ENDDO
       DO it = 1, T + 1
          print*,it,'pi_t = ',pi_t(it),'price_t = ',price_t(it)
       ENDDO
    ELSE
       OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_pi.txt')
       DO it = 1, T
          READ (11, *) pi_t(it)
       ENDDO
       CLOSE(11)
       !Convert to quarterly
       DO it = 1, T
          pi_t(it) = (1.d0 + pi_t(it))**0.25d0 - 1.d0
       ENDDO

       price_t(1) = 1.0d0 + pi_t(1)
       DO it = 2, T
          price_t(it) = price_t(it - 1)*(1.0d0 + pi_t(it))
       ENDDO
    ENDIF

    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_tfp.txt')
    DO it = 1, T
       READ (11, *) tfp_t(it)
    ENDDO
    CLOSE(11)
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_rscost.txt')
    DO it = 1, T
       READ (11, *) rscost_t(it)
    ENDDO
    CLOSE(11)
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_lumpsum.txt')
    DO it = 1, T
       READ (11, *) lumpsum_t(it)
    ENDDO
    CLOSE(11)
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_aggc.txt')
    DO it = 1, T
       READ (11, *) cons_hh_t(it)
    ENDDO
    CLOSE(11)
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_taylorshock.txt')
    DO it = 1, T
       READ (11, *) taylorshock_t(it)
    ENDDO
    CLOSE(11)
    !Read in the guess for the haircut
    OPEN(unit=11,	file='Input'//sep//'Transition'//sep//'transition_haircut.txt')
    DO it = 1, T
       READ (11, *) haircut_t(it)
    ENDDO
    CLOSE(11)




    rscost_t(T + 1)      = rscost_t(T)
    lumpsum_t(T + 1)     = lumpsum_t(T)
    taylorshock_t(T + 1) = taylorshock_t(T)
    phase_t(T + 1)       = phase_t(T)

!    print*,'rscost',rscost
!    print*,'rscost_t',rscost_t

!    print*,'rscost',lumpsum
!    print*,'rscost_t',lumpsum_t

!    print*,'rscost',pi
!    print*,'rscost_t',pi_t

!    print*,'rscost',ph
!    print*,'rscost_t',ph_t


    IF(use_steady_guess .EQ. 1) THEN
       rscost_t = rscost
       lumpsum_t = lumpsum
       pi_t = pi
       ph_t = ph
       price_t = price_steady
!       haircut_t = 0.0d0
!       haircut = 0.0d0
       cons_hh_t = gAggC
    ENDIF



    !Store steady state phase variables
    ALLOCATE(prob_ss_steady(ns,ns))
    prob_ss_steady = prob_ss
    rocost_steady  = rocost
    max_ltv_steady = max_ltv
    min_xs_steady  = min_xs
    curv_ms_steady = curv_ms
    min_xb_steady  = min_xb
    max_xb_steady  = max_xb

  END SUBROUTINE sub_loadtransition


  SUBROUTINE sub_backsolve

    USE mod_globvar

    !Initialize relevant prices, value functions, and policy functions
    qm0				= qm0_t(:,:,:,:,T + 1)
    Jreo			= Jreo_t(:,T + 1)
    default			= default_t(:,:,:,:,T + 1)
    psell			= psell_t(:,:,:,:,T + 1)
    locmown			= locmown_t(:,:,:,:,T + 1)
    locbown_sp1_low = locbown_sp1_low_t(:,:,:,:,T + 1)
    locbown_sp1_upp = locbown_sp1_upp_t(:,:,:,:,T + 1)
    bown_sp1_weight = bown_sp1_weight_t(:,:,:,:,T + 1)
    Wown			= Wown_t(:,:,:,:,T + 1)
    Wrent			= Wrent_t(:,:,T + 1)
    Vown			= Vown_t(:,:,:,:,T + 1)
    Vrent			= Vrent_t(:,:,T + 1)
    PRINT *, ' '
    PRINT *, 'Back solve'
    DO it = 1, T
       it_actual = T - it + 1
       !Set transition phase
       CALL sub_transition_phase(phase_t(it_actual))
       !Compute real rates using the inflation path, Taylor rule, and Fisher equation
       IF (fwdguide == 1) THEN
          IF (it_actual <= T_fwdguide) THEN
             taylor_pi = 0.d0
          ELSE
             taylor_pi = 1.25d0
          ENDIF
       ENDIF
       IF(inf_taylor .eq. 1) THEN
        !  nom_i_t(it_actual + 1) = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi_t(it_actual))/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it_actual)) - 1.d0)
         nom_i_t(it_actual + 1) = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it_actual))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((1.d0 + pi_t(it_actual))/(1.d0 + pi_steady))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it_actual))) - 1.d0)
       ELSE
        !  nom_i_t(it_actual + 1) = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price_t(it_actual))/(1.d0 + price_steady(it_actual)))**taylor_pi*exp(taylorshock_t(it_actual)) - 1.d0)
          nom_i_t(it_actual + 1) = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it_actual))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((price_t(it_actual))/(price_steady(it_actual)))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it_actual))) - 1.d0)
       ENDIF

       rf_t(it_actual + 1)    = (1.d0 + nom_i_t(it_actual + 1))/(1.d0 + pi_t(it_actual + 1)) - 1.d0
       !Calculate wages and aggregate supply using the intermediate goods firm pricing equation
!!!NOTE: need path of labor shock distribution for agglaborsupply in the case with time-varying prob_ss transitions. Right now sum(mudist(:,is)) = stationary(is) always.
       w_t(it_actual)         = wagedetermine(pi_t(it_actual + 1),pi_t(it_actual),rf_t(it_actual + 1),tfp_t(it_actual),aggsupply_t(it_actual + 1),w_t(it_actual + 1),cons_hh_t(it_actual))
       gaggC                  = cons_hh_t(it_actual)
       aggsupply_t(it_actual) = tfp_t(it_actual)*agglaborsupply(w_t(it_actual))
       !Calculate adjustment costs and dividends
       IF(as_if_rigid .EQ. 0) THEN
          adjustcost_t(it_actual) = adjust_theta/2.d0*(pi_t(it_actual) - pi_steady)**2.d0*aggsupply_t(it_actual)
       ELSE
          adjustcost_t(it_actual) = 0.0d0
       ENDIF
       dividends_t(it_actual)  = aggsupply_t(it_actual) - w_t(it_actual)*agglaborsupply(w_t(it_actual)) - adjustcost_t(it_actual) - firmFC
     ENDDO

     DO it = 1, T
        it_actual = T - it + 1
        !Set transition phase
        CALL sub_transition_phase(phase_t(it_actual))

        IF(debt_nominal .eq. 0) THEN
          IF(it_actual .GT. 1) THEN
            lumpsum_t(it_actual) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend - rf_t(it_actual)*qbDebt_Steady
          ELSE
            lumpsum_t(it_actual) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + ((pi_t(1) - nom_i_t(1))/(1.d0 + pi_t(1)))*qbDebt_Steady
          ENDIF
        ELSE
          IF(it_actual .GT. 1) THEN
            !lumpsum_t(it_actual) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + price_steady(it_actual)*qbDebt_Steady/price_t(it_actual) - (1.d0+nom_i_t(it))*price_steady(it_actual - 1)*qbDebt_Steady / price_t(it_actual)
            ! lumpsum_t(it) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + (price_steady(it_actual) - (1.d0+nom_i_t(it))*price_steady(it_actual - 1) / price_t(it_actual))*qbDebt_Steady
            lumpsum_t(it_actual) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + ( pi_steady - nom_i_t(it_actual) )*(price_steady(it_actual - 1) / price_t(it_actual))*qbDebt_Steady
          ELSE
            lumpsum_t(it_actual) = labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + ((pi_steady - nom_i_t(1)) / price_t(1))*qbDebt_Steady
          ENDIF

       ENDIF
       IF(it_actual .GT. 1) THEN
          PRINT*,'t =',it_actual
          PRINT*,'Real T=',labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend - rf_t(it_actual)*qbDebt_Steady
          PRINT*,'Nom T=',labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + ( pi_steady - nom_i_t(it_actual) )*(price_steady(it_actual - 1) / price_t(it_actual))*qbDebt_Steady
          PRINT*,'Nom T=',labortax*w_t(it_actual)*agglaborsupply(w_t(it_actual)) - govspend + price_steady(it_actual)*qbDebt_Steady/price_t(it_actual) -(1.d0+nom_i_t(it_actual))*price_steady(it_actual - 1)*qbDebt_Steady / price_t(it_actual)
       ENDIF
       !Set today's prices and next period's wage and government transfers
        ph         = ph_t(it_actual)
        qb         = 1.d0/(1.d0 + rf_t(it_actual + 1))
        qb_annuity = psurvive*qb
        qm         = 1.d0/(1.d0 + pi_t(MIN(T + 1,it_actual + 2)))*1.d0/(1.d0 + rf_t(MIN(T + 1,it_actual + 2)))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(MIN(T + 1,it_actual + 2)))	!Continuation interest rates (qmprime)
        w          = w_t(it_actual + 1)
        lumpsum    = lumpsum_t(it_actual + 1)
        dividends  = dividends_t(it_actual + 1)
        glabor     = agglaborsupply(w)
        StockMarket_t(it_actual) = dividends_t(it_actual + 1) + StockMarket_t(it_actual + 1)/(1.0d0 + rf_t(MIN(T + 1,it_actual + 2)))
        !Backsolve (calculates today's mortgage/bond prices and REO/household value/policy functions) and calculate today's wage
        pi     = pi_t(it_actual + 1)			!Inflation between today and next period
        rscost = rscost_t(it_actual + 1)		!Mortgage servicing costs
        CALL sub_backsolve_q
        price  = price_t(it_actual)
        pi     = pi_t(it_actual)				!Inflation between last period and today
        qm     = 1.d0/(1.d0 + pi_t(it_actual + 1))*1.d0/(1.d0 + rf_t(it_actual + 1))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(it_actual + 1))
        CALL sub_backsolve_reo
        CALL sub_backsolve_ylb
        CALL sub_backsolve_valfn
        WRITE(*, '(a4, i4, a10, f8.5, a10, f10.4, a12, f15.6)'), ' t =', it_actual, ', diff_q =' ,diff_q, ', diff_j =' ,diff_j, ', diff_val =' ,diff_val
        !Export transition path objects
        qm0_t(:,:,:,:,it_actual)			 = qm0
        Jreo_t(:,it_actual)					 = Jreo
        default_t(:,:,:,:,it_actual)		 = default
        psell_t(:,:,:,:,it_actual)			 = psell
        locmown_t(:,:,:,:,it_actual)		 = locmown
        locbown_sp1_low_t(:,:,:,:,it_actual) = locbown_sp1_low
        locbown_sp1_upp_t(:,:,:,:,it_actual) = locbown_sp1_upp
        bown_sp1_weight_t(:,:,:,:,it_actual) = bown_sp1_weight
        Wown_t(:,:,:,:,it_actual)			 = Wown
        Wrent_t(:,:,it_actual)				 = Wrent
        Vown_t(:,:,:,:,it_actual)			 = Vown
        Vrent_t(:,:,it_actual)				 = Vrent
    ENDDO

  END SUBROUTINE sub_backsolve

  SUBROUTINE sub_forwarditerate

    USE mod_globvar

    INTEGER :: ig, ih, is, isprime


    IF (algorithm == 2) THEN
!!!For each t, find equilibrium ph and pi conditional on the current guesses for t+1, t+2, etc.
      htot_reo      = htot_reo_init
      mu3own_start  = mu3own_init
      mu3rent_start = mu3rent_init
      bown_start    = bown_init
      brent_start   = brent_init
      locmown_start = locmown_init
      mu_s		  = mu_s_init
      it       = 1
      CALL sub_transition_phase(phase_t(it))
      w         = w_t(it)
      gAggC     = cons_hh_t(it)
      glabor    = agglaborsupply(w)
      aggsupply = tfp*glabor
      adjustcost= adjustcost_t(it)
      lumpsum   = lumpsum_t(it)
      dividends = dividends_t(it)
      haircut   = haircut_t(it)
      CALL sub_distribution_sp31
      !Update labor efficiency distribution from last period to today
      mu_s_new = 0.d0
      DO is = 1, ns
         DO isprime = 1, ns
            mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
         ENDDO
      ENDDO
      mu_s  = mu_s_new
      pi = price_t(1) - 1.d0
      ph            = ph_t(it)
      tfp	          = tfp_t(it)
      rscost        = rscost_t(it)
      rf	          = rf_t(it + 1)
      qb	          = 1.d0/(1.d0 + rf_t(it + 1))
      qb_annuity    = psurvive*qb
      qm            = 1.d0/(1.d0 + pi_t(it + 1))*1.d0/(1.d0 + rf_t(it + 1))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(it + 1))
      qm0	          = qm0_t(:,:,:,:,it)
      Jreo_tomorrow = Jreo_t(:,it + 1)
      CALL sub_supplydemand_firstperiod_fast
      IF (genseries == 1) THEN
      !!!Calculate time series
         CALL sub_timeseries
      ENDIF
      !Update end-of-period housing stock
      hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
      !Calculate stock of houses abandoned due to death of owner
      htot_reo_death = 0.d0
      DO ig = 1, ng
         ih = hfun(ig)
         htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
      ENDDO
      !Calculate REO stock at the beginning of next period
      DO ih = 1, nh
         htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
      ENDDO

       DO it = 2, T
          !Transition of households from subperiod 3 in period t - 1 to subperiod 1 in period t
          CALL sub_transition_phase(phase_t(it))
          w         = w_t(it)
          gAggC     = cons_hh_t(it)
          glabor    = agglaborsupply(w)
          aggsupply = tfp*glabor
          adjustcost= adjustcost_t(it)
          lumpsum   = lumpsum_t(it)
          dividends = dividends_t(it)
          haircut   = haircut_t(it)
          CALL sub_distribution_sp31
          mu_s_new = 0.d0
          DO is = 1, ns
             DO isprime = 1, ns
                mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
             ENDDO
          ENDDO
          mu_s = mu_s_new
          !Find ph that clears the housing market
          IF(shoot_pi .eq. 0) THEN
             IF (it == 1) THEN
                pi = price_t(1) - 1.d0
             ELSE
                pi = price_t(it)/price_t(it - 1) - 1.d0
             ENDIF
          ELSE
             pi = pi_t(it)
          ENDIF

          ph            = ph_t(it)
          tfp	          = tfp_t(it)
          rscost        = rscost_t(it)
          rf	          = rf_t(it + 1)
          qb	          = 1.d0/(1.d0 + rf_t(it + 1))
          qb_annuity    = psurvive*qb
          qm            = 1.d0/(1.d0 + pi_t(it + 1))*1.d0/(1.d0 + rf_t(it + 1))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(it + 1))
          qm0	          = qm0_t(:,:,:,:,it)
          Jreo_tomorrow = Jreo_t(:,it + 1)
          CALL sub_supplydemand_fast
          IF (genseries == 1) THEN
!!!Calculate time series
             CALL sub_timeseries
          ENDIF
          !Update end-of-period housing stock
          hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
          !Calculate stock of houses abandoned due to death of owner
          htot_reo_death = 0.d0
          DO ig = 1, ng
             ih = hfun(ig)
             htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
          ENDDO
          !Calculate REO stock at the beginning of next period
          DO ih = 1, nh
             htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
          ENDDO
       ENDDO
    ELSEIF (algorithm == 3) THEN
!!!For each t, find equilibrium ph and pi conditional on the current guesses for t+1, t+2, etc.
        htot_reo      = htot_reo_init
        mu3own_start  = mu3own_init
        mu3rent_start = mu3rent_init
        bown_start    = bown_init
        brent_start   = brent_init
        locmown_start = locmown_init
        mu_s		  = mu_s_init
        it       = 1
        CALL sub_transition_phase(phase_t(it))
        !Update labor efficiency distribution from last period to today
        mu_s_new = 0.d0
        DO is = 1, ns
           DO isprime = 1, ns
              mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
           ENDDO
        ENDDO
        mu_s  = mu_s_new
        IF (it == 1) THEN
           pi = price_t(1) - 1.d0
        ELSE
           pi = price_t(it)/price_t(it - 1) - 1.d0
        ENDIF
        ph      = ph_t(it)
        lumpsum = lumpsum_t(it)
        CALL sub_supplydemand_firstperiod
        IF (genseries == 1) THEN
        !!!Calculate time series
           CALL sub_timeseries
        ENDIF
        !Store policy functions and distributions
        mu3own_start  = mu3own
        mu3rent_start = mu3rent
        bown_start    = bown
        brent_start   = brent
        locmown_start = locmown
        !Update end-of-period housing stock
        hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
        !Calculate stock of houses abandoned due to death of owner
        htot_reo_death = 0.d0
        DO ig = 1, ng
           ih = hfun(ig)
           htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
        ENDDO
        !Calculate REO stock at the beginning of next period
        DO ih = 1, nh
           htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
        ENDDO
       DO it = 2, T
          !Find ph and pi that clear the housing and bond/goods markets
          CALL sub_transition_phase(phase_t(it))
          !Update labor efficiency distribution from last period to today
          mu_s_new = 0.d0
          DO is = 1, ns
             DO isprime = 1, ns
                mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
             ENDDO
          ENDDO
          mu_s  = mu_s_new
          IF(shoot_pi .eq. 0) THEN
             IF (it == 1) THEN
                pi = price_t(1) - 1.d0
             ELSE
                pi = price_t(it)/price_t(it - 1) - 1.d0
             ENDIF
          ELSE
             pi = pi_t(it)
          ENDIF
          ph      = ph_t(it)
          lumpsum = lumpsum_t(it)
          CALL sub_supplydemand
          IF (genseries == 1) THEN
!!!Calculate time series
             CALL sub_timeseries
          ENDIF
          !Store policy functions and distributions
          mu3own_start  = mu3own
          mu3rent_start = mu3rent
          bown_start    = bown
          brent_start   = brent
          locmown_start = locmown
          !Update end-of-period housing stock
          hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
          !Calculate stock of houses abandoned due to death of owner
          htot_reo_death = 0.d0
          DO ig = 1, ng
             ih = hfun(ig)
             htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
          ENDDO
          !Calculate REO stock at the beginning of next period
          DO ih = 1, nh
             htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
          ENDDO
       ENDDO
     ELSEIF (algorithm == 5) THEN
     !!!For each t, find equilibrium ph and pi conditional on the current guesses for t+1, t+2, etc.
       htot_reo      = htot_reo_init
       mu3own_start  = mu3own_init
       mu3rent_start = mu3rent_init
       bown_start    = bown_init
       brent_start   = brent_init
       locmown_start = locmown_init
       mu_s		  = mu_s_init
       it       = 1
       CALL sub_transition_phase(phase_t(it))
       !Update labor efficiency distribution from last period to today
       mu_s_new = 0.d0
       DO is = 1, ns
          DO isprime = 1, ns
             mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
          ENDDO
       ENDDO
       mu_s  = mu_s_new
       IF (it == 1) THEN
          pi = price_t(1) - 1.d0
       ELSE
          pi = price_t(it)/price_t(it - 1) - 1.d0
       ENDIF
       ph      = ph_t(it)
       lumpsum = lumpsum_t(it)
       CALL sub_supplydemand_firstperiod
       IF (genseries == 1) THEN
       !!!Calculate time series
          CALL sub_timeseries
       ENDIF
       !Store policy functions and distributions
       mu3own_start  = mu3own
       mu3rent_start = mu3rent
       bown_start    = bown
       brent_start   = brent
       locmown_start = locmown
       !Update end-of-period housing stock
       hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
       !Calculate stock of houses abandoned due to death of owner
       htot_reo_death = 0.d0
       DO ig = 1, ng
          ih = hfun(ig)
          htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
       ENDDO
       !Calculate REO stock at the beginning of next period
       DO ih = 1, nh
          htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
       ENDDO

        DO it = 2, T
           !Find ph and pi that clear the housing and bond/goods markets
           CALL sub_transition_phase(phase_t(it))
           !Update labor efficiency distribution from last period to today
           mu_s_new = 0.d0
           DO is = 1, ns
              DO isprime = 1, ns
                 mu_s_new(isprime) = mu_s_new(isprime) + mu_s(is)*prob_ss(is,isprime)
              ENDDO
           ENDDO
           mu_s  = mu_s_new
           IF(shoot_pi .eq. 0) THEN
              IF (it == 1) THEN
                 pi = price_t(1) - 1.d0
              ELSE
                 pi = price_t(it)/price_t(it - 1) - 1.d0
              ENDIF
           ELSE
              pi = pi_t(it)
           ENDIF
           ph      = ph_t(it)
           lumpsum = lumpsum_t(it)
           CALL sub_supplydemand_slow
           IF (genseries == 1) THEN
     !!!Calculate time series
              CALL sub_timeseries
           ENDIF
           !Store policy functions and distributions
           mu3own_start  = mu3own
           mu3rent_start = mu3rent
           bown_start    = bown
           brent_start   = brent
           locmown_start = locmown
           !Update end-of-period housing stock
           hstock_prev = hstock_prev*(1.d0 - deprec_h) + yh
           !Calculate stock of houses abandoned due to death of owner
           htot_reo_death = 0.d0
           DO ig = 1, ng
              ih = hfun(ig)
              htot_reo_death(ih) = htot_reo_death(ih) + SUM(mu3own(:,:,ig,:))*(1.d0 - psurvive)*(1.d0 - deprec_h_death)
           ENDDO
           !Calculate REO stock at the beginning of next period
           DO ih = 1, nh
              htot_reo(ih) = (htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))*(1.d0 - deprec_h) + htot_reo_death(ih)
           ENDDO
        ENDDO
    ENDIF
    price_t(T + 1)    = price_steady(T + 1)
    pricenew_t(T + 1) = price_steady(T + 1)
    cons_hhnew_t(T + 1) = cons_hh_t(T + 1)
    ! pricenew_t(T + 1) = pricenew_t(T)*(1.d0 + pinew_t(T + 1))

  END SUBROUTINE sub_forwarditerate

  SUBROUTINE sub_supplydemand_firstperiod

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime

    !************ Variables for transition path ************!!!

    !Government transfer loop variables
    REAL(8) :: toler_lumpsum, dlumpsum, taxrev, lumpsumnew, err_lumpsum, err_cons, toler_cons, toler_haircut, err_haircut
    INTEGER :: max_iterlumpsum, iter_lumpsum, iter_haircut, max_iterhaircut

    !Inflation loop variables
    REAL(8) :: toler_excessdemand_pi, min_dpi, pi_low, pi_upp, supply_pi, supply_pi_upp, supply_pi_low, demand_pi, demand_pi_low, demand_pi_upp, excess_demand_pi1, excess_demand_pi0, excessdemand_pi_upp, &
         excessdemand_pi_low, dpi, diff_pi, piweight, nom_i, cons_hh_today, haircut_low, haircut_upp
    REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: mu3own_pi_low, mu3own_pi_upp
    REAL(8), DIMENSION(:,:), ALLOCATABLE :: mu3rent_pi_low, mu3rent_pi_upp
    REAL(8), DIMENSION(:), ALLOCATABLE :: htot_reo_new_pi_low, htot_reo_new_pi_upp
    INTEGER :: max_iterpi, iter_pi, sign_excessdemand_pi1, sign_excessdemand_pi0, max_itercons

    !Other internal variables
    REAL(8) :: pi_today, w_today, lumpsum_today, dividends_today, pi_tomorrow, w_tomorrow, rscost_tomorrow, qm_tomorrow, lumpsum_tomorrow, dividends_tomorrow, StockMarket_today

    !Allocate variables
    ALLOCATE(mu3own_pi_low(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_low(ny,ns))
    ALLOCATE(mu3own_pi_upp(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_upp(ny,ns))
    ALLOCATE(htot_reo_new_pi_low(nh))
    ALLOCATE(htot_reo_new_pi_upp(nh))

!!!Find ph, pi, and transfers that satisfy the government budget constraint and clear the housing and bond/goods markets
    max_iterpi            = 100
    max_itercons          = 5
    max_iterph            = 100
    IF(algorithm .lt. 5) THEN
      max_iterhaircut       = 1 !20
    ELSE
      max_iterhaircut = 10
    ENDIF
    toler_excessdemand_pi = 1e-8
    toler_excessdemand    = 1e-5
    toler_cons            = 1e-5
    toler_haircut         = 1e-5
    dlumpsum              = 0.9d0
    min_dpi               = 1e-8
    min_dph               = 1e-6
    !Initialize
    IF (fwdguide == 1) THEN
       IF (it <= T_fwdguide) THEN
          taylor_pi = 0.d0
       ELSE
          taylor_pi = 1.25d0
       ENDIF
    ENDIF
    tfp                = tfp_t(it)
    pi_today           = pi
    cons_hh_today      = cons_hh_t(it)
    lumpsum_today      = lumpsum
    qm_tomorrow        = 1.d0/(1.d0 + pi_t(MIN(T + 1,it + 2)))*1.d0/(1.d0 + rf_t(MIN(T + 1,it + 2)))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(MIN(T + 1,it + 2)))
    w_tomorrow		     = w_t(it + 1)
    lumpsum_tomorrow   = lumpsum_t(it + 1)
    dividends_tomorrow = dividends_t(it + 1)
    rscost_tomorrow    = rscost_t(it + 1)
    Jreo_tomorrow      = Jreo_t(:,it + 1)
    !Inflation loop
    excessdemand_pi_low = 999.d0
    excessdemand_pi_upp = 999.d0
    pi_low              = -0.5d0
    pi_upp              = 0.5d0
    haircut = haircut_t(1)

    DO iter_pi = 1, max_iterpi
       !Compute today's real rate using the Taylor rule and Fisher equation
       pi          = pi_today
       gAggC       = cons_hh_today
       IF (it == 1) THEN
          price = 1.d0 + pi
       ELSE
          price = price_t(it - 1)*(1.d0 + pi)
       ENDIF
       !Guess inflation vs guess price level [uncomment next line for price level]
       IF(shoot_pi .eq. 1) THEN
         pi_tomorrow = pi_t(it + 1)
       ELSE
         pi_tomorrow = price_t(it + 1)/price - 1.d0
       ENDIF

      !  IF(inf_taylor .eq. 1) THEN
      !    nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
      !  ELSE
      !    nom_i       = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price)/(1.d0 + price_steady(it)))**taylor_pi*exp(taylorshock_t(it)) - 1.d0)
      !  ENDIF
       IF(inf_taylor .eq. 1) THEN
         IF(rho_taylor .GT. 0.01) THEN
           nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it))) - 1.d0)
         ELSE
           nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
         ENDIF
       ELSE
          nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((price/(price_steady(it)))**taylor_pi))**(1.0d0-rho_taylor)))*EXP(taylorshock_t(it))) - 1.d0)
       ENDIF

      !  nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
       ! nom_i       = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price)/(1.d0 + price_steady(it)))**taylor_pi*exp(taylorshock_t(it)) - 1.d0)
       rf          = (1.d0 + nom_i)/(1.d0 + pi_tomorrow) - 1.d0
       qb          = 1.d0/(1.d0 + rf)
       qb_annuity  = psurvive*qb
       !Compute mortgage prices
       qm0				= qm0_t(:,:,:,:,it + 1)
       Jreo			= Jreo_t(:,it + 1)
       default			= default_t(:,:,:,:,it + 1)
       psell			= psell_t(:,:,:,:,it + 1)
       locmown			= locmown_t(:,:,:,:,it + 1)
       locbown_sp1_low = locbown_sp1_low_t(:,:,:,:,it + 1)
       locbown_sp1_upp = locbown_sp1_upp_t(:,:,:,:,it + 1)
       bown_sp1_weight = bown_sp1_weight_t(:,:,:,:,it + 1)
       qm              = qm_tomorrow
       w               = w_tomorrow
       lumpsum         = lumpsum_tomorrow
       dividends       = dividends_tomorrow
       pi              = pi_tomorrow
       rscost          = rscost_tomorrow
       CALL sub_backsolve_q
       qm = 1.d0/(1.d0 + pi_tomorrow)*1.d0/(1.d0 + rf)*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_tomorrow)
       !Compute wage and dividends using firm pricing equation
       pi              = pi_today
       w               = wagedetermine(pi_tomorrow,pi,rf,tfp,aggsupply_t(it + 1),w_t(it),cons_hh_today)
       w_today         = w
       aggsupply       = tfp*agglaborsupply(w)
       glabor          = agglaborsupply(w)
       taxrev	         = labortax*w*agglaborsupply(w)
       lumpsum_today   = taxrev - (govspend + govdebt*(1.0d0-qb))

       IF(debt_nominal .eq. 0) THEN
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + (pi_today - nom_i_t(it))*qbDebt_Steady/(1.0d0 + pi_today)
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_today - nom_i_t(1))/(1.d0 + pi_today))*qbDebt_Steady
         ENDIF
       ELSE
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + ( pi_steady - nom_i_t(it) )*(price_steady(it - 1) / price)*qbDebt_Steady
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_steady - nom_i_t(1)) / price)*qbDebt_Steady
         ENDIF

       ENDIF



       IF(as_if_rigid .eq. 0) THEN
          adjustcost      = adjust_theta/2.d0*(pi_today - pi_steady)**2.d0*aggsupply
       ELSE
          adjustcost      = 0.0d0
       ENDIF
       dividends       = aggsupply - w*agglaborsupply(w) - adjustcost - firmFC
       dividends_today = dividends
       StockMarket_today = qb*StockMarket_t(it)
       !House price loop
       excessdemand_low = 999.d0
       excessdemand_upp = 999.d0
       ph_low           = 0.1d0
       ph_upp           = 10.d0
       DO iter_ph = 1, max_iterph
          w         = w_t(it + 1)
          dividends = dividends_t(it + 1)
          lumpsum   = lumpsum_t(it + 1)
          Jreo      = Jreo_t(:,it + 1)
          Wown	  = Wown_t(:,:,:,:,it + 1)
          Wrent     = Wrent_t(:,:,it + 1)
          Vown      = Vown_t(:,:,:,:,it + 1)
          Vrent     = Vrent_t(:,:,it + 1)
          CALL sub_backsolve_reo
          CALL sub_backsolve_ylb
          CALL sub_backsolve_valfn

          haircut_low = 0.0d0
          haircut_upp = 0.1d0
          !Here we figure out what the haircut has to be
          !Code updated 28 June to have the haircut be a part of the outer loop
          !so that now we read in the haircut in sub_loadtransition from the file
          !transition_haircut
          haircut_new = haircut
          do iter_haircut = 1, max_iterhaircut
            ! haircut       = (haircut_low + haircut_upp) / 2.0d0
            haircut = (haircut + haircut_new)/2.0d0
            haircut = max(-0.1d0,min(0.1d0,haircut))

            bown_today    = bown
            brent_today   = brent
            locmown_today = locmown
               !Update distribution from last period to today
             w         = w_today
             lumpsum   = lumpsum_today
             dividends = dividends_today
             mu3own    = mu3own_start
             mu3rent   = mu3rent_start
             bown      = bown_start
             brent     = brent_start
             locmown   = locmown_start
             glabor    = agglaborsupply(w)
             CALL sub_distribution_sp31
             bown    = bown_today
             brent   = brent_today
             locmown = locmown_today
             !Stock of new REO housing
             htot_reo_new  = 0.d0
             htot_reo_temp = 0.d0
             DO is = 1, ns
                DO ig = 1, ng
                   ih = hfun(ig)
                   CALL omp_set_num_threads(num_threads_simul)
                   !$omp parallel private(id,iy)
                   !$omp do
                   DO im = 1, nm
                      id = omp_get_thread_num() + 1
                      DO iy = ny_neg + 1, ny
                         !New REOs
                         htot_reo_temp(ih,id) = htot_reo_temp(ih,id) + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                      ENDDO
                   ENDDO
                   !$omp end do
                   !$omp end parallel
                ENDDO
             ENDDO
             DO id = 1, num_threads_simul
                htot_reo_new = htot_reo_new + htot_reo_temp(:,id)
             ENDDO
             !Distribution transition from subperiod 1 to subperiod 2
             CALL sub_distribution_sp12
             !Distribution transition from subperiod 2 to subperiod 3
             CALL sub_distribution_sp23
             !Calculate market clearing
             CALL sub_marketclear_transition

             ! err_b = ((govdebt + MBS)*(1.0d0-ChinaDemand) + StockMarketValue*tfp*aggL - qb_annuity*bonds_hh/qb)
             !Payout is supposed to be psurvive*bonds_hh, i.e. the
             IF(debt_nominal .EQ. 0) THEN
              !  In first period doesn't matter level or fraction, because coming from SS where they are the same
              !  IF(ChinaLevel == 1) THEN
              !    haircut_new = 1.0d0 - (mortrev + govdebt - bankFC - ChinaDemand_steady + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ELSE
                  haircut_new = 1.0d0 - ((mortrev+govdebt-bankFC)*(1.0d0-ChinaDemand) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ENDIF
             ELSE
              !  IF(ChinaLevel == 1) THEN
              !    haircut_new = 1.0d0 - ((mortrev+price_steady(it)*govdebt/price - bankFC - price_steady(it)*ChinaDemand_steady/price) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ELSE
                 haircut_new = 1.0d0 - ((mortrev+ (1.0d0 + nom_i_t(it))*qbDebt_Steady/price-bankFC)*(1.0d0-ChinaDemand) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ENDIF
             ENDIF
             err_haircut = haircut - haircut_new
             WRITE(*, '(a12, i4, a6, f11.7, a9, f11.7, a7, f11.7, a7, f11.7)'), '   Iter_hc =', iter_haircut, ', hc =' ,haircut, ', hcnew =' , haircut_new, ', Rev =' ,mortrev + dividends_today + govdebt + StockMarket_today, ', Lia =' ,psurvive*bonds_hh + ChinaDemand_steady
            !  print*,'Stockmarket = ',StockMarket_today
            !  print*,'dividends_today = ',dividends_today
            !  print*,'mortrev = ',mortrev
            !  print*,'govdebt = ',govdebt

             IF( abs(err_haircut) .LT. toler_haircut) THEN
               EXIT
             ENDIF
            !  IF(err_haircut .LT. 0.0d0) THEN
            !    haircut =

          ENDDO
!          haircut = 0.0d0
          !Compute housing supply and demand
          CALL sub_housingclear
          IF (ABS(dph) < min_dph .OR. (demand > 0.d0 .AND. supply > 0.d0 .AND. diff < toler_excessdemand) .OR. iter_ph == max_iterph) THEN
             EXIT
          ENDIF
          !Exit in the event of nonexistence
          IF ((ph < ph_low .AND. excess_demand1 < 0.d0) .OR. (ph > ph_upp .AND. excess_demand1 > 0.d0)) THEN
             PRINT *, ' '
             PRINT *, 'Nonexistence!'
             STOP
          ENDIF
          ph = ph + dph
       ENDDO
       IF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          phweight     = excessdemand_low/(excessdemand_low - excessdemand_upp)
          ph			 = (1.d0 - phweight)*ph_low + phweight*ph_upp
          demand       = (1.d0 - phweight)*demand_low + phweight*demand_upp
          supply       = (1.d0 - phweight)*supply_low + phweight*supply_upp
          mu3own       = (1.d0 - phweight)*mu3own_low + phweight*mu3own_upp
          mu3rent      = (1.d0 - phweight)*mu3rent_low + phweight*mu3rent_upp
          htot_reo_new = (1.d0 - phweight)*htot_reo_new_low + phweight*htot_reo_new_upp
       ELSEIF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) > 100.d0) THEN
          ph			 = ph_low
          demand       = demand_low
          supply       = supply_low
          mu3own       = mu3own_low
          mu3rent      = mu3rent_low
          htot_reo_new = htot_reo_new_low
       ELSEIF (ABS(excessdemand_low) > 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          ph			 = ph_upp
          demand       = demand_upp
          supply       = supply_upp
          mu3own       = mu3own_upp
          mu3rent      = mu3rent_upp
          htot_reo_new = htot_reo_new_upp
       ENDIF
       !Compute excess demand [Dont think we need this since we do it in the haircut]
      !  CALL sub_marketclear_transition
       !Check stuff

        IF(debt_nominal .eq. 0) THEN

          BondDemand_t(it)  = mortcost+qbDebt_Steady
          IF(ChinaLevel .eq. 0) THEN
            !Note, this is likely not correct
            ! ChinaPayment_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
            ! ChinaShare_t(it)  = ChinaDemand
          ELSE
            ChinaPayment_t(it) = ChinaDemand*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
            ChinaShare_t(it)  = ChinaDemand_steady / BondDemand_t(it)
          ENDIF
       ELSE
          BondDemand_t(it)  = mortcost+qbDebt_Steady*price_steady(it)/price
          IF(ChinaLevel .eq. 0) THEN
            ! ChinaPayment_t(it) = (mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_tomorrow))/price - (mortrev+price_steady(it)*govdebt/price-bankFC))*ChinaDemand
            ! ChinaShare_t(it)  = ChinaDemand
          ELSE
            ChinaPayment_t(it) = ChinaDemand*(mortrev + (1.0d0 + nom_i_t(it))*qbDebt_Steady/price - bankFC) - ChinaDemand_nominal_t(it)/price
            ChinaShare_t(it)  = (ChinaDemand_nominal_t(it)/price) / BondDemand_t(it)
          ENDIF
        ENDIF
        aggdemand = aggdemand + ChinaPayment_t(it)

       demand_pi         = aggdemand !govdebt
       supply_pi         = aggsupply !bonds_hh
       excess_demand_pi1 = demand_pi - supply_pi
       IF (excess_demand_pi1 >= 0.d0) THEN
          sign_excessdemand_pi1 = 1
          IF (pi > pi_low) THEN
             pi_low              = pi
             demand_pi_low       = demand_pi
             supply_pi_low       = supply_pi
             excessdemand_pi_low = excess_demand_pi1
             mu3own_pi_low       = mu3own
             mu3rent_pi_low      = mu3rent
             htot_reo_new_pi_low = htot_reo_new
          ENDIF
       ELSE
          sign_excessdemand_pi1 = -1
          IF (pi < pi_upp) THEN
             pi_upp              = pi
             demand_pi_upp       = demand_pi
             supply_pi_upp       = supply_pi
             excessdemand_pi_upp = excess_demand_pi1
             mu3own_pi_upp       = mu3own
             mu3rent_pi_upp      = mu3rent
             htot_reo_new_pi_upp = htot_reo_new
          ENDIF
       ENDIF
       IF (iter_pi >=2 ) THEN
          IF (excess_demand_pi0 >= 0.d0) THEN
             sign_excessdemand_pi0 = 1
          ELSE
             sign_excessdemand_pi0 = -1
          ENDIF
       ENDIF
       IF (iter_pi == 1) THEN
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          !diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),MIN(0.2d0,diff_pi)*MAX(0.01d0,ABS(pi))*sign_excessdemand_pi1))		!Change min(0.2d0,... to adjust initial increment
          !				dpi     = max(0.5d0*(pi_low - pi),min(0.5d0*(pi_upp - pi),0.1d0*max(0.01d0,abs(pi))*sign_excessdemand_pi1))						!Change min(0.1d0,... to adjust initial increment
       ELSEIF (sign_excessdemand_pi1*sign_excessdemand_pi0 == -1) THEN
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),-0.5d0*dpi))
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
!          diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*ABS(pi)/ABS(dpi)
          ENDIF
       ELSE
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),dpi))
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
!          diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*ABS(pi)) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))/ABS(dpi)
          ENDIF
       ENDIF
       excess_demand_pi0 = excess_demand_pi1
       WRITE(*, '(a12, i4, a6, f11.7, a9, f11.7, a6, f11.7, a6, f11.7, a6, f12.8, a6, f12.8)'), '   Iter_pi =', iter_pi, ', pi =' ,(1.d0 + pi)**4.d0 - 1.d0, ', pinew =' ,(1.d0 + pi + dpi)**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh
       IF (ABS(dpi) < min_dpi .OR. ABS(excess_demand_pi1) < toler_excessdemand_pi .OR. iter_pi == max_iterpi) THEN
          EXIT
       ENDIF
       !Update inflation
       pi            = pi + dpi
       pi_today      = pi
       cons_hh_today = cons_hh


    ENDDO
    IF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       piweight     = excessdemand_pi_low/(excessdemand_pi_low - excessdemand_pi_upp)
       pi			 = (1.d0 - piweight)*pi_low + piweight*pi_upp
       demand_pi    = (1.d0 - piweight)*demand_pi_low + piweight*demand_pi_upp
       supply_pi    = (1.d0 - piweight)*supply_pi_low + piweight*supply_pi_upp
       mu3own       = (1.d0 - piweight)*mu3own_pi_low + piweight*mu3own_pi_upp
       mu3rent      = (1.d0 - piweight)*mu3rent_pi_low + piweight*mu3rent_pi_upp
       htot_reo_new = (1.d0 - piweight)*htot_reo_new_pi_low + piweight*htot_reo_new_pi_upp
    ELSEIF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) > 100.d0) THEN
       pi			 = pi_low
       demand_pi    = demand_pi_low
       supply_pi    = supply_pi_low
       mu3own       = mu3own_pi_low
       mu3rent      = mu3rent_pi_low
       htot_reo_new = htot_reo_new_pi_low
    ELSEIF (ABS(excessdemand_pi_low) > 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       pi			 = pi_upp
       demand_pi    = demand_pi_upp
       supply_pi    = supply_pi_upp
       mu3own       = mu3own_pi_upp
       mu3rent      = mu3rent_pi_upp
       htot_reo_new = htot_reo_new_pi_upp
    ENDIF
    aggdemand = demand_pi
    aggsupply = supply_pi
    !Calculate quantities for aggregate supply and demand
    svcorigcost_t(it)  = svcorigcost
    sellcost_REO_t(it) = sellcost_REO
    sellcost_hh_t(it)  = sellcost_hh
    sellcost_t(it)     = sellcost_hh + sellcost_REO
    buycost_t(it)      = buycost
    bonds_hh_t(it)     = bonds_hh
    cons_hhnew_t(it)   = cons_hh
    cons_rent_t(it)    = cons_rent
    maint_hh_t(it)     = maint_hh
    maint_REO_t(it)    = maint_REO
    maint_t(it)        = maint_hh + maint_REO
    sh_t(it)           = sh
    pl_t(it)           = pl
    adjustcost_t(it)   = adjustcost
    aggdemand_t(it)    = aggdemand
    aggsupply_t(it)    = aggsupply
    haircutnew_t(it)   = haircut_new
    !Update inflation, lump sum transfers, and house prices
    pinew_t(it)      = pi
    IF (it == 1) THEN
       price = 1.d0 + pi
    ELSE
       price = price_t(it - 1)*(1.d0 + pi)
    ENDIF
    pricenew_t(it)   = price
    lumpsumnew_t(it) = lumpsum
    phnew_t(it)      = ph
    WRITE(*, '(a4, i4, a6, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a9, f9.6)'), ' t =', it, ', ph =' ,ph_t(it), ', phnew =' ,phnew_t(it), ', DS =' ,(hdemand - hsupply)/(0.5d0*(hdemand + hsupply)), ', Hd =' ,hactual, ', Hs =' ,hstock_prev*(1.d0 - deprec_h) + yh, ', Tr =' ,lumpsum_t(it), ', Trnew =' ,lumpsumnew_t(it)
    WRITE(*, '(a14, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6)'), '        , pi =' ,(1.d0 + pi_t(it))**4.d0 - 1.d0, ', pinew =', (1.d0 + pinew_t(it))**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh

  END SUBROUTINE sub_supplydemand_firstperiod

  SUBROUTINE sub_supplydemand_firstperiod_fast

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime

    !************ Variables for transition path ************!!!

    !Other internal variables
    REAL(8) :: excessdemand_pi

!!!Find ph that clears the housing market and also calculate lump sum transfers and excess demand in the goods market
    toler_excessdemand = 1e-4
    min_dph            = 1e-5
    max_iterph         = 100
    excessdemand_low   = 999.d0
    excessdemand_upp   = 999.d0
    ph_low             = 0.1d0
    ph_upp             = 10.d0
    DO iter_ph = 1, max_iterph
       w         = w_t(it + 1)
       dividends = dividends_t(it + 1)
       lumpsum   = lumpsum_t(it + 1)
       Jreo      = Jreo_t(:,it + 1)
       Wown	  = Wown_t(:,:,:,:,it + 1)
       Wrent     = Wrent_t(:,:,it + 1)
       Vown      = Vown_t(:,:,:,:,it + 1)
       Vrent     = Vrent_t(:,:,it + 1)
       CALL sub_backsolve_reo
       CALL sub_backsolve_ylb
       CALL sub_backsolve_valfn
       !Stock of new REO housing
       htot_reo_new  = 0.d0
       htot_reo_temp = 0.d0
       DO is = 1, ns
          DO ig = 1, ng
             ih = hfun(ig)
             CALL omp_set_num_threads(num_threads_simul)
             !$omp parallel private(id,iy)
             !$omp do
             DO im = 1, nm
                id = omp_get_thread_num() + 1
                DO iy = ny_neg + 1, ny
                   !New REOs
                   htot_reo_temp(ih,id) = htot_reo_temp(ih,id) + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                ENDDO
             ENDDO
             !$omp end do
             !$omp end parallel
          ENDDO
       ENDDO
       DO id = 1, num_threads_simul
          htot_reo_new = htot_reo_new + htot_reo_temp(:,id)
       ENDDO
       !Distribution transition from subperiod 1 to subperiod 2
       CALL sub_distribution_sp12
       !Distribution transition from subperiod 2 to subperiod 3
       CALL sub_distribution_sp23
       !Compute housing supply and demand
       CALL sub_housingclear

       IF (ABS(dph) < min_dph .OR. (demand > 0.d0 .AND. supply > 0.d0 .AND. diff < toler_excessdemand) .OR. iter_ph == max_iterph) THEN
          EXIT
       ENDIF
       !Exit in the event of nonexistence
       IF ((ph < ph_low .AND. excess_demand1 < 0.d0) .OR. (ph > ph_upp .AND. excess_demand1 > 0.d0)) THEN
          PRINT *, ' '
          PRINT *, 'Nonexistence!'
          STOP
       ENDIF
       ph = ph + dph
    ENDDO
    IF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
       phweight     = excessdemand_low/(excessdemand_low - excessdemand_upp)
       ph			 = (1.d0 - phweight)*ph_low + phweight*ph_upp
       demand       = (1.d0 - phweight)*demand_low + phweight*demand_upp
       supply       = (1.d0 - phweight)*supply_low + phweight*supply_upp
       mu3own       = (1.d0 - phweight)*mu3own_low + phweight*mu3own_upp
       mu3rent      = (1.d0 - phweight)*mu3rent_low + phweight*mu3rent_upp
       htot_reo_new = (1.d0 - phweight)*htot_reo_new_low + phweight*htot_reo_new_upp
    ELSEIF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) > 100.d0) THEN
       ph			 = ph_low
       demand       = demand_low
       supply       = supply_low
       mu3own       = mu3own_low
       mu3rent      = mu3rent_low
       htot_reo_new = htot_reo_new_low
    ELSEIF (ABS(excessdemand_low) > 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
       ph			 = ph_upp
       demand       = demand_upp
       supply       = supply_upp
       mu3own       = mu3own_upp
       mu3rent      = mu3rent_upp
       htot_reo_new = htot_reo_new_upp
    ENDIF
    IF(max_iterph == 1) ph = ph + dph
    !Compute excess demand in the goods market
    w = w_t(it)
    glabor          = agglaborsupply(w)
    CALL sub_marketclear_transition
    IF(debt_nominal .EQ. 0) THEN
         haircut_new = 1.0d0 - ((mortrev+govdebt-bankFC)*(1.0d0-ChinaDemand) + dividends + qb*StockMarket_t(it)) / (psurvive*bonds_hh)
    ELSE
        haircut_new = 1.0d0 - ((mortrev+ (1.0d0 + nom_i_t(it))*qbDebt_Steady/price_t(it)-bankFC)*(1.0d0-ChinaDemand) + dividends + qb*StockMarket_t(it)) / (psurvive*bonds_hh)
    ENDIF


    IF(debt_nominal .eq. 0) THEN

      BondDemand_t(it)  = mortcost+qbDebt_Steady
      IF(ChinaLevel .eq. 0) THEN
        ChinaPayment_t(it) = ChinaPay_Steady
        ChinaShare_t(it)  = ChinaDemand

        !Note, this is likely not correct
        ! ChinaPayment_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
        ! ChinaShare_t(it)  = ChinaDemand
      ELSE
        ChinaPayment_t(it) = ChinaDemand*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
        ChinaShare_t(it)  = ChinaDemand_steady / BondDemand_t(it)
      ENDIF
   ELSE
      BondDemand_t(it)  = mortcost+qbDebt_Steady*price_steady(it)/price_t(it)
      IF(ChinaLevel .eq. 0) THEN
        ChinaPayment_t(it) = ChinaPay_Steady*price_steady(it)/price_t(it)
        ChinaShare_t(it)  = ChinaDemand

        ! ChinaPayment_t(it) = (mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_tomorrow))/price - (mortrev+price_steady(it)*govdebt/price-bankFC))*ChinaDemand
        ! ChinaShare_t(it)  = ChinaDemand
      ELSE
        ChinaPayment_t(it) = ChinaDemand*(mortrev + (1.0d0 + nom_i_t(it))*qbDebt_Steady/price_t(it) - bankFC) - ChinaDemand_nominal_t(it)/price_t(it)
        ChinaShare_t(it)  = (ChinaDemand_nominal_t(it)/price_t(it)) / BondDemand_t(it)
      ENDIF
   ENDIF


  !  PRINT*,'Real CP=', ChinaShare_t(it - 1)*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
  !  PRINT*,'Nom CP=', ChinaShare_t(it - 1)*(mortrev + (1.0d0 + nom_i_t(it))*price_steady(it - 1)*qbDebt_Steady/price_t(it) - bankFC) - ChinaDemand_nominal_t(it)/price_t(it)
  !  PRINT*,'Real CS=', ChinaDemand_steady / BondDemand_t(it)
  !  PRINT*,'Nom CS=',(ChinaDemand_nominal_t(it + 1)/price_t(it)) / BondDemand_t(it)
  !  print*,'Real BD=',mortcost+qb*govdebt
  !  print*,'Nom BD=',mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_t(it + 1)))/price_t(it)
    aggdemand        = aggdemand + ChinaPayment_t(it)

    excessdemand_pi = (aggdemand - aggsupply)/(aggdemand + aggsupply) !govdebt - bonds_hh
    !Update inflation, lump sum transfers, and house prices
    IF(shoot_pi .eq. 1) THEN
       pinew_t(it)    = pi + max(0.0025d0,pi_steady)*10.d0*excessdemand_pi
       pricenew_t(it) = (1.0d0 + pinew_t(it))
    ELSE
      pinew_t(it)    = pi + max(0.0025d0,pi_steady)*10.d0*excessdemand_pi
      pricenew_t(it) = (1.0d0 + pinew_t(it))
      !  pricenew_t(it) = price_t(it)*(1.d0 + excessdemand_pi/2.0d0)
      !  IF (it == 1) THEN
      !     pinew_t(1) = pricenew_t(1) - 1.d0
      !  ELSE
      !     pinew_t(it) = pricenew_t(it)/pricenew_t(it - 1) - 1.d0
      !  ENDIF
    ENDIF

    lumpsumnew_t(it) = lumpsum_t(it)
    phnew_t(it)      = ph





    !Store other quantities for aggregate supply and demand
    svcorigcost_t(it)  = svcorigcost
    sellcost_REO_t(it) = sellcost_REO
    sellcost_hh_t(it)  = sellcost_hh
    sellcost_t(it)     = sellcost_hh + sellcost_REO
    buycost_t(it)      = buycost
    bonds_hh_t(it)     = bonds_hh
    cons_hhnew_t(it)   = cons_hh
    cons_rent_t(it)    = cons_rent
    maint_hh_t(it)     = maint_hh
    maint_REO_t(it)    = maint_REO
    maint_t(it)        = maint_hh + maint_REO
    sh_t(it)           = sh
    pl_t(it)           = pl
    adjustcost_t(it)   = adjustcost
    aggdemand_t(it)    = aggdemand
    aggsupply_t(it)    = aggsupply
    haircutnew_t(it)   = haircut_new
    WRITE(*, '(a4, i4, a6, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a9, f9.6)'), ' t =', it, ', ph =' ,ph_t(it), ', phnew =' ,phnew_t(it), ', DS =' ,(hdemand - hsupply)/(0.5d0*(hdemand + hsupply)), ', Hd =' ,hactual, ', Hs =' ,hstock_prev*(1.d0 - deprec_h) + yh, ', Tr =' ,lumpsum_t(it), ', Trnew =' ,lumpsumnew_t(it)
    WRITE(*,'(a9, f9.6, a9, f9.6)'),'China =  ',ChinaPayment_t(it),'ChinaR =',(mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaShare_t(it - 1)
    WRITE(*, '(a14, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6)'), '        , pi =' ,(1.d0 + pi_t(it))**4.d0 - 1.d0, ', pinew =', (1.d0 + pinew_t(it))**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh

  END SUBROUTINE sub_supplydemand_firstperiod_fast



  SUBROUTINE sub_supplydemand_slow

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime

    !************ Variables for transition path ************!!!

    !Government transfer loop variables
    REAL(8) :: toler_lumpsum, dlumpsum, taxrev, lumpsumnew, err_lumpsum, err_cons, toler_cons, toler_haircut, err_haircut
    INTEGER :: max_iterlumpsum, iter_lumpsum, iter_haircut, max_iterhaircut

    !Inflation loop variables
    REAL(8) :: toler_excessdemand_pi, min_dpi, pi_low, pi_upp, supply_pi, supply_pi_upp, supply_pi_low, demand_pi, demand_pi_low, demand_pi_upp, excess_demand_pi1, excess_demand_pi0, excessdemand_pi_upp, &
         excessdemand_pi_low, dpi, diff_pi, piweight, nom_i, cons_hh_today, haircut_low, haircut_upp
    REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: mu3own_pi_low, mu3own_pi_upp
    REAL(8), DIMENSION(:,:), ALLOCATABLE :: mu3rent_pi_low, mu3rent_pi_upp
    REAL(8), DIMENSION(:), ALLOCATABLE :: htot_reo_new_pi_low, htot_reo_new_pi_upp
    INTEGER :: max_iterpi, iter_pi, sign_excessdemand_pi1, sign_excessdemand_pi0, max_itercons

    !Other internal variables
    REAL(8) :: pi_today, w_today, lumpsum_today, dividends_today, pi_tomorrow, w_tomorrow, rscost_tomorrow, qm_tomorrow, lumpsum_tomorrow, dividends_tomorrow, StockMarket_today

    !Allocate variables
    ALLOCATE(mu3own_pi_low(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_low(ny,ns))
    ALLOCATE(mu3own_pi_upp(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_upp(ny,ns))
    ALLOCATE(htot_reo_new_pi_low(nh))
    ALLOCATE(htot_reo_new_pi_upp(nh))

!!!Find ph, pi, and transfers that satisfy the government budget constraint and clear the housing and bond/goods markets
    max_iterpi            = 100
    max_itercons          = 5
    max_iterph            = 100
    max_iterhaircut       = 10 !20
    toler_excessdemand_pi = 1e-8
    toler_excessdemand    = 1e-5
    toler_cons            = 1e-5
    toler_haircut         = 1e-5
    dlumpsum              = 0.9d0
    min_dpi               = 1e-8
    min_dph               = 1e-6
    !Initialize
    IF (fwdguide == 1) THEN
       IF (it <= T_fwdguide) THEN
          taylor_pi = 0.d0
       ELSE
          taylor_pi = 1.25d0
       ENDIF
    ENDIF
    tfp                = tfp_t(it)
    pi_today           = pi
    cons_hh_today      = cons_hh_t(it)
    lumpsum_today      = lumpsum
    qm_tomorrow        = 1.d0/(1.d0 + pi_t(MIN(T + 1,it + 2)))*1.d0/(1.d0 + rf_t(MIN(T + 1,it + 2)))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(MIN(T + 1,it + 2)))
    w_tomorrow		     = w_t(it + 1)
    lumpsum_tomorrow   = lumpsum_t(it + 1)
    dividends_tomorrow = dividends_t(it + 1)
    rscost_tomorrow    = rscost_t(it + 1)
    Jreo_tomorrow      = Jreo_t(:,it + 1)
    haircut            = haircut_t(it)

    !Inflation loop
    excessdemand_pi_low = 999.d0
    excessdemand_pi_upp = 999.d0
    pi_low              = -0.5d0
    pi_upp              = 0.5d0

    DO iter_pi = 1, max_iterpi
       !Compute today's real rate using the Taylor rule and Fisher equation
       pi          = pi_today
       gAggC       = cons_hh_today
       IF (it == 1) THEN
          price = 1.d0 + pi
       ELSE
          price = price_t(it - 1)*(1.d0 + pi)
       ENDIF
       !Guess inflation vs guess price level [uncomment next line for price level]
       IF(shoot_pi .eq. 1) THEN
         pi_tomorrow = pi_t(it + 1)
       ELSE
         pi_tomorrow = price_t(it + 1)/price - 1.d0
       ENDIF

      !  IF(inf_taylor .eq. 1) THEN
      !    nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
      !  ELSE
      !    nom_i       = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price)/(1.d0 + price_steady(it)))**taylor_pi*exp(taylorshock_t(it)) - 1.d0)
      !  ENDIF
       IF(inf_taylor .eq. 1) THEN
         IF(rho_taylor .GT. 0.01) THEN
           nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it))) - 1.d0)
         ELSE
           nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
         ENDIF
       ELSE
          nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((price/(price_steady(it)))**taylor_pi))**(1.0d0-rho_taylor)))*EXP(taylorshock_t(it))) - 1.d0)
       ENDIF

      !  nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
       ! nom_i       = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price)/(1.d0 + price_steady(it)))**taylor_pi*exp(taylorshock_t(it)) - 1.d0)
       rf          = (1.d0 + nom_i)/(1.d0 + pi_tomorrow) - 1.d0
       qb          = 1.d0/(1.d0 + rf)
       qb_annuity  = psurvive*qb
       !Compute mortgage prices
       qm0				= qm0_t(:,:,:,:,it + 1)
       Jreo			= Jreo_t(:,it + 1)
       default			= default_t(:,:,:,:,it + 1)
       psell			= psell_t(:,:,:,:,it + 1)
       locmown			= locmown_t(:,:,:,:,it + 1)
       locbown_sp1_low = locbown_sp1_low_t(:,:,:,:,it + 1)
       locbown_sp1_upp = locbown_sp1_upp_t(:,:,:,:,it + 1)
       bown_sp1_weight = bown_sp1_weight_t(:,:,:,:,it + 1)
       qm              = qm_tomorrow
       w               = w_tomorrow
       lumpsum         = lumpsum_tomorrow
       dividends       = dividends_tomorrow
       pi              = pi_tomorrow
       rscost          = rscost_tomorrow
       CALL sub_backsolve_q
       qm = 1.d0/(1.d0 + pi_tomorrow)*1.d0/(1.d0 + rf)*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_tomorrow)
       !Compute wage and dividends using firm pricing equation
       pi              = pi_today
       w               = wagedetermine(pi_tomorrow,pi,rf,tfp,aggsupply_t(it + 1),w_t(it),cons_hh_today)
       w_today         = w
       aggsupply       = tfp*agglaborsupply(w)
       glabor          = agglaborsupply(w)
       taxrev	         = labortax*w*agglaborsupply(w)
       lumpsum_today   = taxrev - (govspend + govdebt*(1.0d0-qb))
       IF(debt_nominal .eq. 0) THEN
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + (pi_today - nom_i_t(it))*qbDebt_Steady/(1.0d0 + pi_today)
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_today - nom_i_t(1))/(1.d0 + pi_today))*qbDebt_Steady
         ENDIF
       ELSE
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + ( pi_steady - nom_i_t(it) )*(price_steady(it - 1) / price)*qbDebt_Steady
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_steady - nom_i_t(1)) / price)*qbDebt_Steady
         ENDIF

       ENDIF


       IF(as_if_rigid .eq. 0) THEN
          adjustcost      = adjust_theta/2.d0*(pi_today - pi_steady)**2.d0*aggsupply
       ELSE
          adjustcost      = 0.0d0
       ENDIF
       dividends       = aggsupply - w*agglaborsupply(w) - adjustcost - firmFC
       dividends_today = dividends
       StockMarket_today = qb*StockMarket_t(it)
       !House price loop
       excessdemand_low = 999.d0
       excessdemand_upp = 999.d0
       ph_low           = 0.1d0
       ph_upp           = 10.d0
       DO iter_ph = 1, max_iterph
          w         = w_t(it + 1)
          dividends = dividends_t(it + 1)
          lumpsum   = lumpsum_t(it + 1)
          Jreo      = Jreo_t(:,it + 1)
          Wown	  = Wown_t(:,:,:,:,it + 1)
          Wrent     = Wrent_t(:,:,it + 1)
          Vown      = Vown_t(:,:,:,:,it + 1)
          Vrent     = Vrent_t(:,:,it + 1)
          CALL sub_backsolve_reo
          CALL sub_backsolve_ylb
          CALL sub_backsolve_valfn

          haircut_low = 0.0d0
          haircut_upp = 0.1d0
          !Here we figure out what the haircut has to be
          !Code updated 28 June to have the haircut be a part of the outer loop
          !so that now we read in the haircut in sub_loadtransition from the file
          !transition_haircut
          haircut_new = haircut
          do iter_haircut = 1, max_iterhaircut
            ! haircut       = (haircut_low + haircut_upp) / 2.0d0
            haircut = (haircut + haircut_new)/2.0d0
            haircut = max(-0.1d0,min(0.1d0,haircut))

            bown_today    = bown
            brent_today   = brent
            locmown_today = locmown
               !Update distribution from last period to today
             w         = w_today
             lumpsum   = lumpsum_today
             dividends = dividends_today
             mu3own    = mu3own_start
             mu3rent   = mu3rent_start
             bown      = bown_start
             brent     = brent_start
             locmown   = locmown_start
             glabor    = agglaborsupply(w)
             CALL sub_distribution_sp31
             bown    = bown_today
             brent   = brent_today
             locmown = locmown_today
             !Stock of new REO housing
             htot_reo_new  = 0.d0
             htot_reo_temp = 0.d0
             DO is = 1, ns
                DO ig = 1, ng
                   ih = hfun(ig)
                   CALL omp_set_num_threads(num_threads_simul)
                   !$omp parallel private(id,iy)
                   !$omp do
                   DO im = 1, nm
                      id = omp_get_thread_num() + 1
                      DO iy = ny_neg + 1, ny
                         !New REOs
                         htot_reo_temp(ih,id) = htot_reo_temp(ih,id) + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                      ENDDO
                   ENDDO
                   !$omp end do
                   !$omp end parallel
                ENDDO
             ENDDO
             DO id = 1, num_threads_simul
                htot_reo_new = htot_reo_new + htot_reo_temp(:,id)
             ENDDO
             !Distribution transition from subperiod 1 to subperiod 2
             CALL sub_distribution_sp12
             !Distribution transition from subperiod 2 to subperiod 3
             CALL sub_distribution_sp23
             !Calculate market clearing
             CALL sub_marketclear_transition

             ! err_b = ((govdebt + MBS)*(1.0d0-ChinaDemand) + StockMarketValue*tfp*aggL - qb_annuity*bonds_hh/qb)
             !Payout is supposed to be psurvive*bonds_hh, i.e. the
             IF(debt_nominal .EQ. 0) THEN
              !  IF(ChinaLevel == 1) THEN
              !    haircut_new = 1.0d0 - (mortrev + govdebt - bankFC - ChinaDemand_steady + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ELSE
                  haircut_new = 1.0d0 - ((mortrev+govdebt-bankFC)*(1.0d0-ChinaShare_t(it - 1)) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ENDIF
             ELSE
              !  IF(ChinaLevel == 1) THEN
              !    haircut_new = 1.0d0 - ((mortrev+price_steady(it)*govdebt/price - bankFC - price_steady(it)*ChinaDemand_steady/price) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              !  ELSE
                 haircut_new = 1.0d0 - ((mortrev+price_steady(it)*govdebt/price-bankFC)*(1.0d0-ChinaShare_t(it - 1)) + dividends_today + StockMarket_today) / (psurvive*bonds_hh)
              ENDIF



             err_haircut = haircut - haircut_new
             WRITE(*, '(a12, i4, a6, f11.7, a9, f11.7, a7, f11.7, a7, f11.7)'), '   Iter_hc =', iter_haircut, ', hc =' ,haircut, ', hcnew =' , haircut_new, ', Rev =' ,mortrev + dividends_today + govdebt + StockMarket_today, ', Lia =' ,psurvive*bonds_hh + ChinaDemand_steady
            !  print*,'Stockmarket = ',StockMarket_today
            !  print*,'dividends_today = ',dividends_today
            !  print*,'mortrev = ',mortrev
            !  print*,'govdebt = ',govdebt

             IF( abs(err_haircut) .LT. toler_haircut) THEN
               EXIT
             ENDIF
            !  IF(err_haircut .LT. 0.0d0) THEN
            !    haircut =

          ENDDO
!          haircut = 0.0d0
          !Compute housing supply and demand
          CALL sub_housingclear
          IF (ABS(dph) < min_dph .OR. (demand > 0.d0 .AND. supply > 0.d0 .AND. diff < toler_excessdemand) .OR. iter_ph == max_iterph) THEN
             EXIT
          ENDIF
          !Exit in the event of nonexistence
          IF ((ph < ph_low .AND. excess_demand1 < 0.d0) .OR. (ph > ph_upp .AND. excess_demand1 > 0.d0)) THEN
             PRINT *, ' '
             PRINT *, 'Nonexistence!'
             STOP
          ENDIF
          ph = ph + dph
       ENDDO
       IF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          phweight     = excessdemand_low/(excessdemand_low - excessdemand_upp)
          ph			 = (1.d0 - phweight)*ph_low + phweight*ph_upp
          demand       = (1.d0 - phweight)*demand_low + phweight*demand_upp
          supply       = (1.d0 - phweight)*supply_low + phweight*supply_upp
          mu3own       = (1.d0 - phweight)*mu3own_low + phweight*mu3own_upp
          mu3rent      = (1.d0 - phweight)*mu3rent_low + phweight*mu3rent_upp
          htot_reo_new = (1.d0 - phweight)*htot_reo_new_low + phweight*htot_reo_new_upp
       ELSEIF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) > 100.d0) THEN
          ph			 = ph_low
          demand       = demand_low
          supply       = supply_low
          mu3own       = mu3own_low
          mu3rent      = mu3rent_low
          htot_reo_new = htot_reo_new_low
       ELSEIF (ABS(excessdemand_low) > 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          ph			 = ph_upp
          demand       = demand_upp
          supply       = supply_upp
          mu3own       = mu3own_upp
          mu3rent      = mu3rent_upp
          htot_reo_new = htot_reo_new_upp
       ENDIF
       !Compute excess demand
      !  CALL sub_marketclear_transition
      IF(debt_nominal .eq. 0) THEN

        BondDemand_t(it)  = mortcost+qbDebt_Steady
        IF(ChinaLevel .eq. 0) THEN
          !Note, this is likely not correct
          ! ChinaPayment_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
          ! ChinaShare_t(it)  = ChinaDemand
        ELSE
          ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
          ChinaShare_t(it)  = ChinaDemand_steady / BondDemand_t(it)
        ENDIF
     ELSE
        BondDemand_t(it)  = mortcost+qbDebt_Steady*price_steady(it)/price
        IF(ChinaLevel .eq. 0) THEN
          ! ChinaPayment_t(it) = (mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_tomorrow))/price - (mortrev+price_steady(it)*govdebt/price-bankFC))*ChinaDemand
          ! ChinaShare_t(it)  = ChinaDemand
        ELSE
          ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + (1.0d0 + nom_i_t(it))*price_steady(it - 1)*qbDebt_Steady/price - bankFC) - ChinaDemand_nominal_t(it)/price
          ChinaShare_t(it)  = (ChinaDemand_nominal_t(it + 1)/price) / BondDemand_t(it)
        ENDIF
      ENDIF
        aggdemand        = aggdemand + ChinaPayment_t(it)
       demand_pi         = aggdemand !govdebt
       supply_pi         = aggsupply !bonds_hh
       excess_demand_pi1 = demand_pi - supply_pi
       IF (excess_demand_pi1 >= 0.d0) THEN
          sign_excessdemand_pi1 = 1
          IF (pi > pi_low) THEN
             pi_low              = pi
             demand_pi_low       = demand_pi
             supply_pi_low       = supply_pi
             excessdemand_pi_low = excess_demand_pi1
             mu3own_pi_low       = mu3own
             mu3rent_pi_low      = mu3rent
             htot_reo_new_pi_low = htot_reo_new
          ENDIF
       ELSE
          sign_excessdemand_pi1 = -1
          IF (pi < pi_upp) THEN
             pi_upp              = pi
             demand_pi_upp       = demand_pi
             supply_pi_upp       = supply_pi
             excessdemand_pi_upp = excess_demand_pi1
             mu3own_pi_upp       = mu3own
             mu3rent_pi_upp      = mu3rent
             htot_reo_new_pi_upp = htot_reo_new
          ENDIF
       ENDIF
       IF (iter_pi >=2 ) THEN
          IF (excess_demand_pi0 >= 0.d0) THEN
             sign_excessdemand_pi0 = 1
          ELSE
             sign_excessdemand_pi0 = -1
          ENDIF
       ENDIF
       IF (iter_pi == 1) THEN
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          !diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),MIN(0.2d0,diff_pi)*MAX(0.01d0,ABS(pi))*sign_excessdemand_pi1))		!Change min(0.2d0,... to adjust initial increment
          !				dpi     = max(0.5d0*(pi_low - pi),min(0.5d0*(pi_upp - pi),0.1d0*max(0.01d0,abs(pi))*sign_excessdemand_pi1))						!Change min(0.1d0,... to adjust initial increment
       ELSEIF (sign_excessdemand_pi1*sign_excessdemand_pi0 == -1) THEN
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),-0.5d0*dpi))
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
!          diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*ABS(pi)/ABS(dpi)
          ENDIF
       ELSE
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),dpi))
          diff_pi = 10.d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
!          diff_pi = ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*ABS(pi)) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))/ABS(dpi)
          ENDIF
       ENDIF
       excess_demand_pi0 = excess_demand_pi1
       WRITE(*, '(a12, i4, a6, f11.7, a9, f11.7, a6, f11.7, a6, f11.7, a6, f12.8, a6, f12.8)'), '   Iter_pi =', iter_pi, ', pi =' ,(1.d0 + pi)**4.d0 - 1.d0, ', pinew =' ,(1.d0 + pi + dpi)**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh
       IF (ABS(dpi) < min_dpi .OR. ABS(excess_demand_pi1) < toler_excessdemand_pi .OR. iter_pi == max_iterpi) THEN
          EXIT
       ENDIF
       !Update inflation
       pi            = pi + dpi
       pi_today      = pi
       cons_hh_today = cons_hh
      !  IF(ChinaLevel .eq. 0) THEN
      !    ChinaDemand_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
      !    BondDemand_t(it)  = mortcost+qb*govdebt
      !  ENDIF

    ENDDO


    IF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       piweight     = excessdemand_pi_low/(excessdemand_pi_low - excessdemand_pi_upp)
       pi			 = (1.d0 - piweight)*pi_low + piweight*pi_upp
       demand_pi    = (1.d0 - piweight)*demand_pi_low + piweight*demand_pi_upp
       supply_pi    = (1.d0 - piweight)*supply_pi_low + piweight*supply_pi_upp
       mu3own       = (1.d0 - piweight)*mu3own_pi_low + piweight*mu3own_pi_upp
       mu3rent      = (1.d0 - piweight)*mu3rent_pi_low + piweight*mu3rent_pi_upp
       htot_reo_new = (1.d0 - piweight)*htot_reo_new_pi_low + piweight*htot_reo_new_pi_upp
    ELSEIF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) > 100.d0) THEN
       pi			 = pi_low
       demand_pi    = demand_pi_low
       supply_pi    = supply_pi_low
       mu3own       = mu3own_pi_low
       mu3rent      = mu3rent_pi_low
       htot_reo_new = htot_reo_new_pi_low
    ELSEIF (ABS(excessdemand_pi_low) > 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       pi			 = pi_upp
       demand_pi    = demand_pi_upp
       supply_pi    = supply_pi_upp
       mu3own       = mu3own_pi_upp
       mu3rent      = mu3rent_pi_upp
       htot_reo_new = htot_reo_new_pi_upp
    ENDIF
    aggdemand = demand_pi
    aggsupply = supply_pi
    !Calculate quantities for aggregate supply and demand
    svcorigcost_t(it)  = svcorigcost
    sellcost_REO_t(it) = sellcost_REO
    sellcost_hh_t(it)  = sellcost_hh
    sellcost_t(it)     = sellcost_hh + sellcost_REO
    buycost_t(it)      = buycost
    bonds_hh_t(it)     = bonds_hh
    cons_hhnew_t(it)   = cons_hh
    cons_rent_t(it)    = cons_rent
    maint_hh_t(it)     = maint_hh
    maint_REO_t(it)    = maint_REO
    maint_t(it)        = maint_hh + maint_REO
    sh_t(it)           = sh
    pl_t(it)           = pl
    adjustcost_t(it)   = adjustcost
    aggdemand_t(it)    = aggdemand
    aggsupply_t(it)    = aggsupply
    haircutnew_t(it)   = haircut_new
    !Update inflation, lump sum transfers, and house prices
    pinew_t(it)      = pi
    IF (it == 1) THEN
       price = 1.d0 + pi
    ELSE
       price = price_t(it - 1)*(1.d0 + pi)
    ENDIF
    pricenew_t(it)   = price
    lumpsumnew_t(it) = lumpsum
    phnew_t(it)      = ph
    WRITE(*, '(a4, i4, a6, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a9, f9.6)'), ' t =', it, ', ph =' ,ph_t(it), ', phnew =' ,phnew_t(it), ', DS =' ,(hdemand - hsupply)/(0.5d0*(hdemand + hsupply)), ', Hd =' ,hactual, ', Hs =' ,hstock_prev*(1.d0 - deprec_h) + yh, ', Tr =' ,lumpsum_t(it), ', Trnew =' ,lumpsumnew_t(it)
    WRITE(*, '(a14, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6)'), '        , pi =' ,(1.d0 + pi_t(it))**4.d0 - 1.d0, ', pinew =', (1.d0 + pinew_t(it))**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh

  END SUBROUTINE sub_supplydemand_slow



  SUBROUTINE sub_supplydemand

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime

    !************ Variables for transition path ************!!!

    !Government transfer loop variables
    REAL(8) :: toler_lumpsum, dlumpsum, taxrev, lumpsumnew, err_lumpsum, err_cons, toler_cons
    INTEGER :: max_iterlumpsum, iter_lumpsum

    !Inflation loop variables
    REAL(8) :: toler_excessdemand_pi, min_dpi, pi_low, pi_upp, supply_pi, supply_pi_upp, supply_pi_low, demand_pi, demand_pi_low, demand_pi_upp, excess_demand_pi1, excess_demand_pi0, excessdemand_pi_upp, &
         excessdemand_pi_low, dpi, diff_pi, piweight, nom_i, cons_hh_today
    REAL(8), DIMENSION(:,:,:,:), ALLOCATABLE :: mu3own_pi_low, mu3own_pi_upp
    REAL(8), DIMENSION(:,:), ALLOCATABLE :: mu3rent_pi_low, mu3rent_pi_upp
    REAL(8), DIMENSION(:), ALLOCATABLE :: htot_reo_new_pi_low, htot_reo_new_pi_upp
    INTEGER :: max_iterpi, iter_pi, sign_excessdemand_pi1, sign_excessdemand_pi0, max_itercons

    !Other internal variables
    REAL(8) :: pi_today, w_today, lumpsum_today, dividends_today, pi_tomorrow, w_tomorrow, rscost_tomorrow, qm_tomorrow, lumpsum_tomorrow, dividends_tomorrow

    !Allocate variables
    ALLOCATE(mu3own_pi_low(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_low(ny,ns))
    ALLOCATE(mu3own_pi_upp(ny,nm,ng,ns))
    ALLOCATE(mu3rent_pi_upp(ny,ns))
    ALLOCATE(htot_reo_new_pi_low(nh))
    ALLOCATE(htot_reo_new_pi_upp(nh))

!!!Find ph, pi, and transfers that satisfy the government budget constraint and clear the housing and bond/goods markets
    max_iterpi            = 100
    max_itercons          = 5
    max_iterph            = 100
    toler_excessdemand_pi = 1e-8
    toler_excessdemand    = 1e-5
    toler_cons            = 1e-5
    dlumpsum              = 0.9d0
    min_dpi               = 1e-8
    min_dph               = 1e-6
    !Initialize
    IF (fwdguide == 1) THEN
       IF (it <= T_fwdguide) THEN
          taylor_pi = 0.d0
       ELSE
          taylor_pi = 1.25d0
       ENDIF
    ENDIF
    haircut            = haircut_t(it)
    tfp                = tfp_t(it)
    pi_today           = pi
    cons_hh_today      = cons_hh_t(it)
    lumpsum_today      = lumpsum
    qm_tomorrow        = 1.d0/(1.d0 + pi_t(MIN(T + 1,it + 2)))*1.d0/(1.d0 + rf_t(MIN(T + 1,it + 2)))*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_t(MIN(T + 1,it + 2)))
    w_tomorrow		     = w_t(it + 1)
    lumpsum_tomorrow   = lumpsum_t(it + 1)
    dividends_tomorrow = dividends_t(it + 1)
    rscost_tomorrow    = rscost_t(it + 1)
    Jreo_tomorrow      = Jreo_t(:,it + 1)
    !Inflation loop
    excessdemand_pi_low = 999.d0
    excessdemand_pi_upp = 999.d0
    pi_low              = -0.5d0
    pi_upp              = 0.5d0
    DO iter_pi = 1, max_iterpi
       !Compute today's real rate using the Taylor rule and Fisher equation
       pi          = pi_today
       gAggC       = cons_hh_today
       IF (it == 1) THEN
          price = 1.d0 + pi
       ELSE
          price = price_t(it - 1)*(1.d0 + pi)
       ENDIF


       !Use price level or inflation to determine tomorrow's inflation
      !  pi_tomorrow = price_t(it + 1)/price - 1.d0
      !  pi_tomorrow = pi_t(it + 1)
       IF(shoot_pi .eq. 1) THEN
         pi_tomorrow = pi_t(it + 1)
       ELSE
         pi_tomorrow = price_t(it + 1)/price - 1.d0
       ENDIF
      !  IF(inf_taylor .eq. 1) THEN
      !    nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
      !  ELSE
      !    nom_i       = max(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + price)/(1.d0 + price_steady(it)))**taylor_pi*exp(taylorshock_t(it)) - 1.d0)
      !  ENDIF
      IF(inf_taylor .eq. 1) THEN
        ! nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it))) - 1.d0)
        IF(rho_taylor .GT. 0.01) THEN
          nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi))**(1.0d0-rho_taylor))*EXP(taylorshock_t(it))) - 1.d0)
        ELSE
          nom_i       = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*((1.d0 + pi)/(1.d0 + pi_steady))**taylor_pi*EXP(taylorshock_t(it)) - 1.d0)
        ENDIF

      ELSE
         nom_i = MAX(0.d0,(1.d0 + rf_steady)*(1.d0 + pi_steady)*( (((1.0d0+nom_i_t(it))/(1.0d0 + nom_i_steady))**rho_taylor )*(((((price/(price_steady(it)))**taylor_pi))**(1.0d0-rho_taylor)))*EXP(taylorshock_t(it))) - 1.d0)
      ENDIF

       rf          = (1.d0 + nom_i)/(1.d0 + pi_tomorrow) - 1.d0
       qb          = 1.d0/(1.d0 + rf)
       qb_annuity  = psurvive*qb
       !Compute mortgage prices
       qm0				= qm0_t(:,:,:,:,it + 1)
       Jreo			= Jreo_t(:,it + 1)
       default			= default_t(:,:,:,:,it + 1)
       psell			= psell_t(:,:,:,:,it + 1)
       locmown			= locmown_t(:,:,:,:,it + 1)
       locbown_sp1_low = locbown_sp1_low_t(:,:,:,:,it + 1)
       locbown_sp1_upp = locbown_sp1_upp_t(:,:,:,:,it + 1)
       bown_sp1_weight = bown_sp1_weight_t(:,:,:,:,it + 1)
       qm              = qm_tomorrow
       w               = w_tomorrow
       lumpsum         = lumpsum_tomorrow
       dividends       = dividends_tomorrow
       pi              = pi_tomorrow
       rscost          = rscost_tomorrow
       CALL sub_backsolve_q
       qm = 1.d0/(1.d0 + pi_tomorrow)*1.d0/(1.d0 + rf)*psurvive*(1.d0 - deprec_h_survive)/(1.d0 + rscost_tomorrow)
       !Compute wage and dividends using firm pricing equation
       pi              = pi_today
       w               = wagedetermine(pi_tomorrow,pi,rf,tfp,aggsupply_t(it + 1),w_t(it),cons_hh_today)
       w_today         = w
       glabor          = agglaborsupply(w)
       aggsupply       = tfp*agglaborsupply(w)
       taxrev	         = labortax*w*agglaborsupply(w)

       IF(debt_nominal .eq. 0) THEN
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + (pi_today - nom_i_t(it))*qbDebt_Steady/(1.0d0 + pi_today)
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_today - nom_i_t(1))/(1.d0 + pi_today))*qbDebt_Steady
         ENDIF
       ELSE
         IF(it .GT. 1) THEN
           lumpsum_today = labortax*w_today*glabor - govspend + ( pi_steady - nom_i_t(it) )*(price_steady(it - 1) / price)*qbDebt_Steady
         ELSE
           lumpsum_today = labortax*w_today*glabor - govspend + ((pi_steady - nom_i_t(1)) / price)*qbDebt_Steady
         ENDIF

       ENDIF

       IF(as_if_rigid .EQ. 0) THEN
          adjustcost      = adjust_theta/2.d0*(pi_today - pi_steady)**2.d0*aggsupply
       ELSE
          adjustcost      = 0.0d0
       ENDIF
       dividends       = aggsupply - w*agglaborsupply(w) - adjustcost - firmFC
       dividends_today = dividends
       !House price loop
       excessdemand_low = 999.d0
       excessdemand_upp = 999.d0
       ph_low           = 0.1d0
       ph_upp           = 10.d0
       DO iter_ph = 1, max_iterph
          w         = w_t(it + 1)
          dividends = dividends_t(it + 1)
          lumpsum   = lumpsum_t(it + 1)
          Jreo      = Jreo_t(:,it + 1)
          Wown	  = Wown_t(:,:,:,:,it + 1)
          Wrent     = Wrent_t(:,:,it + 1)
          Vown      = Vown_t(:,:,:,:,it + 1)
          Vrent     = Vrent_t(:,:,it + 1)
          CALL sub_backsolve_reo
          CALL sub_backsolve_ylb
          CALL sub_backsolve_valfn
          bown_today    = bown
          brent_today   = brent
          locmown_today = locmown
             !Update distribution from last period to today
          w         = w_today
          glabor    = agglaborsupply(w)
           lumpsum   = lumpsum_today
           dividends = dividends_today
           mu3own    = mu3own_start
           mu3rent   = mu3rent_start
           bown      = bown_start
           brent     = brent_start
           locmown   = locmown_start
           CALL sub_distribution_sp31
           bown    = bown_today
           brent   = brent_today
           locmown = locmown_today
           !Stock of new REO housing
           htot_reo_new  = 0.d0
           htot_reo_temp = 0.d0
           DO is = 1, ns
              DO ig = 1, ng
                 ih = hfun(ig)
                 CALL omp_set_num_threads(num_threads_simul)
                 !$omp parallel private(id,iy)
                 !$omp do
                 DO im = 1, nm
                    id = omp_get_thread_num() + 1
                    DO iy = ny_neg + 1, ny
                       !New REOs
                       htot_reo_temp(ih,id) = htot_reo_temp(ih,id) + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                    ENDDO
                 ENDDO
                 !$omp end do
                 !$omp end parallel
              ENDDO
           ENDDO
           DO id = 1, num_threads_simul
              htot_reo_new = htot_reo_new + htot_reo_temp(:,id)
           ENDDO
           !Distribution transition from subperiod 1 to subperiod 2
           CALL sub_distribution_sp12
           !Distribution transition from subperiod 2 to subperiod 3
           CALL sub_distribution_sp23
           !Calculate market clearing
           CALL sub_marketclear_transition
          !Compute housing supply and demand
          CALL sub_housingclear
          IF (ABS(dph) < min_dph .OR. (demand > 0.d0 .AND. supply > 0.d0 .AND. diff < toler_excessdemand) .OR. iter_ph == max_iterph) THEN
             EXIT
          ENDIF
          !Exit in the event of nonexistence
          IF ((ph < ph_low .AND. excess_demand1 < 0.d0) .OR. (ph > ph_upp .AND. excess_demand1 > 0.d0)) THEN
             PRINT *, ' '
             PRINT *, 'Nonexistence!'
             STOP
          ENDIF
          ph = ph + dph
       ENDDO
       IF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          phweight     = excessdemand_low/(excessdemand_low - excessdemand_upp)
          ph			 = (1.d0 - phweight)*ph_low + phweight*ph_upp
          demand       = (1.d0 - phweight)*demand_low + phweight*demand_upp
          supply       = (1.d0 - phweight)*supply_low + phweight*supply_upp
          mu3own       = (1.d0 - phweight)*mu3own_low + phweight*mu3own_upp
          mu3rent      = (1.d0 - phweight)*mu3rent_low + phweight*mu3rent_upp
          htot_reo_new = (1.d0 - phweight)*htot_reo_new_low + phweight*htot_reo_new_upp
       ELSEIF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) > 100.d0) THEN
          ph			 = ph_low
          demand       = demand_low
          supply       = supply_low
          mu3own       = mu3own_low
          mu3rent      = mu3rent_low
          htot_reo_new = htot_reo_new_low
       ELSEIF (ABS(excessdemand_low) > 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
          ph			 = ph_upp
          demand       = demand_upp
          supply       = supply_upp
          mu3own       = mu3own_upp
          mu3rent      = mu3rent_upp
          htot_reo_new = htot_reo_new_upp
       ENDIF
       !Compute excess demand
      !  CALL sub_marketclear_transition
      IF(debt_nominal .eq. 0) THEN

        BondDemand_t(it)  = mortcost+qbDebt_Steady
        IF(ChinaLevel .eq. 0) THEN
          !Note, this is likely not correct
          ! ChinaPayment_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
          ! ChinaShare_t(it)  = ChinaDemand
        ELSE
          ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
          ChinaShare_t(it)  = ChinaDemand_steady / BondDemand_t(it)
        ENDIF
     ELSE
        BondDemand_t(it)  = mortcost+qbDebt_Steady*price_steady(it)/price
        IF(ChinaLevel .eq. 0) THEN
          ! ChinaPayment_t(it) = (mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_tomorrow))/price - (mortrev+price_steady(it)*govdebt/price-bankFC))*ChinaDemand
          ! ChinaShare_t(it)  = ChinaDemand
        ELSE
          ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + (1.0d0 + nom_i_t(it))*price_steady(it - 1)*qbDebt_Steady/price - bankFC) - ChinaDemand_nominal_t(it)/price
          ChinaShare_t(it)  = (ChinaDemand_nominal_t(it + 1)/price) / BondDemand_t(it)
        ENDIF
      ENDIF
       aggdemand         = aggdemand + ChinaPayment_t(it)
       demand_pi         = aggdemand !govdebt
       supply_pi         = aggsupply !bonds_hh
       excess_demand_pi1 = demand_pi - supply_pi
       IF (excess_demand_pi1 >= 0.d0) THEN
          sign_excessdemand_pi1 = 1
          IF (pi > pi_low) THEN
             pi_low              = pi
             demand_pi_low       = demand_pi
             supply_pi_low       = supply_pi
             excessdemand_pi_low = excess_demand_pi1
             mu3own_pi_low       = mu3own
             mu3rent_pi_low      = mu3rent
             htot_reo_new_pi_low = htot_reo_new
          ENDIF
       ELSE
          sign_excessdemand_pi1 = -1
          IF (pi < pi_upp) THEN
             pi_upp              = pi
             demand_pi_upp       = demand_pi
             supply_pi_upp       = supply_pi
             excessdemand_pi_upp = excess_demand_pi1
             mu3own_pi_upp       = mu3own
             mu3rent_pi_upp      = mu3rent
             htot_reo_new_pi_upp = htot_reo_new
          ENDIF
       ENDIF
       IF (iter_pi >=2 ) THEN
          IF (excess_demand_pi0 >= 0.d0) THEN
             sign_excessdemand_pi0 = 1
          ELSE
             sign_excessdemand_pi0 = -1
          ENDIF
       ENDIF
       IF (iter_pi == 1) THEN
          diff_pi = 10.0d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),MIN(0.2d0,diff_pi)*MAX(0.01d0,ABS(pi))*sign_excessdemand_pi1))		!Change min(0.2d0,... to adjust initial increment
          !				dpi     = max(0.5d0*(pi_low - pi),min(0.5d0*(pi_upp - pi),0.1d0*max(0.01d0,abs(pi))*sign_excessdemand_pi1))						!Change min(0.1d0,... to adjust initial increment
       ELSEIF (sign_excessdemand_pi1*sign_excessdemand_pi0 == -1) THEN
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),-0.5d0*dpi))
          diff_pi = 10.0d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*ABS(pi)/ABS(dpi)
          ENDIF
       ELSE
          dpi     = MAX(0.5d0*(pi_low - pi),MIN(0.5d0*(pi_upp - pi),dpi))
          diff_pi = 10.0d0*ABS(excess_demand_pi1)/(0.5*demand_pi + 0.5*supply_pi)															!The 10.d0 in front is simply to speed things up
          IF (ABS(dpi) > MIN(1.d0,diff_pi)*ABS(pi)) THEN
             dpi = dpi*MIN(1.d0,diff_pi)*MAX(0.01d0,ABS(pi))/ABS(dpi)
          ENDIF
       ENDIF
       excess_demand_pi0 = excess_demand_pi1
       WRITE(*, '(a12, i4, a6, f11.7, a9, f11.7, a6, f11.7, a6, f11.7, a6, f12.8, a6, f12.8)'), '   Iter_pi =', iter_pi, ', pi =' ,(1.d0 + pi)**4.d0 - 1.d0, ', pinew =' ,(1.d0 + pi + dpi)**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh
       IF (ABS(dpi) < min_dpi .OR. ABS(excess_demand_pi1) < toler_excessdemand_pi .OR. iter_pi == max_iterpi) THEN
          EXIT
       ENDIF
       !Update inflation
       pi            = pi + dpi
       pi_today      = pi
       cons_hh_today = cons_hh



    ENDDO
    IF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       piweight     = excessdemand_pi_low/(excessdemand_pi_low - excessdemand_pi_upp)
       pi			 = (1.d0 - piweight)*pi_low + piweight*pi_upp
       demand_pi    = (1.d0 - piweight)*demand_pi_low + piweight*demand_pi_upp
       supply_pi    = (1.d0 - piweight)*supply_pi_low + piweight*supply_pi_upp
       mu3own       = (1.d0 - piweight)*mu3own_pi_low + piweight*mu3own_pi_upp
       mu3rent      = (1.d0 - piweight)*mu3rent_pi_low + piweight*mu3rent_pi_upp
       htot_reo_new = (1.d0 - piweight)*htot_reo_new_pi_low + piweight*htot_reo_new_pi_upp
    ELSEIF (ABS(excessdemand_pi_low) < 100.d0 .AND. ABS(excessdemand_pi_upp) > 100.d0) THEN
       pi			 = pi_low
       demand_pi    = demand_pi_low
       supply_pi    = supply_pi_low
       mu3own       = mu3own_pi_low
       mu3rent      = mu3rent_pi_low
       htot_reo_new = htot_reo_new_pi_low
    ELSEIF (ABS(excessdemand_pi_low) > 100.d0 .AND. ABS(excessdemand_pi_upp) < 100.d0) THEN
       pi			 = pi_upp
       demand_pi    = demand_pi_upp
       supply_pi    = supply_pi_upp
       mu3own       = mu3own_pi_upp
       mu3rent      = mu3rent_pi_upp
       htot_reo_new = htot_reo_new_pi_upp
    ENDIF
    aggdemand = demand_pi
    aggsupply = supply_pi
    !Calculate quantities for aggregate supply and demand
    svcorigcost_t(it)  = svcorigcost
    sellcost_REO_t(it) = sellcost_REO
    sellcost_hh_t(it)  = sellcost_hh
    sellcost_t(it)     = sellcost_hh + sellcost_REO
    buycost_t(it)      = buycost
    bonds_hh_t(it)     = bonds_hh
    cons_hhnew_t(it)   = cons_hh
    cons_rent_t(it)    = cons_rent
    maint_hh_t(it)     = maint_hh
    maint_REO_t(it)    = maint_REO
    maint_t(it)        = maint_hh + maint_REO
    sh_t(it)           = sh
    pl_t(it)           = pl
    adjustcost_t(it)   = adjustcost
    aggdemand_t(it)    = aggdemand
    aggsupply_t(it)    = aggsupply
    !Update inflation, lump sum transfers, and house prices
    pinew_t(it)      = pi
    pricenew_t(it)   = price
    lumpsumnew_t(it) = lumpsum
    phnew_t(it)      = ph
    haircutnew_t(it) = 0.0d0
    WRITE(*, '(a4, i4, a6, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a9, f9.6)'), ' t =', it, ', ph =' ,ph_t(it), ', phnew =' ,phnew_t(it), ', DS =' ,(hdemand - hsupply)/(0.5d0*(hdemand + hsupply)), ', Hd =' ,hactual, ', Hs =' ,hstock_prev*(1.d0 - deprec_h) + yh, ', Tr =' ,lumpsum_t(it), ', Trnew =' ,lumpsumnew_t(it)
    WRITE(*, '(a14, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6)'), '        , pi =' ,(1.d0 + pi_t(it))**4.d0 - 1.d0, ', pinew =', (1.d0 + pinew_t(it))**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh

  END SUBROUTINE sub_supplydemand

  SUBROUTINE sub_supplydemand_fast

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices for the individual state
    INTEGER :: is, ie, ig, ih, iy, im, imprime

    !************ Variables for transition path ************!!!

    !Other internal variables
    REAL(8) :: excessdemand_pi

!!!Find ph that clears the housing market and also calculate lump sum transfers and excess demand in the goods market
    toler_excessdemand = 1e-4
    min_dph            = 1e-5
    max_iterph         = 100
    excessdemand_low   = 999.d0
    excessdemand_upp   = 999.d0
    ph_low             = 0.1d0
    ph_upp             = 10.d0
    DO iter_ph = 1, max_iterph
       w         = w_t(it + 1)
       dividends = dividends_t(it + 1)
       lumpsum   = lumpsum_t(it + 1)
       Jreo      = Jreo_t(:,it + 1)
       Wown	  = Wown_t(:,:,:,:,it + 1)
       Wrent     = Wrent_t(:,:,it + 1)
       Vown      = Vown_t(:,:,:,:,it + 1)
       Vrent     = Vrent_t(:,:,it + 1)
       CALL sub_backsolve_reo
       CALL sub_backsolve_ylb
       CALL sub_backsolve_valfn
       !Stock of new REO housing
       htot_reo_new  = 0.d0
       htot_reo_temp = 0.d0
       DO is = 1, ns
          DO ig = 1, ng
             ih = hfun(ig)
             CALL omp_set_num_threads(num_threads_simul)
             !$omp parallel private(id,iy)
             !$omp do
             DO im = 1, nm
                id = omp_get_thread_num() + 1
                DO iy = ny_neg + 1, ny
                   !New REOs
                   htot_reo_temp(ih,id) = htot_reo_temp(ih,id) + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                ENDDO
             ENDDO
             !$omp end do
             !$omp end parallel
          ENDDO
       ENDDO
       DO id = 1, num_threads_simul
          htot_reo_new = htot_reo_new + htot_reo_temp(:,id)
       ENDDO
       !Distribution transition from subperiod 1 to subperiod 2
       CALL sub_distribution_sp12
       !Distribution transition from subperiod 2 to subperiod 3
       CALL sub_distribution_sp23
       !Compute housing supply and demand
       CALL sub_housingclear

       IF (ABS(dph) < min_dph .OR. (demand > 0.d0 .AND. supply > 0.d0 .AND. diff < toler_excessdemand) .OR. iter_ph == max_iterph) THEN
          EXIT
       ENDIF
       !Exit in the event of nonexistence
       IF ((ph < ph_low .AND. excess_demand1 < 0.d0) .OR. (ph > ph_upp .AND. excess_demand1 > 0.d0)) THEN
          PRINT *, ' '
          PRINT *, 'Nonexistence!'
          STOP
       ENDIF
       ph = ph + dph
    ENDDO
    IF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
       phweight     = excessdemand_low/(excessdemand_low - excessdemand_upp)
       ph			 = (1.d0 - phweight)*ph_low + phweight*ph_upp
       demand       = (1.d0 - phweight)*demand_low + phweight*demand_upp
       supply       = (1.d0 - phweight)*supply_low + phweight*supply_upp
       mu3own       = (1.d0 - phweight)*mu3own_low + phweight*mu3own_upp
       mu3rent      = (1.d0 - phweight)*mu3rent_low + phweight*mu3rent_upp
       htot_reo_new = (1.d0 - phweight)*htot_reo_new_low + phweight*htot_reo_new_upp
    ELSEIF (ABS(excessdemand_low) < 100.d0 .AND. ABS(excessdemand_upp) > 100.d0) THEN
       ph			 = ph_low
       demand       = demand_low
       supply       = supply_low
       mu3own       = mu3own_low
       mu3rent      = mu3rent_low
       htot_reo_new = htot_reo_new_low
    ELSEIF (ABS(excessdemand_low) > 100.d0 .AND. ABS(excessdemand_upp) < 100.d0) THEN
       ph			 = ph_upp
       demand       = demand_upp
       supply       = supply_upp
       mu3own       = mu3own_upp
       mu3rent      = mu3rent_upp
       htot_reo_new = htot_reo_new_upp
    ENDIF
    IF(max_iterph == 1) ph = ph + dph

    !Compute excess demand in the goods market
    w = w_t(it)
    glabor          = agglaborsupply(w)
    CALL sub_marketclear_transition
    IF(debt_nominal .eq. 0) THEN

      BondDemand_t(it)  = mortcost+qbDebt_Steady
      IF(ChinaLevel .eq. 0) THEN
        !Note, this is likely not correct
        ! ChinaPayment_t(it) = (mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaDemand
        ! ChinaShare_t(it)  = ChinaDemand
        ChinaPayment_t(it) = ChinaPay_Steady
        ChinaShare_t(it)  = ChinaDemand

      ELSE
        ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
        ChinaShare_t(it)  = ChinaDemand_steady / BondDemand_t(it)
      ENDIF
   ELSE
      BondDemand_t(it)  = mortcost+qbDebt_Steady*price_steady(it)/price_t(it)
      IF(ChinaLevel .eq. 0) THEN
        ! ChinaPayment_t(it) = (mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_tomorrow))/price - (mortrev+price_steady(it)*govdebt/price-bankFC))*ChinaDemand
        ! ChinaShare_t(it)  = ChinaDemand
        ChinaPayment_t(it) = ChinaPay_Steady*price_steady(it)/price_t(it)
        ChinaShare_t(it)  = ChinaDemand

      ELSE
        ChinaPayment_t(it) = ChinaShare_t(it - 1)*(mortrev + (1.0d0 + nom_i_t(it))*price_steady(it - 1)*qbDebt_Steady/price_t(it) - bankFC) - ChinaDemand_nominal_t(it)/price_t(it)
        ChinaShare_t(it)  = (ChinaDemand_nominal_t(it)/price_t(it)) / BondDemand_t(it)
      ENDIF
   ENDIF
  !  PRINT*,'Real CP=', ChinaShare_t(it - 1)*(mortrev + qbDebt_Steady*(1.d0 + nom_i_t(it))/(1.d0 + pi_t(it)) - bankFC) - ChinaDemand_steady
  !  PRINT*,'Nom CP=', ChinaShare_t(it - 1)*(mortrev + (1.0d0 + nom_i_t(it))*price_steady(it - 1)*qbDebt_Steady/price_t(it) - bankFC) - ChinaDemand_nominal_t(it)/price_t(it)
  !  PRINT*,'Real CS=', ChinaDemand_steady / BondDemand_t(it)
  !  PRINT*,'Nom CS=',(ChinaDemand_nominal_t(it + 1)/price_t(it)) / BondDemand_t(it)
  !  print*,'Real BD=',mortcost+qb*govdebt
  !  print*,'Nom BD=',mortcost+qb*(nominal_debt_t(it + 1)/(1.0d0 + pi_t(it + 1)))/price_t(it)
    aggdemand        = aggdemand + ChinaPayment_t(it)

    excessdemand_pi = (aggdemand - aggsupply)/(aggdemand + aggsupply) !govdebt - bonds_hh
    !Update inflation, lump sum transfers, and house prices
    IF(shoot_pi .eq. 1) THEN
       pinew_t(it)    = pi + max(0.0025d0,pi_steady)*2.d0*excessdemand_pi
       pricenew_t(it) = price_t(it - 1) * (1.0d0 + pinew_t(it))
    ELSE
      pinew_t(it)    = pi + max(0.0025d0,pi_steady)*2.d0*excessdemand_pi
      pricenew_t(it) = price_t(it - 1) * (1.0d0 + pinew_t(it))

      !  pricenew_t(it) = price_t(it)*(1.d0 + excessdemand_pi/5.0d0)
      !  IF (it == 1) THEN
      !     pinew_t(1) = pricenew_t(1) - 1.d0
      !  ELSE
      !     pinew_t(it) = pricenew_t(it)/pricenew_t(it - 1) - 1.d0
      !  ENDIF
    ENDIF
    ! gAggC            = cons_hh_t(it)
!     IF(debt_nominal .eq. 0) THEN
!       lumpsumnew_t(it) = labortax*w*agglaborsupply(w) - (govspend + govdebt*(1.0d0-qb))
!     ELSE
!       ! lumpsumnew_t(it) = labortax*w*agglaborsupply(w) - (govspend + govdebt*(1.0d0-qb)*price_steady(it)/price_t(it))
! !      lumpsumnew_t(it) = labortax*w*agglaborsupply(w) - (govspend + (nominal_debt_t(it) - nominal_debt_t(it + 1)/(1.0d0 + nom_i_t(it + 1)))/price_t(it))
!       lumpsumnew_t(it) = labortax*w*agglaborsupply(w) - (govspend + (nominal_debt_t(it) - qb*nominal_debt_t(it + 1)/(1.0d0 + pi_t(it + 1)))/price_t(it))
!
!       !lumpsumnew_t(it) = labortax*w*agglaborsupply(w) - (govspend + govdebt*(1.0d0-qb*(1.0d0+pi_steady))*price_steady(it)/price_t(it))
!       !lumpsumnew_t(it)   = labortax*w*agglaborsupply(w) - (govspend + govdebt*(pi_steady - nom_i_t(it + 1))*price_steady(it)/price)
!       !lumpsumnew_t(it)   = labortax*w*agglaborsupply(w) - (govspend + govdebt*(pi_steady - nom_i_t(it + 1))*price_steady(it)/price/(1.0d0 + nom_i_t(it + 1)))
!     ENDIF
    lumpsumnew_t(it) = lumpsum_t(it)
    phnew_t(it)      = ph





    !Store other quantities for aggregate supply and demand
    svcorigcost_t(it)  = svcorigcost
    sellcost_REO_t(it) = sellcost_REO
    sellcost_hh_t(it)  = sellcost_hh
    sellcost_t(it)     = sellcost_hh + sellcost_REO
    buycost_t(it)      = buycost
    bonds_hh_t(it)     = bonds_hh
    cons_hhnew_t(it)   = cons_hh
    cons_rent_t(it)    = cons_rent
    maint_hh_t(it)     = maint_hh
    maint_REO_t(it)    = maint_REO
    maint_t(it)        = maint_hh + maint_REO
    sh_t(it)           = sh
    pl_t(it)           = pl
    adjustcost_t(it)   = adjustcost
    aggdemand_t(it)    = aggdemand
    aggsupply_t(it)    = aggsupply
    haircutnew_t(it)   = 0.0d0
    WRITE(*, '(a4, i4, a6, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a9, f9.6)'), ' t =', it, ', ph =' ,ph_t(it), ', phnew =' ,phnew_t(it), ', DS =' ,(hdemand - hsupply)/(0.5d0*(hdemand + hsupply)), ', Hd =' ,hactual, ', Hs =' ,hstock_prev*(1.d0 - deprec_h) + yh, ', Tr =' ,lumpsum_t(it), ', Trnew =' ,lumpsumnew_t(it)
    WRITE(*,'(a9, f9.6, a9, f9.6)'),'China =  ',ChinaPayment_t(it),'ChinaR =',(mortcost + qb*govdebt - (mortrev + govdebt - bankFC))*ChinaShare_t(it - 1)
    WRITE(*, '(a14, f9.6, a9, f9.6, a6, f9.6, a6, f9.6, a6, f9.6, a6, f9.6)'), '        , pi =' ,(1.d0 + pi_t(it))**4.d0 - 1.d0, ', pinew =', (1.d0 + pinew_t(it))**4.d0 - 1.d0, ', AD =' ,aggdemand, ', AS =' ,aggsupply, ', BD =', BondDemand_t(it)*(1.0d0 - ChinaShare_t(it)) + qb*Stockmarket_t(it), ', BS =', qb_annuity*bonds_hh

  END SUBROUTINE sub_supplydemand_fast

  SUBROUTINE sub_housingclear

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    !Indices
    INTEGER :: is, ig, ih, im, iy

    !Other internal variables
    REAL(8) :: supply_owner, supply_REO, hsupply_2, hdemand_2
    REAL(8), DIMENSION(num_threads_simul) :: hsupply_temp

!!!Compute housing supply - version 1
    yh           = zh*((1.d0 - landshare)*ph*zh)**((1.d0 - landshare)/landshare)
    hsupply		 = yh
    hsupply_temp = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          ih = hfun(ig)
          CALL omp_set_num_threads(num_threads_simul)
          !$omp parallel private(id,iy)
          !$omp do
          DO im = 1, nm
             id = omp_get_thread_num() + 1
             DO iy = ny_neg + 1, ny
                !Owners who sell voluntarily
                hsupply_temp(id) = hsupply_temp(id) + psell(iy,im,ig,is)*hh(ih)*mu1own(iy,im,ig,is)
             ENDDO
          ENDDO
          !$omp end do
          !$omp end parallel
       ENDDO
    ENDDO
    DO id = 1, num_threads_simul
       hsupply = hsupply + hsupply_temp(id)
    ENDDO
    !Supply of housing from REOs
    DO ih = 1, nh
       hsupply = hsupply + psell_reo(ih)*hh(ih)*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
!!!!Compute housing supply - version 2
    supply_owner = 0.d0
    supply_REO   = 0.d0
    DO ig = 1, ng
       ih = hfun(ig)
       supply_owner = supply_owner + hh(ih)*(SUM(mu1own(:,:,ig,:)) - SUM(mu2own(:,:,ig,:)))
    ENDDO
    DO ih = 1, nh
       supply_owner = supply_owner - hh(ih)*htot_reo_new(ih)
       supply_REO   = supply_REO + hh(ih)*(htot_reo_new(ih) + htot_reo(ih))*psell_reo(ih)
    ENDDO
    hsupply_2 = yh + supply_owner + supply_REO
!!!Compute housing demand - version 1
    hdemand = 0.d0
    DO is = 1, ns
       DO iy = ny_neg + 1, ny
          hdemand = hdemand + pbuy(iy,is)*hbuy(iy,is)*mu2rent(iy,is)
       ENDDO
    ENDDO
!!!Compute housing demand - version 2
    hdemand_2 = 0.d0
    DO ig = 1, ng
       ih = hfun(ig)
       hdemand_2 = hdemand_2 + hh(ih)*(SUM(mu3own(:,:,ig,:)) - SUM(mu2own(:,:,ig,:)))
    ENDDO
!!!Compute end of period housing stock
    hactual = 0.d0
    !End-of-period owner-occupied housing
    DO ig = 1, ng
       ih      = hfun(ig)
       hactual = hactual + hh(ih)*SUM(mu3own(:,:,ig,:))
    ENDDO
    !End-of-period unsold REO inventories
    DO ih = 1, nh
       hactual = hactual + (1.d0 - psell_reo(ih))*hh(ih)*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
    !Compute excess demand and update ph
    !		hdemand = hdemand_2			!Comment this line out to use version 1 of hdemand
    !		hsupply = hsupply_2			!Comment this line out to use version 1 of hsupply
    IF (algorithm_ph == 1) THEN
       demand = hdemand
       supply = hsupply
    ELSE
       demand = hactual
       supply = hstock_prev*(1.d0 - deprec_h) + yh
    ENDIF
    excess_demand1 = demand - supply
    IF (excess_demand1 >= 0.d0) THEN
       sign_excessdemand1 = 1
       IF (ph > ph_low) THEN
          ph_low           = ph
          demand_low       = demand
          supply_low       = supply
          excessdemand_low = excess_demand1
          mu3own_low       = mu3own
          mu3rent_low      = mu3rent
          htot_reo_new_low = htot_reo_new
       ENDIF
    ELSE
       sign_excessdemand1 = -1
       IF (ph < ph_upp) THEN
          ph_upp			 = ph
          demand_upp       = demand
          supply_upp       = supply
          excessdemand_upp = excess_demand1
          mu3own_upp       = mu3own
          mu3rent_upp      = mu3rent
          htot_reo_new_upp = htot_reo_new
       ENDIF
    ENDIF
    IF (iter_ph >=2 ) THEN
       IF (excess_demand0 >= 0.d0) THEN
          sign_excessdemand0 = 1
       ELSE
          sign_excessdemand0 = -1
       ENDIF
    ENDIF
    IF (iter_ph == 1) THEN
       IF (algorithm_ph == 1) THEN
          diff = ABS(excess_demand1)/(0.5*demand + 0.5*supply)
       ELSE
          diff = ABS(excess_demand1)
       ENDIF
       dph = MAX(0.5d0*(ph_low - ph),MIN(0.5d0*(ph_upp - ph),MIN(0.05d0,diff)*ph*sign_excessdemand1))		!Change min(0.05d0,... to adjust initial increment
    ELSEIF (sign_excessdemand1*sign_excessdemand0 == -1) THEN
       dph = MAX(0.5d0*(ph_low - ph),MIN(0.5d0*(ph_upp - ph),-0.5d0*dph))
       IF (algorithm_ph == 1) THEN
          diff = ABS(excess_demand1)/(0.5*demand + 0.5*supply)
       ELSE
          diff = ABS(excess_demand1)
       ENDIF
       IF (ABS(dph) > MIN(1.d0,diff)*ph) THEN
          dph = dph*MIN(1.d0,diff)*ph/ABS(dph)
       ENDIF
    ELSE
       dph = MAX(0.5d0*(ph_low - ph),MIN(0.5d0*(ph_upp - ph),dph))
       IF (algorithm_ph == 1) THEN
          diff = ABS(excess_demand1)/(0.5*demand + 0.5*supply)
       ELSE
          diff = ABS(excess_demand1)
       ENDIF
       IF (ABS(dph) > MIN(1.d0,diff)*ph) THEN
          dph = dph*MIN(1.d0,diff)*ph/ABS(dph)
       ENDIF
    ENDIF
    excess_demand0 = excess_demand1
    WRITE(*, '(a15, i4, a6, f9.5, a7, f9.5, a6, f9.6, a6, f9.6, a6, f8.5, a6, f8.5, a7, f8.5)'), '      Iter_ph =' ,iter_ph, ', ph =' ,ph, ', dph =' ,dph, ', DH =', hdemand, ', SH =', hsupply, ', Hd =', hactual, ', Hs =', hstock, ', own =', SUM(mu3own)

  END SUBROUTINE sub_housingclear

  SUBROUTINE sub_marketclear_transition

    USE mod_globvar

    IMPLICIT NONE

    INTEGER :: is, ig, ih, iy, im, imprime, isprime, ieprime, ie
    REAL(8) :: qm0price

    !******* External functions to be called *******
    REAL(8), EXTERNAL :: agglaborsupply

    !Housing stock
    hownstock = 0.d0
    DO ig = 1, ng
       ih = hfun(ig)
       hownstock = hownstock + hh(ih)*SUM(mu3own(:,:,ig,:))
    ENDDO

    !Mortgage revenues, costs, and net payments
    mortrev     = 0.d0
    mortcost    = 0.d0
    svcorigcost = 0.d0
    netmortpay  = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          ih = hfun(ig)
          !Mortgage Costs (Subperiod 3)
          DO iy = 1, ny
             !New mortgages
             imprime     = locmown(iy,1,ig,is)
             qm0price    = (1.d0 - bown_weight(iy,1,ig,is))*qm0(imprime,locbown_low(iy,1,ig,is),ig,is) + bown_weight(iy,1,ig,is)*qm0(imprime,locbown_upp(iy,1,ig,is),ig,is)
             mortcost    = mortcost + (1.d0 + rocost)*(1.d0 + rscost)*qm0price*mm(imprime)*hh(ih)*mu3own(iy,1,ig,is)
             netmortpay  = netmortpay - qm0price*mm(imprime)*hh(ih)*mu3own(iy,1,ig,is)
             svcorigcost = svcorigcost + ((1.d0 + rocost)*(1.d0 + rscost) - 1.d0)*qm0price*mm(imprime)*hh(ih)*mu3own(iy,1,ig,is)
          ENDDO
          DO im = 2, nm
             DO iy = ny_neg + 1, ny
                !Vintage mortgages
                qm0price = (1.d0 - bown_weight(iy,im,ig,is))*qm0(im,locbown_low(iy,im,ig,is),ig,is) + bown_weight(iy,im,ig,is)*qm0(im,locbown_upp(iy,im,ig,is),ig,is)
                mortcost = mortcost + (1.d0 + rocost)*(1.d0 + rscost)*qm0price*mm(im)*hh(ih)*mu3own(iy,im,ig,is)
             ENDDO
          ENDDO
          !Mortgage Revenues (Subperiod 1 next period)
          DO im = 2, nm
             DO iy = ny_neg + 1, ny
                !Homeowners that sell
                mortrev     = mortrev + 1.d0/(1.d0 + pi)*mm(im)*hh(ih)*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                netmortpay  = netmortpay + 1.d0/(1.d0 + pi)*mm(im)*hh(ih)*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                !Homeowners that enter foreclosure (see below for REOs)
                mortrev     = mortrev - MAX(0.d0,Jreo(ih) - 1.d0/(1.d0 + pi)*mm(im)*hh(ih))*(1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
                !Homeowners that default but stay out of foreclosure
                qm0price    = (1.d0 - bown_weight(iy,im,ig,is))*qm0(im,locbown_low(iy,im,ig,is),ig,is) + bown_weight(iy,im,ig,is)*qm0(im,locbown_upp(iy,im,ig,is),ig,is)
                mortrev     = mortrev + (1.d0 + rocost)*qm0price*1.d0/(1.d0 + pi)*mm(im)*hh(ih)*(1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*(1.d0 - pforeclose)*mu1own(iy,im,ig,is)
                !Homeowners that make a payment (payment + value of continuation contract)
                qm0price    = (1.d0 - bown_sp1_weight(iy,im,ig,is))*qm0(locmown(iy,im,ig,is),locbown_sp1_low(iy,im,ig,is),ig,is) + bown_sp1_weight(iy,im,ig,is)*qm0(locmown(iy,im,ig,is),locbown_sp1_upp(iy,im,ig,is),ig,is)
                mortrev     = mortrev + (1.d0/(1.d0 + pi)*mm(im)*hh(ih) - qm*mown(iy,im,ig,is)*hh(ih) + (1.d0 + rocost)*qm0price*mown(iy,im,ig,is)*hh(ih))*(1.d0 - psell(iy,im,ig,is))*(1.d0 - default(iy,im,ig,is))*mu1own(iy,im,ig,is)
                netmortpay  = netmortpay + (1.d0/(1.d0 + pi)*mm(im)*hh(ih) - qm*mown(iy,im,ig,is)*hh(ih))*(1.d0 - psell(iy,im,ig,is))*(1.d0 - default(iy,im,ig,is))*mu1own(iy,im,ig,is)
                svcorigcost = svcorigcost + rscost*(1.d0 + rocost)*qm0price*mown(iy,im,ig,is)*hh(ih)*(1.d0 - psell(iy,im,ig,is))*(1.d0 - default(iy,im,ig,is))*mu1own(iy,im,ig,is)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    !Purchases and sales of REOs
    DO ih = 1, nh
       mortrev  = mortrev + Jreo(ih)*(htot_reo(ih) + htot_reo_new(ih))
       mortcost = mortcost + (1.d0 - deprec_h)*qb*Jreo_tomorrow(ih)*(htot_reo(ih) + htot_reo_new(ih))*(1.d0 - psell_reo(ih))
       !mortcost = mortcost + (Jreo(ih) - psell_reo(ih)*xs_reo(ih) + (1.d0 - psell_reo(ih))*usercost*ph*hh(ih))*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO

    !Household bonds
    bonds_hh = 0.d0
    DO is = 1, ns
       DO iy = ny_neg + 1, ny
          !Renters
          bonds_hh = bonds_hh + brent(iy,is)*mu3rent(iy,is)
       ENDDO
       DO ig = 1, ng
          !Homeowners without a mortgage
          DO iy = 1, ny
             bonds_hh = bonds_hh + bown(iy,1,ig,is)*mu3own(iy,1,ig,is)
          ENDDO
          DO im = 2, nm
             !Homeowners with mortgages
             DO iy = ny_neg + 1, ny
                bonds_hh = bonds_hh + bown(iy,im,ig,is)*mu3own(iy,im,ig,is)
             ENDDO
          ENDDO
       ENDDO
    ENDDO

!!!Check Walras' law!!!
    aggsupply = tfp*glabor

    !Household consumption
    cons_hh   = 0.d0
    cons_rent = 0.d0
    DO is = 1, ns
       DO iy = ny_neg + 1, ny
          !Renters
          cons_hh = cons_hh + crent(iy,is)*mu3rent(iy,is)
          cons_rent = cons_rent + hcrent(iy,is)*mu3rent(iy,is)
       ENDDO
       DO ig = 1, ng
          !Homeowners without a mortgage
          DO iy = 1, ny
             cons_hh = cons_hh + cown(iy,1,ig,is)*mu3own(iy,1,ig,is)
          ENDDO
          DO im = 2, nm
             !Homeowners with mortgages
             DO iy = ny_neg + 1, ny
                cons_hh = cons_hh + cown(iy,im,ig,is)*mu3own(iy,im,ig,is)
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    !Transaction costs from buying houses
    buycost     = 0.d0
    nethousebuy = 0.d0
    DO is = 1, ns
       DO iy = ny_neg + 1, ny
          buycost     = buycost + (xbuy(iy,is) - ph*hbuy(iy,is))*pbuy(iy,is)*mu2rent(iy,is)
          nethousebuy = nethousebuy + xbuy(iy,is)*pbuy(iy,is)*mu2rent(iy,is)
       ENDDO
    ENDDO

    !Transaction costs from selling houses
    sellcost_hh  = 0.d0
    sellcost_REO = 0.d0
    REOsales     = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          ih = hfun(ig)
          DO im = 1, nm
             DO iy = ny_neg + 1, ny
                !Homeowners who sell voluntarily
                sellcost_hh = sellcost_hh + (ph*hh(ih) - xsell(iy,im,ig,is))*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                nethousebuy = nethousebuy - xsell(iy,im,ig,is)*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    !REO sales
    DO ih = 1, nh
       sellcost_REO = sellcost_REO + (ph*hh(ih) - xs_reo(ih))*psell_reo(ih)*(htot_reo(ih) + htot_reo_new(ih))
       REOsales     = REOsales + hh(ih)*psell_reo(ih)*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
    sellcost = sellcost_hh + sellcost_REO

    !Housing maintenance expenses
    maint_hh  = 0.d0
    maint_REO = 0.d0
    DO ig = 1, ng
       ih       = hfun(ig)
       maint_hh = maint_hh + usercost*ph*hh(ih)*SUM(mu3own(:,1,ig,:))
       DO iy = ny_neg + 1, ny
          maint_hh = maint_hh + MIN(yy(iy),usercost*ph*hh(ih))*SUM(mu3own(iy,2:nm,ig,:))
       ENDDO
    ENDDO
    !REOs
    DO ih = 1, nh
       maint_REO = maint_REO + usercost*ph*hh(ih)*(1.d0 - psell_reo(ih))*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
    maint = maint_hh + maint_REO

    !Resources spent on construction
    sh = ((1 - landshare)*ph*zh)**(1.d0/landshare)
    pl = ph*zh*landshare*sh**(1.d0 - landshare)

    !Aggregate demand (doesn't include pricing adjustment costs because they are zero in steady state)
!!!NOTE: maint vs. maint_hh?
    aggdemand = cons_hh + pa*cons_rent + govspend + svcorigcost + maint + sellcost + buycost + sh + pl + adjustcost + firmFC + bankFC

  END SUBROUTINE sub_marketclear_transition

  SUBROUTINE sub_timeseries

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices
    INTEGER :: is, ig, ih, iu, im, iy

    !Foreclosures
    REAL(8) :: foreclosures, frate
    !Housing sales, inventories, and time on market
    REAL(8) :: existsales, inventories, months_supply, reosaleshare, htot, salesrate, existsales_noreo, inventories_noreo, months_supply_noreo, mean_tomsell, meas_sell, mean_tomsell_prefshocks, &
         meas_sell_prefshocks, mean_TOMREO

    !Foreclosure rate
    foreclosures = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          DO im = 2, nm
             DO iy = ny_neg + 1, ny
                foreclosures = foreclosures + (1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*mu1own(iy,im,ig,is)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    IF (SUM(mu1own(:,2:nm,:,:)) > 0.d0) THEN
       frate = foreclosures/SUM(mu1own(:,2:nm,:,:))
    ELSE
       frate = 0.d0
    ENDIF

    !Sales volume, inventories, months of supply, REO sales share, and sales rate (sales is the measure of housing traded, inventories is the measure of housing put up for sale but not sold, months of supply is annualized
    !inventories/sales, REO sales share is REO sales volume as a share of existing home sales, and sales rate is existing sales as a share of the housing stock
    existsales  = 0.d0
    inventories = 0.d0
    reosales    = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          ih = hfun(ig)
          !Homeowners
          DO im = 1, nm
             DO iy = ny_neg + 1, ny
                !Owners who voluntarily sell
                IF (psell(iy,im,ig,is) > 0.d0) THEN
                   existsales  = existsales + psell(iy,im,ig,is)*hh(ih)*mu1own(iy,im,ig,is)
                   inventories = inventories + (1.d0 - psell(iy,im,ig,is))*hh(ih)*mu1own(iy,im,ig,is)
                ENDIF
                !Owners who default
                IF (psell_reo(ih) > 0.d0) THEN
                   !Sales volume and inventories from new REOs
                   existsales  = existsales + hh(ih)*(1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*psell_reo(ih)*mu1own(iy,im,ig,is)
                   reosales    = reosales + hh(ih)*(1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*psell_reo(ih)*mu1own(iy,im,ig,is)
                   inventories = inventories + hh(ih)*(1.d0 - psell(iy,im,ig,is))*default(iy,im,ig,is)*pforeclose*(1.d0 - psell_reo(ih))*mu1own(iy,im,ig,is)
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    !Sales volume and inventories from existing REOs
    DO ih = 1, nh
       existsales  = existsales + hh(ih)*psell_reo(ih)*(htot_reo(ih) + htot_reo_new(ih))
       reosales    = reosales + hh(ih)*psell_reo(ih)*(htot_reo(ih) + htot_reo_new(ih))
       inventories = inventories + hh(ih)*(1.d0 - psell_reo(ih))*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
    months_supply = inventories/(4.0*existsales)*12.0
    reosaleshare  = reosales/existsales
    !Compute beginning of period housing stock and sales rate
    htot = 0.d0
    DO ig = 1, ng
       ih   = hfun(ig)
       htot = htot + hh(ih)*SUM(mu1own(:,:,ig,:))
    ENDDO
    DO ih = 1, nh
       htot = htot + hh(ih)*(htot_reo(ih) + htot_reo_new(ih))
    ENDDO
    salesrate = 4.0*existsales/htot

    !Months supply excluding REOs
    existsales_noREO  = 0.d0
    inventories_noREO = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          ih = hfun(ig)
          !Homeowners
          DO im = 1, nm
             DO iy = ny_neg + 1, ny
                !Owners who voluntarily sell
                IF (psell(iy,im,ig,is) > 0.d0) THEN
                   existsales_noREO  = existsales_noREO + psell(iy,im,ig,is)*hh(ih)*mu1own(iy,im,ig,is)
                   inventories_noREO = inventories_noREO + (1.d0 - psell(iy,im,ig,is))*hh(ih)*mu1own(iy,im,ig,is)
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    months_supply_noREO = inventories_noREO/(4.0*existsales_noREO)*12.0

    !Mean seller search duration
    mean_TOMsell            = 0.d0
    meas_sell               = 0.d0
    mean_TOMsell_prefshocks = 0.d0
    meas_sell_prefshocks    = 0.d0
    DO is = 1, ns
       DO ig = 1, ng
          iu = ufun(ig)
          !Sellers
          DO im = 1, nm
             DO iy =  ny_neg + 1, ny
                IF (psell(iy,im,ig,is) > 0.d0) THEN
                   mean_TOMsell = mean_TOMsell + ((1.d0 - psell(iy,im,ig,is))/psell(iy,im,ig,is) + 0.5)*12.0*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                   meas_sell    = meas_sell + psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                   IF (iu == 2) THEN
                      mean_TOMsell_prefshocks = mean_TOMsell_prefshocks + ((1.d0 - psell(iy,im,ig,is))/psell(iy,im,ig,is) + 0.5)*12.0*psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                      meas_sell_prefshocks    = meas_sell_prefshocks + psell(iy,im,ig,is)*mu1own(iy,im,ig,is)
                   ENDIF
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    IF (meas_sell > 1e-6) THEN
       mean_TOMsell = mean_TOMsell/meas_sell
    ELSE
       mean_TOMsell = 0.d0
    ENDIF
    IF (meas_sell_prefshocks > 1e-6) THEN
       mean_TOMsell_prefshocks = mean_TOMsell_prefshocks/meas_sell_prefshocks
    ELSE
       mean_TOMsell_prefshocks = 0.d0
    ENDIF

    !Compute REO sale discount and search duration
    IF (SUM(htot_reo) + SUM(htot_reo_new) > 0.d0) THEN
       DO ih = 1, nh
          IF (psell_reo(ih) > 0.d0) THEN
             mean_TOMREO = mean_TOMREO + ((1.d0 - psell_reo(ih))/psell_reo(ih) + 0.5)*12.0*(htot_reo(ih) + htot_reo_new(ih))/(SUM(htot_reo) + SUM(htot_reo_new))
          ENDIF
       ENDDO
    ENDIF

    series_full(it + 1,1)  = ph
    series_full(it + 1,2)  = (1.d0 + rf)**4.d0 - 1.d0
    series_full(it + 1,3)  = (1.d0 + pi)**4.d0 - 1.d0
    series_full(it + 1,4)  = w
    series_full(it + 1,5)  = lumpsum
    series_full(it + 1,6)  = aggsupply
    series_full(it + 1,7)  = aggdemand
    series_full(it + 1,8)  = cons_hh
    series_full(it + 1,9)  = pa*cons_rent
    series_full(it + 1,10) = govspend
    series_full(it + 1,11) = sh
    series_full(it + 1,12) = pl
    series_full(it + 1,13) = buycost
    series_full(it + 1,14) = sellcost
    series_full(it + 1,15) = svcorigcost
    series_full(it + 1,16) = maint
    series_full(it + 1,17) = adjustcost
    series_full(it + 1,18) = frate
    series_full(it + 1,19) = salesrate
    series_full(it + 1,20) = reosaleshare
    series_full(it + 1,21) = months_supply
    series_full(it + 1,22) = months_supply_noREO
    series_full(it + 1,23) = mean_TOMsell
    series_full(it + 1,24) = mean_TOMsell_prefshocks
    series_full(it + 1,25) = mean_TOMREO

  END SUBROUTINE sub_timeseries

  SUBROUTINE sub_initialize_panel

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !********* Variables for panel simulation *********!!!
    INTEGER :: ig, iy, im, is, ie
    REAL(8), DIMENSION(:), ALLOCATABLE :: rand_mu_own, rand_mu_s_own, rand_mu_y_srent, rand_mu_g_sown, rand_mu_m_gsown, rand_mu_y_mgsown

    !Allocate variables for panel
    ALLOCATE(rand_psell_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_foreclose_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_pbuy_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_hseparate_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_hdeprec_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_survive_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_weight_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_s_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_e_nt(num_id,T_init + T + 1))
    ALLOCATE(rand_all_nt(num_id*(T_init + T + 1),10))
    ALLOCATE(iy_nt(num_id,T_init + T + 1))
    ALLOCATE(im_nt(num_id,T_init + T + 1))
    ALLOCATE(ig_nt(num_id,T_init + T + 1))
    ALLOCATE(is_nt(num_id,T_init + T + 1))
    ALLOCATE(ie_nt(num_id,T_init + T + 1))
    !Panel output
    ALLOCATE(cons_nt(num_id,T_init + T + 1))
    ALLOCATE(earn_nt(num_id,T_init + T + 1))
    ALLOCATE(nonlaborincome_nt(num_id,T_init + T + 1))
    ALLOCATE(bonds_nt(num_id,T_init + T + 1))
    ALLOCATE(xsell_nt(num_id,T_init + T + 1))
    ALLOCATE(psell_nt(num_id,T_init + T + 1))
    ALLOCATE(sell_nt(num_id,T_init + T + 1))
    ALLOCATE(default_nt(num_id,T_init + T + 1))
    ALLOCATE(foreclose_nt(num_id,T_init + T + 1))
    ALLOCATE(phsim_nt(num_id,T_init + T + 1))
    ALLOCATE(rfsim_nt(num_id,T_init + T + 1))
    ALLOCATE(rmsim_nt(num_id,T_init + T + 1))
    ALLOCATE(pisim_nt(num_id,T_init + T + 1))
    ALLOCATE(wsim_nt(num_id,T_init + T + 1))
    !Random number generators
    ALLOCATE(rand_mu_own(num_id))
    ALLOCATE(rand_mu_s_own(num_id))
    ALLOCATE(rand_mu_y_srent(num_id))
    ALLOCATE(rand_mu_g_sown(num_id))
    ALLOCATE(rand_mu_m_gsown(num_id))
    ALLOCATE(rand_mu_y_mgsown(num_id))
    CALL random_seed

    !Initialize panel
    IF (load_panel == 1) THEN
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_initialdraw.txt')
       DO in = 1, num_id
          READ (11, '(i5, i5, i5, i5)'), iy_nt(in,1), im_nt(in,1), ig_nt(in,1), is_nt(in,1)
       ENDDO
       CLOSE(11)
    ELSE
       PRINT *, ' '
       PRINT *, 'Initialize panel'
       PRINT *, ' '
       CALL omp_set_num_threads(num_threads)
       !$omp parallel private(is,iy,ig,im)
       !$omp do
       DO in = 1, num_id
          CALL RANDOM_NUMBER(rand_mu_own(in))
          CALL RANDOM_NUMBER(rand_mu_s_own(in))
          !Determine homeownership status
          IF (rand_mu_own(in) <= SUM(mu1rent)/(SUM(mu1rent) + SUM(mu1own))) THEN
             CALL RANDOM_NUMBER(rand_mu_y_srent(in))
             !Renter
             ig_nt(in,1) = 0
             im_nt(in,1) = 1
             !Labor shock
             DO is = 1, ns
                IF (rand_mu_s_own(in) <= SUM(mu1rent(:,1:is))/SUM(mu1rent)) THEN
                   is_nt(in,1) = is
                   !Determine cash at hand
                   DO iy = ny_neg + 1, ny
                      IF (rand_mu_y_srent(in) <= SUM(mu1rent(1:iy,is))/SUM(mu1rent(:,is))) THEN
                         iy_nt(in,1) = iy
                         EXIT
                      ENDIF
                   ENDDO
                   EXIT
                ENDIF
             ENDDO
          ELSE
             !Homeowner
             CALL RANDOM_NUMBER(rand_mu_g_sown(in))
             CALL RANDOM_NUMBER(rand_mu_m_gsown(in))
             CALL RANDOM_NUMBER(rand_mu_y_mgsown(in))
             !Labor shock
             DO is = 1, ns
                IF (rand_mu_s_own(in) <= SUM(mu1own(:,:,:,1:is))/SUM(mu1own)) THEN
                   is_nt(in,1) = is
                   !Determine house size
                   DO ig = 1, ng
                      IF (rand_mu_g_sown(in) <= SUM(mu1own(:,:,1:ig,is))/SUM(mu1own(:,:,:,is))) THEN
                         ig_nt(in,1) = ig
                         !Determine mortgage debt
                         DO im = 1, nm
                            IF (rand_mu_m_gsown(in) <= SUM(mu1own(:,1:im,ig,is))/SUM(mu1own(:,:,ig,is))) THEN
                               im_nt(in,1) = im
                               !Determine cash at hand
                               DO iy = ny_neg + 1, ny
                                  IF (rand_mu_y_mgsown(in) <= SUM(mu1own(1:iy,im,ig,is))/SUM(mu1own(:,im,ig,is))) THEN
                                     iy_nt(in,1) = iy
                                     EXIT
                                  ENDIF
                               ENDDO
                               EXIT
                            ENDIF
                         ENDDO
                         EXIT
                      ENDIF
                   ENDDO
                   EXIT
                ENDIF
             ENDDO
          ENDIF
       ENDDO
       !$omp end do
       !$omp end parallel
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_initialdraw.txt', status = 'replace')
       DO in = 1, num_id
          WRITE(11, '(i5, i5, i5, i5)'), iy_nt(in,1), im_nt(in,1), ig_nt(in,1), is_nt(in,1)
       ENDDO
       CLOSE(11)
    ENDIF

    !Generate stochastic realizations for panel simulation
    IF (load_stochastic == 1) THEN
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_psell_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_psell_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_foreclose_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_foreclose_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_pbuy_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_pbuy_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_hseparate_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_hseparate_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_hdeprec_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_hdeprec_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_survive_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_survive_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_weight_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_weight_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_s_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_s_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_e_nt.txt')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             READ (11, *) rand_e_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
    ELSE
       CALL random_seed
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             CALL RANDOM_NUMBER(rand_psell_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_foreclose_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_pbuy_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_hseparate_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_hdeprec_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_survive_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_weight_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_s_nt(in,itsim))
             CALL RANDOM_NUMBER(rand_e_nt(in,itsim))
          ENDDO
       ENDDO
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_psell_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_psell_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_foreclose_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_foreclose_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_pbuy_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_pbuy_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_hseparate_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_hseparate_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_hdeprec_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_hdeprec_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_survive_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_survive_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_weight_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_weight_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_s_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_s_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
       OPEN(unit=11,	file='Input'//sep//'Panel'//sep//'panel_rand_e_nt.txt', status = 'replace')
       DO in = 1, num_id
          DO itsim = 1, T_init + T + 1
             WRITE (11, '(1f12.8)'), rand_e_nt(in,itsim)
          ENDDO
       ENDDO
       CLOSE(11)
    ENDIF
    DO itsim = 1, T_init + T + 1
       it = MAX(itsim - T_init,0)
       IF (it == 0) THEN
          CALL sub_transition_phase(0)
       ELSE
          CALL sub_transition_phase(phase_t(it))
       ENDIF
       DO in = 1, num_id
          DO ie = 1, ne
             IF (rand_e_nt(in,itsim) <= SUM(emeas(1:ie))) THEN
                ie_nt(in,itsim) = ie
                EXIT
             ENDIF
          ENDDO
          IF (itsim >= 2) THEN
             DO is = 1, ns
                IF (rand_s_nt(in,itsim) <= SUM(prob_ss(is_nt(in,itsim - 1),1:is))) THEN
                   is_nt(in,itsim) = is
                   EXIT
                ENDIF
             ENDDO
          ENDIF
       ENDDO
    ENDDO

    !Simulate panel during initial stationary equilibrium
    CALL sub_transition_phase(0)
    IF (T_init > 1) THEN
       DO itsim = 1, T_init - 1
          CALL sub_simulate
          WRITE(*, '(a7, i4)'), ' tsim =', itsim
       ENDDO
    ENDIF

    !Initialize panel at period 1 of the transition
    itsim = T_init
    !Need to change timing of simulation? Subperiod 3 from t - 1 followed by subperiods 1 and 2 in period t?
    CALL sub_simulate
    WRITE(*, '(a7, i4)'), ' tsim =', itsim

  END SUBROUTINE sub_initialize_panel

  SUBROUTINE sub_simulate

    USE mod_globvar
    USE mod_interpolate
    USE omp_lib

    IMPLICIT NONE

    !Indices
    INTEGER :: in, is, iy_sp1, im_sp1, ig_sp1, ih_sp1, iy_sp2, im_sp2, ig_sp2, ih_sp2, iy_sp3, im_sp3, ig_sp3, ih_sp3, isprime, ieprime

    !Other internal variables
    REAL(8), DIMENSION(ns) :: laborchoice, disutility
    REAL(8) :: rand_weight, rand_psell, rand_foreclose, ysell, yweight, ydefault, ypay, rand_pbuy, ybuy, rand_hseparate, rand_hdeprec, rand_survive, yprime
    INTEGER :: depreciate, iylow, iyupp

    !******* External functions to be called *******
    REAL(8), EXTERNAL :: laborsupply, labordisutil

    !Labor supply and disutility
    DO is = 1, ns
       laborchoice(is) = laborsupply(w,ss(is))
       disutility(is)  = labordisutil(laborchoice(is),ss(is))
    ENDDO

    CALL omp_set_num_threads(num_threads)
    !$omp parallel private(is,iy_sp1,im_sp1,ig_sp1,ih_sp1,iy_sp2,im_sp2,ig_sp2,ih_sp2,iy_sp3,im_sp3,ig_sp3,ih_sp3,isprime,ieprime,depreciate,iylow,iyupp,rand_weight,rand_psell,rand_foreclose,ysell,yweight,ydefault,ypay,rand_pbuy,ybuy,rand_hseparate,rand_hdeprec,rand_survive,yprime)
    !$omp do
    DO in = 1, num_id
       rand_weight = 0.5d0		!rand_weight_nt(in,itsim)		!Comment out this line to do stochastic interpolation of cash at hand in *all* subperiods
       !Output
       phsim_nt(in,itsim) = ph
       rfsim_nt(in,itsim) = (1.d0 + rf)**4 - 1.d0
       rmsim_nt(in,itsim) = (1.d0/qm)**4 - 1.d0
       pisim_nt(in,itsim) = (1.d0 + pi)**4 - 1.d0
!!!Subperiod 1
       !Individual states
       iy_sp1 = iy_nt(in,itsim)
       im_sp1 = im_nt(in,itsim)
       ig_sp1 = ig_nt(in,itsim)
       IF (ig_sp1 == 0) THEN
          ih_sp1 = 0
       ELSE
          ih_sp1 = hfun(ig_sp1)
       ENDIF
       is = is_nt(in,itsim)
       !Stochastic outcomes
       rand_psell     = rand_psell_nt(in,itsim)
       rand_foreclose = rand_foreclose_nt(in,itsim)
       IF (ig_sp1 >= 1) THEN
          !Homeowners
          IF (im_sp1 == 1) THEN
             psell_nt(in,itsim)     = psell(iy_sp1,1,ig_sp1,is)
             xsell_nt(in,itsim)     = xsratio(iy_sp1,1,ig_sp1,is)
             default_nt(in,itsim)   = 0.d0
             foreclose_nt(in,itsim) = 0.d0
             !Homeowners with no mortgage
             IF (rand_psell < psell(iy_sp1,1,ig_sp1,is) .OR. psell(iy_sp1,1,ig_sp1,is) > 1.d0 - 1e-4) THEN
                sell_nt(in,itsim) = 1.d0
                !Successful sale
                im_sp2  = 1
                ig_sp2  = 0
                ysell   = yy(iy_sp1) + xsell(iy_sp1,1,ig_sp1,is)
                iylow   = ly(ysell)
                iyupp   = uy(ysell)
                yweight = wy(ysell)
                IF (yweight <= rand_weight) THEN
                   iy_sp2 = iylow
                ELSE
                   iy_sp2 = iyupp
                ENDIF
             ELSE
                sell_nt(in,itsim) = 0
                !No sale
                iy_sp2 = iy_sp1
                im_sp2 = 1
                ig_sp2 = ig_sp1
             ENDIF
          ELSE
             psell_nt(in,itsim)   = psell(iy_sp1,im_sp1,ig_sp1,is)
             xsell_nt(in,itsim)   = xsratio(iy_sp1,im_sp1,ig_sp1,is)
             !Homeowners with a mortgage
             IF (rand_psell < psell(iy_sp1,im_sp1,ig_sp1,is) .OR. psell(iy_sp1,im_sp1,ig_sp1,is) > 1.d0 - 1e-4) THEN
                !Successful sale
                sell_nt(in,itsim)      = 1.d0
                default_nt(in,itsim)   = 0.d0
                foreclose_nt(in,itsim) = 0.d0
                im_sp2	= 1
                ig_sp2  = 0
                ysell   = yy(iy_sp1) - 1.d0/(1.d0 + pi)*mm(im_sp1)*hh(ih_sp1) + xsell(iy_sp1,im_sp1,ig_sp1,is)
                iylow   = ly(ysell)
                iyupp   = uy(ysell)
                yweight = wy(ysell)
                IF (yweight <= rand_weight) THEN
                   iy_sp2 = iylow
                ELSE
                   iy_sp2 = iyupp
                ENDIF
             ELSEIF (rand_psell > psell(iy_sp1,im_sp1,ig_sp1,is) .AND. default(iy_sp1,im_sp1,ig_sp1,is) > 0.d0) THEN
                !Default
                sell_nt(in,itsim)	 = 0.d0
                default_nt(in,itsim) = 1.d0
                IF (rand_foreclose < pforeclose .OR. pforeclose > 1.d0 - 1e-4) THEN
                   !Bank forecloses
                   foreclose_nt(in,itsim) = 1.d0
                   im_sp2    = 1
                   ig_sp2    = 0
                   ydefault  = yy(iy_sp1) - MAX(0.d0,MIN(0.d0,1.d0/(1.d0 + pi)*mm(im_sp1)*hh(ih_sp1) - Jreo(ih_sp1)))
                   iylow     = ly(ydefault)
                   iyupp     = uy(ydefault)
                   yweight   = wy(ydefault)
                   IF (yweight <= rand_weight) THEN
                      iy_sp2 = iylow
                   ELSE
                      iy_sp2 = iyupp
                   ENDIF
                ELSE
                   !Bank does not foreclose
                   foreclose_nt(in,itsim) = 0.d0
                   !Bank does not foreclose
                   iy_sp2   = iy_sp1
                   im_sp2   = im_sp1
                   ig_sp2   = ig_sp1
                ENDIF
             ELSE
                sell_nt(in,itsim)      = 0.d0
                default_nt(in,itsim)   = 0.d0
                foreclose_nt(in,itsim) = 0.d0
                !Make a regular payment
                ig_sp2 = ig_sp1
                im_sp2 = locmown(iy_sp1,im_sp1,ig_sp1,is)
                ypay   = yy(iy_sp1) - (1.d0/(1.d0 + pi)*mm(im_sp1) - qm*mown(iy_sp1,im_sp1,ig_sp1,is))*hh(ih_sp1)
                IF (ypay < 0.d0) THEN
                   ypay = ypay/ABS(ylb(ig_sp1,is))
                ENDIF
                iylow   = ly(ypay)
                iyupp   = uy(ypay)
                yweight = wy(ypay)
                IF (yweight <= rand_weight) THEN
                   iy_sp2 = iylow
                ELSE
                   iy_sp2 = iyupp
                ENDIF
             ENDIF
          ENDIF
       ELSE
          psell_nt(in,itsim)     = 0.d0
          xsell_nt(in,itsim)     = 0.d0
          sell_nt(in,itsim)      = 0.d0
          default_nt(in,itsim)   = 0.d0
          foreclose_nt(in,itsim) = 0.d0
          !Renters
          iy_sp2 = iy_sp1
          im_sp2 = 1
          ig_sp2 = 0
       ENDIF
!!!Subperiod 2
       !Stochastic outcomes
       rand_pbuy = rand_pbuy_nt(in,itsim)
       IF (ig_sp2 >= 1) THEN
          !Homeowners
          iy_sp3 = iy_sp2
          im_sp3 = im_sp2
          ig_sp3 = ig_sp2
       ELSE
          !Renters
          IF (rand_pbuy < pbuy(iy_sp2,is) .OR. pbuy(iy_sp2,is) > 1.d0 - 1e-4) THEN
             !Successful purchase
             im_sp3 = 1
             ih_sp3 = lochbuy(iy_sp2,is)
             ig_sp3 = gfun(ih_sp3,1)
             ybuy   = yy(iy_sp2) - xbuy(iy_sp2,is)
             IF (ybuy < 0.d0) THEN
                ybuy = ybuy/ABS(ylb(ig_sp3,is))
             ENDIF
             iylow   = ly(ybuy)
             iyupp   = uy(ybuy)
             yweight = wy(ybuy)
             IF (yweight <= rand_weight) THEN
                iy_sp3 = iylow
             ELSE
                iy_sp3 = iyupp
             ENDIF
          ELSE
             iy_sp3 = iy_sp2
             im_sp3 = 1
             ig_sp3 = 0
          ENDIF
       ENDIF
!!!Subperiod 3
       rand_weight    = rand_weight_nt(in,itsim + 1)
       rand_hseparate = rand_hseparate_nt(in,itsim + 1)
       rand_hdeprec   = rand_hdeprec_nt(in,itsim + 1)
       rand_survive   = rand_survive_nt(in,itsim + 1)
       isprime        = is_nt(in,itsim + 1)
       ieprime        = ie_nt(in,itsim + 1)
       !Output
       earn_nt(in,itsim + 1)           = labor(isprime,ieprime)*laborchoice(isprime)*w
       nonlaborincome_nt(in,itsim + 1) = lumpsum + dividends
       wsim_nt(in,itsim + 1)           = w
       IF (ig_sp3 >= 1) THEN
          !Homeowners
          IF (im_sp3 == 1) THEN
             !Homeowners without a mortgage
             cons_nt(in,itsim)  = cown(iy_sp3,1,ig_sp3,is)
             bonds_nt(in,itsim) = bown(iy_sp3,1,ig_sp3,is)
             IF (rand_survive < psurvive) THEN
                !Households who survive
                IF (rand_hdeprec < deprec_h_survive) THEN
                   depreciate = 1
                   !House fully depreciates
                   yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends + bown(iy_sp3,1,ig_sp3,is))
                   iylow   = ly(yprime)
                   iyupp   = uy(yprime)
                   yweight = wy(yprime)
                   IF (yweight <= rand_weight) THEN
                      iy_nt(in,itsim + 1) = iylow
                   ELSE
                      iy_nt(in,itsim + 1) = iyupp
                   ENDIF
                   im_nt(in,itsim + 1) = 1
                   ig_nt(in,itsim + 1) = 0
                ELSE
                   depreciate = 0
                   !House does not depreciate
                   yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends + bown(iy_sp3,1,ig_sp3,is))
                   iylow   = ly(yprime)
                   iyupp   = uy(yprime)
                   yweight = wy(yprime)
                   IF (yweight <= rand_weight) THEN
                      iy_nt(in,itsim + 1) = iylow
                   ELSE
                      iy_nt(in,itsim + 1) = iyupp
                   ENDIF
                   im_nt(in,itsim + 1) = locmown(iy_sp3,1,ig_sp3,is)
                   IF (ufun(ig_sp3) == 1 .AND. nu == 2 .AND. rand_hseparate < pseparate) THEN
                      !Housing shock
                      ig_nt(in,itsim + 1) = gfun(hfun(ig_sp3),nu)
                   ELSE
                      !No housing shock
                      ig_nt(in,itsim + 1) = ig_sp3
                   ENDIF
                ENDIF
             ELSE
                !Households who die and are replaced by newborns
                yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends)
                iylow   = ly(yprime)
                iyupp   = uy(yprime)
                yweight = wy(yprime)
                IF (yweight <= rand_weight) THEN
                   iy_nt(in,itsim + 1) = iylow
                ELSE
                   iy_nt(in,itsim + 1) = iyupp
                ENDIF
                im_nt(in,itsim + 1) = 1
                ig_nt(in,itsim + 1) = 0
             ENDIF
          ELSE
             !Homeowners with a mortgage
             cons_nt(in,itsim)  = cown(iy_sp3,im_sp3,ig_sp3,is)
             bonds_nt(in,itsim) = bown(iy_sp3,im_sp3,ig_sp3,is)
             IF (rand_survive < psurvive) THEN
                !Households who survive
                IF (rand_hdeprec < deprec_h_survive) THEN
                   depreciate = 1
                   !House fully depreciates
                   yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends + bown(iy_sp3,im_sp3,ig_sp3,is))
                   iylow   = ly(yprime)
                   iyupp   = uy(yprime)
                   yweight = wy(yprime)
                   IF (yweight <= rand_weight) THEN
                      iy_nt(in,itsim + 1) = iylow
                   ELSE
                      iy_nt(in,itsim + 1) = iyupp
                   ENDIF
                   im_nt(in,itsim + 1) = 1
                   ig_nt(in,itsim + 1) = 0
                ELSE
                   depreciate = 0
                   !House does not depreciate
                   yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends + bown(iy_sp3,im_sp3,ig_sp3,is))
                   iylow   = ly(yprime)
                   iyupp   = uy(yprime)
                   yweight = wy(yprime)
                   IF (yweight <= rand_weight) THEN
                      iy_nt(in,itsim + 1) = iylow
                   ELSE
                      iy_nt(in,itsim + 1) = iyupp
                   ENDIF
                   im_nt(in,itsim + 1) = im_sp3
                   IF (ufun(ig_sp3) == 1 .AND. nu == 2 .AND. rand_hseparate < pseparate) THEN
                      !Housing shock
                      ig_nt(in,itsim + 1) = gfun(hfun(ig_sp3),nu)
                   ELSE
                      !No housing shock
                      ig_nt(in,itsim + 1) = ig_sp3
                   ENDIF
                ENDIF
             ELSE
                !Households who die and are replaced by newborns
                yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends)
                iylow   = ly(yprime)
                iyupp   = uy(yprime)
                yweight = wy(yprime)
                IF (yweight <= rand_weight) THEN
                   iy_nt(in,itsim + 1) = iylow
                ELSE
                   iy_nt(in,itsim + 1) = iyupp
                ENDIF
                im_nt(in,itsim + 1) = 1
                ig_nt(in,itsim + 1) = 0
             ENDIF
          ENDIF
       ELSE
          depreciate = 0
          !Renters
          cons_nt(in,itsim)  = crent(iy_sp3,is)
          bonds_nt(in,itsim) = brent(iy_sp3,is)
          IF (rand_survive < psurvive) THEN
             !Households who survive
             yprime             = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends + brent(iy_sp3,is))
             iylow              = ly(yprime)
             iyupp              = uy(yprime)
             yweight            = wy(yprime)
             IF (yweight <= rand_weight) THEN
                iy_nt(in,itsim + 1) = iylow
             ELSE
                iy_nt(in,itsim + 1) = iyupp
             ENDIF
             im_nt(in,itsim + 1) = 1
             ig_nt(in,itsim + 1) = 0
          ELSE
             !Households who die and are replaced by newborns
             yprime  = MAX(0.d0,(1.d0 - labortax)*labor(isprime,ieprime)*laborchoice(isprime)*w + lumpsum + dividends)
             iylow   = ly(yprime)
             iyupp   = uy(yprime)
             yweight = wy(yprime)
             IF (yweight <= rand_weight) THEN
                iy_nt(in,itsim + 1) = iylow
             ELSE
                iy_nt(in,itsim + 1) = iyupp
             ENDIF
             im_nt(in,itsim + 1) = 1
             ig_nt(in,itsim + 1) = 0
          ENDIF
       ENDIF
    ENDDO
    !$omp end do
    !$omp end parallel

  END SUBROUTINE sub_simulate

END MODULE mod_transition
