{'fraction_non_alphanumeric': 0.0328465046470011, 'fraction_numerical': 0.013392599434278127, 'mean_word_length': 4.801741460147355, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatial Correlation in Weather Forecast Accuracy: A Functional Time Series Approach

A functional time series approach is proposed for investigating spatial correlation in daily maximum temperature forecast errors for 111 cities spread across the U.S. The modelling of spatial correlation is most fruitful for longer forecast horizons, and becomes less relevant as the forecast horizon shrinks towards zero. For 6-day-ahead forecasts, the functional approach uncovers interpretable regional spatial effects, and captures the higher variance observed in inland cities versus coastal cities, as well as the higher variance observed in mountain and midwest states. The functional approach also naturally handles missing data through modelling a continuum, and can be implemented efficiently by exploiting the sparsity induced by a B-spline basis. The temporal dependence in the data is modeled through temporal dependence in functional basis coefficients. Independent first order autoregressions with generalized autoregressive conditional heteroskedasticity [AR(1)+GARCH(1,1)] and Student-t innovations work well to capture the persistence of basis coefficients over time and the seasonal heteroskedasticity reflecting higher variance in winter. Through exploiting autocorrelation in the basis coefficients, the functional time series approach also yields a method for improving weather forecasts and uncertainty quantification. The resulting method corrects for bias in the weather forecasts, while reducing the error variance.


Introduction
A functional time series is a time-indexed sequence of stochastic processes {f t (τ )} ∞ t=1 where each f t (·) is a random function on the domain T . By unifying functional data analysis with time series analysis, it presents an approach to modelling randomness on curves, surfaces, and other phenomena varying over a spatial continuum where these functional data are observed regularly over time and exhibit serial dependence. Ramsay and Silverman (2005) and Tsay (2010) provide background for functional data analysis and time series analysis respectively, while Hörmann and Kokoszka (2012) provides background for functional time series. Hyndman and Shahid Ullah (2007) and Hyndman and Booth (2008) propose a forecasting approach for functional time series and apply it to demographic data. Hyndman and Shang (2009) proposes a weighted functional approach which assigns more weight to recent observations and yields an improvement in forecast accuracy. Further developments in functional forecasting with applications to demographic data are found in Shang et al. (2011), Hyndman et al. (2013), and Dokumentov et al. (2018). Aue et al. (2017)  The most common examples of functional data are samples of 1-dimensional curves, and a contribution of this paper is the analysis of a model for 2-dimensional surface data. In this application, the spatial domain T is a rectangle in R 2 containing the range of longitudes and latitudes covering the lower 48 states. The data set is taken from the 2018 American Statistical Association Data Expo 3 , consisting of daily maximum temperature forecasts from the National Weather Service for 111 cities spread across the US (excluding Alaska and Hawaii) over the period from July 2014 to September 2017. The locations of the cities are illustrated in Figure 1. Forecasts range from same-day to six-days-ahead and are compared to actual temperature recorded at city airports.
A functional time series approach is applied to investigate and extract the structure of the spatial correlation in forecast errors. Through modelling the entire continuum instead of individual points, the proposed functional data approach also naturally handles missing data. This is a vital benefit, as data records are frequently incomplete and forecasts are not always available at every location. This paper also extends the methodology of Hyndman and Shahid Ullah (2007) with the incorporation of a heteroskedastic time series model to account for the higher unpredictability of weather in winters compared to summers.
Through exploiting both the spatial and temporal dependence in weather forecast errors, one can estimate the next day's expected forecast error and correct the forecast accordingly.
In this paper, a four-step procedure is used to facilitate this estimation. Section 2.1 describes the initial pre-processing step, where we first construct continuous surface data from the discrete observations using manually selected B-splines. The resolution of the chosen B-splines are set so that cities are able to be grouped into regions of interest. Once these surfaces are built, the second step is to reduce dimension to a smaller spatial basis which describes the most important modes of variation, and this is described in Section 2.2. After the spatial basis is chosen, the next step is to model the temporal dependence in the random coefficients which scale the basis functions, and this is described in Section 2.3. The time series model facilitates the final step, the prediction of forecast errors, and this is described in Section 3.1.

Methodology
We first fix a forecast horizon h of interest, in this case between 0 to 6 days ahead. For the given forecast horizon, let Y t (τ ) be the forecast error of maximum temperature on day t for the city located at spatial coordinates τ ≡ (longitude, latitude). Spatial correlation is captured through the following spatio-temporal random effect model: (1) Here, µ(τ ) represents the mean forecast error and is assumed fixed over time. Estimation of the mean function is described in Section 2.1. The spatial basis functions ϕ 1 (τ ), . . . , ϕ K (τ ) describe the main modes of variation in the forecast errors, capturing spatial dependence across different regions of the U.S. The construction of the basis functions is described in Sections 2.1 and 2.2. The random coefficients β kt capture temporal correlation. They are modelled using independent AR(1)+GARCH(1,1) processes, as described in Section 2.3. Lastly, ε t (τ ) is a white noise process independent of the random coefficients, and is assumed i.i.d. N(0, σ 2 ) for all t and τ . It is used to capture, including any measurement errors, the remaining variation not explained by the spatiotemporal random effect.

Pre-Processing with Spline Smoothing
A two-step procedure is employed to construct the spatial basis functions ϕ k . As an initial preprocessing step, we first need to construct continuous surface data from the discretely observed values using manually selected B-splines, where the knot sequence is at the right level of tightness to capture regional groupings of the cities. Once these surfaces are built, the singular value decomposition (SVD) is used in Section 2.2 to reduce dimension and identify the ϕ k which describe the most important modes of variation.
For simplicity, we consider the spatial domain as a subset of R 2 . Specifically, we define a 2- For this particular dataset, knot sequences with 13 equally spaced interior knots were able to capture interesting regional groupings of the cities, and the knots are visualized against the cities by the grey grid lines in Figure [24, 24, 24, 24, 25.79, 27.57, . . . , 49, 49, 49, 49], and For each day t = 1, . . . , m, we fit splines such that the coefficientsĉ t,1 , . . . ,ĉ t,289 solve the following optimization problem: where τ t,1 , . . . , τ t,nt are the observation locations available on day t. This handles (moderate amounts of) missing data naturally since missing observations Y t (τ ) are simply omitted from the objective function. Also, since B-splines have compact support, the resulting system is sparse and can be solved efficiently.
As the number of coefficients exceeds the number of observation locations, a truncated singular value decomposition is used to solve for the coefficients. The truncated SVD acts as a form of regularization similar to ridge regression, and further details can be found in Hansen (1987). In Figure 1, there is a lack of observations in regions of the rectangle outside the borders of the United States. The truncated SVD yields stable results by imposing zero coefficients for the splines that lie in these data-free regions.
We denote the resulting coefficient matrix by C = [ĉ t,j ]. We use the column meansC = c 1 , . . . ,c 289 to estimate the mean function µ(τ ) as: is the sample mean of all observed forecast errors.

Constructing Spatial Basis Functions ϕ k
As each spline function captures variation specific to a sub-rectangle of the domain, the amount of information in each region of the domain is encoded through the coefficients of its corresponding splines. Through studying the covariance of the B-spline coefficients, we are able to piece together variation across multiple sub-rectangles to study the important regional effects in the weather forecast errors. A principal component analysis yields a ranked order of groupings through the eigendecomposition of the covariance matrix, or equivalently, the singular value decomposition of the mean-centered data matrix. The eigenvectors of the covariance matrix identify weightings of the spline coefficients which contribute the most variation, with said variation quantified by the size of the eigenvalue.
After mean-centering the columns of C, the singular value decomposition C −C = UΣV T provides principal component loadings as the columns of V , allowing for dimension reduction that is mean-square optimal in the coefficient space Col(C −C). Refer to Jolliffe (2002) for further background on principal component analysis and singular value decomposition. Assuming the SVD is written in descending order of singular values, the spatial basis functions ϕ k (τ ) are built using the first K columns of V , representing the K most important principal components: The first four basis functions are shown in Figure 2 for the 6-day-ahead forecast horizon. The first basis function represents an inland versus coastal effect, as cities further inland had greater variance in their forecast errors compared to coastal cities. The second basis function represents an east versus west effect, with the opposing signs of the regions allowing for a differentiation between the regions. The third and fourth basis functions represent mountain state and midwest state effects respectively, as these regions have the most unpredictable weather.
This structure of spatial correlation is most prevalent in 6-day-ahead forecasts, but vanishes as the forecast horizon shrinks to zero. For example, Figure 3 shows the first four basis functions for same-day forecasts, and these basis functions lack any coherent spatial structure.
For the 6-day-ahead forecasts, a cutoff of K = 20 basis functions is used. Each of the 20 selected basis functions ranged from 9% to 1.4% of explained variance, while the basis functions after the 20th each contribute less than 1.4% of explained variance and revealed no interesting spatial structure. The results of Table 2 were stable above K = 20 basis functions, justified by Figure 7.

Modelling Random Coefficients
Once the spatial basis functions are defined as above, coefficients β 1t , . . . , β Kt for each day t are estimated on this reduced K-dimensional basis as solutions to the following optimization problem: Temporal dependence in forecast errors is then modelled through these coefficients. Empirically, an AR(1)+GARCH(1,1) model with Student-t innovations was found to provide a good description of the observed coefficients. Specifically, where F t−1 is the information set up to time t − 1 (the σ-field generated by u k,t−1 , u k,t−2 , . . . for all k). The u 1t , . . . , u Kt are assumed conditionally independent given F t−1 . The proposed model implies the following conditional mean and covariance functions: The resulting innovations for the first basis function u 1t are shown in Figure 4 on the top left. Notably, the innovations exhibit a seasonal heteroskedasticity with winter weather being the most unpredictable. The GARCH process characterizes the heteroskedasticity well, as the standardized innovations u 1t /η 1t exhibit approximately constant variance, and the squared standardized innovations show no significant autocorrelation.
For the 6-day-ahead weather forecasts, the AR+GARCH parameter estimates for the first four basis functions are shown in Table 1

Empirical Performance
For each pair of cities (i, j), sample spatial correlation for the 6-day-ahead forecasts is computed before and after accounting for the spatial basis functions. More specifically, given a pair of cities located at τ i and τ j , the top of Figure 5 shows and the bottom of Figure 5 shows for the residualsε t (τ ) = Y t (τ ) −μ(τ ) − K k=1 β kt ϕ k (τ ). Above, the sums are over the days t with no missing observations, andȲ · (τ i ) andε · (τ i ) are sample averages of Y t (τ i ) andε t (τ ) respectively over such t.
The cities are numbered from 1 to 111 (shown on the xand y-axes), and cities are ordered from east to west, resulting in a concentration of high correlation along the main diagonal in the first figure. After accounting for K = 20 spatial basis functions, the second figure indicates the lack of spatial correlation in the residuals and provides evidence that the proposed model provides an adequate approximation of the observed spatial correlation in forecast errors.
The correlograms before and after accounting for the basis functions are shown in Figure 6, visualizing the sample spatial correlations against city distance. The cutoff of K = 20 is justified in Figure 7, where the squared Frobenius norm of the residual spatial correlation matrix ( i,j (ρ after i,j ) 2 ) is computed for different values of K. The Frobenius norm stops significantly decreasing after K = 20 indicating that additional basis functions do not explain any more spatial correlation, and capture individual city variation rather than regional variation.   i,j ) 2 (blue) for a changing number of basis functions K. After K = 20 basis functions, the Frobenius norm of the residual spatial correlation matrix stops decreasing, indicating the remaining basis functions capture local variation rather than regional variation.

