\section{Introduction} \label{sec::introduction}
	
	
	Estimating continuous treatment effects in panel data is essential for many important questions. For example, environmental  researchers study how health or economic outcomes change with exposure to high temperatures, air pollution, or other continuous environmental factors.  In economics, applied work commonly studies the impact of a continuous variable such as price (e.g. for demand models in industrial organization) or distance (e.g. for evaluating programs in urban economics) in panel settings.  
	Panel data lets researchers control for time-constant factors that are not recorded in the data, by using observations of the same unit over several time periods. In the widely used fixed effects approach, researchers model time-constant factors as indicator variables for each unit in the data (\cite{wooldridge2005violating}). 
	Most applications using the fixed effects approach for panel analysis rely on linear models and do not attempt to flexibly model high dimensional covariates and their interactions.
	This can can cause significant biases if the linear model is not correctly specified, or if there is heterogeneity in effects. 
	
	Machine learning (ML) approaches like Lasso provide tools to estimate flexible functions of high-dimensional data, but standard models can introduce bias and do not incorporate fixed effects terms. Through regularization and overfitting, ML can result in biased estimates and invalid confidence intervals (\cite{chernozhukov2017double}). 
	There are a growing number of double/debiased machine learning (DML) approaches that address these biases and can construct valid confidence intervals (\cite{chernozhukov2018biased,chernozhukov2022automatic,rothenhausler2019incremental}).
	Current DML approaches for continuous treatments do not account for common features of panel data, such as the fixed effects term. 
	
	This paper fills this gap by introducing a DML estimator for continuous treatment effects in panel settings with additive fixed effects. 
	We present a new panel data estimator of average derivatives that allows for high dimensionality and heterogeneity. After introducing the estimator, we show that it performs well in simulation trials and prove that it is asymptotically normal.
	In simulations, we show that our estimator maintains low bias and valid confidence intervals, while standard Lasso does not. 
	We then use the estimator to study an important question in agricultural economics, the impact of extreme heat exposure on crop yields. Our work relates to the panel, DML, continuous treatment effect, and climate economics literature. First, we discuss the applied and theoretical motivation for the estimator. We then discuss our application to measure damages to U.S. agriculture from extreme heat exposure. 
	We show that in both cases, our DML approach allows greater model flexibility while preserving low standard errors.
	
	
	Flexible continuous effects in panel data are under-explored in the theoretical literature relative to binary treatment effects. While there are many interesting flexible panel methods for binary treatment effects, including factor models and synthetic controls, there are no such methods for continuous effects. The paper studying binary treatment effects in panel flexibly that is most related to ours is \citet{belloni2016inference}. This paper uses a fixed-effects approach, as we do, and allows for the outcome model to be a general function of covariates estimated with ML, but does not allow for general treatment effect heterogeneity as we do. \citet{chernozhukov2022automatic} and \citet{chernozhukov2017orthogonal} allow for general treatment effect heterogeneity and continuous treatment variables. However, \citet{chernozhukov2017orthogonal} considers the dynamic panel and imposes the potentially strong assumption of sequential exogeneity, and \citet{chernozhukov2022automatic} does not consider panels with few observations per unit. Both these papers use a random correlated effects approach, which imposes a restriction on the relationship between the covariates and the unobserved fixed effects. Our fixed-effects approach does not impose that restriction. 
	
	This paper focuses on the \emph{average derivative}, although there are many other potentially interesting continuous treatment effects. For example,
	\citet{colangelo2020double} and \citet{klosin2021automatic} study the dose-response curve, and \citet{chernozhukov2022automatic} considers average policy effects, average treatment effects, and the average equivalent variation bound. 
	The methods in our paper could be adapted to other objects, such as average policy effects or conditional average derivatives. We focus on the average derivative because it can be estimated at parametric root $n$ rates and because there are many applications of the average derivative in practice.
	Many papers estimate an elasticity, which is equivalent to estimating an average derivative when the outcome variable is the logarithm of a raw value. 
	The average derivative is also studied by, among others, \citet{rothenhausler2019incremental} and \citet{imbens2009identification}.
	
	
	We also contribute to the DML literature more generally by introducing a new optimization procedure to estimate the de-biasing component. We show that the optimization procedure improves performance in Monte Carlo simulations. Furthermore, though our paper is the first to use DML to estimate average derivatives with a fixed effects panel data structure, we contribute more generally to the literature on average derivatives by introducing a new analytical method for estimating average derivatives that is very stable in simulations.
	
	We apply our estimator to study the impact of extreme heat on United States (U.S.) corn yields. 
	We estimate the elasticity of corn yield with respect to a marginal increase in extreme heat after controlling for a set of weather variation compiled by \citet{abatzoglou2013development}. 
	Our estimator is well suited to this setting because of the importance of fixed effects modeling to isolate the impact of weather shocks (\cite{Deschenes2007}) and the nonlinear relationships between continuous environmental factors such as temperature and precipitation (\cite{Schlenker2009}).
	
	This application is most closely related to \citet{Schlenker2009}, which estimates this same elasticity. 
	They introduce a parsimonious model that captures the response of crop yields to temperature and precipitation. Our approach differs because we consider flexible functions and interactions between these terms. 
	Using our approach, we project economic damages from extreme heat by the year 2050 that are $\$1.18$ billion more than projections using the weather variation and linear model from \citet{Schlenker2009}. 
	We estimate the value of this elasticity over time, and cannot reject that the elasticity remains constant over time. 
	This corroborates findings from \citet{Schlenker2009} and \citet{Burke2016} who conclude that there has been limited adaptation to extreme heat in U.S. agriculture. 
	
	The application is also related to a literature in environmental economics on studying the economic impacts of climate change. 
	\citet{rode2021estimating} estimates elasticities of energy use with respect to temperature, to construct a global estimate of the social cost of carbon from energy use.
	\citet{Burke2015} and \citet{Hsiang2017} estimate economic damages of climate change globally and within the U.S., respectively, by combining estimates of elasticities of multiple economic sectors with respect to temperature.
	These approaches model the impacts of temperature by separating heat exposure into bins and estimating a coefficient for each bin, a model advocated by \citet{Dell2014} and \citet{hsiang2016climate}. 
	These approaches can have limitations for high-dimensional weather variation, as the number of variables grows exponentially if the researcher wishes to include interactions between weather factors. Our project is the first, to our knowledge, to estimate these elasticities using a machine learning approach, allowing us to flexibly model higher-dimensional variation. 
	
	
	Our project is also related to a growing field applying machine learning techniques in environmental economics. 
	Many researchers have used machine learning for predictive properties, such as forecasting crop yields (\cite{Crane-Droesch2018}), filling missing data to track global poverty (\cite{Jean2016}), or more effectively assign treatment (\cite{knittel2019using}).
	Recent work has begun using machine learners for measurement tasks. 
	\citet{Deryugina2019} uses a machine learning approach to measure the costs of air pollution.
	\citet{stetter2022using} use a de-biased machine learning approach to measure effectiveness of an agricultural intervention. 
	Our paper is the first application in environmental economics, to our knowledge, to use  double machine learning to estimate continuous treatment effects. 
	
	
	The paper is structured in the following way. Section \ref{sec::estimation} sets up the framework of the paper,  introduces the parameter of interest, and presents our estimator. Simulation design and results are given in Section \ref{sec::simulations}. Section \ref{sec::application} covers our application. Section \ref{sec::conclusion} concludes. 
	
	
	
	
	
	\section{Estimation} \label{sec::estimation}
	
	
	\subsection{Notation and Definitions} \label{sec::notation}
	
	We work in a panel data setting with $n$ individuals and $T$ time periods. 
	As is often the case in economic data, we assume that $n$ is large but $T$ is small.
	We assume we have independent and identically distributed data $(W_1, \cdots, W_n)$ where the $W_i = \{(X_i, D_i, Y_i)\}_{t=1}^T$ are copies of a random variable $W$ with support $\{ \mathcal{W} = \mathcal{X} \times \mathcal{D} \times \mathcal{Y}\}_{t=1}^T$, with a cumulative distribution function (cdf) $F_{YDX}(Y,D,X)$. We use capital letters to denote random variables and lowercase letters to denote their possible values. For each unit in a population $X_{i,t} \in {\rm I}\kern-0.18em{\rm R}^{h}$ denotes a vector of covariates, with $h$ potentially large, and $D_{i,t} \in {\rm I}\kern-0.18em{\rm R}$ denotes the treatment variable. 
	
	For a given variable $X$, we use the notation $\Delta {X}_{i,t} :=  X_{i,t} - X_{i,t-1}$ for the first difference transformation. For the first difference transformation of a function $f$ of a variable $X$, we apply the function $f$ before taking the difference: $\Delta f(X_{i,t}) :=  f(X_{i,t}) - f(X_{i,t-1})$. 
	
	Define $|\cdot|_1$ as the $\ell_1$ norm; that is, $|\beta|_1 = \sum_{j=1}^p|\beta_j|$ where $\beta_j$ is the $j^\text{th}$ component of $\beta$ and $p$ is the length of $\beta$.

	
	
	\subsection{Parameter of Interest} \label{sec::param of interest}
	
	We are estimating a general additive fixed effects panel model. 
	\begin{equation}
		\label{eqn:basic_regression}
		Y_{i,t} = a_i + \gamma_0(D_{i,t}, X_{i,t}) + \epsilon_{i,t} \quad \quad E[\epsilon_{i,t}| a_i, X_{i,1}, \cdots, X_{i,T}, D_{i,1}, \cdots, D_{i,T} ] = 0
	\end{equation}
	
	Here $a_i$ represents individual fixed effects, and $\gamma$ is a flexible high dimensional function of treatment, covariates, interactions, and higher order terms. 
	We assume that $\gamma$ is constant throughout time and that it can be estimated well with Lasso. 
	This assumption implies a form of sparsity on covariates. However we do not assume that the fixed effects $a_i$ are sparse.
	
	We apply Lasso after applying a set of basis functions to transform $\{D_{i,t}, X_{i,t}\}$ into a high-dimensional set of covariates. 
	We define a $p \times 1$ dictionary of basis functions $b$ that transforms our original vector of covariates, so that $b(D_{i,t}, X_{i,t}) \in {\rm I}\kern-0.18em{\rm R}^p$. 
	Basis functions can include any desired transformations of the covariates, such as polynomial terms or interactions between variables.
	The assumption that Lasso estimates $\gamma$ implies that there exists a sparse parameter vector $\beta_0 \in {\rm I}\kern-0.18em{\rm R}^p$ such that $\gamma_0(D_{i,t}, X_{i,t}) = b(D_{i,t}, X_{i,t})'\beta_0$. 
	
	We choose these modeling assumptions because they match those commonly used in applied work, while relaxing functional form assumptions on $\gamma_0$. 
	We assume the additive error term $\epsilon_{i,t}$ is mean zero conditional on the history of covariates, an assumption called strict exogeneity that is frequently used in applied work (\cite{wooldridge2010econometric}).  
	
	Our estimation target is the average of a continuous treatment effect, the moment function $m$ given in \eqref{eq:main_target}. In our leading application, we consider the average derivative:
	\begin{equation}
		\begin{aligned}
			\label{eq:main_target}
			\tau_0 &=  {\rm I}\kern-0.18em{\rm E}[ m(W_{i,t}, \gamma_0) ] = {\rm I}\kern-0.18em{\rm E}\left[  \frac{\partial \gamma_0(D_{i,t}, X_{i,t})}{\partial D_{i,t}} \right] 
		\end{aligned}
	\end{equation}

	Note here that the expectation is over both $D_{i,t}$ and $X_{i,t}$, and that the target derivative is only a function of $D_{i,t}$ and not its lagged value. The causal interpretation of \ref{eq:main_target} is the average causal effect of a marginal increase in
	treatment. The average derivative has been studied by, among others, \citet{imbens2009identification} and \citet{rothenhausler2019incremental}.
	When $y_{i,t}$ is in log scale, this parameter captures the elasticity of $y$ with respect to a marginal change in $D$. 
	

	To account for the fixed effect term $a_i$, we introduce a first-differenced version of \ref{eqn:basic_regression}.
	In short panels, it is not possible to consistently estimate $a_i$. By taking a first difference, we can remove the time-invariant factor $a_i$ and consistently estimate $\gamma_0$ (\cite{wooldridge2010econometric}). 
	\begin{equation}
		\label{eq:first_difference}
		Y_{i,t} - Y_{i,t-1} =  \gamma_0(D_{i,t}, X_{i,t}) - \gamma_0(D_{i,t-1}, X_{i,t-1}) + \epsilon_{i,t} - \epsilon_{i,t-1}
	\end{equation}
	\begin{equation}
		\label{eq:delta_first_difference}
		\Delta Y_{i,t}  =  \Delta \gamma_0(D_{i,t}, X_{i,t}) + \Delta \epsilon_{i,t} 
	\end{equation}
	
	Because we apply the first difference transformation after taking the function $\gamma_0$ of our data, our linear representations of $\gamma_0$ and $\Delta \gamma_0$ share the same parameter vector $\beta_0$. That is, if $\gamma_0(D_{i,t}, X_{i,t}) := b(D_{i,t}, X_{i,t})'\beta_0$, then $\Delta \gamma_0(D_{i,t}, X_{i,t}) := \Delta b(D_{i,t}, X_{i,t})'\beta_0$.
	
	Note that we can express the estimation target in terms of the average derivative of $\Delta \gamma_0$:

	
	\begin{equation}
		\label{eqn:average_derivative}
		\tau_0 = {\rm I}\kern-0.18em{\rm E}\left[\frac{\partial \gamma_0(D_{i,t},X_{i,t})}{\partial D_{i,t}}\right] = {\rm I}\kern-0.18em{\rm E}\left[\frac{\partial \Delta \gamma_0(D_{i,t}, X_{i,t})}{\partial D_{i,t}}\right]
	\end{equation}
	
	
	
	
	\subsection{Estimator}
	We construct a de-biased estimator of the average derivative using two Lasso estimated high dimensional functions,  $\Delta \hat{\gamma}$ and $\hat{\alpha}$. The first, $\Delta \hat{\gamma}$, is an estimate of the regression function $\Delta \hat{\gamma}$ in equation \ref{eq:delta_first_difference}, how we estimate this function is explained in detail in Section \ref{sec:lasso details}. The second, $\hat{\alpha}$ is a de-biasing term. We describe $\hat{\alpha}$ in more detail and explain how we construct it in Section \ref{sec:Riesz Representer details} after introducing the full estimator. In short one can consider $\Delta \hat{\gamma}$ to be the initial ML estimate of the effect, and the $\hat{\alpha}$ to be an add on that de-biases this initial estimate. 
	
	
	We use a cross-folds procedure to reduce the risk of overfitting. 
	First the researcher chooses the number of splits $L$ ($L = 5$ is commonly used). Then each unit's indices are randomly partitioned into the $L$ equally sized groups. We use $\ell$ to denote these groups, $\ell = 1, \dots, L$. Denote observations in group $\ell$ by $W_{\ell}$. All observations from a single panel unit are placed in the same fold.
	This is important because the data in different folds should be independent, and there is a dependence within the observations of a single unit.
	Our functions $\hat{\gamma}_\ell$ and $\hat{\alpha}_\ell$ are trained using observations not in group $\ell$. 
	The full estimator takes these two functions $\Delta \hat{\gamma}_\ell$ and $\hat{\alpha}_\ell$ from each fold and uses them to construct our estimate $\hat{\tau}$:
	
	
	\begin{equation}
		\label{eq:estimator}
		\begin{aligned}
			\hat{\tau} &= \frac{1}{n (T - 1)} \sum_{\ell=1}^{L} \sum_{i \in \ell} \sum_{t = 2}^{T} \hat{\tau}_{\ell; i,t} \\
			\hat{\tau}_{\ell; i,t} = \frac{\partial \Delta \hat{\gamma}_{\ell}(D_{i,t}, X_{i,t})}{\partial D_{i,t}} 
			&+  \hat{\alpha_{\ell}}(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i, t- 1} )(\Delta Y_{i,t} - \Delta \hat{\gamma}_{\ell}(D_{i,t}, X_{i,t})) 
		\end{aligned}
	\end{equation}
	
	To compute the asymptotic variance of our estimator, it is necessary to account for correlation of our target within panel units. We assume that errors have a constant correlation within a panel unit but are uncorrelated between panel units.
	Let  $ \hat{\tau}_{\ell; i} = 1 / (T - 1) \sum_{t = 2}^T \hat{\tau}_{\ell; i,t}$. 
	Then the asymptotic variance is: 
	
	\begin{equation}
		\hat{V} = \frac{1}{n (T -1)}   \sum_{\ell=1}^{L} \sum_{i \in \ell} \left\{ \sum_{t = 2}^{T} (\hat{\tau}_{\ell; i,t} - \hat{\tau})^2 +   2 \sum_{t = 2}^{T - 1}  \sum_{t' = t + 1}^{T} (\hat{\tau}_{\ell; i,t} - \hat{\tau}_{\ell; i})(\hat{\tau}_{\ell; it'} - \hat{\tau}_{\ell; i}) \right\}
	\end{equation}
	
	
	Assumptions and the proof for asymptotic normality of our estimator are given in Appendix \ref{section:asymptotic_normality}.
	We work with $(T - 1)$ time periods rather than $T$ because we removed one time period by first-differencing the data.
	
	In the following subsections, we describe how we compute the derivative of the regression function and how we estimate the bias correction term $\alpha$. 
	
	\subsubsection{Derivative of Regression Function}
	\label{sec:lasso details}
	There are two general steps to find the derivative of the regression function ,$\frac{\partial \Delta \hat{\gamma}_{\ell}(D_{i,t}, X_{i,t})}{\partial D_{i,t}} $      ,  in \ref{eq:estimator}.
	\begin{enumerate}
		\item Estimate $\Delta \hat{\gamma}_{\ell}$ for each fold.
		
		\begin{enumerate}
			\item Transform the covariates $\{D_{i,t}, X_{i,t}\}$ using a flexible specification. We do so by using polynomial basis functions of terms and interactions, although other approaches like kernel functions or splines could be used as long as the derivatives are bounded. Let $b(D_{i,t}, X_{i,t})$ denote the resulting $p \times 1 $ dictionary of functions \footnote{For example, when we estimate $\Delta \gamma_{\ell}$ we set $b(D_{i,t}, X_{i,t})$ to be a third order polynomial set of the covariate
				variables and interactions between $D$ and each covariate in $X$.}. Then let $\Delta b(D_{i,t}, X_{i,t}) := b(D_{i,t}, X_{i,t}) - b(D_{i,t-1}, X_{i,t -1})$. 
			We set each function in the dictionary $b$ to have mean $0$ and variance $1$; Appendix \ref{sec:normalization} discusses this procedure. 
			
			\item Find a vector of coefficients $\hat{\beta}$  for our dictionary such that $\Delta \hat{\gamma}(D_{i,t},X_{i,t}) := \Delta b(D_{i,t},X_{i,t})' \hat{\beta}$ is a sparse linear approximation of $\Delta \gamma_0(D_{i,t},X_{i,t})$.
			We do so by solving the following Lasso problem:
			\begin{equation}
				\hat{\beta} = \mathop{\mathrm{argmin}}_{\beta} \left\{\frac{1}{n (T - 1)} \sum_{i = 1}^{n} \sum_{t = 1}^{T} (\Delta Y_{i,t} - \Delta b(D_{i,t}, X_{i,t})' \beta)^2 + r_L|\beta|_1 \right\}
			\end{equation}
			This procedure depends on the regularization weight $r_L$, which we determine by finding values that minimize test-set error in a cross-folds procedure. This procedure is described in Appendix \ref{sec:tuning alpha and gamma}.
			
		\end{enumerate}
		\item Calculate the derivative. 
		
		
		We calculate the derivative analytically.
		Many DML approaches use numeric differentiation; we discuss this alternative in Appendix \ref{sec:numerical differentiation}. Our procedure uses the estimate of $\hat{\beta}_{\ell}$ from the previous step to compute the derivative and derivatives of each function in our dictionary of basis functions. 
		
		\begin{enumerate}
			\item Construct the dictionary $b_D$, a $p \times 1$ dictionary of derivatives of each basis function in $b$.
			For each basis function $b^{j}$ for $j = 1, \dots, p$ in our dictionary of basis functions, define its derivative as follows: 
			\begin{equation}
				b_D^{j}(D, X) = \frac{\partial b^{j}(D,X)}{\partial D}
			\end{equation}
			\item Estimate the average derivative as: 
			\begin{equation}
				{\rm I}\kern-0.18em{\rm E}\left[{\frac{\partial \Delta \hat{\gamma}(D_{i,t}, X_{i,t})}{\partial D_{i,t}}}\right] = {\rm I}\kern-0.18em{\rm E}[b_D(D_{i,t},X_{i,t})'\hat{\beta}]
			\end{equation}
		\end{enumerate}
	\end{enumerate}
	
	\begin{example}
		
		Consider a simple setting where $X_{i,t} \in  {\rm I}\kern-0.18em{\rm R}$, and where $\gamma_0(D_{i,t}, X_{i,t}) = D^2_{i,t} X_{i,t}$. 
		Our basis function dictionary is $b(D_{i,t},X_{i,t}) = \{D_{i,t}, X_{i,t}, D^2_{i,t} X_{i,t} \}$.
		In our linear representation, $\beta_0 = \{0, 0, 1\}$. 
		
		In step 1, we obtain an estimate $\hat{\beta}$ using Lasso. 
		In step 2, we first define the derivative of the basis functions. 
		Here, $b_D(D_{i,t},X_{i,t}) = \{1, 0, 2D_{i,t} X_{i,t} \}$. 
		The estimated average derivative is then: $\hat{\beta}_1 + {\rm I}\kern-0.18em{\rm E}[2D_{i,t}X_{i,t}]\hat{\beta}_3$, where $\hat{\beta}_j$ is the $j^\text{th}$ component of $\hat{\beta}$.
	\end{example}
	
	\subsubsection{De-biasing term}
	\label{sec:Riesz Representer details}
	Now the next term of equation \eqref{eq:estimator} is the de-biasing term, $\hat{\alpha_{\ell}}(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i, t- 1} )$.  Our $\alpha$ is based on the methods of \citet{chernozhukov2022automatic}, which use the Riesz Representation theorem. In short this theorem states that for a given linear functional $m$, there exists a function $\alpha$ such that the following holds for any function $b$:
	
	\begin{equation}
		\label{eq:rr_theorem}
		{\rm I}\kern-0.18em{\rm E}[m(W, b)] = {\rm I}\kern-0.18em{\rm E}[\alpha(W)b(W)]
	\end{equation}
	
	
	
	
	In our case using the definition of the moment and our derivative bases, we have $m(W, b) = \Delta b_D(D_{i,t}, X_{i,t})$.
	These basis functions are the same standardized basis functions as in Section \ref{sec:lasso details}. 
	Plugging this into the above yields: 
	
	\begin{equation}
		\label{eq:rr_theorem_plugged_in}
		{\rm I}\kern-0.18em{\rm E}[\Delta b_D(D_{i,t}, X_{i,t}) ] = {\rm I}\kern-0.18em{\rm E}[\alpha(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i,t-1})\Delta b(D_{i,t}, X_{i,t})]
	\end{equation}
	
	Because this equality holds regardless of the function $b$, it is possible to estimate $\alpha$ from data independently of estimating the function $\gamma$.
	In order to construct an estimate $\hat{\alpha}$, we first assume that $\alpha_0$ has a sparse linear form: $\alpha_0(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i,t-1}) = \Delta b(D_{i,t}, X_{i,t})' \rho_0$. 
	Our estimate is then $\hat{\alpha}(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i,t-1}) = \Delta b(D_{i,t}, X_{i,t})' \hat{\rho}$. 
	We follow \citet{chernozhukov2022automatic} and find $\hat{\rho}$ to minimize the squared loss between $\alpha_0$ and $\hat{\alpha}$:
	
	\begin{equation}
		\hat{\rho} = \mathop{\mathrm{argmin}}_\rho \left\{ {\rm I}\kern-0.18em{\rm E}[(\alpha_0(D_{i,t}, X_{i,t}, D_{i,t-1}, X_{i,t-1}) - \Delta b(D_{i,t}, X_{i,t})' \rho)^2] + r_\alpha|\rho|_1 \right\} 
	\end{equation}
	Additional details about this solution are given in Appendix \ref{sec::origin alpha}. 
	\citet{chernozhukov2022automatic} provides an iterative procedure to solve this problem; we also introduce an alternate approach using an optimization package to solve the minimization problem.
	
	Our implementation using an optimization package guarantees that we find an optimal solution to this minimization problem.
	Iterative approaches may fail to converge to the true value of the parameter, with finite training time. Given the convex nature of the minimization problem, we are able to find an exact solution using optimization software.
	We use the Python package CvxPy (\cite{diamond2016cvxpy}) to set up the problem and the convex optimization software Mosek (\cite{mosek}) to solve the problem. 
	In simulation trials, we compare the performance of iterative and optimizer approaches for determining $\hat{\alpha}$. 
	
	With an estimate $\hat{\alpha}$, we can estimate the average derivative using only the $Y_{i,t}$ observations: $\hat{\tau} = {\rm I}\kern-0.18em{\rm E}[\hat{\alpha}(D_{i,t},X_{i,t},D_{i,t-1}, X_{i,t-1})\Delta Y_{i,t}]$.
	
	We now provide an example to illustrate the role of the Riesz Representer function in a simple setting. 
	
	\begin{example}
		Consider the average derivative of a function $\gamma$ of $X$ where $X \sim N(0,1)$. Via integration by parts after expanding the expectation,  $\displaystyle {\rm I}\kern-0.18em{\rm E}\left[ \frac{\partial \gamma(X)}{\partial X}\right] = {\rm I}\kern-0.18em{\rm E}[X \gamma(X)]$.
	
	
		
		That is, the Riesz Representer is $\alpha_0(X) = X$ when $X \sim N(0,1)$ and $m$ is the derivative operator.
		The Riesz Representer depends on the distribution of data and the operator $m$, but does not depend on the choice of function $\gamma$.
		
		If our basis function is the identity (i.e. $b(X_i) = X_i$), then in our linear representation above, $\rho_0 = \{1\}$. 
		We construct an estimate $\hat{\rho}$ using our Lasso procedure, and then estimate the average derivative as: $\hat{\tau} = {\rm I}\kern-0.18em{\rm E}[ X_i \hat{\rho} Y_i]$
	\end{example}
	
	
	
	\section{Simulations} \label{sec::simulations}
	
	
	We demonstrate the performance of our estimator through a Monte Carlo simulation exercise. We use the data generating process (DGP) defined below to create 1000 different datasets. We report the true value of the derivative, our estimate, the average bias of our estimates, and mean squared error (MSE) between true and estimated values. 
	We compare the performance of our estimator with OLS on the untransformed set of covariates (OLS Linear), OLS on the basis function transformation of covariates (OLS Poly), Lasso without a bias correction term, and DML using an iterative estimation procedure (DML Iterative).
	
	
	\paragraph{DGP}: $N = 1000$ number of individuals, $T = 2$ number of time periods, $h = 20$ original number of $X$ covariates. 
	Our basis function transformation takes third order polynomials of each variable, then adds interactions between each $\{D, X^{(j)}\}$ pair and their polynomials. We have a total of $p=244$ covariates after applying the basis function transformation.
	
	We generate outcome variables according to the following function:
	\begin{equation}
		Y_{i,t} = a_i + D_{i,t} + D_{i,t}^2 + D_{i,t}^3 + D_{i,t} X^{(1)}_{i,t}+ .1 \theta \mathbf{X}_{i,t}  + \epsilon_{i,t}
	\end{equation}
	
	To match real-world panel settings, we impose a correlation between $a_i, D_{i,t} $ and $\mathbf{X}_{i,t}$. 
	We set $\theta$ so that the $j^\text{th}$ element is $\theta_j = 1/j^2$.
	Fixed effects $a_i$, covariates $\mathbf{X}_{i,t}$, and random noise $\epsilon_{i,t}$ are drawn from Gaussian distributions: $a_i \sim N(1,1)$, $X^{(j)}_{i,t} \sim N(a_i, 1) \ j = 1,\dots,h$, and $\epsilon_{i,t} \sim N(0,1)$, while the treatment is correlated with $\mathbf{X}_{i,t}$ but includes simulation draws from a Beta distribution:  $D_{i,t} \sim .1 \theta \mathbf{X}_{i,t} + Beta(1, 7) $
	
	We show the results of this simulation trial by summarizing the bias and MSE in Table \ref{tab:simulation results} and plotting the distribution of the error in Figure \ref{fig:ridge plot of sims }. 
	
	\begin{table}[ht]
		\centering
		\input{tables/derivative_estimates_p20t2}
		\caption{Summary of derivative estimates from 1000 bootstrap trials of our simulation procedure. Bias is the average of the estimated value of the derivative minus the true value of the derivative in each simulation draw. ``MSE $\tau$'' is the mean squared error between the true average derivative and the estimated average derivative in each simulation draw, while ``MSE  $\gamma$ in sample'' and ``MSE $\gamma$ cross folds'' refer to the mean squared error of regression from own-sample and out-of-sample estimation. }
		\label{tab:simulation results}
	\end{table}
	
	\begin{figure}[ht]
		\centering
		\includegraphics[width=\textwidth]{images/ridge_plot_p20t2.png}
		\caption{Visualization of distributions of error per simulation trial, for the five methods above. Error is defined as estimated value minus true value. }
		\label{fig:ridge plot of sims }
	\end{figure}
	
	Table \ref{tab:simulation results} demonstrate the value of our approach for applied research. If an applied researcher were faced with our simulated data with 20 covariates and they wanted to use classic OLS methods they could run either a classic linear regression with 20 covaraiates (``OLS Linear'')  or they could run a more flexible linear regression with polynomial terms of these 20 base covaraites leading to 244 covariates (``OLS Poly''). Note that in order to recover the average derivative using OLS Poly, the researcher must use either derivatives of the basis functions, as described in Section \ref{sec:lasso details}, or numerical differentiation. 
	
	Researchers may prefer OLS Linear to OLS Poly because the large number of regressors in OLS Poly can can lead to large standard errors, even though OLS Linear may be biased due to misspecification. The bias and variance trade off of OLS Linear and OLS Poly is exactly what we see in the table.  OLS Linear has smaller standard errors than OLS Poly, which is why applied researchers may prefer it, however OLS Linear has a huge bias relative to OLS Poly. 
	
	Our DML method provides an alternative that preserves the low standard errors of OLS Linear and the low bias of flexible modeling such as OLS Poly. Due to the regularization of Lasso, our estimates have considerably lower standard errors as OLS Poly even though they use the same set of basis functions. 
	Figure \ref{fig:ridge plot of sims } visualizes this benefit - the error from DML and OLS Poly are both approximately centered around 0, but the distribution of errors using OLS Poly is much wider. 
	In this simulation trial, our method has lower bias than OLS Poly and lower standard errors than OLS Linear.
	
	Reassuringly, our DML procedure results in an approximately de-biased estimate of the average derivative. 
	This is especially clear relative to Lasso without any correction. The bias of base Lasso is 100 times larger than the bias of our proposed method.  
	Our estimates also have substantially lower bias than OLS using the untransformed covariates (``OLS Linear'') and lower bias than OLS using the basis function transformation  (``OLS Poly''). 
	We expect these results, as misspecification of OLS Linear can induce a bias in estimating the average derivative, and OLS Poly may be overfitting the data (as shown by the high error in MSE cross folds).
	Our optimization-based DML approach results in lower bias and standard deviation in the result than DML Iterative, although both dramatically improve over Lasso without a correction. This comparison was done by creating our estimator and then using two different algorithms to solve for the de-baising term. 
	
	We use the mean squared error of estimating the true parameter (MSE $\tau$) to compare results, incorporating differences in both  bias and variance. 
	The MSE $\tau$ using our estimation procedure is the lowest among all models considered, roughly half the magnitude from Lasso or OLS Linear and one quarter the magnitude from OLS Poly. 
	The MSE $\tau$ is closest to the result using DML Iterative, and shows that our general optimization solution to de-biasing is leading to lower bias. 
	
	\section{Application} \label{sec::application}
	\label{sec:application}
	In our application, we measure the elasticity of corn yields with respect to extreme heat exposure and use this to project the damages from expected changes in the distribution of extreme heat following climate change.
	The treatment variable $D_{i,t}$ is annual aggregate exposure to temperatures above 29 $^{\circ}$ C; for simplicity we refer to this variable as extreme heat. 
	This variable roughly captures the amount of heat stress a plant experiences.
	As \citet{Schlenker2009} demonstrate, crop yields are generally decreasing in extreme heat, while increasing in heat exposure below 29 $^{\circ}$ C.
	The covariates $X_{i,t}$ include other weather features, such as heat exposure below this temperature threshold, precipitation, humidity, and solar radiation. 
	Heat exposure is measured in Growing Degree Days (GDD), the amount of time a crop is exposed to temperatures between an upper and lower bound during the March-August growing season.
	
	\subsection{Data Description}
	
	We work with two datasets for the main estimation: corn yield from the USDA's Survey of Agriculture, and weather records from \citet{abatzoglou2013development}. 
	To form projections about the impacts under climate change, we use data about projected future impacts across the U.S. that are generously shared by Burke and Emerick. 
	
	For corn yields, we use county-level reports from the USDA Survey of Agriculture, supplemented with additional data generously provided by \citet{Burke2016}.
	We use the metric bushels per acre, and focus on U.S. counties east of the 100$^{\circ}$ West meridian. 
	This region is commonly used in studies of U.S. agriculture because the region west of this meridian relies on heavily subsidized irrigation. 
	In our sample period, over 90\% of corn in the U.S. was grown in this region. 
	This dataset also includes the planted area of corn per county in each year. 
	
	For weather observations, we use the GridMET weather dataset from \citet{abatzoglou2013development}. 
	This dataset contains a rich set of daily weather variables since 1979 at high spatial resolution (4 km) across the United States. In our main results we follow \citet{Schlenker2009} and include only precipitation, beneficial heat exposure, and damaging heat exposure. In Appendix \ref{section:new_weather_vars}, we include additional results with the dataset's 9 raw variables\footnote{GridMET also includes several derived variables, such as Energy Release Component, for specific ecological applications.}: specific humidity, precipitation, minimum relative humidity, maximum relative humidity, surface downwelling shortwave flux in air (a measure of solar radiation), minimum air temperature, maximum air temperature, wind speed, and wind  direction.
	We generate a county-level daily dataset by taking the average of daily weather from all grid observations within each U.S. county. 
	
	To construct a dataset we can use to study the affects of weather on annual crop yields, we aggregate daily weather observations up to the March-August growing season. 
	We take averages of the variables  specific humidity,  minimum relative humidity, maximum relative humidity, surface downwelling shortwave flux in air (a measure of solar radiation), wind speed, and wind  direction. 
	We take the sum of precipitation. We then construct the variables beneficial and damaging heat exposure, total growing season heat exposure below and above (respectively) 29 $^{\circ}$ C. 
	This gives us a county-level dataset with rich weather variation for our empirical application.
	
	
	\subsection{Empirical Results}
	
	In this section we present our primary empirical results. In our primary results the set of covariates we use comes from the current literature (\cite{Schlenker2009,Burke2016}) and use total growing season precipitation and heat exposure below and above 29 $^{\circ}$ C. The specification is our preferred one as it allows for the cleanest comparison between DML and the current literature given that it only uses the variables currently used by applied researchers.  In another set of empirical results in Appendix \ref{section:new_weather_vars} we use a second set of covariates. This second set of covariates includes go beyond the data used in applied work and include all  weather variables at the total growing season level. Though this gives us more variables to work with, but may also make the interpretation of our effect harder to understand. We include results using this second set as they may be of interest and because the difference between estimates becomes smaller. In section we are comparing results using the first common set of weather variables and focus on the changes from DML relative to OLS Linear and not the changes from including additional weather features.  
	
	We estimate the model using the standardization, differencing, cross folds, analytical derivative, and de-biasing procedures described in Section \ref{sec::estimation}.
	We use third order polynomial basis functions.
	We consider interactions up to depth 2 of all variables. 
	Estimation of Lasso for both the regression problem and the Riesz representer uses the optimization package Mosek (\cite{mosek}).
	Table \ref{tab:ag results table} shows the results from the main analysis.
	
	We weight estimates by the area of corn planted when computing coefficients or average derivatives. This choice is common in the applied literature, as the object of interest is the elasticity of yield in an acre of corn with a marginal increase in damaging heat exposure. 
	We follow \citet{Burke2016} and weight each observation by the area of corn planted in that county in that year.
	Appendix \ref{sec:: ag results no weighting} includes results without this weighting choice. Without weighting the difference between OLS and DML estimates become even larger.  
	
	\begin{table}[p]
		\small  
		\centering 
		\caption{Temperature and Precipitation}
		\input{tables/weighted_derivative_estimates_simple.tex}
		\label{tab:simplest table}
		\caption[Elasticity of Corn Yields using DML]{\footnotesize Estimates of elasticity of corn yields with respect to increase in growing season exposure to extreme heat, using two sets of weather covariates.
			Standard errors (in parentheses) are clustered at the county level. 
			See text for estimation details.
		
		
		}
		\label{tab:ag results table}
	\end{table}
	
	
	\begin{figure}[p]
		\centering
	
		\centering
		\includegraphics[scale=.65]{images/weighted_climate_extrapolation_simple.png}
		\caption{Temperature and Precipitation}
		\label{fig:simple climate extrapolation}
		\caption[Climate extrapolation]{ \footnotesize Extrapolating impacts of extreme heat to crop yields by the year 2050, using elasticities from OLS Linear and DML. 
			Each dot represents a central estimate from a model, and the error bar represents the 95\% confidence interval. 
			Dotted line represents the median value across climate models.
			See text for estimation details. 
		}
		\label{fig:climate extrapolation}
	\end{figure}
	
	The average derivative estimate can be interpreted as the elasticity of corn yield with respect to additional exposure to extreme heat. 
	That is, each estimate is the percentage by which yields change with an additional growing degree day of heat exposure above 29 $^{\circ}$ C over the growing season. 
	The magnitude is relatively large - this suggests that an increase of a single growing degree day is associated with corn yields declining by 0.51\%-0.56\%. 
	This finding is line with findings form other analyses of damages to crop yields from extreme heat.\footnote{\cite{Burke2016} estimate this same elasticity in a range from -0.0036 to -0.0062, depending on the specification.}
	
	
	
	
	MSE from both in-sample and out-of-sample prediction is slightly higher for OLS Linear than any of the flexible methods. There is not a large difference between the MSE of the OLS Poly and Lasso models.
	This indicates that, for these sets of weather variables, OLS Poly does not over-fit the data. This is likely because we use a large data-set. 
	
	To illustrate the significance of these estimates, Figure \ref{fig:climate extrapolation} shows the projected overall change in crop yield by 2050 due solely to an increase in damaging heat exposure. 
	We apply the estimated elasticity to climate projections from a range of scenarios. 
	To generate these projections, we find the weighted average degree of expected warming from 2015 to 2050 under various climate scenarios, where we weigh the degree of warming by the average area of corn cultivated per county.
	Climate scenarios are derived from 18 global climate
	models running the A1B emissions scenario. 
	The per-county degree of expected warming under each climate scenario are generously shared by \citet{Burke2016}. 
	We compare the OLS Linear and DML estimates.
	
	These results illustrate that modeling assumptions significantly impact overall damage projections with the simple set of weather variables. 
	There are large sources of uncertainty resulting from the varying climate projections, though the difference becomes smaller when using all weather covariates as seen in Appendix \ref{section:new_weather_vars}. 
	The projected damages are quite significant - in the median emissions scenario, log yields are 0.37-0.41 lower than a world that does not experience climate change. 
	To put this in dollar value, the 2017 Census of Agriculture reported the total value of sales of corn in the United States as \$ 51.2 billion . 
	That range of damage estimates translates to a dollar value of \$16.1
	-\$17.3 billion (in 2017 dollars).
	Our preferred estimate, DML with the widely used simple set of weather covariates, finds an estimated damage of \$17.3 billion.
	
	The difference between OLS Linear and DML estimates is  statistically significant and economically meaningful in our primary empirical results. 
	OLS Poly and DML find similar results, although the standard error from OLS Poly is significantly higher than that of DML.
	The difference between OLS Linear and DML is significant with a p value of 0.00012, while the difference between OLS Linear and OLS Poly is significant with a p value of 0.0055. 
	Median damage estimates using DML instead of OLS correspond to an additional \$1.18 billion. 
	
	These values overstate the degree of damages from climate change for several reasons, but still provide a valuable illustration of the economic significance of using this measurement technique. 
	This is not a complete projection of climate change damages, but only highlights the contribution of the elasticity estimated above using various modeling assumptions.
	Also note that an elasticity is the marginal impact at the observed distribution of weather characteristics; as the distribution of weather patterns changes with climate change, this elasticity will likely change. 
	Additionally, this estimate does not attempt to account for adaptation to climate change. 
	
	We use our estimation procedure to study how the elasticity has changed over time, using 2-period panels from 1980 through 2019. 
	A key advantage of our DML procedure is that it allows de-biased estimation of flexible functions, even in relatively short panels. This allows us to study how the elasticity has changed over time. 
	This research design is similar to \citet{Barreca2016}, who study the mortality response from temperature during separate decades across the 20th century to conclude that there has been significant adaptation to extreme heat exposure. 
	
	\begin{figure}[t]
		\centering
		{
			\captionsetup{justification=centering}
		
			\begin{subfigure}[b]{0.35\textwidth}
				\centering
				\includegraphics[width=\textwidth]{images/weighted_line_plot_OLS_Linear_simple_2.png}
				\caption{OLS Linear \\ Trend 3.43e-05, p-val 0.433}
			\end{subfigure}
			\begin{subfigure}[b]{0.35\textwidth}
				\centering
				\includegraphics[width=\textwidth]{images/weighted_line_plot_DML_simple_2.png}
				\caption{DML  \\ Trend 7.56e-05, p-val 0.538 }
			\end{subfigure}
			
		}
		\caption[Elasticity over time]{ \footnotesize Estimating elasticity of corn yield with respect to extreme heat over time. We use our estimation procedure on 2-period samples from 1980 through 2019. Line shows central estimate, and grey band shows 95\% confidence interval.    }
		\label{fig:elasticities over time}
	\end{figure}
	
	
	Our estimates suggest that the degree of damages from extreme heat has remained approximately stable from 1980 to 2019. 
	Figure \ref{fig:elasticities over time} shows these estimates, as well as the value and statistical significance of the trend in elasticity. 
	The trend and standard errors are computed using weighted least squares, weighted by the inverse of the sum of the variance of the elasticity estimate and the residual from the regression of elasticity on the year.      
	Using both methods with both sets of weather variables, we are unable to reject the null hypothesis at the $p=0.05$ level, suggesting little adaptation to climate change. 
	This finding is in line with the literature on adaptation in U.S. agriculture. \citet{Schlenker2009} and  \citet{Burke2016} also find limited evidence of adaptation to extreme heat. 
	
	
	\section{Conclusion} \label{sec::conclusion}
	
	
	
	In this paper, we presented the first DML based estimator for continuous treatment effects in the classic fixed effect panel setting. It is also the first DML estimator in this setting to allow for general unrestricted heterogeneity in treatment effects. Our estimator is one of the first methods for flexibly estimating continuous treatment effects in high-dimensional panel data, and our estimation approach improves upon the current state-of-the-art methods in simulations. We gave clear guidance on how to implement our methods so they can be of use to an applied audience. 
	
	Our estimator is of broad applied interest - it allows for fixed effects in panel data, unrestricted heterogeneity in treatment effects, and ML tools, and keeps standard errors small even with very flexible models.
	This allows the researcher to flexibly model high-dimensional data while addressing unobservable fixed effects and machine learning bias. 
	
	This is one of the first applications of DML in environmental economics.  We applied the estimator to study how extreme heat impacts crop yields. In our application, our DML estimator leads to estimates of treatment effects that are significantly different from the estimates based on OLS. Using a widely used set of weather variation, the difference is significant at the $p = 0.001$ level and economically meaningful. Our method predicts an effect of warmer temperature that is \$1.18 billion larger than current OLS based estimates.
	
	\newpage
	\printbibliography
	
	\newpage 
	










