The Arctic is warming at least twice as fast as the global average (Osborne et al., 2018). This phenomenon of Arctic amplification in surface air temperature is closely connected to a dramatic multi-decade reduction in Northern sea ice. Indeed, since accurate satellite measurements began in 1978, the extent of Arctic summer sea ice has shrunk by about 40 percent, a loss in area comparable to the western continental United States. This drop in sea ice is one of the most conspicuous warning signs of ongoing climate change, but the reduction of sea ice will also play an important role in determining the pace of future global climate change and has wide-ranging implications for the polar region and the rest of the world. At a regional level, diminishing sea ice alters polar ecosystems and habitats and introduces major economic opportunities and risks. For example, new deposits of natural gas, petroleum, and other natural resources will become accessible for extraction, emission, and possible spillage (Petrick et al., 2017). Also, reduced ice coverage facilitates tourism and the use of Arctic shipping lanes, which are shorter than traditional passages via the Suez or Panama Canals. These new routes reduce sailing times but increase Arctic environmental risks from, for example, discharges, spills, and soot deposits (Bekkers et al., 2016). Finally, melting sea ice will have geopolitical consequences for Arctic sea-lane control (Ebinger and Zambetakis, 2009). Although these proximal Arctic effects are important, the far-reaching implications of diminished Arctic sea ice for regulating global climate and weather are even more consequential. Less sea ice and more open water diminishes the reflectivity (or albedo) of the Arctic region, so that, over time, a greater share of solar heat is absorbed by the earth, which leads to further Arctic amplification and increased temperatures worldwide (Hudson, 2011). These higher Arctic temperatures promote thawing and erosion of the polar permafrost, which can result in the release of large amounts of carbon dioxide and methane and provide a significant impetus to further global warming (Tanski et al., 2019). Increasing Arctic temperatures also hasten the melting of the Greenland ice sheet, further pushing up sea levels (Trusel et al., 2018). Finally, a warming Arctic and loss of sea ice cover can alter the global dynamics of ocean and air streams, and this effect already appears to be changing weather patterns at sub-polar latitudes (Petoukhov and Semenov, 2010) and weakening thermohaline ocean circulation including the current that warms Europe (Liu and Fedorov, 2019). In brief, the loss of Arctic sea ice is not only a stark signal of a changing climate, but it also plays an integral role in the timing and intensity of further global climate change. Not surprisingly then, the downward trend in Arctic sea ice has been the subject of hundreds of research studies. The forecasting literature alone is voluminous and impressive in both methodology and substance. 1 Although the importance of accurate polar prediction is hard to overstate, substantial uncertainty still remains about the future evolution of sea ice. Indeed, obtaining a deeper understanding of Arctic sea ice loss has been called a "grand challenge of climate science" (Kattsov et al., 2010). Much forward-looking sea ice analysis has been based on large-scale climate models, which represent of the fundamental physical, chemical, and biological drivers of the earth's climate. These models attempt to capture the dynamics of the oceans, atmosphere, and cryosphere at a high frequency and a granular level of geographic and spatial detail. Such structural physical models are invaluable for understanding climate variation, determining event and trend attribution, and assessing alternative scenarios. However, from a forecasting perspective, climate models have generally underestimated the amount of lost sea ice in recent decades (Stroeve et al., 2007;Stroeve et al., 2012b;Jahn et al., 2016;and Rosenblum and Eisenman, 2017). In addition, long-range sea ice projections can differ widely across climate models (Stroeve and Notz, 2015). Given the global significance of Arctic conditions and the progress yet to be made on specifying structural global climate models, purely statistical projections of sea ice are an obvious complementary approach. As a practical matter -across many disciplines -parsimonious statistical representations often produce forecasts that are at least as accurate as detailed structural models. Specifically for forecasting Arctic sea ice, there is already some evidence that small-scale statistical models with no explicitly embedded physical science can have some success (Guemas et al., 2016;Wang et al., 2016). Therefore, we provide a statistical analysis of the long-run future evolution of Arctic sea ice. Our work is distinguished by its use of intrinsically stochastic "unobserved components" models, with detailed attention to trend, seasonality, and serial correlation. Based on several decades of satellite data, we provide statistical forecasts of the future loss of Arctic sea ice. 2 Importantly, these forecasts provide probability assessments for a range of long-run outcomes and quantify both model parameter uncertainty and intrinsic uncertainty. Of particular interest are probability assessments of the timing of an ice-free Arctic, an outcome with vital economic and climate consequences (Massonnet et al., 2012;Snape and Forster, 2014;and Jahn et al., 2016). For this analysis, we also introduce a novel statistical modeling mechanism -a shadow ice interpretation -that allows us to readily account for the zero lower bound on the extent of Arctic sea ice in our model. Our resulting distributional forecasts suggest an ice-free Arctic summer is more likely than not within two decades -much sooner than the projections from many large-scale climate models. In particular, we contrast our statistical forecasts with projections from the ensemble of model simulations conducted for the fifth Coupled Model Intercomparison Project (CMIP5) -a highly-regarded central source for international global climate model projections. On average, these climate models envisage ice-free Arctic conditions close to the end of the century (assuming a range of business-as-usual carbon emissions paths). Thus, besides their relevance for environmental and economic planning, our probability assessments may also provide a useful benchmark for assessing or calibrating global climate models going forward. We proceed as follows.
In section 2, we introduce a linear statistical trend model and use it to produce long-range sea ice point forecasts. In section 3, we introduce the "shadow ice" concept to account for the zero-ice lower bound. In section 4, we generalize to a nonlinear (quadratic) statistical model and to interval forecasts that incorporate several forms of uncertainty. In section 5, we compare our statistical model forecasts to global climate model forecasts with particular attention to hard versus soft landings at zero ice. In section 6, we make probabilistic assessments of several sea ice scenarios. We conclude in section 7. Arctic sea ice has been continuously monitored since 1978 using satellite-based passive microwave sensing, which is unaffected by cloud cover or a lack of sunlight. For a polar region divided into a grid of individual cells, the satellite data provide cell-by-cell brightness readings, which can be converted into fractional ice surface coverage estimates for each cell. Sea ice extent, SIE -a very commonly-used measure of total ice area -is the total area of all cells with at least 15 percent ice surface coverage. That is, SIE rounds down cells with measured coverage of less than 15 percent to zero and rounds up cells that pass the 15 percent threshold to full coverage. The up-rounding in SIE is effectively a useful bias correction, as melting pools on summer ice surfaces can be mistaken for ice-free open water. For this reason, we follow common practice and use monthly average SIE data from November 1978 through October 2019 from the National Snow and Ice Data Center (NSIDC). 3 The NSIDC data use the NASA team algorithm to convert the satellite microwave brightness readings into measured ice coverage (Fetterer et al., 2017). 4,5  Figure 1 plots the time series of Arctic SIE -each monthly average observation is a dot -with an overall estimated linear trend superimposed. The clear downward trend is accompanied by obvious seasonality. A more subtle feature is the possible time variation in the seasonal effects, which may be trending at different rates and possibly nonlinearly. These effects turn out to be of interest in a complete statistical representation of sea-ice dynamics. A simple initial representation to capture this variation is a linear statistical model with twelve intercepts, one for each month, each of which may be differently trending, and potentially serially correlated stochastic shocks:

where the D i 's are monthly dummy variables (D it =1 in month i and 0 otherwise, i=1, ..., 12) and T IM E is a time dummy (T IM E t =t). 6,7 We estimate model (1) -and other versions below -by maximizing the Gaussian likelihood. Detailed regression results for model (1) are in column (6) of Table A1 in Appendix A. Figure 2 shows the resulting linear trends for all twelve months, highlighting March in blue and September in red. All of the monthly trends slope downward -an indication of a warming climate -and are highly significant. The slopes of the linear trends also differ across months (Serreze and Meier, 2019;Cavalieri and Parkinson, 2012). In particular, the seasonally lowice months of July through October have notably steeper downward sloping trends than the seasonally high-ice months of December through May. The estimated September trend, for example, is twice as steep as the March trend, and the difference is highly statistically significant. These linear trends are also extrapolated out of sample (shaded gray) through the end of the century. For example, September sea ice extent is projected to reach zero just after 2072. Such linear point forecasts are a useful first step, but they can be improved by allowing for nonlinearity in the trends and by quantifying forecast uncertainty -as described in section 4. First, however, we elucidate a "shadow ice" modeling approach that takes into account the fact that the measured amount of sea ice is bounded below by zero. One consideration for downward trending statistical models for Arctic sea ice is that the measured amount of sea ice will always be non-negative. By contrast, extrapolations of simple trend models will eventually push into negative territory. There are various functional forms that can be used to model such bounded time series, and the appropriate representation depends very much on the details of the real-world phenomenon under examination. 8 Some bounds act like reflecting barriers, so the variable of interest spends very little time at the constraint. Other bounds are absorbing states, and once reached, they may be sustained for some time. With positive amounts of sea ice, fluctuations in SIE can serve as a rough approximation for changes in the amount of thermal energy in the Arctic; that is, hotter and colder ocean surface temperatures are reflected in less or more ice, respectively. 9 However, this connection 8 One modeling approach is to rescale the bounded time series data to the real line using, say, a log-ratio transformation (Wallis, 1987). Alternatively, a time series can be modeled in the original bounded sample space using, for example, the beta autoregressive model of Rocha and Cribari-Neto (2008). By contrast, some physical science analyses simply terminate model simulations when zero ice conditions are reached, as in Obryk et al. (2019). 9 The amount of sea ice depends on a variety of factors including ocean and air temperature, ocean salinity, cloud cover, and wind, current, and wave action. However, several studies have identified ocean heat content as a major driver of sea ice coverage, notably, Arthun et al. (2012), Schlichtholz (2011), andSelyuzhenok et al. (2020). breaks down when the ice disappears: While SIE is fixed at zero, the surface temperature of the Arctic ocean can continue to warm. Furthermore, the warmer the ocean becomes, the less likely there will be a quick return of sea ice, which is indicative of a partially-absorbing state. 10To account for this effect, the negative values of sea ice produced by a statistical model can be viewed as a rough expression of ocean temperature.
Thus, we redefine the left-hand side variable of the unconstrained model as a shadow surface ice extent, SIE * . We view SIE * as a notional variable that equals measured surface ice when positive, but that may also go negative to represent further increases in ocean thermal energy more broadly. Formally, to translate negative model-based sea ice values into nonnegative sea ice observations, our shadow ice model modifies the unconstrained model (1) to respect the zero lower bound for ice:

That is, we now interpret our earlier unconstrained linear model of surface ice as a model of shadow ice, SIE * , so that the observed extent of sea ice, SIE, is the maximum of SIE * and zero.11

The shadow ice model respects the nonlinearity of observed ice at the zero lower bound but retains tractability. It also allows us to translate the long-range forecasts from downward trending models like model (1) -including distributional projections -into observed data that are always non-negative. As a matter of physical interpretation, very negative values of shadow ice extent, SIE * , represent environments in which the thermal content of the Arctic ocean is high enough that an immediate return to a positive SIE is unlikely. This shadow ice structure provides an intuitive and simple approximation of the thermodynamics of the Arctic ocean transition between sea ice and open water and serves as a useful modeling tool for observed SIE dynamics.12

A downward linear trend is a common representation of the secular decline in Arctic sea ice, but linearity is not assured by the physical science. There are a variety of climate feedback mechanisms that could hasten or retard the pace of sea ice loss. The well-known ice albedo effect occurs as sea ice cover is reduced, and the resulting darker ocean surface absorbs more energy, which in turn further reduces sea ice (Stroeve et al., 2012b). This feedback effect amplifies sea ice seasonality and may progressively steepen the downward trend in SIE over time (Schroeder et al., 2014). The geography of the perimeter of the Arctic Ocean, which is partially constrained by land that blocks expanding winter ice, becomes less relevant as sea ice shrinks, and relaxing this constraint may allow greater seasonal variation and a steeper downward sea ice trend (Serreze and Meier, 2019).13 However, there are offsetting negative feedback mechanisms -associated, for example, with increased cloud cover -that could slow the rate of sea ice loss over time (Portner et al., 2019). Indeed, Stroeve and Notz (2015) argue against trend amplification in favor of trend constancy (linearity) or trend attenuation. Moreover, as described in the next section, long-range SIE projections from large-scale global climate models appear dominated by feedback mechanisms that slow the rate of September sea ice loss over time. 14The lack of a complete understanding of the drivers of Arctic sea ice recommends consideration of a flexible empirical SIE specification, so we generalize from linear to quadratic trends: v t ∼ iid(0, σ 2 )

