\section{Introduction}
\input{introduction}

\section{Related Work}
\input{related_work}
\label{sec:related_work}



\section{Sub-seasonal Climate Forecasting}
\label{sec:ssf_difficulty}
\input{problem_statement}


\section{Methods}
\input{method}
\label{sec:models}

\section{Data and Experimental Setup}
\input{experiment_set_up}
\label{sec:exp_setup}

\section{Experimental Results}
\input{results}
\label{sec:results}


\section{Conclusion}
\input{conclusion}
\label{sec:conclusion}


\section{Broader Impact}
\input{broader_impact}

\section*{Acknowledgement}
The research was supported by NSF grants OAC-1934634, IIS-1908104, IIS-1563950, IIS-1447566, IIS-1447574, IIS-1422557, CCF-1451986. The authors would like to acknowledge the computing support from the Minnesota Supercomputing Institute (MSI) at the University of Minnesota.

\bibliographystyle{plain}

\section{Difficulty of the Problem}

\label{sec:app_ssf_difficulty}

\subsection{Dependence between historical data and forecasting target}
In section \ref{sec:ssf_difficulty}, the dependence between the most recent historical data (the normalized average temperature of week -2 \& -1) and the forecasting target (the normalized average temperature of week 3 \& 4) is measured by maximum information coefficient (MIC). Here we show the results measured by Pearson correlation coefficient \cite{wasserman2013all}, and Spearman's rank correlation coefficient \cite{wasserman2013all} (Figure \ref{fig:corr_margin}). Small values ($\le$0.2) of Pearson correlation and Spearman's rank correlation at a majority of locations, which verify that there is little information shared between the most recent date and the forecasting target, once again, demonstrate how difficult SSF is.


\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{mi_marginal.png}
\caption{Pearson correlation, Spearman's rank correlation and MIC between 2m temperature of week -2 \& -1 and week 3 \& 4. Small values ($\le$0.2) of Pearson correlation and Spearman's rank correlation at a majority of locations verify the fact, as we illustrate in the main paper using MIC, that there is little information shared between the most recent date and the forecasting target.}
\label{fig:corr_margin}
\end{figure}






\subsection{Relative R$\mathbf{^2}$}

In the main paper, we introduce cosine similarity, which is widely used in weather prediction evaluation, as an evaluation metric. Here we formally define the other evaluation metric, namely relative $R^2$ as

\begin{equation}
\text{Relative } R^2=1-\text{Relative MSE}=
1-\frac{\sum_{i=1}^n(\mathbf{y}^*_i-\hat{\mathbf{y}}_i)^2}{\sum_{i=1}^n(\mathbf{y}^*_i-\bar{\mathbf{y}}_{\text{train}})^2}~,
\end{equation}

where $\hat{\mathbf{y}}$ denotes a vector of predicted values, and $\mathbf{y}^*$ be the corresponding ground truth. We use relative $R^2$ to evaluate the relative predictive skill of a given prediction $\hat{y}$ compared to the best constant predictor $\bar{y}_{\text{train}}$, the long-term average of target variable at each date and each target location computed from training set. A model which achieves a positive relative $R^2$ is, at least, able to predict the sign of $y^*$ accurately. The results of temporal and spatial relative $R^2$ over the US mainland of ML models discussed in section \ref{sec:models} are shown in Table \ref{tab:r2_all} and Figure \ref{fig:spatial_result_r2} respectively.


