bayesTFR: Bayesian Fertility Projection

get.obs.estimate.diff <- function(mcmc)
{
  mcmc$eps_unc <- list()
  for(country in 1:min(mcmc$meta$nr_countries, length(mcmc$meta$raw.data)))
  {
    if (is.null(mcmc$meta$raw.data[[country]])) next
    tmp.year <- mcmc$meta$raw.data[[country]]$year
    tmp.tfr <- mcmc$meta$raw.data[[country]]$tfr
    if (mcmc$meta$annual.simulation)
    {
      tmp.ind <- (tmp.year - mcmc$meta$start.year + 1)
      ref.tfr <- mcmc$meta$tfr_matrix_all[, country]
      interpolated.tfr <- ref.tfr[round(tmp.ind + 0.01)]
    }
    else
    {
      tmp.ind <- (tmp.year - mcmc$meta$start.year + 2) / 5
      left.distance <- tmp.ind - floor(tmp.ind)
      ref.tfr <- mcmc$meta$tfr_matrix_all[, country]
      interpolated.tfr <- ref.tfr[pmax(floor(tmp.ind), 1)] * (5 - left.distance) / 5 + 
        ref.tfr[pmin(pmax(floor(tmp.ind), 1) + 1, length(ref.tfr))] * left.distance / 5
    }
    delta <- tmp.tfr - interpolated.tfr
    mcmc$eps_unc[[country]] <- delta
  }
  return(mcmc)
}

get.obs.estimate.diff.original <- function(mcmc)
{
  mcmc$meta$raw_data.original$eps <- NA
  if (mcmc$meta$annual.simulation)
  {
    for(year in mcmc$meta$start.year:mcmc$meta$present.year)
    {
      year.ind <- year - mcmc$meta$start.year + 1
      if (length(mcmc$meta$ind.by.year[[year.ind]]) == 0) next
      ref.tfr <- mcmc$meta$tfr_matrix_all[year.ind, mcmc$meta$country.ind.by.year[[year.ind]]]
      idx <- mcmc$meta$ind.by.year[[year.ind]]
      mcmc$meta$raw_data.original$eps[idx] <- mcmc$meta$raw_data.original$tfr[idx] - ref.tfr
    }
  }
  else 
  {
    years <- as.numeric(rownames(mcmc$meta$tfr_matrix_all))
    count <- 0
    for (year in c(years, rev(years)[1]+5))
    {
      count <- count + 1
      if (length(mcmc$meta$ind.by.year[[count]]) == 0) next
      idx <- mcmc$meta$ind.by.year[[count]]
      tfr.left <- mcmc$meta$tfr_matrix_all[max(count-1, 1), mcmc$meta$country.ind.by.year[[count]]]
      tfr.right <- mcmc$meta$tfr_matrix_all[min(count, nrow(mcmc$meta$tfr_matrix_all)), mcmc$meta$country.ind.by.year[[count]]]
      left.distance <- mcmc$meta$left.distance[[count]]
      ref.tfr <- (tfr.left * (5 - left.distance) + tfr.right * left.distance) / 5
      mcmc$meta$raw_data.original$eps[idx] <- mcmc$meta$raw_data.original$tfr[idx] - ref.tfr
    }
  }
  return(mcmc)
}

update.obs.estimate.diff <- function(mcmc, country, year)
{
  if (mcmc$meta$annual.simulation)
  {
    indices <- which(round(mcmc$meta$raw_data_extra[[country]]$year - mcmc$meta$start.year + 1.01) == year)
  }
  else
  {
    indices <- which((mcmc$meta$raw_data_extra[[country]]$year < (year * 5 + mcmc$meta$start.year + 3)) & 
                       (mcmc$meta$raw_data_extra[[country]]$year > (year * 5 + mcmc$meta$start.year -7)))
  }
  if (length(indices) > 0)
  {
    tmp.tfr <- mcmc$meta$raw_data_extra[[country]]$tfr[indices]
    tmp.year <- mcmc$meta$raw_data_extra[[country]]$year[indices]
    if (mcmc$meta$annual.simulation)
    {
      tmp.ind <- (tmp.year - mcmc$meta$start.year + 1)
      ref.tfr <- mcmc$meta$tfr_all[, country]
      interpolated.tfr <- ref.tfr[round(tmp.ind +0.01)]
    }
    else
    {
      tmp.ind <- (tmp.year - mcmc$meta$start.year + 2) / 5
      left.distance <- tmp.ind - floor(tmp.ind)
      ref.tfr <- mcmc$meta$tfr_all[, country]
      interpolated.tfr <- ref.tfr[pmax(floor(tmp.ind), 1)] * (5 - left.distance) / 5 + 
        ref.tfr[pmin(floor(tmp.ind) + 1, length(ref.tfr))] * left.distance / 5
    }
    return (tmp.tfr - interpolated.tfr)
  }
  else
  {
    return (numeric(0))
  }
}

estimate.bias.sd.raw <- function(mcmc)
{
  for(country in 1:min(mcmc$meta$nr_countries, length(mcmc$meta$raw.data)))
  {
    if (is.null(mcmc$meta$raw.data[[country]])) next
    source <- factor(as.character(mcmc$meta$raw.data[[country]]$source))
    estimation.method <- factor(as.character(mcmc$meta$raw.data[[country]]$method))
    diff <- mcmc$eps_unc[[country]]
    regressor <- '1'
    if (length(levels(source)) > 1)
    {
      regressor <- paste0(regressor, ' + source')
    }
    if (length(levels(estimation.method)) > 1)
    {
      regressor <- paste0(regressor, ' + estimation.method')
    }
    m1 <- lm(as.formula(paste0('diff ~ ', regressor)))
    bias <- predict(m1)
    abs.residual <- abs(residuals(m1))
    m2 <- lm(as.formula(paste0('abs.residual ~ ', regressor)))
    std <- predict(m2) * sqrt(pi/2)
    std[std < 1e-6] <- 0.1
    std[std < (abs(bias) / 2)] <- abs(bias[std < (abs(bias) / 2)])
    mcmc$meta$raw.data[[country]]$bias <- bias
    mcmc$meta$raw.data[[country]]$std <- std
  }
  
  return(mcmc)
}

