\section{Introduction}\label{intro}

By definition, a \emph{matroid} on a finite set $E$ is given by a rank function
$ \rk: 2^E \to \Z_{\geq 0} $
satisfying the following:
\begin{enumerate}[]\itemsep 5pt
\item[] (Submodularity) For any $A_1,A_2 \subseteq E$, we have $\rk(A_1 \cup A_2) + \rk(A_1 \cap A_2) \leq \rk(A_1) + \rk(A_2)$.
\item[] (Monotonicity) For any $A_1 \subseteq A_2 \subseteq E$, we have $\rk(A_1) \leq \rk(A_2)$.
\item[] (Boundedness) For any $A \subseteq E$, we have $\rk(A) \leq |A|$.
\item[] (Normalization) The rank of the empty subset is zero.
\end{enumerate}
Sans ``boundedness'', the axioms above define a \textit{polymatroid}.
Throughout this paper, we assume that the polymatroid is \emph{loopless}:
\begin{enumerate}[]\itemsep 5pt
\item[] (Looplessness) The rank of any nonempty subset is nonzero.
\end{enumerate}
If $P$ is a polymatroid on $E$, then its \emph{rank} is $\rk(P) \coloneqq \rk(E)$.
A \textit{flat} of $P$ is a subset $F \subseteq E$ that is maximal among sets of its rank.
Ordered by inclusion, the flats of $P$ form a lattice $\L_P$.\footnote{Unlike in the case of matroids, the lattice of flats of a polymatroid can fail to be graded or atomic.}
The intersection of two flats is a flat, so any subset $A$ of $E$ is contained in a unique minimal flat $\cl_P(A)$, called the \emph{closure} of $A$ in $P$, which is obtained by intersecting all flats that contain $A$.

Matroids can be viewed as combinatorial abstractions of hyperplane arrangements.
More generally, polymatroids can be viewed as combinatorial abstractions of subspace arrangements.

\begin{example}\label{subspace}
  Let $V_1, \ldots, V_n$ be linear subspaces of a vector space $V$ over a field $\mathbb{F}$.
  The rank function
  \[ \rk(A) \coloneqq \codim_V(\cap_{i \in A} V_i) \]
   defines a polymatroid $P$ on the set of indices $E=\{1,\ldots,n\}$, which is a matroid if and only if every $V_i$ is a hyperplane.
   The map
   $F \mapsto \cap_{i \in F} V_i$
   is a bijection between the flats of $P$ and the subspaces of $V$ obtained by intersecting some of the $V_i$'s.
   A polymatroid arising in this way is said to be \word{realizable} over $\mathbb{F}$, and the subspace arrangement is called a \emph{realization} of $P$ over $\mathbb{F}$.
\end{example}

Much of a hyperplane arrangement's combinatorial data is captured by intersection theory on its wonderful compactification \cite{dCP95}.
To extend this from hyperplane arrangements to non-realizable matroids, one must replace the wonderful compactification with a combinatorial object, the Bergman fan of a matroid.
The purpose of the present paper is to introduce the \emph{Bergman fan of a polymatroid}, a combinatorial model for the wonderful compactification of a subspace arrangement.
As in the case of matroids, the Bergman fan of a polymatroid is a tropical variety of degree one.
In Section \ref{sec:chow}, we show that the Chow ring of a polymatroid satisfies the K\"{a}hler package with respect to any strictly convex piecewise linear function on its Bergman fan,
 recovering and strengthening a result of Pagaria and Pezzoli \cite[Theorems 4.7 and 4.21]{PP21}.

Our construction of the Bergman fan is inspired by a geometric observation: over an infinite field, the wonderful compactification of any subspace arrangement can be realized as the wonderful compactification of a hyperplane arrangement, taken with respect to an appropriate building set (\cref{rmk:geomOfBS}).
The construction immediately reveals that the Bergman fan of a polymatroid and the Bergman fan of the associated matroid have the same support.
Thus, the K\"{a}hler package for the polymatroid follows from that of the associated matroid \cite{AHK18}
and the general fact that the validity of the K\"{a}hler package for the Chow ring of a fan depends only on the support of the fan \cite{ADH20}.

\subsection{The Bergman fan of a Boolean polymatroid}\label{intro:boolean}
An important special case is that of Boolean polymatroids.
Let $\pi:\widetilde E \to E$ be a surjective map between finite sets.
The \textit{Boolean polymatroid} $B(\pi)$ is the polymatroid on $E$ defined by the rank function
\[
\rk_{B(\pi)}(A) = |\pi^{-1}(A)| \ \ \text{for $A \subseteq E$.}
\]
We write $N_{\widetilde E}$ for $\mathbb{Z}^{\widetilde E} / \mathbb{Z}(1,1,\ldots,1)$,
and, for a subset $S$ of $\widetilde E$, write $\e_S$ for the vector
$\sum_{i \in S} \e_i$ in $N_{\widetilde E} \otimes \R$.

\begin{dfn}\label{dfn:booleanBergman}
 The \word{Bergman fan} $\Sigma_{B(\pi)}$ of the Boolean polymatroid $B(\pi)$ is the fan in $N_{\widetilde E}\otimes \R$ with cones
  \[
    \sigma_{\mathscr F, A} \coloneqq \cone( \e_{\pi^{-1}(F_1)}, \ldots, \e_{\pi^{-1}(F_k)}) + \cone(\e_i)_{i \in S},
  \]
  for every chain $\mathscr{F} = \{\emptyset \subsetneq F_1 \subsetneq F_2 \subsetneq \cdots \subsetneq F_k \subsetneq E \}$
  and subset $S$ of $\widetilde E$ not containing any fiber of $\pi$.
\end{dfn}

Throughout the paper, we write $n$ for the cardinality of $E$.

\begin{dfn}
  An \emph{ordered transversal} of $\pi$ is a sequence $s_1,\ldots,s_n$ of elements of $\tilde{E}$ such that each fiber of $\pi$ contains exactly one element of the sequence.
  The \textit{polypermutohedron}  $Q(\pi)$ is the convex hull of the vectors
  $\sum_{i =1}^n i \e_{s_i}$ in $\R^{\widetilde E}$,
    where $s_1, \ldots, s_n$ range over all ordered transversals of $\pi$.
\end{dfn}

In \cref{appendix:polytope}, we show $\Sigma_{B(\pi)}$ is the inner normal fan of the polypermutohedron $Q(\pi)$. It follows immediately from Definition \ref{dfn:booleanBergman} that $\Sigma_{B(\pi)}$ is a complete unimodular fan in $N_{\widetilde E}$.

\begin{example}
When $\pi$ is a bijection, an ordered traversal of $\pi$ is a permutation of $\tilde{E}$, and $Q(\pi)$ is the standard permutohedron in $\R^{\widetilde E}$. This recovers a familiar statement: the Bergman fan of the Boolean matroid is the normal fan of the standard permutohedron.
\end{example}
\begin{example}
When $E$ is a singleton, an ordered transversal of $\pi$ is an element $\tilde{E}$, and $Q(\pi)$ is the standard simplex in $\R^{\widetilde E}$.
Thus, the Bergman fan of a Boolean polymatroid on a singleton is the normal fan of the standard simplex.
\end{example}

\subsection{The Bergman fan of a polymatroid}

Let $P$ be a polymatroid on $E$, and let $\pi:\widetilde E \to E$ be a surjective map satisfying $\rk_P(i)=|\pi^{-1}(i)|$ for every $i$ in $E$.

