# A Machine Learning Approach to Measuring Climate Adaptation

Max Vilgalys111<EMAIL_ADDRESS>Eastern Research Group, Inc.

###### Abstract

I measure adaptation to climate change by comparing elasticities from short-
run and long-run changes in damaging weather. I propose a debiased machine
learning approach to flexibly measure these elasticities in panel settings. In
a simulation exercise, I show that debiased machine learning has considerable
benefits relative to standard machine learning or ordinary least squares,
particularly in high-dimensional settings. I then measure adaptation to
damaging heat exposure in United States corn and soy production. Using rich
sets of temperature and precipitation variation, I find evidence that short-
run impacts from damaging heat are significantly offset in the long run. I
show that this is because the impacts of long-run changes in heat exposure do
not follow the same functional form as short-run shocks to heat exposure.

##### Keywords:

Climate change, directional derivative, machine learning, panel data

##### JEL Classification:

C14, C33, C55, Q51, Q54

## 1 Introduction

Measuring adaptation to recent climate change is important for informing
climate policy and projecting damages of future climate change. The extent of
recent adaptation can signal when policy interventions are needed and give a
sense of how much adaptation will be possible to projected changes.
Researchers can study this by modeling how weather shocks impact economic
outcomes, and examining how impacts change with more exposure to the weather
shock (e.g. [5]) or over time (e.g. [2]). These studies depend on accurate
models of the relationship between weather and economic outcomes.

Machine learning (ML) can help model these relationships when researchers have
rich data, but lack expert guidance to suggest a functional form. When domain
experts provide a model of how weather shocks impact economic outcomes,
researchers can use this to accurately measure and compare impacts. Otherwise
the researcher must learn this relationship from data. Economists typically
use classic statistical tools to model weather-economic relationships without
imposing strong functional form assumptions [18]. Such tools work well with
low-dimensional weather variation, but can lead to high variance or
inconsistent estimates as the dimensionality of covariates increases. Even a
model with only temperature and precipitation can become high-dimensional if a
researcher flexibly models each variable and the interactions between them. ML
is well suited for flexible modeling in such settings [27], and can help
researchers take full advantage of high-dimensional weather variation in
modern panel datasets.

I introduce an ML approach to study adaptation to damaging heat exposure in
United States (U.S.) corn and soy production. High temperatures are generally
damaging for crop growth, although adaptation may offset some of these damages
in the long run. [33] introduce a parsimonious model of how an annual shock of
heat exposure impacts crop yield. They model crop yields as a piecewise-linear
function of heat exposure below and above a crop-specific temperature
threshold, where heat exposure below the threshold is beneficial and heat
exposure above the threshold is damaging. By estimating the parameter on
damaging heat using different sources of variation, [32], [5], and [24] argue
that the observed degree of adaptation to damaging heat exposure will not be
sufficient to offset projected climate damages.

I measure the degree of adaptation after learning the crop yield-weather
relationship from data. I compare impacts from short-run and long-run changes
in damaging heat exposure. With a low-dimensional linear model, this is
equivalent to the method from [5]. The approach allows me to incorporate more
flexible models and high-dimensional weather variation, making the method
suitable for applications with rich weather variation and little expert
guidance on how that variation impacts economic outcomes. I use ML to model
the crop yield-weather relationship, specifically Least Absolute Shrinkage and
Selection Operator (Lasso) and a neural network (NNet). My ML models account
for additive fixed effects while flexibly modeling temperature, precipitation,
and interactions between them.

To estimate the degree of adaptation, I estimate the elasticity of crop yields
with respect to damaging heat exposure. This elasticity summarizes how
damaging heat exposure impacts crop yields, and can be computed for each
model. There is no single parameter I can compare across models, as in [5].
Instead, I find the elasticity by fitting a regression function of log crop
yields on temperature and precipitation and computing the average directional
derivative in the direction of a marginal increase in damaging heat exposure.
I estimate these regression functions via ordinary least squares (OLS) and ML.

I implement a debiasing procedure to address bias from standard ML models. ML
approaches can induce bias from overfitting or regularization, but there are
approaches to reduce this bias when estimating causal parameters or statistics
based on regression functions [7, 8]. I adapt the estimator from [8]. This
approach uses double machine learning (DML), where standard ML estimates are
corrected using a second ML algorithm. I adjust the second ML algorithm to
suit the panel setting. I then apply this DML procedure to debias the
elasticity estimates.

I compare elasticities to estimate the degree that short-run impacts from
damaging heat exposure are offset in the longer term. I construct panel
datasets of long-run and short-run variation in crop yield and weather
variables from 1990-2019. I then estimate the elasticity of crop yield with
respect to damaging heat exposure in each dataset. I compute this elasticity
using OLS, DML, and ML without bias correction. Comparing these estimates, I
find the extent that short-run impacts from damaging heat exposure are offset
in the longer term.

Before taking this approach to the data, I conduct a simulation exercise to
evaluate DML estimates of an elasticity relative to OLS and standard ML
approaches. Each simulation trial uses the empirical distribution of
temperature and precipitation from U.S. counties, but simulates outcome
variables. I examine the performance of my DML procedure relative to OLS and
ML without bias correction, comparing the bias and variance of recovering the
true elasticity. I use three potential sets of weather variables: a simple,
commonly used set of annual temperature and precipitation variables, a set
that includes richer variation in annual temperature, and a set that includes
monthly observations of precipitation and the rich variation in temperature.

This simulation exercise clearly highlights the benefits of using debiased
machine learning for high-dimensional settings. Both NNet and Lasso have
advantages over OLS in high-dimensional cases, although NNet has lower bias.
With the rich set of annual temperature variation, DML estimates have
significantly lower variance than OLS and lower bias than standard ML
estimates. With the set of monthly temperature and precipitation, DML
estimates have significantly less bias and variance than OLS, and lower bias
than standard ML.

I then apply the approach to a dataset of U.S. crop yields, and find evidence
that adaptation is offsetting impacts from damaging heat exposure. I implement
Lasso, NNet, and OLS estimators on a county-level dataset of corn and soy
yields from 1990 to 2019. I compare estimates of the elasticity using short-
run and long-run variation. I take 500 bootstrap trials of the estimated
elasticity, using the three sets of weather variation as in the simulation
exercise. Using the simple annual set of weather variables as in [32] and [5],
I find that there has been little to no significant adaptation to climate
change in corn or soy production. This confirms the results from [5].

Using more flexible sets of weather variation, I find that a large share of
short-run impacts from damaging heat exposure are offset in the long run. With
short-run variation, I find statistically and economically significant
declines in yield from a marginal increase in damaging heat exposure. However,
I do not find evidence of such declines when using long-run variation. These
results hold for both corn and soy. The primary difference comes from using a
richer set of weather variables. I make this same conclusion using OLS with
the more flexible set of weather variation, although the DML approach results
in smaller confidence intervals. This shows that a substantial degree of the
short-run impacts from damaging heat exposure are offset in the long run,
suggesting substantial adaptation to this heat exposure.

This result differs dramatically from the conclusions by [5] and other
analyses [32, 24]. This is likely explained by model misspecification for the
impact of a long-run shift in heat exposure on crop yields. I show that in the
panel with long-run variation, the simple model from [33] does not adequately
summarize the flexible role of temperature variables. While damaging heat
exposure is correlated with declines in crop yield, other temperature
variation is better able to explain these declines. This suggests that there
is limited adaptation to some damaging feature of climate change, but not to
the specific feature of marginal increase in damaging heat exposure.

This paper is related to several literatures. First is a literature on
estimating the degree of adaptation to climate change. Measuring adaptation to
climate change requires understanding how weather influences economic
outcomes. For a review of economics literature on measuring the economic
impacts of the weather, see [12]. [18] provides an overview of econometric
approaches to measuring these impacts. Much of this literature focuses on
agriculture, as this sector is directly exposed to weather and hence is
particularly vulnerable to climate change [35]. The first approach to studying
impacts of climate change used the Ricardian approach, where researchers
compare the value of agricultural land in cross sections. [26] forecast the
impacts of climate change by regressing average temperature and agricultural
property value in a cross section of U.S. counties. This approach is
susceptible to omitted variable bias, and subsequent work has focused on
addressing specific omitted variables such as endogenous changes in farmer
technology [23] or nonfarm income [30]. Other approaches to estimate the
potential for adaptation in agriculture involve economy-wide simulations [10],
production changes in historical migrations [37, 28], natural experiments [16,
15], or panel approaches.

Panel approaches address omitted variable bias by identifying adaptation from
annual or long-term variation within a panel dataset. [32] uses panel data to
estimate the elasticity of crop yields with respect to extreme heat exposure,
and conclude that there is limited potential to adapt to climate change
because these damages are similar in the southern and northern U.S. despite
climatic differences. [2] use a flexible model of temperature exposure to
document how the mortality consequences of extreme heat declined over the
20${}^{\text{th}}$ century. [5] show that panel variation can be used to
estimate adaptation to recent climate change by using separate sources of
variation to identify the impacts of weather shocks and shifts in average
temperature. [24] provides an alternate approach that partially identifies the
degree of possible adaptation by considering the role of ex-ante and ex-post
adaptation to heat exposure shocks.

My paper is most closely related to [5]. Like their paper, I estimate the
degree that damages to corn and soy yields from short-run changes in weather
are offset over longer exposures. I also use crop and weather data from U.S.
agriculture. My approach differs because I consider richer sets of weather
variables, and use DML to model learn the relationship between these data and
crop yields. I conclude that there has been a higher degree of adaptation to
damaging heat exposure.