estimate.bias.sd.original <- function(mcmc, iso.unbiased=NULL, covariates=c('source', 'method'), 
                                      cont_covariates=NULL, source.col.name="source", countries = NULL)
{
  mcmc$meta$raw_data.original$bias <- NA
  mcmc$meta$raw_data.original$std <- NA
  if (is.null(mcmc$meta$bias_model))
  {
    mcmc$meta$bias_model <- list()
    mcmc$meta$std_model <- list()
  }
  if(is.null(countries)) countries <- 1:mcmc$meta$nr_countries
  for(country in countries)
  {
    ISO.code <- get.country.object(country, meta=mcmc$meta, index=TRUE)
    if (is.na(ISO.code$code)) next
    index.by.country <- which(mcmc$meta$raw_data.original$country_code == ISO.code$code)
    if (length(index.by.country) == 0) next
    regressor <- '1'
    diff <- mcmc$meta$raw_data.original$eps[index.by.country]
    if (length(covariates) > 0)
    {
      for (i in 1:length(covariates))
      {
        covariate <- covariates[i]
        if(!covariate %in% colnames(mcmc$meta$raw_data.original)){
            warning("Covariate ", covariate, " not available in the data.")
            next
        }
        assign(paste0('covariate_', i), factor(as.character(mcmc$meta$raw_data.original[index.by.country, covariate])))
        if (length(levels(get(paste0('covariate_', i)))) > 1)
          regressor <- paste0(regressor, ' + covariate_', i)
      }
    }
    if (length(cont_covariates) > 0)
    {
      for (i in 1:length(cont_covariates))
      {
        covariate <- cont_covariates[i]
        if(!covariate %in% colnames(mcmc$meta$raw_data.original)){
            warning("Cont. covariate ", covariate, " not available in the data.")
            next
        }
        assign(paste0('cont_covariate_', i), mcmc$meta$raw_data.original[index.by.country, covariate])
        if (max(get(paste0('cont_covariate_', i))) - min(get(paste0('cont_covariate_', i))) > 1e-6)
        {
          regressor <- paste0(regressor, ' + cont_covariate_', i)
        }
      }
    }
    
    m1 <- lm(as.formula(paste0('diff ~ ', regressor)), model = FALSE)
    bias <- predict(m1)
    abs.residual <- abs(residuals(m1))
    m2 <- lm(as.formula(paste0('abs.residual ~ ', regressor)), model = FALSE)
    std <- predict(m2) * sqrt(pi/2)
    std[std < 1e-6] <- 0.1
    std[std < (abs(bias) / 2)] <- abs(bias[std < (abs(bias) / 2)]) 
    mcmc$meta$raw_data.original$bias[index.by.country] <- bias
    mcmc$meta$raw_data.original$std[index.by.country] <- std
    ## Optional
    if (ISO.code$code %in% iso.unbiased)
    {
      if (source.col.name %in% covariates)
      {
        index.by.country.vr.estimate <- which((mcmc$meta$raw_data.original$country_code == ISO.code$code) & 
                                                (mcmc$meta$raw_data.original[, source.col.name] %in% c("VR")))
        mcmc$meta$raw_data.original$bias[index.by.country.vr.estimate] <- 0
        mcmc$meta$raw_data.original$std[index.by.country.vr.estimate] <- 0.016
      }
      else
      {
        warning("Covariate source not found. iso.unbiased is not used.\n")
      }
    }
    # to save space
    attr(m1$terms, ".Environment") <- NULL
    attr(m2$terms, ".Environment") <- NULL
    #m1$fitted.values <- NULL # cannot be deleted as the summary function would not work
    #m2$fitted.values <- NULL
    
    mcmc$meta$bias_model[[country]] <- m1
    mcmc$meta$std_model[[country]] <- m2
  }
  return(mcmc)
}

get.eps.all.phases <- function(Dlpar, mcmc, country)
{
  eps_return <- numeric(length = mcmc$meta$T_end-1)
  if (mcmc$meta$lambda_c[country] > mcmc$meta$start_c[country]){
    id2 <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country] - 1)
    dl <- - Dlpar[5]/(1 + exp(- 2*log(9)/Dlpar[1] *(mcmc$meta$tfr_all[id2, country] - Dlpar[1] - Dlpar[2] -Dlpar[3] - Dlpar[4] + 0.5*Dlpar[1]))) + 
      Dlpar[5]/(1 + exp(- 2*log(9)/Dlpar[3] *(mcmc$meta$tfr_all[id2, country] - Dlpar[4] - 0.5*Dlpar[3])))
    dl <- dl * ifelse(mcmc$meta$annual.simulation, 1, 5)
    dl[mcmc$meta$tfr_all[id2, country] < 1] <- 0
    dl[dl<0] <- 0
    eps_return[id2] <- mcmc$meta$tfr_all[id2 + 1, country] - mcmc$meta$tfr_all[id2, country] + dl
  }
  if (mcmc$meta$start_c[country] > 1)
  {
    eps_return[1:(mcmc$meta$start_c[country] - 1)] <- mcmc$meta$tfr_all[2:mcmc$meta$start_c[country], country] - 
      mcmc$meta$tfr_all[1:(mcmc$meta$start_c[country]-1), country]
  }
  if (mcmc$meta$lambda_c[country] < mcmc$meta$T_end_c[country])
  {
    id3 <- which(mcmc$meta$id_phase3 == country)
    eps_return[mcmc$meta$lambda_c[country]:(mcmc$meta$T_end_c[country] - 1)] <- mcmc$meta$tfr_all[(mcmc$meta$lambda_c[country]+1):mcmc$meta$T_end_c[country], country] - 
      mcmc$meta$tfr_all[mcmc$meta$lambda_c[country]:(mcmc$meta$T_end_c[country]-1), country] * mcmc$rho.c[id3] - (1-mcmc$rho.c[id3]) * mcmc$mu.c[id3]
  }
  return (eps_return)
}

get.sd.all.phases <- function(mcmc, country)
{
  sd_return <- mcmc$sd_Tc[, country]
  if (mcmc$meta$start_c[country] > 1)
  {
    sd_return[1:(mcmc$meta$start_c[country] - 1)] <- mcmc$sd_eps_tau
  }
  if (mcmc$meta$lambda_c[country] < mcmc$meta$T_end_c[country])
  {
    sd_return[mcmc$meta$lambda_c[country]:(mcmc$meta$T_end_c[country] - 1)] <- mcmc$sigma.eps
  }
  return (sd_return)
}

