R/functions.R

In ABRSQOL: Quality-of-Life Solver for "Measuring Quality of Life under Spatial Frictions"

#' ABRSQOL numerical solution algorithm to invert a quality of life measure
#'
#' @description
#' This toolkit implements a numerical solution algorithm
#' to invert a quality of life (QoL) from observed data
#' in various programming languages. The QoL measure is
#' based on Ahlfeldt, Bald, Roth, Seidel (2024):
#' Measuring quality of life under spatial frictions.
#' Unlike the traditional Rosen-Roback measure, this measure
#' accounts for mobility frictions—generated by idiosyncratic
#' tastes and local ties—and trade frictions—generated by
#' trade costs and non-tradable services, thereby reducing
#' non-classical measurement error.
#' When using this programme or the toolkit in your work, please cite the paper.
#'
#' @details
#' Notice that quality of life is identified up to a constant.
#' Therefore, the inverted QoL measures measure has a relative
#' interpretation only. We normalize the QoL relative to the first
#' observation in the data set. It is straightforward to rescale
#' the QoL measure to any other location or any other value (such
#' as the mean or median in the distribution of QoL across locations).
#'
#' @param df input data containing variables refenced by
#' following arguments: data.frame or matrix
#' @param w wage index variable name or column index:
#' character or integer (or list), default='w'
#' @param p_H floor_space_price variable name or column index:
#' character or integer, default='p_H'
#' @param P_t tradable_goods_price variable name or column index:
#' character or integer, default='P_t'
#' @param p_n local_services_price variable name or column index:
#' character or integer, default='p_n'
#' @param L residence_population variable name or column index:
#' character or integer, default='L'
#' @param L_b hometown_population variable name or column index:
#' character or integer, default='L_b'
#' @param alpha Income share on non-housing consumtpion:
#' double, default=0.7
#' @param beta Share of tradable goods in non-housing consumption:
#' double, default=0.5
#' @param gamma Idiosyncratic taste dispersion (inverse labour
#' supply elasticity):
#' double, default=3
#' @param xi Valuation of local ties: double, default=5
#' @param conv Convergence parameter (Hgher value increases spead of,
#' convergence and risk of bouncing): double, default=0.5
#' @param tolerance Value used in stopping rule (The mean absolute error (MAE).
#' Smaller values imply greater precision and longer convergence):
#' double, default=1e-10
#' @param maxiter Maximum number of iterations after which the algorithm
#' is forced to stop: integer, default=10000
#'
#' @return inverted quality of life measure as Numeric vector
#' (identified up to a constant)
#'
#' @examples
#' # Example 1:
#' # load testdata,
#' # run QoL inversion with default parameters and append result
#' data(ABRSQOL_testdata)
#' my_dataframe <- ABRSQOL_testdata
#' my_dataframe$qol1 <- ABRSQOL(df=ABRSQOL_testdata)
#' my_dataframe
#'
#' # Example 2: load your data from csv, run inversion, save result as csv
#' # my_dataframe <- read.csv("path/to/your/csv_filename.csv")
#' # my_dataframe$qol2 <- ABRSQOL(
#' #  # supply your dataset as a dataframe
#' #  df=my_dataframe,
#' #  # and specify the corresponding variable name for your dataset
#' #  w = 'wage',
#' #  p_H = 'floor_space_price',
#' #  P_t = 'tradable_goods_price',
#' #  p_n = 'local_services_price',
#' #  L = 'residence_pop',
#' #  L_b = 'L_b',
#' #  # freely adjust remaining parameters
#' #  alpha = 0.7,
#' #  beta = 0.5,
#' #  gamma = 3,
#' #  xi = 5.5,
#' #  conv = 0.3,
#' #  tolerance = 1e-11,
#' #  maxiter = 50000
#' #)
#' #write.csv(my_dataframe, 'qol_of_my_data.csv'
#'
#' # Example 3: Reference variables in your dataset by using the column index
#' my_dataframe$qol3 <- ABRSQOL(
#'   df=my_dataframe,
#'   w = 1,
#'   p_H = 3,
#'   P_t = 4,
#'   p_n = 2,
#'   L = 6,
#'   L_b = 5
#' )
#'
#' # Example 4: Having named the variables in your data according to the
#' # default parameters, you can ommit specifying variable names
#' my_dataframe$qol4 <- ABRSQOL(
#'   df=my_dataframe,
#'   alpha = 0.7,
#'   beta = 0.5,
#'   gamma = 3,
#'   xi = 5.5,
#'   conv = 0.5
#' )
#'
#' @export
ABRSQOL <- function(
  df, # data.frame or matrix containing dataset

  # SPECIFY VARIABLES NAMES OR COLUMN INDEX
  w = "w", # 2: Wage index
  p_H = "p_H", # 3: Floor space price index
  P_t = "P_t", # 4: Tradable goods price index
  p_n = "p_n", # 5: Local services price index
  L = "L", # 6: Residence population
  L_b = "L_b", # 7: Hometown population

  # DEFINE PARAMETER VALUES
  alpha = 0.7, #income share on non-housing; 1-alpha expenditure on housing
  beta = 0.5, # share of alpha that is spent on tradable good
  gamma = 3, # Own calculations
  xi = 5.5, # Own calculations
  # CONVERGENCE AND STOPPING PARAMETERS
  conv = 0.5, # convergence parameter
  tolerance = 1e-10, # Tolerance level for loop
  maxiter = 10000
) {

  if (maxiter < 1) {
    stop("Maximum iterations must at least 1.
    Larger values like default of 10000 are recommended.")
  }

  # Extract key variables from input dataframe/matrix
  L_b <- matrix(data=as.vector(unlist(df[L_b])), ncol=length(L_b), byrow=FALSE)
  L <- matrix(data = as.vector(unlist(df[L])), ncol = length(L), byrow = FALSE)
  w <- matrix(data = as.vector(unlist(df[w])), ncol = length(w), byrow = FALSE)
  P_t <- df[[P_t]]
  p_H <- df[[p_H]]
  p_n <- df[[p_n]]


  # if there are unequal number of rows (n_obs) among variables throw error
  if(length(unique(
    c(length(L_b), length(L), length(w), length(P_t) ,length(p_H), length(p_n))
  )) > 1){
    error_msg <- paste(
      "\nVariables do not have the same length:",
      "\nL_b: ", length(L_b), "\nL: ", length(L), "\nw: ", length(w),
      "\nP_t: ", length(P_t), "\np_H: ", length(p_H), "\np_n: ", length(p_n)
    )
    stop(error_msg)
  }
  # else save units of observation as J
  J <- length(L_b)

  # if there are unequal number of dimensions throw an error
  if(length(unique(c(dim(L_b)[2],dim(L)[2],dim(w)[2]))) > 1){
    error_msg <- paste(
      "\nVariables do not have the same dimension:",
      "\nL_b: ", dim(L_b)[2], "\nL: ", dim(L)[2], "\nw:", dim(w)[2]
    )
    stop(error_msg)
  }
  # else assign theta as the number of dimensions (mostly will be 1)
  Theta <- dim(L_b)[2]


  ## Inversion

  # Adjust L_b to have same sum as L
  L_bar <- sum(L) # total number of workers in dataset
  L_b_adjust <- L_bar / apply(L_b, 2, function(x) sum(x))
  L_b <- t(t(L_b) * L_b_adjust)

  # Express all variables in relative differences
  # Calculate relative employment, L_hat
  L_hat <- sweep(L, 2, L[1,], `/`)
  # Calculate relative wages, w_hat
  w_hat <- sweep(w, 2, w[1,], `/`)
  # Calculate relative price levels
  P_t_hat <- P_t / P_t[1]
  p_H_hat <- p_H / p_H[1]
  p_n_hat <- p_n / p_n[1]

  # Calculate aggregate price level
  P_hat<-(P_t_hat^(alpha * beta))*(p_n_hat^(alpha*(1-beta)))*(p_H_hat^(1-alpha))
  P    <-(P_t    ^(alpha * beta))*(p_n    ^(alpha*(1-beta)))*(p_H   ^(1-alpha))

  # Relative Quality of life (A_hat)
  # Guess values relative QoL
  A_hat <- matrix(1, J, Theta) # First guess: all locations have the same QoL
  A <- A_hat

  O_total <- Inf # Starting value for loop
  iter_count <- 1 # Counts the number of iterations

  while (O_total > tolerance && iter_count <= maxiter){

    # (1) Calculate model-consistent aggregation shares, Psi_b
    nom <- (as.vector(A) * w  * as.vector(1/P)) ^(gamma)
    Psi_b<-(sweep((exp(xi)-1)*nom,2,apply(nom,2,function(x) sum(x)),`/`)+1)^(-1)

    # (2) Calculate mathcal_L
    mathcal_L <- apply(L_b * Psi_b, 2, function(x) sum(x))+L_b*Psi_b*(exp(xi)-1)

    # (3) Calculate relative mathcal_L
    mathcal_L_hat <- sweep(mathcal_L, 2, mathcal_L[1, ], `/`)

    # (4) Calculate relative QoL, A_hat, according to equation (17)
    A_hat_up <- as.vector(P_hat) * (1/w_hat) * (L_hat/mathcal_L_hat) ^ (1/gamma)

    # (5) Calculate deviations from inital guesses for QoL levels
    O_total <- sum(abs(A_hat_up - A_hat)) / J

    # Update QoL levels for next iteration of loop
    A_hat <- conv * A_hat_up + (1 - conv) * A_hat
    A <- A_hat

    message(
      "\rSolving for QoL: ",
      "Iteration: ", iter_count, "/",maxiter,
      ". Objective function: ", O_total, " > ", tolerance,
      appendLF = FALSE
    )

    # Next iteration
    iter_count <- iter_count + 1
  }

  if (O_total > tolerance) {
    warning(
      "WARNING: Exceeded maximum number of iterations: ", iter_count, "/",
      maxiter, ". Objective function: ", O_total, " >= ", tolerance
    )
  } else {
    message(
      "\rSUCCESS. Solved for QoL at iteration: ", iter_count, "/", maxiter,
      ". Objective function: ", O_total, " < ", tolerance
    )
  }

  return(as.vector(A))

}