12  Merging: annual data

Objective

As in Chapters 5, we merge the different datasets to produce the data table used in the local projections that uses annual data instead of monthl or quarterly data.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(cli)

12.1 Load Intermediate Files

The (annual) weather data (Chapter 1) can be loaded:

load("../data/output/weather/weather_annual_regions_df.rda")

The agricultural data (Chapter 2):

load("../data/output/minagri/dataset_agri_2001_2015.rda")

The macroeconomic data (Chapter 3):

load("../data/output/macro/df_macro.rda")

The share of natural regions and the El Niño–Southern Oscillation (Chapter 4):

load("../data/output/natural_region_dep.rda")
load("../data/output/weather/ONI_temp.rda")

12.2 Merge the Files

We add ENSO data to the weather dataset:

Weather <- weather_annual_regions_df |> 
  # Add ENSO data
  left_join(
    ONI_temp |> 
      mutate(year = as.numeric(Year)) |> 
      group_by(year) |> 
      summarise(ONI = mean(ONI)), 
    by = c(
      "year" = "year"
    )
  ) |>
  group_by(IDDPTO) |> 
  mutate( 
    temp_min_dev_ENSO   = temp_min - mean(temp_min),
    temp_max_dev_ENSO   = temp_max - mean(temp_max),
    temp_mean_dev_ENSO  = temp_mean - mean(temp_mean),
    precip_sum_dev_ENSO = precip_sum - mean(precip_sum))|> 
  ungroup() |> 
  labelled::set_variable_labels(
    temp_min_dev_ENSO   = "Deviation of Min. Temperature from ENSO Normals",
    temp_max_dev_ENSO   = "Deviation of Max. Temperature from ENSO Normals",
    temp_mean_dev_ENSO  = "Deviation of Mean Temperature from ENSO Normals",
    precip_sum_dev_ENSO = "Deviation of Total Rainfall from ENSO Normals",
  )

For international prices:

library(readxl)
int_prices <- read_excel(
  path = "../data/raw/Macro/IMF_DATA.xls",
  sheet = "SELECTED",
  col_types = "text") |> 
  mutate(date = lubridate::ymd(str_c(YEAR,  "-", MONTH, "-01"))) |> 
  rename(
    "FPI"           = "FPI_PFOOD", 
    "FERTILIZER"    = "FERTILIZER_PFERT",   
    "IndexOIL"      = "IndexOIL_POILAPSP",
    "PriceOIL"      = "PriceOIL_POILAPSP", 
    "CORN"          = "CORN_PMAIZMT", 
    "ARROZ CÁSCARA" = "RICE_PRICENPQ", 
    "TRIGO"         = "WHEAT_PWHEAMT"
  ) |> 
  mutate(CORN2 = CORN, FPI2 = FPI) |> 
  select(-c(MONTH, YEAR, CPI_Peru_IMF, CPI_US_IMF, IndexOIL, PriceOIL)) |> 
  pivot_longer(cols = !date, names_to = "product", values_to = "price_int") |> 
  mutate(
    product = case_when(
      product == "FPI"  ~ "PAPA", 
      product == "FPI2" ~ "YUCA", 
      product == "CORN" ~ "MAÍZ AMILÁCEO", 
      product == "CORN2"~ "MAÍZ AMARILLO DURO", 
      TRUE ~ product
    ), 
    product_eng = case_when(
      product == "PAPA" ~ "Potato",
      product == "YUCA" ~ "Cassava",
      product == "ARROZ CÁSCARA" ~ "Rice",
      product == "MAÍZ AMARILLO DURO" ~ "Dent corn"
    ),
    price_int = as.numeric(price_int)
  ) |> 
  arrange(product, date) |> 
  mutate(year = year(date)) |>  
  group_by(product, year) |> 
  mutate(price_int = mean(price_int)) |>  
  ungroup() |> 
  select(-date) |> 
  unique() |> 
  group_by(product) |> 
  arrange(year) |> 
  mutate(int_price_inf = (price_int/lag(price_int) - 1) * 100) |> 
  ungroup() |> 
  arrange(product, year) |> 
  filter(! year == 2000)