\section{Data and experimental setup}
\label{sec:app_exp_setup}
\subsection{Data sources}
The data described in Table \ref{table: data} were downloaded from the following sources:
\begin{itemize}[leftmargin=*]
    \item Temperature (tmp2m):
    \url{https://www.esrl.noaa.gov/psd/data/gridded/data.cpc.globaltemp.html}
    \item Soil moisture (sm): \url{https://www.esrl.noaa.gov/psd/data/gridded/data.cpcsoil.html}
    \item Sea surface temperature (sst): \url{https://www.ncdc.noaa.gov/oisst}
    \item Relative humidity (rhum), sea level pressure (slp), and geopotential height (hgt): \url{ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/}
    \item Multivariate ENSO index (MEI): \url{https://psl.noaa.gov/enso/mei/}
    \item Ni\~no indices: \url{https://www.cpc.ncep.noaa.gov/data/indices/wksst8110.for}
    \item North Atlantic Oscillation (NAO) index: \url{ftp://ftp.cpc.ncep.noaa.gov/cwlinks/norm.daily.nao.index.b500101.current.ascii}
    \item Madden Julian Oscillation (MJO) phase \& amplitude: \url{http://www.bom.gov.au/climate/mjo/graphics/rmm.74toRealtime.txt}
\end{itemize}


\subsection{PCA prepossessing}

As mentioned in section \ref{sec:exp_setup} of the main paper, one way for feature extraction is to apply PCA to spatial-temporal variables. To do so, let's consider sst of Pacific ocean as an example. Daily sst of Pacific ocean is originally stored in a matrix, of which each element represents the sea surface temperature at each grid point of Pacific ocean. The covariance matrix can be computed by flattening each matrix into a 1-D vector, viewing each element in the matrix as a feature and each date as one observation. Such covariance matrix captures spatial connection among grid points of Pacific ocean. By considering all dates from 1986 to 2016, we can extract the top 10 principal components (PCs) as features based on PC loadings computed from the corresponding covariance.


\begin{figure}[t]
    \centering
    \subfigure[Sequential feature set for Mar 1, 2018]{\includegraphics[width=0.45\columnwidth]{feature_new.png}}
    \hspace{5mm}
    \subfigure[Evaluation pipeline]{\includegraphics[width=0.35\columnwidth]{model_selection.png}}
    \vspace{-2mm}
    \caption{(a) Sequential feature set: to construct feature set at Mar. 1, 2018, we concatenate covariates from Mar. 1 in 2018, 2017, and 2016, their corresponding $7^{th}$, $14^{th}$, and $28^{th}$ days in the past, and $7^{th}$, $14^{th}$, and $28^{th}$ ``future'' days in 2017 and 2016. (b) Evaluation pipeline: to test SSF in Jan 2017,
    the training set covers historical 30 year ends at Dec 4, 2016 (the last available date). 5 validation sets include dates from each Jan between 2012 to 2016, with the corresponding training sets generated by applying a moving window of 10 years and a stride of 365 days on data start from 2000.}
     \label{fig:model_selection}
    
\end{figure}


 \subsection{Feature set construction}
 
To better utilize historical information, we construct a sequential feature set by including not only covariates of the target date, but also covariates of the $7^{th}$, $14^{th}$, and $28^{th}$ day previous from the target date, as well as the day of the year of the target date in the past 2 years and both the historical past and future dates around the day of the year of the target date in the past 2 years. Such selection of historical dates mainly bases on the temporal correlation. Figure \ref{fig:model_selection}(a) provides a detailed example on how to construct feature set for Mar 1, 2018: we concatenate covariates from Mar. 1 in 2018, 2017, and 2016, their corresponding $7^{th}$, $14^{th}$, and $28^{th}$ days in the past, and $7^{th}$, $14^{th}$, and $28^{th}$ ``future'' days in 2017 and 2016. In total, we include $H=18$ historical days in our feature set for each date.


\subsection{Evaluation pipeline}


Predictive models are created independently for each month in 2017 and 2018. To mimic a live forecasting system, we generate 105 test dates during 2017-2018, one for each week, and group them into 24 test sets by their month of the year. Given a test set, our evaluation pipeline consists of two parts (Figure~\ref{fig:model_selection}(b)):

\begin{itemize}[leftmargin=*]
\setlength\itemsep{0em}
\item  ``5-fold'' training-validation pairs for hyper-parameter tuning, based on a ``sliding-window'' strategy designed for time-series data. Each validation set uses the data from the same month of the year as the test set. For instance, if the test set is Jan 2017, the corresponding 5 validation sets are Jan 2012, Jan 2013, Jan 2014, Jan 2015, and Jan 2016 respectively. Each validation set corresponds to a training set containing 10 years of data and ending 28 days before the first date in the validation set. Specifically, if the validation set starting from Jan 1, 2016,
the training set is from Dec 4, 2005 to Dec 4, 2015. Such construction is equivalent to apply a sliding-window of 10-year with a stride of 365 days on data from 2002.

\item The training-test pair, where the training set, including 30-year data in the past, ends 28 days before the first date in the test set. For example, to test SSF in Jan 2017, i.e., Jan 1, Jan 8, Jan 15, Jan 22, and Jan 29, the training set starts from Dec 4, 1986 and ends at Dec 4, 2016, which is the $28^{th}$ day before Jan 1, and the last date we have the ground truth for the target variable.
\end{itemize}




\section{Additional Results}
\label{sec:app_results}


 
\begin{table}
  \caption{Comparison of relative $R^2$ of tmp2m forecasting for test sets over 2017-2018. A positive relative $R^2$ indicates a model predicting the sign of the target variable correctly. XGBoost achieves the highest relative $R^2$. 
  }

  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
  \hline
  \multicolumn{5}{c}{\textbf{Temporally Global Set}}\\
\hline
    \textbf{XGBoost - one day }&\textbf{0.0760(0.03)} &\textbf{0.0974(0.03)}  &\textbf{-0.0449(0.03)}&\textbf{ 0.2434(0.03)}\\
  Lasso - one day &  0.0552(0.02) &0.0321(0.02)& -0.0309(0.01)& 0.1295(0.02)\\
    \hline
     \textbf{Encoder (LSTM)-Decoder (FNN)}&  \textbf{-0.0353 (0.05)} &\textbf{0.0596(0.05)}&\textbf{-0.2409 (0.06)} &\textbf{ 0.2426 (0.05)}\\
   
    
     FNN & -0.5777(0.29)& -0.0183(0.15)&  -0.0794(0.13)& 0.0213(0.13)\\
    CNN & -0.0564(0.03)& 0.0284(0.02) & -0.0266(0.02)& 0.0570(0.02)\\
     CNN-LSTM &  -0.1164(0.05) &0.0263(0.03) & -0.0862(0.03)& 0.0698(0.03)\\
    \hline
\textbf{LS with NAO \& all nino  - daily }& \textbf{0.0418(0.01)} 
&\textbf{0.0535(0.01)}  &\textbf{-0.0078(0.01)}  & \textbf{0.0949(0.01)} \\
Damped persistence & 0.0266(0.01) &0.0414(0.02) & -0.0542(0.02) & 0.1354(0.02)\\
    \hline 
    MultiLLR  &-0.0571 (0.02) & 0.0034 (0.02)
&-0.1156 (0.03) & 0.0797 (0.02)\\
AutoKNN  &0.0181 (0.01) & 0.0260 (0.02)&-0.0531 (0.02) & 0.1041 (0.01)\\
    \hline
\multicolumn{5}{c}{\textbf{Temporally Local Set}}\\
\hline
XGBoost - one day &-0.0337(0.03) &0.0396(0.03)&  -0.1310(0.04) &0.1873(0.03)\\
Lasso - one day &    -0.0028(0.02)& 0.0327(0.02) & -0.0613(0.02)& 0.0996(0.02)\\
Encoder (LSTM)-Decoder (FNN)&-0.2333 (0.06) & -0.1116 (0.06)
&-0.4694 (0.09) & 0.1808 (0.06)\\
\hline
\end{tabular}}
\label{tab:r2_all}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{spatial_result_r2.png}
\caption{Temporal relative $R^2$ over the US mainland of ML models discussed in section \ref{sec:models} for temperature prediction over 2017-2018. Large positive values (green) closer to 1 indicates better predictive skills.}
\label{fig:spatial_result_r2}
\end{figure}



\subsection{Temporal and spatial results of relative R$\mathbf{^2}$}

Table \ref{tab:r2_all} lists the mean, the median, the 0.25 quantile, the 0.75 quantile, and their corresponding standard errors of relative $R^2$ for all models.  A positive relative $R^2$ indicates a model can at least predict the sign of the target variable correctly. Again, XGBoost achieves the highest predictive skill in terms of both the mean and the median, demonstrating its predictive power. Linear regression, like Lasso, with a proper feature set has good predictive performance. Both XGBoost and Lasso have larger positive relative $R^2$ in terms of the mean, and can still outperform climatology and two climate baseline models, i.e., LS with NAO \& Ni\~no, and damped persistence. Even though Encoder (LSTM)-Decoder (FNN) has a slightly negative mean relative $R^2$, it has the second largest median and 0.75 quantile among all models, showing its potential for further improvement. 

Figure \ref{fig:spatial_result_r2} shows the spatial relative $R^2$ of all methods. XGBoost and Lasso are able to achieve positive relative $R^2$ for most of the target locations. Encoder (LSTM)-Decoder (FNN) shows better predictive skill over the southern US compared to other regions. MultiLLR and AutoKNN manages to obtain non-negative relative $R^2$ for the coastal area in the western US but their predictive performance drops in the rest of locations. All other baseline methods struggle to reach positive relative $R^2$ for most of the target locations. 



\subsection{Analysis on feature importance}

To emphasis the importance of the land-based covariates, e.g., soil moisture and the ocean-based covariates, e.g., NAO and Ni\~no indices, we compare the prediction performance among (1) the model trained with all covariates, (2) the model trained without soil moisture, and (3) the model trained without NAO and Ni\~no indices (Table~\ref{tab:cos_without} and Table~\ref{tab:r2_without}). Most models experience a performance deterioration when we exclude certain ``important'' covariates.




\subsection{The influence of feature sequence length} We compare the prediction performance under 3 different settings, referred to as ``one day'', ``four days", and ``all days'' respectively. For feature set construction, ``one day'' includes covariates at the target date only, ``four days'' also covers the  $7^{th}$, $14^{th}$, and $28^{th}$ days previous to the target date, and ``all days'' uses the exact feature sequence we use for LSTM-based models. Comparison of predictive skills under each setting, measured by both cosine similarity and relative $R^2$, can be found in Table~\ref{tab:cos_seq_len} and Table~\ref{tab:r2_seq_len}. Both XGBoost and Lasso enjoy a performance boost using ``one day'' values. Especially for XGBoost, the performance of ``one day'' is approximately 50\% better than using ``all days''. A possible explanation for such performance degradation as we increase the feature sequence length is that both models weight covariates from different dates exactly the same without considering temporal information, thus more noise has been introduced. 

\begin{table}[t]
  \caption{Comparison of cosine similarity of tmp2m forecasting for test sets over 2017-2018 using different feature set. Excluding soil moisture or climate indices (NAO \& Ni\~no) leads to a deterioration in the predictive performance.
  }

  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
\hline
\textbf{XGBoost - one day}&\textbf{ 0.3044(0.03)}& \textbf{0.3447(0.05)}&  \textbf{0.0252(0.05)}&\textbf{ 0.5905(0.04})\\
XGBoost - one day (w/o soil moisture)&0.2685(0.03)& 0.2797(0.05) & 0.0703(0.04)& 0.5492(0.05)\\
XGBoost - one day (w/o nao \& all nino) & 0.2081(0.03)& 0.1640(0.05) & -0.0588(0.04) &0.5246(0.05)\\
\hline
 Lasso - one day &  0.2499(0.04)& 0.2554(0.06) & -0.0224(0.05) &0.5604(0.06)\\
 Lasso - one day (w/o soil moisture)&0.2638(0.04)& 0.2912(0.05) & 0.0032(0.06)& 0.5655(0.05)\\
 Lasso - one day (w/o nao \& all nino)& 0.1956(0.04)& 0.2573(0.07) & -0.1657(0.06) &0.5533(0.05)\\
\hline
Encoder (LSTM)-Decoder (FNN) & 0.2616 (0.04) & 0.2995 (0.07)
&{-0.0719 (0.06)} & {0.6310 (0.05)}\\
Encoder (LSTM)-Decoder (FNN)(w/o soil moisture)&0.2157 (0.04) & 0.2909 (0.07)
&-0.1106 (0.07) & 0.5443 (0.07)\\

Encoder (LSTM)-Decoder (FNN)(w/o nao \& all nino)
&0.2236 (0.04) & 0.2395 (0.06)
&-0.1527 (0.07) & 0.5989 (0.06)\\
\hline 
\end{tabular}}
\label{tab:cos_without}
\end{table}


\begin{table}[th]
  \caption{Comparison of relative $R^2$ of tmp2m forecasting for test sets over 2017-2018. Excluding soil moisture or climate indices (NAO \& Ni\~no) leads to a smaller or even negative relative $R^2$, showing that it becomes harder for the model to predict the sign of the target variable correctly.}
 
  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
\hline
\textbf{XGBoost - one day }&\textbf{0.0760(0.03)} &\textbf{0.0974(0.03)}  &\textbf{-0.0449(0.03)}&\textbf{ 0.2434(0.03)}\\
XGBoost - one day (w/o soil moisture)&0.0370(0.03)& 0.0322(0.03) & -0.0564(0.03) &0.2225(0.03)\\
XGBoost - one day (w/o nao \& all nino)&-0.0161(0.03)& -0.0079(0.04) & -0.1618(0.03)& 0.2426(0.04)\\
\hline
Lasso - one day &  0.0552(0.02) &0.0321(0.02)& -0.0309(0.01)& 0.1295(0.02)\\
Lasso - one day (w/o soil moisture)&-0.0161(0.03)& -0.0079(0.04) & -0.1618(0.03) &0.2426(0.04)\\
Lasso - one day (w/o nao \& all nino)&0.0003(0.02)& 0.0457(0.02) & -0.1113(0.03)& 0.1641(0.02)\\
\hline
Encoder (LSTM)-Decoder (FNN) &  -0.0353 (0.05) &0.0596(0.05)&-0.2409 (0.06) &0.2426 (0.05)\\
Encoder (LSTM)-Decoder (FNN)(w/o soil moisture)
&-0.1083 (0.05) & 0.0314 (0.05)
&-0.3022 (0.08) & 0.2252 (0.05)\\

Encoder (LSTM)-Decoder (FNN)(w/o nao \& all nino)
&-0.0802 (0.04) & 0.0124 (0.05)
&-0.3032 (0.06) & 0.2446 (0.05)\\
\hline 
\end{tabular}}
\label{tab:r2_without}
\end{table}




\begin{table}[t]
   \caption{Comparison of spatial cosine similarity for tmp2m forecasting over 2017-2018 using various length of feature sequence. Including longer historical sequence leads to a deterioration in the predictive performance of XGBoost and Lasso.}
 
   \centering
   \resizebox{0.8\textwidth}{!}{
   \begin{tabular}{c|c c c c}
   \hline
   Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
     \hline
     XGBoost - all days  &0.2080(0.03) &0.1582(0.05) & -0.0466(0.05)& 0.5383(0.05)\\
     XGBoost - four days  &0.2433(0.03)&0.2203(0.05) & 0.0561(0.04)& 0.5168(0.06)\\
    \textbf{XGBoost - one day }&\textbf{ 0.3044(0.03)}& \textbf{0.3447(0.05)}&  \textbf{0.0252(0.05)}&\textbf{ 0.5905(0.04})\\
     \hline
     Lasso - all days  &0.2160(0.04)&0.2258(0.07)&-0.1381(0.06)&0.5384(0.06)\\
     Lasso - four days &0.2247(0.04)& 0.1952(0.07)&0.0572(0.06)&-0.5700(0.06)\\
    \textbf{ Lasso - one day }&  \textbf{0.2499(0.04)}& \textbf{0.2554(0.06)} & \textbf{-0.0224(0.05}) &\textbf{0.5604(0.06)}\\
     \hline
     \end{tabular}
     }
     \label{tab:cos_seq_len}
\end{table}



\begin{table}[t]
  \caption{Comparison of relative $R^2$ (with training set mean) for tmp2m prediction for test set over 2017-2019 using different length of feature sequence. Including longer historical sequence leads to a smaller or even negative relative $R^2$ for both XGBoost and Lasso.}

   \centering
   \resizebox{0.8\textwidth}{!}{
   \begin{tabular}{c|c c c c}
   \hline
  Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
     \hline
     XGBoost - all days  &-0.0200(0.03)& -0.0010(0.04)&  -0.1499(0.04)& 0.2304(0.04)\\
     XGBoost - four days  &0.0242(0.03)& 0.0193(0.03)&  -0.0786(0.03)& 0.1882(0.04)\\
     \textbf{XGBoost - one day }&\textbf{0.0760(0.03)} &\textbf{0.0974(0.03)}  &\textbf{-0.0449(0.03)}&\textbf{ 0.2434(0.03)}\\
     \hline
     Lasso - all days  &-0.0167(0.03)& 0.0367(0.03)&  -0.0639(0.02)& 0.1588(0.03)\\
     Lasso - four days &0.0518(0.02) &0.0266(0.02) & -0.0542(0.02)& 0.1653(0.03)\\
     \textbf{Lasso - one day} &  \textbf{0.0552(0.02}) &\textbf{0.0321(0.02)}&\textbf{ -0.0309(0.01)}& \textbf{0.1295(0.02)}\\
     \hline
     \end{tabular}
     }
     \label{tab:r2_seq_len}
\end{table}


\subsection{Discussion on deep learning models}

{\bf Results of DL models.} Table \ref{tab:cos_deep} and Table \ref{tab:r2_deep} compare the predictive skills of 5 DL models discussed in section~\ref{sec:results}, measured by both cosine similarity and relative $R^2$. Significant improvements can been observed as we evolve from the standard Encoder (LSTM)-Decoder (LSTM), to Encoder (LSTM)-Decoder (FNN)-last step, where ``last step'' indicates that FNN Decoder only uses the last step of the output sequence from LSTM Encoder, and finally to Encoder (LSTM)-Decoder (FNN) with FNN Decoder uses every step of the output sequence from LSTM Encoder.

One issue with Encoder (LSTM)-Decoder (FNN) is that the input features are shared by all target locations, which requires the model to identify the useful information for each locations without any help from the input.


{\bf Autoregressive (AR) component.}  Currently, the Encoder(LSTM)-Decoder(FNN) clearly considers climate covariates on a global scale, which are shared by all target locations. Nevertheless, SSF depends on not only global climate condition but also local weather change. Therefore, we seek a way to improve the model by adding an autoregressive (AR) component to capture the ``local'' information from historical data. We consider two variants of Encoder (LSTM)-Decoder (FNN). The first variant contains an AR component with the input as historical temperature at each target location, denoted as Encoder (LSTM)-Decoder (FNN)+AR. The second one includes both historical temperature and historical temporal climate variables, i.e., climate indices, as input features, denoted as Encoder (LSTM)-Decoder (FNN)+AR (CI). For both models, the final forecast is computed as a linear combination of the prediction from Encoder (LSTM)-Decoder (FNN) and the prediction from AR component for each location. Unexpectedly, as shown in Table \ref{tab:cos_deep} and Table \ref{tab:r2_deep}, simply adding the AR component to our Encoder(LSTM)-Decoder(FNN) does not help the model to perform better. However, we believe there is a better way to involve local information, and such modification is a promising direction that worth investigation in the future.


\begin{table}[t]
  \caption{Comparison of cosine similarity of tmp2m forecasting for test sets over 2017-2018 using different deep learning architectures. 
  }
  
  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
\hline
Encoder (LSTM)-Decoder (LSTM) & 0.0740(0.03) &0.0358(0.04) & -0.1569(0.03) &0.2584(0.04)\\
Encoder (LSTM)-Decoder (FNN)-last step &0.1614 (0.05)
& 0.2061 (0.08) &-0.2590 (0.08) & 0.5720 (0.08)\\
\textbf{Encoder (LSTM)-Decoder (FNN)} &\textbf{ 0.2616 (0.04}) & \textbf{0.2995 (0.07)}
&\textbf{-0.0719 (0.06)} & \textbf{0.6310 (0.05)}\\
Encoder (LSTM)-Decoder (FNN)+AR &0.1733 (0.04) & 0.1922 (0.06) &-0.0863 (0.07) & 0.5225 (0.06)\\
Encoder (LSTM)-Decoder (FNN)+AR (CI) & 0.1852 (0.04) & 0.1986 (0.05) &-0.0838 (0.06) & 0.5164 (0.05)\\
\hline
\end{tabular}}
\label{tab:cos_deep}
\end{table}


\begin{table}[t]
  \caption{Comparison of relative $R^2$ of tmp2m forecasting for test sets over 2017-2018. A positive relative $R^2$ indicates a model predicting the sign of the target variable correctly.  
  }

  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
\hline
Encoder (LSTM)-Decoder (LSTM) & -0.3947(0.05)& -0.2999(0.05) & -0.6606(0.08) &-0.0537(0.05)\\
Encoder (LSTM)-Decoder (FNN)-last step &-0.1709 (0.06) & 0.0217 (0.06)
&-0.4569 (0.11) & 0.2278 (0.06)\\
\textbf{Encoder (LSTM)-Decoder (FNN)}&  \textbf{-0.0353 (0.05)} &\textbf{0.0596(0.05)}&\textbf{-0.2409 (0.06)} &\textbf{ 0.2426 (0.05)}\\
Encoder (LSTM)-Decoder (FNN)+AR &-0.0414 (0.04) & -0.0041 (0.05) &-0.3027 (0.07) & 0.2309 (0.05))\\
Encoder (LSTM)-Decoder (FNN)+AR (CI) &-0.0563 (0.03) & -0.0380 (0.05) &-0.2365 (0.05) & 0.1951 (0.04)\\
\hline
\end{tabular}}
\label{tab:r2_deep}
\end{table}








\subsection{Results of all methods} 
\vspace{-2mm}
\textbf{Temporal results.} Table \ref{tab:cos_all} lists the mean, the median, the 0.25 quantile, the 0.75 quantile, and their corresponding standard errors of spatial cosine similarity of all methods. Additional results based on relative $R^2$ can be found in Appendix \ref{sec:app_results}. XGBoost, Encoder (LSTM)+Decoder (FNN) and Lasso accomplish higher predictive skills than other presented methods, and can outperform climatology and two climate baseline models, i.e., LS with NAO \& Ni\~no, and damped persistence. Overall, XGBoost achieves the highest predictive skill in terms of both the mean and the median, demonstrating its predictive power. Surprisingly, linear regression with a proper feature set has good predictive performance. Even though DL models are not the obvious winner, with careful architectural selections, they still show promising results. 



\begin{table}[t]
  \caption{Comparison of spatial cosine similarity of tmp2m forecasting for test sets over 2017-2018. XGBoost and Encoder (LSTM)-Decoder (FNN) have the best performance. Models achieve better performance using temporally global set compared to temporally local set. 
 }
  \resizebox{\textwidth}{!}{
  \begin{tabular}{c|c c c c}
  \hline
    Model & Mean(se) & Median (se) & 0.25 quantile (se) & 0.75 quantile (se)\\
    \hline
 \multicolumn{5}{c}{\textbf{Temporally Global Dataset}}\\
    \hline
    \textbf{XGBoost - one day}&\textbf{ 0.3044(0.03)}& \textbf{0.3447(0.05)}&  \textbf{0.0252(0.05)}&\textbf{ 0.5905(0.04})\\
   
  Lasso - one day &  0.2499(0.04)& 0.2554(0.06) & -0.0224(0.05) &0.5604(0.06)\\
    \hline
     \textbf{Encoder (LSTM)-Decoder (FNN)} &\textbf{ 0.2616 (0.04}) & \textbf{0.2995 (0.07)}
&\textbf{-0.0719 (0.06)} & \textbf{0.6310 (0.05)}\\

 
    
     FNN & 0.0792(0.01) & 0.0920(0.02) & 0.0085(0.02) &0.1655(0.02) \\
     CNN &0.1688(0.04) &0.2324(0.06)& -0.0662(0.06)& 0.4768(0.04)\\
     CNN-LSTM &  0.1743(0.04) &0.2867(0.06) & -0.1225(0.07) &0.5148(0.04)\\
  
    \hline
\textbf{LS with NAO \& Ni\~no }&\textbf{ 0.2415(0.03)} & \textbf{0.3169(0.04)}& \textbf{0.0454(0.05)} &\textbf{0.4624(0.03)}\\
 Damped persistence &0.2009(0.04)  &0.2310(0.06)  & -0.0884(0.06) & 0.5335(0.05))\\
    \hline 
    MultiLLR &0.0684 (0.03) & 0.1046 (0.05)
&-0.1764 (0.06) & 0.3156 (0.04)\\
AutoKNN &0.1457 (0.03) & 0.1744 (0.05)
&-0.1018 (0.06) & 0.4000 (0.04)\\
    \hline
\multicolumn{5}{c}{\textbf{Temporally Local Dataset}}\\
\hline
XGBoost - one day &0.1965(0.04)& 0.2345(0.05)&  -0.0636(0.06)& 0.5178(0.05)\\
 Lasso - one day & 0.1631(0.04) &0.2087(0.06) & -0.1178(0.05)& 0.5059(0.05)\\ 
Encoder (LSTM)-Decoder (FNN) & 0.1277 (0.04) & 0.1272 (0.06)
&-0.1558 (0.06) & 0.4971 (0.06)\\
 \hline
\end{tabular}}
\label{tab:cos_all}
\end{table}