\begin{dfn}\label{dfn:maximalBergman}
   The \word{Bergman fan} $\Sigma_P$ of the polymatroid $P$ is the subfan of $\Sigma_{B(\pi)}$ with cones
  \[
    \sigma_{\mathscr F, S} \coloneqq \cone( \e_{\pi^{-1}(F_1)}, \ldots, \e_{\pi^{-1}(F_k)}) + \cone(\e_i)_{i \in S},
  \]
  one for every chain of flats $\mathscr F = \{\emptyset = F_0 \subsetneq F_1 \subsetneq F_2 \subsetneq \cdots \subsetneq F_k \subsetneq E \}$ of $P$ and a subset $S$ of $\widetilde E$ such that
\[
\rk(F \cup \pi(T)) > \rk(F) + |T| \ \ \text{for every proper flat $F$ in $\mathscr{F}$ and every nonempty subset $T$ of $S\setminus \pi^{-1}(F)$.}
\]
\end{dfn}

The Bergman fan $\Sigma_P$ is unimodular with respect to $N_{\widetilde E}$, defining a smooth toric variety $X_P$ over $\mathbb{C}$. We write $A(\Sigma_P)$ for the Chow ring of $X_P$.
We relate $A(\Sigma_P)$ to the \emph{Chow ring} of the polymatroid $P$ (Definition \ref{dfn:ppchow}), denoted $\DP(P)$, introduced in \cite[Section 4]{PP21}. Our main result states the following.

\begin{manual}{\cref{thm:isomorphism}}
 There is a natural isomorphism of graded rings $\DP(P) \cong A(\Sigma_P)$.
\end{manual}

In \cref{cor:gb}, we use \cref{thm:isomorphism} to recover a Gr\"obner basis for $\DP(P)$ found in \cite{PP21}.
In \cref{cor:hodge}, we prove the K\"{a}hler package for $A(\Sigma_P)$ with respect to the cone of strictly convex piecewise linear function on $\Sigma_P$, extending the K\"{a}hler package for $\DP(P)$ with respect to the $\sigma$-cone in  \cite{PP21}. See Remark \ref{ConeRemark} for a comparison of the two cones.

\subsection{Building sets}
In fact, our results hold for polymatroids $P$ equipped with a \word{geometric building set} $\G$ (\cref{sec:buildingsets}). The statements above are specializations of our results to the case when $\G$ consists of all nonempty flats of $P$.
In maximal generality, we define the \word{Bergman fan} of $(P, \G)$, denoted $\Sigma_{P,\G}$ (\cref{dfn:bergmanfan}).
The Chow ring $\DP(P, \G)$ associated to $(P, \G)$ was introduced in \cite{PP21}, and it is isomorphic to
$A(\Sigma_{P, \G})$ (\cref{thm:isomorphism}).
All corollaries continue to hold, including the K\"{a}hler package for $\DP(P, \G)$.

\subsection*{Organization}
In \cref{sec:multisymmetric}, we develop the combinatorics of multisymmetric matroids and lifts, a key tool throughout this paper.
We use lifts to define the Bergman fan of a polymatroid (with respect to a geometric building set) in \cref{sec:buildingsets}.
Finally, in \cref{sec:chow}, we show that the Chow ring of the Bergman fan agrees with the polymatroid Chow ring of \cite{PP21} and derive consequences.
Examples \ref{subspace} and \ref{rmk:why-geometric}, and Remarks \ref{rmk:geomOfBS}, \ref{tropicalization}, and \ref{rmk:geoProof} explain the geometry underlying this work.

\subsection*{Acknowledgements}
We thank Spencer Backman for comments and conversations. The second author is partially supported by a Simons Investigator Grant and NSF Grant DMS-2053308, the third is supported by an NDSEG fellowship, and the last is supported by a Sloan fellowship.

\iffalse
\textcolor{blue}{(Please move this explicit description to the main body.) The Bergman fan has the associated \word{Chow ring}
\[
  A(\Sigma_P) \coloneqq \Z[z_F : F \text{ nonempty proper flat of $P$}] \otimes \R[z_i : i \in \widetilde E] / I,
\]
where $I$ is an ideal generated by the following polynomials:
\begin{align*}
  z_{F_1} &z_{F_2},  \qquad \text{$F_1$ and $F_2$ are incomparable nonempty proper flats of $P$,} \\
  z_F &\prod_{i \in S} z_i,  \quad \text{$S$ is a subset of $\widetilde E$ satisfying $\rk(F \cup \pi(S)) \leq \rk(F) + |S \setminus \pi^{-1}(F)|$,} \\
  \sum_{i \in F} z_F &- \sum_{j \in G} z_G,  \quad \text{$i$ and $j$ are elements of $\widetilde E$.}
\end{align*}}
\fi

\section{Multisymmetric matroids}\label{sec:multisymmetric}
Many proofs in \cref{sec:chow} reduce statements about polymatroids to known statements about matroids.
The key tool for this reduction is multisymmetric matroids, a new cryptomorphic formulation of polymatroids.
Let $\sym_{\widetilde{E}}$ denote the symmetric group on a finite set $\widetilde{E}$. We continue to assume that all  (poly)matroids are loopless.
\begin{dfn}\label{dfn:multisymmetric}
  A \emph{multisymmetric matroid} is a matroid $M$ on $\widetilde{E}$ equipped with a partition
  $\widetilde{E} = \widetilde{E}_1 \sqcup \cdots \sqcup \widetilde{E}_n$
  such that the action of $\Gamma = \sym_{\widetilde{E}_1} \times \cdots \times \sym_{\widetilde{E}_n}$ on $\widetilde{E}$ takes flats to flats.
The \word{geometric part} of a subset $S \subseteq \widetilde{E}$ is $S^{\hspace{0.05ex}\mathrm{geo}} \coloneqq \cap_{\gamma \in \Gamma} (\gamma \cdot S)$. We call a subset $S\subseteq \widetilde{E}$ \word{geometric} if $S = S^{\hspace{0.05ex}\mathrm{geo}}$.
\end{dfn}

For a multisymmetric matroid $M$, we write $\L_M^\Gamma$ for its poset of geometric flats.

\begin{example}
  Any matroid $M$ on $\widetilde{E}$ can be given the \word{trivial} multisymmetric structure by setting $\Gamma = \prod_{e \in \widetilde{E}} \sym_{\{e\}}$. In this case, $\L_M^\Gamma = \L_M$.
\end{example}
\begin{example}\label{ex:uniform}
  If $M$ is multisymmetric on $\widetilde{E} = \widetilde{E}_1$, then $M$ is a uniform matroid, and $\L_M^\Gamma$ is $\{\emptyset, \widetilde{E}\}$.
\end{example}
\begin{example}
  Let $\widetilde{E}$ be the set of edges of the complete graph $K_4$, and let $M$ be the graphic matroid of $K_4$.
  No transposition of $\sym_{\widetilde{E}}$ preserves the flats of $M$, so $M$ has no non-trivial multisymmetric structures.
\end{example}

Closure in a multisymmetric matroid is restricted by the group action.
\begin{lem}\label{cor:union}
  If $M$ is multisymmetric and $S \subseteq \widetilde{E}$ is geometric, then $\cl_M(S)$ is also geometric.
\end{lem}
\begin{proof}
For any $S \subseteq \widetilde{E}$ geometric and $\gamma \in \Gamma$, we have $\gamma \cdot \cl_M(S) = \cl_M(\gamma\cdot S) = \cl_M(S)$.
\end{proof}

\begin{cor}\label{cor:sublattice}
  If $M$ is multisymmetric, then the geometric flats of $M$ form a sublattice of $\L_M$.
\end{cor}
\begin{proof}
  If $F$ and $G$ are two geometric flats, then $\cl_M(F \cup G)$ is geometric by \cref{cor:union}.
  The intersection of two geometric flats is also geometric.
  In other words, the set of geometric flats is closed under both join and meet, and therefore forms a sublattice of $\L_M$.