Let us merge all these datasets in a single one:

data_total <- 
  data_total |> 
  # Add macroeconomic data
  left_join(
    df_macro |> rename(gdp = y),
    by = "date"
  ) |> 
  dplyr::select(
    product_eng, region,region_id, product, year, Value_prod,
    rer_hp, r_hp, pi, ind_prod) |> 
  group_by(region, product, year) |> 
  mutate(
    Value_prod     = sum(Value_prod),   
    rer_hp         = mean(rer_hp),
    r_hp           = mean(r_hp),
    pi             = mean(pi),
    ind_prod       = mean(ind_prod)
  ) |>  
  ungroup() |> 
  unique() |> 
  # Add commodity prices data
  left_join(
    int_prices,
    by =  c(
      "product", "product_eng", "year")
  ) |> 
  group_by(product, region) |> 
  mutate(
    int_price_inf = mFilter::hpfilter(
      int_price_inf, freq = 1600, type = "lambda")[["cycle"]]
  ) |> 
  # Add weather data and ENSO 
  left_join(
    Weather |> 
      dplyr::select(-IDDPTO),
    by = c(
      "year" = "year",
      "region" = "DEPARTAMEN"
    )
  )

save(data_total, file = "../data/output/dataset_2001_2015-annual.rda")

Here are the first rows of that tibble:

data_total

# A tibble: 2,160 × 45
# Groups:   product, region [144]
   product_eng     region   region_id product   year Value_prod   rer_hp    r_hp
   <chr>           <chr>        <int> <chr>    <dbl>      <dbl>    <dbl>   <dbl>
 1 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2001     10110  -0.0828   0.986 
 2 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2002     11493  -0.218   -1.41  
 3 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2003      5850   0.444    0.220 
 4 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2004      7315.  0.00578 -0.109 
 5 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2005      6126. -0.818   -0.207 
 6 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2006      4471.  0.00305  0.269 
 7 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2007      5240.  1.57     0.0279
 8 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2008      4909. -0.427    1.17  
 9 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2009      5660. -0.344   -0.456 
10 Amylaceous corn AMAZONAS         1 MAÍZ AM…  2010      4946. -1.04    -1.26  
# ℹ 2,150 more rows
# ℹ 37 more variables: pi <dbl>, ind_prod <dbl>, price_int <dbl>,
#   int_price_inf <dbl>, temp_min <dbl>, temp_max <dbl>, temp_mean <dbl>,
#   precip_sum <dbl>, precip_piscop_sum <dbl>, perc_gamma_precip <dbl>,
#   perc_gamma_precip_piscop <dbl>, gdd_rice <dbl>, gdd_maize <dbl>,
#   gdd_potato <dbl>, gdd_cassava <dbl>, hdd_rice <dbl>, hdd_maize <dbl>,
#   hdd_potato <dbl>, hdd_cassava <dbl>, temp_min_dev <dbl>, …

12.3 Dataset for the Local Projections

Now, let us create the dataset specifically used to estimate the models.

Let us make sure that the region data are encoded as a factor.

data_total <- 
  data_total |> 
  mutate(region_id = factor(region_id))

The crops we focus on:

crops <- c("Rice", "Dent corn", "Potato", "Cassava")

The number of observation in each region, for each crop:

data_total |> 
  group_by(product_eng, region_id) |> 
  summarise(n = sum(Value_prod <= 0)) |> 
  arrange(desc(n))

`summarise()` has grouped output by 'product_eng'. You can override using the
`.groups` argument.

# A tibble: 144 × 3
# Groups:   product_eng [6]
   product_eng     region_id     n
   <chr>           <fct>     <int>
 1 Amylaceous corn 15           15
 2 Amylaceous corn 21           15
 3 Amylaceous corn 23           15
 4 Amylaceous corn 24           15
 5 Potato          15           15
 6 Potato          21           15
 7 Potato          23           15
 8 Potato          24           15
 9 Rice            3            15