mcmc.update.tfr <- function(country, mcmc)
{
  Dlpar <- c((mcmc$U_c[country]-mcmc$Triangle_c4[country])*
               exp(mcmc$gamma_ci[country,])/
               sum(exp(mcmc$gamma_ci[country,])), 
             mcmc$Triangle_c4[country], 
             mcmc$d_c[country])
  if (mcmc$meta$lambda_c[country] > mcmc$meta$start_c[country]){
    epsT.idx <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country]-1)
  }
  else{
    epsT.idx <- c()
  }
  eps_tfr_prev <- get.eps.all.phases(Dlpar, mcmc, country)
  sd_tfr_prev <- get.sd.all.phases(mcmc, country)
  eps_tfr_prop <- eps_tfr_prev
  sd_tfr_prop <- sd_tfr_prev
  if (country %in% mcmc$meta$id_phase3) {id3 <- which(mcmc$meta$id_phase3 == country)}
  for (year in 1:mcmc$meta$T_end_c[country])
  {
    tmp.years <- max(1, year-1):min(year, mcmc$meta$T_end_c[country] -1)
    if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      tmp.extra.year <- max(2, year-1):min(year+1, mcmc$meta$T_end_c[country] -1)
      tmp.extra.year <- tmp.extra.year[tmp.extra.year > mcmc$meta$start_c[country] & tmp.extra.year < mcmc$meta$lambda_c[country]]
      tmp.more.years <- tmp.years[!(tmp.years %in% tmp.extra.year)]
    }
    sd_f <- ifelse(mcmc$finished.iter>=100, 2.4 * max(mcmc$meta$tfr_sd_all[year, country], 0.01), 0.2)
    prop_tfr <- rnorm(1, mcmc$meta$tfr_all[year, country], sd_f)
    while (prop_tfr < 0.5) prop_tfr <- rnorm(1, mcmc$meta$tfr_all[year, country], sd_f)
    if (mcmc$meta$annual.simulation)
    {
      indices <- which(round(mcmc$meta$raw_data_extra[[country]]$year - mcmc$meta$start.year + 1.01) == year)
    }
    else
    {
      indices <- which((mcmc$meta$raw_data_extra[[country]]$year < (year * 5 + mcmc$meta$start.year + 3)) & 
                         (mcmc$meta$raw_data_extra[[country]]$year > (year * 5 + mcmc$meta$start.year - 7)))
    }
    tfr_old <- mcmc$meta$tfr_all[year, country]
    eps_prev <- mcmc$eps_unc[[country]][indices]
    loglik_prev <- sum(dnorm(eps_prev, mean=mcmc$meta$raw_data_extra[[country]]$bias[indices], sd=mcmc$meta$raw_data_extra[[country]]$std[indices], log=TRUE))
    if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      loglik_prev <- loglik_prev + sum(dnorm(eps_tfr_prev[tmp.extra.year], mean = mcmc$rho.phase2 * eps_tfr_prev[tmp.extra.year - 1], 
                                             sd = sd_tfr_prev[tmp.extra.year], log=TRUE))
      loglik_prev <- loglik_prev + sum(dnorm(eps_tfr_prev[tmp.more.years], mean = 0, sd = sd_tfr_prev[tmp.more.years], log=TRUE))
    }
    else loglik_prev <- loglik_prev + sum(dnorm(eps_tfr_prev[tmp.years], mean = 0, sd = sd_tfr_prev[tmp.years], log=TRUE))
    mcmc$meta$tfr_all[year, country] <- prop_tfr
    eps_new <- update.obs.estimate.diff(mcmc, country, year)
    loglik_new <- sum(dnorm(eps_new, mean=mcmc$meta$raw_data_extra[[country]]$bias[indices], sd=mcmc$meta$raw_data_extra[[country]]$std[indices], log=TRUE))
    if (year < mcmc$meta$start_c[country])
    {
      eps_tfr_prop[tmp.years] <- mcmc$meta$tfr_all[tmp.years + 1, country] - mcmc$meta$tfr_all[tmp.years, country]
    }
    else if (year > mcmc$meta$lambda_c[country])
    {
      eps_tfr_prop[tmp.years] <- mcmc$meta$tfr_all[tmp.years + 1, country] - mcmc$meta$tfr_all[tmp.years, country] * 
        mcmc$rho.c[id3] - (1-mcmc$rho.c[id3]) * mcmc$mu.c[id3]
    }
    else 
      eps_tfr_prop <- get.eps.all.phases(Dlpar, mcmc, country)
    
    if ((year >= mcmc$meta$start_c[country]) && (year < mcmc$meta$lambda_c[country]))
    {
      tmp <- (prop_tfr - mcmc$S_sd) * ifelse(prop_tfr > mcmc$S_sd, -mcmc$a_sd, mcmc$b_sd)
      sd_tfr_prop[year] <- ifelse(mcmc$const_sd_dummie_Tc[year, country] == 1, mcmc$const_sd, 1) * 
        ifelse((mcmc$sigma0 + tmp > mcmc$meta$sigma0.min), mcmc$sigma0 + tmp, mcmc$meta$sigma0.min)
    }
    if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      loglik_new <- loglik_new + sum(dnorm(eps_tfr_prop[tmp.extra.year], mean = mcmc$rho.phase2 * eps_tfr_prop[tmp.extra.year - 1], 
                                             sd = sd_tfr_prop[tmp.extra.year], log=TRUE))
      loglik_new <- loglik_new + sum(dnorm(eps_tfr_prop[tmp.more.years], mean = 0, sd = sd_tfr_prop[tmp.more.years], log=TRUE))
    }
    else loglik_new <- loglik_new + sum(dnorm(eps_tfr_prop[tmp.years], mean = 0, sd = sd_tfr_prop[tmp.years], log=TRUE))
    
    # loglik_new <- loglik_new + sum(dnorm(eps_tfr_prop[tmp.years], mean=0, sd=sd_tfr_prop[tmp.years], log=TRUE))
    accept.prob <- exp(loglik_new - loglik_prev)
    if (runif(1) < accept.prob)
    {
      if (year %in% epsT.idx)
      {
        if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2 && year > mcmc$meta$start_c[country]) 
          mcmc$eps_Tc[year, country] <- eps_tfr_prop[year] - mcmc$rho.phase2 * eps_tfr_prop[year - 1]
        else mcmc$eps_Tc[year, country] <- eps_tfr_prop[year]
        mcmc$sd_Tc[year, country] <- sd_tfr_prop[year]
        mcmc$add_to_sd_Tc[year, country] <- (prop_tfr - mcmc$S_sd) * ifelse(prop_tfr > mcmc$S_sd, -mcmc$a_sd, mcmc$b_sd)
      }
      if ((year - 1) %in% epsT.idx)
      {
        if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2 && year > mcmc$meta$start_c[country] + 1) 
          mcmc$eps_Tc[year - 1, country] <- eps_tfr_prop[year - 1] - mcmc$rho.phase2 * eps_tfr_prop[year - 2]
        else mcmc$eps_Tc[year - 1, country] <- eps_tfr_prop[year - 1]
        mcmc$sd_Tc[year - 1, country] <- sd_tfr_prop[year - 1]
      }
      if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2 && (year + 1) %in% epsT.idx && year > mcmc$meta$start_c[country] - 1)
        mcmc$eps_Tc[year + 1, country] <- eps_tfr_prop[year + 1] - mcmc$rho.phase2 * eps_tfr_prop[year]
      
      mcmc$eps_unc[[country]][indices] <- eps_new
      eps_tfr_prev <- eps_tfr_prop
      sd_tfr_prev <- sd_tfr_prop
      
    }
    else
    {
      mcmc$meta$tfr_all[year, country] <- tfr_old
      sd_tfr_prop <- sd_tfr_prev
    }
    mu0 <- mcmc$meta$tfr_mu_all[year, country]
    sigma0 <- mcmc$meta$tfr_sd_all[year, country]
    mcmc$meta$tfr_mu_all[year, country] <- (mcmc$finished.iter * mu0 + mcmc$meta$tfr_all[year, country]) / (mcmc$finished.iter + 1)
    mu <- mcmc$meta$tfr_mu_all[year, country]
    mcmc$meta$tfr_sd_all[year, country] <- sqrt(((mcmc$finished.iter-1) * sigma0 ** 2 + 
                                                  mcmc$finished.iter * (mu0 - mu) ** 2 + 
                                                  (mcmc$meta$tfr_all[year, country] - mu) ** 2)/mcmc$finished.iter)
    
  }
  if (country %in% mcmc$meta3$id_phase3) mcmc$observations[[id3]] <- mcmc$meta$tfr_all[mcmc$meta$lambda_c[country]:mcmc$meta$T_end_c[country], country]
  if (country %in% mcmc$meta$id_DL)
    mcmc$data.list[[country]][epsT.idx] <- mcmc$meta$tfr_all[epsT.idx, country]
}