We label model (3) as the "general" quadratic model as no constraints are imposed on the twelve quadratic (α k ) parameters. The linear model (2) of course emerges under the constraint α 1 =...=α 12 =0. Figure 3 shows the SIE estimation fits and forecasts from the general quadratic model, by month, with March and September again highlighted in blue and red, respectively. The trend curvatures for all months are strikingly similar: the trends for all months decrease at an increasing rate. That is, the estimated α k coefficients are negative for every month, indicating that SIE is diminishing at an increasing rate. (Detailed estimation results of model (3) appear in column (1) of Table A1 in Appendix A.) The size of the estimated negative coefficients on the quadratic trend terms and their statistical significance vary by month. The most negative and significant α k coefficients are in the seasonally low-ice "summer" months of August, September, and October, and these months show the greatest trend rates of decline. 15 An F test of the joint hypothesis that α 8 =α 9 =α 10 =0 produces a p-value of 0.00. 16 Relative to a linear trend model, the nonlinear trend model forecasts lower sea ice at long horizons. February-April SIE point forecasts nevertheless remain well above zero through the century, but August-October point forecasts approach zero much more quickly. Indeed the quadratic September point forecast hits zero in 2045. Table 1 summarizes the results from investigation of the summer and non-summer differences in quadratic trend curvature using two standard model selection criteria: the Akaike information criterion (AIC ) and the Bayesian information criterion (BIC ). The AIC and BIC are estimates of out-of-sample forecasting performance (mean-squared error), formed by penalizing estimates of in-sample forecasting performance for degrees of freedom used in model fitting and differing only in the precise penalty applied (Diebold, 2007). The table reports these model selection criteria for six versions of the quadratic trend model (3) with various equality constraints imposed on the quadratic coefficients α 1 , ..., α 12 . Corresponding estimation results appear in Table A1 in Appendix A. The models favored by AIC and BIC -that is, those with smaller values -are very similar and involve summer and non-summer restrictions. The AIC selects a model with equal α k 's for the nine non-summer months (November-July) and unconstrained α k 's for the three summer months (August-October). The BIC, which penalizes degrees of freedom more harshly, selects a slightly more constrained model, with the non-summer α k 's again constrained to be equal and the three summer α k 's also constrained to be equal. Compared to the unconstrained version of equation ( 3) -the general model -both the AIC and BIC prefer specifications with some summer and non-summer equality constraints imposed, because their imposition saves degrees of freedom without substantially degrading fit.
All told, the results of Table 1 suggest a "simplified" quadratic model, namely model (3) with both summer and non-summer quadratic coefficients constrained separately to equality: α 8 =α 9 =α 10 and α 11 =α 12 =α 1 =...=α 7 . The regression diagnostics reported in Appendix A suggest good performance of the simplified quadratic model: R 2 is almost perfect; the residual graph and Durbin-Watson statistic (DW) are consistent with random noise; and residual skewness, kurtosis, and density estimates are consistent with normality. All told, the model appears to do a fine job of capturing SIE signal, reducing its complicated dynamics to Gaussian white