10 Rice            8            15
# ℹ 134 more rows

12.3.1 Definition of the Variable of Interest

Warning

We compute percentage deviation of production from the regional average over the period of interest, but we will actually not use those values in the subsequent estimations. We will, however, use the demeaned values computed in this section.

This section outlines a two-step procedure for expressing agricultural production data at the annual regional level for a specific crop as a percentage deviation from the regional crop-specific average over the period of interest. The procedure involves handling missing values.

• Step 1: Handling Missing Values

  In the first step, we address missing values by linear interpolation. This approach helps us estimate the missing values by considering the neighboring data points.

  • Step 1.1: Imputing missing values with linear interpolation.

    The missing values get replaced by linear interpolation. However, if there are more than two consecutive missing values, they are not replaced with interpolated values. Instead, the series for the specific crop in the given region is split based on the locations of the missing values. The split with the highest number of consecutive non-missing values is retained, while the other splits are discarded.

  • Step 1.2: Dropping Series with Remaining Missing Values

    After imputing missing values using the moving median, we check if any missing values still remain in the dataset. If there are any remaining missing values for a particular series, we choose to exclude that series from further analysis. By doing so, we ensure that the subsequent detrending process is performed only on reliable and complete data.

• Step 2: Normalized Agricultural Production

  For each month ( m ), region ( i ), and crop ( c ), we calculate the average production over the entire period (January 2001 to December 2015): \[\overline{y}_{c,i,m} = \frac{1}{n_{T_c}} \sum_{t=1}^{T_c} y_{c,i,m,t}^{\text{raw}} \] Then, we express agricultural production relative to the average: \[y_{c,i,m,t} = \begin{cases} \frac{y_{c,i,m,t}^{\text{raw}}}{\overline{y}_{c,i,m}}, & \overline{y}_{c,i,m} > 0\\ 0, & \overline{y}_{c,i,m} = 0 \end{cases}\] Values of \(y_{c,i,m,t}>1\) means that the production for crop \(c\) in region \(i\) during month \(m\) of year \(t\) is higher than the average monthly production for that crop and region over the period 2001 to 2015. For example, a value of 1.5 means that the production is 50% higher than average.

• Step 2 (alternative version): Deviation from regional monthly average, in percent (this step is useless in the new version of the analysis: it lead to discard too many observations)

  Once we have addressed the missing values, we proceed to the second step, which consists in computing the deviation of production from the regional average. First, we compute the average production of each crop \(c\) in each region \(i\): \[\overline{y}_{c,i} = \frac{1}{n_{T_c}} \sum_{t=1}^{T_c} y_{c,i,t}^{raw}\] Then, we compute the percentage deviation from this average at each date (i.e., year) \(t\): \[y_{c,i,t} = \frac{y_{c,i,t}^{raw} - \overline{y}_{c,i}}{\overline{y}_{c,i}}\]

Let us implement this process in R. First, we need to define two functions to handle the missing values:

• The get_index_longest_non_na() function retrieves the indices of the longest consecutive sequence without missing values from a given input vector. It helps us identify the positions of elements in that sequence.

• The keep_values_longest_non_na() function uses the obtained indices to create a logical vector. Each element of this vector indicates whether the corresponding element in the input vector belongs to the longest consecutive sequence of non-missing values. This allows us to filter the data and retain only the values from the longest consecutive sequence without missing values.

These two functions combined help us handle missing data in the weather series and ensure that we work with the most complete sequences for each region and crop.

The first function:

#' Returns the index of the longest sequence of non NA values in a vector
#'
#' @param y vector of numerical values
#' @export
get_index_longest_non_na <- function(y) {
  split_indices <- which(is.na(y))
  nb_obs <- length(y)

  if (length(split_indices) == 0) {
    res <- seq_len(nb_obs)
  } else {
    idx_beg <- c(1, split_indices)
    if (idx_beg[length(idx_beg)] != nb_obs) {
      idx_beg <- c(idx_beg, nb_obs)
    }
    lengths <- diff(idx_beg)
    ind_max <- which.max(lengths)
    index_beginning <- idx_beg[ind_max]
    if(!index_beginning == 1 | is.na(y[index_beginning])) {
      index_beginning <- index_beginning + 1
    }
    index_end <- idx_beg[ind_max] + lengths[ind_max]
    if(is.na(y[index_end])) {
      index_end <- index_end - 1
    }
    res <- seq(index_beginning, index_end)
  }
  res
}

