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)). 1There 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 Figure 1. Of course, the Arctic SIE measures all trend down in every month, with steeper trends for the low-ice "summer" months (e.g., August, September, October). (Note the different axis scales for different months.) There are also systematic differences across indicators. SIE G , for example, tends to be high, and SIE J tends to be low, while SIE S and SIE B are intermediate. But the deviations between various pairs of measures are not rigid; that is, they are not simply parallel translations of each other. Instead, there are sizable time-varying differences among the various measures. 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. Table 1 summarizes the operative dates of the various satellites and sensors relevant for Arctic sea ice measurement. The first multi-frequency sensor equipped on a satellite was the Scanning Multichannel Microwave Radiometer (SMMR) launched in 1978 (Cavalieri et al. (1996)). Starting in 1987, later sensors -the Special Sensor Microwave Imager (SSM/I) and the Special Sensor  Microwave Imager/Sounder (SSMIS) -offered higher resolution images. 4 In 2002 and 2012, respectively, the Advanced Microwave Scanner Radiometer for EOS (AMSR-E) and Advanced Microwave Scanner Radiometer 2 (AMSR2) sensors were launched and provided further improvements in resolution (Comiso et al. (2017)). 5 Given the inclinations of the satellite orbits and the spherical shape of the earth, all of the satellites share an inability to observe the Arctic "pole hole" -a circular region at the very top of the world. The size of the pole hole varies across sensors, but historically, there is full confidence that the area covered by the pole hole fulfills the 15 percent SIC requirement (Meier and Stewart (2019)). Table 1 also describes the underlying source data for our four SIE indicators. These measures use algorithms to transform the raw satellite brightness data into SIC and SIE values. We now turn to a more detailed discussion of these algorithms to illuminate the differences across SIE indicators.
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 Goddard Space Flight Center -and (2) the Near-Real-Time DMSP SSMIS Daily Polar Gridded Sea Ice Concentrations (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 Goddard Space Flight Center, 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)). The Japan Aerospace Exploration Agency (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 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

12 See https://seaice.uni-bremen.de/sea-ice-concentration-amsr-eamsr2/time-series/.
13 Dynamic factor analysis is closely related to principal components analysis, but the dynamic factor model provides a fully-specified probabilistic modeling framework, in which estimation and factor extraction via the Kalman filter are statistically efficient. Under conditions the two approaches coincide in large samples, but those conditions include a very large number of indicators. Those conditions are violated in our case as we have only four indicators, so dynamic factor analysis is preferable. For a much more complete account, see Stock and Watson (2011). 14 The approach parallels Aruoba et al. (2016), who extract latent "true" U.S. GDP from noisy expenditure-side and income-side estimates. with

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 Figure 1. SIE B loads with an estimated coefficient marginally different from 1 but significant at the 5% level. Of course, the SIE S loading is 1 by construction. While Bremen's level offset (with respect to SII) is arguably negligible, those of Jaxa and especially Goddard are sizable. Hence, the estimation results accord with Figure 1, with SIE S and SIE B more in the center of the range and SIE J and SIE G more extreme. 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  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 ( âi , b j , and ĉk ) as reported in the λ S =1 columns of Table 2. 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 ( âi , b j , and ĉk ) as reported in the λ G =1 columns of Table 2. The λ S =1 and λ G =1 results are very similar. 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)). For a broad and insightful overview of the evolution and causes of reduced Arctic sea ice cover, seeShalina et al. (2020). See the notes for their Figure

3.3, page 3-13. Parkinson and Cavalieri (2008) discuss reasons for using a 15 percent cutoff. For detailed discussion of sensor characteristics see https://nsidc.org/ancillary-pages/ smmr-ssmi-ssmis-sensors. Early in the sample, operational problems prevented data delivery for several days during 1986 and between December 1987 and January 1988(Comiso (2017)). For more recent technical difficulties, see https://www.nrl. navy.mil/WindSat/Description.php. See https://doi.org/10.7265/N5K072F8. Maslanik and Stroeve (1999) takes the data from the National Oceanic and Atmospheric Administration Comprehensive Large Array-data Stewardship System (NOAA CLASS);(Cavalieri et al. (1996)) uses data processed at the NASA Goddard Space Flight Center. 8  For detailed algorithm description seeComiso et al. (2017b). For the data, see https://doi.org/10.5067/ 7Q8HCCWS4I0R. See also https://nsidc.org/support/faq/nasa-team-vs-bootstrap-algorithm. For description and data see https://kuroshio.eorc.jaxa.jp/JASMES/climate/index.html. Monthly data from https://seaice.uni-bremen.de/data/amsr2/today/extent_n_ 19720101-20181231_amsr2.txt. We emphasize that our model is meant to be a simple benchmark, and that many potentially-important variations and extensions are possible. For example, one could alternatively entertain stochastic as opposed to deterministic trend and seasonality. A simple approach would be separate month-by month modeling so that there is no seasonality, whether with one unit root, as in (for month m) T REND m,t = d m + T REND m,t-1 + u m,t , or two unit roots, as in T REND m,t = d mt + T REND m,t-1 + u m,t , where d m,t = d m,t-1 + v m,t . Note that we do not assume Gaussian shocks, and that it is not necessary to assume Gaussian shocks, as we can still maximize the Gaussian likelihood even if the shocks are not truly Gaussian, and the resulting (pseudo-)MLE still has the good properties of consistency, asymptotic normality, etc. Indeed when making Figure2we added a tiny constant (0.1) to SIE S to make it easier to distinguish SIE S from SIE * (λ S =1).