\textbf{Spatial results.} Figure~\ref{fig:spatial_result_all} shows the temporal cosine similarity of all methods evaluated on test sets described in Section \ref{sec:exp_setup}. Among all methods, XGBoost and the Encoder (LSTM)-Decoder (FNN) achieve the overall best performance, regarding the number of locations with positive temporal cosine similarity. Qualitatively, coastal and south regions, in general, are easier to predict compared to inland regions (e.g., Midwest). Such a phenomenon might be explained by the influence of the slow-moving component, i.e., Pacific and Atlantic Ocean. Such component exhibits inertia or memory, in which anomalous condition can take relatively long period of time to decay. However, each model has its own favorable and disadvantageous regions. For example, XGBoost and Lasso do poorly in Montana, Wyoming, and Idaho, while Encoder (LSTM)-Decoder (FNN) performs much better on those regions. The observations naturally imply that the ensemble of multiple models is a promising future direction.




\textbf{Comparison with the state-of-the-art methods.} MultiLLR and AutoKNN are two state-of-the-art methods designed for SSF on western US \cite{hwang2019improving}. Both methods have shown good forecasting performance on the original target region. However, over the inland region (Midwest), Northeast, and South region, the methods do not perform so well (Figure~\ref{fig:spatial_result_all}). To be fair, even though a similar set of climate variables have been used in our work compared to the original paper \cite{hwang2019improving}, how we prepossess the data and construct the feature set are slightly different. Such differences may lead to relatively poor performance for these two methods, especially for MultiLLR.