The second one:

#' Returns a logical vector that identifies the longest sequence of non NA
#' values within the input vector
#' 
#' @param y numeric vector
keep_values_longest_non_na <- function(y) {
  ids_to_keep <- get_index_longest_non_na(y)
  ids <- seq(1, length(y))
  ids %in% ids_to_keep
}

Note

Those two functions are defined in weatherperu/R/utils.R.

We define a function, pct_prod_production(), that takes the data frame of observations as input, as well as a crop name and a region ID. It returns a tibble with the following variables:

• product_eng: the English name of the crop
• region_id: the ID of the region
• year: year
• y_new_normalized (our variable of interest in Chapter 14): the production demeaned by the average for the crop of interest in the region of interest
• y_new: the production (in tons) where missing values were imputed, if possible
• y_dev_pct: the production expressed as the percentage deviation from the average (for the crop of interest, in the region of interest)
• y: same as y_dev_pct but without an estimated quadratic trend estimated by OLS
• t: trend.

#' Computes the percentage deviation of production from annual regional average
#'
#' @param df data
#' @param crop_name name of the crop
#' @param region_id id of the region
#'
#' @returns data frame with the product, the region id, the date, the production
#' with imputed missing values (`y_new`), the production demeaned (`y_new_normalized`), 
#' the percentage deviation from mean production (`y_dev_pct`), the percentage 
#' deviation from mean production minus an estimated quadratic trend (estimated
#'  by OLS) (`y`), and, a trend (`t`)
#' @export
#' @importFrom dplyr filter arrange mutate select row_number group_by
#' @importFrom tidyr nest unnest
#' @importFrom purrr map
#' @importFrom imputeTS na_interpolation
#' @importFrom stats lm predict residuals
pct_prod_production <- function(df,
                                crop_name,
                                region_id) {
  # The current data
  df_current <-
    df |>
    filter(
      product_eng == !!crop_name,
      region_id == !!region_id
    ) |>
    arrange(year)

  ## Dealing with missing values ----
  # Look for negative production values
  df_current <-
    df_current |>
    mutate(
      y_new = ifelse(Value_prod < 0, NA, Value_prod)
    )

  if (any(is.na(df_current$y_new))) {

    # Replacing NAs by interpolation
    # If there are more than two contiguous NAs, they are not replaced
    df_current <-
      df_current |>
      mutate(
        y_new = imputeTS::na_interpolation(y_new, maxgap = 3)
      )

    # Removing obs at the beginning/end if they are still missing
    df_current <-
      df_current |>
      mutate(
        row_to_keep = !(is.na(y_new) & row_number() %in% c(1:2, (n()-1):(n())))
      ) |>
      filter(row_to_keep) |>
      select(-row_to_keep)

    # Keeping the longest series of continuous non-NA values
    df_current <-
      df_current |>
      mutate(
        row_to_keep = keep_values_longest_non_na(y_new)
      ) |>
      filter(row_to_keep) |>
      select(-row_to_keep)
  }

  rle_y_new <- rle(df_current$y_new)
  check_contiguous_zeros <- rle_y_new$lengths[rle_y_new$values==0]


  ## Percent deviation from regional average over the period
  resul <-
    df_current |>
    ungroup() |>
    mutate(
      y_new_normalized = case_when(
        mean(y_new) == 0~ 0,
        TRUE ~ y_new / mean(y_new)
      ),
      y_dev_pct = case_when(
        mean(y_new) == 0 ~ 0,
        TRUE ~ (y_new - mean(y_new)) / mean(y_new)
      )
    ) |>
    arrange(year) |>
    mutate(t = row_number()) |>
    nest(.by = c(product_eng, region_id)) |>
    mutate(
      ols_fit   = map(data, ~ lm(y_new_normalized ~ -1 + t + I(t^2), data = .x)),
      resid     = map(ols_fit, residuals),
      fitted    = map(ols_fit, predict)
    ) |>
    unnest(cols = c(data, resid, fitted)) |>
    mutate(
      y = resid
    ) |>
    select(
      product_eng, region_id, year,
      y_new, y_dev_pct, y_new_normalized, y, t
    ) |>
    ungroup() |>
    arrange(year)

  resul
}

