ASFR solver: Given target fertility parameters determine suitable ASFRs

asfr_solver.R

# ASFR Solver
# Given target fertility parameters determine suitable ASFRs
#
# Jonas Schöley <[email protected]>
#
# This script determines a set of age specific fertility rates under
# the following constraints:
# a) a total fertility rate
# b) a mean age at first birth
# c) a set of proportional deviations from initial (reference) rates
#
# Weights can be given to the relative importance of each of the
# constraints.

# input parameters
input <- list(
  # midpoints of age groups
  ages = c(15, 20, 25, 30, 35, 40, 45)+2.5,
  # width of last age group
  wlast = 5,
  # initial ASFRs
  initial_asfrs = c(0.0031, 0.0269, 0.0671, 0.0888, 0.0512, 0.0134, 0.0012),
  # target TFR
  target_tfr = 1.51,
  # target mean age at first child
  target_mafb = 32.9,
  # target ratio to initial ASFRs
  target_ratio = c(0.697, 0.842, 0.980, 1.499, 1.888, 1.489, 1.917),
  # weights for the 3 constraints 
  weights = c(tfr=10,mafb=1,ratios=1)
)

FindMatchingASFRs <- function (input) {
  
  # constants
  n = length(input$ages)
  w = c(diff(input$ages), input$wlast)
  
  # the loss function of the optimization problem
  Loss <- function (delta, input, n, w) {
    prop_delta <- exp(delta)
    optim_asfrs <- input$initial_asfrs*prop_delta
    mafb <- t(input$ages)%*%prop.table(optim_asfrs)
    tfr <- t(w)%*%c(optim_asfrs)
    loss <-
      input$weights[1]*(tfr-input$target_tfr)^2 +
      input$weights[2]*(mafb-input$target_mafb)^2 +
      sum(delta^2) +
      input$weights[3]*sum((input$target_ratio-prop_delta)^2)
    return(c(loss))
  }
  
  # optimize
  fit <- optim(rep(1, n), Loss, gr = NULL, input, n, w,
               control = list(maxit = 1e4))
  
  # results
  ratio_optim <- exp(fit$par)
  asfr_optim <- input$initial_asfrs*ratio_optim
  mafb_optim <- c(t(input$ages)%*%prop.table(asfr_optim))
  tfr_optim <- c(t(w)%*%c(asfr_optim))
  result <- list(
    asfr = asfr_optim, tfr = tfr_optim, mafb = mafb_optim,
    ratio = ratio_optim, fit = fit
  )
  
  return(result)
  
}

FindMatchingASFRs(input)