\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{spatial_result_all.png}
\vspace{-5mm}
\caption{Temporal cosine similarity over the US mainland of ML models discussed in Section \ref{sec:models} for temperature prediction over 2017-2018. Large positive values (green) closer to 1 indicates better predictive skills. Overall, XGBoost and Encoder (LSTM)-Decoder (FNN) perform the best. Qualitatively, coastal and south regions are easier to predict than inland regions (e.g., Midwest).}
\label{fig:spatial_result_all}
\end{figure}

\subsection{Analysis and exploration}

We analyze and explore several important aspects that could influence the performance of ML models.


\textbf{Temporally ``local'' vs. ``global'' dataset.} Our training set consists of
all calendar months over the past 30 years, which we refer to as the temporally ``global'' dataset. Another way to construct the training set is to consider calendar months only within the temporal neighborhood of the test date. For instance, to build a predictive model to forecast June in 2017, the training set can only contains dates in June (from earlier years), and months that are close to June, e.g., April, May, July, and August, over the past 30 years. Such a construction account for the seasonal dependence of predictive relations, for example summer predictions are not trained with winter data. We name such dataset as a temporally ``local'' dataset. A comparison between the ``global'' and ``local'' datasets has been listed in Table \ref{tab:cos_all} where a significant drop in cosine similarity can be noticed when using ``local'' dataset for all of our best predictors, including XGBoost, Lasso, and Encoder (LSTM)-Decoder (FNN). We suspect such performance drop from ``global'' to ``local'' dataset may come from the reduction in the number of effective samples.

