anchor.diff:  !------------------------------------------------------------------------------
   
       ! Construct auxiliary objects
  +    num_obs = SIZE(data_array, 1)
  +    ALLOCATE(crit_val(num_obs))
  +
       shocks_cov = MATMUL(shocks_cholesky, TRANSPOSE(shocks_cholesky))
       is_deterministic = ALL(shocks_cov .EQ. zero_dble)
   
       ! Initialize container for likelihood contributions
       crit_val = zero_dble
   
  -    j = 1
  -
  -    DO i = 0, num_agents_est - 1
  -
  -        DO period = 0, num_periods -1
  -
  -            ! Extract observable components of state space as well as agent decision.
  -            exp_a = INT(data_array(j, 5))
  -            exp_b = INT(data_array(j, 6))
  -            edu = INT(data_array(j, 7))
  -            edu_lagged = INT(data_array(j, 8))
  -
  -            choice = INT(data_array(j, 3))
  -            is_working = (choice == 1) .OR. (choice == 2)
  -
  -            ! Transform total years of education to additional years of education and create an index from the choice.
  -            edu = edu - edu_start
  -
  -            ! This is only done for alignment
  -            idx = choice
  -
  -            ! Get state indicator to obtain the systematic component of the agents payoffs. These feed into the simulation of choice probabilities.
  -            k = mapping_state_idx(period + 1, exp_a + 1, exp_b + 1, edu + 1, edu_lagged + 1)
  -            payoffs_systematic = periods_payoffs_systematic(period + 1, k + 1, :)
  -
  -            ! Extract relevant deviates from standard normal distribution.
  -            draws_prob_raw = periods_draws_prob(period + 1, :, :)
  -
  -            ! Prepare to calculate product of likelihood contributions.
  -            crit_val_contrib = 1.0
  -
  -            ! If an agent is observed working, then the the labor market shocks are observed and the conditional distribution is used to determine the choice probabilities.
  -            dist = zero_dble
  -            IF (is_working) THEN
  -
  -                ! Calculate the disturbance, which follows a normal distribution.
  -                dist = clip_value(LOG(data_array(j, 4)), -HUGE_FLOAT, HUGE_FLOAT) - clip_value(LOG(payoffs_systematic(idx)), -HUGE_FLOAT, HUGE_FLOAT) 
  -
  -                ! If there is no random variation in payoffs, then the observed wages need to be identical their systematic components. The discrepancy between the observed wages and their systematic components might be small due to the reading in of the dataset.
  -                IF (is_deterministic) THEN
  -                    IF (dist .GT. SMALL_FLOAT) THEN
  -                        rslt = zero_dble
  +    DO j = 1, num_obs
  +
  +        period = INT(data_array(j, 2))  
  +
  +        ! Extract observable components of state space as well as agent decision.
  +        exp_a = INT(data_array(j, 5))
  +        exp_b = INT(data_array(j, 6))
  +        edu = INT(data_array(j, 7))
  +        edu_lagged = INT(data_array(j, 8))
  +
  +        choice = INT(data_array(j, 3))
  +        is_working = (choice == 1) .OR. (choice == 2)
  +
  +        ! Transform total years of education to additional years of education and create an index from the choice.
  +        edu = edu - edu_start
  +
  +        ! This is only done for alignment
  +        idx = choice
  +
  +        ! Get state indicator to obtain the systematic component of the agents payoffs. These feed into the simulation of choice probabilities.
  +        k = mapping_state_idx(period + 1, exp_a + 1, exp_b + 1, edu + 1, edu_lagged + 1)
  +        payoffs_systematic = periods_payoffs_systematic(period + 1, k + 1, :)
  +
  +        ! Extract relevant deviates from standard normal distribution.
  +        draws_prob_raw = periods_draws_prob(period + 1, :, :)
  +
  +        ! Prepare to calculate product of likelihood contributions.
  +        crit_val_contrib = 1.0
  +
  +        ! If an agent is observed working, then the the labor market shocks are observed and the conditional distribution is used to determine the choice probabilities.
  +        dist = zero_dble
  +        IF (is_working) THEN
  +
  +            ! Calculate the disturbance, which follows a normal distribution.
  +            dist = clip_value(LOG(data_array(j, 4)), -HUGE_FLOAT, HUGE_FLOAT) - clip_value(LOG(payoffs_systematic(idx)), -HUGE_FLOAT, HUGE_FLOAT) 
  +
  +            ! If there is no random variation in payoffs, then the observed wages need to be identical their systematic components. The discrepancy between the observed wages and their systematic components might be small due to the reading in of the dataset.
  +            IF (is_deterministic) THEN
  +                IF (dist .GT. SMALL_FLOAT) THEN
  +                    rslt = zero_dble
                           RETURN
  -                    END IF
  -                END IF
  -
  -
  -            END IF
  -
  -            ! Simulate the conditional distribution of alternative-specific value functions and determine the choice probabilities.
  -            counts = zero_int
  -            prob_obs = zero_dble
  -
  -            DO s = 1, num_draws_prob
  -                
  -                ! Extract deviates from (un-)conditional normal distributions.
  -                draws_stan = draws_prob_raw(s, :)
  -
  -                ! Construct independent normal draws implied by the agents state experience. This is need to maintain the correlation structure of the disturbances. Special care is needed in case of a deterministic model, as otherwise a zero division error occurs. 
  -                IF (is_working) THEN 
  -                    
  -                    IF (is_deterministic) THEN
  -
  -                        prob_wage = HUGE_FLOAT
  -
  -                    ELSE
  -
  -                        IF (choice == 1) THEN
  -                            draws_stan(idx) = dist / shocks_cholesky(idx, idx)
  -                        ELSE
  -                            draws_stan(idx) = (dist - shocks_cholesky(idx, 1) * draws_stan(1)) / shocks_cholesky(idx, idx)
  -                        END IF
  -
  -                        prob_wage = normal_pdf(draws_stan(idx), zero_dble, one_dble) / SQRT(shocks_cov(idx, idx))
  -                    
  -                    END IF
  -
  -                ELSE 
  -                    prob_wage = one_dble
  -                END IF
  -
  -                ! As deviates are aligned with the state experiences, create the conditional draws. Note, that the realization of the random component of wages align withe their observed counterpart in the data.
  -                draws_cond = MATMUL(draws_stan, TRANSPOSE(shocks_cholesky))
  -
  -                ! Extract deviates from (un-)conditional normal distributions and transform labor market shocks.
  -                draws = draws_cond
  -                draws(1) = clip_value(EXP(draws(1)), zero_dble, HUGE_FLOAT)
  -                draws(2) = clip_value(EXP(draws(2)), zero_dble, HUGE_FLOAT)
  -
  -                ! Calculate total payoff.
  -                CALL get_total_value(total_payoffs, period, payoffs_systematic, draws, mapping_state_idx, periods_emax, k, states_all)
  -
  -                ! Record optimal choices
  -                counts(MAXLOC(total_payoffs)) = counts(MAXLOC(total_payoffs)) + 1
  -
  -                ! Get the smoothed choice probability
  -                prob_choice = get_smoothed_probability(total_payoffs, idx)
  -                prob_obs = prob_obs + prob_choice * prob_wage
  -
  -            END DO
  -
  -            ! Determine relative shares
  -            prob_obs = prob_obs / num_draws_prob
  -
  -            ! If there is no random variation in payoffs, then this implies a unique optimal choice.
  -            IF (is_deterministic) THEN
  -                IF  ((counts(idx) .EQ. num_draws_prob) .EQV. .FALSE.) THEN
  -                    rslt = zero_dble
  -                    RETURN
                   END IF
               END IF
   
  -            ! Adjust  and record likelihood contribution
  -            crit_val(j) = prob_obs
  -
  -            j = j + 1
  +        END IF
  +
  +        ! Simulate the conditional distribution of alternative-specific value functions and determine the choice probabilities.
  +        counts = zero_int
  +        prob_obs = zero_dble
  +
  +        DO s = 1, num_draws_prob
  +                
  +            ! Extract deviates from (un-)conditional normal distributions.
  +            draws_stan = draws_prob_raw(s, :)
  +
  +            ! Construct independent normal draws implied by the agents state experience. This is need to maintain the correlation structure of the disturbances. Special care is needed in case of a deterministic model, as otherwise a zero division error occurs. 
  +            IF (is_working) THEN 
  +                    
  +                IF (is_deterministic) THEN
  +
  +                    prob_wage = HUGE_FLOAT
  +
  +                ELSE
  +
  +                    IF (choice == 1) THEN
  +                        draws_stan(idx) = dist / shocks_cholesky(idx, idx)
  +                    ELSE
  +                        draws_stan(idx) = (dist - shocks_cholesky(idx, 1) * draws_stan(1)) / shocks_cholesky(idx, idx)
  +                    END IF
  +
  +                    prob_wage = normal_pdf(draws_stan(idx), zero_dble, one_dble) / SQRT(shocks_cov(idx, idx))
  +                    
  +                END IF
  +
  +            ELSE 
  +                prob_wage = one_dble
  +            END IF
  +
  +            ! As deviates are aligned with the state experiences, create the conditional draws. Note, that the realization of the random component of wages align withe their observed counterpart in the data.
  +            draws_cond = MATMUL(draws_stan, TRANSPOSE(shocks_cholesky))
  +
  +            ! Extract deviates from (un-)conditional normal distributions and transform labor market shocks.
  +            draws = draws_cond
  +            draws(1) = clip_value(EXP(draws(1)), zero_dble, HUGE_FLOAT)
  +            draws(2) = clip_value(EXP(draws(2)), zero_dble, HUGE_FLOAT)
  +
  +            ! Calculate total payoff.
  +            CALL get_total_value(total_payoffs, period, payoffs_systematic, draws, mapping_state_idx, periods_emax, k, states_all)
  +
  +            ! Record optimal choices
  +            counts(MAXLOC(total_payoffs)) = counts(MAXLOC(total_payoffs)) + 1
  +
  +            ! Get the smoothed choice probability
  +            prob_choice = get_smoothed_probability(total_payoffs, idx)
  +            prob_obs = prob_obs + prob_choice * prob_wage
   
           END DO
  +
  +        ! Determine relative shares
  +        prob_obs = prob_obs / num_draws_prob
  +
  +        ! If there is no random variation in payoffs, then this implies a unique optimal choice.
  +        IF (is_deterministic) THEN
  +            IF  ((counts(idx) .EQ. num_draws_prob) .EQV. .FALSE.) THEN
  +                rslt = zero_dble
  +                RETURN
  +            END IF
  +        END IF
  +
  +        ! Adjust  and record likelihood contribution
  +        crit_val(j) = prob_obs
   
       END DO
   
       ! Scaling
       crit_val = clip_value(LOG(crit_val), -HUGE_FLOAT, HUGE_FLOAT)
  -    rslt = -SUM(crit_val) / (num_agents_est * num_periods)
  +    rslt = -SUM(crit_val) / (num_obs * num_periods)
   
       ! If there is no random variation in payoffs and no agent violated the implications of observed wages and choices, then the evaluation return a value of one.
       IF (is_deterministic) THEN
  
