Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
71,898 characters · 18 sections · 106 citation commands
Arctic Amplification of Anthropogenic Forcing: A Vector Autoregressive Analysis
\qquad \qquad Maximilian G\"obel}
\newgeometry{left=1.3cm, right = 1.3cm, bottom = 2.0cm, top = 2.0cm} \thispagestyle{empty}
\restoregeometry
\pagenumbering{arabic}
With 3.74 million square kilometers on September 15$^\text{th}$ 2020, Arctic sea ice extent ranked second-to-lowest in history, after the record minimum in 2012. A persistent retreat of SIE may further accelerate global warming and threaten the composition of the Arctic's ecosystem (ScSim10). The Coupled Model Intercomparison Project (CMIP) assembles estimates of long-run projections of Arctic sea ice from many climate models. These models try to reproduce the geophysical dynamics and interrelations among various variables, influencing the evolution of global climate.
With CMIP being in its 6$^{\text{th}}$ phase (CMIP6), climate models now provide more realistic forecasts of the Arctic's sea ice cover compared to its predecessor CMIP5 (see StroeveEtAl12, Notzetal2020). The majority of contributors to CMIP6 see the Arctic's September mean sea ice to retreat below the $1 \times 10^6$ $km^2$ mark before the year 2050. Despite following the hitherto accepted physical laws of our climate, its chaotic nature, i.e. the still obscure interplay of various climate variables, imposes a major burden on climate models. Repeated initialization with differing starting conditions is intended to reduce uncertainty and biases surrounding initial parameters. The byproduct is a wide range of projections of key climate variables Notzetal2020. In addition to such tuning, these simulations require large amounts of input data and a coupling of various sub-models TaylorEtAl09.
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 depicts the interaction of key variables describing the state of the Arctic. In such a setup, the downward SIE path will be an implication of a complete dynamic system based on the observed climate record. We provide a formal statistical assessment of different hypotheses about the historical path of SIE and outline the implications for the future. The effects on Arctic sea ice arising from various physical processes -- and the uncertainty surrounding their estimation -- can both be quantified without resorting to use a climate model. \vskip 0.2cm
{\sc Feedback Loops.} Feedbacks are dynamics initially triggered by an external shock to the system. Such a disturbance can be of radiative nature or not.\footnote{In contrast, internal variability, another source of climactic variation, describes fluctuations emerging from within the climate system KayEtAl2015.} Our analysis aims at better understanding how feedback loops -- and their interaction with anthropogenic carbon dioxide ($CO_2$) forcing -- shape the response of key Arctic variables, and most notably, sea ice.\footnote{A detailed description of various feedbacks, which the VARCTIC is capable of modeling, can be found in GoosseEtAl2018.} $CO_2$ forcing is already widely suspected to be the main driver behind long-run SIE evolution (see Meier2014, notz2017arctic). Feedback loops are well documented in the literature (see Park13, Winton13, StueEtAl18, McGraw19) and their understanding is crucial for enhancing the predictability of the Arctic's sea ice cover (Wang16, Notzetal2020). Only an approach that considers the interaction of many 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. CMIP6 models consider many variables, but high variation in sea ice projections (see Notzetal2020) suggests (among other things) widespread uncertainty around how strongly feedback loops may amplify external forcing. To shed more compelling statistical light on the matter, we borrow a methodology from economics.
\vskip 0.2cm
{\sc The Varctic.} Our analysis focuses on the evolution of the long-term trajectory of SIE and the interdependent processes behind it. The modeling approach, which we propose, achieves a desirable balance between purely statistical and theoretical/structural approaches. In many fields, statistical approaches often have a better forecasting record than theory-based models.\footnote{ When it comes to September Arctic sea ice, statistical approaches supplanted dynamical models for at least the last three years, as per the Sea Ice Prediction Network's Sea Ice Outlook post-season reports SIOAug2019. Statistical models showed much less disparity than theory-based alternatives and, most importantly, consistently provided a median forecast closer to the realized value.} An obvious drawback is that the successful statistical model may provide little to no explanation of the underlying physical processes.
A Vector Autoregression (VAR) lives in a useful middle ground. It is a statistical model that yet generates forecasts by iterating a complete system of difference equations in multiple endogenous variables. These interactions can be analyzed and provide an explanation for the resulting forecasts. Considering all this, we propose the VAR for the Arctic (VARCTIC), a statistical approach that (i) can generate long-run forecasts, (ii) can explain them as the result of feedback loops and external forcing (iii) allows us to analyze how the Arctic responds to exogenous impulses/anomalies.
\vskip 0.2cm
{\sc Roadmap.} We first discuss the data and its transformation in section (ref). We proceed with discussing the VAR model, its identification and Bayesian estimation in section (ref). Section (ref) contains the empirical results which comprise (i) a long-run forecast of SIE, (ii) impulse response functions of the VAR, (iii) an exploration of the transmission mechanism (feedback loops), and (iv) a conditional forecasting analysis. We conclude and propose directions for future research in section (ref).
Our data set comprises eighteen time series, proxying the Arctic's climate system, and accounting for potential feedback loops among the different constituents. The sample covers monthly observations from 1980 through 2018. We rely on standard data providers (see stroeve2018), which are listed in Table (ref) in the appendix. We combine eight variables, which importance has been highlighted by the existing literature (Meier2014), into VARCTIC 8, our benchmark specification. Fortunately, variables can easily be added/removed from a VAR. Bayesian shrinkage ensures that a larger model will not overfit -- the latter aspect is further explained in section (ref). Therefore, the appendix contains VARCTIC 18 which includes an additional 10 series from the reanalysis product MERRA2 (MERRA2) as a robustness check. To summarize compactly, the two specifications considered in this paper are:
A comprehensive overview of all variables (including those of VARCTIC 18), their acronyms, and links to data providers can be found in the appendix in Table (ref).\footnote{The primary goal was to assemble empirical data on key climate variables. To capture the most prominent feedbacks on SIE (see Meier2014), we augmented the observed series for $CO_2$, SIE, PR, and the assimilated PIOMAS product SIT, with data from model output. Our choice of data series is conditioned on whether they are (i) operated by well-established climate science institutions (ii) cited in the literature. } We want the VARCTIC to be a credible approximation of a completely endogenous system, where local processes jointly determine each other, without significant external dependencies outside of forcing.\footnote{This is precisely what allows us to iterate the system forward (see section (ref)(ref)) in order to obtain statistical forecasts based on a dynamic system.} Thus, we restrict the spatial coverage to a regional rather than a global scale. In line with the literature notz2016observed, all variables (except $CO_2$, {which is measured globally, and SST, which is measured over the Northern Hemisphere HorvathEtAl2020}) are monthly means over all grid-cells between 30$^{\circ}$N and 90$^{\circ}$N latitude. {This region matches the spatial coverage of the Sea Ice Index and is in the neighborhood of the lower bound of 40$^{\circ}$N latitude applied in HorvathEtAl2020.\footnote{Previous studies have emphasized the interdependencies between weather effects in the midlatitudes and the Arctic (mcgrawbarnes2018,ScreenEtAl2015).} It is a legitimate concern that averaging over too large of a region could wrongfully blend together mid-latitude events with others more specific to the Arctic circle. Fortunately, all key findings remain unchanged when restricting the gridded variables of TCC, PR, AT, and Albedo to the 60$^{\circ}$N-90$^{\circ}$N domain. An interesting avenue for future research is to consider a (larger) VAR where means over various latitude ranges are included -- so to study their interactions and relationship with SIE.} Further, we follow OelkeEtAl2003 and use AT measured at a sigma-level of 0.995, i.e. at 0.995 $times$ each grid-cell's surface-level pressure. For its part, the important choice of VARCTIC 8's variables themselves (and additions in VARCTIC 18) will be motivated extensively in section (ref).
The raw data is highly seasonal --- but the feedback loops we wish to estimate and extrapolate, reside in the (stochastic) trend components and short-run anomalies. Hence, we proceed to transform the data so that the resulting VARCTIC is fitted on deviations from seasonal means. For our benchmark analysis, we use a simple and transparent transformation: we de-seasonalize our data by regressing a particular variable $y^{raw}$ on a set of monthly dummies. That is, for each variable we run the regression
with $y_t$ being defined as $y_t \equiv y_{t}^{raw} - \sum_{m=1}^{12}\hat{\alpha}_m D_m$. $D_m$ is an indicator that is 1 if date $t$ is in month $m$ and 0 otherwise. The estimates of $\alpha_m$, $\hat{\alpha}_m$, are obtained by ordinary least squares. This is exactly equivalent to de-meaning each data series month by month and is a more flexible approach to modeling seasonality than using Fourier series.\footnote{Of course, if we were using higher-frequency data -- like daily observations, then the Fourier approach would be much more parsimonious and potentially preferable hyndman2010forecasting. The dummies approach to taking out seasonality only requires 12 coefficients with monthly, but 365 with hypothetical daily data.} Finally, we keep our filtered data $y$ in levels. We do not want to employ first differences or growth rate transformations to make the data stationary. Such an action would suppress long-run relationships which are an important object of interest. Figure (ref) in the Appendix shows the data after being filtered with monthly dummies.\footnote{Note that $CO_2$ is available without seasonality from the data provider (NOAA-ESRL) and thus was not passed through the dummies filter.}
\textcolor{black}{Pre-processing the data can influence results. Moreover, DieRu19 and Meier2014 document seasonal variability in SIE trends. As a natural robustness check, we also consider a very different approach to eliminate seasonality. In appendix (ref), we reproduce our results with a data set of stochastically de-seasonalized variables obtained from the approach of structural time series (Harvey90 and harvey2014). In short, this extension allows for seasonality to evolve (slowly) over time, which could be a feature of some Arctic time series.}
In this section, we review the VAR: the model; its identification; its Bayesian estimation. Furthermore, we discuss the construction of the long-run forecasts and impulse response functions as tools to understand the VARCTICs' results.
Vector Autoregressions are dynamic simultaneous systems of equations. They can characterize a linear dynamic system in discrete time. The methodology was introduced to macroeconomics by sims1980 and is now so widely used that it almost became a field of its own (see kilian2017svar). It is a multivariate model in the sense that $\boldsymbol{y}_t$ in
is an $M$ by 1 vector. This means that the dynamic system incorporates $M$ variables. $\Psi_p$'s parameterizes how each of these variables is predicted by its own lags and lags of the $M-1$ remaining variables. The matrix $A$ characterizes how the $M$ different variables interact contemporaneously --- e.g., how AT affects SIE within the same month (a time unit $t$ in our setup). Finally, the disturbances are mutually uncorrelated with mean zero: \[ \boldsymbol{\varepsilon}_t = \left[ \varepsilon_{1,t},\enskip ... \enskip , \varepsilon_{M,t} \right] ~ \sim ~ \, N \left (0 , ~ I_M \right ). \] Equation ((ref)) is the so-called structural form of the VAR, which cannot be estimated because $A$ is not identified by the data. For clarity, the elements of $A$ are not plain regression coefficients, but structural model parameters. Attempting to estimate those directly via least squares would be plagued by simultaneity bias kilian2017svar. Rather, structural VAR estimation proceeds in two steps. First, an estimable "reduced-form" VAR is fitted to the data. That is, we run
where $\boldsymbol{c} = A^{-1} \Psi_0$ and $\Phi_p = A^{-1} \Psi_p$ are both regression coefficients. $\boldsymbol{u}_t$ are now regression residuals \[ \boldsymbol{u}_t = \left[ u_{1,t},\enskip ... \enskip , u_{M,t} \right] ~ \sim ~ \, N \left (0 , ~ \Sigma_u \right ) \] which are allowed to be cross-correlated. By construction, $\Sigma_u=A^{-1'}A^{-1}$. This last relationship is key to the second, so-called "identification", step. In words, this means the covariance matrix of regression residuals from the first step ($\Sigma_u$) can be used as raw material to retrieve the "structural" $A$ --- the latter which, as we stressed earlier, cannot be estimated directly. The process for obtaining $A$ by decomposing $\Sigma_u$ is addressed on its own in section (ref).
The methodology has many advantages over simple autoregressive distributed lags (ARDL) regression that have gained some popularity in the econometric and climate literature. For instance, in mcgrawbarnes2018, the argument for inclusion of lags of the dependent variable can be interpreted as one for completeness of the modeled dynamic system, as guaranteed by an adequately specified VAR.
The symmetry of the VAR allows for it to generate forecasts by simply iterating the model.\footnote{Further, forecasting does not rely on the matrix $A$.} Assuming the chosen variables to characterize the system completely, we can forecast its future state by iterating a particular mapping. To do so, we use a representation that exploits the fact that any VAR($P$) (that is, with $P$ lags) can be rewritten as a VAR(1), using the so-called companion matrix.\footnote{In short, any VAR($P$) in $M$ variables can be rewritten as a VAR(1) in $M \times P$ variables, such that the theoretical analysis can be carried out with the less burdensome VAR(1) kilian2017svar. $\boldsymbol{Y}_t$ are stacked $\boldsymbol{y}_{t-p}$'s for $p=1,...,P$.} Thus, obtaining forecasts amounts to iterate
where $F$ is the one-month ahead forecasting function, while $\boldsymbol{\kappa}$ and $\Phi$ are the companion-form analogs of $c$ and $\Phi_p$'s in (ref). This equation provides forecasts of all variables, $h$ periods from time $t$. An obvious $t$ to consider is $T$, the end of the sample. The fact that we can obtain predictions by simply iterating the system, is of interest to generate scenarios for the Arctic. First, the prediction will rely on an explainable mechanism -- potentially mixing external forcing and internal feedback loops -- rather than a purely statistical relationship. Second, our forecast does not rely in any way on external data or forecasts made exogenously by some other entity, which would rely on assumptions implicitly incompatible with ours. Nevertheless, in some cases, it may be desirable to mix some external forecasts/scenarios of certain variables (like $CO_2$) that may be less successfully characterized by the VAR. We consider just that in section (ref).
While conditional and unconditional forecasting are important byproducts of the VARCTIC, another objective of our analysis is to understand -- from a statistical standpoint -- the underlying process driving interactions between key Arctic variables. For instance, by forecasting SIE conditional on various emission scenarios, we will later show that anthropogenic $CO_2$ forcing is the main driving force behind the long-run forecast --- cutting emissions dramatically would prevent SIE from going to 0.\footnote{In contrast, an unconditional forecast lets the VARCTIC generate internally future paths for all variables (including $CO_2$) without relying on externally developed scenarios.} This important result rests solely on the reduced-form VAR. However, to uncover and interpret the mechanism that amplifies $CO_2$'s effect on SIE, we need an identification scheme for instantaneous relationships. In time series analysis, the identification problem originates from simultaneity in the data. Multivariate time series data can tell us whether $X_{t-1} \rightarrow Y_{t} \quad \text{or} \quad Y_{t-1} \rightarrow X_{t}$ is more plausible. This is predictive causality in the sense of granger1969. However, the data by itself cannot distinguish $X_{t} \rightarrow Y_{t} \quad \text{from} \quad X_{t} \leftarrow Y_{t}.$ In words, a single correlation between $X_t$ and $Y_t$ can be generated by two different causal structures. Within the VAR, the problem boils down to the need for identifying $A$ in equation ((ref)). Yet, the data only procures us with the variance-covariance matrix of the residuals $\hat{\Sigma}$. The identification problem emerges from the fact that $A$ is not the only matrix satisfying $\hat{\Sigma}_u=A^{-1'}A^{-1}$. Fortunately, there exist many ways to pin down a single $A$ matrix without having to delve into too much theory, which is partially responsible for the popularity of VARs among applied economists. The strategy we opt for is the traditional Choleski decomposition of $\hat{\Sigma}_u$. Mechanically, it provides a lower-triangular matrix $C$, satisfying $\hat{\Sigma}_u=C'C$ (where $C\equiv A^{-1}$ for convenience). Its purpose is to transform regressions residuals $u_t$ (equation (ref)) into uncorrelated structural shocks $\boldsymbol{\varepsilon}_t$ (equation (ref)). This is done by reversing the relationship $\boldsymbol{u}_t=C\boldsymbol{\varepsilon}_t$. Uncorrelatedness is essential (as further discussed in section (ref)) to study how the VARCTIC responds to a given impulse, keeping everything else constant. Such a causal claim would be impossible when considering an impulse from correlated residuals $u_t$ as those always co-move. In other words, studying $u_t$ assuming everything else stays constant is generally inconsistent with the data. In sum, the Choleski decompostion is one way to transform the observed (but practically useless) $u_t$ into the very useful (but originally unobserved) fundamental shocks $\boldsymbol{\varepsilon}_t$.
The assumption underlying such an approach to orthogonalization is a causal ordering of shocks. First, it is worth cataloging the relationships, i.e., which get restricted by the ordering choice and which do not. The dynamics (lead-lag relationships as characterized by $\Psi$) of the VAR are exempted as they are already completely identified by the data itself. Rather, the ordering restricts how variables interact together within the same month, conditional on the previous state of the system. This is done by making an explicit assumption about the composition of (reduced-form) deviations of Arctic variables from their predicted values (i.e., the anomalies). Precisely, the lower-triangular structure of (ref) implies that residuals of a variable at position $i$ are only constituted of structural shocks $\boldsymbol{\varepsilon}_t$ of variables ordered before it. To make that explicit, we report $\boldsymbol{u}_t=C\boldsymbol{\varepsilon_t}$ in full:
Only if variable $i$ is ordered below variable $j$, will a "fundamental" shock to $j$ affect variable $i$ contemporaneously. Otherwise, variable $i$ will experience the effect of that shock with a lag of at least one month (which corresponds to one time unit in the application). For example, the $CO_2$ anomalies (which means, unpredictable by the past behavior of any of the eight variables) are assumed to be composed of structural $CO_2$ shocks only. This implies that the effect of other variables on $CO_2$ take at least a month (but perhaps more) to set in. In contrast, SIE or Albedo anomalies can be composed of a variety of fundamental shocks. Those restrictions are not without cost as the ordering of the variables may influence our understanding of the mechanism uncovered by the VARCTIC. This is why the ordering must be motivated based on the application at hand.\footnote{Moreover, when possible, the robustness of results to some reasonable alterations of the ordering should be assessed.} \vskip 0.2cm
{\sc Motivating The Ordering.} It is well established that the melting SIE and the state of the Arctic environment are both results of exogenous (to other Arctic variables) human action (DaiLuoSong19, notz2016observed). We view the Arctic system as being subject to feedback loops that 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 unprecedented increase in $CO_2$ emissions and subsequent rise in global temperature, none of these mechanisms would have been so evident in effect (Ams2010, MelRichYo14).\footnote{ Meier2014 give an in-depth description of the various internal factors, their mutual interaction and their response to carbon dioxide.} Consequently, we order $CO_2$ first. The implication is that shocks to any of the other variables can impact $CO_2$ with a minimal delay of one month. In contrast, $CO_2$ can impact any variable in the system either contemporaneously, in the short/medium/long run, or both.
In the spirit of many medium to large BVAR applications to macroeconomic data (bernanke2005measuring, christiano1999monetary, stock2005implications and banbura2010large), we classify the variables, describing the internal climate variability, into fast-moving and slow-moving ones. TCC, PR and AT are classified as fast-moving. Absorbing short- and longwave radiation, clouds have a significant impact on the earth's energy balance and thus its overall heat content Carslaw2002. But clouds eventually carry precipitation with not unambiguously determined effects on SIE (Park13, Meier2014). We order both variables before the temperature variables AT and SST. Besides AT, also SST, especially warmer water from the Atlantic Ocean, contributed to shaping the historically unprecedented decline of SIE over the last four decades Meier2014. Here we follow Park13 who state that besides the cooling effects of a melting ice cover, SST is highly influenced by currents and winds, transferring warmer energy from lower to higher latitudes. We therefore place SST at the boundary of fast- and slow-moving variables.
The last block of variables comprises, SIE, SIT, and Albedo. SIT is an underestimated determinant of how SIE reacts to both external forcing and internal variability (Meier2014, Park13). Thicker layers make the ice more resilient and increase Albedo, while thin ice is more easily advected by winds, making SIE more sensitive to extreme events Meier2014. We order SIT -- and Albedo -- after SIE because we hypothesize that the effect of shocks of the former can only influence the latter with a certain delay. For instance, shocks to SIT via increased water precipitation or strong winds will immediately reduce SIT but SIE only with a certain lag. Lastly, we regard Albedo as being driven contemporaneously by all other factors.
To wrap up, it is worth re-emphasizing that identification, via the described ordering, is necessary to interpret and understand the mechanisms captured by the VARCTIC. However, ordering choices do not alter forecasts. Mechanically, this happens because the potentially contentious matrix $A$ does not enter the forecasting equation (ref).
{\sc On Potentially Excluded Mechanisms.} We consider VARCTIC 18 in part to confirm that the key channels are already accounted for in VARCTIC 8. For example, studies have emphasized the role of incoming long- and shortwave radiation and their interactions with SIE and SIT (see BurtEtAl2016, DaiLuoSong19). The impact of downwelling longwave radiation (DLW) on SIE is not direct, but transmitted via DLW's influence on AT. Here, thickness of sea ice is crucial, as thinner ice is more susceptible to DLW than thicker layers ParkEtAl2015. As we will show later (like in figure (ref)), accounting for both short- and longwave radiation in VARCTIC 18, the forecast of an ice-free Arctic deviates only marginally from the ice-free date projected by VARCTIC 8. This result suggests that short- and longwave radiation does not impact SIE directly, but rather affects the evolution of the Arctic's sea ice cover via other variables (e.g., AT and SIT), which VARCTIC 8 already accounts for. In a similar line of thought, upper-ocean heat-content may also influence to the evolution of SIE. Studies have found that anomalies in the temperature of the upper-ocean layers and anomalies in SST do coincide ParkEtAl2015, making an extension of both VARCTIC models dispensable.
However, it is not excluded that some non-local processes (e.g., poleward atmospheric heat transport) do contribute to sea-ice loss through channels not represented in both VARCTICs. As stated earlier, we opted for including local processes only (in addition to $CO_2$) because this makes the VARCTIC a complete system where all $M$ variables are internally modeled and forecasted jointly. Adding non-local processes raises the additional question of how to model their external dependence, a complication left for future research.\footnote{The literature on Global VARs could provide a natural place to start pesaran2009forecasting.}
Since sims1980, the dominant approach for studying the properties of the VAR around its deterministic path has been impulse response functions (IRFs) to structural shocks. Thanks to the orthogonalization strategy discussed in (ref), we converted plain regression residuals into orthogonal shocks.\footnote{Mathematically, we took a linear combination of the VAR residuals (an unpredictable change in a variable of interest, $\boldsymbol{u}_t$) such that $\boldsymbol{u}_t=C\boldsymbol{\varepsilon_t}$.} The dynamic effect of these specific disturbances (the impulse) can be analyzed as that of a randomly assigned treatment.\footnote{Of course, one could look at how the system responds to an impulse from a residual $u_{m,t}$, but the interpretation will be rather weak because those are correlated across equations.} Uncorrelatedness of $\varepsilon_{m,t}$ implies the "keeping everything else constant" interpretation -- hence, a causal meaning for IRFs -- is guaranteed by construction.
It is natural to wonder about the meaning of uncorrelated shocks in a physical system. Mechanically, these shocks are the difference between the realized state of a variable and its predicted value as per the previous state of the dynamic system. These unpredictable anomalies, which emerge from outside a well-specified VARCTIC, are the key to understanding the dynamic properties of the model. A now obvious example of a shock will be that of $CO_2$ emissions reduction in 2020: it is inevitable that the observed emissions will be lower than what was predicted by the endogenous system since the latter excludes "pandemics". Any model that is partially incomplete will be subject to external shocks. The study of such exogenous impulses may be alien-sounding, especially when contrasted with the deterministic environment of a climate model. Nonetheless, understanding the properties of a climate model by conditioning on a particular RCP scenario is equivalent to conditioning on a series of shocks. Hence, one can understand the VARCTIC and its IRFs as expanding the number of potentially exogenous sources of forcing. Of course, our later focus on $CO_2$ shocks is expressively motivated by the fact that the latter is a well-accepted source of exogenous forcing in climate systems.
The impulse response function of a variable $m$ to a one standard deviation shock of $\varepsilon_{\tilde{m},t}$ is defined as
Thus, it is the expected difference, $h$ months after "impact", between an Arctic system that responded to an unexpected $CO_2$ increase, and the same system where no such increase occurred. In a linear VAR with one lag ($P=1$), the IRF of all variables can easily be computed from the original estimates using the formula
where $e_{\tilde{m}}$ is a vector with $\sigma_{\varepsilon_{\tilde{m}}}$ in position $\tilde{m}$ and zero elsewhere. This just means that we are looking at the individual effect of $\varepsilon_{\tilde{m}}$ while all other structural disturbances are shut down.\footnote{In the case of a linear VAR with $P>1$ lags, we must use the companion matrix form. \textcolor{black}{The relevant formula (equation (ref)) can be found in the discussion of appendix (ref)}. }
The latter discussion focused on analyzing how our dynamic system responds to an external/unforeseeable 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 propagation of shocks when they do occur -- the IRF transmission mechanism. In section (ref), we focus our attention on $CO_2$ and AT shocks and quantify the amplification effect of different channels.
We use a Bayesian VAR in the tradition of litterman1980. There are two crucial advantages of doing so. First, Bayesian inference does not depend on whether the VAR system is stationary or not (FanWen92). We are effectively modeling variables in levels and expecting at least one explosive root. Frequentist inference is notoriously complicated in such setups (almosteverythingaboutUR) and even standard approaches for non-stationary data have well-known robustness problems (elliott1998). From a practical point of view, using non-stationary data means that standard test statistics (like Granger Causality tests) will be undermined by faulty standard errors, potentially leading to erroneous conclusions.
Second, for us to consider a system of many variables estimated with a relatively small number of observations, Bayesian shrinkage can be beneficial to out-of-sample forecasting performance and help in reducing estimation uncertainty (like those of IRFs). In fact, VARs are known to suffer from the curse of dimensionality as the number of parameters scales up very fast with the number of endogenous variables.\footnote{Such a situation motivates McGraw19 to use the LASSO.} Via informative priors, Bayesian inference provides a natural way to impose soft/stochastic constraints (that is, constraints are not imposed to bind) and yet keep inference possible banbura2010large.\footnote{For instance, running a VAR with LASSO would induce some form of shrinkage but inference is far from easy.} Furthermore, we are interested in transformations (forecasting paths, impulse response functions) of the parameters rather than the parameters themselves. Inference for such objects can easily be obtained by transforming draws from the posterior distribution. All these procedures are well established in the macroeconometrics community and packages are available in most statistical programming software BearMatlab. An extended discussion of the prior, its motivation for time series data and details on the exact values of (data-driven) hyperparameters used, can all be found in section (ref).
Finally, the maximal lag order of the VAR, $P$ in equations (ref) and (ref), must be chosen.\footnote{To re-emphasize, this means $y_{t-p}$ for $p=1,...,P$ are included.} Its selection is yet another incarnation of the bias-variance trade-off. We fix the number of lags in VARCTIC 8 to $P=12$ and to $P=3$ in VARCTIC 18 respectively. That choice is based on the Deviance Information Criterion (DIC) -- the Bayesian analog to popular information criteria used for model selection. Accordingly, the superior VARCTIC 8 would set $P=3$ (DIC=-6988)\footnote{The lower, the better.}, a choice which only provides a marginal improvement with respect to $P=12$ (DIC=-6894).\footnote{Additionally, the reported DIC for $P=12$ is superior to other natural candidates such as $P=1$ and $P=24$.} Since structural analysis is an essential part of this paper, we err on the side of having slightly higher variance, but potentially richer dynamics for IRFs. In the large VARCTIC 18, the need for shrinkage is magnified and $P=3$ is the obvious more reasonable choice.
An extraneous question, which can benefit from verification by DIC, is whether trends should be included. We hypothesize that VARCTIC 8 is a complete, divergent system which can endogenously explain the trending behavior of all its variables by the joint action of $CO_2$ forcing and feedback loops. If that were not to be true, including linear trends would noticeably improve model fit, and lower the DIC even further. Backing our claim that the VARCTIC needs no additional (and hardly climatically-explainable) statistical crutch, the DIC from including trends is worse (now DIC=-6817 for VARCTIC 8) than that of the original model.
A VAR contains many coefficients -- there are $8 \times (8 \times 12 +1)= 776$ in the baseline VARCTIC.\footnote{The same arithmetic gives a total of 990 parameters in VARCTIC 18.} 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 study the properties of the VARCTIC by looking at its implied forecasts and its IRFs.
We report here the unconditional forecast of our main VARs. VARCTIC 8 suggests SIE to hit the zero lower bound around 2060 (see Figure (ref)), whereas VARCTIC 18 projects the Arctic to be ice-free at about the same time (see Figure (ref)).\footnote{We include in the graph the in-sample deterministic component of the VAR (as discussed in PLR, which is essentially a long-run forecast, starting from 1980 (the same sort of which we are doing right now for the next decades) using the VAR estimates of 12 lags.} The shaded area shows 90% of all the potential paths of the respective VARCTIC. That is, VARCTIC 8 dates the Arctic to be totally ice-free for the first time somewhere between 2052 and 2073 with a probability of 90%. VARCTIC 18 slightly extends that time frame to the year 2079.
For the two models, the median scenario has SIE being less than 1 $times$ $10^6$ km$^2$ by 2054 and 2060 respectively. The 1 $times$ $10^6$ km$^2$ is more likely an interesting quantity since the "regions north of Greenland/Canada will retain some sea ice in the future even though the Arctic can be considered as 'nearly sea ice free' at the end of summer." (wang2009sea). The corresponding credible regions mark the period 2047-2065 for VARCTIC 8 and 2047-2069 for VARCTIC 18 respectively. These dates and time spans range in the close neighborhood of previous climate model simulations (see JahnEtAl2016). For both VARCTICs, less than 5% of the simulated paths hit 0 before 2050, making it an unlikely scenario according to our calculations. In essence, the two models suggest SIE melting at a rate that is slower than DieRu19's results, but much faster than most CMIP5 models StroeveEtAl12, and in line with the latest CMIP6 calculations Notzetal2020.\footnote{Note that augmenting VARCTIC 8 with other greenhouse gases such as methane (CH$_4$) procures near-identical results (i.e., forecasting and forthcoming IRFs). This reinforces the view that $CO_2$ plays a distinct and important role in determining the fate of SIE.}
Nonetheless, it is natural to ask how much we can trust a forecast made 40 years ahead, based on 40 years of data behind. To a large extent, answering this amounts to catalog what types of uncertainty the 90% credible regions incorporate, and those they do not. These reflect both forecasting uncertainty (the presence of shocks) and parameter uncertainty. The latter means the 90% regions reflect what happens to the spread of forecast paths when small disturbances are incorporated in the (estimated) coefficient matrix. In other words, those bands conveniently (and correctly) quantify prediction uncertainty accounting for the fact that we are iterating something that is estimated. All things considered, uncertainty is correctly calibrated as long as the model is correctly specified. As it is the case with any statistical approach, the VARCTIC necessarily assumes that the physical reactions estimated on previous decades' data remain valid for those to come. Thus, if the future holds unprecedented nonlinear mechanisms or previously undetectable relationships\footnote{notz2016observed find that in nearly all CMIP5 models the negative relationship between $CO_2$ and SIE was not prevalent until the second half of the twentieth century.}, the VARCTIC can hardly accommodate for that. In contrast, any intensification of phenomena characterized by our 8 key variables (like Albedo feedback) should be successfully captured out-of-sample. With VARCTIC results being in accord with the recent CMIP6 consensus, our specification seems to capture the main drivers and dynamics of the Arctic ecosystem.\footnote{Though, we acknowledge recent research, which stresses the role of ozone depleting substances (ODS) -- another form of anthropogenic greenhouse gases -- in the warming of the Arctic region over the last decades PolvaniEtAl2020.} Finally, future $CO_2$ emissions are an uncontested source of uncertainty for long-run SIE forecasts. Section (ref) studies how those (and their credible regions) behave under standard forcing scenarios.
Figure (ref) displays impulse response functions -- the response of SIE to a positive shock of one standard deviation to any of the model's M variables. To reflect parameter uncertainty, we additionally report the 90% credible region for each IRF. This means the gray bands contain 90% of the posterior draws from VARCTIC 8. Those are crucial to determine whether the attached IRFs describe significant physical phenomena or not. Particularly, when the credible region extends to both positive and negative sides, the IRF characterizes a phenomenon of negligible importance. In those instances (e.g., many IRFs at horizon $h>24$ months), the posterior mean's (black line) difference from 0 could merely be due to parameter uncertainty, and can be thought of as approximately 0.
The resulting impact of $CO_2$ 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_2$ shocks emerge to have a lasting downward effect on SIE. The relevance of the $CO_2$/SIE relation is not a surprise (notz2016observed). Moreover, this behavior is distinct from other shocks that rather have a significant short-run effect but no significant effect after more than roughly six months. More precisely, the effect of $CO_2$ impulses more than a year to settle in (not significant for approximately 15 months) but ends up having a continuing downward effect on trend SIE of approximately -0.005 $10^6$ km$^2$. This mechanically implies that a one-off $CO_2$ deviation from its predicted value/trend leads to a cumulative impact that is ever increasing in absolute terms (as displayed later in Figure (ref)(b)). It is important to remember that this is the effect of an unexpected increase in $CO_2$ which is to be contrasted with the systematic effect that will be studied later. However, in the framework of this section -- where $CO_2$ is allowed to endogenously respond to Arctic variables -- this is as close as one can get to obtain an experimental/exogenous variation needed to evaluate a dynamic causal effect. -0.005 $10^6$ is roughly 0.1% of the last deterministic trend value of SIE. $CO_2$ shocks, by construction of our linear VAR, have mean 0 and there are approximately as many positive and negative shocks in-sample. The linearity and symmetry of the VAR imply that these durable effects are present for both upward and downward deviations from the deterministic trend.
Other shocks have sizeable impacts that eventually vanish, which is the traditional IRF shape one would expect to see from a VAR on macroeconomic data. For instance, AT and Albedo IRFs clearly have the expected sign. However, they do not have the striking lasting impact of $CO_2$ perturbations. To rationalize the short-lived $IRF(\text{AT}\rightarrow \text{SIE})$, it is worth re-emphasizing what is meant by an AT shock. It is an AT anomaly that is not explicable by (i) the previous state of the system and (ii) other structural shocks ordered before it ($CO_2$, TTC, PR). As an example, one could think of the 2007 record low SIE (at that time) being attributed to an abnormally high “atmospheric flow of warm and humid air” from lower latitudes into the Arctic region graversen2011. As we will see in section (ref), a $CO_2$ shock triggers (with a significant delay) a persistent increase in AT, which eventually impacts SIE downward through the systematic part of the VAR. Thus, the short-lived response of SIE to AT shocks does not rule out a lasting impact of AT on SIE. Rather, it means that when it occurs, the origin of the anomaly is not AT itself, but likely $CO_2$.
Similar to an unforeseeable AT shock, a one-time Albedo shock will not have a lasting effect on SIE neither. This does not preclude Albedo to amplify other shocks as we will see in the next section.\footnote{For a discussion of VARCTIC 18 results, see section (ref).} Finally, a rightful concern is whether IRFs remain unaltered upon sensibly altering the ordering of section (ref). To a large extent, they do. For instance, placing SST and AT before TCC and PR brings no noticeable change. So does moving Albedo from last to second \textcolor{black}{(see section (ref))}.
The melting of SIE is happening much faster than many other phenomena that are also believed to be set in motion by the steady increase of $CO_2$ emissions. Many recent papers (notz2012observations, wang2012sea, serreze2015arctic, notz2017arctic) argue with theory/climate models or correlations that external $CO_2$ forcing is responsible for the long-run trajectory of SIE. Some of these findings led notz2016observed to conclude that climate models severely underestimate the impact of $CO_2$ on SIE.
A rather consensual view is that the very nature of the Arctic system leads to the amplification of such external forcing shocks. An understanding -- from observational data -- on how the Arctic may amplify -- or not -- certain external forces is still pending. Fortunately, a VAR can quantify the contribution of different variables in explaining how a dynamic system responds to an external impulse.
A potential approach that has a long history in econometrics is the use of Granger Causality (GC) tests. Those consist of evaluating predictive causal statements (such as $X_{t-1} \rightarrow Y_{t}$ and $Y_{t-1} \rightarrow X_{t}$) via significance tests in time series regressions. They have been recently advocated for climate applications by mcgrawbarnes2018. Nevertheless, those tests often fall short of answering questions of interest. First, the meaning of the test is not obvious when more than two variables are included and/or if one is interested in multi-horizon impacts. Second, in the wake of a GC test rejection, the block of reduced-form coefficients\footnote{Precisely, we mean coefficients on lags of $X_t$ in a regression of those on $Y_t$ (including lags of $Y_t$ as well).}, which we know to be of some statistical importance, are very hard to interpret. In other words, we know some channel matters, but we have little idea how it matters.\footnote{Similar concerns led us to discard liang2014's quantitative causality since the currently available formula only applies to bivariate systems. Further, it does not allow for contemporaneous relationships which are clearly present in our application (and a feature of almost any discretely sampled multivariate time series).}
In light of the above, we rather opt for IRF Decomposition. As the name suggests, the physical reaction characterized by IRFs will be decomposed as a sum of transmission channels, which contributions' magnitudes and signs are directly informative. Less abstractly, the consequential negative response of SIE to $CO_2$ shocks is likely composed of a direct effect and many entangled indirect effects (e.g., that of AT and Albedo). Understanding those in the dynamic setup of a VAR is much more intricate than in a static regression setting. This is so because IRFs -- for horizons greater than one -- are obtained by iterating predictions, which means $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 transmission mechanism in VARs. It consists, rather simply, of shutting down "channels" and plotting what the response to a shock would be, given that this very channel had been shut (sims2006does, bernanke1997systematic, bachmann2012confidence). We can deploy this methodology to find and quantify the most important channels through which $CO_2$ and temperature shocks impact SIE.
For VARCTIC 8, figures (ref)(a) and (ref)(b) show the responses of SIE to an unexpected increase in one standard-deviation of $CO_2$. The blue line pictures the case of the baseline VARCTIC 8 with 90% credible region. The remaining six lines show the response of SIE to the same shock but shutting 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_2$ shock.\footnote{See bachmann2012confidence for details.}
{The top panel of figure (ref)} reveals -- without great surprise -- the importance of temperature (especially AT) in translating $CO_2$ anomalies into decreasing SIE. That is, we observe that shutting down these channels leads to a smaller absolute response which means that those variables can be considered as amplification channels. Given the atypical shape of the $CO_2$ IRF, the scale of figure (ref)(a) makes less visible the action of channels that only alter the longer-run effect. Since those effects are durably negative (at different levels), their cumulative effect will more clearly reveal their relative importance. Thus, figure (ref)(b) displays the cumulative impact of selected (more important) channels. The two temperature channels are responsible for approximately one fourth of the cumulative effect of $CO_2$ on SIE after 3 years. More precisely, restricting temperature variables to not respond to a positive $CO_2$ shock, decreases (in absolute terms) the after-3-years impact from -0.13 $10^6$ km$^2$ to -0.1 $10^6$ km$^2$. Of course, it was expected that temperature should be a major conductor of such shocks. We also observe similar quantitative effects for both SIT and Albedo in isolation. Most strikingly, we find that the conjunction of the Albedo and SIT amplification channels is responsible for amplifying the effect of $CO_2$ shocks by a non-negligible 50%.
The Albedo-amplification matches evidence reported in several studies (see Pero2012, Bjoerk2013, Park13) using various different methodologies. In contrast, our results for SIT-amplification contribute to a literature where a consensus has yet to emerge. The ice-growth-SIT feedback describes the observation that a thinning of the sea ice cover induces more rapid ice formation during winter, a compensating effect which slows down melting bitz2004mechanism,GoosseEtAl2018. Other studies have emphasized the positive feedback between $SIT$ and $SIE$, where a thinning ice cover further accelerates melting by being less resilient to climate forcing Park13,Kwok2018. Our results unequivocally support the latter to be most empirically prevalent.\footnote{It is plausible that the ice-growth-SIT feedback explains why both $IRF(\text{SIT}\rightarrow \text{SIE})$ (figure (ref)) and SIT's influence on $IRF(CO_2\rightarrow \text{SIE})$ (figure (ref)) take over 6 months to completely settle in --- its seasonal character dampens the (early) positive feedback effect.} Nearly identical results are obtained when using AT, Albedo, PR, and TCC averaged between 60$^{\circ}$N and 90$^{\circ}$N latitude, suggesting most (if not all) of the action comes from higher latitudes -- hence our focus on local processes when explaining those results.
This section focused on how and why SIE responds to $CO_2$ shocks. In section (ref), we rather look at the effect of the systematic increase of $CO_2$ level.
AT-shocks are movements in AT that are not predictable given the past state of the system and are orthogonal to other shocks in the system, most notably $CO_2$. In other words, we are looking at the effect of unexpected higher/lower AT that is uncorrelated with other shocks in the system. As we saw in figure (ref), such AT anomalies have a pronounced short-run effect on trend SIE for no longer than ten months after the shocks. This means that unlike $CO_2$, the cumulative effect of AT disturbances stabilizes about 1.5 years after the event.
{In figure (ref)(c)}, we clearly observe (again) an important role for the thinning of ice and the Albedo effect amplifying the response of SIE to AT shocks. In fact, we see in figure {(ref)(d)} that without them, the long-run impact is the same as the instantaneous one. Thus, this is evidence to suggest that the AT shock's long-run cumulative impact of -0.24 $10^6$ km$^2$ is mostly a result of the action of feedback loops.
If $CO_2$'s trend is mostly or solely affected by factors outside of those considered in the VAR, the forecast of SIE can be improved by treating $CO_2$ forcing as exogenous and using an external forecast rather than the one internally generated by the VAR. Conditional forecasting can be achieved in VARs following the approach of waggoner1999conditional. As we will see, this will markedly sharpen the bands around our forecasts, suggesting that a great amount of uncertainty is related to the future path of $CO_2$ emissions. Additionally, this brings the VARCTIC conceptually closer to standard analyses on the future of the Arctic (StroeveEtAl12, stroeve2018, Notzetal2020).
{In the spirit of SigmondEtAl2018, who constrain the levels of AT in their climate model, we will look at $CO_2$ emissions under three different representative concentration pathways (RCP) and investigate their impact on the evolution of Arctic sea ice.} Figure (ref) shows a steady increase in $CO_2$ emissions over the last three decades, but several RCPs paint different pictures for the trajectory of carbon emissions until the end of the century. Figure (ref) shows the different paths of $CO_2$ under RCP 2.6, RCP 6, and RCP 8.5, as well as the projected path following VARCTIC 8. Most interestingly, we find our completely endogenous and unconditional forecast of $CO_2$ to lay somewhere between the "very bad" RCP 8.5 scenario and the "business-as-usual" RCP 6 one.
Figure (ref) shows VARCTIC 8's projection of Arctic SIE when conditioning the out-of-sample path of $CO_2$ on the three different RCP scenarios. For reference, the figure also includes projected SIE with the future path of $CO_2$ endogenously determined within the model, as discussed above. The pictured effect is dramatic: if emissions were reduced as to follow the RCP 2.6 scenario, whose $CO_2$ emissions are still at the higher boundary of what the Paris Agreement demands, the Arctic would be far from blue and even 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 unconditional VARCTIC 8, but would still be completely gone during the 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 stroeve2018 reported using a very different methodology (extrapolating a linear relationship). Their bivariate (SIE and $CO_2$) analysis suggests the Arctic summer months to be ice-free by 2050. However, in contrast, our results 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 based solely on the reduced form.\footnote{Important to note is the fact that the very last in-sample observations for $CO_2$ even range above the RCP 8.5 values, which generates the slight upward jump in case of the latter scenario.} Overall, these results reinforce the view that anthropogenic $CO_2$ is the main driver behind the current melting of SIE as well as the main source of uncertainty around the future SIE path. Furthermore, the optimistic \emph{RCP 2.6} results suggest that internal variability by itself cannot lead to the complete melting of SIE, even when starting from today's level. Overall, the VARCTIC yields similar conclusions about the importance of $CO_2$ to that of DaiLuoSong19 and notz2016observed. It is reassuring to see that climate models' conclusions can be corroborated by a transparent approach that relies solely and directly on the multivariate time series properties of the observational record.
The previous section documented the evolution of SIE conditional on several $CO_2$ trajectories, treating the latter as an exogenous driver. This section seeks to quantify the importance of feedback effects when it comes to translating an RCP path into SIE loss. That is, we attempt to quantify to which extent the Albedo- and SIT-effects can be held responsible for amplifying the impact of $CO_2$ forcing and thus accelerating the melting of SIE.
Following the findings of section (ref), in which we identified SIT and Albedo to carry potential for mitigating the adverse influence of $CO_2$ on SIE, we ask the question about how SIE would evolve, if SIT and Albedo were to remain constant at a certain level over the forecasting period. In particular, we repeat the forecasting exercise of the previous section for all three RCP scenarios, but keep SIT and Albedo constant until the end of the forecasting 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 artificial shocks to both SIT and Albedo in each forecasting step, which off-set their response to the external forcing variable. As we are modeling a dynamic and interconnected system, these shocks do affect all the other variables (except for $CO_2$ on which we condition our forecast).
Figure (ref) documents the corresponding results for RCP 8.5, RCP 6 and RCP 2.6. For each scenario, we show three different cases: (i) the projection of SIE under the respective RCP; (ii) the evolution of SIE under the respective RCP while keeping Albedo constant at its last in-sample deterministic value; (iii) the projection of SIE while keeping both Albedo and SIT constant at their last respective deterministic value. The latter are shown to be undeniable accelerants. First, fixing Albedo to its 2019 value and thus shutting down this particular long-run amplification effect postpones the date of reaching $1 \times 10^6$ $km^2$ by a bit less than a decade under both RCP 8.5 and RCP 6. Arctic sea ice thickness plays a major role for the reaction and resilience of SIE to anthropogenic forcing. Figure (ref) re-enforces this view by showing that preventing both SIT and Albedo from further decay could potentially postpone the zero-SIE event to the next century under RCP 6. Under \emph{RCP 8.5}, shutting down both amplification channels starting from 2020 leads to SIE crossing the bar of $1 \times 10^6$ $km^2$ about a decade later.\footnote{The graphs are cut at the $1 \times 10^6$ $km^2$ bar as keeping SIT constant (which the thought experiment suggests) is untenable as SIE approaches 0: SIT cannot be constrained to be positive if SIE is 0.} This feeds into the pictured non-linearity and acceleration of SIE loss and provides a potential justification for the finding in DieRu19 that a quadratic trend is a preferable approximation of long-run summer months' SIE evolution.
We proposed the VARCTIC as a middle ground alternative to purely theoretical or statistical modeling. It generates long-run forecasts that embody the interaction of many key variables without the inevitable opacity of climate models. First, we focus our attention on how the Arctic system responds to exogenous impulses and propagates them. Our results show that $CO_2$ anomalies have an unusually lasting effect on SIE which takes about a year to settle in. It is the only impulse that has the property of durably affecting SIE. Albedo and SIT are shown to play an important role in amplifying the response of SIE to $CO_2$ and AT shocks. In both cases, the conjunction of the two effects can double the cumulative impact of such shocks after two years.
Second, we focus on the systematic/deterministic part of the VARCTIC and conduct conditional forecasting experiments that again seek to quantify the effect of anthropogenic $CO_2$ and how feedback loops can amplify it. We condition on the future path of $CO_2$ and show that, within the context of our model, it is the prime source of uncertainty for the long-run forecast of SIE. RCP 8.5 implies 0 September SIE around 2054, RCP 6 says so around 2075 and finally, RCP 2.6 ($\sim$ Paris Accord) implies that such an event would never happen. We conclude the analysis by evaluating to which extent internal knock-on effects can amplify the long-run effect of $CO_2$ forcing. Our results provide statistical backing for the view that $CO_2$ shocks trigger feedbacks of other climate variables (as characterized here by Albedo and SIT), which substantially accelerate the speed at which SIE is headed toward 0.
There are many methodological extensions within the VAR paradigm that could be of interest for future cryosphere research. For instance, Smooth-Transition VARs (with a popular application in AG2012) could be used to accommodate for dynamics evolving over the seasonal cycle. Additionally, ScreenDeser2019 remark the importance of changing weather phenomena that transition through decadal cycles, such as the pacific oscillation. Time-varying parameters VARs could evaluate the quantitative relevance of such phenomena. Finally, some recent attention (chavas2019dynamic) has been given to the potentially non-linear relationship between $CO_2$ and SIE. Methods that blend time series econometrics and Machine Learning of the like in GC_MRF could reveal interesting insights on complex/time-varying relationships in the Arctic.