\begin{figure}[t]
    \centering
    \subfigure[XGBoost]{\includegraphics[width=0.45\columnwidth]{feature_importance_xgboost_daily_nao_nino_updated.png}}
    \hspace{3mm}
    \subfigure[Lasso]{\includegraphics[width=0.43\columnwidth]{feature_importance_lasso_daily_nao_nino_updated.png}}
    \vspace{-2mm}
    \caption{Feature importance scores computed from (a) XGBoost and (b) Lasso. Darker color means a covariate is of the higher importance. The first 8 rows contains the top 10 principal components (PCs) extracted from 8 spatial-temporal covariates respectively, and the last row includes all the temporal indices. Land component, e.g., soil moisture ($3^{rd}$ row from the top) and ocean components, e.g., sst (Pacific and Atlantic) and some climate indices are the most commonly selected covariates.}
     \label{fig:feature_importance}
    
\end{figure}


\textbf{Feature importance.} At sub-seasonal time scales,
it is widely believed \cite{predictability_s2s,delsole2017predictability} that land and ocean are the important sources of predictability, while the impact of atmosphere is limited (Figure~\ref{fig:corr_tmp2m} (a)). We study which covariate(s) are important, considered by ML models, based on the feature importance score. In particular, we compute the feature importance score from 2 ML models, XGBoost and Lasso (Figure~\ref{fig:feature_importance}). For XGBoost, the importance score is computed using the average information gain across all tree nodes a feature/covariate splits, while for Lasso, we simply count the non-zero coefficients of each model. The reported feature importance score is the average over 24 models (one per month in 2017-2018). In addition, we also provide measurement of feature importance based on Shapley value  (Figure~\ref{fig:shapley_values}), a concept from game theory~\cite{lipovetsky2001shapley}. To determine the Shapley value of a given feature, we compute the prediction difference between a model trained with and without that feature. Since the effect of suppressing a feature also depends on other features, we have to consider all possible subsets of other features, and compute the Shapley values as a weighted average of all possible differences. As shown in Figure \ref{fig:feature_importance} and \ref{fig:shapley_values}, among all covariates, soil moisture ($3^{rd}$ row from the top) is the variable that has been constantly selected by both models. Another set of important covariates is the family of Ni\~no indices. A LS model using those indices alone as predictors performs fairly well (Table \ref{tab:cos_all}). Sea surface temperatures (of Pacific and Atlantic) are also commonly selected. Such observations indicate that ML models pick up ocean-based covariates, some land-based covariates, and almost entirely ignore the atmosphere-related covariates, which are well aligned with domain knowledge ~\cite{predictability_s2s,delsole2017predictability}.