# get.log.lik.year <- function(year.ind, mcmc, Dlpar, phase3par, id_phase1, id_phase2, id_phase3, tfr=NULL, prev=TRUE)
# {
# ## Work for computing log-likelihood from year.ind -> year.ind+1
#   if (is.null(tfr))
#   {
#     tfr_0 <- mcmc$meta$tfr_all[year.ind, ]
#     tfr_1 <- mcmc$meta$tfr_all[year.ind + 1, ]
#   }
#   else if (prev)
#   {
#     tfr_0 <- tfr
#     tfr_1 <- mcmc$meta$tfr_all[year.ind + 1, ]
#   }
#   else
#   {
#     tfr_0 <- mcmc$meta$tfr_all[year.ind, ]
#     tfr_1 <- tfr
#   }
#   
#   eps <- numeric(mcmc$meta$nr_countries)
#   eps[id_phase1] <- tfr_1[id_phase1] - tfr_0[id_phase1]
#   if (!is.null(tfr))
#   {
#     dl <- - Dlpar[id_phase2,5]/(1 + exp(- 2*log(9)/Dlpar[id_phase2, 1] * 
#                                           (tfr_0[id_phase2] - 0.5 * Dlpar[id_phase2, 1] - Dlpar[id_phase2, 2] -
#                                              Dlpar[id_phase2, 3] - Dlpar[id_phase2, 4]))) + 
#       Dlpar[id_phase2, 5]/(1 + exp(- 2*log(9)/Dlpar[id_phase2, 3] *
#                                      (tfr_0[id_phase2] - Dlpar[id_phase2, 4] - 0.5*Dlpar[id_phase2,3])))
#     dl <- dl * ifelse(mcmc$meta$annual.simulation, 1, 5)
#     dl[tfr_0[id_phase2] < 1] <- 0
#     dl[dl<0] <- 0
#     eps[id_phase2] <- tfr_1[id_phase2] - (tfr_0[id_phase2] - dl)
#   }
#   else
#   {
#     eps[id_phase2] <- mcmc$eps_Tc[year.ind, id_phase2]
#   }
#   eps[id_phase3] <- tfr_1[id_phase3] - phase3par[id_phase3, 1] - 
#     phase3par[id_phase3, 2] * (tfr_0[id_phase3] - phase3par[id_phase3, 1])
#   
#   std <- numeric(mcmc$meta$nr_countries)
#   std[id_phase1] <- mcmc$sd_eps_tau
#   std[id_phase3] <- mcmc$sigma.eps
#   if (!is.null(tfr))
#   {
#     const_dummy <- ifelse(mcmc$const_sd_dummie_Tc[year.ind, 1], mcmc$const_sd, 1)
#     tmp_add <- tfr_0[id_phase2] - mcmc$S_sd
#     tmp_add[tfr_0[id_phase2] > mcmc$S_sd] <- - tmp_add[tfr_0[id_phase2] > mcmc$S_sd] * mcmc$a_sd
#     tmp_add[tfr_0[id_phase2] <= mcmc$S_sd] <- tmp_add[tfr_0[id_phase2] <= mcmc$S_sd] * mcmc$b_sd
#     std[id_phase2] <- tmp_add + mcmc$sigma0
#     std[std <= mcmc$meta$sigma0.min] <- mcmc$meta$sigma0.min
#   }
#   else std[id_phase2] <- mcmc$sd_Tc[year.ind, id_phase2]
#   
#   log.lik <- dnorm(eps, mean=0, sd=std, log=TRUE)
#   output <- list(log.lik = log.lik)
#   if (!is.null(tfr))
#   {
#     output$eps <- eps
#     output$std <- std
#   }
#   return (output)
# }