Predicting Forecast Errors to Improve Forecast Accuracy
Through exploiting autocorrelation in the coefficients, the AR+GARCH parameter estimates can be used to predict the next day's basis coefficients, as β k,t+1 |F t ∼ t ν k (ψ k β kt , η 2 k,t+1 ), and (13) Substituting the parameter estimatesψ k ,ω k ,α k ,γ k , andν k in place of the true parameters yields an approximate distribution which can be used for prediction and uncertainty quantification.
By setting the predicted coefficient toβ k,t+1 =ψ k β kt , this can then be used to predict next day's weather forecast errors by settinĝ These predicted errorsŶ t+1 (τ ) can then be used to adjust the next day's weather forecast (of the same horizon). More explicitly, for forecast horizon h and location τ , if F t,t+h (τ ) is the forecast on day t of the maximum temperature on day t + h and A t+h (τ ) is the actual temperature on day t + h so that Y t (τ ) = F t,t+h (τ ) − A t+h (τ ), we can define an adjusted forecast for day t + 1 as The improvement on weather forecasts can then be assessed by comparing the adjusted errors to the unadjusted errors Y t+1 (τ ).
The mean and variance of Y and Z across all t and τ is shown in Table 2.
The change in mean error indicates that the proposed method eliminates bias in the weather forecasts, while the reduction of standard deviation by 12% indicates a significant average improvement in forecast accuracy. The shift in error distribution is illustrated in Figure 8, with the adjusted distribution centered at 0 and with smaller variance.

