<jupyter_start><jupyter_text>Forecasting comparison<jupyter_code>suppressPackageStartupMessages({
    library(tidyverse)
    library(lubridate)
    library(tsibble)
    library(fable)
    library(gtools)
    library(feasts)
})
source("../functions/utils.R")
proj_lc_full <- readRDS("../data/projections/county_final.rds")
acs_1_full <- readRDS("../data/prepped_data/counties_1_acs.rds")
proj_lc <- proj_lc_full %>%
    mutate(lig = log(ig_count_imptd)) %>%
    select(countyid, year, lprd_4, lu90_4, ll90_4, lig)
acs_1 <- acs_1_full %>%
    mutate(lcens = log(chh_)) %>%
    select(countyid, year, lcens)
head(acs_1_full)
proj_data <- proj_lc %>%
    inner_join(acs_1)
print(dim(proj_data))
head(proj_data)
train <- proj_data %>%
    filter(year != 2018)
test <- proj_data %>% 
    filter(year == 2018)
train <- train %>% 
    as_tsibble(index = year, key = countyid)<jupyter_output><empty_output><jupyter_text>Fit TSLM models<jupyter_code>fit <- train %>%
    model(lm = TSLM(lcens ~ trend()))
fc <- fit %>% 
    forecast(h = 1) %>%
    hilo(level = 90) %>%
    unpack_hilo("90%") %>%
    select(-lcens, -year)
head(fc)
res <- test %>%
    left_join(fc)
head(res)<jupyter_output>Joining, by = c("countyid", "year")<jupyter_text>Fit fixed effect autoregressive model<jupyter_code>train <- train %>%
    group_by(countyid) %>%
    mutate(
        lcens_l1 = lag(lcens),
        lcens_l2 = lag(lcens, 2),
        lcens_l3 = lag(lcens, 3),
        lcens_l4 = lag(lcens, 4),
        lcens_l5 = lag(lcens, 5),
        lcens_l6 = lag(lcens, 6)
    ) %>%
    ungroup()
library(fixest)
f1 <- formula(lcens ~ lcens_l1 + lcens_l2 + lcens_l3 
              + lcens_l4 + lcens_l5 + lcens_l6 | countyid)
f2 <- formula(lcens ~ lcens_l1 + lcens_l2 + lcens_l3 
              + lcens_l4 + lcens_l5 + lcens_l6 + lprd_4 | countyid)
m_ar <- feols(f1 , data = train)
m_ar_ig <- feols(f2, data = train)
summary(m_ar)
summary(m_ar_ig)
test_2 <- proj_data %>%
    group_by(countyid) %>%
    mutate(
        lcens_l1 = lag(lcens),
        lcens_l2 = lag(lcens, 2),
        lcens_l3 = lag(lcens, 3),
        lcens_l4 = lag(lcens, 4),
        lcens_l5 = lag(lcens, 5),
        lcens_l6 = lag(lcens, 6)
    ) %>%
    ungroup() %>%
    filter(year == 2018)
head(test_2)
res$m_ar <- predict(m_ar, test_2)
res$m_ar_ig <- predict(m_ar_ig, test_2)
head(res)
rbind(get_errors(res$lcens, res$)
      get_errors(res$lcens, res$.mean, "LSTM projection", round = TRUE),
      get_errors(res$lcens, res$lig, "Raw IG estimate", round = TRUE),
      get_errors(res$lcens, res$lprd_4, "Adjusted IG estimate", round = TRUE),
      get_errors(res$lcens, res$m_ar, "FE model w/ 6 years of lags", round = TRUE),
      get_errors(res$lcens, res$m_ar_ig, "Previous, but add current infogroup", round = TRUE))
options(repr.plot.width=12, repr.plot.height=12)

# split the data by population decile
res <- res %>%
    mutate(acs_quartile = quantcut(lcens, q = 4, na.rm = TRUE))
levels(res$acs_quartile) <- c(1:4)

countyids <- res %>% 
    filter(acs_quartile == 1) %>%
    sample_n(9) %>%
    .$countyid

res_samp <- res %>%
    filter(countyid %in% countyids)
tr_samp <- train %>%
    filter(countyid %in% countyids)

ggplot() +
    geom_line(
        data = tr_samp,
        aes(x = year, y = lcens),
        size = 2
    ) +
    geom_smooth(
        data = tr_samp,
        aes(x = year, y = lcens),
        method = "lm",
        fill = NA, 
        color = 'red',
        linetype = 'dashed',
        size = 2
    ) +
    geom_point(
        data = res_samp,
        aes(x = year, y = lcens),
        size = 2
    ) +
    geom_point(
        data = res_samp,
        aes(x = year, y = .mean),
        color = 'red',
        size = 2,
        position = position_nudge(x = .25)
    ) +
    geom_point(
        data = res_samp,
        aes(x = year, y = lprd_4),
        color = 'darkgreen',
        size = 2,
        position = position_nudge(x = .5)
    ) +
    geom_point(
        data = res_samp,
        aes(x = year, y = lig),
        color = 'blue',
        size = 2,
        position = position_nudge(x = .75)
    ) +
    geom_errorbar(
        data = res_samp,
        aes(x = year, ymin = `90%_upper`, ymax = `90%_lower`),
        color = 'red',
        width = 0.25,
        position = position_nudge(x = .25)
    ) + 
    geom_errorbar(
        data = res_samp,
        aes(x = year, ymin = ll90_4, ymax = lu90_4),
        color = 'darkgreen',
        width = 0.25,
        position = position_nudge(x = .5)
    ) +
    facet_wrap(. ~ countyid) +
    theme_classic()<jupyter_output>`geom_smooth()` using formula 'y ~ x'
{'fraction_non_alphanumeric': 0.1509269356597601, 'fraction_numerical': 0.025081788440567066, 'mean_word_length': 2.599686028257457, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '112237', 'license': 'CC0-1.0', 'n_tokens_mistral': 1971, 'n_tokens_neox': 1796, 'n_words': 353, 'path': 'notebooks/TSLM_projection.ipynb', 'repo_name': 'aridf/infogroup_pop_adj', 'size': '4585'}