Second is a growing literature on applying ML methods in economics. [21]
discuss applications of predictive machine learning in economics, and [38] and
[27] provide a practical guide to algorithms. Several recent papers have used
ML to measure important outcomes in environmental economics. [11] proposes a
semi-parametric NNet and uses it to study the impact of climate change on corn
yields. [13] uses a ML approach to measure the costs of air pollution. [6] use
ML to refine estimates of energy efficiency improvements. [36] use a DML
approach to measure effectiveness of an agricultural intervention, and [22]
introduce a DML approach to measure elasticities in a panel setting. There are
also numerous applications within agriculture; for a review, [25].

My paper is most related to [11] and [22]. [11] estimates a NNet that accounts
for unobservable county-level fixed effects, and uses the model to predict
yield under counterfactual climate change scenarios. I use the same technique
to address county-level fixed effects, although I modify the algorithm in
order to recover derivatives from the network and to ensure that the network
is differentiable. [11] studies the impacts of future climate change, while I
use this tool to understand adaptation to recent climate change. Like [22], I
estimate the elasticity of crop yield with respect to an increase in damaging
heat exposure. I consider more flexible representations of temperature
variation, and I apply the estimator to measure adaptation to recent climate
change.

Within the literature on machine learning, this paper applies results from the
emerging field of DML. [7] prove that sample splitting and constructing
Neyman-orthogonal moment conditions can yield approximately debiased machine
learning estimates in certain settings. [34] extend the Neyman-orthogonal
moment condition approach to several other statistical targets, including
structural derivatives. In an alternate approach, [8, 9] give an approximately
debiased estimator for a more general class of linear functionals based on the
Riesz representation theorem. [22] adopt this approach in panel settings and
prove asymptotic normality of the estimator for the average derivative. Like
[22], my approach applies the result from [8] in panel settings; I use a
different approach to address fixed effects and consider a more flexible
representation of temperature variation.

The rest of this paper proceeds as follows. In Section 2, I describe the data
used for this project. I illustrate the degree of climate change and describe
the transformations required to generate the growing degree days. In Section
3, I describe the methods used; this includes details on how to compute
average derivatives and an explanation of the debiased machine learning
estimation approach. In Section 4, I give details and results of the
simulation exercise. In Section 5, I present and discuss results from using
the estimation procedure to measure the degree of adaptation to climate
change. Section 6 concludes.

## 2 Data

For the empirical application, I use weather and crop data from U.S. corn and
soy production from 1990-2019. I consider counties east of the 100°West
meridian, which defines an agricultural region of significant corn and soy
cultivation. From 1990-2019, this region produced over 93% of the nation’s
corn and over 99% of the nation’s soy. Crop data are annual yield (bushels per
acre) of corn and soy, from the U.S. Department of Agriculture’s Survey of
Agriculture. These data also include the area planted (acres) in each county.

Weather data are generously shared by [32], who provide a gridded dataset of
daily temperature and precipitation from a network of consistently reporting
weather stations. As in [32] and [29], I consider weather during the March-
August growing season. I aggregate the gridded dataset to a county-level
dataset of daily maximum temperature, minimum temperature, and precipitation.
I then transform daily temperature exposure into growing degree days (GDD) at
a monthly level and for the growing season. Figure 1 illustrates the
transformation from daily temperature observations to GDD.

(a)

(b)

(c)

(d) Monthly Flexible

(e) Yearly Flexible

(f) Yearly Linear

Figure 1: Transformations from daily minimum and maximum temperature records
into the weather variables used in the analysis. (a) shows daily temperature
within one day. In (b), this temperature is translated into growing degree
days (GDD). (c) illustrates how daily GDD observations are aggregated into
cumulative GDD exposure. (d) is the Monthly Flexible variable, where daily
temperature is aggregated into GDD in each temperature bin for each month. In
(e), this is aggregated into total growing season heat exposure in each
temperature bin. In (f), this is further aggregated into total growing season
heat exposure above and below the dotted line, 29°C.

Throughout this paper, I consider three sets of weather variables per county:
total growing season heat exposure above and below 29°C plus total growing
season precipitation (Yearly Linear), total growing season heat exposure in
each 1°C temperature bin plus total growing season precipitation (Yearly
Flexible), and monthly heat exposure in each 1°C bin and monthly precipitation
for each month of the growing season (Monthly Flexible). Figure 1 illustrates
the transformations of temperature variables. The Yearly Linear
transformation, shown in Figure 1(f), is widely used in economic analysis.
[32] demonstrate that, in a panel setting, regression based on heat exposure
above and below a crop-specific damaging threshold explains as much variation
in yields as more flexible models. [5] identify 29°C as the threshold for
damaging heat exposure in the long run, for both corn and soy. I use the 29°C
threshold throughout the analysis, and refer to heat exposure above 29°C as
damaging heat exposure.

(a) Corn yields

(b) Soy yields

(c) GDD above 29°C

Figure 2: Differences in county-level values of corn yields, soy yields, and
growing-season GDD above 29°C, between 1990-1999 averages and 2010-2019
averages. Maps include counties east of the 100 °W meridian. In (a) and (b),
the first color bin includes all counties with yield declines (roughly 4% of
counties) and the remaining distribution is divided evenly among the remaining
bins. In (c), ‘X’ (‘O’) hatching indicates counties with an increase
(decrease) in heat exposure above 29°C.

Variation in the degree of climate change between counties is used to estimate
the long-run elasticities in our sample. Figure 2(c) shows illustrates the
difference in GDD exposure above 29°C per US county, from the period 1990-1999
to the period from 2010-2019. Heat exposure increases in 80% of counties in
the sample, but there is substantial variation in the degree of warming even
within states. There are many pairs of neighboring counties where one
experienced an increase in damaging heat exposure and the other experienced a
decrease. This variation is plausibly uncorrelated with unobservable factors
such as soil quality or per-state fixed effects. See [5] for a detailed
argument that this type of variation can identify the effects of climate
change.

While many counties experienced damaging heat exposure increases and
decreases, only a handful had average yields decline between these time
periods. Figure 2(a) and Figure 2(b) show changes in average yields of corn
and soy over this time period. Only 4.03% of counties saw corn yields decline
and only 3.96% saw soy yields decline, among counties that grew corn or soy in
both 1990-1999 and 2010-2019. This suggests that improved farming technology
increased yields throughout the sample, and that these improvements exceeded
damages from increased damaging heat exposure.

## 3 Methods

In this section, I discuss the methods used to estimate the degree of
adaptation. I first define adaptation in terms of average directional
derivatives, and describe how to measure these average directional
derivatives. I then introduce the ordinary least squares (OLS), Lasso and
neural network (NNet) procedures, including details on how to account for
fixed effect terms in each model and the cross-folds training procedures. I
also describe the double machine learning (DML) approach used to adjust for
bias in the standard machine learning (ML) estimates.

### 3.1 Adaptation

I define adaptation as the amount of short-term impact to yield from damaging
heat exposure that is offset in the longer term. This definition encompasses
all adaptation behaviors an agent makes to their production technology,
management practices, or variety choice within each crop as they are exposed
to climate change. It does not capture some other important margins of
adaptation such as crop switching or exit from agriculture. [5] provide
evidence that these margins of adaptation are limited. This definition also
does not capture changes that alter the impact of damaging heat independently
from an individual’s exposure to climate change, such as the decline in heat-
related mortality as studied by [2]. First I define the estimation target, as
introduced by [5]. I generalize this definition in terms of elasticities, and
discuss how to estimate those elasticities using flexible functional forms
and/or richer sets of weather variation.

Following [5], I measure adaptation by comparing the elasticity of crop yields
with respect to extreme heat using short-run and long-run variation. Short-run
variation comes from year-to-year changes in an annual panel of weather
observations and log crop yields. To capture long-run variation, I first
average weather and crop yield data over a long period (I use 10-year periods)
and then construct a two-period panel. The elasticity computed using short run
variation captures the extent that a weather shock of extreme heat in a single
year impacts yields, while the elasticity computed using long run variation
captures the extent that exposure to a long period of increased heat exposure
will impact average crop yields. Estimates using long-run variation arguably
capture the extent of damages from climate change, because there is a change
in average heat exposure over a relatively long period where farmers have time
to adjust to those changes. Figure 2 illustrates long-run variation in the
sample, comparing the differences in long-run average crop yields and damaging
temperature exposure over the decades I study. Comparing impacts using these
sources of variation, I can conclude whether the short-term damages are offset
in the long term.

It is straightforward to recover these elasticities in a model where log yield
is linear in total growing season heat exposure. [5] use such a model, which I
adapt below:

$y_{it}=a_{i}+\beta_{1}lower_{it}+\beta_{2}higher_{it}+g(prec_{it})+\varepsilon_{it}$
(1)

where $higher_{it}$ ($lower_{it}$) is total growing season heat exposure above
(below) the damaging temperature threshold, $g$ is some function of
precipitation, $a_{i}$ is an additive per-county fixed effect term, and
$\varepsilon_{it}$ is an additive error term. [5] estimate this equation with
OLS, after using the within transformation to remove the fixed effect
term222[5] advocate using first differences to remove the fixed effect in the
long-run variation panel; I use a within transformation approach as this is
numerically equivalent to taking first differences in a two-period panel. .
The key parameter here is $\beta_{2}$, which captures the extent that log crop
yield changes with marginal increase in damaging heat exposure. This is
equivalent to the elasticity of crop yield with respect to damaging heat
exposure. Note that $\beta_{2}$ is expected to be negative in the short run,
as temperature shocks in this range are damaging to crop growth [33, 32]. The
estimate of the share of short-run damages that are offset in the longer term
is therefore
$(\hat{\beta}_{2}^{SR}-\hat{\beta}_{2}^{LR})/\hat{\beta}_{2}^{SR}=1-\hat{\beta}_{2}^{LR}/\hat{\beta}_{2}^{SR}$,
where $\hat{\beta}_{2}^{SR}$ ($\hat{\beta}_{2}^{LR}$) is the estimate using
the short-run (long-run) variation. [5] take bootstrap samples of this ratio
and fail to reject the null hypothesis that the ratio is different from 0.