Figure 4 shows the simplified quadratic model fits and forecasts. The simplified model point forecasts are very similar to those of the general quadratic model in Figure 3, with a zero-ice September also reached in 2045, but the rank ordering of the months in terms of SIE is better preserved going forward. Going beyond these point forecasts, an important advantage of a formal statistical approach is that it can quantify the amount of future uncertainty. Figure 4 supplements the simplified quadratic trend point forecasts with interval or probability density forecasts. When making long-horizon interval forecasts, it is crucial to account for parameter estimation error, because its deleterious effects grow with the forecast horizon. For example, although parameter estimation error may have small effects on 6month-ahead intervals, it will be greatly compounded for 600-month-ahead intervals. An estimate of the time-t standard deviation of the forecast error, which accounts for parameter estimation error, is δ t =s 1 + x t (X X) -1 x t , where s is the standard error of the regression, x t is a 36×1 column vector of time-t right-hand-side variables, X is a T ×36 matrix whose columns contain the regression's right-hand-side variables over time, and T is sample size. 17We use δ t to produce the ± 2 δ t pointwise prediction intervals of Figure 4. Under normality of the shocks underlying the simplified quadratic model, the ± 2 δ t intervals are approximate 95 percent confidence intervals. 18 The intervals widen rapidly; indeed the September interval starts to include zero before 2040. 19Finally, in Figure 5 we zoom in on the 2060-2062 (36-month) segment of the simplified quadratic model shadow ice forecast. The point forecast trends down and is seasonally below zero by then. (Indeed, as already discussed, the point forecast is seasonally below zero well before then.) Moreover, the ± 2 s.e. intervals widen over time, and entire intervals are seasonally below zero by then. Hence it appears that, with near certainty, summer SIE will vanish by 2060. This result is also clear from our earlier Figure 4, but Figure 5 highlights it in a different and complementary way. Of the many analyses of the long-term future evolution of Arctic ice, most have focused on projections from large-scale climate models. Such models are based on the underlying physical, chemical, and biological processes that govern the dynamics of weather and climate across ocean, air, ice, and land. The models fit an immense number of variables at a high temporal frequency and a granular spatial scale (e.g., a 30-minute time interval. a 100km worldwide horizontal resolution, and 40 vertical levels in the oceans and atmosphere). Dozens of scientific groups around the world have constructed and currently maintain such models. Occasionally, these groups conduct concurrent simulations as part of the Coupled Model Intercomparison Project (CMIP) that involve common sets of inputs including carbon emissions scenarios. The most recently completed iteration or phase of this project is the fifth one, denoted CMIP5 (Taylor et al., 2012). The CMIP5 model comparison study was the main source of climate projections included by the International Panel on Climate Change (IPCC) in its landmark Fifth Assessment Report. Arctic SIE is a key variable projected by climate models, and the models included in CMIP5 are generally judged to provide a better fit to Arctic sea ice than earlier CMIP itera-tions. 20 Figure 6 shows three projections for September SIE constructed as averages across sets of CMIP5 global climate models. These multi-model mean projections are constructed under three different scenarios, or Representative Concentration Pathways (RCPs), for future greenhouse gas concentrations. The brown, yellow, and blue lines are averages of the climate models, respectively, under a high level of carbon emissions (RCP8.5), a medium level (RCP6.0, yellow), and a low level (RCP4.5, blue). 21 The first two higher emissions scenarios are viewed as more likely business-as-usual outcomes and are the most relevant to compare to statistical projections that extend the historical sample of past data and assume a continuation of the world economy's current population and development trajectories. 22  The solid red line in Figure 6 shows the September trend from the simplified quadratic model estimated on the full sample from November 1978 to October 2019. The dotted lines bracketing this forecast provide an approximate 95 percent confidence interval. However, the CMIP5 climate model projections are only based on data through 2005 and do not include the past dozen or so years of sea ice observations. For comparability to these climate model projections, we re-estimated the simplified quadratic model using data from 1978 through 2005, and the September trend from this pre-2006 model is shown as the dashed red line. An interesting first result is that from the close conjunction of the two simplified model statistical projections (the solid and dashed red lines), the addition of data from 2006 to 2019 does not lead to a significant revision in the statistical trend model. The very modest differences between the pre-2006 estimated quadratic trend and the full-sample version is an indication of the stability and suitability of the simplified statistical model. Comparing the statistical and climate model projections in Figure 6 reveals two salient features. First, throughout the forecast period, the full-sample statistical model projection is significantly lower than any of the climate multi-model mean projections.
Indeed, the climate model mean projections are well outside the 95 percent confidence intervals. Relative to the statistical model estimated on data before 2006 or the full-sample version, the climate model means are overestimating the 2019 SIE trend by about 1.0 million km 2 . This wedge is 20 Stroeve et al. (2007) describes the poor sea ice fit of the CMIP3 models, while the somewhat improved fit of the CMIP5 models is noted by in Stroeve et al. (2012b). 21 The climate model data are described in Serreze and Meier (2019) and Stroeve et al. (2012b) and were kindly provided by Andrew Barrett at the National Snow and Ice Data Center. The model sets averaged for each scenario are not identical, with the RCP4.5, RCP6.0, and RCP8.5 scenarios based on 26, 8, and 25 climate models, respectively. 22 In the RCP8.5 scenario, continuing increases in greenhouse gas emissions through the end of the century raise the 2100 global average temperature by about 4.0-6.0 • C above pre-industrial levels (Reidmiller et al., 2018). In the RCP4.5 scenario, greenhouse gas emissions level off before mid-century, and the 2100 global average temperature is approximately 2.0-3.0 • C above pre-industrial levels. projected to increase dramatically over the next two decades. Indeed, the pre-2006 and fullsample statistical trend models project zero ice in 2042 and 2044, respectively, but none of the climate model mean projections reach a completely ice-free Arctic in this century. Notably, even assuming a high level of emissions (RCP8.5) -which is a scenario with continuing increases in average global surface temperatures throughout this century -the multi-model mean projection never reaches zero Arctic summer sea ice. Second, in contrast to the statistical projection, the climate model mean projections show a decreasing rate of ice loss over time -that is, a concave rather than convex structure. Specifically, the climate model projections all display a roughly linear decline for the first couple of forecast decades (through around 2040) and then start to level out. For the lower emissions RCP4.5 scenario, the leveling out and deceleration of sea ice partly reflects a slowdown in the pace of global temperature increases. This effect is not at work in the RCP8.5 scenario, which has global temperatures steadily climbing through 2100. However, very close to zero sea ice extent, the leveling out of SIE in RCP8.5 appears reflect the hypothesized difficulty of melting the thick sea ice clinging near northern coastlines -notably in Greenland and Canada. Climate models generally assume that these coastal regions will retain landfast sea ice for a time even after the open Arctic Sea is free of ice. Therefore, a common definition of "ice-free" or "nearly ice-free" in the literature is a threshold of 1.0 million km 2 rather than zero SIE (Wang and Overland, 2009). Still, even with this higher threshold, the three climate model mean projections only reach a nearly ice-free Arctic in 2068, 2089, and >2100 for successively lower emissions scenarios, respectively. By contrast, the pre-2006 statistical projection reaches the higher 1.0 million km 2 threshold in 2037, and the full-sample statistical model reaches that level in 2039. Some have argued that climate models generally do well in representing the large-scale evolution of Arctic sea ice (Stroeve and Notz, 2015), but a number of studies have noted that the CMIP5 global climate models overpredicted the amount of Arctic sea ice (Massonnet et al., 2012;Stroeve et al., 2012a;Serreze and Meier, 2019). That overprediction continues, and its source is not well understood. One proposed correction to this overprediction has been to focus on the models that fit the historical SIE data better according to certain metrics (Wang and Overland, 2012). However, there is no agreed upon model selection criterion. Also, Rosenblum and Eisenman (2017) discount such model selection because models with more accurate sea ice readings also tend to overpredict global temperatures, so the selected climate models may be getting sea ice loss right for the wrong reason. Finally, it should be noted that from a statistical viewpoint, focusing on a simple average forecast from many models has been shown to be a robust prediction strategy (Diebold and Shin, 2019). In some respects, the wide differences between the statistical and climate model projections are not too surprising. In climate models, the monthly observations on total Arctic SIE are a high-level output from complex, nonlinear, granular representations of the relevant underlying science. Obtaining good SIE predictions from these models requires correctly specifying a host of detailed subsidiary processes. Such a bottom-up modeling procedure has important advantages in structural interpretation and counterfactual scenario analysis. However, in a variety of disciplines, a bottom-up procedure, which carries the possibility that small misspecifications can accumulate and affect high-level aggregates, has not been found to improve prediction relative to a top-line procedure that directly models the object of interest (Diebold, 2007). Thus, based on broad previous experience, we believe that direct statistical projections of Arctic SIE are likely to be relatively accurate. An advantage of a formal statistical model is its ability to make density forecasts for a range of possible reduced Arctic sea ice scenarios. Of particular interest are probability assessments of an ice-free or nearly ice-free Arctic, that is, the probability that SIE t equals zero or is less than or equal to some threshold γ, respectively. Formally, such event probabilities can be denoted as P (SIE t ≤γ), which represents the probability that sea ice extent is less than or equal to γ in month t. We estimate these scenario probability distributions using the simplified quadratic model and a stochastic simulation procedure that accounts for parameter estimation uncertainty and allows for potentially non-Gaussian serially correlated stochastic shocks. From a given set of simulated paths, we estimate the event probabilities of interest as the proportion of simulated paths in which the event occurs out of the total number of paths. 23An event that has attracted much attention in the literature is the initial occurrence of an ice-free or nearly ice-free September.
We calculate the probability for each September with date t 0 so that SIE t 0 ≤γ and SIE t >γ for all t<t 0 . Specifically, for a given γ and simulation i, we determine the year in which September SIE t first reaches γ, and then we cumulate across all simulations to build a distributional estimate. The red lines in Figure 7 provide the resulting probability distributions for an initial "ice-free" Arctic September for γ=0, γ=1, and γ=2, that is, for progressively more lenient definitions of "ice-free." As noted above, the middle value of γ, which represents Arctic sea ice of less than 1 million km 2 , is a popular benchmark in the literature (Wang and Overland, 2012). For that particular threshold, the statistical model produces a distribution centered at 2039. 24 Of particular interest is the probability distribution of dates taking into account model parameter uncertainty and stochastic shock uncertainty, and this distribution shows about a 60 percent probability of an effectively ice-free September Arctic occurring in the 2030s. The distribution with γ=1 is bracketed on either side by distributions that use the higher and lower thresholds. The γ=0 distribution has a median date of 2044 and a 95 percent range from 2039 to 2053. The γ=2 distribution has a median date of 2033 and an earlier and slightly narrower 95 percent range from 20230 to 2039. The climate modeling literature has pointed to several factors that underpin the uncertainty in the timing of an initial September ice-free Arctic including natural climate variability, emissions path uncertainty, and uncertain sea ice dynamics (Serreze and Meier, 2019). These factors are at least partially accounted for in our analysis. The multi-model mean CMIP5 climate projections described above are well outside the above-described distributional ranges. For a threshold of 1 million km 2 , the mean projections reach a nearly ice-free Arctic in 2068 and 2089 under the RCP8.5 and RCP6.0 scenarios, respectively. The former date is denoted by a vertical brown bar in Figure 7. The range of dates across individual models is also extremely wide, stretching well over a century (Jahn et al., 2016). To narrow this range, researchers have omitted models with poor performance using a variety of sea ice metrics. With such model selection procedures, the range of dates for a first nearly ice-free September is narrowed greatly to a 20-year span that runs from the 2040s to 2060s (Massonnet et al., 2012;Thackeray and Hall, 2019). 25 Even such a carefully circumscribed span is roughly a decade later than the simplified quadratic trend model produces. Moreover, the climate model simulation exercises are not designed to yield formal error estimates or measures of uncertainty as the spread of climate model forecasts in an ensemble is insufficient to completely characterize forecast uncertainty. Finally, we note the availability of density forecasts for a variety of richer joint scenarios of interest. As one example, Jahn et al. (2016) describe other definitions of "ice-free" that involve, for example, 5-year running means. Alternatively, "ice-free" may require no ice for several consecutive months, for example, to accommodate meaningful freight shipping, tourism, mining, and commercial fishing (Yevgeny et al., 2017). For example, the strong Autumn demand for international freight shipping to satisfy year-end Western holiday consumer demand could make a multi-month ice-free Arctic shipping lane of interest. In this case, the probability distribution of the initial occurrence of an "ice-free" summer -a joint ice-free August, September, and October -could be relevant. Figure 7 shows this distribution in black assuming an ice-free threshold of γ=1. This density is notably shifted rightand more right-skewed -compared to the September scenarios. A rapidly warming Arctic is an ominous sign of the broader climate change caused by human activity, but declining Arctic sea ice also has an important influence on the pace and intensity of future climate change. Using statistical models, we have provided probabilistic projections of 21st-century Arctic sea ice that account for both intrinsic uncertainty and parameter estimation uncertainty. These projections indicate that summertime Arctic sea ice will quickly diminish and disappear -with about a 60 percent probability of an effectively ice-free September Arctic occurring within two decades. By contrast, the average projection from leading climate models implies an initial seasonal ice-free Arctic several decades later -even assuming a business-as-usual emissions path. The slow and decreasing pace at which large-scale climate models reach an ice-free Arctic may be a serious shortcoming, and such conservative projections of sea ice loss could be a misleading guide for global climate policy. 26There are numerous examples, across many disciplines, showing that parsimonious statistical representations can provide forecasts that are at least as accurate as the ones from detailed structural models. 27 However, rather than treat statistical models as just forecast competitors to structural climate models, there is likely to be scope to use them as complementary representations going forward. The mechanisms governing Arctic sea ice loss and connecting that loss to atmospheric, oceanic, and permafrost responses are not fully captured in climate models. Statistical models may be able to help assist in bridging such gaps until a more complete understanding is available (Parkinson, 2010). In addition, statistical models may also be able to provide a benchmark for model performance that can be used for climate model evaluation, estimation/calibration, and selection. In particular, our statistical approach summarizes a key aspect of the data, namely its historical dynamics and the associated forecast path (call it P AT H D , where D stands for "Data"). A climate model with calibrated parameters θ produces a corresponding summary for the model (call it P AT H M (θ), where M stands for "Model"). So, at least conceptually, we can consider In practice such formal procedures are infeasible for massive global climate models, which attempt to address all aspects of the high-dimensional climate state, of which sea ice is but one small part.
Hence, for example, the "right" window for complete climate model estimation by indirect inference is much more complicated than just a sea ice path. Less formal procedures directly motivated by the above considerations can nevertheless be highly revealing, as in our key Figure 6.