get.log.lik.year <- function(year.ind, mcmc, Dlpar, phase3par, tfr=NULL, prev=TRUE, prev_two=FALSE)
{
  ## Work for computing log-likelihood from year.ind -> year.ind+1
  id_phase1 <- mcmc$meta$id_phase1_by_year[[year.ind]]
  id_phase2 <- mcmc$meta$id_phase2_by_year[[year.ind]]
  id_phase3 <- mcmc$meta$id_phase3_by_year[[year.ind]]
  if (year.ind > 1) id_phase2_prev <- mcmc$meta$id_phase2_by_year[[year.ind-1]]
  
  if (is.null(tfr) || prev_two)
  {
    tfr_0 <- mcmc$meta$tfr_all[year.ind, ]
    tfr_1 <- mcmc$meta$tfr_all[year.ind + 1, ]
  }
  else if (prev)
  {
    tfr_0 <- tfr
    tfr_1 <- mcmc$meta$tfr_all[year.ind + 1, ]
  }
  else
  {
    tfr_0 <- mcmc$meta$tfr_all[year.ind, ]
    tfr_1 <- tfr
  }
  
  eps <- numeric(mcmc$meta$nr_countries)
  eps[id_phase1] <- tfr_1[id_phase1] - tfr_0[id_phase1]
  if (!is.null(tfr))
  {
    dl <- - Dlpar[id_phase2,5]/(1 + exp(- 2*log(9)/Dlpar[id_phase2, 1] * 
                                          (tfr_0[id_phase2] - 0.5 * Dlpar[id_phase2, 1] - Dlpar[id_phase2, 2] -
                                             Dlpar[id_phase2, 3] - Dlpar[id_phase2, 4]))) + 
      Dlpar[id_phase2, 5]/(1 + exp(- 2*log(9)/Dlpar[id_phase2, 3] *
                                     (tfr_0[id_phase2] - Dlpar[id_phase2, 4] - 0.5*Dlpar[id_phase2,3])))
    dl <- dl * ifelse(mcmc$meta$annual.simulation, 1, 5)
    dl[tfr_0[id_phase2] < 1] <- 0
    dl[dl<0] <- 0
    eps[id_phase2] <- tfr_1[id_phase2] - (tfr_0[id_phase2] - dl)
    if ((year.ind > 1) && !(is.null(mcmc$meta$ar.phase2)) && mcmc$meta$ar.phase2)
    {
      joint_phase2 <- intersect(id_phase2, id_phase2_prev)
      if (prev_two) tfr_minus_1 <- tfr
      else tfr_minus_1 <- mcmc$meta$tfr_all[year.ind - 1, ]
      dl_prev <- - Dlpar[joint_phase2,5]/(1 + exp(- 2*log(9)/Dlpar[joint_phase2, 1] * 
                                            (tfr_minus_1[joint_phase2] - 0.5 * Dlpar[joint_phase2, 1] - Dlpar[joint_phase2, 2] -
                                               Dlpar[joint_phase2, 3] - Dlpar[joint_phase2, 4]))) + 
        Dlpar[joint_phase2, 5]/(1 + exp(- 2*log(9)/Dlpar[joint_phase2, 3] *
                                       (tfr_minus_1[joint_phase2] - Dlpar[joint_phase2, 4] - 0.5*Dlpar[joint_phase2,3])))
      dl_prev <- dl_prev * ifelse(mcmc$meta$annual.simulation, 1, 5)
      dl_prev[tfr_minus_1[joint_phase2] < 1] <- 0
      dl_prev[dl_prev<0] <- 0
      eps_prev <- tfr_0[joint_phase2] - (tfr_minus_1[joint_phase2] - dl_prev)
      eps[joint_phase2] <- eps[joint_phase2] - mcmc$rho.phase2 * eps_prev
    }
  }
  else
  {
    eps[id_phase2] <- mcmc$eps_Tc[year.ind, id_phase2]
  }
  eps[id_phase3] <- tfr_1[id_phase3] - phase3par[id_phase3, 1] - 
    phase3par[id_phase3, 2] * (tfr_0[id_phase3] - phase3par[id_phase3, 1])
  
  std <- numeric(mcmc$meta$nr_countries)
  std[id_phase1] <- mcmc$sd_eps_tau
  std[id_phase3] <- mcmc$sigma.eps
  if (!is.null(tfr))
  {
    const_dummy <- ifelse(mcmc$const_sd_dummie_Tc[year.ind, 1], mcmc$const_sd, 1)
    tmp_add <- tfr_0[id_phase2] - mcmc$S_sd
    tmp_add[tfr_0[id_phase2] > mcmc$S_sd] <- - tmp_add[tfr_0[id_phase2] > mcmc$S_sd] * mcmc$a_sd
    tmp_add[tfr_0[id_phase2] <= mcmc$S_sd] <- tmp_add[tfr_0[id_phase2] <= mcmc$S_sd] * mcmc$b_sd
    std[id_phase2] <- tmp_add + mcmc$sigma0
    std[std <= mcmc$meta$sigma0.min] <- mcmc$meta$sigma0.min
  }
  else std[id_phase2] <- mcmc$sd_Tc[year.ind, id_phase2]
  
  log.lik <- dnorm(eps, mean=0, sd=std, log=TRUE)
  output <- list(log.lik = log.lik)
  if (!is.null(tfr))
  {
    output$eps <- eps
    output$std <- std
  }
  return (output)
}