\begin{figure}[h]
    \centering
    \subfigure[XGBoost]{\includegraphics[width=0.46\columnwidth]{shap_xgboost_daily_nao_nino_updated.png}}
    \hspace{3mm}
    \subfigure[Lasso]{\includegraphics[width=0.45\columnwidth]{shap_lasso_daily_nao_nino_updated.png}}
    \vspace{-2mm}
    \caption{Shapley values computed from (a) XGBoost and (b) Lasso. Darker color means a covariate is of the higher importance. The observation that the most commonly selected covariates are from land and ocean are consistent with Figure \ref{fig:feature_importance}.}
    \label{fig:shapley_values}
    
\end{figure}


\textbf{The influence of feature sequence length.} To adapt the usage of LSTM, we construct a sequential feature set, which consists not only the target date, but also 17 other dates preceding the target date. However, other ML models, e.g., XGBoost and Lasso, which are not designed to handle sequential data, experience a drastic performance drop when we include more information from the past. More precisely, by including covariates from the full historical sequence, the performance of XGBoost drops approximately 50\% compared to the XGBoost model using covariates from the most recent date only. A possible explanation for such performance degradation, as we increase the feature sequence length, is that both models weight covariates from different dates exactly the same without considering temporal information, thus some irrelevant historical information might mislead the model. In Appendix~\ref{sec:app_results}, we provide a more detailed comparison among results obtained from various sequence lengths.