\end{proof}
\begin{cor}\label{cor:polymat}
  If $M$ is multisymmetric on $\widetilde{E} = \widetilde{E}_1 \sqcup \cdots \sqcup \widetilde{E}_n$, then $\L_M^\Gamma$ is the lattice of flats of the polymatroid $P$ on the set of indices $E = \{1, \ldots, n\}$ defined by the rank function
    $\rk_P(A) \coloneqq \rk_M( \cup_{i \in A} \widetilde{E}_i).$

\end{cor}
\begin{proof}
  Define $\pi: \widetilde{E} \to E$ by setting $\pi^{-1}(i) = \widetilde{E}_i$.
If $F$ is a flat of $P$, then for all $F \subsetneq A \subseteq \{1, \dotsc, n\}$,
\[ \rk_M(\pi^{-1}(F)) = \rk_P(F) < \rk_P(A) = \rk_M(\pi^{-1}(A)). \]
By \cref{cor:union}, $\cl_M(\pi^{-1}(F))$ is geometric, so we conclude that $\cl_M(\pi^{-1}(F)) = \pi^{-1}(F)$.
In other words, $\pi^{-1}(F)$ is a flat of $M$.
Conversely, if $F$ is not a flat of $P$, then $\rk_P(F) = \rk_P(F \cup i)$ for some $i$ not in $F$.
This implies $\cl_M(\pi^{-1}(F)) \supseteq \pi^{-1}(F \cup i)$, so $\pi^{-1}(F)$ is not a flat of $M$.
Therefore, $F$ is a flat of $P$ if and only if $\pi^{-1}(F)$ is a geometric flat of $M$, and we have an isomorphism of lattices
\[ \L_P \longrightarrow \L_M^\Gamma, \qquad F \longmapsto \pi^{-1}(F). \qedhere\]
\end{proof}

\begin{lem}\label{lem:closure}
  Let $M$ be a multisymmetric matroid on $\widetilde{E} = \widetilde{E}_1 \sqcup \cdots \sqcup \widetilde{E}_n$.
\begin{enumerate}\itemsep 5pt
  \item If $S \subseteq \widetilde{E}$, then either $\cl_M(S) \cap \widetilde{E}_i = \widetilde{E}_i$ or $\cl_M(S) \cap \widetilde{E}_i = S \cap \widetilde{E}_i$.\label{lem:extraclosure}
  \item If $F$ is a flat of $M$, then $\rk_M(F) = \rk_M(F^{\hspace{0.05ex}\mathrm{geo}}) + |F \setminus F^{\hspace{0.05ex}\mathrm{geo}}|$. \label{lem:closurerank}
\end{enumerate}
\end{lem}

\begin{proof}
    A permutation of $\widetilde{E}_i \setminus S$ induces an automorphism of $M$ that fixes $S$. Any such automorphism also fixes $\cl_M(S)$ because automorphisms commute with closure.
    Hence, if $(\widetilde{E}_i \setminus S) \cap \cl_M(S)$ is nonempty, then $\widetilde{E}_i \subseteq \cl_M(S)$. This proves the first part.

      For the second part, let $\{s_1, \ldots, s_k\} = F \setminus F^{\hspace{0.05ex}\mathrm{geo}}$. If \cref{lem:closurerank} fails then for some $1 \leq i \leq k,$
       \[
         \rk_M(F^{\hspace{0.05ex}\mathrm{geo}}\cup \{s_1, \ldots, s_i\}) = \rk_M(F^{\hspace{0.05ex}\mathrm{geo}}\cup \{s_1, \ldots, s_{i+1}\}).
       \]
       Therefore if $s_{i+1} \in \widetilde{E}_j$ then $\widetilde{E}_j \subseteq \cl_M(F^{\hspace{0.05ex}\mathrm{geo}}\cup \{s_1, \ldots, s_i\}) \subseteq F$ by the first part.
        Consequently, $s_{i+1} \in F^{\hspace{0.05ex}\mathrm{geo}}$, a contradiction.
\end{proof}

\subsection{Lifts}
  Let $P$ be a polymatroid on $E = \{1, \ldots, n\}$ and
   $M$ a multisymmetric matroid on $\widetilde{E}~=~\widetilde{E}_1~\sqcup~\cdots~\sqcup~\widetilde{E_n}$. If $P$ is the polymatroid given by Corollary \ref{cor:polymat}, then we say that $M$ is a \word{multisymmetric lift} of $P$.
  If $\rk_M(\widetilde{E}_i) = |\widetilde{E_i}|$ for all $1 \leq i \leq n$, then we say that $M$ is a \word{minimal} multisymmetric lift of $P$.
\begin{thm}\label{lift}
A polymatroid $P$ has a unique minimal multisymmetric lift $\widetilde P$ constructed as follows.
  For $1 \leq i \leq n$, let $\widetilde E_i = \{1, \ldots, \rk_P(i)\}$, and
$ \widetilde E = \widetilde E_1 \sqcup \cdots \sqcup \widetilde E_n$. Define the projection $\pi: \widetilde E \to E$ by $\pi^{-1}(i) = \widetilde E_i$.

  The minimal multisymmetric lift of $P$ is the matroid $\widetilde P$ on $\widetilde E$ with rank function
  \[
    \rk_{\widetilde P}(S) = \min \{ \rk_P(A) + |S \setminus \pi^{-1}(A)| : A \subseteq E \}.
\]
\end{thm}
\begin{example}\label{rmk:why-geometric}
  Suppose $P$ is the polymatroid realized by an arrangement of subspaces $V_1,\ldots, V_n$ in $V$ as in \cref{subspace}.
  The minimal multisymmetric lift $\widetilde P$ is realized by any hyperplane arrangement
  \[ \{V_{i,j} : 1 \leq i \leq n, \; 1 \leq j \leq \rk_P(i)\}, \]
  where $V_{i,1}, \ldots, V_{i, \rk_P(i)}$ are generic hyperplanes containing $V_i$.
  Under the correspondence between flats of $\widetilde P$ and intersections of the $V_{i,j}$'s, the geometric flats correspond to the subspaces of $V$ arise in for any choice of $\{V_{i,j}\}_{i,j}$.
\end{example}
  Versions of the construction in \cref{lift} make many independent appearances in the literature, for example, in \cite[\S2]{H72} and \cite[Propositions 3.1 and 3.2]{L77}. The most complete treatment we are aware of is \cite[\S2]{N86}, whose terminology differs from ours.
\subsubsection*{Notation.} We continue to use the notations of \cref{lift} in the remainder of this section. For visual clarity, we often write $M = \widetilde P$.
As usual, $\Gamma$ stands for the product of symmetric groups that acts on $M$.

\medskip
\begin{prop}\label{liftwelldefined}
The minimal lift $M = \widetilde P$ is a multisymmetric matroid. Explicitly,
\begin{enumerate}\itemsep 5pt
\item $\rk_{M}$ is a matroid rank function, and
\item the action of $\Gamma$  preserves flats.
\end{enumerate}
\end{prop}
\begin{proof}
 The proof of the first part is reproduced from \cite[Proof of Theorem 11.1.9]{O11}. We need to check that $\rk_M$ is non-negative, increasing, submodular, and satisfies $\rk_M(S)  \leq |S|$ for $S \subseteq \widetilde E$.

It is clear that values of $\rk_M$ are non-negative.
Let $S \subseteq \widetilde E$ and $s \in \widetilde E \setminus S$. For all $A \subseteq E$,
\[
\rk_P(A) + |S \setminus \pi^{-1}( A)| \leq \rk_P(A) + |S\cup a \setminus \pi^{-1}( A)|,
\]
so $\rk_M$ is increasing. Moreover, by induction on $|S|$,
\begin{align*}
  \rk_M(S \cup s)
  =& \min\{\rk_P(A) + |S \cup s \setminus \pi^{-1}(A)| : A \subseteq E\}\\
  \leq& \min\{\rk_P(A) + |S \setminus \pi^{-1}(A)| : A \subseteq E\} + 1 = \rk_M(S) + 1\leq |S|+1 = |S \cup s|,