get.log.lik.raw <- function(year.ind, mcmc, tfr)
{
  ## Work for computing log-likelihood from year.ind -> year.ind+1
  if ((length(mcmc$meta$ind.by.year[[year.ind]]) == 0) && mcmc$meta$annual.simulation) return (numeric(mcmc$meta$nr_countries))
  if (!mcmc$meta$annual.simulation && (length(mcmc$meta$ind.by.year[[year.ind]]) == 0 && length(mcmc$meta$ind.by.year[[year.ind + 1]]) == 0)) 
    return (numeric(mcmc$meta$nr_countries))
  output.loglik <- numeric(mcmc$meta$nr_countries)
  
  if (length(mcmc$meta$ind.by.year[[year.ind]]) > 0)
  {
    bias <- mcmc$meta$raw_data.original$bias[mcmc$meta$ind.by.year[[year.ind]]]
    std <- mcmc$meta$raw_data.original$std[mcmc$meta$ind.by.year[[year.ind]]]
    if (mcmc$meta$annual.simulation) tfr.ref <- tfr[mcmc$meta$country.ind.by.year[[year.ind]]]
    else
    {
      tfr.ref <- tfr[mcmc$meta$country.ind.by.year[[year.ind]]] * mcmc$meta$left.distance[[year.ind]] / 5
      tfr.ref <- tfr.ref + mcmc$meta$tfr_matrix_all[max(1, year.ind - 1), mcmc$meta$country.ind.by.year[[year.ind]]] * 
        (5 - mcmc$meta$left.distance[[year.ind]]) / 5
    }
    delta <- mcmc$meta$raw_data.original$tfr[mcmc$meta$ind.by.year[[year.ind]]] - tfr.ref
    loglik.delta <- dnorm(delta, mean=bias, sd=std, log=TRUE)
    df <- data.frame(idx = mcmc$meta$country.ind.by.year[[year.ind]], loglik = loglik.delta)
    res <- aggregate(df$loglik, by=list(ind=df$idx), sum)
    output.loglik[res$ind] <- res$x
  }
  if (!mcmc$meta$annual.simulation && length(mcmc$meta$ind.by.year[[year.ind + 1]]) > 0)
  {
    bias <- mcmc$meta$raw_data.original$bias[mcmc$meta$ind.by.year[[year.ind + 1]]]
    std <- mcmc$meta$raw_data.original$std[mcmc$meta$ind.by.year[[year.ind + 1]]]
    tfr.ref <- tfr[mcmc$meta$country.ind.by.year[[year.ind + 1]]] * (5 - mcmc$meta$left.distance[[year.ind + 1]]) / 5
    tfr.ref <- tfr.ref + mcmc$meta$tfr_matrix_all[min(mcmc$meta$T_end, year.ind + 1), mcmc$meta$country.ind.by.year[[year.ind + 1]]] * 
      mcmc$meta$left.distance[[year.ind + 1]] / 5
    delta <- mcmc$meta$raw_data.original$tfr[mcmc$meta$ind.by.year[[year.ind + 1]]] - tfr.ref
    loglik.delta <- dnorm(delta, mean=bias, sd=std, log=TRUE)
    df <- data.frame(idx = mcmc$meta$country.ind.by.year[[year.ind + 1]], loglik = loglik.delta)
    res <- aggregate(df$loglik, by=list(ind=df$idx), sum)
    output.loglik[res$ind] <- output.loglik[res$ind] + res$x
  }
  
  return (output.loglik)
}

