The diminishing extent of Arctic sea ice is a key indicator of climate change as well as an accelerant for future global warming. Since 1978, Arctic sea ice has been measured using satellite-based microwave sensing; however, different measures of Arctic sea ice extent have been made available based on differing algorithmic transformations of the raw satellite data. We propose and estimate a dynamic factor model that combines four of these measures in an optimal way that accounts for their differing volatility and cross-correlations. We then use the Kalman smoother to extract an optimal combined measure of Arctic sea ice extent. It turns out that almost all weight is put on the NSIDC Sea Ice Index, confirming and enhancing confidence in the Sea Ice Index and the NASA Team algorithm on which it is based. Climate change is among the most pressing issues of our time, with many severe economic, environmental, and geopolitical consequences. Recently, the application of time series analytical methods to this topic -and, more broadly, a "climate econometrics" -has emerged as a vibrant research literature, as highlighted, for example, in Hillebrand et al. (2020) and the references therein. One important issue that these methods can address is the loss of Arctic sea ice. The loss of Arctic sea ice is a vital focus point of climate study. It is both an ongoing conspicuous effect of climate change and a cause of additional climate change via feedback loops. In particular, reduced Arctic sea ice boosts solar energy absorption via decreased albedo due to darkening color (e.g., Stroeve et al. (2012), Pistone et al. (2019), Diebold and Rudebusch (2019)) and increased methane release due to melting permafrost (e.g., Vaks et al. (2020)). 1 There are, however, several alternative measures of Arctic sea ice extent based on different processing methodologies of the underlying satellite-based microwave measurement data, and the choice among these measures is not clear-cut (Bunzel et al. (2016)). In this paper, we study four such sea ice extent (SIE) measures, which we denote as Sea Ice Index (SIE S ), Goddard Bootstrap (SIE G ), JAXA (SIE J ), and Bremen (SIE B ). The top panel of Figure 1 provides time series plots of these four measures of Arctic SIE for the satellite measurement era, which started in 1978. The four measures appear almost identical, because their scale is dominated by large seasonal swings. However, the effects of seasonality can be removed by plotting each month separately for the four series, as done in the lower twelve panels of All of this suggests treating the various measures as noisy indicators of latent true sea ice extent, which in turn suggests the possibility of blending them into a single combined indicator with less measurement error. Indeed some prominent studies have used simple equally-weighted averages of competing indicators, with precisely that goal. For example, a recent report on the state of the cryosphere (IPCC (2019)) uses a simple average of three indicators. 2 Simple averages, however, are often sub-optimal. Optimality generally requires use of weighted averages giving, for example, less weight to noisier indicators. Motivated by these considerations, in this paper, we propose and explore a dynamic factor state-space model that combines the various published indicators into an optimal measure of sea ice extent, which we extract using the Kalman smoother. We proceed as follows. In section 2, we describe the four leading Arctic sea ice extent indicators that we study, and the satellites, sensors, and algorithms used to produce them. In section 3, we propose a basic dynamic-factor state-space model for sea ice extent and use it to obtain optimal extractions of latent extent. We conclude in section 4. Sea ice extent (SIE) indicators are constructed from satellite measurements of the earth's surface using passive microwave sensing, which is unaffected by cloud cover or a lack of sunlight. Several steps are necessary to convert raw reflectivity observations into final SIE measurements. First, for a polar region divided into a grid of individual cells, various sensors record a brightness reading or "brightness temperature" for each cell. An algorithm then transforms these brightness readings into fractional surface coverage estimates -sea ice concentration (SIC) values -for each grid cell. Finally, SIE is calculated by summing the area of all cells with at least 15 percent ice surface coverage. 3 This up-rounding in SIE is effectively a bias correction, as determining the edge between ice and water can be especially difficult in the summer, when, for example, melting pools on summer ice surfaces can be mistaken for ice-free open water (Meier and Stewart (2019)). Different algorithms for processing the raw measurements importantly shape the final SIE estimates. In addition, the SIE series are not based on identical raw data because they use somewhat different satellites and sensors (Comiso et al. (2017), Comiso (2007)). In this section, we review some aspects of the satellites, sensors, and algorithms that underlie the SIE measures. Once brightness data have been recorded by satellite sensors, an algorithm converts the measurements into estimates of SIC. Here, we discuss the algorithms and other details of the various SIE indicators. Updated on a daily basis and distributed by the National Snow and Ice Data Center (NSIDC), the Sea Ice Index (SII or SIE S ) combines two separate Sea Ice indicators: (1) the Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data (NASA Team) (Cavalieri et al. (1996)) -produced at the -and (2) the Near-Real-Time DMSP SSMIS (NRTSI) (Maslanik and Stroeve (1999)) -produced by the NSIDC itself. 6 A time-lag of about one year between the SIE estimates by NASA Team and its publication in the NSIDC database requires the NRTSI to complement the SII. The NRTSI follows the NASA Team algorithm as closely as possible, but inconsistencies between the two series cannot be ruled out entirely (Fetterer et al. (2017)). In particular, the two sub-indicators use brightness temperatures from different providers.
7 These raw readings can be distorted by weather effects, making open water look like sea ice cover. Therefore, post-calculation quality checks apply land and ocean masks, to remove errorneous and implausible ice covers. However, NASA Team and NRTSI do not apply the exact same filters (Fetterer et al. (2017)). The former algorithm additionally screens the data manually for falsely detected ice formation (Cavalieri et al. (1996)), which can enhance accuracy but also reduce transparency of the final measurements. As Table 1 shows, the SII obtains the raw data from different generations of satellites and sensors. To make the data comparable, a linear least-squares model on the brightness temperatures, as reported by the two distinct sensors for an overlapping period of operation, is intended to adjust the reference points of 100 percent sea ice and 100 percent open water. These tie points then remain fixed over the lifetime of the new system (Cavalieri et al. (2012)). Another sea ice indicator, distributed by the NSIDC, relies on SIC estimates from the Goddard Bootstrap algorithm (SIE G ). 8 Despite the NASA Team and the Goddard Bootstrap algorithms having both been developed at the NASA there are some differences between the two approaches. These arise mostly from the calibration of tie points: While the NASA Team adjusts these reference points for 100 percent open water and 100 percent ice only when a new satellite or sensor becomes operational, the Goddard Bootstrap algorithm adjusts these reference points on a daily basis to account for varying weather conditions (Comiso et al. (2017)). Differing weather filters and sensitivities to varying physical temperature also lead to differences in the final measurements (Comiso et al. (1997)). In contrast to the NASA Team, the strength of the Goddard Bootstrap algorithm is the identification of melting sea ice. Therefore, the Goddard Bootstrap algorithm provides more accurate estimates of the edge of the ice cover (Goldstein et al. (2018)). Although differences between SIE S and SIE G are generally assessed to be small, they cannot necessarily be neglected (Goldstein et al. (2018)). 9 The differences between the NASA Team and Goddard Bootstrap algorithms occur especially during the melting season, when the former generally reports larger deviations from ship or radar observations. However, the relative accuracy of the two algorithms is not clear-cut, as the Goddard Bootstrap algorithm is highly sensitive to physical temperature and underestimates SIC during winter times in the higher latitudes of the Arctic region. As listed in Table 1, both SIE S and SIE G rely on the same set of instruments, which have been criticized for their low spatial resolution (Goldstein et al. (2018)). (JAXA) sea ice measure (SIE J ) uses an adapted version of the Goddard Bootstrap algorithm to derive SIC measures from satellite readings with higher spatial resolution (Comiso et al. (2017b)). 10 However, readings from these high resolution satellites are only available since 2000, so their data must merged with observations from older sensors to extend the data coverage to 1978 (Comiso et al. (2017)). SIE J also distinguishes itself from SIE S and SIE G by using a 5-day moving average of observations to compensate for potentially missing data. Using observations delivered by the high-resolution AMSR-E sensor, a group of researchers at the University of Bremen developed the ARTIST Sea Ice (ASI) algorithm (Spreen et al. (2008)) to estimate daily SIC. 11 The time-series (SIE B ) uses different algorithms for different sensors. Until the launch of the AMSR-E sensor in 2003, SIE B used the NASA Team algorithm to transform 9 See also brightness readings into SIC values. From then on, the NASA Team algorithm is replaced by the ASI. 12 3 Optimal Extraction of Latent Extent The four sea ice indicators discussed above differ in terms of the raw data sources and the algorithms used to process the raw data. They can be viewed as distinct indicators of an unobserved or latent "true" sea ice extent, SIE * . Blending such a set of noisy indicators can produce a single series with less measurement error. Here we formalize this intuition in a state-space dynamic-factor model, from which we extract an optimal composite estimate of SIE * from the four component indicators. 13 We work in a state-space environment, modeling each of the four indicators (SIE S , SIE J , SIE B , and SIE G ) as driven by latent true sea ice extent, SIE * , with an additive measurement error. 14 As discussed previously, some indicators present level shifts with respect to one another, most often resulting from how they respectively deal with tie points. It is thus preferable to constrain both the trend and dynamics to follow a common factor, while leaving level offsets unconstrained. The measurement equation is where Moreover, c i = 0 when we normalize  i = 1 for i  {S, J, B, G}. The transition equation is where  t  iid(0,  2  ) is orthogonal to  t at all leads and lags. Various modeling approaches are distinguished by their treatment of T REND t and SEASONAL t . We follow Diebold and Rudebusch (2019) and allow for 12 monthly deterministic seasonal effects, each of which is endowed with (possible) deterministic quadratic trend. 15 This results in a blended deterministic "trend/seasonal" given by where D i indicates month i and T IME indicates time. Hence the full transition equation is The model is already in state-space form, and one pass of the Kalman filter, initialized with the unconditional state mean and covariance matrix, provides the 1-step prediction errors necessary to construct the Gaussian pseudo-likelihood, which we maximize using the EM algorithm and calculate standard errors from the analytic Hessian matrix. 16 Following estimation, we use the Kalman smoother to obtain the best linear unbiased extraction of SIE * from the estimated model. The smoother averages across indicators, but it desirably produces optimally weighted averages rather than simple averages.
The smoother also averages over time, using data both before and after time t to estimate SIE * t , which is also necessary for optimal extraction, due to the serial correlation in SIE * . For details see Harvey (1989). One or more restrictions are necessary for identification. The standard approach is to normalize a factor loading, which amounts to an unbiasedness assumption. Normalizing  S =1, for example, amounts to an assumption that SIE S is unbiased for SIE * . Whether there truly exists such an unbiased indicator (and if so, which) is of course an open question -one can never know for sure. SIE S and SIE G are the most widely used indicators (Meier et al. (2014), Peng et al. (2013)), so it is natural to consider normalizing on  S or  G . We explore both. The estimated measurement equation ( 1), normalized with  S =1, is where standard errors appear beneath each estimated loading. All indicators are estimated to load heavily on SIE * , with all  s very close to 1. SIE J and SIE G load least heavily ( J =0.950, G =0.961)), in accord with their generally less abrupt trend in Alternatively, the estimated measurement equation normalized with The estimated loadings in equations ( 7) and ( 8), corresponding to  S =1 and  G =1 respectively, are effectively identical up to the normalization. Now consider the associated measurement error covariance matrix (3). The estimate for the with implied estimated correlation matrix Note that  2 GG is much higher than any of  2 SS ,  2 JJ , and  2 BB , potentially due to different indicators using different methods to determine tie points, i.e., reference points of brightness for 100% sea ice and 100% open water. The choice is crucial for accurate measurement of SIC within grid cells. Tie points, moreover, need not be constant, as brightness readings are sensitive to weather effects and atmospheric forcings (Ivanova et al. (2015)). Dynamic tie-point calibration is potentially desirable because it can decrease the bias of SIC measurements (Comiso et al. (2017)). The latest version of the Goddard Bootstrap algorithm, in particular, calibrates tie points daily. One would expect, however, that the bias reduction from dynamic tie-point calibration may come at the cost of potential discontinuities that increase measurement error variance. Our results confirm that conjecture. Our estimate of  2 GG is about 13 times that of  2 BB (which uses constant tie points). Alternatively, the estimated measurement error covariance matrix for the  G =1 normalization is 10 0.0361 with implied estimated correlation matrix Now let us move to the transition equation ( 6). Using the  S =1 normalization we obtain =0.704 [0.041] and trend/seasonal parameter estimates ( ai , b j , and ck ) as reported in the  S =1 columns of Trends for all months are highly significant and downward sloping. The trends for summer months (August-November) display a notable negative, and generally statistically significant, curvature, whereas for non-summer months the quadratic trend terms are generally small and statistically insignificant. For the  G =1 normalization we get =0.719 [0.042] and trend/seasonal parameter estimates ( ai , b j , and ck ) as reported in the  G =1 columns of In The fact that SIE * ( S =1) and SIE * ( G =1) are different, both in terms of level and variation around the level, limits their usefulness for research focusing on level, as the level depends entirely on identifying assumptions. However, and crucially, in an important sense SIE * ( S =1) and SIE * ( G =1) are highly similar: The model and identification scheme makes SIE * ( S =1) and SIE * ( G =1) identical up to a linear transformation. This is clear in Figure 4, which plots the two competing extracted factors, month-by-month. Regressions of SIE * ( G =1) on SIE * ( S =1) yield highly-significant intercepts not far from zero, highly-significant slopes near 1.0, and R 2 values above 0.999, for each month. Because SIE * ( S =1) and SIE * ( G =1) are identical up to a linear transformation, it makes no difference which SIE * we use for research focused on linear relationships between SIE * and other aspects of climate (e.g., various radiative forcings). The obvious choice, then, is SIE * ( S =1), which is just SIE S itself, dispensing with the need to estimate the factor model. We propose a dynamic factor model for four leading Arctic sea ice extent indicators. We estimate the model and use it in conjunction with the Kalman smoother to produce a statistically-optimal combination of the individual indicators, effectively "averaging out" the individual measurement errors. We explore two identification strategies corresponding to two different factor loading normalizations. The corresponding two extracted combined measures (latent factors) are identical up to a linear transformation, so either one can be used to explore relationships of sea ice extent and other variables. Interestingly, however, the extracted factor for one of the normalizations puts all weight on the Sea Ice Index. Hence the Sea Ice Index alone is a statistically optimal "combination" and one can simply use it alone with no loss, dispensing with the need to estimate the factor model. That is, there is no gain from combining the Sea Ice Index with other indicators, confirming and enhancing confidence in the Sea Ice Index and the NASA Team algorithm on which it is based, and similarly lending credibility -in a competition against very sophisticated opponents -to the NSIDC's claim that the Sea Ice Index is the "final authoritative SMMR, SSM/I, and SSMIS passive microwave sea ice concentration record" (Fetterer et al. (2017)).