Arctic Amplification of Anthropogenic Forcing: A Vector Autoregressive Analysis
aa r X i v : . [ ec on . E M ] M a y Arctic Amplification of Anthropogenic Forcing:A Vector Autoregressive Analysis
Philippe Goulet Coulombe Maximilian Göbel University of Pennsylvania ISEG - Universidade de Lisboa
First Draft: December 11, 2019This Draft: May 7, 2020Latest Draft Here
Abstract
Arctic sea ice extent (SIE) in September 2019 ranked second-to-lowest in history and is trend-ing downward. The understanding of how internal variability amplifies the effects of external CO forcing is still limited. We propose the VARCTIC, which is a Vector Autoregression (VAR)designed to capture and extrapolate Arctic feedback loops. VARs are dynamic simultaneous sys-tems of equations, routinely estimated to predict and understand the interactions of multiplemacroeconomic time series. Hence, the VARCTIC is a parsimonious compromise between full-blown climate models and purely statistical approaches that usually offer little explanation ofthe underlying mechanism. Our "business as usual" completely unconditional forecast has SIEhitting 0 in September by the 2060’s. Impulse response functions reveal that anthropogenic CO emission shocks have a permanent effect on SIE – a property shared by no other shock. Further,we find Albedo- and Thickness-based feedbacks to be the main amplification channels throughwhich CO anomalies impact SIE in the short/medium run. Conditional forecast analyses revealthat the future path of SIE crucially depends on the evolution of CO emissions, with outcomesranging from recovering SIE to it reaching 0 in the 2050’s. Finally, Albedo and Thickness feed-backs are shown to play an important role in accelerating the speed at which predicted SIE isheading towards 0. Department of Economics, [email protected]. ISEG, Universidade de Lisboa.
Introduction
In 2019, the minimum extent of Arctic sea ice ranked second-to-lowest in history, with thelowest occurring in September 2012. A persistent retreat of SIE may further accelerate globalwarming and threaten the composition of the Arctic’s ecosystem (Screen and Simmonds (2010)).While depriving many people from their traditional livelihood, a retreating ice cover has al-ready offered new shipping routes and oil exploration projects over recent years (Meier et al.,2014), increasing business activity in the region. M OTIVATION . The Coupled Model Intercomparison Project (CMIP), assembles estimates oflong-run projections of Arctic sea ice from a large number of climate models. These modelstry to reproduce the geophysical dynamics and interrelations among various variables, influ-encing the evolution of global climate. With CMIP now being in its 6 th phase (CMIP6), climatemodels provide more realistic forecasts of the Arctic’s sea ice cover compared to its predeces-sor CMIP5 (see Stroeve et al. (2012), Notz et al. (2020)). The majority of contributors to CMIP6see the Arctic’s sea ice to retreat below the 1 × km mark before the year 2050. Despitefollowing the hitherto accepted physical laws of our climate, the chaotic nature of the latter,i.e. the still obscure interplay of various climate variables, imposes a major burden on cli-mate models. In trying to replicate the observed behavior of our climate, each model is re-runseveral times with differing initial conditions resulting in a wide range of projections of keyclimate variables (Notz et al., 2020). In addition to this tuning, these simulations require largeamounts of input data and a coupling of various sub-models (Taylor et al., 2012).The above raises the question whether an approach that is statistical and yet multivariate can paint a more conciliating picture. This means estimating a statistical system that depictsthe interaction of key variables describing the state of the Arctic. In such a setup, the down-ward SIE path will be an implication of a complete dynamic system based on the observedclimate record. We can provide a formal statistical assessment of different hypotheses aboutthe historical path of SIE and their implications for the future. We can quantify a specific effect– and the uncertainty surrounding its estimation – without resorting to use a climate model. F EEDBACK L OOPS . Our analysis focuses on the interaction between anthropogenic carbondioxide ( CO ) forcing and internal variability. The former is already widely suspected tobe the main driver behind long-run SIE evolution (see Meier et al. (2014), Notz (2017)). In-ternal variability consists of feedback loops that are well documented in the literature (seeParkinson and Comiso (2013), Winton (2013), Stuecker et al. (2018), McGraw and Barnes (2020))and are crucial for further enhancing the predictability of the Arctic’s sea ice cover (Wang et al.(2016), Notz et al. (2020)). It is clear that only an approach that considers the interaction ofmany variables in a flexible way – and thus numerous potential sources for feedback loops– has a chance to depict a reliable statistical portrait of the Arctic. The still high variation inCMIP6 Arctic sea ice projections (see Notz et al. (2020)) suggests that there is widespread un-ertainty around the question to what extent internal variability can amplify external forcing.In order to quantify this, we present a methodology from economics to achieve just that. T HE V ARCTIC . Our analysis focuses on the evolution of the long-term trajectory of SIE andthe interconnected forces behind it. The modeling approach we propose achieves a desirablebalance between purely statistical and theoretical/structural approaches. In many fields sta-tistical approaches have a much better forecasting record than theory-based models. An obvi-ous drawback is that the successful model provides little to no explanation of the underlyingmechanism. A Vector Autoregression (VAR) lives in a useful middle ground. It is a statistical model,but yet generates forecasts by iterating a complete system of difference equations in multipleendogenous variables. These interactions can be analyzed and provide an explanation for theresulting forecasts. When estimated with Bayesian techniques, VARs are known to providecompetitive forecasts – very often as good as black-box models (Giannone et al. (2015)). Inthe light of all these considerations, we propose the VAR for the Arctic (VARCTIC) statisticalapproach that can (i) generate long-run forecasts, (ii) explain them as the result of feedbackloops and external forcing (iii) allows us to analyze how the Arctic responds to exogenousimpulses/anomalies. P REVIEW OF R ESULTS . Our "business as usual" forecast has SIE reaching 0 in September ofthe 2060s and predicts SIE to be below one million km by the mid 2050s. By studying theimpulse response functions of a Bayesian VAR, we report that CO shocks have the uniqueproperty of an everlasting impact on SIE. Additionally, we document that the correspondingresponses of Albedo, Air Temperature and Thickness largely amplify the middle-run impactof CO anomalies on SIE. Further, we use conditional forecast analysis to evaluate the long-runeffect of the systematic CO increase. We consider three different CO emission scenarios andshow that abiding by the Paris Accord could eventually bring back SIE to 2010s levels. Thetwo other standard CO paths lead SIE to 0 in the 2050s or the 2070s. We find that internalvariability as characterized by Albedo and Thickness feedbacks, while not the original sourceof the decay, amplify external forcing greatly. In the worst-case scenario for the CO path,canceling these forces starting from 2020 would postpone reaching 1 × km by a decade. R OADMAP . We first discuss the data and its transformation in section 2. Secondly, we discussthe VAR model, its identification and Bayesian estimation in section 3. Section 4 containsthe empirical results which comprise of (i) a long-run forecast of SIE, (ii) impulse responsefunctions of the VAR, (iii) exploring the transmission mechanism (feedback loops) and (iv)conditional forecasting analysis. This is also a well-known concern in Machine Learning, which generated an ongoing interpretability litera-ture that aims at opening the so-called black box. Via informative priors, Bayesian shrinkage can help in reducing the variance of the densely parameterizedVAR estimates – the same way LASSO or Ridge regularization helps in Machine Learning. Data
Our data set comprises eighteen time series, proxying the Arctic’s climate system, and account-ing for potential feedback loops among different components. The sample covers monthlyobservations from 1980 through 2018. Even if the accuracy of estimates of various variableshas improved over the last decades, uncertainties and measurement errors still remain. Inparticular, providing data for long-term analysis often requires the merging of data sources ofdifferent quality and reliability (Meier et al., 2014). Nevertheless, our data is derived fromwell accepted sources (see Stroeve and Notz (2018)), such as the
National Snow & Ice DataCenter’s (NSIDC)
Sea Ice Index series (Fetterer et al. (2017)), which proxies SIE, NASA’s
Mod-ern Era Retrospectiveanalysis for Research and Applications (MERRA-2) for atmospheric variables(Gelaro et al., 2017), or the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS)(Zhang and Rothrock, 2003).We combine 8 variables, which importance has been highlighted by the existing literature(Meier et al. (2014)), into VARCTIC 8, our benchmark specification. Fortunately, variables caneasily be added/removed from a VAR and Bayesian shrinkage ensures that a larger modelwill not overfit – the latter aspect is further explained in section 3.5. Therefore, we considerin the appendix a VARCTIC 18 which includes an additional 10 series from the reanalysis prod-uct MERRA2 (Gelaro et al. (2017)) as a robustness check. To summarize compactly, the twospecifications considered in this paper are:I
VARCTIC 8 : CO , TotCloudCover, PrecipCMAP, AirTemp , SeaSurfTemp, SIE, Thick-ness, Albedo;II. VARCTIC 18 : SWGNT, SWTNT, CO , LWGNT, TotCloudCover, TAUTOT, PrecipCMAP,TS, AirTemp, SeaSurfTemp, LWGAB, LWTUP, LWGEM, SIE, Age, Thickness, EMIS, Albedo.A comprehensive overview of all variables (including those of VARCTIC 18), their acronymsand links to data providers can be found in the appendix in Table 1. The justification for thesubset of these variables included in VARCTIC 8 and the extensions of VARCTIC 18 will bediscussed extensively in section 3.3.The raw data is highly seasonal as displayed in Figure 15. However, the feedback loopswe wish to estimate and extrapolate reside in the (stochastic) trend components and short-runanomalies, which is a limited part of the SIE’s variance. That is, we wish our VAR to explainlong-run fluctuations and short-run variability rather than seasonality. Hence, we proceedto transform the data so that the resulting VARCTIC is fitted on deviations from seasonalmeans.For our benchmark analysis, we use a simple and transparent transformation: we de-seasonalized our data by regressing a particular variable y raw on a set of monthly dummies.3hat is, for each variable we run the regression y rawt = ∑ m = α m D m + residual t and y t is defined as y t ≡ y rawt − ∑ m = ˆ α m D m . D m is an indicator that is 1 if the date t is in month m and 0 otherwise. The estimates of α m ’s, ˆ α m ’s, are obtained by ordinary least squares. Thisis exactly equivalent to de-meaning each data series month by month and is a more flexibleapproach to modeling seasonality than using Fourier series. Finally, we keep our filtereddata y in levels. We do not want to employ first differences or growth rate transformationsto make the data stationary. Such an action would suppress long-run relationships whichare an important object of interest. Figure 1 shows the data after being filtered with monthlydummies. Figure 1: Deseasonalized Series: 8 Variables
CO2
TotCloudCover
PrecipCMAP
AirTemp
SeaSurfTemp
SIE
Thickness
ALBEDO
Pre-processing the data can influence results. Moreover, Diebold and Rudebusch (2019)and Meier et al. (2014) document seasonal variability in SIE trends. As a natural robustnesscheck, we also consider a very different approach to eliminate seasonality. In appendix A.6,we reproduce our results with a data set of stochastically de-seasonalized variables obtainedfrom the approach of structural time series (Harvey (1990) and Harvey and Koopman (2014)).In short, this extension allows for seasonality to evolve (slowly) over time, which could be afeature of some Arctic time series. Of course, if we were using higher-frequency data – like daily observations, then the Fourier approach wouldbe much more parsimonious and potentially preferable (Hyndman, 2010). The dummies approach to taking outseasonality only requires 12 coefficients with monthly but 365 with hypothetical daily data. CO was not de-seasonalized. The VARCTIC
Reenforcing feedback loops have been the subject of countless climate studies. However, re-cently proposed statistical approaches to model them – as advocated in McGraw and Barnes(2020) and McGraw and Barnes (2018) – are yet incomplete in their quest to fully unlock thepotential of macroeconometrics. In this section, we review the VAR: the model; its identifica-tion; its Bayesian estimation. Furthermore, we discuss the construction of the long-forecastsand impulse response functions as tools to understand the VARCTICs’ results.
Vector Autoregressions are dynamic simultaneous systems of equations. They can character-ize a linear dynamic system in discrete time. The methodology was introduced to macroeco-nomics by Sims (1980) and is now so widely used that it almost became a field of its own (seeKilian and Lütkepohl (2017)). It is a multivariate model in the sense that yyy t in Ayyy t = Ψ + P ∑ p = Ψ p yyy t − p + εεε t , (1)is an M by 1 vector. This means the dynamic system incorporates M variables. Each of thesevariables is predicted by its own lags and lags of the M − A characterizes how the M different variables interact contemporaneously. Finally, thedisturbances are mutually uncorrelated disturbances with mean zero: εεε t = [ ε t , ... , ε M , t ] ∼ N ( I M ) .Equation (1) is the so-called structural form of the VAR which cannot be estimated because A is not identified by the data. An estimable version is the reduced-form VAR yyy t = ccc + P ∑ p = Φ p yyy t − p + uuu t , (2)where ccc = A − Ψ , Φ p = A − Ψ p and uuu t are now regression residuals uuu t = [ u t , ... , u M , t ] ∼ N ( Σ u ) with Σ u = A − ′ A − by construction. While standard, obtaining an estimate of A by decom-posing Σ u is more complicated than running a simple regression is addressed on its own insection 3.3.The methodology has many advantages over simple autoregressive distributed lags (ARDL)regression that have gained some popularity in the econometric and climate literature. In test-5ng for the importance of blocks of lags, McGraw and Barnes (2018) argue in favor of doingGranger causality tests rather than running ARDL regressions that exclude lags of the depen-dent variable. In essence, this is a well-known omitted variable bias that led eventually to theadoption of VARs by the macroeconometric profession. For instance, if the true data gener-ating process is a VAR for each of the M equations, excluding lags of y m , t in equation m willlead to spurious results. The advantages the authors described for Granger causality tests arethus de facto included in a VAR analysis. Their argument for inclusion of lags of the depen-dent variable can be interpreted as one for completeness of the modeled dynamic system, asguaranteed by an adequately specified VAR.
The symmetry of the VAR allows for it to generate forecasts by simply iterating the model. Assuming the chosen variables to characterize the system completely, we can forecast its futurestate by iterating a particular mapping. To do so, we use a representation that exploits thefact that any VAR( P ) (that is, with P lags) can be rewritten as a VAR(1) using the so-calledcompanion matrix. Thus, obtaining forecasts amounts to iterateˆ
YYY t + = F ( ˆ YYY t ) ≡ κκκ + Φ ˆ YYY t , to obtain ˆ YYY t + h = F h ( YYY t ) . (3)This equation provides forecasts of all variables, h periods from time t . An obvious t to con-sider is T , the end of the sample. The fact that we can obtain predictions by simply iteratingthe system, is of particular interest for generating scenarios for the Arctic. First, the predictionwill rely on an explainable mechanism – potentially mixing external forcing and internal feed-back loops – rather than a purely statistical relationship. Second, our forecast does not rely inany way on external data or forecasts made exogenously by some other entity which wouldrely on assumptions implicitly incompatible with ours. Nevertheless, in some cases, it may bedesirable to mix some external forecasts/scenarios of certain variables (like CO ) that may beless successfully characterized by the VAR and we do just that in section 4.4.The canonical macroeconomic VAR analysis seeks to explain business cycle fluctuations(expansions, recessions). Thus, many applications focus on modeling growth rates or de-viations from trend rather than levels. In such setups, information about the deterministiclong-run component has been suppressed, which implies that the VAR prediction (3) usuallyconverges to a constant for each variable as h → ∞ . This behavior is in line with standardmacroeconomic theory: an equilibrium path implies market forces balancing each other until This connection will be further discussed in sections 3.4 and 4.3. Further, it does not rely on the matrix A . In short, any VAR( P ) in M variables can be rewritten as a VAR(1) in M × P variables so the theoreticalanalysis can be carried out with the less burdensome VAR(1) (Kilian and Lütkepohl (2017)). YYY t are stacked yyy t − p ’sfor p =
1, ..., P . Complete lay out of Φ in equation (3) can be found in the appendix in equation A.3. While conditional and unconditional forecasting are important byproducts of our VARCTIC,another important objective of our analysis is to understand the underlying mechanism froma statistical standpoint. For instance, we will later show that anthropogenic CO is a keydriver of the long-run forecast (cutting emissions dramatically would prevent SIE from goingto 0). This important result rests solely on the reduced-form VAR. However, to uncover andinterpret the mechanism that amplifies CO ’s effect on SIE, we need an identification schemefor instantaneous relationships.Fundamentally, the typical time series identification problem originates from simultaneityin the data. That is, multivariate time series data can tell us whether X t − → Y t or Y t − → X t is more plausible. This is predictive causality in the sense of Granger (1969). However, thedata by itself cannot distinguish X t → Y t from X t ← Y t .This is a simple example of the simultaneity problem: one correlation between X t and Y t can begenerated by two different causal structures. Fortunately, the vast VAR literature has providedmany tools to address the identification problem.If one seeks to do structural analysis with the VAR, the above problem boils down to theneed for A in equation (1)). Yet, the data only procures us with the variance-covariance matrixof the residuals ˆ Σ . The identification problem emerges from the fact that A is not the only ma-trix that can satisfy ˆ Σ u = A − ′ A − . For instance, for any orthogonal rotation matrix H (whichhas the property of H = H − ) we can have ˜ A = H A and thus ˆ Σ u = ˜ A − ′ ˜ A − is also satisfied byconstruction. Mercifully, there exist many ways to pin down a single A matrix without havingto delve into too much theory, which is partially responsible for the popularity of VARs amongapplied economists. Furthermore, if the cross-equation correlation of residuals (off-diagonalelements of ˆ Σ u ) is small, then the impact on results of the chosen scheme will be rather limited. Much more on this can be found in Kilian and Lütkepohl (2017) Σ u , which implies acausal ordering of the data. It implies that residuals of a variable at position i are only consti-tuted of structural shocks εεε t (the unpredictable and uncorrelated disturbances of equation (1))of variables situated before it. Numerous more sophisticated alternatives have been proposedover the years (Kilian and Lütkepohl (2017)). The prime motivation for this is the occurrenceof "puzzles" (IRFs with counter-intuitive signs) especially when VARs are used to estimate theeffect of monetary policy on inflation. As figure 3 will show later, there are no such puzzlesin our application. Even better, we will find that by shuffling the ordering around (withinreasonable bounds), our results do not change in any significant way (see section A.4). T HE O RDERING . When using Cholesky decomposition to identify a VAR, the ordering of thevariables may influence the effect and transmission of shocks. Only if variable i is orderedbelow variable j , will a shock to j affect variable i contemporaneously. Otherwise, variable i will experience the effect of that very shock only with a lag.It is a relatively well established fact that the melting SIE and the responsible Arctic en-vironment are both results of exogenous (to other Arctic variable) human action (Dai et al.(2019), Notz and Stroeve (2016)). We view the Arctic system as being subject to feedback loopsthat may amplify the effect of exogenous shocks way beyond their original impact. However,the original stimulus is very likely to be anthropogenic, given that without the unprecedentedincrease in CO emissions and subsequent rise in global temperature, none of these mecha-nisms would have been so evident in effect (Amstrup et al. (2010), Melillo et al. (2014)). In2018, carbon dioxide accounted for 81.3% of greenhouse gases emitted in the United States.These greenhouse gases absorb infrared radiation, preventing the inherent heat content frombeing transmitted into outer space, triggering a response of global temperature (EPA, 2020)and initiating a cycle of knock-on effects on other factors of internal climate variability. Inour benchmark VARCTIC, CO is the only representative of the group of external forcing vari-ables. Therefore, we order CO first, meaning that shocks to any of the other variables mayonly impact CO with a one-period lag. In the spirit of many medium to large BVAR applications to macroeconomic data (Bernanke et al.(2005), Christiano et al. (1999), Stock and Watson (2005) and Ba ´nbura et al. (2010)), we classifythe variables, describing the internal climate variability, into fast-moving and slow-movingones. Cloud Cover, Precipitation and Air Temperature are classified as fast-moving. Absorb-ing short- and longwave radiation, clouds have a significant impact on theearth’s energy bal-ance and thus its overall heat content (Carslaw et al., 2002). But clouds, or the accumulation ofcondensed water and dust particles, eventually carry precipitation with not unambiguouslydeterminded effects on SIE (Parkinson and Comiso (2013), Meier et al. (2014)). We order bothvariables, Cloud Cover and Precipitation, before the temperature variables, i.e. Air Tempera- Meier et al. (2014) give an in-depth description of the various internal factors, their mutual interaction andtheir response to carbon dioxide. O N E XCLUDED P OTENTIAL M ECHANISMS . We consider VARCTIC 18 in part as a way toconfirm that key phenomena are already accounted for in VARCTIC 8. For example, studieshave emphasized the role of incoming long- and shortwave radiation and their interactionswith SIE and Thickness (see Burt et al. (2016), Dai et al. (2019)). The impact of downwellinglongwave radiation (DLW) on SIE is not direct, but transmitted via DLW’s influence on AT.Here, the thickness of sea ice is crucial, as thinner sea ice is more susceptible to DLW thanthicker layers (Park et al., 2015). As we will show later (like in figure 12), accounting for bothshort- and longwave radiation in VARCTIC 18, the forecast of an ice-free Arctic deviates onlymarginally from the ice-free date projected by the VARCTIC 8. This result suggests that short-and longwave radiation do not have a direct impact on SIE, but rather affect the evolutionof the Arctic’s sea ice cover via other variables (e.g. AT and Thickness), which VARCTIC 8already accounts for.In a similar line of thought, upper-ocean-heat content may also contribute to the evolutionof SIE. Examining the Barents Sea, Årthun et al. (2012) attribute the increased sea ice loss toan elevated influx of ocean heat. Studies have, however, shown that anomalies in the temper-ature of the upper-ocean layers and anomalies in SST do coincide during winter and springespecially over the Barents Sea (Park et al., 2015). We regard this as evidence for making anextension of both VARCTIC models dispensable.9 .4 Impulse Response Functions
Since Sims (1980), the dominant approach for studying the properties of the VAR around itsdeterministic path has been impulse response functions (IRFs) to structural shocks. Thanks tothe orthogonalization strategy discussed in 3.3, we converted plain regression residuals intoorthogonal shocks. The dynamic effect of these specific disturbances (the impulse) can beanalyzed as that of a randomly assigned treatment. IRFs take this structural meaning asa response to a fundamental shock when we rather look at uncorrelated impulses from ε m , t .Uncorrelatedness of the latter implies the "keeping everything else constant" interpretation –hence, a causal meaning for IRFs – is guaranteed by construction.As we just argued, a natural way of understanding the VAR is to look at its response toplausibly exogenous impulses called shocks . It is natural to wonder what is the meaning ofsuch shocks in a physical system. Mechanically, these shocks are the difference between therealized state of a variable and its predicted value as per the previous state of the dynamicsystem. These unpredictable anomalies, which emerge from outside a well-specified VARC-TIC, are the key to understand the dynamic properties of the model. A now obvious exampleof a shock will be that of CO emissions reduction in 2020: it is inevitable that the observedemissions will be lower than what was predicted by the endogenous system since the lat-ter excludes "pandemics". Any model that is partially incomplete will be subject to externalshocks. The study of such exogenous impulses may be alien-sounding, especially when con-trasted with the deterministic environment of a climate model. Nevertheless, understandingthe properties of a climate model by conditioning on a particular RCP scenario is equivalentto conditioning on a series of shocks. Hence, one can understand the VARCTIC and its IRFs asexpanding the number of potentially exogenous sources of forcing. Of course, our later focuson CO shocks is expressively motivated by the fact that the latter is a well-accepted source ofexogenous forcing in climate systems. The impulse response function of a variable m to a one standard deviation shock of ε ˜ m , t isdefined as I RF ( ˜ m → m , h ) = E ( y m , t | yyy t , ε t , ˜ m = σ ε ˜ m ) − E ( y m , t | yyy t , ε t , ˜ m = ) . Mathematically, we took a linear combination of the VAR residuals (an unpredictable change in a variable ofinterest, uuu t ) such that uuu t = A εεε ttt . Of course, one could look at how the system respond to an impulse from a residual u m , t , but the interpretationwill be rather weak because those are cross-correlated across equations. For instance, air temperature (AT) is largely determined by the previous values of variables in the system,but not completely. For macroeconomists, the idea of a single shock driving an otherwise completely endogenous system echoesback to the foundational work of Kydland and Prescott (1982) in which productivity shocks are considered as thesole driver of all aggregate macroeconomic fluctuations. That is, productivity is not predictable by the cyclicalsystem and driven by exogenous forces. Hence, in that setup, one would need an external productivity growthscenario to construct long-run forecasts the same way RCPs are needed in climate science. Lastly, it is worthmentioning that the VAR paradigm moves away from the original Kydland and Prescott (1982) setup by (amongmany other things) allowing for more than one shock – a necessary feature for economic models to be reconciledwith reality. CO increase and the same system where no such increase occurred. In a linearVAR with one lag ( P = all variables can easily be computed from the originalestimates using the formula I RF ( ˜ m → mmm , h ) = Ψ h A − e ˜ m where e ˜ m is vector with σ ε ˜ m in position ˜ m and zero elsewhere. This just means that we arelooking at the individual effect of ε ˜ m while all other structural disturbances are shut down. The latter discussion focused on analyzing how our dynamic system responds to an exter-nal/unpredictable impulse, which is a standard way of interpreting VAR systems. Of course,we are also interested in the "systematic" part of the VAR that is responsible for the propa-gation of shocks when they do occur – the IRF transmission mechanism. In section 4.3, wefocus our attention on CO and Air Temperature shocks and quantify the amplification effectof different channels. Over the years, many extensions to Sims (1980) original work have been proposed and somehave specific advantages that make them more adequate for our application. Precisely, we usea Bayesian VAR in the tradition of Litterman (1980). There are two crucial advantages of doingso. First, Bayesian inference does not change whether the VAR system is stationary or not(Fanchon and Wendel (1992)). We are effectively modeling variables in levels and expecting atleast one explosive root. Frequentist inference is notoriously complicated in such setups (Choi(2015)) and even standard approaches for non-stationary data have well-known robustnessproblems (Elliott (1998)). From a practical point of view, using non-stationary data (as we thinkis necessary here) means that standard test statistics (like popular Granger Causality tests) willbe undermined by faulty standard errors, potentially leading to erroneous conclusions.Second, for us to consider a system of many variables estimated with a relatively smallnumber of observations, Bayesian shrinkage can be beneficial to out-of-sample forecastingperformance and help in reducing estimation uncertainty (like those of IRFs). In fact, VARsare known to suffer from the curse of dimensionality as the number of parameters scales upvery fast with the number of endogenous variables. Via informative priors, Bayesian infer-ence provides a natural way to impose soft/stochastic constraints (that is, constraints are notimposed to bind) and yet keep inference (Ba ´nbura et al. (2010)). Furthermore, we are in-terested in transformations (forecasting paths, impulse response functions) of the parametersrather than the parameters themselves. Inference for such objects can easily and naturally be In the case of a linear VAR with P > Such a situation motivates McGraw and Barnes (2020) use of the LASSO. For instance, doing a VAR with LASSO would induce some form of shrinkage but inference is far from easy. The prior structure, however, must be chosen. We estimate our benchmarkBayesian VARCTIC with a standard Minnesota prior. In this simple framework, Σ , the variance-covariance matrix of the VAR residuals, is treated as known. Thus, the remaining param-eters of the model reduce to the vectorized matrix β = vec (cid:16)(cid:2) Φ · · · Φ ppp c (cid:3) ⊤ (cid:17) of dimension ( M p + M ) ×
1. The posterior distribution of β , π ( β | y ) , is obtained by the product of thelikelihood function of the data f ( y | β ) , and the prior distribution of β , π ( β ) . Hence, by sam-pling from the posterior distribution π ( β | y ) ∝ f ( y | β ) π ( β ) we can quantify both the uncertainty around β , but also more interesting transformationsof it, such as IRFs and forecasts. The prior distribution for β is the multivariate normaldistribution π ( β ) ∼ N ( β , Σ ) . The Minnesota prior is a specific structure for values of both β and Σ . An extended discussion the prior, its motivation for time series data and detailson the exact values of hyperparameters used can all be found in section A.3. Finally, we fixthe number of lags in the VARCTIC 8 to P =
12 and to P = A VAR contains many coefficients – there are 8 × ( × + ) =
776 in the baseline VARCTIC8. Staring at them directly is unproductive and a single coefficient (or even a specific block)carries little meaning by itself. As it is common with VARs in macroeconomics, we rather studythe properties of the VARCTIC by looking at its implied forecasts and its impulse responsefunctions. Setting priors’ tightness in such a way can be understood as analogous (at a philosophical level) to settingtuning parameters using cross-validation in Machine Learning. This choice is motivated by the fact that it facilitates the optimization of hyperparameters. As it turns out, op-timizing tuning parameters has more impact on resulting IRFs and their respective credible regions than treating Σ as unknown, when using for instance an Independent Normal Wishart (with Gibbs sampling). In the latter case, the credible region will naturally comprehend the uncertainty from the act of recursiveforecasting itself, but also the fact that it relies on unknown parameters that must be estimated. As a reference, a Ridge regression would imply β =
000 and Σ being a diagonal matrix with identical diagonalelements. The same arithmetic gives a total of 990 parameters in the VARCTIC 18. .1 The "Business as Usual" Forecast We report here the unconditional forecast of our main VAR. The VARCTIC 8 suggests SIE to hitthe zero lower bound around 2060 (see Figure 2), whereas in the VARCTIC 18 specification, theArctic would be ice-free about the same time (see Figure 12). The shaded area shows 90% ofall the potential paths of the respective VARCTIC. That is, the VARCTIC 8 dates the Arctic to betotally ice-free for the first time somewhere between 2052 and 2073 with a probability of 90%.The VARCTIC 18 slightly extends that time frame to the year 2079. For the two models, themedian scenario has SIE being less than one 10 km by 2054 and 2060 respectively. The one 10 km is more likely an interesting quantity since the "regions north of Greenland/Canada willretain some sea ice in the future even though the Arctic can be considered as ’nearly sea ice free’at the end of summer." (Wang and Overland (2009)). The corresponding credible regions markthe period 2047-2065 for the VARCTIC 8 and 2047-2069 for the VARCTIC 18 respectively. Thesedates and time spans range in the close neighborhood of previous climate model simulations(see Jahn et al. (2016)). For both VARCTICs, less than 5% of the simulated paths hit 0 before2050, making it an unlikely scenario according to our calculations. In essence, the two modelssuggest SIE melting at a rate that is slower than Diebold and Rudebusch (2019)’s results, butmuch faster than most CMIP5 models (Stroeve et al. (2012)) and in line with the latest CMIP6calculations (Notz et al. (2020)). Figure 3 shows the impulse response functions with the 90% credible region. Thus, the blueband reflects parameters’ uncertainty and contains 90% of the posterior draws from VARCTIC8. We display the response of SIE to a positive shock of one standard deviation to any of themodel’s M variables. Hence, when the credible region contains the 0 line, it means that morethan 5% of the posterior draws produce an IRF which sign is opposite to that of the poste-rior mean (the dark blue line). This implies that such an IRF does not describe a significantphenomenon, implying that the posterior probability of observing the opposite-signed effectis non-negligible.The resulting impact of CO anomalies on SIE is sizable and most importantly, durable.While the sign of the response is highly uncertain and weak for more than a year, CO shocksemerge to have a permanent downward effect on SIE. The relevance of the CO /SIE relation is We include in the graph the in-sample determinisic component of the VAR (as discussed in Giannone et al.(2019), which is essentially a long-run forecast starting from 1980 (the same sort of which we are doing right nowfor the next decades) using the VAR estimates of 12 lags. Diebold and Rudebusch (2019) stress the point that quadratic trends are likely to differ across months (es-pecially summer vs non-summer months). We accommodate for that using a refinement of the stochasticallyde-seasonalized series in section A.6. An interesting but more sophisticated extension – that would howeverbe beyond the scope of this paper – is to estimate a smooth-transition VARCTIC where dynamics are allowedto vary across seasons. In fact, many other such non-linearities/state-dependencies could be investigated andtested against our benchmark linear VARCTIC. CO impulses takes almost a year to settle in (not significant forapproximately 10 periods) but ends up having a permanent downward effect on trend SIEof approximately -0.005 10 km month after month. This mechanically implies that a one-off CO deviation from its predicted value/trend leads to a cumulative impact that is everincreasing in absolute terms (as displayed later in Figure 4b). It is important to remember thatthis is the effect of an unexpected increase in CO which is to be contrasted with the systematiceffect that will be studied later. However, in the framework of this section – where CO isallowed to endogenously respond to Arctic variables, this is as close as one can get to obtainan experimental/exogenous variation needed to evaluate a dynamic causal effect. -0.005 10 is roughly 0.1% of the last deterministic trend value of SIE, which is about the size of theGreat Salt Lake. CO shocks, by construction of our linear VAR, have mean 0 and there areapproximately as many positive and negative shocks in-sample. The linearity and symmetryof the VAR imply that these permanent effects are present for both upward and downwarddeviations from the deterministic trend.Other shocks have important impacts that eventually vanish, which is the traditional IRFshape one would expect to see from a VAR on macroeconomic data. For instance, Air Temper-ature and Albedo IRFs clearly have the expected sign. However, they do not have the strikinglasting impact of CO perturbations. In other words, a one-time Albedo shock will not have alasting effect on SIE as reported in Figure 2. This does not preclude Albedo to amplify othershocks as we will see in the next section.By looking at the yearly autocorrelation of SIE, Notz and Marotzke (2012) note that an ex-ceptionally low SIE is usually followed by a higher SIE the following year. Notz (2017) sur-14igure 3: IRFs: Response of Sea Ice Extentveys three main sources of negative feedback within a given year. The response of SIE to aSIE anomaly suggest indeed a small negative feedback – usually 6 months later. For instance,as Notz (2017) mentions, this could be the result of thinner ice that forms during winter beingable to grow much faster than thicker ice (that itself did not melt in the summer). Nevertheless,this does not preclude Thickness shocks from having an expected and mechanical negative ef-fect in the short run. Furthermore, as we will see later, the response of Thickness to CO shockswill itself contribute to the accelerated decay of SIE in the middle run.For a discussion of VARCTIC 18 results, see section A.5. We now turn to assess the validityof such an hypothesis by opening the black-box of the VAR transmission mechanism. We doso with IRF decompositions. CO and Temperature Shocks by Feedback Loops The melting of SIE is happening much faster than many other phenomena that are also be-lieved to be set in motion by the steady increase of CO emissions. Many recent papers(Notz and Marotzke (2012), Wang and Overland (2012), Serreze and Stroeve (2015), Notz (2017))argue with theory/climate models or correlations that external CO forcing is responsible forthe long-run trajectory of SIE. Some of these findings led Notz and Stroeve (2016) to concludethat climate models severely underestimate the impact of CO on SIE.15 rather consensual view is that the very nature of the Arctic system leads to the amplifi-cation of such external forcing shocks. For instance, the Albedo effect (less ice to reflect heat,more heat, less ice, repeat) has received a lot of attention in the literature. Another hypoth-esis is that of the thickness channel. As ice melts, older and thicker ice is replaced by newerand thinner ice. The latter is prone to melt faster when the next summer comes around. Thisthickness channel also has received some recent attention (Meier et al. (2014)). For instance,Parkinson and Comiso (2013) points out that the thinning of Arctic sea ice increases its sensi-tivity to abnormally adverse weather phenomena. However, things may not be that straight-forward as summarized in Notz (2017). The sum of feedbacks at the yearly frequency is infact likely to be negative . Univariate autocorrelation properties of SIE reveal an annual self-correction mechanism: a low September SIE in year t is usually followed by higher SIE thesame month next year. Hence, an understanding – from observational data – on how theArctic may amplify or not certain external forces is still pending. As Stroeve and Notz (2018)puts it:
It is difficult however, to robustly assess the contribution of internal variability to the ob-served loss, as this is only possible with climate models, which differ widely in their esti-mated magnitude of internal variability of the Arctic sea-ice cover.
In general, the VAR can quantify the contribution of different variables in explaining howa dynamic system responds to an external impulse. The VARCTIC encompasses differentamplification hypotheses can quantify which channels empirically matter and which ones donot. The only question left is how to convincingly extract that information from our model.
A potential approach that has a long history in econometrics is the use of Granger Causality(GC) tests. They have been recently advocated for climate applications by McGraw and Barnes(2018). Nevertheless, those tests only carry very limited information that quite often fall shortof answering questions of interest. First, the meaning of the test is not obvious when morethan two variables are included and if the researcher is interested in multi-horizon impacts, asdiscussed in Dufour and Renault (1998). Second, even if for some reason, rejecting the null ofa GC test is meaningful somehow, the block of reduced-form coefficients that we know to be ofsome statistical importance are very hard to interpret. For instance, when there are more than2 variables, the significant coefficients (by themselves) reveal close to nothing about indirecteffects. In sum, in the wake of a GC test rejection, we know some channel matters but we havelittle to no idea how it matters. In the light of all this, we choose another route that we believecould have wide applicability in empirical climate research, beyond the VARCTIC.When we observe that SIE is negatively affected by CO shocks, that response can be com-posed of a direct effect and many complicated indirect effects. Understanding indirect effects The negative autocorrelations are usually quite significant, so it is not a "regression to the mean" illusion.
16n the dynamic setup of a VAR is much more intricate than in a static regression setting. This isso because IRFs – for horizons greater than one – are obtained by iterating predictions, whichmeans X can impact Y through Z , but also through any of its lags.We employ a strategy that has been used in macroeconomics to better understand the trans-mission mechanism in VARs. It consists, rather simply, of shutting down "channels" and plot-ting what the response to a shock would be, given that channel had been shut. It is the VARequivalent of the partial (rather than total) derivative interpretation of ordinary least squaresregression coefficients. To further clarify, we proceed with examples from the literature thatfollowed Sims and Zha (2006)’s contribution. Bernanke et al. (1997) study the effect of oil priceshocks on the US economy, assuming the monetary policy authority had not responded in theusual way it has. This helps understanding whether oil prices themselves or subsequent (andsystematic) interest rate tightening is the real cause for a sequential decrease in economic activ-ity. Bachmann and Sims (2012) studies how an unexpected fiscal stimulus leads to increasedeconomic activity. They document that it is mainly through increasing "confidence" of eco-nomic actors that fiscal policy boosts GDP, especially in highly uncertain times like recessions.Running back to our Arctic concerns, we can deploy the same methodology to find and quan-tify the most important channels through which CO and temperature shocks impacts SIE. Inthe absence of a consensus for the name of this procedure, we will refer to it as IRF Decomposi-tion . CO Shocks
For the VARCTIC 8, Figure 4 shows the responses of SIE to an unexpected increase in onestandard-deviation of CO . The blue line pictures the case of the baseline VARCTIC 8 with90% credible region. The remaining 6 lines show the response of SIE to the same shock, butshutting down key transmission channels. In terms of implementation, it consists of imposing hypothetical shocks to one of the other variables which off-sets their own response to a CO shock. Figure 4a reveals – without great surprise – the importance of Temperature (especially AirTemperature) in translating CO anomalies into decreasing SIE. That is, we observe that shut-ting down these channels leads to a smaller absolute response which means that those vari-ables can be considered as amplification channels . Given the atypical shape of the CO IRF, thescale of Figure 4 makes less visible the action of channels that only alter the longer-run ef-fect. Since those effects are permanently negative (at different levels), their cumulative effectwill more clearly reveal their relative importance. Thus, Figure 4b displays the cumulativeimpact of selected (more important) channels. The two temperature channels are responsi-ble for approximately one fourth of the cumulative effect of CO on SIE after 3 years. Moreprecisely, restricting temperature variables to not respond to a positive CO shock, decreases See Bachmann and Sims (2012) for details. CO (a) IRF Decomposition (b) Cumulative IRF Decomposition (in absolute terms) the after-3-years impact from -0.13 10 km to -0.1 10 km . Of course, itwas expected that temperature should be a major conductor of such shocks. We also observesimilar quantitative effects for both thickness and the albedo effect in isolation. Most strik-ingly, we find out that the conjunction of the Albedo and Thickness amplification channels isresponsible for amplifying the effect of CO shocks by a non-negligible 50%. In the case ofthickness, this goes in line with the view that thinner ice has a harder time withstanding atmo-spheric forcings or abnormally warm ocean currents (Parkinson and Comiso (2013)), whichcan create a positive feedback loop. Furthermore, the amplification by Albedo also matchesevidence reported in several studies (see Perovich and Polashenski (2012), Björk et al. (2013),Parkinson and Comiso (2013)) that use very different methodologies.This IRF Decomposition exercise provides statistical evidence to suggest that (i) anthro-pogenic CO anomalies are a key driver behind the trajectory of SIE and (ii) Arctic feedbackloops (Albedo effect & the thinning of ice) eventually doubles CO initial impact on SIE. Thesefindings, obtained from a complete statistical model, broadly go in line with the emerging con-sensus that the Arctic evolution is driven by a combination of anthropogenic forcing variablesand the inherent dynamics of the Arctic itself (Stroeve and Notz (2018)). This section focusedon how and why SIE responds to CO shocks . In section 4.4, we rather look at the effect of the systematic increase of CO level. Air Temperature (AT) shocks are movements in AT that are not predictable given the past stateof the system and are orthogonal to other shocks in the system, most notably CO . In otherwords, we are looking at the effect of unexpected higher/lower AT that are uncorrelated withother shocks in the system. As we saw in Figure 3, such AT anomalies have a pronounced18hort-run effect on trend SIE for about a year after the shocks. This means that unlike CO , thecumulative effect of AT disturbances stabilizes about 1.5 year after the event.Figure 5: Transmission Mechanism Analysis - Shock to Air Temperature (a) IRF Decomposition (b) Cumulative IRF Decomposition In Figure 5a, we clearly observe (again) an important role for the thinning of ice and theAlbedo effect amplifying the response of SIE to AT shocks. In fact, we see in Figure 5b thatwithout them, the long-run impact is the same as the instantaneous one. Thus, this is evidenceto suggest that the AT shock’s long-run cumulative impact of -0.24 10 km is mostly a resultof the action of feedback loops. CO Emissions Scenarios If CO ’s trend is mostly or solely affected by factors outside of those considered in the VAR, theforecast of SIE can be improved by treating CO forcing as exogenous and using an externalforecast rather than the one internally generated by the VAR. Conditional forecasting can beachieved in VARs following the approach of Waggoner and Zha (1999). As we will see, thiswill markedly sharpen the bands around our forecasts, suggesting that a great amount ofuncertainty is related to the future path of CO emissions. It is also more common practice inclimate science to provide conditional rather than purely unconditional forecasts (Stroeve et al.(2012), Stroeve and Notz (2018)).The benefits of such an analysis are in fact twofold. The first, as we have seen, is obtainingpotentially better forecasts. The second is that when a specific scenario corresponds to a policychoice, we can evaluate the effect of such policies. In the spirit of Sigmond et al. (2018) whoconstrain the levels of AT in their climate model, we will look at CO emissions under three For the sake of clarity, the VARCTIC considered in this section treats CO as an exogenous variable for whichthe out-of-sample path is known. This is a natural approach given that long-run CO increase is undoubtlyanthropogenic. CO emissionsscenarios. For instance, Stroeve and Notz (2018) considers the evolution of SIE conditional onthose using a simple linear relationship between the two variables in levels and find that CO emissions can be a decisive factor between having the Arctic ice-free in the next 50 years ornot. Hence, it is of interest to see if the more complete VARCTIC 8 produces results in linewith their simpler statistical analysis. Additionally, in our case, we will be able to decomposesuch projections and assess the impact of internal variability in section 4.5.Figure 1 shows a steady increase in CO emissions over the last three decades, but sev-eral RCPs paint different pictures for the trajectory of carbon emissions until the end of thecentury. Therefore, we now consider the evolution of SIE under three different representa-tive concentration pathways: the RCP 2.6 , the
RCP 6 and the
RCP 8.5 . The
RCP 2.6 represents alow-mitigation scenario, in which CO emissions peak around mid-century (van Vuuren et al.,2011). Following this trajectory would be necessary to comply with the Paris Agreement(UNFCCC, 2015). The second scenario projects CO emissions – measured in gigatones ofcarbon per year – to peak during mid-century and taper off thereafter. This is very much in linewith levels of CO , projected by models in the absence of any climate-policies (van Vuuren et al.,2011). RCP 8.5 serves as the most pessimistic scenario. Figure 6a shows the different paths of CO under the three different RCP scenarios, as well as the projected path following VARCTIC8. Most interestingly, we find our completely endogenous and unconditional forecast of CO to lay somewhere between the "very bad" RCP 8.5 scenario and the "business-as-usual"
RCP 6 one.Figure 6b shows the variation of Arctic SIE in the VARCTIC when conditioning the out-of-sample path of CO on the three different RCP scenarios, as well as under the scenario of lettingthe model determine the future path of CO endogenously. The pictured effect is dramatic. Ifemissions were reduced as to follow the RCP 2.6 scenario, whose CO emissions are still at thehigher boundary of what the Paris Agreement demands, the Arctic would be far from blue andeven recover earlier losses by the end of the century. If emissions follow the more likely RCP 6 ,SIE would vanish later than projected by the VARCTIC 8, but still be completely gone duringthe 2070’s. In the worst case scenario,
RCP 8.5 , we obtain an ice-free September by the mid-2050’s. Interestingly, this result is very close to what Stroeve and Notz (2018) reported using avery different methodology (extrapolating a linear relationship). Their bivariate (SIE and CO )analysis suggests the Arctic summer months to be ice-free by 2050. However, in contrast, ourresults are much more optimistic than theirs in terms of SIE conditional on the (rather unlikely) RCP 2.6 scenario. Such analysis is not conditional on the identification scheme since it is based20igure 6: VARCTIC Projections & Different RCPs (a) Evolution of CO emissions until the End of the Century under different Scenarios(b) Evolution of SIE under different Scenarios of CO solely on the reduced form. Overall, these results reinforce the view that anthropogenic CO is the main driver behind the current melting of SIE as well as the main source of uncertaintyaround the future SIE path. Furthermore, the optimistic RCP 2.6 results suggest that internalvariability by itself cannot lead to the complete melting of SIE, even when starting from today’slevel. Overall, it is interesting to see that the VARCTIC yields similar conclusions about theimportance of CO to that of Dai et al. (2019) and Notz and Stroeve (2016). It is reassuring tosee that climate models conclusions can be corroborated by a transparent approach that relies Important to note is the fact that the very last in-sample observations for CO even range above the RCP 8.5 values, which generates the slight upward jump in case of the latter scenario. CO in Figure 1, one could worry that itmerely acts as a proxy for an omitted linear trend. We view the use of trends as undesirablein our multivariate setup as it would undermine the capacity of the VARCTIC to be a "com-plete" model. Including a trend would make it rely on a unknown/unexplained latent force– which is at odds with the main goal of our modeling strategy. Nevertheless, for the sakeof completeness, we estimated such models to find out that the inclusion of an exogenoustime trend is in fact not preferred by the data according to the Deviance Information Criterion(DIC, a generalization of the well-known AIC). The VARCTIC 8 has a DIC of -6894.35 andthe VARCTIC 8 with the exogenous trend has -6817.32. The smallest value being preferred,this justifies on a data-driven basis the exclusion of the trend. While seemingly technical ofpoint, this carries great meaning: the specification of the VARCTIC 8 system, based solely ondynamic relationships of observable data, can generate/simulate the observed SIE downwardpath. Additionally, the role CO as a central driving (downward) force is unmistakable inFigure 6b. The previous section documented the evolution of SIE conditional on several CO trajectories,treating it as an exogenous driver. This section seeks to quantify the importance of inter-nal variability when it comes to translating a RCP path into SIE loss. That is, we attempt toquantify to which extent the albedo effect and thickness effects can be held responsible foramplifying the effect of CO forcing and thus accelerating the melting of SIE.Following the findings of section 4.3, in which we identified Thickness and Albedo to showpotential for mitigating the adverse influence of CO on SIE, we ask the question about howSIE would evolve, if Thickness and Albedo were to remain constant at a certain level over theforecasting period. In particular, we repeat the forecasting exercise of the previous section forall three RCP scenarios, but keep Thickness and Albedo constant until the end of the forecast-ing period. For both variables we set the level equal to the value, which is given by the series’deterministic component at the end of the sample period. By doing so, we create artificialshocks to both Thickness and Albedo in each forecasting step, which off-set their response tothe external forcing variable. As we are modeling a dynamic and interconnected system, theseshocks do affect all the other variables (except for CO on which we condition our forecast).Figure 7 documents the corresponding results for RCP 8.5 , RCP 6 and
RCP 2.6 . For each sce-nario we show three different cases: (i) the projection of SIE under the respective RCP; (ii) theevolution of SIE under the respective RCP while keeping albedo constant at its last in-sample This approach makes more sense in macroeconomics where the trend is usually seen as an exogenous in-creasing productivity process which is the object of a different field of study (Growth). In the other words, thereis a strong belief that the trend and deviations from it arise from two very different models. This is not the casein our application to the Arctic where the trend is not a nuisance, but rather part of the essence. We view thisexercise more as of a way to illustrate the dynamic properties of the Arctic cryosphere asestimated by our VAR, in opposition to a definitive quantitative assessment of a potentiallyimplementable counterfactual.The role of internal variability as described by both Albedo and Thickness is undeniable.First, fixing Albedo to its 2019 value and thus shutting down this particular long-run amplifi-cation effect postpones the zero-SIE date by a bit less than a decade in both
RCP 8.5 and
RCP6 . Thickness of Arctic sea ice plays a major role for the reaction and resilience of SIE to anthro-pogenic forcing. Figure 7 re-enforces this view by showing that preventing both Thicknessand Albedo from further decay could in fact postpone the zero-SIE event to the next centuryunder
RCP 6 . Under
RCP 8.5 , shutting down both amplification channels starting from 2020leads to SIE crossing the bar of 1 × km about a decade later. Overall, this suggest amoderate contribution of both sources of feedback: shutting them down does not prevent SIEfrom heading towards 0 quickly. Nevertheless, this feeds into the pictured non-linearity andacceleration of SIE loss. In fact, for the two more realistic RCPs, shutting down both channelsimmediately makes the trend look more linear. Of course, those eventually accelerate, in ac-cordance with CO , but it is much later. Under the RCP 6 , we would still see a blue Arcticbefore the turn of the century, but the decrease flattens out in the very long-run. This providesa potential justification for the finding in Diebold and Rudebusch (2019) that a quadratic trendis a preferable approximation of long-run summer months’ SIE evolution. Finally, by symme-try of the VAR, we expect feedback loops that accelerate melting to also accelerate the reversephenomena after the 2050’s under
RCP 2.6 . However, this happens with a long delay, whichexplains why the curves constraining thickness and albedo are in fact above the green line forall the forecasting period. In that line of thought, fixing Thickness while letting SIE decrease is more likely to necessitate shocks of asize that has not been observed in our sample. The graphs are cut at the 1 × km bar as keeping Thickness constant (which the thought experimentsuggests) is untenable as SIE approaches 0: Thickness cannot be constrained to be positive if SIE is 0. (a) RCP 8.5(b) RCP 6(c) RCP 2.6 Conclusion and Directions for Future Research
We proposed the VARCTIC as a middle ground alternative to purely theoretical or statisticalmodeling. It generates long-run forecasts that embody the interaction of many key variableswithout the inevitable opacity of climate models. First, we focus our attention on exogenousimpulses which in the context of a structural VAR, have a meaning of quasi-experimental dis-turbances. Our results show that CO anomalies have a permanent effect on SIE which takesabout a year to settle in. It is the only impulse that has the property of permanently affectingSIE. Other impulses usually die out after a year and a half. We show that Albedo and Thick-ness play an important role in amplifying the response of SIE to CO and Air Temperatureshocks. In both cases, the conjunction of the two effects can double the cumulative impact ofsuch shocks after two years.Second, we focus on the systematic/deterministic part of the VARCTIC and conduct condi-tional forecasting experiments that again seek to quantify the effect of anthropogenic CO andhow internal variability can amplify it. We condition on the future path of CO and show that,within the context of our model, it is the prime source of uncertainty for the long-run forecastof SIE. RCP 8.5 implies 0 September SIE around 2054,
RCP 6 says so around 2075 and finally,
RCP 2.6 ( ∼ Paris Accord) implies that such an event would never happen. We conclude theanalysis by evaluating to which extent internal variability can amplify the long-run effect of CO forcing. Overall, our results provide statistical backing for the view that internal variabil-ity (as characterized here by Albedo and Thickness) can indeed transform relatively linearlytrending CO emissions into a non-linearly melting SIE.Regarding future research on the Arctic ecosystem, there are many methodological exten-sions available within the VAR paradigm that could be of some interest for climate scientists.For instance, there is a class of models called Smooth-Transition VAR (with a popular appli-cation in Auerbach and Gorodnichenko (2012)) that could be put to use to accommodate fordynamics (read VAR coefficients) evolving over the seasonal cycle. Moreover, there is an abun-dant literature that allows for structural breaks (immediate discrete change of parameters)and time-varying parameters (slow/smooth change). For instance, Screen and Deser (2019)remark the importance of changing weather phenomena that transition through decadal cy-cles, such as the pacific oscillation. Extensions that would allow parameters to evolve throughtime could evaluate the quantitative relevance of such phenomena. Furthermore, some recentattention (Chavas and Grainger (2019)) has been given to the potentially non-linear relation-ship between CO and SIE. Methods that blend time series econometrics, Machine Learningand abundant data of the like in Goulet Coulombe (2020) could reveal interesting insights oncomplex/time-varying relationships in the Arctic. In sum, the VAR methodology and timeseries econometrics still offer a rather unexploited potential for research on the Arctic.25 eferences Amstrup, S., DeWeaver, E., Douglas, D., Marcot, B., Durner, G., Bitz, C., and Bailey, D. (2010).Greenhouse gas mitigation can reduce sea-ice loss and increase polar bear persistence.
Na-ture , 468(37326):955–958.Auerbach, A. J. and Gorodnichenko, Y. (2012). Measuring the output responses to fiscal policy.
American Economic Journal: Economic Policy , 4(2):1–27.Bachmann, R. and Sims, E. R. (2012). Confidence and the transmission of government spend-ing shocks.
Journal of Monetary Economics , 59(3):235–249.Ba ´nbura, M., Giannone, D., and Reichlin, L. (2010). Large bayesian vector auto regressions.
Journal of applied Econometrics , 25(1):71–92.Bernanke, B. S., Boivin, J., and Eliasz, P. (2005). Measuring the effects of monetary policy: afactor-augmented vector autoregressive (favar) approach.
The Quarterly journal of economics ,120(1):387–422.Bernanke, B. S., Gertler, M., Watson, M., Sims, C. A., and Friedman, B. M. (1997). Systematicmonetary policy and the effects of oil price shocks.
Brookings papers on economic activity ,1997(1):91–157.Björk, G., Stranne, C., and Borenäs, K. (2013). The sensitivity of the arctic ocean sea icethickness and its dependence on the surface albedo parameterization.
Journal of Climate ,26(4):1355–1370.Burt, M., Randall, D., and Branson, M. (2016). Dark warming.
Journal of Climate , 29(2):705–719.Carslaw, K., Harrison, R., and Kirkby, J. (2002). Cosmic rays, clouds, and climate.298(5599):1732–1737.Chavas, J.-P. and Grainger, C. (2019). On the dynamic instability of arctic sea ice. npj Climateand Atmospheric Science , 2(1):1–7.Choi, I. (2015).
Almost all about unit roots: Foundations, developments, and applications . CambridgeUniversity Press.Christiano, L. J., Eichenbaum, M., and Evans, C. L. (1999). Monetary policy shocks: What havewe learned and to what end?
Handbook of macroeconomics , 1:65–148.Dai, A., Luo, D., Song, M., and Liu, J. (2019). Arctic amplification is caused by sea-ice lossunder increasing co2.
Nature Communications , 10(121).26iebold, F. X. and Rudebusch, G. (2019). Probability assessments of an ice-free arctic: Com-paring statistical and climate model projections.Dieppe, A., Legrand, R., and van Roye, B. (2016). The bear toolbox.
Working Paper Series ,(1934).Dufour, J.-M. and Renault, E. (1998). Short run and long run causality in time series: theory.
Econometrica , pages 1099–1125.Elliott, G. (1998). On the robustness of cointegration methods when regressors almost haveunit roots.
Econometrica , 66(1):149–158.EPA (2020). Draft inventory of u.s. greenhouse gas emissions and sinks 1990-2018.Fanchon, P. and Wendel, J. (1992). Estimating var models under non-stationarity and cointe-gration: alternative approaches for forecasting cattle prices.
Applied Economics , 24(2):207–217.Fetterer, F., Knowles, K., Meier, W. N., Savoie, M., and Windnagel, A. K. (2017). Sea ice index,version 3. https://doi.org/10.7265/N5K072F8 . Accessed: 2019-08-29.Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C., Dar-menov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C.,Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi,R., Merkova, D., Nielsen, J., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert,S. D., Sienkiewicz, M., and Zhao, B. (2017). The modern-era retrospective analysis for re-search and applications, version 2 (merra-2).
Journal of Climate , 30(14):5419–5454.Giannone, D., Lenza, M., and Primiceri, G. E. (2015). Prior selection for vector autoregressions.
Review of Economics and Statistics , 97(2):436–451.Giannone, D., Lenza, M., and Primiceri, G. E. (2019). Priors for the long run.
Journal of theAmerican Statistical Association , 114(526):565–580.Goulet Coulombe, P. (2020). The macroeconomy as a random forest.Granger, C. W. (1969). Investigating causal relations by econometric models and cross-spectralmethods.
Econometrica: Journal of the Econometric Society , pages 424–438.Harvey, A. (1990).
Forecasting, Structural Time Series Models, and the Kalman Filter , CambridgeUniversity Press.Harvey, A. and Koopman, S. (2014). Structural time series models.
Wiley StatsRef: StatisticsReference Online . 27arvey, A. C. and Todd, P. (1983). Forecasting economic time series with structural and box-jenkins models: A case study.
Journal of Business & Economic Statistics , 1(4):299–307.Hyndman, R. J. (2010). Forecasting with long seasonal periods.
Hyndsight blog .Jahn, A., Kay, J., Holland, M., and Hall, D. (2016). How predictable is the timing of a summerice-free arctic?
Geophysical Research Letters , 43(17):9113–9120.Kilian, L. and Lütkepohl, H. (2017).
Structural vector autoregressive analysis . Cambridge Uni-versity Press.Kydland, F. E. and Prescott, E. C. (1982). Time to build and aggregate fluctuations.
Economet-rica: Journal of the Econometric Society , pages 1345–1370.Litterman, R. B. (1980).
A Bayesian procedure for forecasting with vector autoregressions . MIT.McGraw, M. C. and Barnes, E. A. (2018). Memory matters: a case for granger causality inclimate variability studies.
Journal of Climate , 31(8):3289–3300.McGraw, M. C. and Barnes, E. A. (2020). New insights on subseasonal arctic–midlatitudecausal connections from a regularized regression model.
Journal of Climate , 33(1):213–228.Meier, W., Hovelsrud, G., van Oort, B., Key, J., Kovacs, K., Michel, C., Haas, C., Granskog, M.,Gerland, S., Perovich, D., Makshtas, A., and Reist, J. (2014). Arctic sea ice in transformation:A review of recent observed changes and impacts on biology and human activity.
Reviewsof Geophysics , 52(3):185–217.Melillo, J., Richmond, T., and Yohe, G. (2014). Climate change impacts in the united states: Thethird national climate assessment.
U.S. Global Change Research Program , (841).Notz, D. (2017). Arctic sea ice seasonal-to-decadal variability and long-term change.
PastGlobal Changes Magazine , 25:14–19.Notz, D., Dörr, J., Bailey, D., Blockley, E., Bushuk, M., Debernard, J., Dekker, E., DeRepentigny,P., Docquier, D., Fuˇckar, N., Fyfe, J., Jahn, A., Holland, M., Hunke, E., Iovino, D., Khos-ravi, N., Massonnet, F., Madec, G., O’Farrell, S., Petty, A., Rana, A., Roach, L., Rosenblum,E., Rousset, C., Semmler, T., Stroeve, J., Tremblay, B., Toyoda, T., Tsujino, H., and Vancop-penolle, M. (2020). Arctic sea ice in cmip6.
Geophysical Research Letters , page e2019GL086749.Notz, D. and Marotzke, J. (2012). Observations reveal external driver for arctic sea-ice retreat.
Geophysical Research Letters , 39(8).Notz, D. and Stroeve, J. (2016). Observed arctic sea-ice loss directly follows anthropogenic co2emission.
Science , 354(6313):747–750. 28ark, H.-S., Lee, S., Kosaka, Y., Son, S.-W., and Kim, S.-W. (2015). The impact of arctic winterinfrared radiation on early summer sea ice.
Journal of Climate , 28(15):6281–6296.Parkinson, C. and Comiso, J. (2013). On the 2012 record low arctic sea ice cover: Combinedimpact of preconditioning and an august storm.
Geophysical Research Letters , 40(7):1356–1361.Perovich, D. and Polashenski, C. (2012). Albedo evolution of seasonal arctic sea ice.
GeophysicalResearch Letters , 39(8).Årthun, M., Eldevik, T., Smedsrud, L., Skagseth, Ø., and Ingvaldsen, R. (2012). Quantifyingthe influence of atlantic heat on barents sea ice variability and retreat.
Journal of Climate ,25(13):4736–4743.Screen, J. and Deser, C. (2019). Pacific ocean variability influences the time of emergence of aseasonally ice-free arctic ocean.
Geophysical Research Letters , 46(4):2222–2231.Screen, J. A. and Simmonds, I. (2010). The central role of diminishing sea ice in recent arctictemperature amplification.
Nature , 464:1334–1337.Serreze, M. C. and Stroeve, J. (2015). Arctic sea ice trends, variability and implications forseasonal ice forecasting.
Philosophical Transactions of the Royal Society A: Mathematical, Physicaland Engineering Sciences , 373(2045):20140159.Sigmond, M., Fyfe, J., and Swart, N. (2018). Ice-free arctic projections under the paris agree-ment.
Nature Climate Change , 8(5):404–408.Sims, C. (2012). Statistical modeling of monetary policy and its effects.
American EconomicReview , 102(4):1187–1205.Sims, C. A. (1980). Macroeconomics and reality.
Econometrica: journal of the Econometric Society ,pages 1–48.Sims, C. A. and Zha, T. (2006). Does monetary policy generate recessions?
MacroeconomicDynamics , 10(2):231–272.Stock, J. H. and Watson, M. W. (2005). Implications of dynamic factor models for var analysis.Technical report, National Bureau of Economic Research.Stroeve, J., Kattsov, V., Barrett, A., Serreze, M., Pavlova, T., Holland, M., and Meier, W. (2012).Trends in arctic sea ice extent from cmip5, cmip3 and observations.
Geophysical ResearchLetters , 39(16).Stroeve, J. and Notz, D. (2018). Changing state of arctic sea ice across all seasons.
EnvironmentalResearch Letters , 13(10):103001. 29tuecker, M., Bitz, C., Armour, K., Proistosescu, C., Kang, S., Xie, S.-P., Kim, D., McGregor, S.,Zhang, W., Zhao, S., Cai, W., Dong, Y., and Jin, F.-F. (2018). Polar amplification dominatedby local forcing and feedbacks.
Nature Climate Change , 8:1076–1081.Taylor, K., Stouffer, R., and Meehl, G. (2012). An overview of cmip5 and the experiment design.
Bulletin of the American Meteorological Society , 93(4):485–498.UNFCCC (2015). Paris agreement.
United Nations Framework Convention on Climate Change ,pages 1–25.van Vuuren, D., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G.,Kram, T., Krey, V., Lamarque, J.-F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith,S. J., and Rose, S. (2011). The representative concentration pathways: An overview.
ClimaticChange , 109(1):5–31.Waggoner, D. and Zha, T. (1999). Conditional forecasts in dynamic multivariate models.
Reviewof Economics and Statistics , 81(4):639–651.Wang, L., Yuan, X., Ting, M., and Li, C. (2016). Predicting summer arctic sea ice concentrationintraseasonal variability using a vector autoregressive model.
Journal of Climate , 29(4):1529–1543.Wang, M. and Overland, J. E. (2009). A sea ice free summer arctic within 30 years?
Geophysicalresearch letters , 36(7).Wang, M. and Overland, J. E. (2012). A sea ice free summer arctic within 30 years: An updatefrom cmip5 models.
Geophysical Research Letters , 39(18).Winton, M. (2013).
Sea Ice–Albedo Feedback and Nonlinear Arctic Climate Change , pages 111–131.American Geophysical Union (AGU).Zhang, J. and Rothrock, D. (2003). Modeling global sea ice with a thickness and enthalpy dis-tribution model in generalized curvilinear coordinates.
Monthly Weather Review , 131(5):845–861. 30
Appendix
A.1 Data Details
Table 1: List of Variables
Abbreviation Description Data Source
Age Gridded monthly mean ofSea Ice Age nsidc.orgAirTemp Gridded monthly mean ofAir Temperature esrl.noaa.govAlbedo Gridded monthly mean ofSurface Albedo disc.gsfc.nasa.gov CO Global monthly mean of CO esrl.noaa.govLWGAB Gridded monthly mean ofSurface Absorbed Longwave Radiation disc.gsfc.nasa.govLWGEM Gridded monthly mean ofLongwave Flux Emitted from Surface disc.gsfc.nasa.govLWGNT Gridded monthly mean ofSurface Net Downward Longwave Flux disc.gsfc.nasa.govLWTUP Gridded monthly mean ofUpwelling Longwave Flux at TOA disc.gsfc.nasa.govPrecipCMAP Gridded monthly mean ofPrecipitation esrl.noaa.govSeaSurfTemp Median northern-hemispheric mean Sea-SurfaceTemperature anomaly (relative to 1961-1990) metoffice.gov.ukSIE Gridded monthly mean ofSea Ice Extent nsidc.orgSWGNT Gridded monthly mean ofSurface Net Downward Shortwave Flux disc.gsfc.nasa.govSWTNT Gridded monthly mean ofTOA Net Downward Shortwave Flux disc.gsfc.nasa.govTAUTOT Gridded monthly mean ofIn-Cloud Optical Thickness of All Clouds disc.gsfc.nasa.govThickness Gridded monthly mean ofSea Ice Thickness psc.apl.uw.eduTotCloudCover Gridded monthly mean ofTotal Cloud Cover esrl.noaa.govTS Gridded monthly mean ofSurface Skin Temperature disc.gsfc.nasa.gov A.2 Transmission Mechanism Analysis for a Shock to SIE
The purpose of the TMA analysis is to assess how the response of variable i to a shock onvariable j changes, if a third variable z were immune to the shock generated by variable j .Here we follow Sims (2012) by differentiating between the direct and indirect effect. The formeris variable i ’s own response to the shock hitting variable j . However, the shock also affectsvariable z , which itself transmits the shock further to variable i . This channel is the indirect effect of a shock to variable j on the response of variable i . Hence, it is the latter that willexplain the role of variable z within the transmission channel of a shock to j on i . To do so,31ims (2012) introduce artificial shocks to variable z , which offset its own response to a shock to j . These artificial shocks have two effects: (i) the IRF of variable z will be zeroed over the wholeIRF horizon; (ii) the indirect channel transmits the artificial shock onto variable i and allows toidentify the direct effect of j on i .This procedure requires the transformation of the structural VAR, given in equation (1) intothe reduced form VAR of equation (2), which reads as follows: yyy t = ccc + P ∑ p = A − Ψ p yyy t − p + A − εεε t , (A.1)where A − is the Choleski decomposition of matrix A in equation (1). This imposes thenecessary restrictions in order to identify the contemporaneous relationships of the variables.In particular, it assumes higher ordered variables to have an immediate effect on variables thatare ranked below, but not vice versa. As CO is ordered first in all of our models, an exogenousshock to carbon dioxide in period t will have an immediate effect on all of the other variables.The companion form of equation (A.1) is YYY t = ccc + Φ YYY t − + A − εεε t , (A.2)where YYY t = (cid:2) y t y t − · · · y t − p − (cid:3) ⊤ and the corresponding companion matrix is Φ = A − Ψ A − Ψ · · · · · · A − Ψ p I · · · I · · · · · · · · · I . (A.3)An equivalent way (to what laid out in section 3.4) of constructing IRFs, i.e. the responseof variable i to a structural shock on variable j over all horizons h =
0, ..., H , is to proceediteratively. Hence, for a given period h , the response of i to a shock hitting j is given by I RF ( j → i , h ) = e i Φ h A − • , j (A.4)where e i is a selection vector of dimension 1 × M with 1 at entry i and 0 otherwise. A − • , j elicitsthe j th column of A − . Following Sims (2012), switching off the indirect effect of a shock tovariable j on i via variable z amounts to I RF ( j → z , h ) = ∀ h =
0, ..., H . That requires the artificial shocks, ε z , h , to be calibrated such that the response of variable z to a shock to variable j is zero over the whole IRF period. Hence, at h = artificial shock ε z ,1 is32 z ,1 = − A − j ,1 A − z ,1 . (A.5)As these shocks are transmitted through time, the artficial response ε z , h has to account for allthe past shocks, ε z , h − , for any periods beyond h = ε z , h = − I RF ( j → z , h ) + ∑ h − h ′ = e z Φ h − h ′ A − • , j ε z , h ′ e z A − z . (A.6)The altered IRFs (that omits the transmission channel z ) for all the variables in the model to ashock to j is I RF − z ( j → iii , h ) = I RF ( j → iii , h ) + h ∑ h ′ = e z Φ h − h ′ A − • , j ε z , h ′ . (A.7)So far, we have reviewed how IRF decomposition works when one is interested in shuttingdown a single channel at a time. In contrast to Sims (2012), our VAR comprises more thanthree variables. Therefore, in some cases, it is desirable to shut-down not only one, but a group Z ∈ M \ { i , j } of indirect channels. To do so, equations (A.5) and (A.6) need to be generalized.At impact, the artificial response of variable z to a shock to j does not only have to offset the direct effect of j , but also the indirect effect of a shock j via the indirect effect of all the other artificial responses ( εεε + z + ,1 ) of those variables in Z which are ordered above z . This amounts tothe following extension of equation (A.5): ε z ,1 = − A − j ,1 A − z ,1 + ∑ m ∈ z + ε m ,1 A − z ,1 ! . (A.8)Also equation (A.6) has to be adjusted accordingly. However, at horizons h > artificial re-sponse ε z , h will not only have to offset the contemporaneous effects of z + , but also compensatefor the artificial responses of all other variables in Z over the period h ′ = · · · h − ε z , h = − I RF ( j → z , h ) + ∑ h − h ′ = e z Φ h − h ′ A − j ε z , h ′ + ∑ h − h ′ = ∑ n ∈ Z ε n , h ′ + ∑ m ∈ z + ε m ,1 e z A − z . (A.9)Equation (A.7) for the modified IRF ( I RF − z ( j → iii , h ) ) remains intact. z + denotes all those variables in Z which are ordered above z . .3 Bayesian Estimation Details The Minnesota prior is a specific structure for values of both β and Σ . In words, it allowsconcisely to parameterize heterogeneity in both the prior mean and variance. It consists ofthree major elements: the first one is about β and the last two concern Σ .1. For any equation y m , t with m =
1, ... M – where M is the total number of observed vari-ables in the VAR – all parameters are shrunk to 0 except for its first own lag y m , t − . The lat-ter is usually shrunk to a value b AR between 0.5 and 1. This can be interpreted as shrink-ing each VAR equation to the much simpler and parsimonious AR(1) process. Giventhe persistent nature of time-series data, this structure for β is much more appropriatethan that of Ridge regression (or LASSO), which shrinks all coefficients homogeneouslytowards 0.2. It is often observed in multivariate time series models that y m , t − → y m , t will be waystronger than almost any of the y ˜ m , t − → y m , t relationships. λ therefore calibrates theintensity of shrinking dynamic cross-correlations rather than auto-correlations.3. Distant partial lag relationships (say y m , t − → y m , t ) are expected to be of smaller magni-tude than close ones like y m , t − → y m , t , and y m , t − → y m , t . λ is in charge of determiningthe intensity of shrinking the coefficients of more distant lags.The overall tightness of the whole prior apparatus is determined by λ . Since we wishto be as agnostic as possible when tuning our model, we optimize/estimate hyperparameters { b AR , λ , λ , λ } within some grid to optimally balance bias and variance. To do so, we use themethodology developed in Giannone et al. (2015). This data-driven way of setting the priors’tightness can be understood as analogous (at a philosophical level) to setting tuning parame-ters using cross-validation in Machine Learning. Details on the exact values of hyperparame-ters used can be found in section A.3. Finally, we fix the number of lags in the VARCTIC 8 to P =
12 and to P = • Autoregressive Coefficient: = 0.9; • Overall tightness is λ = 0.3; • Cross-variable weighting is λ = 0.5; • Lag decay is λ = 1.5; • Exogenous variable tightness: λ = 100. For further details, explicit mathematical formulation of the prior and additional discussion on priors forVARs, the reader is referred to (Dieppe et al., 2016).
34t is worth remembering that the different λ ’s are prior variances . A larger value of λ and λ implies a looser prior, whereas a higher λ assigns less importance to lagged values. To draw aparallel to penalized regression (like Ridge and LASSO), a small λ in a Bayesian VAR increaseregularization in way analogous to how increasing the λ RIDGE – that is, by pushing the BVARestimate ˆ Φ away from ˆ Φ OLS .Put shortly, λ guides the overall level of shrinkage in the model. λ is indicative of howmuch the cross-variable dynamic relationship are shrunk to zero relative to own lags. It is oftenobserved in multivariate time series models that Y t − → Y t will be way stronger than muchof the X t − → Y t relationships. Thus, λ defines how we a priori think that autocorrelationsshould be more important to explain Y t than dynamic cross-correlations. Finally, λ is yetanother hyperparameter in charge of determining the tightness of a reasonable prior: far awaylag relationships (say Y t − → Y t ) are expected to be less important than close ones like Y t − → Y t and Y t − → Y t .In the following subsections, we report basic results – namely, the long run forecast andIRFs – with tighter and looser priors. In the latter case the point estimates approach whatwould have been obtained by classical Maximum Likelihood. Results remain both quali-tatively and quantitatively unchanged. Additionally, we experiment with alternative priorspecifications and again similar results. A.3.1 Altering the Priors
The hyperparamter λ , is tuned to be more relaxed than the ones of the benchmark VARCTIC8 in section 4, whereas the lag decay is strengthened and the AR coefficient is slightly reduced.Lags remain at 12. The present specification reads as follows: • Autoregressive Coefficient: = 0.8; • Overall tightness is λ = 1; • Cross-variable weighting is λ = 1; • Lag decay is λ = 3; • Exogenous variable tightness: λ = 100.In this specification, SIE is projected to hit the zero-lower bound in 2061. The DIC of -6828.57 is higher than its counterpart of the VARCTIC 8 (-6894.35).35igure 8: Trend Sea Ice ExtentFigure 9: IRFs: Response of Sea Ice Extent36 .4 Different Ordering In this section, we check the sensitivity of the responses of SIE to a shock of any of the othervariables when varying the ordering of variables compared to the benchmark VARCTIC 8in section 4. The priors and lags remain unaltered to the specification outlined in sectionA.3. The ordering now reads: CO , AirTemp, SeaSurfTemp, TotCloudCover, PrecipCMAP,SIE, Thickness, ALBEDO.Figure 10: IRFs: Response of Sea Ice ExtentA comparison of the responses of the benchmark VARCTIC 8 in Figure 3 and the IRFs afterreordering the model (Figure 10) documents the robustness of results to different identificationschemes. A second – more radical – variation in the model set-up locates Albedo at positiontwo. Hence, a shock to Albedo will contemporaneously affect all the other variables except CO : CO , ALBEDO, TotCloudCover, PrecipCMAP, AirTemp, SeaSurfTemp, SIE, Thickness.37igure 11: IRFs: Response of Sea Ice ExtentFor most of the effects, the shapes remain robust in comparison with Figure 3. Only theresponse to air temperatures deviates visibly with the statistically significant impact in theshort-run now vanishing. A.5 Results of VARCTIC 18
VARCTIC 18, including all the variables displayed in Figure 15, tests the robustness of theVARCTIC 8 projection of SIE. The ordering of variables in the VARCTIC 18 reads as follows:SWGNT, SWTNT, CO , LWGNT, TotCloudCover, TAUTOT, PrecipCMAP, TS, AirTemp, Sea-SurfTemp, LWGAB, LWTUP, LWGEM, SIE, Age, Thickness, EMIS and Albedo. Due to the in-creased number of variables, the lags were reduced to 3 and the estimation period starts in 1984due to some series unavailability. With more parameters to estimate, the prior-specificationslightly tightens and reads as follows: • Autoregressive Coefficient: = 0.8; • Overall tightness is λ = 0.5; • Cross-variable weighting is λ = 0.5; • Lag decay is λ = 3; 38 Exogenous variable tightness: λ = 100.The VARCTIC 18 predicts Arctic sea ice to reach the zero-lower bound by the year 2062,which is in the very neighborhood of the VARCTIC 8 (see Figure 12). This result suggeststhe VARCTIC 8 to comprise the key variables for a proper and robust projection of Arctic seaice. With the key-mechanisms for forecasting SIE aparently being captured by the benchmarkspecification, we conduct the investigation of the main feedback-channels and amplificationmechansims of sections 4.3 through 4.5 by using the VARCTIC 8.Figure 12: Trend Sea Ice ExtentFor completeness, the impulse response of SIE to a shock to any of the other variables isshown in Figure 13. IRFs of key variables remains roughly unchanged in VARCTIC 18. Mostinterestingly, in the VARCTIC 18, not only the CO shock has the effect of triggering a per-manently decreasing SIE, but also LWGAB, which measures the longwave radiation absorbedby the surface. Great many other shocks have statistically significant impacts in the shortrun. However, none has the lasting (and damaging) impact of CO and LWGAB. These lattertwo shocks have the unique property of permanently pushing the system out of the formerequilibrium.The forecast of SIE under the specification of the VARCTIC 18 is shown in Figure 14. Theprojected ice-free dates under the RCP 6 and
RCP 8.5 scenarios are consistent with the resultsreported by the VARCTIC 8 in Figure 6a. The trajectory of SIE under
RCP 2.6 , however, slightlychanges and seems to stabilize rather than recover by the end of the century.
A.6 Stochastic De-seasonalization
As discussed in section 2, the raw data (see 15) is highly seasonal except those that have beenpre-treated by the data provider. 39igure 13: IRFs: Response of Sea Ice ExtentFigure 14: Evolution of CO emissions until the End of the Century Different ScenariosVARCTIC 18
As a robustness check, we verify that our main results hold if we employ a radically differ-ent technique to take out seasonality. In this subsection, we adopt the approach of structural40igure 15: Raw Data: 18 Variables
CO2
TotCloudCover
TAUTOT
PrecipCMAP
SWGNT
SWTNT
LWGNT
AirTemp TS SeaSurfTemp
LWGAB
LWTUP
LWGEM
SIE
Thickness
Age
EMIS
ALBEDO time series (Harvey (1990) and Harvey and Koopman (2014)) where y raw is split in three some-what intuitive parts: y rawt = µ t + γ t + η t a trend component µ t ; a seasonality component γ t and a (possibly autocorrelated) noise com-ponent η t . Each of them are stochastic and have their own law of motion. The structure andlaw of motions we use follow the well-established Harvey Basic Structural Model (Harvey and Todd(1983)). The model reads as follows. µ t = µ t − + β t + u t β t = β t − + v t γ t = − ∑ m = γ t − m + w t ( η t , u t , v t , w t ) ∼ iid N ( Σ ) = σ ηη σ uu σ vv
00 0 0 σ ww The law of motion is that of Harvey and Todd (1983) and fits in a traditional state spacemodel. The trend µ t is a random walk with a stochastic drift. The drift β t is itself evolvingaccording to a random walk. For instance, this means that µ SIE , t , the trend of SIE, is trendingdown stochastically at a rate β SIE , t . That (negative) growth rate is itself allowed to evolve. Aquick look at a flexibly modeled trend of SIE suggests that allowing for a time-varying growthrate is necessary given the acceleration and deceleration of the melting of SIE in the 2000’s.Figure 16 shows the complete set of stochastic trends resulting from the BSM.Figure 16: Basic Structural Model: 8 Variables Extracted Trends adjusted for average September-Seasonality
CO2
TotCloudCover
PrecipCMAP
AirTemp
SeaSurfTemp
SIE
Thickness
ALBEDO
The extraction of trends as a first step and their subsequent modeling as a second stepis analogous to standard practice in macroeconomics, but not similar. In macroeconomics, itis customary in a strand of empirical work to filter the data as a pre-processing step. TheVAR is then estimated on the extracted cycles, which is simply the difference of the raw dataand the estimated trend. Here, we are indeed doing the filtering step first but using trendscomponents – rather than seasonality and short-run noise – for the second step. However, ourtrend components µ t are rather stochastic with respect to what is usually seen in economics. A.6.1 The Benchmark Specification and Results
Following Giannone et al. (2015), we obtain the optimal hyperparameters: • Autoregressive Coefficient: = 1; 42
Overall tightness is λ = 0.3; • Cross-variable weighting is λ = 0.5; • Lag decay is λ = 1.51; • Exogenous variable tightness: λ = 100;The date of the zero-lower bound of the stochastic de-seasonalized version remains in theneighborhood of the benchmark model. In this specification, the Arctic would be ice-free bythe year 2062. Figure 17: Trend Sea Ice Extent Stochastic De-seasonalization
Stochastic De-seasonalization
Figure 19: Evolution of SIE under different Scenarios of CO Different ScenariosStochastic De-seasonalization - Extracted trend adjusted for mean September-seasonality