mcmc.update.tfr.year <- function(mcmc, countries = NULL)
{
  Dlpar <- cbind((mcmc$U_c - mcmc$Triangle_c4) * 
                   exp(mcmc$gamma_ci)/ apply(exp(mcmc$gamma_ci), 1, sum), 
                 mcmc$Triangle_c4, mcmc$d_c)
  nr_countries <- mcmc$meta$nr_countries
  phase3par <- matrix(nrow = nr_countries, ncol=2)
  phase3par[mcmc$meta$id_phase3,] <- cbind(mcmc$mu.c, mcmc$rho.c)
  ## _r represents phases from year -> year + 1
  id_phase1_r <- mcmc$meta$id_phase1_by_year[[1]]
  id_phase3_r <- mcmc$meta$id_phase3_by_year[[1]]
  id_phase2_r <- mcmc$meta$id_phase2_by_year[[1]]
  if (!is.null(countries))
  {
    id_phase1_r <- intersect(id_phase1_r, countries)
    id_phase3_r <- intersect(id_phase3_r, countries)
    id_phase2_r <- intersect(id_phase2_r, countries)
  }
  indices_outliers_m1 <- mcmc$meta$yearly.outliers[[as.character(1)]]
  if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
  {
    indices_outliers_0 <- sort(unique(c(mcmc$meta$yearly.outliers[[as.character(1)]], 
                                        mcmc$meta$yearly.outliers[[as.character(2)]])))
    
    indices_outliers_p1 <- sort(unique(c(mcmc$meta$yearly.outliers[[as.character(1)]], 
                                         mcmc$meta$yearly.outliers[[as.character(2)]])))
  }
  else
  {
    indices_outliers_0 <- mcmc$meta$yearly.outliers[[as.character(1)]]
    indices_outliers_p1 <- mcmc$meta$yearly.outliers[[as.character(2)]]
  }
  for (year in 1:mcmc$meta$T_end)
  {
    if (year > 1)
    {
      id_phase1_l <- id_phase1_r
      id_phase2_l <- id_phase2_r
      id_phase3_l <- id_phase3_r
      id_phase1_r <- mcmc$meta$id_phase1_by_year[[year]]
      id_phase3_r <- mcmc$meta$id_phase3_by_year[[year]]
      id_phase2_r <- mcmc$meta$id_phase2_by_year[[year]]
      if (!is.null(countries))
      {
        id_phase1_r <- intersect(id_phase1_r, countries)
        id_phase3_r <- intersect(id_phase3_r, countries)
        id_phase2_r <- intersect(id_phase2_r, countries)
      }
    }
    loglik_orig <- numeric(nr_countries)
    loglik_proposed <- numeric(nr_countries)
    tfr_orig <- mcmc$meta$tfr_all[year, ]
    if (mcmc$finished.iter <= 100) sd_f <- rep(0.2, nr_countries) else sd_f <- pmax(2.4 * mcmc$meta$tfr_sd_all[year, ], 0.01)
    tfr_proposed <- rep(1, nr_countries)
    if (is.null(countries))
      tfr_proposed <- rnorm(nr_countries, tfr_orig, sd_f)
    else
      tfr_proposed[countries] <- rnorm(length(countries), tfr_orig[countries], sd_f[countries])
    while(any(tfr_proposed < 0.5)) 
    {
      idx <- which(tfr_proposed < 0.5)
      if (mcmc$finished.iter <= 100) sd_f_std <- 0.2 else sd_f_std <- sd_f[idx]
      tfr_proposed[idx] <- rnorm(length(idx), tfr_orig[idx], sd_f_std)
    }
    
    ## Compute Log likelihood
    ## For outliers keep log likelihood to be frozen to 0 for both 
    ## loglik_orig annd loglik_proposed
    if ((year < mcmc$meta$T_end - 1) && !is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      loglik_next <- get.log.lik.year(year + 1, mcmc, Dlpar, phase3par)$log.lik
      loglik_next_prop <- get.log.lik.year(year + 1, mcmc, Dlpar, phase3par, tfr_proposed, prev_two = TRUE)
      
      loglik_next[indices_outliers_p1] <- 0.0
      loglik_next_prop$log.lik[indices_outliers_p1] <- 0.0
      loglik_orig <- loglik_orig + loglik_next
      loglik_proposed <- loglik_proposed + loglik_next_prop$log.lik
    }
    if (year < mcmc$meta$T_end)
    {
      loglik_mid <- get.log.lik.year(year, mcmc, Dlpar, phase3par)$log.lik
      loglik_mid_prop <- get.log.lik.year(year, mcmc, Dlpar, phase3par, tfr_proposed)
      loglik_mid[indices_outliers_0] <- 0.0
      loglik_mid_prop$log.lik[indices_outliers_0] <- 0.0
      
      loglik_orig <- loglik_orig + loglik_mid
      
      loglik_proposed <- loglik_proposed + loglik_mid_prop$log.lik
    }
    if (year > 1)
    {
      loglik_prev_prop <- get.log.lik.year(year-1, mcmc, Dlpar, phase3par, tfr_proposed, prev = FALSE)
      loglik_prev[indices_outliers_m1] <- 0.0
      loglik_prev_prop$log.lik[indices_outliers_m1] <- 0.0
      
      loglik_orig <- loglik_orig + loglik_prev
      loglik_proposed <- loglik_proposed + loglik_prev_prop$log.lik
    }
    # Here should be fine.
    loglik_orig <- loglik_orig + get.log.lik.raw(year, mcmc, mcmc$meta$tfr_all[year, ])
    loglik_proposed <- loglik_proposed + get.log.lik.raw(year, mcmc, tfr_proposed)
    
    # Accept those with likelihood not dropping so much
    accept.prob <- exp(loglik_proposed - loglik_orig)
    rv_unif <- runif(nr_countries)
    idx_accept <- which(rv_unif < accept.prob)
    mcmc$meta$tfr_all[year, idx_accept] <- tfr_proposed[idx_accept]
    if (year < mcmc$meta$T_end)
    {
      idx_update <- intersect(idx_accept, id_phase2_r)
      mcmc$eps_Tc[year, idx_update] <- loglik_mid_prop$eps[idx_update]
      mcmc$sd_Tc[year, idx_update] <- loglik_mid_prop$std[idx_update]
      tmp <- tfr_proposed[idx_update] - mcmc$S_sd
      idx <- which(tmp < 0)
      tmp[tmp > 0] <- tmp[tmp > 0] * (-mcmc$a_sd)
      tmp[idx] <- tmp[idx] * mcmc$b_sd
      mcmc$add_to_sd_Tc[year, idx_update] <- tmp
    }
    if (year > 1)
    {
      idx_update <- intersect(idx_accept, id_phase2_l)
      mcmc$eps_Tc[year - 1, idx_update] <- loglik_prev_prop$eps[idx_update]
    }
    if ((year < mcmc$meta$T_end - 1) && !is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      idx_update <- intersect(idx_accept, mcmc$meta$id_phase2_by_year[[year + 1]])
      mcmc$eps_Tc[year + 1, idx_update] <- loglik_next_prop$eps[idx_update]
    }
    
    loglik_prev <- loglik_mid
    loglik_prev[idx_accept] <- loglik_mid_prop$log.lik[idx_accept]
    
    mu0 <- mcmc$meta$tfr_mu_all[year, ]
    sigma0 <- mcmc$meta$tfr_sd_all[year, ]
    mcmc$meta$tfr_mu_all[year, ] <- (mcmc$finished.iter * mu0 + mcmc$meta$tfr_all[year, ]) / (mcmc$finished.iter + 1)
    mu <- mcmc$meta$tfr_mu_all[year, ]
    mcmc$meta$tfr_sd_all[year, ] <- sqrt(((mcmc$finished.iter-1) * sigma0 ** 2 + 
                                            mcmc$finished.iter * (mu0 - mu) ** 2 + 
                                            (mcmc$meta$tfr_all[year, ] - mu) ** 2)/mcmc$finished.iter)
    
    indices_outliers_m1 <- indices_outliers_0
    indices_outliers_0 <- indices_outliers_p1
    if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2)
    {
      indices_outliers_p1 <- sort(unique(c(mcmc$meta$yearly.outliers[[as.character(year + 1)]], 
                                           mcmc$meta$yearly.outliers[[as.character(year + 2)]])))
    }
    else
    {
      indices_outliers_p1 <- mcmc$meta$yearly.outliers[[as.character(year + 2)]]
    }
  }
  for (country in 1:nr_countries)
  {
    if (!is.null(countries) && !(country %in% countries)) next
    if (country %in% mcmc$meta3$id_phase3) 
    {
      id3 <- which(mcmc$meta$id_phase3 == country)  
      mcmc$observations[[id3]] <- mcmc$meta$tfr_all[mcmc$meta$lambda_c[country]:mcmc$meta$T_end_c[country], country]
    }
    if (country %in% mcmc$meta$id_DL){
        idx2 <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country] - 1)
        mcmc$data.list[[country]][idx2] <- mcmc$meta$tfr_all[idx2, country]
    }
  }
}