\subsection{What happened with DL models?}


As discussed in Section \ref{sec:ssf_difficulty}, applying black-box DL models naively does not work well for SSF. The improvement (Table \ref{tab:cos_all}), as we evolve from FNN to CNN-LSTM, and finally to Encoder (LSTM)-Decoder (FNN), demonstrates how the network architecture plays an important role, and leaves us plenty of space for further advancement. Below we focus on the discussion of feature representation, and the architecture design for sequence modeling. More discussion is included in Appendix~\ref{sec:app_results}.


\textbf{Feature representation: CNN vs. PCA.} Since SSF can be considered as a spatial-temporal prediction problem, CNN~\cite{lecun1998gradient} is a natural choice to handle the spatial aspect of each climate covariate by viewing it as a map, and can be applied as a ``supervised'' way for learning feature representation. However, our results imply that models with CNN has limited predictive skills regarding both spatial and temporal cosine similarity. CNN, while doing convolution using a small kernel, mainly focus on spatially localized regions. However, the strong spatial correlation of climate variables restricts the effectiveness of CNN kernels on feature extraction. Meanwhile, PCA, termed Empirical Orthogonal Functions (EOF)~\cite{lorenz1956empirical,von2001statistical} in climate science, is a commonly used ``unsupervised'' feature representation method, which focuses on low-rank modeling of spatial covariance structure revealing spatial connection. Our results (Table \ref{tab:cos_all}) illustrate that PCA-based models have higher predictive skills than CNN-based models, verifying that PCA is a better technique for feature extraction in SSF.

\textbf{Sequential modeling: Encoder-Decoder.} With features extracted by PCA,  we formulate SSF as a sequential modeling problem \cite{sutskever2014sequence}, where the input is the covariates sequence described in Section \ref{sec:exp_setup}, and the output is the target variable. Due to its immense success in sequential modeling \cite{srivastava2015unsupervised,venugopalan2015sequence}, the standard Encoder-Decoder, where both Encoder and Decoder are LSTM \cite{LSTM}, is the first architecture we investigate. Unfortunately, the model does not perform well and suffers from over-fitting, possibly caused by overly complex architecture. To reduce the model complexity, we replace the LSTM Decoder with an FNN Decoder which takes only the last step of the output sequence from the Encoder. Such change leads to an immediate boost of predictive performance. However, the input of the FNN Decoder mainly contains information encoded from the latest day in the input sequence and can only embed limited amount of historical information owing to the recurrent architecture of LSTM.  To further improve the performance, we adjust the connection between Encoder and Decoder, such that FNN Decoder takes every step of the output sequence from LSTM Encoder, which makes a better use of historical information. Eventually, this architecture achieves the best performance among all investigated Encoder-Decoder variants (see a detailed comparisons in Appendix~\ref{sec:app_results}). 