For example, for potatoes in region with id 1:

pct_prod_production(
  df = data_total, 
  crop_name = "Potato", 
  region_id = 1
)

# A tibble: 15 × 8
   product_eng region_id  year  y_new y_dev_pct y_new_normalized        y     t
   <chr>       <fct>     <dbl>  <dbl>     <dbl>            <dbl>    <dbl> <int>
 1 Potato      1          2001 50077    -0.188             0.812  0.566       1
 2 Potato      1          2002 57096    -0.0747            0.925  0.462       2
 3 Potato      1          2003 57630    -0.0660            0.934  0.281       3
 4 Potato      1          2004 50169.   -0.187             0.813 -0.00355     4
 5 Potato      1          2005 47543.   -0.230             0.770 -0.182       5
 6 Potato      1          2006 53307.   -0.136             0.864 -0.197       6
 7 Potato      1          2007 74439.    0.206             1.21   0.0643      7
 8 Potato      1          2008 85019.    0.378             1.38   0.182       8
 9 Potato      1          2009 80966.    0.312             1.31   0.0899      9
10 Potato      1          2010 66088.    0.0710            1.07  -0.150      10
11 Potato      1          2011 51551.   -0.165             0.835 -0.358      11
12 Potato      1          2012 59051.   -0.0430            0.957 -0.181      12
13 Potato      1          2013 66423.    0.0765            1.08   0.0217     13
14 Potato      1          2014 59995.   -0.0277            0.972  0.0278     14
15 Potato      1          2015 66222.    0.0732            1.07   0.266      15

We can apply this function to all crops of interest, in each region. Let us define a table that contains all the possible values for the combination of crops and regions:

product_and_regions <- 
  data_total |> 
  filter(product_eng %in% crops) |> 
  dplyr::select(product_eng, region_id, region) |> 
  unique()

Adding missing grouping variables: `product`

Then we apply the pct_prod_production() function to all these different cases, and store the results in a list named df_pct_pred_production:

df_pct_pred_production <- vector(mode = "list", length = nrow(product_and_regions))
cli_progress_bar(total = nrow(product_and_regions))
for(i in 1:nrow(product_and_regions)){
  df_pct_pred_production[[i]] <- pct_prod_production(
    df = data_total, 
    crop_name = product_and_regions$product_eng[i], 
    region_id = product_and_regions$region_id[i]
  )
  cli_progress_update(set = i)
}

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

The elements of the list are all tibbles, with the same column names. We can merge them in a single tibble.

df_pct_pred_production <- bind_rows(df_pct_pred_production)

We can have a look at the number regions with 0 values for the agricultural production.

df_pct_pred_production |> 
  group_by(product_eng, region_id) |> 
  summarise(nb_0 = sum(y_new == 0)) |> 
  arrange(desc(nb_0))

`summarise()` has grouped output by 'product_eng'. You can override using the
`.groups` argument.

# A tibble: 96 × 3
# Groups:   product_eng [4]
   product_eng region_id  nb_0
   <chr>       <fct>     <int>
 1 Cassava     8            15
 2 Cassava     22           15
 3 Potato      15           15
 4 Potato      16           15
 5 Potato      21           15
 6 Potato      23           15
 7 Potato      24           15
 8 Rice        3            15
 9 Rice        8            15
10 Rice        10           15
# ℹ 86 more rows

Now, let us add the other columns to the tibble that contains the new data:

df <- df_pct_pred_production |> 
  left_join(
    data_total,
    join_by(product_eng, region_id, year)
  )