one.step.mcmc3.sampling <- function(mcmc)
{
  meta <- mcmc$meta
  niter <- mcmc$iter
  nr.countries <- meta$nr.countries
  countries.index <- meta$id_phase3
  ardata <- mcmc$observations
  Ts <- mcmc$length_obs
  gamma.mu.low <- 1/(meta$sigma.mu.prior.range[2])^2
  gamma.mu.up <- if (meta$sigma.mu.prior.range[1] == 0) NA else 1/(meta$sigma.mu.prior.range[1])^2
  gamma.rho.low <- 1/(meta$sigma.rho.prior.range[2])^2
  gamma.rho.up <- if (meta$sigma.rho.prior.range[1] == 0) NA else 1/(meta$sigma.rho.prior.range[1])^2
  gamma.eps.low <- 1/(meta$sigma.eps.prior.range[2])^2
  gamma.eps.up <- if (meta$sigma.eps.prior.range[1] == 0) NA else 1/(meta$sigma.eps.prior.range[1])^2

  #### Start MCMC

  mu.integral.to.mC <- (mu.rho.integral(mcmc[['mu']], mcmc[['sigma.mu']], low=0))^(-nr.countries)
  prop.mu <- proposal.mu.rho(mcmc[['mu.c']], mcmc[['sigma.mu']], nr.countries,
                             meta[['mu.prior.range']][1], meta[['mu.prior.range']][2])
  accept.prob <- min(((mu.rho.integral(prop.mu, mcmc[['sigma.mu']], low=0))^(-nr.countries))/mu.integral.to.mC, 1)
  if (runif(1) < accept.prob) {
    mcmc[['mu']] <- prop.mu
    mcmc[['recompute.mu.integral']] <- TRUE
  } else mcmc[['recompute.mu.integral']] <- FALSE

  # Metropolis-Hastings for sigma_mu=1/sqrt(lambda_mu)
  S <- sum((mcmc[['mu.c']]-mcmc[['mu']])^2)
  if(mcmc[['recompute.mu.integral']]) mu.integral.to.mC <- mu.rho.integral(mcmc[['mu']], mcmc[['sigma.mu']], low=0)^(-nr.countries)
  prop.lambda.mu <- rgamma.trunc((nr.countries-1)/2, S/2, low=gamma.mu.low, high=gamma.mu.up)
  accept.prob <- min(((mu.rho.integral(mcmc[['mu']], 1/prop.lambda.mu, low=0))^(-nr.countries))/mu.integral.to.mC, 1)
  if (runif(1) < accept.prob) {
    mcmc[['sigma.mu']] <- 1/sqrt(prop.lambda.mu)
    recompute.mu.integral <- TRUE
  } else recompute.mu.integral <- FALSE

  # Slice sampling for rho
  mcmc[['rho']] <- slice.sampling(mcmc[['rho']], logdensity.mu.rho, 1,
                                  low=meta[['rho.prior.range']][1], up=meta[['rho.prior.range']][2],
                                  par.c=mcmc[['rho.c']], sd=mcmc[['sigma.rho']],
                                  c.low=0, c.up=1)
  # Metropolis-Hastings for rho
  # rho.integral.to.mC <- (mu.rho.integral(mcmc[['rho']], mcmc[['sigma.rho']], low=0, up=1))^(-nr.countries)
  # prop.rho <- proposal.mu.rho(mcmc[['rho.c']], mcmc[['sigma.rho']], nr.countries,
  # meta[['rho.prior.range']][1], meta[['rho.prior.range']][2])
  # accept.prob <- min(((mu.rho.integral(prop.rho, mcmc[['sigma.rho']], low=0, up=1))^(-nr.countries))/rho.integral.to.mC, 1)
  # if (runif(1) < accept.prob) {
  # mcmc[['rho']] <- prop.rho
  # recompute.rho.integral <- TRUE
  # } else recompute.rho.integral <- FALSE

  mcmc[['sigma.rho']] <- slice.sampling(mcmc[['sigma.rho']], logdensity.sigma.mu.rho, 1,
                                        low=meta$sigma.rho.prior.range[1], up=meta$sigma.rho.prior.range[2],
                                        par.c=mcmc[['rho.c']], mean=mcmc[['rho']],
                                        c.low=0, c.up=1)

  # Metropolis-Hastings for sigma_rho=1/sqrt(lambda_rho)
  # S <- sum((mcmc[['rho.c']]-mcmc[['rho']])^2)
  # if(recompute.rho.integral)
  # rho.integral.to.mC <- (mu.rho.integral(mcmc[['rho']], mcmc[['sigma.rho']], low=0, up=1))^(-nr.countries)
  # prop.lambda.rho <- rgamma.trunc((nr.countries-1)/2, S/2, low=gamma.rho.low, high=gamma.rho.up)
  # accept.prob <- min(((mu.rho.integral(mcmc[['rho']], 1/prop.lambda.rho, low=0, up=1))^(-nr.countries))/rho.integral.to.mC, 1)
  # if (runif(1) < accept.prob) {
  # mcmc[['sigma.rho']] <- 1/sqrt(prop.lambda.rho)
  # recompute.rho.integral <- TRUE
  # } else recompute.rho.integral <- FALSE

  sigma.eps.sq <- mcmc$sigma.eps^2
  sigma.mu.sq <- mcmc$sigma.mu^2
  sigma.mu.sq.inv <- 1/sigma.mu.sq
  mu.over.sigma.sq <- mcmc$mu/sigma.mu.sq
  sigma.rho.sq <- mcmc$sigma.rho^2
  sigma.rho.sq.inv <- 1/sigma.rho.sq
  rho.over.sigma.sq <- mcmc$rho/sigma.rho.sq
  S.eps <- STn <- 0
  one.minus.rho <- 1-mcmc$rho.c
  one.minus.rho.sq <- one.minus.rho^2
  W <- one.minus.rho.sq/sigma.eps.sq

  # country-specific parameters - Gibbs sampler
  for(country in 1:nr.countries) {
    f.ct <- ardata[[country]][2:Ts[country]]
    f.ctm1 <- ardata[[country]][1:(Ts[country]-1)]
    # mu.c
    s <- sum((f.ct - mcmc$rho.c[country]*f.ctm1)/one.minus.rho[country])
    nomin <- W[country] * s + mu.over.sigma.sq
    denom <- (Ts[country]-1) * W[country] + sigma.mu.sq.inv
    mcmc$mu.c[country] <- rnorm.trunc(mean=nomin/denom, sd=1/sqrt(denom), low=0)

    # rho.c
    d1 <- f.ctm1 - mcmc$mu.c[country]
    a <- sum(d1^2)/sigma.eps.sq
    b <- sum(d1*(f.ct - mcmc$mu.c[country]))/sigma.eps.sq
    nomin <- b + rho.over.sigma.sq
    denom <- a + sigma.rho.sq.inv
    mcmc$rho.c[country] <- rnorm.trunc(mean=nomin/denom, sd=1/sqrt(denom), low=0, #high=1-10*.Machine$double.xmin
                                       high=0.999999)
    S.eps <- S.eps + sum((f.ct - mcmc$mu.c[country] - mcmc$rho.c[country]*d1)^2)
    STn <- STn + Ts[country]-1
  }
  # Gibbs for sigma.eps
  mcmc$sigma.eps <- 1/sqrt(rgamma.trunc((STn-1)/2, S.eps/2, low=gamma.eps.low, high=gamma.eps.up))

}