\end{align*}
so it only remains to check that $\rk_M$ is submodular.
Let $S_1, S_2 \subseteq \widetilde E$ and
  $A_1, A_2 \subseteq E$ such that $\rk_M(S_i) = \rk_{P}(A_i) +
  |S_i \setminus \pi^{-1}(A_i)|$. Then
\begin{align*}
\rk_M&(S_1) + \rk_M(S_2)
= \rk_{P}(A_1) + |S_1 \setminus \pi^{-1}(A_1)| + \rk_{P}(A_2) + |S_2 \setminus \pi^{-1}(A_2)|  \\
\geq& \rk_{P}(A_1 \cup A_2) + \rk_{P}(A_1 \cap A_2) + |S_1 \setminus \pi^{-1}(A_1)| + |S_2 \setminus \pi^{-1}(A_2)| \quad \text{by submodularity of $\rk_{P}$} \\
\geq& \rk_{P}(A_1 \cup A_2) + \rk_{P}(A_1 \cap A_2) + |(S_1 \cup S_2) \setminus \pi^{-1}(A_1 \cup A_2)| + |(S_1 \cap S_2) \setminus \pi^{-1}(A_1 \cap A_2)| \\
\geq& \rk_M(S_1 \cup S_2)  + \rk_M(S_1 \cap S_2).
\end{align*}
For the second part, let $\gamma \in \Gamma$ and $S \subseteq \widetilde E$.
  For any $A \subseteq E$, we have
   $|S \setminus \pi^{-1}(A)| = |(\gamma \cdot S) \setminus \pi^{-1}(A)|$, and hence
   $\rk_M(S) = \rk_M(\gamma \cdot S)$.
\end{proof}
The following lemma implies $\widetilde P$ is a multisymmetric lift of $P$.
\begin{lem}\label{Gfixed}
  Fix notation as in \cref{lift}. If $S \subseteq \widetilde E$ is stable under $\Gamma$, then
    $\rk_M(S) = \rk_P(\pi(S)).$
\end{lem}
\begin{proof}
  Since $S$ is stable under the action of $\Gamma$, $S$ is a union of
  fibers of $\pi$. Hence,
    \begin{align*}
      \rk_{M}(S) =& \min \{ \rk_P(A) + |S \setminus \pi^{-1}(A)| : A \subseteq E \}\quad \text{by the definition of $\rk_{M}$}\\
      =& \min \{ \rk_P(A) + |S \setminus \pi^{-1}(A)| : A \subseteq \pi(S) \} \quad \text{because $\rk_P$ is increasing}\\
      =& \min \bigg\{ \rk_P(A) + \sum_{i\in \pi(S)\setminus
      A}|\pi^{-1}(i)| : A \subseteq \pi(S) \bigg\}\quad \text{because
      $S$ is a union of fibers of $\pi$}\\
      = & \min \bigg\{ \rk_P(A) + \sum_{i\in \pi(S)\setminus A}\rk_P(i)
      : A \subseteq \pi(S) \bigg\} \quad\text{by the definition of $\pi$
      in \cref{lift}}\\
      \geq& \rk_P(\pi(S))\quad \text{by submodularity of $\rk_P$.}
    \end{align*}
    On the other hand, taking $A=\pi(S)$ in the definition of $\rk_M(S)$, we have $\rk_{M}(S)\leq \rk_P(\pi(S))$.