To recreate this ratio with a more flexible functional form, I replace the
parameter estimates with average directional derivatives. Consider a more
general functional form:

$y_{it}=a_{i}+\gamma(X_{it})+\varepsilon_{it}$ (2)

Here, $X_{it}$ is a collection of weather variables, $\gamma$ is a general
function, $a_{i}$ is an additive fixed effect term, and $\varepsilon_{it}$ is
an additive error term. I use three different collections of weather variables
for $X_{it}$, as described in Figure 1. As $y_{it}$ is the log of crop yields,
the elasticity of crop yield with respect to some variable is equivalent to
the average directional derivative of $\gamma$ with respect to that variable.

The analogue to $\beta_{2}$ from Equation 1 is therefore the average
directional derivative of $\gamma$ with respect to $higher_{it}$. Let
$\theta^{SR}$ ($\theta^{LR}$) be this average directional derivative using
short-run (long-run) variation. I then take bootstrap samples of the ratio
$1-\hat{\theta}_{2}^{LR}/\hat{\theta}_{2}^{SR}$ to test the hypothesis that
damages from short-run heat exposure are offset in the longer run.

This average directional derivative is equivalent to an average partial
derivative when using the Yearly Linear set of weather variables.
Specifically, take
$\hat{\theta}=\mathbb{E}[\partial\hat{\gamma}(X_{it})/\partial higher_{it}]$.
When using a linear specification, this average derivative is equivalent to
the parameter estimate from OLS. When using NNet, this can be recovered after
training the network (see Section 3.3.b). When the specification involves
basis functions (such as OLS with polynomial basis functions or Lasso), I
recover this derivative by taking the average of the dot product of the
gradient of the basis functions and the estimated coefficients. I describe
this procedure in Section 3.2.

When I use another set of weather variation, it is also necessary to account
for the extent that each weather variable contributes to the total growing
season heat exposure. I find this using the chain rule:

$\hat{\theta}=\sum_{X_{it}\in
X_{it}^{higher}}\frac{\partial\hat{\gamma}(X_{it})}{\partial
X_{it}}\frac{\partial X_{it}}{\partial higher_{it}}$ (3)

Where $X_{it}^{higher}$ is the set of weather variables that are summed to
reach total growing season heat exposure above the temperature threshold. I
set $\partial X_{it}/\partial higher_{it}=X_{it}/higher_{it}$; this captures
the assumption that additional marginal heat will be distributed
proportionally to the empirical heat exposure distribution.

This allows me to measure adaptation to climate change from recent panel
variation while using a flexible model of high-dimensional weather variation.
In the following sections, I describe how I estimate the regression function
and average directional derivative using various estimation techniques.

### 3.2 Ordinary Least Squares methods

Ordinary Least Squares (OLS) is a classical statistics approach to estimating
this elasticity. I use methods with a linear functional form (OLS Linear), and
after applying a basis function transformation of polynomial functions and
interactions (OLS Poly). In both cases, I use the within transformation to
remove a county-level fixed effect term, and then include yearly fixed effect
terms via dummy variables.

I use basis functions that include interactions and flexible functional forms
of the original data. I specify polynomial basis functions of all terms, as
well as interactions between these polynomial expansions. To produce a
tractable model, I limit the space of potential interactions to interactions
between heat exposure and precipitation variables within the same time period.
For example, in the Monthly Flexible specification, I consider cumulative GDD
in July between 28 and 29 C interacted with precipitation in July, as well as
squared values of each term and the interactions between those polynomial
expansions, but do not consider interactions of that variable with cumulative
GDD in any other temperature bin, or precipitation in any other month. I use
3rd order polynomials for the Yearly Linear and Yearly Flexible variable sets,
and 2nd order polynomials for the Monthly Flexible variable set. I then scale
each flexible basis function so that it has mean zero and variance 1.

For OLS Linear, I use the identity set of basis functions; that is,
$b(X_{it})=X_{it}$. For Yearly Linear, $X_{it}$ has 3 covariates; for Yearly
Flexible $X_{it}$ has 41; and for Monthly Flexible $X_{it}$ has 246. For OLS
Poly, I use the set of basis functions described above. For Yearly Linear,
$b(X_{it})$ has 30 covariates; for Yearly Flexible $b(X_{it})$ has 486; and
for Monthly Flexible $b(X_{it})$ has 1464.

OLS assumes that the following is a true model of the relationship:

$y_{it}=a_{i}+b(X_{it})\beta_{0}+\varepsilon_{it}$ (4)

As is common in economics, I assume that we have relatively short panels where
it is not possible to consistently estimate $a_{i}$ by including dummy
variables. I therefore use the within transformation to remove county-level
fixed effect terms:

$\ddot{y}_{it}=\ddot{b}(X_{it})\beta_{0}+\varepsilon_{it}$ (5)