Code
For reproducing our analysis and figures, MATLAB code is available at: https://github.com/pjang23/dataexpo2018-functionaltimeseries To accomodate those without the MATLAB Econometrics Toolbox™, the MATLAB package 'ARMAX-GARCH-K-SK' by Gabrielsen (2021) is included as an alternative for fitting the AR+GARCH models.

Conclusions
We have introduced a functional time series approach to investigating spatial correlation in weather forecast accuracy. The modelling of spatial correlation is most fruitful for the longer forecast horizons, and becomes less relevant as the forecast horizon shrinks towards zero. For 6-day-ahead weather forecasts, the functional approach uncovers interpretable regional spatial effects, and captures the higher variance observed in inland cities versus coastal cities, as well as the higher variance observed in mountain and midwest states. The functional approach also naturally handles missing data and can be implemented efficiently by exploiting the sparsity induced by using a B-spline basis.
Independent first order autoregressions with generalized autoregressive conditional heteroskedasticity [AR(1)+GARCH(1,1)] and Student-t innovations worked well to capture the persistence of coefficients over time and the seasonal heteroskedasticity reflecting higher variance in winter. Autocorrelation in the basis coefficients can further be exploited to improve weather forecasts, especially at longer horizons, and the resulting approach eliminates bias while reducing error variance.