Let us also impute missing values for the weather variables.

weather_variables <- 
  weather_annual_regions_df |> 
  select(where(is.numeric)) |> 
  select(-year) |> 
  colnames()

The current number of missing values:

df |> 
  summarise(
    across(
      .cols = !!weather_variables,
      .fns = ~ sum(is.na(.x)),
      .names = "{.col}_nb_na"
    )
  ) |> 
  unlist()

                temp_min_nb_na                 temp_max_nb_na 
                             0                              0 
               temp_mean_nb_na               precip_sum_nb_na 
                             0                              0 
       precip_piscop_sum_nb_na        perc_gamma_precip_nb_na 
                             0                              0 
perc_gamma_precip_piscop_nb_na                 gdd_rice_nb_na 
                             0                              0 
               gdd_maize_nb_na               gdd_potato_nb_na 
                             0                              0 
             gdd_cassava_nb_na                 hdd_rice_nb_na 
                             0                              0 
               hdd_maize_nb_na               hdd_potato_nb_na 
                             0                              0 
             hdd_cassava_nb_na             temp_min_dev_nb_na 
                             0                              0 
            temp_max_dev_nb_na            temp_mean_dev_nb_na 
                             0                              0 
          precip_sum_dev_nb_na    precip_piscop_sum_dev_nb_na 
                             0                              0 
            gdd_rice_dev_nb_na            gdd_maize_dev_nb_na 
                             0                              0 
          gdd_potato_dev_nb_na          gdd_cassava_dev_nb_na 
                             0                              0 
            hdd_rice_dev_nb_na            hdd_maize_dev_nb_na 
                             0                              0 
          hdd_potato_dev_nb_na          hdd_cassava_dev_nb_na 
                             0                              0 

In case of missing values, we use linear interpolation to replace them:

df <- 
  df |> 
  mutate(
    across(
      .cols = !!weather_variables,
      .fns = ~ imputeTS::na_interpolation(.x, maxgap = 3)
    )
  )

The number of remaining missing values:

df |> 
  summarise(
    across(
      .cols = !!weather_variables,
      .fns = ~ sum(is.na(.x)),
      .names = "{.col}_nb_na"
    )
  ) |> 
  unlist()

                temp_min_nb_na                 temp_max_nb_na 
                             0                              0 
               temp_mean_nb_na               precip_sum_nb_na 
                             0                              0 
       precip_piscop_sum_nb_na        perc_gamma_precip_nb_na 
                             0                              0 
perc_gamma_precip_piscop_nb_na                 gdd_rice_nb_na 
                             0                              0 
               gdd_maize_nb_na               gdd_potato_nb_na 
                             0                              0 
             gdd_cassava_nb_na                 hdd_rice_nb_na 
                             0                              0 
               hdd_maize_nb_na               hdd_potato_nb_na 
                             0                              0 
             hdd_cassava_nb_na             temp_min_dev_nb_na 
                             0                              0 
            temp_max_dev_nb_na            temp_mean_dev_nb_na 
                             0                              0 
          precip_sum_dev_nb_na    precip_piscop_sum_dev_nb_na 
                             0                              0 
            gdd_rice_dev_nb_na            gdd_maize_dev_nb_na 
                             0                              0 
          gdd_potato_dev_nb_na          gdd_cassava_dev_nb_na 
                             0                              0 
            hdd_rice_dev_nb_na            hdd_maize_dev_nb_na 
                             0                              0 
          hdd_potato_dev_nb_na          hdd_cassava_dev_nb_na 
                             0                              0 

We add labels to the new columns:

df <- 
  df |> 
  labelled::set_variable_labels(
    y_new = "Quarterly Agricultural Production (tons)",
    y_dev_pct = "Quarterly Agricultural Production (pct. deviation from regional quarterly mean)"
  )

12.3.2 Saving the file

The dataset that can be used to estimate the impact of weather shocks on agricultural production can be saved in the data output folder:

save(df, file = "../data/output/df_lp_annual.rda")