positive.code: 
      REAL(our_dble), INTENT(IN)      :: periods_payoffs_systematic(num_periods, max_states_period, 4)
      REAL(our_dble), INTENT(IN)      :: periods_draws_prob(num_periods, num_draws_prob, 4)
      REAL(our_dble), INTENT(IN)      :: data_array(num_agents_est * num_periods, 8)
      REAL(our_dble), INTENT(IN)      :: shocks_cholesky(4, 4)
      REAL(our_dble), INTENT(IN)      :: periods_emax(num_periods, max_states_period)
  
      INTEGER(our_int), INTENT(IN)    :: mapping_state_idx(num_periods, num_periods, num_periods, min_idx, 2)
      INTEGER(our_int), INTENT(IN)    :: states_all(num_periods, max_states_period, 4)
  
      !/* internal objects        */
  
      REAL(our_dble)                  :: crit_val(num_agents_est * num_periods)
      REAL(our_dble)                  :: draws_prob_raw(num_draws_prob, 4)
      REAL(our_dble)                  :: payoffs_systematic(4)
      REAL(our_dble)                  :: crit_val_contrib
  
negative.code: !******************************************************************************
  !******************************************************************************
  SUBROUTINE fort_evaluate(rslt, periods_payoffs_systematic, mapping_state_idx, periods_emax, states_all, shocks_cholesky, data_array, periods_draws_prob)
  
      ! DEVELOPMENT NOTES
      !   
      !   This has assumed-shape arguments as it is used when the work is distributed across the slaves as well.
      !
  
      !/* external objects        */
  
  