Wildfire Smoke and Transportation
November 2, 2026
This instructor-authored guide uses public smoke and transportation data to illustrate a reproducible study.
An illustrative scope
Begin with daily flight disruptions at Northeast airports. A county-quarter analysis of transportation employment and wages can be a separate extension.
Code snippets are templates. Confirm the source schema, sample definition, and required inputs before running them.
Activity measures
Limits on interpretation
These outcomes measure activity, service disruptions, or labor-market conditions. They do not directly measure firm revenue or profit.
An illustrative question
How does exposure to Canadian wildfire smoke relate to flight cancellations and delays at Northeast airports?
Define the comparison before choosing an estimator:
| Element | Example decision |
|---|---|
| Unit | One origin airport on one local calendar day |
| Outcome | Cancelled flights divided by scheduled flights |
| Exposure | County surface-smoke PM2.5 |
| Comparison | Within-airport variation, with a defended adjustment set |
| Claim | A descriptive association, prediction, or causal effect |
Comparisons across modes require compatible scales and observation periods.
Potential policy uses
Additional evidence needed
Connect each proposed policy use to a result that could inform a specific decision.
Organize the literature by mechanism, outcome, and research design. The course bibliography supplies starting points:
| Debate | Start with | Research-design use |
|---|---|---|
| Smoke and labor markets | Borgschulte et al. (2024) | A comparison for a QCEW extension |
| Demand channel = avoidance | Burke et al. (2022); Neidell (2009) | People change behavior when smoke or alerts arrive |
| Flight data in an environmental design | Schlenker and Walker (2016) | How airport operations enter a credible design |
| Why the Northeast may differ | Heutel et al. (2021) | Harm per hot day is smaller where heat is common; smoke may hurt more where it is rare |
The course bibliography contains these references. Support any proposed literature gap with a documented search and study-by-study comparison.
| Borgschulte et al. (2024) | Transportation example | |
|---|---|---|
| Smoke | HMS smoke days, nationwide | Canadian-dominated smoke, Northeast |
| Outcome | County labor-market outcomes (earnings, employment) | Transportation employment and wages, NAICS 48–49 |
| Timing | Smoke days per quarter | Match it, so the two estimates can be compared |
Read the published data and specification sections before constructing a comparable panel.
Short window, COVID period removed
5 years
2018, 2019, 2023, 2024, and 2025 remain. Check how concentrated exposure is in the extreme smoke seasons.
2006–2025, COVID removed
17 years
BTS on-time data, QCEW, and the smoke data all cover the full period. Canadian smoke reached the East before 2018: check 2010 in canadian_dominated_summary.csv.
Excluding March 2020–December 2022 also drops the 2021 episode. State that cost; a robustness check can bring 2021 back.
Record for each source
Check before combining sources
Use the course data guide for bundle definitions and cite the original providers.
1 · The Canadian share is a proxy
Share of same-day fire radiative power within 1,500 km that is in Canada. It ignores winds and smoke height.
Transported smoke takes days to arrive, so compare share_ca_lag03 (days \(d-3\) to \(d\)) with the same-day share. The bundle flag ca_dominated_sameday uses 50%.
2 · Sparse shares are not a complete smoke series
The county-day fire-share table includes only days with positive unrounded CSU smoke.
Read smoke from the full county-day series. A missing share row can accompany zero or negative smoke; it does not justify replacing smoke with zero. Keep unavailable fire shares as NA.
3 · Verify each product’s day definition
Check the BTS flight-date convention, the exact CSU smoke release, and the weather product used in the analysis. For ERA5, verify how hourly UTC observations were aggregated to a daily value.
A local calendar day and a UTC day cover different hours. Document the alignment rule before joining daily observations.
| Code | Cause | Possible design use |
|---|---|---|
| A | Carrier | — |
| B | Weather | Mechanism evidence |
| C | National Aviation System: non-extreme weather, air-traffic-control limits | Mechanism evidence: FAA flow restrictions in smoke |
| D | Security | Placebo outcome: smoke should not move it |
Use all cancellations and delays as the outcome
Smoke-related cancellations need not be coded B; carriers may record them under C. Use the codes as evidence on the mechanism, not to define the outcome.
Verify the reporting-carrier threshold for each sample year: the set of carriers that report can change over time.
01_bts_download.R
Monthly PREZIP files; keep Northeast flights as you read
→
02_bts_clean.R
Airport-day panel; grain check
→
03_merge_smoke.R
Airport → county; CSU smoke, Canadian share, weather
→
04_audit.R
Rows, keys, duplicates, missing values, ranges → a table in the report
library(tidyverse)
ne_states <- c("CT", "ME", "MA", "NH", "RI", "VT", "NJ", "NY", "PA")
# One unzipped monthly PREZIP file; column names as in TranStats
read_bts_month <- function(path) {
read_csv(path, show_col_types = FALSE,
col_select = c(FlightDate, Origin, OriginState,
Cancelled, CancellationCode)) |>
filter(OriginState %in% ne_states)
}
airport_day <- list.files("data/raw/bts", pattern = "\\.csv$",
full.names = TRUE) |>
map(read_bts_month) |>
list_rbind() |>
summarise(n_sched = n(),
n_cancel = sum(Cancelled),
n_wx = sum(CancellationCode %in% "B"),
n_nas = sum(CancellationCode %in% "C"),
n_sec = sum(CancellationCode %in% "D"),
.by = c(Origin, FlightDate)) |>
mutate(cancel_rate = n_cancel / n_sched)
stopifnot(!anyDuplicated(airport_day[c("Origin", "FlightDate")])) # one row per airport-day%in% treats the blank cancellation code of an operated flight as “no match,” where == would return NA.
library(sf)
ne_fips <- c("09", "23", "25", "33", "44", "50", "34", "36", "42")
counties <- st_read("wildfire_data/clean/counties_2020.gpkg", quiet = TRUE) |>
filter(census_region == "Northeast")
# airports: one row per Origin with lon and lat (BTS airport master coordinates)
airports_sf <- st_as_sf(airports, coords = c("lon", "lat"), crs = 4326)
airport_county <- airports_sf |>
st_join(select(counties, GEOID), join = st_nearest_feature) |>
st_drop_geometry()
read_csu <- function(f) {
read_csv(f, col_types = cols(GEOID = col_character())) |>
filter(str_sub(GEOID, 1, 2) %in% ne_fips) |>
select(GEOID, date, smoke_pm25)
}
csu <- list.files("wildfire_data/clean/csu_smoke_county_day", full.names = TRUE) |>
map(read_csu) |>
list_rbind()
panel <- airport_day |>
left_join(airport_county, by = "Origin") |>
left_join(csu, by = join_by(GEOID, FlightDate == date))st_nearest_feature catches coastal airports that fall just outside the shoreline-clipped county outlines.
Create an audit table from each input actually read:
audit <- function(df, key, date_col) {
tibble(rows = nrow(df),
dup_keys = sum(duplicated(df[key])),
rows_w_na = sum(!complete.cases(df)),
first = min(df[[date_col]]),
last = max(df[[date_col]]))
}
bind_rows(
airport_day = audit(airport_day, c("Origin", "FlightDate"), "FlightDate"),
csu = audit(csu, c("GEOID", "date"), "date"),
.id = "file"
)COVID is a design decision, not a check
Changes in transportation demand during the pandemic can affect the comparison. Pre-specify whether and how to exclude a pandemic period, explain the rationale, and compare sensitivity to that choice.
\[ \text{smoke}_{ct} = \text{HMS flag}_{ct} \times \left(\text{total PM2.5}_{ct} - \text{background}_{ct}\right) \]
Begin with an episode plot and a clearly defined binned comparison. Removing weather or calendar effects is a modeling choice that needs explanation.
Figure A · One episode
If the visibility mechanism is real, it should show here first.
Figure B · All airport-days
Treat these as illustrative bins. Report the observations in each and check sensitivity to the cut points.
library(fixest)
fig_b <- panel |>
filter(!between(FlightDate, as.Date("2020-03-01"), as.Date("2022-12-31"))) |>
mutate(month = month(FlightDate),
wday = wday(FlightDate, label = TRUE),
smoke_bin = cut(smoke_pm25, c(-Inf, 0, 5, 15, 35, Inf),
labels = c("<= 0", "0-5", "5-15", "15-35", "> 35")))
# Remove airport-by-month and weekday means
fit0 <- feols(cancel_rate ~ 1 | Origin^month + wday, data = fig_b)
fig_b$resid <- resid(fit0, na.rm = FALSE)
fig_b |>
filter(!is.na(smoke_bin)) |>
summarise(mean_resid = mean(resid, na.rm = TRUE), n = n(), .by = smoke_bin) |>
ggplot(aes(smoke_bin, mean_resid)) +
geom_col() +
labs(x = "Smoke PM2.5 (µg/m³)",
y = "Cancellation rate relative to\nairport-month and weekday means")Report n for each bin in the caption: the top bins will hold few airport-days, and readers need to know that.
Northeast counties shaded by Canadian-dominated days above 15 µg/m³ (2006–2025), with airports sized by flights. Classwork 5 covers the joins.
ca_days <- read_csv("wildfire_data/clean/county_day_fire_share_2006_2025.csv.gz",
col_select = c(GEOID, smoke_pm25, ca_dominated_lag03),
col_types = cols(GEOID = col_character())) |>
filter(str_sub(GEOID, 1, 2) %in% ne_fips, smoke_pm25 > 15, ca_dominated_lag03) |>
count(GEOID, name = "ca_days_gt15")
flights_by_airport <- summarise(airport_day, flights = sum(n_sched), .by = Origin)
counties |>
left_join(ca_days, by = "GEOID") |>
mutate(ca_days_gt15 = replace_na(ca_days_gt15, 0)) |> # here, no row = zero days
ggplot() +
geom_sf(aes(fill = ca_days_gt15), color = "white", linewidth = 0.1) +
geom_sf(data = inner_join(airports_sf, flights_by_airport, by = "Origin"),
aes(size = flights), shape = 21, fill = "white") +
scale_fill_viridis_c(name = "Canadian-dominated\ndays above 15 µg/m³") +
theme_void()\[ \text{CancelRate}_{at} = \beta\, \text{Smoke}_{c(a),t} + W_{c(a),t}'\gamma + \alpha_{a \times m(t)} + \delta_{\text{dow}(t)} + \lambda_{y(t)} + \varepsilon_{at} \]
\(a\): airport · \(t\): day · \(c(a)\): the airport’s county · \(W\): precipitation, temperature, wind · \(\alpha_{a \times m}\): airport-by-month-of-year effects · \(\delta\): day-of-week (and holiday) effects · \(\lambda\): year effects
Unit and outcome
fepois(), offset = ~ log(n_sched))Sample and extension
Within an airport, over time
Lecture 7’s problem: shocks persist within a unit → cluster by airport
Across airports, on the same day
Smoke is a regional shock: on June 7, 2023, nearly every Northeast airport was treated at once → cluster by date
# panel now also holds weather (tmax, prcp, wind) and month, wday, year
m1 <- feols(cancel_rate ~ smoke_pm25 + tmax + prcp + wind |
Origin^month + wday + year,
data = panel, cluster = ~ Origin + FlightDate)
m2 <- fepois(n_cancel ~ smoke_pm25 + tmax + prcp + wind |
Origin^month + wday + year,
offset = ~ log(n_sched),
data = panel, cluster = ~ Origin + FlightDate)
etable(m1, m2)Clustering by state leaves only nine clusters, too few for reliable standard errors. fixest runs the same within estimator as plm in Lecture 7, but handles many fixed effects and two-way clustering.
Candidate controls: precipitation, temperature, and wind. Explain their proposed roles in the causal model.
For a total-effect estimand: avoid conditioning on a variable that mediates smoke’s effect. Justify the role of visibility and distinguish it from other weather controls.
smoke_episodes_2006_2025.csv already lists every regional episode with its dates and a canadian_dominated flag. Days −7 to −2 illustrate a lead test; days 0 to 7 describe the post-episode response.
episodes <- read_csv("wildfire_data/clean/smoke_episodes_2006_2025.csv") |>
filter(region_type == "Census region", region == "Northeast",
canadian_dominated)
windows <- episodes |>
reframe(FlightDate = start_date + (-7:14), rel_day = -7:14,
.by = episode_id)
es_panel <- panel |>
inner_join(windows, by = "FlightDate", relationship = "many-to-many")
es <- feols(cancel_rate ~ i(rel_day, ref = -1) | Origin^episode_id + wday,
data = es_panel, cluster = ~ Origin + FlightDate)
iplot(es)Two episodes less than three weeks apart produce overlapping windows: trim them, or keep each day with its nearest episode, and say which.
Effect heterogeneity
Prediction
| Check | Purpose | Limitation to examine |
|---|---|---|
| Leads before an episode | Look for changes preceding measured exposure | A one-day lead may still lie inside a persistent smoke event |
| A justified placebo outcome | Probe a mechanism or design assumption | Explain why smoke should not affect the selected outcome |
| Hold out a state | Check geographic sensitivity | Exposure and flights may be concentrated in a few states |
| Hold out a year | Check temporal sensitivity | Extreme seasons may differ from the training years |
Pre-specified checks
State what result would challenge the interpretation before inspecting the estimate.
Exposure and outcomes
Coverage and sample
Phase 1
BTS only. Download, filter to Northeast airports, build the airport-day panel, merge smoke and the Canadian share. Figures A and B; Map 1.
Phase 2
Baseline. Two-way-clustered regressions; the event study on the episode list; lead and placebo-outcome tests.
Phase 3
QCEW. NAICS 481, 484, 485 by county-quarter, with BMZ (2024) timing. Robustness with HMS and Childs smoke.
Phase 4
Extend and write. Heterogeneity interactions, the ML prediction comparison, a first draft of results and limitations.
Complete and check one primary outcome before adding a new mode, frequency, or data source. The phases are examples, not deadlines.
Scope
Which single outcome can the available data measure and the proposed design support?
Mechanism
What additional evidence could distinguish a visibility mechanism from changes in demand?
Feasibility
Which input, definition, or join is the main feasibility risk for the first panel?