representations. 27 A classic example from economics is Nelson (1972), who showed that simple statistical models forecast the economy as well as large-scale structural models based on economic theory. 3), with various equality constraints imposed on the quadratic coefficients α 1 , ..., α 12 . "NONE" denotes no constraints, which corresponds to the general quadratic model (3). "Seq" ("Summer equal") denotes August-October equal (α 8 =α 9 =α 10 ). "NSeq" ("Non-Summer equal") denotes November-July equal (α 11 =α 12 =α 1 =...=α 7 ). "Seq+NSeq" denotes summer months equal and non-summer months (separately) equal, which corresponds to the simplified quadratic model. "ALLeq" denotes all months equal (α 1 =...=α 12 ). "ALL0" denotes all months 0 (α 1 =...=α 12 =0), which corresponds to the linear model (1). In each column, we show constrained parameter estimates in boldface. "*" denotes significance at the ten percent level. We estimate scenario probability distributions using the simplified quadratic model (that is, the quadratic model (3) subject to the constraints α 8 =α 9 =α 10 and α 11 =α 12 =α 1 =...=α 7 ) and a simulation procedure that accounts for parameter estimation uncertainty and allows for potentially non-Gaussian serially correlated stochastic shocks. To describe this simulation procedure, it is useful to re-write the quadratic model (3) in a more concise notation: ) and β =(δ 1 , ..., δ 12 , γ 1 , ..., γ 12 , α 1 , ..., α 12 ). We estimate the model by maximum likelihood using the historical data sample t = 1, ..., T , as discussed earlier, and then we simulate i = 1, ..., 10, 000 future paths based on the estimated model. Simulation i proceeds as follows:

1. Draw β (i) from the sampling distribution of the estimated parameter vector β, and form x t β (i) , t = T +1, ..., 2099M 12.

by sampling with replacement from the observed vt 's (with equal weights 1/T ), t = T +1, ..., 2099M 12. 3. Draw ρ (i) from the sampling distribution of ρ, and form ε

t , t = T +1, ..., 2099M 12, using the last observed historical residual εT as the initial condition. Finally, we estimate event probabilities of interest as the relative frequency of occurrence across the simulated paths. Another measure of Arctic ice is sea ice area (SIA), which adds together the measured fractions of ice-covered areas of all cells that pass the 15 percent threshold. For brevity, we report results only for SIE, but results for SIA are qualitatively identical. We interpolate the missing December 1987 and January 1988 observations with fitted values from a regression on trend and monthly dummies estimated using the full data sample. For a comparison of algorithms, and evidence of strong performance of the NASA team algorithm, seeDiebold et al. (2020). Earlier work in this tradition includesPeng et al. (2018). Note that although we allow for seasonal and trending intercepts, we assume a constant disturbance variance across months and time. There are so many interesting directions that we have reserved volatility modeling for future work. For example, high-frequency (daily) SIE data are available, from which one could accurately estimate, and then model, lower-frequency (monthly) realized SIE volatility. See Stroeve and Notz (2018)  for a description of the phenomenon by which the ocean must release heat back into the atmosphere before sea ice can form again in winter. A similar framework has been successfully applied in finance to model nominal interest rates near their zero lower bound(Christensen and Rudebusch, 2014;Bauer and Rudebusch, 2016). Wang et al. (2016) take a different approach by simply constraining the model by the lower bound (for sea ice concentration in their case), so negative predicted values are set to zero. In a dynamic model with lagged sea ice, this procedure will result in a representation of a physical bound that is much closer to a reflecting barrier. Currently, the southern limit of winter ice is bounded by land except in the Bering Sea, Sea of Okhotsk, East Greenland Sea, Barents Sea, and Baffin Bay. More extreme forms of nonlinearity -such as discontinuous breaks, tipping points, and thresholds -are possible but viewed as less likely(Stroeve and Notz, 2015).Goldstein et al. (2018) argue that Arctic sea ice is best modeled by step-like shifting means at fitted breakpoints. However, modified statistical information criteria that properly account for the implicit flexibility of such breakpoints -followingHall et al. (2013) do not favor such a shifting mean models relative to a linear trend. Minimum ice coverage occurs at the end of summer in September because sea ice coverage continues to shrink for any above-freezing temperatures. For convenience, we refer to the low-ice months of August, September, October as "summer" even though they span late summer through early autumn. Because T IM E 2 and T IM E are correlated, an insignificant α k coefficient would not necessarily imply that nonlinearity is unimportant. See for exampleJohnston and DiNardo (1972), pp. 153-155, for derivation of this canonical result. Normality of the v shocks does not appear unreasonable, as the simplified quadratic model residuals have skewness and kurtosis of -0.17 and 3.73, respectively. We also obtained similar results without assuming normality via bootstrap simulation, which we discuss in section 6. By contrast, interval forecasts that fail to account for parameter estimation uncertainty quickly approach (by about 12 months ahead) the fixed-width interval ± 2 δ, where δ is the estimated unconditional standard deviation of the AR(1) disturbance in equation (3), σ2 /(1-ρ 2 ) ≈ 0.32 and fail to widen with forecast horizon. We provide details on this simulation procedure in Appendix B.
For further discussion of the econometrics of threshold event probabilities, seeBauer and Rudebusch (2016). Both the mean and median round to 2039 despite a slightly skewed distribution toward longer times. Again, the caveats ofRosenblum and Eisenman (2017) regarding such model selection are relevant. Our statistical results are in line with the concerns ofStroeve et al. (2007) regarding the possibility that the slow projected decline in Arctic sea ice by climate models suggests they are underestimating the effects of greenhouse gases. More recently,Guarino et al. (2020) also emphasize deficiencies in climate models' Arctic