\end{proof}
\begin{proof}[Proof of \cref{lift}]
  Let $P$ be a polymatroid, and let $M = \widetilde P$ on $\widetilde E = \widetilde E_1 \sqcup \cdots \sqcup \widetilde E_n$, acted upon by $\Gamma$, be as in the statement of \cref{lift}.
  By \cref{liftwelldefined} and \cref{Gfixed}, $M$ is a multisymmetric lift of $P$. It is minimal because
  \[
    \rk_M(\widetilde E_i) = \rk_P(\pi(\widetilde E_i)) = \rk_P(i) = |\widetilde E_i|.
   \]
   It remains to check the uniqueness statement.

   Suppose $M'$ is a second minimal multisymmetric lift of $P$. Then there is an equivariant bijection between their ground sets; hence, we may think of both $M$ and $M'$ as matroids on $\widetilde E$.
   By \cref{Gfixed}, $\rk_M(A) = \rk_{M'}(A)$ for any $A \subseteq \widetilde E$ stable under the $\Gamma$-action.
   By \cref{lem:closure}\cref{lem:closurerank}, a multisymmetric matroid is determined by the ranks of its $\Gamma$-stable sets, so $M \cong M'$.
\end{proof}

\subsection{Operations}
The formation of minimal multisymmetric lifts commutes with some polymatroid operations.
Let $F$ be a flat of a polymatroid $P$ on $E$. The \word{restriction} of $P$ to $F$, denoted $P|_F$, is the polymatroid on $F$ with rank function
\[
  \rk_{P|_F}(A) \coloneqq \rk_P(A), \quad A \subseteq F.
\]
There is a lattice isomorphism
\[ \L_{P|_F} \to \{G \in \L_P: G \leq F\},\quad H \mapsto H. \]

If $P_1$ and $P_2$ are polymatroids on $E^1$ and $E^2$, respectively, then their \textit{direct sum} is the polymatroid $P_1 \oplus P_2$ on $E^1 \sqcup E^2$ with rank function
\[
  \rk_{P_1 \oplus P_2}(S) = \rk_{P_1}(S \cap E^1) + \rk_{P_2}(S \cap E^2).
\]
There is a lattice isomorphism
\[
  \L_{P_1} \times \L_{P_2} \to \L_{P_1 \oplus P_2}, \quad (F, G) \mapsto F \sqcup G.
  \]

\begin{lem}\label{lem:sum}
  If $P_1$ and $P_2$ are polymatroids, then
    $\widetilde{P_1 \oplus P_2} = \widetilde P_1 \oplus \widetilde P_2.$
\end{lem}
\begin{proof}
  From Corollary \ref{cor:polymat} and the definition of $P_1 \oplus P_2$, it follows that $\widetilde P_1 \oplus \widetilde P_2$ lifts $P_1 \oplus P_2$.
  If $\Gamma_1$ and $\Gamma_2$ are the groups acting on $\widetilde P_1$ and $\widetilde P_2$, then $\Gamma_1 \times \Gamma_2$ acts on $\widetilde P_1 \oplus \widetilde P_2$, so it is multisymmetric.
  It is minimal by minimality of $\widetilde P_1$ and $\widetilde P_2$, so the lemma holds by the uniqueness statement of \cref{lift}.
\end{proof}

\begin{lem}\label{lem:restriction}
  If $F$ is a flat of a polymatroid $P$, then
  $\widetilde{P|_F} = \widetilde P|_{\pi^{-1}(F)}.$
\end{lem}
\begin{proof}
  The rank functions of both sides are obtained by restriction, so $\widetilde P|_{\pi^{-1}(F)}$ lifts $\widetilde{P|_F}$.
  If $\Gamma$ acts on $\widetilde P$, then a subset of $\Gamma$'s factors acts on $\widetilde{P}|_{\pi^{-1}(F)}$, so it is multisymmetric.
  It is minimal because $\widetilde P$ is, so \cref{lift} implies the lemma.
\end{proof}

\section{Building sets and Bergman fans}\label{sec:buildingsets}
Here, we recall the combinatorics of geometric building sets. We then generalize the definition of the Bergman fan of a polymatroid given in the introduction by associating a fan to each polymatroid equipped with a geometric building set.
\subsection{Geometric building sets}
Let $P$ be a polymatroid.
If $\mathscr{G} \subseteq \L_P$ and $F$ is a flat, let
\[ \mathscr{G}_{\leq F} \coloneqq \{ G \in \mathscr{G} : G \leq F\} \]
and write $\max \mathscr{G}$ for the set of maximal elements of $\mathscr{G}$.
\begin{dfn}
  A \word{geometric building set}\footnote{Our terminology is that of \cite[Definition 4.4]{FK04}. In \cite{PP21}, geometric building sets are called ``combinatorial building sets'', but \cite{FK04} uses this term for collections that satisfy only the isomorphism condition.
  The ``combinatorial building sets'' of \cite{FK04} are the same as the ``building sets'' of \cite{FY04}. Postnikov considers only building sets for the Boolean lattice in \cite[Definition 7.1]{P09}; all such building sets are geometric.} of a polymatroid $P$ is a collection $\G$ of nonempty flats such that for all $F \in \L_P \setminus \{\emptyset\}$, the map
  \begin{equation*}
    \prod_{G \in \max \G_{\leq F}} \L_{P|_G} \to \L_{P|_F} \label{factorization}
    \end{equation*}
  is an isomorphism, and
  \begin{equation*}
 \sum_{G \in \max \G_{\leq F}} \rk(G) = \rk(F). \label{combobuildingset}
\end{equation*}
  If $E \in \G$, then a \word{nested set} of $P$ with respect to $\G$ is a subset $\N \subseteq  \G$ such that for all $\{F_1, \ldots, F_k\} \subseteq \N$ pairwise incomparable with $k \ge 2$, we have that
  \[ \cl_P(F_1 \cup \cdots \cup F_k)  \not \in \G.\]
\end{dfn}
With respect to a fixed geometric building set, a subset of a nested set is nested, so nested sets form a simplicial complex.
All building sets are henceforth assumed to contain $E$\footnote{In the realizable case, this assumption guarantees the associated wonderful compactification is smooth and can be described as an iterated blow-up \cite[\S4.1]{dCP95}.  Combinatorially, we lose nothing by this assumption \cite[Remark 4.1]{PP21}.}.
\begin{example}
  The maximal geometric building set of $P$ is the collection of all nonempty flats.
  With respect to this building set, the nested sets are flags of nonempty flats.
\end{example}

\begin{lem}\label{liftedbs}
  Let $P$ be a polymatroid, $\widetilde P$ its minimal multisymmetric lift, and $\pi$ as in \cref{lift}.
  If $\mathscr G$ is a geometric building set for $P$, then
\[ \widetilde{\mathscr G} = \{ \pi^{-1}(G) : G \in \mathscr G\} \cup \{\text{atoms of $\mathscr L_{\widetilde P}$}\} \]
is a geometric building set for $\widetilde P$.
\end{lem}
\begin{proof}
Let $F$ be a flat of $M = \tilde P$.
By \cref{lem:sum} and \cref{lem:restriction}, the map
\[ \prod_{G \in \max \widetilde \G_{\leq F^{\hspace{0.05ex}\mathrm{geo}}}} \L_{M|_G} \to \L_{M|_F} \]
  factors into a chain of isomorphisms
  \[
    \prod_{G \in \max \widetilde \G_{\leq F^{\hspace{0.05ex}\mathrm{geo}}}} \L_{M|_G}
    \cong \prod_{H \in \max \G_{\leq \pi(F^{\hspace{0.05ex}\mathrm{geo}})}} \L_{\widetilde{P|_H}}
    \cong \L_{\widetilde{P|_{\pi(F^{\hspace{0.05ex}\mathrm{geo}})}}} \cong \L_{M|_{F^{\hspace{0.05ex}\mathrm{geo}}}},
\]
and by \cref{Gfixed},
\[
  \sum_{G \in \max \widetilde \G_{\leq F^{\hspace{0.05ex}\mathrm{geo}}}} \rk_{M}(G)
  = \sum_{G \in \max \widetilde \G_{\leq F^{\hspace{0.05ex}\mathrm{geo}}}} \rk_{P}(\pi(G))
  = \sum_{H \in \max \G_{\leq \pi(F^{\hspace{0.05ex}\mathrm{geo}})}} \rk_{P}(H)
  = \rk_P(\pi(F^{\hspace{0.05ex}\mathrm{geo}})) = \rk_M(F^{\hspace{0.05ex}\mathrm{geo}}).
\]
Consequently, by \cref{lem:closure}\cref{lem:closurerank},
\[
  \prod_{G \in \max \widetilde\G_{\leq F}} \L_{M|_G}
   = \prod_{G \in \max \widetilde\G_{\leq F^{\hspace{0.05ex}\mathrm{geo}}}}  \L_{M|_G} \times \prod_{i \in F \setminus F^{\hspace{0.05ex}\mathrm{geo}}}  \L_{M|_{\{i\}}}
  \cong  \L_{M|_{F^{\hspace{0.05ex}\mathrm{geo}}}} \times \prod_{i \in F \setminus F^{\hspace{0.05ex}\mathrm{geo}}} \L_{M|_{\{i\}}} \cong \L_{M|_F},
\]
and
\[
  \sum_{G \in \max \widetilde{\mathscr{G}}_{\leq F}} \rk_{M}(G)
  = \rk_M(F^{\hspace{0.05ex}\mathrm{geo}}) + \sum_{i \in F \setminus F^{\hspace{0.05ex}\mathrm{geo}}} \rk_{M}(i)
  = \rk_M(F^{\hspace{0.05ex}\mathrm{geo}})+|F \setminus F^{\hspace{0.05ex}\mathrm{geo}}|
  = \rk_M(F),
\]
as desired.
\end{proof}
\begin{lem}\label{lem:nested}
  If $\G$ is a geometric building set for $P$, then $\N \subseteq \G$ is $\G$-nested if and only if $\widetilde \N = \{ \pi^{-1}(F) : F \in \N\}$ is $\widetilde \G$-nested.
\end{lem}
\begin{proof}
  Set $M = \widetilde P$.
  Let $\{F_1, \ldots, F_k\}\subseteq \mathscr N$ be a set of pairwise incomparable flats.
  By \cref{cor:sublattice}, the geometric flats of $M$ form a sublattice of $\L_M$, isomorphic to $\L_P$ by \cref{lift}. Hence,
  \[
    \cl_M(\pi^{-1}(F_1) \cup \cdots \pi^{-1}(F_k)) \in \widetilde \G  \iff  \cl_P(F_1 \cup \cdots \cup F_k) \in \G.  \qedhere
  \]
\end{proof}

\begin{rmk}[Geometry of building sets]\label{rmk:geomOfBS}
  Let $\{V_i\}_i$ be a subspace arrangement in $V$, defining a polymatroid $P$.
  Let $\G$ be a building set for $P$.
  By \cite[\S1.6]{dCP95}, the \textit{wonderful compactification} of $\{V_i\}_i$ with respect to $\G$ can be constructed by blowing up $\P(V)$ along all subspaces $\cap_{i \in F} \P(V_i)$ with $F \in \G$, first blowing up those of dimension 0, then those of dimension 1, and so on.
Let $\{V_{ij}\}_{ij}$ be a hyperplane arrangement realizing $\widetilde P$, as in \cref{rmk:why-geometric}.
  Blowing up a codimension 1 subvariety is an isomorphism, so \cref{lem:nested} implies that the wonderful compactification of $\{V_i\}_i$ with respect to a building set $\G$ is isomorphic to the wonderful compactification of $\{V_{ij}\}_{ij}$ with respect to $\widetilde \G$.
\end{rmk}

\subsection{Bergman fans}
Let $P$ be a polymatroid on $E$, with minimal lift $\widetilde{P}$ on $\widetilde E$.
Let $\R^{\widetilde E}$ be the vector space spanned by $\e_i$ for $i \in \widetilde E$, and write $\e_S \coloneqq \sum_{i \in S} \e_i$ for $S \subseteq \widetilde E$.
If $\mathscr S \subseteq 2^{\widetilde{E}}$ is a collection of subsets, write
\[
  \sigma_{\mathscr S} \coloneqq \cone(\e_S : S \in \mathscr S) \subseteq \R^{\widetilde E} / \R (1,1,\ldots,1).
\]
\begin{dfn}\label{dfn:bergmanfan}
  Let $P$ be a polymatroid on $E$, and $\G \subseteq \L_P$ a geometric building set of $P$.
  The \word{Bergman fan} associated to $(P, \G)$ is
  $\Sigma_{P, \G} \coloneqq \{\sigma_{\N}\}_\N$,
  where $\N$ ranges over all $\widetilde\G$-nested sets of $\widetilde P$ such that $\widetilde E \not \in \N$.
\end{dfn}
If $P$ is a matroid and $\G$ is the maximal geometric building set of $P$, then $\widetilde P = P$ and $\widetilde \G = \G$. In this case, $\Sigma_{P,\G}$ coincides with the Bergman fan of \cite[Definition 3.2]{AHK18}.

\begin{lem}
  If $P$ is a polymatroid and $\G$ is the maximal geometric building set, then $\Sigma_{P, \G}$ coincides with the Bergman fan $\Sigma_P$ defined in the introduction (\cref{dfn:maximalBergman}).
\end{lem}

\begin{proof}
  Fix a proper flat $F$ of $P$ and a subset $S \subseteq \widetilde{E}$. It suffices to show that the closure of $\pi^{-1}(F) \cup T$ is not geometric for all nonempty $T \subseteq S \setminus \pi^{-1}(F)$ if and only if the following inequality holds for all nonempty $T \subseteq S \setminus \pi^{-1}(F)$,
  \[
  \rk_P(F \cup \pi(T)) >\rk_P(F) + |T|.
  \]
  First, assume the displayed inequality holds for all such $T$. The inequality's left-hand side is the rank of the smallest geometric flat containing $\pi^{-1}(F) \cup T$. Because the right-hand side is at least $\rk_P(\pi^{-1}(F) \cup T)$, the closure of $\pi^{-1}(F) \cup T$ is not geometric.

  Conversely, assume there is a nonempty $T\subseteq S \setminus \pi^{-1}(F)$ for which the inequality fails. We claim there is a nonempty $T'\subseteq S \setminus \pi^{-1}(F)$ such that the closure of $\pi^{-1}(F)\cup T'$ is geometric. If $\rk_P(F \cup \pi(T)) = \rk_P(F) + |T|$, then take $T'=T$. Otherwise, if $\rk_P(F \cup \pi(T)) < \rk_P(F) + |T|$, then there is a circuit $C \subseteq \pi^{-1}(F) \cup T$ such that $C \setminus \pi^{-1}(F)$ is nonempty. Writing $C \setminus \pi^{-1}(F) = \{c_1,\ldots,c_k\}$, we have $c_i \in \cl_{\widetilde{P}}(\pi^{-1}(F) \cup C \setminus c_i)$ for all $i$ (see \cite[Lemma 1.4.2]{O11}). Therefore by symmetry, $\pi^{-1}(\pi(c_i)) \subseteq \cl_{\widetilde{P}}(\pi^{-1}(F) \cup C)$ for all $i$, so the closure of $\pi^{-1}(F) \cup C$ is geometric. Now we can take $T'=C\setminus \pi^{-1}(F)$.
\end{proof}

\medskip
\begin{lem}\label{polytope}
  If $P$ is a polymatroid and $\G$ is a geometric building set of $P$, then $\Sigma_{P, \G}$ is a subfan of the normal fan of a convex polytope.
\end{lem}
\begin{proof}
  Let $E$ be the ground set of $P$ and $B$ be the Boolean polymatroid with $\rk_B(i) = \rk_P(i)$ for all $i \in E$.
  Suppose $\G$ is a building set for $P$. The definition of geometric building set implies for any flat $F$ of $P$, $\max \G_{\leq F}$ is a partition of $F$. Hence, $\widetilde \G$ is a building set of both $\widetilde B$ and $\widetilde P$, and there is an inclusion of fans
  $\Sigma_{P, \G} = \Sigma_{\widetilde P, \widetilde \G} \subseteq \Sigma_{\widetilde B, \widetilde \G}$ (In fact, if all single-element subsets of $E$ are flats of $P$, then $\G$ is a geometric building set for $B$, and we obtain $\Sigma_{P,\G} \subset \Sigma_{\widetilde B, \widetilde \G} = \Sigma_{B, \G}$).
  The lattice of flats of $\widetilde B$ is isomorphic to the lattice of subsets of $\widetilde E$,  so $\Sigma_{\widetilde B, \widetilde \G}$ is the normal fan of a convex polytope by \cite[Theorem 7.4]{P09}.
\end{proof}

By \cite{Bri96}, for any polymatroid $P$ and a geometric building set $\G$, $A(\Sigma_{P,\G})$ has the following presentation.
\begin{prop}\label{fanchow}
  The \word{Chow ring} of $\Sigma_{P,\G}$ satisfies
  \[
    A(\Sigma_{P,\G}) = \Z[z_G : G \in \widetilde \G \setminus \{\widetilde E\}] / I_{P, \G}
  \]
  where $I_{P, \G}$ is the ideal generated by
  \begin{align*}
      z_{G_1} \cdots z_{G_k} \quad &\text{for any not $\widetilde \G$-nested collection $\{G_1, \ldots, G_k\}$},\\
      \sum_{i \in G} z_G - \sum_{j \in F} z_F \quad &\text{for any } i,j \in \widetilde E.
  \end{align*}
In particular, the Chow ring of $\Sigma_P$ satisfies
\[
  A(\Sigma_P) = \Z[z_F : F \text{ nonempty proper flat of $P$}] \otimes \Z[z_i : i \in \widetilde E] / I_P,
\]
where $I_P$ is an ideal generated by the following polynomials where $z_\emptyset$ is replaced by 1 wherever it appears:
\begin{align*}
  z_{F_1} z_{F_2},  \quad &\text{$F_1$ and $F_2$ are incomparable proper flats of $P$,} \\
  z_F \prod_{i \in T} z_i,  \quad &\text{$F$ is a proper flat and  $T\subseteq \widetilde E\setminus \pi^{-1}(F)$ is nonempty satisfying $\rk(F \cup \pi(T)) \leq \rk(F) + |T|$,} \\
  \sum_{i \in F} z_F - \sum_{j \in G} z_G,  \quad &\text{$i$ and $j$ are elements of $\widetilde E$.}
\end{align*}

\end{prop}

We identify $A^1(\Sigma_{P,\G})$ with the space of piecewise linear functions on the support of $\Sigma_{P, \G}$, modulo global linear functions. Explicitly, a piecewise linear function $\ell$ is a representative of
\[
  \sum_{\mathbf{u}_F} \ell(\mathbf{u}_F) z_F \in A^1(\Sigma_{P, \G}),
\]
where the sum is over all primitive ray generators $\mathbf{u}_F$ of rays of $\Sigma_{P, \G}$.

\begin{rmk}\label{fypresentation}
  A slightly different presentation of $A(\Sigma_{P, \G})$ is used by \cite{FY04}.
  For $F \in \widetilde \G$, set
    \[ y_F = \begin{cases} -\sum_{i \in G} z_G, &F = \widetilde E \\ z_F, &\text{otherwise.} \end{cases} \]
    In terms of the $y_F$'s, $A(\Sigma_{P, \G})$ is defined by the ideal $I_{FY}$ generated by
  \begin{align*}
    y_{G_1} \cdots &y_{G_k}, \quad \text{$\{G_1, \ldots, G_k\}$ not $\widetilde \G$-nested and }\\
    \sum_{i \in G}& y_G, \quad i \in \widetilde E.
  \end{align*}
\end{rmk}
\begin{rmk}[Tropicalization and the Bergman fan]\label{tropicalization}
  Suppose $P$ is the matroid realized by an essential hyperplane arrangement $V_1,\ldots,V_n$ in $V$, defined by linear functionals $\ell_1,\ldots,\ell_n$.
  The inclusion
  \[ \P(V \setminus \cup_i V_i) \hookrightarrow \P((\mathbb{F}^*)^n), \quad v \mapsto [\ell_1(v): \ldots: \ell_n(v)] \]
  shows that $\P(V \setminus \cup_i V_i)$ is a very affine variety. Its tropicalization is the support of $\Sigma_P$, by \cite[\S9.3]{S02} and \cite[\S3]{AK06} (see also \cite[Theorem 4.1]{FS05}).
  The corresponding statement for realizable polymatroids does not make sense: the complement of a subspace arrangement may not be very affine, so tropicalization  \cite[Definition 3.2.1]{MS15} is not defined. Nevertheless, if $\{V_i\}_i$ is a subspace arrangement (so $P$ is a polymatroid), then generic hyperplanes $\{V_{ij}\}_{ij}$ realizing $\widetilde P$ as in \cref{rmk:why-geometric} define a subtorus of $\P(\prod_i V/V_i)$. Tropicalizing this torus's intersection with $\P(V \setminus \cup_i V_i)$
  gives the support of $\Sigma_{\widetilde P}$, which coincides with that of $\Sigma_P$.

\end{rmk}

\section{Chow rings of polymatroids}\label{sec:chow}
If $\cS$ is a collection of sets, write $\cup \cS$ for the union of the elements of $\cS$.
In \cite[Section 4]{PP21}, Pagaria and Pezzoli define the Chow ring of a polymatroid as follows\footnote{Unlike \cite{PP21}, we use $\Z$-coefficients.}.

\begin{dfn}\label{dfn:ppchow}
 Let $P$ be a polymatroid and $\G \subseteq \L_P$ a geometric building set.
 The \word{Chow ring} of $(P, \G)$ is
 \[
   \DP(P, \G) \coloneqq \Z[x_F : F \in \G] / I,
 \]
 where $I$ is the ideal generated by
\[
x_{G_1} \cdots x_{G_k} \Big(\sum_{H \geq G} x_H \Big)^b
\]
for $G \in \G$, $\mathscr S = \{G_1, \ldots, G_k\} \subseteq \G$, and
$b \geq \rk_P(G) - \rk_P( \cup \mathscr S_{<G} ).$
\end{dfn}
When $P$ is realizable by an arrangement of subspaces, $\DP(P, \G)$ is the cohomology ring of De Concini \& Procesi's wonderful compactification of the arrangement complement.
In this section, we make use of the theory of Gr\"obner bases. For background on this subject, see \cite[Chapter 15]{E95}.

\begin{thm}\label{thm:isomorphism}
  Let $P$ be a polymatroid on $E$, $\widetilde E$ the ground set of $\widetilde P$, and $\pi: \widetilde E \to E$ the projection.
  There is an isomorphism $\operatorname{DP}(P, \mathscr G) \cong A(\Sigma_{P, \G})$ sending $x_F$ to $y_{\pi^{-1}(F)}$.
\end{thm}
\begin{proof}
    Let $I_{DP} \subseteq \Z[x_F : F \in \G]$ be the defining ideal of
  $\DP(P, \G)$, and let $I_{FY} \subseteq \Z[y_F : F \in \widetilde \G]$ be the defining ideal of $A(\Sigma_{P,\G})$ as in \cite{FY04} (See
  \cref{fypresentation}). We define the following map on polynomial rings.
  \[
   \varphi: \Z[x_F : F \in \G] \to \Z[y_F : F \in \widetilde \G], \quad x_F \mapsto y_{\pi^{-1}(F)}
  \]
    First we show that $\varphi(I_{DP}) \subseteq I_{FY}$. Write
  $f \in I_{DP}$ for one of the defining relations of $I_{DP}$:
\[
  f = \Big(\prod_{F \in \mathscr S} x_F \Big)\Big(\sum_{\mathscr G \ni
  H \geq G} x_H \Big)^b.
  \]
  By \cite[Theorems 1 and 3]{FY04}, $I_{FY}$ contains the following two types of polynomials:
  \begin{align*}
\prod_{F \in \mathscr S}& y_F, \qquad \mathscr S \text{ not $\widetilde{\mathscr{G}}$-nested,} \\
    \prod_{F \in \N} y_F \Big( \sum_{H \geq G}& y_G \Big)^d, \quad \text{$\N$ a nested antichain, $\cup \N < G$, and $d = \rk(G) - \rk(\cup N)$.}
\end{align*}
If $\mathscr S$ is not $\mathscr G$-nested, then $\tilde{\mathscr{S}} \coloneqq \{ \pi^{-1}(F) : F \in \mathscr S\}$ is not $\tilde{\mathscr G}$-nested by \cref{lem:nested}. Hence, $\varphi(f)$ is divisible by a relation of the first type.
Otherwise, $\mathscr S$ is $\mathscr{G}$-nested, so $\tilde {\mathscr {S}}$ is $\widetilde{\mathscr{G}}$-nested.
In this case, $\varphi(f)$ is divisible by a relation of the second type because
\begin{align*}
 b &\geq
\rk_P(G) - \rk_P(\cup \mathscr S_{<G})
  = \rk_{\tilde P}(\pi^{-1}(G)) - \rk_{\tilde P}(\cup \mathscr{\tilde S}_{< \pi^{-1}(G)})
  = \rk_{\tilde P}(\pi^{-1}(G)) - \rk_{\tilde P}(\cup \max \mathscr{\tilde S}_{<\pi^{-1}(G)}).
\end{align*}
This proves that $\varphi(I_{DP}) \subseteq I_{FY}$, so $\varphi$ descends to
$\bar \varphi: \DP(P,\G) \to A(\Sigma_{P,\G})$.

If $F \in \widetilde \G$ is a flat of rank greater than 1, then $y_{F}$ is in the image of $\bar \varphi$.
By the linear relation $\sum_{i \in G} y_G = 0$, it follows that $y_i$ is also in the image of $\bar\varphi$. Therefore, $\bar\varphi$ is surjective.
It remains to show that $\bar \varphi$ is injective.
By \cite[Theorem 2]{FY04}, the generators of $I_{FY}$ in the previous paragraph
are a Gr\"obner basis with respect to any lexicographic monomial order $<$ in which $F_1 \subseteq F_2$ implies $y_{F_1} > y_{F_2}$.
Any such order is an elimination order with respect to $\{y_i: i \in \widetilde E\}$.
By \cite[Proposition 15.29]{E95}\footnote{Eisenbud's proof of this statement works over $\Z$ because all leading coefficients in our Gr\"obner basis are 1.\label{zgbnote}},
the generators of $I_{FY}$ in the previous paragraph
that do not involve any $y_i$, $i \in \widetilde E$, are a Gr\"obner basis for $\im(\varphi) \cap I_{FY}$.
Any such polynomial is the image of a generator of $I_{DP}$, so $\varphi^{-1}(I_{FY}) = I_{DP}$.
This implies $\bar \varphi$ is an isomorphism.
\end{proof}
The proof of \cref{thm:isomorphism} also shows the following, originally obtained from \cite[Corollary 2.8]{PP21}.
\begin{cor}\label{cor:gb}
  The defining relations of $\DP(P, \G)$ in
  \cref{dfn:ppchow} form a Gr\"obner basis with respect to any lexicographic $<$ such that $x_{F_1} < x_{F_2}$ whenever $F_1 \supsetneq F_2$.
\end{cor}
This recovers the monomial basis of \cite{PP21}.
\begin{cor}\label{cor:basis}
  The following monomials are a $\Z$-basis of $\DP(P,\G)$:
  \[
    x_{G_1}^{a_1} \cdots x_{G_k}^{a_k}
  \]
  where $\N = \{G_1, \ldots, G_k\}$ is a nested set of $\G$, and
  \[ 1 \leq a_i < \rk(F_i) - \rk(\cup \N_{<F_i})\quad \text{for all $1 \leq i \leq k$.} \]

\end{cor}
\begin{proof}
Immediate from \cref{thm:isomorphism} and \cite[Corollary 1]{FY04} (or \cref{cor:gb} and \cite[Theorem 15.3]{E95}).
\end{proof}
\begin{rmk}
  In \cite[Corollary 2.8]{PP21}, degree-lexicographic order (also called ``homogeneous lexicographic'' or ``graded lexicographic'' order) is used.
  Since the defining relations of $\DP(P,\G)$ in
  \cref{dfn:ppchow} are all homogeneous, their initial terms with respect to the lex or degree-lex orders are the the same.
  Consequently, they are a Gr\"obner basis with respect to one order if and only if they are with respect to the other.
\end{rmk}
\begin{rmk}\label{rmk:geoProof}
      Suppose $P$ is a polymatroid with building set $\G$, realized by an arrangement $\{V_i\}_i$ in $V$. Let $\{V_{ij}\}_{ij}$ realize $\widetilde P$ as in \cref{rmk:why-geometric}. In this case, an alternate proof of \cref{thm:isomorphism} is possible.
  By \cite[\S4]{FY04}, $A(\Sigma_{P, \G})$ is isomorphic to the Chow ring of the wonderful compactification of $\{V_{ij}\}_{ij}$ with respect to $\widetilde \G$.
  By \cref{rmk:geomOfBS}, the wonderful compactifications of $\{V_{ij}\}_{ij}$ with respect to $\widetilde \G$ and $\{V_i\}_i$ with respect to $\G$ are isomorphic. Hence, their Chow rings are isomorphic. The Chow ring of the latter space is isomorphic to $\DP(P, \G)$ by a comparison of the presentations in \cref{dfn:ppchow} and \cite[Theorem 5.2]{dCP95}.
\end{rmk}
In the remainder of this section, we recover \cite[Theorem 4.7]{PP21} and generalize \cite[Theorem 4.21]{PP21} using \cref{thm:isomorphism} and the tropical Hodge theory of \cite{ADH20}.
If $R$ is a $\Z$-algebra, define $R_\Q \coloneqq R \otimes_\Z \Q$ and $R_\R$ likewise.
\begin{cor}\label{cor:hodge}
  \newcommand{\mathbb{K}}{\mathbb{K}}
  Let $\mathbb{K} = \Q, \R$.
  Let $P$ be a polymatroid of rank $r$ and $\G$ a geometric building set.
  Let $\ell$ be any $\mathbb{K}$-valued strictly convex piecewise linear function on $\Sigma_{P,\G}$, viewed as an element of $A^1(\Sigma_{P,\G})_\mathbb{K} \cong \DP^1(P,\G)_\mathbb{K}$.
  \begin{enumerate}\itemsep 5pt
  \item \label{poincare} (Poincar\'e duality) There is an isomorphism
    \[ \deg: \DP(P,\G)^{r-1} \to \Z, \]
    and for all $0 \leq k < r/2$, the pairing
      \[ \DP^k(P,\G) \times \DP^{r-k-1}(P, \G) \to \Z, \quad (a,b) \mapsto \deg(ab)  \]
    is non-degenerate.
\item \label{lefschetz} (Hard Lefschetz) For every $0 \leq k < r/2$, the multiplication map
  \[ \DP^k(P,\G)_\mathbb{K} \to \DP^{r-k-1}(P,\G)_\mathbb{K}, \quad a \mapsto \ell^{r - 2k-1} a \]
  is an isomorphism.
\item \label{hodge-riemann} (Hodge-Riemann) For every $0 \leq k < r/2$, the bilinear form
  \[
    \DP^k(P,\G)_\mathbb{K} \times \DP^k(P,\G)_\mathbb{K} \to \mathbb{K}, \quad (a,b) \mapsto (-1)^k \deg(\ell^{r-2k-1} ab)
  \]
  is positive definite on the kernel of multiplication by $\ell^{r-2k}$.
  \end{enumerate}
\end{cor}
\begin{proof}
  Let $M = \widetilde P$.
  By \cite[Theorem 4]{FY04}, $\Sigma_M$ is a refinement of $\Sigma_{P,\G}$. In particular, the two fans have the same support.
  By \cite[Proposition 2.4]{AHK18} and \cref{polytope}, both fans are subfans of the normal fans of convex polytopes, so both fans support strictly convex piecewise linear functions. By \cite[Proposition 2]{FY04}, $\Sigma_{P, \mathscr{G}}$ is a smooth fan.

  The desired statements now follow by applying \cite[Theorem 1.6]{ADH20}, \cite[Proposition 5.2]{AHK18}, and \cite[Theorems 6.19 and 8.8]{AHK18}.
  Loosely, \cite[Theorem 1.6]{ADH20} says properties \cref{poincare}, \cref{lefschetz}, and \cref{hodge-riemann} hold for $A(\Sigma_{P,\G})_\R \cong \DP(P,\G)_\R$ if and only if they hold for $A(\Sigma_M)_\R$, and \cite{AHK18} verifies them for $A(\Sigma_M)_\R$.

  Properties \cref{lefschetz} and \cref{hodge-riemann} for $\DP(P, \G)_\Q$ follow immediately from those for $\DP(P, \G)_\R$.
  For \cref{poincare}, note that \cite{ADH20}'s Poincar\'e duality arguments go through over $\Z$. (Explicitly, one must check statements 6.6--6.9, and Propositions 6.16, 6.17 of \cite{ADH20}.)
\end{proof}
\begin{rmk}\label{ConeRemark}
  In \cite{PP21}, \cref{cor:hodge} is proved for $\ell$ in the \word{$\sigma$-cone} \cite[Definition 4.15]{PP21}, the positive span of
  \[
    - \sum_{G \in \G_{\geq F}} x_G, \quad F \in \G
  \]
  in $\DP(P, \G)$. The $\sigma$-cone is generally a proper subset of the cone of strictly convex piecewise linear functions on $\Sigma_{P, \G}$.
  For example, if $M$ is a loopless matroid on $E$ and $\G$ is its maximal building set, then for any $i \in E$,
  \[
    \beta \coloneqq \sum_{i \not \in F} z_F = - \sum_{|G| > 1} \big(|G| - 1\big) y_G
  \]
  is in the closure of the cone of strictly convex piecewise linear functions on $\Sigma_M$ \cite[Proposition 4.3, Lemma 9.7]{AHK18}.
  However, $\beta$ may not be in the closure of the $\sigma$-cone, e.g. when $M$ is Boolean of rank at least 3.
  For comparison of the $\sigma$-cone to the ample cone of the wonderful compactification, see \cite[Remark~4.22]{PP21}.
\end{rmk}


\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