where the double dot denotes the within transformation, i.e.
$\ddot{y}_{it}:=y_{it}-\text{mean}(\\{y_{it}\forall t\\})$ and
$\ddot{b}(X_{it}):=\text{mean}(\\{b(X_{it}\forall t\\})$. Note that to
construct $\ddot{b}(X_{it})$, the mean of all observations in that panel unit
is subtracted after applying the basis function transformation. This ensures
that $\beta_{0}$ is the same parameter vector between models. I use OLS to
estimate $\hat{\beta}$.

I then compute the average directional derivative by projecting my estimate of
$\beta_{0}$ on partial derivatives of the basis functions. Define the gradient
of the basis function as $b_{higher}$; this is a 1 by $p$ dictionary of the
derivative of each basis function with respect to $higher_{it}$. The true
average directional derivative is then
$\theta_{0}=\mathbb{E}[b_{higher}(X_{it})\beta_{0}]$ and its estimate is
$\hat{\theta}=\mathbb{E}[b_{higher}(X_{it})\hat{\beta}]$.

In specifications in my main analysis, I include per-year fixed effects as
dummy variables. I assume that there are enough observations per time period
to consistently estimate these variables separately. The derivative of each
per-year fixed effect is zero, so including these terms does not change how I
estimate the average directional derivative.

### 3.3 Machine Learning Methods

Machine Learning (ML) methods include a range of estimation techniques that
allow consistent function approximation in high-dimensional settings, when the
number of covariates is large relative to the number of observations. Such
settings poses challenges for classical statistical methods such as OLS,
binning, or kernel regression. I use Lasso and Neural Networks (NNets),
although the same procedure could be used for another algorithm such as random
forests, support vector machines, or other methods. I focus on these two
machine learning algorithms because each enables the researcher to incorporate
linear fixed effects and to evaluate derivatives without numerical
differentiation. Standard ML methods can induce bias in regression analysis; I
overcome this bias by using a procedure from [8].

The average derivative is computed from a ML regression of the output variable
on weather inputs. Write $\gamma(\cdot;\lambda)$ to denote the flexible
machine learning function, emphasizing the dependence on the hyperparameter
$\lambda$. The hyperparameter is a researcher-specified value that influences
the behavior of the model, such as the regularization penalty in Lasso or the
network width in NNet. The machine learner is estimated as a regression
function:
$\mathbb{E}[\ddot{y}_{it}|\ddot{X}_{it}]=\hat{\ddot{\gamma}}(X_{it};\lambda)$.
Section 3.3.a and Section 3.3.b give details on each estimation procedure. Let
$m(\gamma,X_{it};\lambda)$ denote the directional derivative of
$\gamma(\cdot;\lambda)$ evaluated on observation $X_{it}$.

I use the double machine learning (DML) procedure from [8] to find an
approximately debiased estimate of the average directional derivative. I
discuss this procedure in more detail in Section 3.3.c. Briefly, I estimate a
second machine learner and use this to construct an estimate of the average
directional derivative that is robust to errors in estimating either the first
or second machine learner. Let $\alpha(X;\kappa)$ denote this second machine
learner, emphasizing the dependence on hyperparameter $\kappa$. [8] show that
the expression
$\mathbb{E}[m(\hat{\gamma},X_{it};\lambda)+\hat{\alpha}(X_{it};\kappa)\hat{\varepsilon}_{it}]$
is an approximately unbiased estimate of the true average directional
derivative, where $\hat{\varepsilon}_{it}$ is the residual from estimating
$y_{it}$. I modify the form of the doubly robust estimator in [8] to account
for the panel structure of the data; details are in Section 3.3.c.

I use a data-driven process to determine the value of the hyperparameters.
First split each panel unit (a U.S. county) into one of the $k$ folds for
cross validation. Grouping observations from each panel unit into the same
fold reduces correlation between training and test data, as counties share
unobservable characteristics that likely influence the distribution of weather
and crop yields. Let $\mathcal{I}_{\ell}$ denote the set of indices in fold
$\ell$, for $\ell\in\\{1,2,\dots,k\\}$. For each of the $k$ folds, I train the
ML and DML algorithm on data not in fold $\ell$, and evaluate the algorithm
only on indices in fold $\ell$. Let $\hat{\gamma}_{\ell}$ and
$\hat{\alpha}_{\ell}$ denote the ML and DML estimators trained on the set of
indices not in fold $\ell$. Then select a hyperparameter value by searching
over a grid of potential values, and selecting the value that minimizes a loss
function. Let $\mathcal{L}_{\gamma}(\gamma_{\ell},\mathcal{I}_{\ell};\lambda)$
be the mean squared error of the function $\gamma_{\ell}$ with the
hyperparameter $\lambda$ on the data in $\mathcal{I}_{\ell}$. Let
$\mathcal{L}_{\alpha}(\alpha_{\ell},\mathcal{I}_{\ell};\kappa)$ be the loss
function of the function $\alpha_{\ell}$ with the hyperparameter $\kappa$ on
the data in $\mathcal{I}_{\ell}$; I describe this loss function in Appendix
A.3.

Once the hyperparameter is selected, I evaluate the (debiased) score on the
test sets using the same folds defined above. The full estimation procedure
for the debiased score is below. To use this procedure without the debiasing
correction, omit the bias correction term
$\hat{\alpha}_{\ell}(X_{it};\hat{\kappa})(\ddot{y}_{it}-\hat{\ddot{\gamma}}_{\ell}(X_{it};\hat{\lambda}))$
from each step.

  1. 1.

Select hyperparameter $\hat{\lambda}$ that minimize test-set mean squared
error of the regression:

$\hat{\lambda}=\operatorname*{arg\,min}_{\lambda}\sum_{\ell=1}^{k}\mathcal{L}_{\gamma}(\hat{\gamma}_{\ell},\mathcal{I}_{\ell};\lambda)$

  2. 2.

Select hyperparameter $\hat{\kappa}$ that minimize test-set loss of the double
machine learner:

$\hat{\kappa}=\operatorname*{arg\,min}_{\kappa}\sum_{\ell=1}^{k}\mathcal{L}_{\alpha}(\hat{\alpha}_{\ell},\mathcal{I}_{\ell};\kappa)$

  3. 3.

Evaluate the debiased score using these hyperparameters:

$\hat{\theta}=\frac{1}{N}\sum_{\ell=1}^{k}\sum_{it\in\mathcal{I}_{\ell}}m(\hat{\gamma}_{\ell},X_{it};\hat{\lambda})+\hat{\alpha}_{\ell}(X_{it};\hat{\kappa})(\ddot{y}_{it}-\hat{\ddot{\gamma}}_{\ell}(X_{it};\hat{\lambda}))$

  4. 4.

Find the asymptotic variance of the estimator, after adjusting for within-
panel-unit correlations. Let
$\hat{\theta}_{\ell;it}:=m(\hat{\gamma}_{\ell},X_{it};\hat{\lambda})+\hat{\alpha}_{\ell}(X_{it};\hat{\kappa})(\ddot{y}_{it}-\hat{\ddot{\gamma}}_{\ell}(X_{it};\hat{\lambda}))$
and $\hat{\theta}_{\ell;i}:=1/T\sum_{t=1}^{T}\hat{\theta}_{\ell;it}$. Then the
asymptotic variance is:

$\hat{V}=\frac{1}{N}\sum_{\ell=1}^{k}\sum_{i\in\mathcal{I}_{\ell}}\left\\{\sum_{t=1}^{T}(\hat{\theta}_{\ell;it}-\hat{\theta})^{2}+2\sum_{t=1}^{T-1}\sum_{t^{\prime}=t+1}^{T}(\hat{\theta}_{\ell;it}-\hat{\theta}_{\ell;i})(\hat{\theta}_{\ell;it^{\prime}}-\hat{\theta}_{\ell;i})\right\\}$

In the following subsections, I describe how to train and evaluate the Lasso
and NNet estimators, and introduce the DML procedure.

#### 3.3.a Lasso

Least absolute shrinkage and selection operator (Lasso) is a regression
procedure that selects a sparse linear model from a researcher-specified set
of basis functions. The procedure finds this sparse linear combination by
minimizing squared error of estimation, while penalizing more complex models
via regularization. The procedure is similar to OLS Poly, but estimates can
differ greatly because of this penalization. The hyperparameter $\lambda$ is
the magnitude of the regularization term, which I determine via the cross
fitting procedure outlined above.

As with OLS, I assume that there is a true linear model of the form Equation 4
and take within transformations to remove county-level fixed effects to result
in Equation 5. I use the same flexible set of basis functions used in OLS
Poly, as defined in Section 3.2. Unlike in OLS Poly, I assume that the true
parameter vector is sparse and find an estimate by solving a regularized
optimization problem.

For each cross-validation fold, I find the estimate $\hat{\beta}_{\ell}$ via
the following minimization problem:

$\hat{\beta}_{\ell}=\operatorname*{arg\,min}_{\beta}\sum_{i\in\mathcal{I}_{\ell}}\sum_{t=1}^{T}(\ddot{y}_{it}-\ddot{b}(X_{it})\beta)^{2}+\lambda|\beta|_{1}$
(6)

where $|\beta|_{1}$ is $\ell 1$ penalty or the sum of the absolute value of
each component of $\beta$. To incorporate yearly fixed effects, I include
dummy variables and do not apply the regularization penalty to the
coefficients on those dummy variables. This reflects the assumption that while
coefficients in $\gamma$ are sparse, coefficients in the fixed effects terms
are not [4]. Additional details on the Lasso procedure are included in
Appendix A.1.

Once I have estimated $\hat{\beta}_{\ell}$, I compute the directional
derivative by projecting my estimate of $\beta_{0}$ on partial derivatives of
the basis functions. This is the same procedure to compute the derivative
using OLS, as described in Section 3.2.

#### 3.3.b Neural Network

A neural network (NNet) learns a relationship between inputs and outputs
through an iterative process, where the best fitting model is selected from a
large space of flexible transformations of all potential interactions of input
features. NNets have several advantages for my application. First, it is
straightforward to compute a gradient of the output of the entire network with
respect to each input feature. Second, NNets can incorporate fixed effects in
panel models, as demonstrated by [11]. Third, estimating NNets does not
require the researcher to specify basis functions.

I compute derivatives of the NNet by extracting gradients computed during
training the network. A NNet is a weighted composition of activation functions
(user-specified transformations) applied to an input vector. After a random
initialization of internal parameters, the model goes through many iterations
of predicting the output variable, calculating a loss from training data, and
updating the parameter values. I use the mean squared error for this loss
function. While training or evaluating the network, the algorithm computes the
derivative of the output with respect to each input. This is commonly used to
adjust parameter values during the training procedure, in a process known as
back propagation. I also use these automatic derivatives to compute the
partial derivative of the prediction with respect to each weather input.

I follow [11] to estimate the model with nonparametric treatment of weather
inputs and additive linear fixed effects. Training this network involves an
iterative procedure that can be interpreted as selecting basis functions from
a large space of candidate functions. After this training procedure is
completed, I use standard econometric techniques to account for linear terms,
treating the nonlinear transformation of inputs as a feature in a linear
model. Specifically, I use the within transformation to remove individual
fixed effects and use OLS to find other linear terms such as yearly fixed
effects. The iterative procedure is performed on the set of training data, and
the OLS step is performed on the test set. This requires an architecture with
a top layer that is linear in all fixed effect components and the output of
the nonlinear network transformation. See [11] for more details on such
networks and their performance relative to linear models or fully
nonparametric models in estimating agricultural yields. Additional details of
the NNet architecture and training procedure are included in Appendix A.2.

#### 3.3.c Double Machine Learning

Double machine learning (DML) is an approach to remove bias from standard
machine learning algorithms. I use an approach introduced by [8], relying on
the statistical theory of the Riesz representer.

There are two related, yet orthogonal, problems in estimating the average
derivative. First is the regression function, and second is the derivative
operator on a regression function. [8] show how to conduct this second step
without estimating a regression function, and that these two methods can be
combined to construct an approximately debiased estimate of the original
target derivative.

This second function can be estimated from data using the Riesz representation
theorem. I denote this second function $\alpha_{0}$, and its estimate
$\hat{\alpha}$. The Riesz representation implies that
$\mathbb{E}[\alpha_{0}(X_{i,t})\ddot{\gamma}_{0}(X_{i,t})]=m(\ddot{\gamma}_{0},X_{i,t})$.
[8] show how to use this fact to estimate $\alpha_{0}$ from data. In a panel
setting, I use the within-transformed function $\ddot{\gamma}$; I use
$m(\ddot{\gamma}_{0},X_{i,t})=m(\gamma_{0},X_{i,t})$ because the functional
$m$ is a directional derivative and the derivative of the mean value is 0. I
assume that $\alpha_{0}$ is linear in the set of basis functions $b$. The true
Riesz representer is then $\alpha_{0}(X_{it}):=\ddot{b}(X_{it})\rho_{0}$ and
its estimate is $\hat{\alpha}(X_{it}):=\ddot{b}(X_{it})\hat{\rho}$. I estimate
$\hat{\rho}$ using an optimization package, based on the moment conditions
identified by [8]. [22] demonstrate that this procedure is effective in panel
settings. Appendix A.3 includes more details on the motivation of this Riesz
representer and the estimation procedure.

This estimate is then used to construct the following doubly robust score:

$\hat{\theta}=\frac{1}{N}\sum_{\ell=1}^{k}\sum_{it\in\mathcal{I}_{\ell}}m(\hat{\gamma}_{\ell},X_{it};\hat{\lambda})+(\ddot{y}_{it}-\hat{\ddot{\gamma}}_{\ell}(X_{it}))\hat{\alpha}_{\ell}(X_{it};\hat{\kappa})$

## 4 Simulation Exercise

I conduct a simulation exercise to compare the performance of the estimation
procedure to OLS while varying the set of weather variables used. To match the
setting as closely as possible, I use the empirical distribution of weather
covariates in all trials. I then apply a simulated production function of the
the piecewise-linear functional form from [32].

The simulation exercise focuses on comparing ML, DML, and OLS in a regression
setting with high-dimensional variation. OLS is correctly specified in all
trials. In cases where OLS is not correctly specified, ML and DML would likely
perform better because they are able to represent a richer set of flexible
functional forms. Simulation exercises by [7] demonstrate this property.

(a) Yearly Linear

(b) Yearly Flexible

(c) Monthly Flexible

Figure 3: Violin plots of distribution of parameter estimates from 1,000 Monte
Carlo simulations. The subfigures indicate which set of weather variables are
used, using the notation from Figure 1. Each subfigure has a separate violin
plot for each method used. The plot shows the density of the parameter at the
value on the y axis, based on the simulation trials. The dotted lines within
each plot show the median and upper/lower quartiles of the distribution. The
dotted horizontal line through the plots shows the true value of the
parameter, -0.05.

I conduct 1,000 Monte Carlo simulation trials of this estimation procedure. In
each trial, I first randomly draw a 1,000 county sample, and randomly select a
2-year panel from these counties. I use the weather observations from this
sample, guaranteeing that there are realistic correlations between weather
variables. Then, I generate a $y$ variable according to the following
function:
$y_{it}=a_{i}+\beta_{1}lower_{it}+\beta_{2}higher_{it}+\beta_{3}prec_{it}+\varepsilon_{it}$,
where $lower$ ($higher$) is the total growing season accumulated GDD below
(above) 29°C and $prec$ is the total growing season precipitation. This
functional form closely matches the parsimonious form suggested by [32]. Note
that in this functional form, OLS Linear is correctly specified for all sets
of weather variables I consider. I set
$\beta_{1}=0.02;\beta_{2}=-0.05;\beta_{3}=0.001$, and take $a_{i}\sim N(1,1)$
and $\varepsilon_{it}\sim N(0,1)$.

Table 1: Summary of results from 1000 Monte Carlo simulation draws of the
estimation procedure. The true value of the average derivative in all trials
is -0.05; $\theta$ is the average of the estimated average derivative from all
trials. Each column represents a separate method used for the regression
function. Standard errors are in parentheses, and are computed from the
distribution of bootstrap values.

The results of this simulation exercise are summarized visually in Figure 3
and numerically in Table 1. I estimate $\beta_{2}$, the parameter of interest,
using OLS, ML, and DML for weather variables in the Yearly Linear, Yearly
Flexible, and Monthly Flexible sets (as illustrated in Figure 1). For OLS, I
OLS Linear and OLS Poly as described in Section 3.2. For ML, I consider Lasso
and NNet estimators as described in Section 3.3.a and Section 3.3.b. For DML,
I adjust each ML result with the DML procedure in Section 3.3.c. I use Lasso
DML and NNet DML to describe the Lasso and NNet estimates with the DML
correction. I use the cross-folds and sample splitting procedure described in
Section 3.

With Yearly Linear weather inputs, OLS performs best among all models. OLS
Poly and OLS Linear have the lowest bias, and OLS Linear has the lowest
variance. The machine learning models performed reasonably well, especially
with the DML correction, although they have greater variance and bias than the
OLS results. The central estimate from both DML approaches lie within 7% of
the true elasticity, and a 95% confidence interval contains the true
elasticity.

With Yearly Flexible weather variables, ML estimates have considerably lower
variance than OLS. The estimates using OLS Linear have a standard deviation
approximately three times as large as the DML estimates, while estimates from
OLS Poly have a standard deviation approximately ten times as large. There is
substantial bias from the Lasso estimates, although the DML correction greatly
reduces this bias. The NNet and NNet DML estimates have very little bias;
central estimates from both are within 0.4% of the true value and have less
bias than OLS Linear. The 95% confidence interval from both DML procedures
contains the true elasticity, although the 95% confidence interval from Lasso
without DML does not.

With Monthly Flexible variables, the benefits of using DML are even more
dramatic. OLS Linear and OLS Poly have extremely high variance, and both have
greater bias than the DML methods. The high variance is especially limiting,
as confidence intervals using these approaches are uninformative. As Figure
3(c) shows, the interquartile range of the bootstrap estimates is not visible
on a plot whose range is 8 times the magnitude of the true target elasticity.
OLS Linear estimates a mean value of 2.359, relative to the true value of
-0.05. OLS Poly has a bias of .018, which is greater than bias from Lasso DML
(0.012) or NNet DML (0.0087).

Note that while mean squared error (MSE) can help suggest a preferred model, a
straightforward comparison of MSE does not select the model with lowest bias.
MSE appears higher among ML models than OLS because the MSE reported by ML is
the MSE of the model evaluated on a test set. When MSE from an ML model is
low, this indicates that ML is providing a better fit because the model
extrapolates well to unseen data. When MSE from an OLS approach is low, this
could indicate overfitting. OLS Poly has the lowest MSE in all trials, but is
overfitting the data in the Yearly Flexible and Monthly Flexible cases. In the
data generating process (DGP), the additive error term has variance 1. The MSE
from OLS Poly is much lower than the true squared error from the DGP,
indicating that OLS Poly is fitting the noise instead of the desired pattern
in the data. The MSE from OLS Linear with Monthly Flexible weather terms is
significantly lower than the true squared error, indicating that OLS Linear
may be overfitting the data in that setting. This suggests that a comparison
of MSE alone should not be used to select the preferred model, but the
researcher can consider test-set MSE and performance in simulation trials.

This simulation exercise shows that when a researcher wishes to measure a
relationship with a high-dimensional set of weather variables, DML can provide
results with lower bias and less variance than OLS. With Yearly Flexible or
Monthly Flexible weather, NNet DML has the lowest bias among all models, and
the DML methods have significantly lower variance. OLS performs best in the
Yearly Linear case, which is expected because OLS is correctly specified and
the weather variation is low-dimensional.

## 5 Results

I implement the above estimation procedure to study adaptation to damaging
heat exposure in U.S. corn and soy production. When using the Yearly Linear
set of weather variables (as in [5]), I find little to no evidence of
adaptation. However, when using a richer set of weather variation, I find
evidence that a considerable share of the short-run damages from extreme heat
are offset with greater exposure to those temperatures. To help explain this
discrepancy, I visualize OLS results with Yearly Flexible temperature
variation to show that the simple linear model may not accurately describe the
long-run role of extreme heat.

I run 500 bootstrap trials of the procedure described in Section 3 to measure
the elasticity of crop yields with respect to extreme heat. As described in
Section 2, I study corn and soy production from 1990-2019 in counties east of
the 100 °W meridian. I use a panel of the full data to capture short-run
variation, and capture long-run variation by comparing average yield and
weather from 1990-1999 to 2010-2019. In each bootstrap trial, I draw a random
subsamples of 80% of counties. This resampling scheme addresses the
intertemporal correlation as suggested by [19], while avoiding the risk of
own-sample bias from machine learning methods by including the same
observations in train and test data.

I estimate the elasticity using the three sets of weather variation described
in Figure 1 and all estimation approaches described in Section 3. Each method
uses the within transformation to remove individual fixed effects, and
includes additive yearly fixed effects. For estimation methods that use a set
of basis functions, I use polynomial expansions of all temperature and
precipitation variables, and interactions of the polynomial expansions of each
precipitation variable with the polynomial expansions of each temperature
variable. I use third-order polynomial expansions for Yearly Linear and Yearly
Flexible weather variables, and second-order polynomial expansions for Monthly
Flexible. The machine learning models use a 5-fold cross validation procedure.

I only report results from OLS Linear, Lasso DML, and NNet DML in this
section. The simulation results showed that the debiasing procedures are
effective at reducing bias from naive machine learning methods, and that
estimates using OLS Poly can have extremely high variance. Results from the
machine learning models without the debiasing procedure and from OLS with
polynomial basis functions are included in Appendix B. Generally, estimates
using machine learning without bias correction are similar to their debiased
counterparts, and estimates using OLS Poly have unacceptably high variance
when using richer sets of weather variation.

(a) Corn, Yearly Linear

(b) Corn, Yearly Flexible

(c) Corn, Monthly Flexible

(d) Soy, Yearly Linear

(e) Soy, Yearly Flexible

(f) Soy, Monthly Flexible

Figure 4: Box plots of 500 bootstrap samples of ratio
$1-\hat{\theta}^{LR}/\hat{\theta}^{SR}$, from both corn and soy production.
The subfigures indicate which set of weather variables are used, using the
notation from Figure 1. The number in parentheses at the bottom of each box
plot is the $p$ value for the one-sided test that the ratio is greater than 0.
Each subfigure has a separate box plot for each method used. The line inside
each box is the median, the edges of the box are the upper and lower quartile
of the distribution, and the whiskers extending from each box are 1.5 times
the interquartile range.

Table 2: Summary of results from 500 bootstrap draws of the estimation
procedure. $\theta$ is the mean of the average derivative over bootstrap
trials. Standard errors are in parentheses, and are computed from the
distribution of bootstrap values. N is the average number of samples per
simulation trial. Stars indicate significance at the
$p=0.05(^{*}),0.01(^{**}),\text{ and }0.001(^{***})$ levels, based on a
Z-score from the mean and standard error of bootstrap trials.

Reassuringly, the results using Yearly Linear weather variation are similar to
those of [5]. Table 2 summarizes the results from the bootstrap trials. Long-
run and short-run datasets both find that additional extreme heat decreases
crop yields, with comparable magnitudes for both corn and soy cultivation. The
magnitude of these elasticities are economically significant – an estimate of
-0.005 implies that crop yields decline by 0.5% for each additional day crops
are exposed to temperatures above 29°C. The findings are similar to the
results reported by [5]. They find that the elasticity of corn yield ranges
from $-0.0037$ to $-0.0062$.

As in [5], with Yearly Linear weather variation I find little to no evidence
that damages from the short run are offset in the long run. Figure 4 shows the
results of taking bootstrap samples of the ratio
$1-\hat{\theta}^{LR}/\hat{\theta}^{SR}$. As described in Section 3.1,
$\hat{\theta}^{SR}\ (\hat{\theta}^{LR})$ is the estimated elasticity with
short-run (long-run) variation. For corn cultivation, I fail to reject the
null hypothesis that this ratio is different from 0 for any estimation method.
I test at the $p=0.05$ level, with a Bonferroni correction to account for
taking three hypothesis tests. For soy cultivation, I find mixed results: for
both OLS Linear and NNet DML I reject the null hypothesis at this level, while
with Lasso DML I fail to reject the null hypothesis. For these soy estimates,
the 95% confidence interval is $[0.1278,0.3488]$ using OLS Linear and
$[0.07153,0.3649]$ using NNet DML. Confidence intervals are computed using the
sample mean and standard deviation among bootstrap trials.

Using Yearly Flexible weather variables, I find significant evidence of
adaptation. The panel with short-run variation finds statistically and
economically significant damages from a marginal increase in extreme heat
exposure. Central estimates range from $-0.01057$ to $-0.01248$ for corn, and
$-0.006983$ to $-0.00788$ for soy. However, I fail to reject the null
hypothesis that the long-run elasticity is different from zero, for both crops
using all specifications. Figure 4 confirms this finding. For both crops and
using all estimation methods, I reject the null hypothesis that the degree of
short-run damages that are offset in the long run is equal to zero.

Using this set of weather variables, my preferred estimator is NNet DML.
Simulation results showed the this estimator performed best with Yearly
Flexible weather variables. The MSE also suggests that the NNet Double is
performing well in this case. The MSE of NNet DML decreases for long-run
regressions for both corn and soy. As the MSE of NNet DML is the test-set MSE,
an improvement relative to the MSE from the Yearly Linear weather variables
reflects am improvement in modeling the true functional form instead of
overfitting. Using this preferred estimate, I find a 95% confidence interval
that the ratio is within $[0.3023,1.222]$ for corn and $[0.6331,2.001]$ for
soy.

The results from using the Monthly Flexible set of weather variables are
similar to those using Yearly Flexible weather variables, although with
greater variance. This greater variance is not surprising, as the simulation
results demonstrated that DML and OLS have higher bias and variance using this
set of weather variables. With short-run variation using this set of weather
variables, there is an economically and statistically significant decline in
yields associated with a marginal increase in heat exposure. With long-run
variation, there is no evidence of significant declines in yields from a
marginal increase in heat exposure. Reassuringly, this evidence is consistent
with the result using Yearly Flexible weather variation. Given the greater
variance of these estimates, I fail to reject the null hypothesis that the
ratio of these elasticities is different from zero except using DML methods.
Using OLS Linear, I find that a marginal increase in heat exposure
significantly increases crop yields.

These results show that, when the temperature-crop yield relationship is
modeled flexibly, there are not significant average declines in yields
associated with a marginal change in long-run damaging heat exposure. I do not
estimate a significant decline from long-run exposure to temperatures above
29°C because other variation is better able to explain differences in average
long-run yields. This is not the case when using short-run variation, where
flexible modeling confirms that a marginal increase in exposure to
temperatures above 29°C is significantly associated with a decline in crop
yields. This finding is novel, and suggests significant adaptation to damaging
heat exposure.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 5: Comparison of relationship between temperature and yields with
short-run and long-run variation. Each plot summarizes the coefficient
estimates from results with the Yearly Flexible set of weather variables,
using NNet DML. The bar plot shows the expected change in yields from exposure
to a single degree day at the given temperature level. Error bars show one
standard deviation above/below each estimated coefficient. The line shows the
best piecewise linear fit to the bar plot, with one piece from 0-29°C and
another from 29-40 °C.

To better explain this result, I examine the estimates of temperature
coefficients from each weather bin of the Yearly Flexible weather variables. I
use my NNet DML procedure to estimate these plots. Appendix Appendix C
includes results using alternate estimation procedures. Due to the
computational cost of estimating each average derivative using the DML
procedure, I do not replicate the DML procedure for each of these temperature
coefficients; the results are similar. Figure 5 visualizes the results of
these regressions. Figure 5(b) and Figure 5(d) show the results of this
estimate using short-run variation. These results are similar to the results
from [32], who noted that the piecewise linear fit (the dotted line in these
figures) explains most of the variation from the flexible modeling approach
(the bar chart in these figures). Here, the piecewise linear model seems to be
a good fit: before around 29°C, additional heat exposure is generally
associated with an increase in yields. After this point, additional heat
exposure is associated with a decline in yields.

The results using long-run variation do not share this pattern. Figure 5(a)
and Figure 5(c) show the results of this estimate using long-run variation.
The coefficients do not demonstrate a clear piecewise linear pattern. This
suggests that the piecewise linear model relating temperature exposure and
crop yield may not be appropriate in long-run models. In Appendix C, I
replicate this figure for different time periods, starting from 1950-1979.
These figures confirm that the piecewise linear model is a better fit when
using short-run variation than when using long run variation. While the
results from the Yearly Linear analysis show that long-run damaging heat
exposure is correlated with a decline in crop yields, these results show that
other temperature variation is better able to explain the long-run changes in
crop yields. This shows the importance of flexible modeling to understanding
the role of a marginal increase in damaging heat exposure.

## 6 Discussion

I introduced a DML procedure that can estimate average directional derivatives
in high-dimensional settings, and demonstrated the benefits of this procedure
over OLS estimation in a simulation exercise. The results show that this
procedure can estimate average directional derivatives with lower bias and
variance than OLS. The simulation results also show that my method is less
reliable than having access to a true, parsimonious model of the underlying
function. In settings where the researcher has high-dimensional weather
variation and does not have a strong prior about the true functional form,
this DML approach can be used to estimate elasticities and the degree of
adaptation to a weather feature.

Applying this estimator to panel of U.S. corn and soy yields with a rich set
of temperature variation, I conclude that there has been significant
adaptation to damaging heat exposure. Using this flexible method, a panel of
short-run damages finds evidence that a marginal increase in heat exposure
above 29°C is damaging for both corn and soy yields. However, by constructing
a dataset of average changes to capture long-run variation, I cannot reject
the null hypothesis that a marginal increase in long-run damaging heat
exposure is unrelated to long-run crop yields. OLS estimates using a Yearly
Flexible set of weather variables support this finding. I demonstrate that
long-run exposure to temperatures above 29°C is not clearly associated with
yield declines, as is the case with short-run exposures. This implies that
there has been significant adaptation to climate change in this setting, as
the short-run damages are significantly offset in the long run.

This approach does not offer insight into what form this adaptation may take.
Farmers have many possible adaptation mechanisms, such as investing in
irrigation, purchasing improved seeds, or adjusting planting times. It is
important to learn which mechanisms may be responsible for offsetting damages
in the long run, in order to study the long-term effectiveness of these
mechanisms or the ability to use them in other agricultural contexts. More
research is needed to understand how these damages are offset.

The apparent contradiction between this result and prior literature suggests
that while there is adaptation to a marginal increase in damaging heat
exposure, there is limited adaptation to some other damaging feature of
climate change. Prior literature, such as [5], did not find evidence of
substantial adaptation to damaging heat exposure. My analysis shows that while
a long-run increase in damaging heat exposure is correlated with a decline in
crop yields, variation in heat exposure at other temperature levels is better
able to explain this pattern. This is evidence of omitted variable bias in
long-run estimates using only beneficial and damaging heat exposure. While
prior results indicate limited adaptation to some change in the temperature
distribution, my detailed analysis finds that there has been substantial
adaptation to damaging heat exposure.

## References

  * [1] MOSEK ApS  “MOSEK Optimizer API for Python 9.3.6”, 2021  URL: https://docs.mosek.com/latest/pythonapi/index.html
  * [2] Alan Barreca, Karen Clay, Olivier Deschenes, Michael Greenstone and Joseph S. Shapiro  “Adapting to climate change: The remarkable decline in the US temperature-mortality relationship over the Twentieth Century”  In _Journal of Political Economy_ 124.1 University of Chicago PressChicago, IL, 2016, pp. 105–159  DOI: 10.1086/684582
  * [3] Jonathan T. Barron  “Continuously Differentiable Exponential Linear Units”  In _arXiv preprint_ , 2017  arXiv: http://arxiv.org/abs/1704.07483
  * [4] Alexandre Belloni, Victor Chernozhukov, Christian Hansen and Damian Kozbur  “Inference in high-dimensional panel models with an application to gun control”  In _Journal of Business & Economic Statistics_ 34.4 Taylor & Francis, 2016, pp. 590–605
  * [5] Marshall Burke and Kyle Emerick  “Adaptation to climate change: Evidence from US agriculture”  In _American Economic Journal: Economic Policy_ 8.3 American Economic Association, 2016, pp. 106–140  DOI: 10.1257/pol.20130025
  * [6] Fiona Burlig, Christopher Knittel, David Rapson, Mar Reguant and Catherine Wolfram  “Machine learning from schools about energy efficiency”  In _Journal of the Association of Environmental and Resource Economists_ 7.6 The University of Chicago Press Chicago, IL, 2020, pp. 1181–1217
  * [7] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney K Newey and James Robins  “Double/debiased machine learning for treatment and structural parameters”  In _Econometrics Journal_ 21.1 Oxford Academic, 2018, pp. C1–C68  DOI: 10.1111/ectj.12097
  * [8] Victor Chernozhukov, Whitney K Newey and Rahul Singh  “Automatic debiased machine learning of causal and structural effects”  In _Econometrica_ 90.3 Wiley Online Library, 2022, pp. 967–1027
  * [9] Victor Chernozhukov, Whitney K Newey and Rahul Singh  “Debiased machine learning of global and local parameters using regularized Riesz representers”  In _The Econometrics Journal_ , 2022
  * [10] Arnaud Costinot, Dave Donaldson and Cory Smith  “Evolving Comparative Advantage and the Impact of Climate Change in Agricultural Markets: Evidence from 1.7 Million Fields around the World”  In _Journal of Political Economy_ 124.1, 2016  URL: http://www.journals.uchicago.edu/t-and-c
  * [11] Andrew Crane-Droesch  “Machine learning methods for crop yield prediction and climate change impact assessment in agriculture”  In _Environ. Res. Lett_ 13, 2018, pp. 114003  DOI: 10.1088/1748-9326/aae159
  * [12] Melissa Dell, Benjamin F. Jones and Benjamin A. Olken  “What Do We Learn from the Weather? The New Climate-Economy Literature”  In _Journal of Economic Literature_ 52.3, 2014, pp. 740–798  DOI: 10.1257/jel.52.3.740
  * [13] Tatyana Deryugina, Garth Heutel, Nolan H. Miller, David Molitor and Julian Reif  “The mortality and medical costs of air pollution: Evidence from changes in wind direction”  In _American Economic Review_ 109.12 American Economic Association, 2019, pp. 4178–4219  DOI: 10.1257/aer.20180279
  * [14] Steven Diamond and Stephen Boyd  “CVXPY: A Python-embedded modeling language for convex optimization”  In _Journal of Machine Learning Research_ 17.83, 2016, pp. 1–5
  * [15] Nick Hagerty  “Adaptation to surface water scarcity in irrigated agriculture”  In _Working paper_ , 2021
  * [16] Richard Hornbeck  “The enduring impact of the American Dust Bowl: Short- and long-run adjustments to environmental catastrophe”  In _American Economic Review_ 102.4, 2012, pp. 1477–1507  DOI: 10.1257/aer.102.4.1477
  * [17] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks”  In _Neural Networks_ 3.5 Pergamon, 1990, pp. 551–560  DOI: 10.1016/0893-6080(90)90005-6
  * [18] Solomon Hsiang  “Climate econometrics”  In _Annual Review of Resource Economics_ 8.1, 2016, pp. 43–75  DOI: 10.1146/annurev-resource-100815-095343
  * [19] George Kapetanios  “A bootstrap procedure for panel data sets with many cross-sectional units”  In _The Econometrics Journal_ 11.2 Oxford University Press Oxford, UK, 2008, pp. 377–395
  * [20] Diederik P. Kingma and Jimmy Lei Ba  “Adam: A method for stochastic optimization”  In _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ International Conference on Learning Representations, ICLR, 2015  DOI: 10.48550/arxiv.1412.6980
  * [21] Jon Kleinberg, Jens Ludwig, Sendhil Mullainathan and Ziad Obermeyer  “Prediction policy problems”  In _American Economic Review_ 105.5, 2015, pp. 491–495  DOI: 10.1257/aer.p20151023
  * [22] Sylvia Klosin and Max Vilgalys  “Estimating Continuous Treatment Effects in Panel Data using Machine Learning with an Agricultural Application”  In _arXiv preprint arXiv:2207.08789_ , 2022
  * [23] Pradeep Kurukulasuriya, Namrata Kala and Robert Mendelsohn  “Adaptation and climate change impacts: a structural Ricardian model of irrigation and farm income in Africa”  In _Climate Change Economics_ 2.02 World Scientific, 2011, pp. 149–174
  * [24] Derek Lemoine  “Estimating the Consequences of Climate Change from Variation in Weather”  In _Manuscript. National Bureau of Economic Research_ , 2018  DOI: 10.3386/w25008
  * [25] Konstantinos G. Liakos, Patrizia Busato, Dimitrios Moshou, Simon Pearson and Dionysis Bochtis  “Machine learning in agriculture: A review”  In _Sensors (Switzerland)_ 18.8 Multidisciplinary Digital Publishing Institute, 2018, pp. 2674  DOI: 10.3390/s18082674
  * [26] Robert Mendelsohn, William D. Nordhaus and Daigee Shaw  “The Impact of Global Warming on Agriculture: A Ricardian Analysis”  In _The American Economic Review_ 84 American Economic Association, 1994, pp. 753–771  DOI: 10.2307/2118029
  * [27] Sendhil Mullainathan and Jann Spiess  “Machine learning: An applied econometric approach”  In _Journal of Economic Perspectives_ 31.2, 2017, pp. 87–106  DOI: 10.1257/jep.31.2.87
  * [28] Alan L Olmstead and Paul W Rhode  “Adapting North American wheat production to climatic challenges, 1839-2009.”  In _Proceedings of the National Academy of Sciences of the United States of America_ 108.2 National Academy of Sciences, 2011, pp. 480–485  DOI: 10.1073/pnas.1008279108
  * [29] Ariel Ortiz-Bobea  “Is Weather Really Additive in Agricultural Production? Implications for Climate Change Impacts”  In _SSRN Electronic Journal_ 1, 2013, pp. 3–41  DOI: 10.2139/ssrn.2374506
  * [30] Ariel Ortiz-Bobea  “The role of nonfarm influences in Ricardian estimates of climate change impacts on US agriculture”  In _American Journal of Agricultural Economics_ 102.3 Wiley Online Library, 2020, pp. 934–959
  * [31] J. Park and I.. Sandberg  “Universal Approximation Using Radial-Basis-Function Networks”  In _Neural Computation_ 3.2 MIT Press - Journals, 1991, pp. 246–257  DOI: 10.1162/neco.1991.3.2.246
  * [32] Wolfram Schlenker and Michael J Roberts  “Nonlinear temperature effects indicate severe damages to U.S. crop yields under climate change.”  In _Proceedings of the National Academy of Sciences of the United States of America_ 106.37 National Academy of Sciences, 2009, pp. 15594–8  DOI: 10.1073/pnas.0906865106
  * [33] Wolfram Schlenker and Michael J. Roberts  “Nonlinear Effects of Weather on Corn Yields*”  In _Review of Agricultural Economics_ 28.3, 2006, pp. 391–398  DOI: 10.1111/j.1467-9353.2006.00304.x
  * [34] Vira Semenova and Victor Chernozhukov  “Debiased machine learning of conditional average treatment effects and other causal functions”  In _Econometrics Journal_ 24.2 Oxford University Press, 2021, pp. 264–289  DOI: 10.1093/ectj/utaa027
  * [35] P R Shukla et al.  “IPCC, 2019: Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems”  Intergovernmental Panel on Climate Change (IPCC), 2019
  * [36] Christian Stetter, Philipp Mennig and Johannes Sauer  “Using Machine Learning to Identify Heterogeneous Impacts of Agri-Environment Schemes in the EU: A Case Study”  In _European Review of Agricultural Economics_ , 2022
  * [37] Richard Sutch  “The Impact of the 1936 Corn Belt Drought on American Farmers’ Adoption of Hybrid Corn”  In _Economics of Research and Innovation in Agriculture_ University of Chicago Press, 2011, pp. 195–223  URL: https://www.nber.org/chapters/c11987
  * [38] Hal R Varian  “Big data: New tricks for econometrics”  In _Journal of Economic Perspectives_ 28.2, 2014, pp. 3–28

##  Appendix Appendix A Machine Learning Estimation Details

This appendix additional details on the ML estimation procedures.

### Appendix A.1 Lasso

I vary the Lasso hyperparameter (the regularization penalty) on a grid of 15
values, evenly distributed (in log space) from $10^{-10}$ to $10^{0}$. The
cross folds validation procedure generally selects a hyperparameter from the
interior of this grid, suggesting that the range is appropriate.

Recall the form of the optimization problem:

$\hat{\beta}_{\ell}=\operatorname*{arg\,min}_{\beta}\sum_{i\in\mathcal{I}_{\ell}}\sum_{t=1}^{T}w_{i}(\ddot{y}_{it}-\ddot{b}(X_{it})\beta)^{2}+\lambda|\beta|_{1}$
(7)

As this is a (weakly) convex optimization problem, I solve for
$\hat{\beta}_{\ell}$ using the optimization package CVXPY [14] with the
optimizer Mosek [1]. In the optimization package implementation, the optimizer
can occasionally fail to find a unique optimal solution. In these cases, I
introduce an $\ell 2$ of $10^{-20}$ to find a unique solution. Using the
optimization package ensures that the Lasso procedure converges, although this
procedure is more memory intensive than coordinate descent and will raise an
error if it fails to converge. I find that the model fails to converge in
roughly 5% of simulation trials.

In the main specification, I include yearly fixed effects terms but do not
apply a regularization to these terms. I accomplish this by selectively
applying the regularization factor to terms in $\beta$ not in those yearly
fixed effects terms. The additional fixed effects are retuned before
evaluating the estimator on the test set; this does not impact the estimate of
the gradient, but does improve the mean squared error of the estimator.

I use the optimization package for two reasons. When the solver does not
return an error, the solution is guaranteed to be optimal, unlike in iterative
methods such as stochastic gradient descent or coordinate descent. Some
packages contain this functionality for iterative methods, but I am not aware
of an implementation in Python.

### Appendix A.2 Neural Network

I use a relatively simple network configuration, and vary the width of the
neural network (NNet) in each simulation trial. Three key researcher degrees
of freedom when using a NNet are the depth (the number of layers), the width
(the number of nodes per hidden layer), and the activation function (the
nonlinear function applied to the outputs of each layer). I use a network with
one hidden layer, the minimal depth for which a NNet can approximate an
arbitrary function [17, 31]. In each trial run, I use cross-folds validation
to select the width of the network. To ensure that the network is
differentiable, I use the Continuously Differentiable Exponential Linear Unit
[3] as the activation function in the network.

I select the hyperparameter from a grid space of 2 to 256, evenly spaced in
log terms. The cross folds procedure generally selects a width from the
interior of this set, suggesting that the range is appropriate. To train the
network, I use batch normalization, the Adam optimizer [20], a learning rate
of $0.01$, and $1000$ epochs of training. I retune the parametric components
(the county and yearly fixed effects terms) before evaluating on the test set.
This does not impact the estimate of the gradient, but does improve the mean
squared error of the estimator.

Simulation results show that this simple configuration performs well, although
it is possible that a more complex network configuration would perform better.
Given the computational cost of exploring a wide range of potential network
configurations, I leave such exploration for future work.

### Appendix A.3 Automatic Double Machine Learning

The theoretical basis for the double machine learning (DML) procedure comes
from the Riesz representation theorem. This theorem states that, for the
linear functional333A functional is a scalar summary of a function. The
functional in this setting is the directional derivative. $m$, there exists a
function $\alpha_{0}$ such that for any function $h$:

$\mathbb{E}[m(h,X_{i})]=\mathbb{E}[h(X_{i})\alpha_{0}(X_{i})]$ (8)

Because this relationship holds regardless of function $h$, one can set known
functions for which the target $m(h,X_{i})$ is known and estimate
$\hat{\alpha}$ from the empirical distribution of $X$. Alternate procedures
involve the researcher solving for the functional form of $\alpha$ given the
functional $m$, and often estimating densities or derivatives of densities. I
use this procedure because it does not require the researcher specifying the
form of $\alpha$, but instead learns $\hat{\alpha}$ from data.

I follow [8] to estimate $\hat{\alpha}$ using Lasso and the basis functions
used to estimate the regression function. The goal is to find the estimator
that minimizes the mean squared error between $\hat{\alpha}$ and $\alpha_{0}$,
where $\alpha_{0}$ is the true Riesz representer:

$\hat{\alpha}=\operatorname*{arg\,min}_{\alpha}\mathbb{E}[(\alpha_{0}(X_{i,t})-\alpha(X_{i,t}))^{2}]$

Using the Lasso functional form, I have
$\hat{\alpha}=\ddot{b}(X_{i,t})\hat{\rho}$. I find $\hat{\rho}$ by solving the
following regularized problem:

$\hat{\rho}=\operatorname*{arg\,min}_{\rho}\mathbb{E}\left[(\alpha_{0}(X_{i,t})-\ddot{b}(X_{i,t})\rho)^{2}\right]+\kappa|\rho|_{1}$

Expanding the polynomial and applying the definition of the Riesz representer,
the expression simplifies:

$\hat{\rho}=\operatorname*{arg\,min}_{\rho}-2\mathbb{E}[m(\ddot{b},X_{i,t})]\rho+\rho^{\prime}\mathbb{E}[\ddot{b}(X_{i,t})^{\prime}\ddot{b}(X_{i,t})]\rho+\kappa|\rho|_{1}$

[8] include an additional iterative procedure to scale each component’s
regularization term by the inverse of its variance; I implement their
suggested procedure without modification.

As this is a (weakly) convex optimization problem, I use the optimization
package Mosek [1] to find the optimal value of parameter vector $\hat{\rho}$
for a given value of regularization $\kappa$. [22] present simulation results
comparing this optimization method to the iterative approach from [8]; they
find that the optimization approach results in lower mean squared error in
estimating the true average derivative.

I select the value of the hyperparameter $\kappa$ by minimizing the above loss
function through the cross-folds procedure described in Section 3.3. That is:

$\mathcal{L}_{\alpha}(\hat{\alpha}_{\ell},X_{i,t};\kappa)=-2m(\ddot{b},X_{i,t})\hat{\rho}_{\ell}+\hat{\rho}_{\ell}^{\prime}\ddot{b}(X_{i,t})^{\prime}\ddot{b}(X_{i,t})\hat{\rho}_{\ell}$

I search the hyperparameter over a grid of values suggested by [8]:
$\kappa=c(N-N_{\ell})^{-1/2}\Phi^{-1}(1-.05/p)$, for
$c\in\\{5/4,1,3/4,5/8,9/16,1/2\\}$, where $\Phi^{-1}$ is the inverse of the
standard normal density function and $N_{\ell}$ is the size of the set of
observations in fold $\ell$.

##  Appendix Appendix B Additional Results

This appendix includes additional results from the bootstrap trials described
in Section 5. Table 3 shows the full results from the estimation of corn
yields, and Table 4 shows the full results from the estimation of soy yields.
As expected, the OLS Poly results have higher variance than other methods,
particularly in higher-dimensional settings. The results from machine learning
approaches without bias correction have very low variance, but are likely
biased relative to the debiased estimates. After applying the bias correction,
the variance typically increases. Figure 6 show the bootstrapped values of the
ratio $1-\hat{\theta}^{LR}/\hat{\theta}^{SR}$ from all models for corn yields,
and Figure 7 show the bootstrapped values of this ratio from soy yields. Note
that the variance from the OLS Poly estimates dominates. These standard errors
are large and in charge, particularly in Figure 6(c) and Figure 7(c).

Table 3: Summary of results from 500 bootstrap draws of the estimation
procedure, for corn production. $\theta$ is the mean of the average derivative
over bootstrap trials. Standard errors are in parentheses, and are computed
from the distribution of bootstrap values. N is the average number of samples
per simulation trial. Stars indicate significance at the
$p=0.05(^{*}),0.01(^{**}),\text{ and }0.001(^{***})$ levels, based on a
Z-score from the mean and standard error of bootstrap trials.

Table 4: Summary of results from 500 bootstrap draws of the estimation
procedure, for soy production. $\theta$ is the mean of the average derivative
over bootstrap trials. Standard errors are in parentheses, and are computed
from the distribution of bootstrap values. N is the average number of samples
per simulation trial. Stars indicate significance at the
$p=0.05(^{*}),0.01(^{**}),\text{ and }0.001(^{***})$ levels, based on a
Z-score from the mean and standard error of bootstrap trials.

(a) Corn, Yearly Linear

(b) Corn, Yearly Flexible

(c) Corn, Monthly Flexible

Figure 6: Box plots of 500 bootstrap samples of ratio
$1-\hat{\theta}^{LR}/\hat{\theta}^{SR}$, from corn production. The subfigures
indicate which set of weather variables are used, using the notation from
Figure 1. The number in parentheses at the bottom is the $p$ value for the
one-sided test that the ratio is greater than 0. Each subfigure has a separate
box plot for each method used. The line inside each box is the median, the
edges of the box are the upper and lower quartile of the distribution, and the
whiskers extending from each box are 1.5 times the interquartile range.

(a) Soy, Yearly Linear

(b) Soy, Yearly Flexible

(c) Soy, Monthly Flexible

Figure 7: Box plots of 500 bootstrap samples of ratio
$1-\hat{\theta}^{LR}/\hat{\theta}^{SR}$, from soy production. The subfigures
indicate which set of weather variables are used, using the notation from
Figure 1. The number in parentheses at the bottom is the $p$ value for the
one-sided test that the ratio is greater than 0. Each subfigure has a separate
box plot for each method used. The line inside each box is the median, the
edges of the box are the upper and lower quartile of the distribution, and the
whiskers extending from each box are 1.5 times the interquartile range.

##  Appendix Appendix C Short-Run and Long-Run Estimates of Yearly
Coefficients

I provide additional results to show the difference between short-run and
long-run estimates of yearly coefficient bins. These results confirm that the
functional form of long-run damages does not follow the well-established
pattern of short-run damages. I demonstrate that this conclusions is robust to
different modeling choices and time periods.

### Appendix C.1 Comparisons using alternate estimators

The following plots show results from the comparison between short-run and
long-run estimates for each yearly coefficient bin. All plots use the sample
from the main estimation, corn and soy yields from 1990-2019. Each plot
summarizes the coefficient estimates from results with the Yearly Flexible set
of weather variables, using either Lasso DML, OLS Poly, or OLS Linear. The bar
plot shows the expected change in yields from exposure to a single degree day
at the given temperature level. Error bars show one standard deviation
above/below each estimated coefficient. The line shows the best piecewise
linear fit to the bar plot, with one piece from 0-29°C and another from 29-40
°C.

(a) Corn, Long-Run Variation, Lasso DML

(b) Corn, Short-Run Variation, Lasso DML

(c) Soy, Long-Run Variation, Lasso DML

(d) Soy, Short-Run Variation, Lasso DML

Figure 8: Comparison of relationship between temperature and yields with
short-run and long-run variation, using Lasso DML.

(a) Corn, Long-Run Variation, OLS Poly

(b) Corn, Short-Run Variation, Lasso DML

(c) Soy, Long-Run Variation, OLS Poly

(d) Soy, Short-Run Variation, OLS Poly

Figure 9: Comparison of relationship between temperature and yields with
short-run and long-run variation, using OLS Poly.

(a) Corn, Long-Run Variation, OLS Linear

(b) Corn, Short-Run Variation, Lasso DML

(c) Soy, Long-Run Variation, OLS Linear

(d) Soy, Short-Run Variation, OLS Linear

Figure 10: Comparison of relationship between temperature and yields with
short-run and long-run variation, using OLS Linear.

### Appendix C.2 Comparisons Over Time

In this section, I provide OLS estimates of Yearly Flexible weather variables
using short-run and long-run variation from different historical periods. In
each figure, I compare short-run estimates using a panel dataset all years in
the 30-year range, and long-run estimates comparing average weather and crop
yield values from the first and last decades in that range.

Each figure contains four separate plots. Each plot summarizes the coefficient
estimates from OLS results with the Yearly Flexible set of weather variables.
The bar plot shows the expected change in yields from exposure to a single
degree day at the given temperature level. Error bars show one standard
deviation above/below each estimated coefficient. The line shows the best
piecewise linear fit to the bar plot, with one piece from 0-29°C and another
from 29-40 °C. Regression is taken after the within transformation, and also
includes yearly fixed effects terms and precipitation; these terms are omitted
from the figures.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 11: Comparison of relationship between temperature and yields with
short-run and long-run variation, from 1950-1979.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 12: Comparison of relationship between temperature and yields with
short-run and long-run variation, from 1960-1989.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 13: Comparison of relationship between temperature and yields with
short-run and long-run variation, from 1970-1999.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 14: Comparison of relationship between temperature and yields with
short-run and long-run variation, from 1980-2009.

(a) Corn, Long-Run Variation

(b) Corn, Short-Run Variation

(c) Soy, Long-Run Variation

(d) Soy, Short-Run Variation

Figure 15: Comparison of relationship between temperature and yields with
short-run and long-run variation, from 1990-2019.