Robust Estimation and Forecasting of Climate Change Using Score-Driven Ice-Age Models

: We use data on the following climate variables for the period of the last 798 thousand years: global ice volume (Ice t ), atmospheric carbon dioxide level (CO 2, t ), and Antarctic land surface temperature (Temp t ). Those variables are cyclical and are driven by the following strongly exogenous orbital variables: eccentricity of the Earth’s orbit, obliquity, and precession of the equinox. We introduce score-driven ice-age models which use robust ﬁlters of the conditional mean and variance, gener-alizing the updating mechanism and solving the misspeciﬁcation of a recent climate–econometric model (benchmark ice-age model). The score-driven models control for omitted exogenous variables and extreme events, using more general dynamic structures and heteroskedasticity. We ﬁnd that the score-driven models improve the performance of the benchmark ice-age model. We provide out-of-sample forecasts of the climate variables for the last 100 thousand years. We show that during the last 10–15 thousand years of the forecasting period, for which humanity inﬂuenced the Earth’s climate, (i) the forecasts of Ice t are above the observed Ice t , (ii) the forecasts of CO 2, t level are below the observed CO 2, t , and (iii) the forecasts of Temp t are below the observed Temp t . The forecasts for the benchmark ice-age model are reinforced by the score-driven models. forecasting performances of the ice-age model of Castle and Hendry (2020) with those of our score-driven ice-age models. We also report impulse responses among global ice volume, atmospheric CO 2 level, and Antarctic land surface temperature, which are robust for the ice-age model of Castle and Hendry (2020). Likelihood-based model performance metrics, diagnostic tests, and forecasting results indicate that the model of Castle and Hendry (2020) is improved.


Introduction
According to Intergovernmental Panel on Climate Change (2021), compared to the period of 1850-1900, the Earth's global surface temperature for the period of 2081-2100 will very likely rise by 3.3 to 5.7°C under the worst-case scenario, which implies dramatic consequences on the nature and wildlife in terrestrial, wetland, and ocean ecosystems, and on humanity with respect to food and water security, migration, health, higher risk of conflict worldwide, reduction of global economic product, and a possible collapse of the current societal organization. Climate change is partly due to the influence of humanity, which started approximately 10-15 thousand years ago, by commencing agricultural activities such as cultivating plants and livestock (Ruddiman 2005). That influence has significantly increased since the industrial revolution (1769-1840), when human activities such as burning fossil fuels (e.g., coal and oil) have importantly increased, and the influence of humanity has further increased with an accelerating growth rate since then. The Earth's population rose from 1 billion in 1800 to more than 7.8 billion in 2021, which was associated with a significant global-scale economic expansion. One of the consequences is the rising greenhouse gas emissions (e.g., carbon dioxide, CO 2 , nitrous oxide, N 2 O, and methane, CH 4 ), which are directly related to global warming.
The atmospheric CO 2 levels and land surface temperature are related to melting glaciers and sea ice. We name the climate-econometric models of those variables as ice-age models, as per the work of Castle and Hendry (2020). During the 4.5 billion-year history of the Earth, the ice volume, atmospheric CO 2 , and land surface temperature simultaneously Using data for the first 698 thousand years of the sample (for which humanity did not influence the Earth's climate), we present out-of-sample forecasts of the global ice volume, atmospheric CO 2 level, and Antarctic land surface temperature for the period of the last 100 thousand years of the sample. We find that the forecasting results of Castle and Hendry (2020) are robust. For the last 10-15 thousand years of the forecasting period, we find that: (i) the forecasts of global ice volume are above the observed global ice volume, (ii) the forecasts of the atmospheric CO 2 level are below the observed CO 2 level, and (iii) the forecasts of Antarctic land surface temperature are below the observed Antarctic land surface temperature.
The remainder of this paper is organized as follows: Section 2 presents the econometric methods. Section 3 presents the data and the empirical results. Section 4 presents the conclusions.

Benchmark Ice-Age Model
Model specification-In the work of Castle and Hendry (2020, chp. 6), estimation and forecasting results are presented for a general unrestricted model (GUM), named the iceage model. That model specification is the benchmark model of the present paper. The dependent variables y t (3 × 1) of the ice-age model are y t = (Ice t , CO t , Temp t ) , where 'Ice' denotes global ice volume, 'CO 2 ' denotes atmospheric carbon dioxide level, and 'Temp' denotes Antarctic-based land surface temperature. The order of the variables in y t is defined in the work of Castle and Hendry (2020), which we use for all models of the present paper. Correspondingly, the ice-age model is specified as follows: where µ t (3 × 1) is the conditional mean of y t given F t−1 ≡ σ(y 1 , . . . , y t−1 , z 1 , . . . , z t ), the reduced-form error term v t ∼ N 3 (0 3×1 , Σ) has a multivariate i.i.d. normal distribution, where the covariance matrix is Σ ≡ ΩΩ (3 × 3), for which Ω (3 × 3) is a lower triangular matrix with positive elements in the diagonal, and z t (9 × 1) includes strongly exogenous explanatory variables. We initialize µ t using the start values of the dependent variables y 1 . We assume that the maximum modulus of the eigenvalues of Γ 1 , denoted as C µ , is less than one. The elements of the vector of explanatory variables z t are three interacting orbital changes over time affecting solar radiation: z t = (Ec t , Ob t , Pr t , Ec t × Ob t , Ec t × Pr t , Ob t × Pr t , Ec 2 t , Ob 2 t , Pr 2 t ) where 'Ec' measures the eccentricity (i.e., non-circularity) of the Earth's orbit, 'Ob' is obliquity measuring the tilt of the Earth's rotational axis relative to the ecliptic, and 'Pr' is a measure of the precession of the equinox (i.e., circular rotation of the rotational axis itself). The conditional mean µ t in Equation (1) includes the vector of constant parameters γ 0 (3 × 1), and parameter matrices Γ 1 (3 × 3), Γ 2 (3 × 9), and Γ 3 (3 × 9). For the benchmark iceage model, we report estimation and forecasting results for Equation (1) which excludes the outlier-dummies. The same outlier-dummies are also excluded from the score-driven iceage models of the present paper, because those models are robust to extreme observations (Harvey 2013).
Impulse responses-We estimate the dynamic effects of the i.i.d. structural-form error term t ≡ Ω −1 v t ∼ N 3 (0 3×1 , I 3 ). The corresponding IRFs are defined as follows: The IRFs are identified using the sign restrictions for the contemporaneous effects among the elements of v t , which are based on simulations of matrix Ω, according to the following procedure (Rubio-Ramirez et al. 2010): (i) A K × K matrixK of independent N(0, 1) numbers is simulated. (ii) The QR decomposition ofK is performed, and the resulting matrices are denoted asQ (orthogonal matrix) andR (upper triangular matrix). (iii) We defineΩ ≡ Ω ×Q , for which the estimates of Ω are used. For the IRFs, 3 million simulations ofK are generated, and only those simulations are used that satisfy the sign restrictions of Table 1, and for each simulation Ω is replaced byΩ in the IRF formulas. For the simulations that satisfy the sign restrictions, we report the mean ± one standard deviation estimates of the IRFs. The sign restrictions of Table 1 are motivated as follows: (i) In the work of Castle and Hendry (2020) the correlation coefficient estimates among the residuals of the ice-age model are reported, which motivate the sign restrictions of Table 1. (ii) For the negative effects of CO 2 on Ice t , and for the positive effects of Ice t on Temp t of Table 1, we refer to the work of Qin and Buehler (2012). For the positive interaction effects between CO 2 and Temp t of Table 1, we also refer to the works of Jouzel et al. (2007) and Lüthi et al. (2008). (iii) For the negative effects of Ice t on CO 2 of Table 1, we refer to the work of Wadham et al. (2019). (iv) For the positive effects of Temp t on CO 2 of Table 1, we refer to the work of Archer et al. (2004). (v) For the negative interaction effects between Temp t and Ice t of Table 1, we refer to the work of Bronselaer et al. (2018).

Score-Driven Ice-Age Models
Castle and Hendry (2020, p. 102) recognize that the benchmark ice-age model suffers from residual autocorrelation. This may be due to omitted exogenous variables, such as the Sun's radiation output or changes in the Earth's magnetic poles. This may also be due to the presence of heteroskedasticity generated by the omitted variables, which might cause fat-tails in the error distribution. There are alternative ways to address those misspecification issues from an econometric perspective. In this paper, we address those issues using the flexible and robust class of score-driven models.

Score-Driven Homoskedastic Ice-Age Model
Model specification-The score-driven ice-age model of this paper uses the score-driven model specification in the work of Harvey (2013, p. 56). The model is specified as follows: where µ t (3 × 1) is the conditional mean of y t given F t−1 ≡ σ(y 1 , . . . , y t−1 , z 1 , . . . , z t , µ 1 ), v t is the multivariate i.i.d. reduced-form error term, z t (9 × 1) is the vector of strongly exogenous explanatory variables, and u t (3 × 1) is the vector of scaled score functions (Harvey 2013). The assumption of strict exogeneity of z t is from the work of Harvey (2013, p. 56), which is supported for z t of Equation (2) in the work of Castle and Hendry (2020, p. 95). We initialize µ t using the start values of the dependent variables y 1 (Harvey 2013). The conditional mean µ t includes the vector of constant parameters γ 0 (3 × 1), and the parameter matrices Γ 1 (3 × 3), Γ 2 (3 × 9), Γ 3 (3 × 9), and Ψ (3 × 3). We assume that the maximum modulus of the eigenvalues of Γ 1 , denoted as C µ , is less than one. For Γ 1 , Γ 2 , and Γ 3 , we use the restrictions of Castle and Hendry (2020), a decision which is motivated by the following general-to-specific model selection procedure. First, for the parameter estimation of the score-driven ice-age model, we start with the aforementioned restrictions of Γ 2 and Γ 3 , and we use unrestricted parameter matrices for Γ 1 and Ψ. Second, we restrict those parameters of Γ 1 and Ψ to zero which are not significant at the 1% level. The same elements of Γ 1 and Ψ are restricted for the score-driven models as for matrix Γ 1 in the work of Castle and Hendry (2020). Hence, the following elements of Ψ are not restricted to zero: Ψ 1,1 , Ψ 1,3 , Ψ 2,2 , Ψ 2,3 , Ψ 3,2 , and Ψ 3,3 .
The reduced-form error term v t ∼ t 3 (0, Σ, ν) has a multivariate i.i.d. t-distribution, where the scale matrix is Σ ≡ ΩΩ (3 × 3), for which Ω (3 × 3) is a lower triangular squared matrix with positive elements in the diagonal, and ν > 2 is the degrees of freedom parameter (the restriction on the parameter space ν > 2 ensures that the covariance matrix of v t is well-defined).
The scaled score function u t is defined as follows. The log of the conditional density of y t is: is the vector of time-invariant parameters, which includes the elements of γ 0 , Γ 1 , Γ 2 , Γ 3 , Ψ, Ω, and ν. The partial derivative of the log conditional density ln f (y t |F t−1 ; Θ) with respect to µ t is (Harvey 2013): The scaled score function u t is defined in the second equality of Equation (7), where v t is multiplied by [1 + (v t Σ −1 v t )/ν] −1 = ν/(ν + v t Σ −1 v t ) ∈ (0, 1). Therefore, the scaled score function is bounded by the reduced-form error term: |u t | < |v t |. All elements of u t are bounded functions of v t for ν < ∞ (Harvey 2013), hence all moments of u t are well-defined. In the work of Harvey (2013), it is shown that u t is multivariate i.i.d. with mean zero and a covariance matrix: We also note that if ν → ∞, then u t → p v t . In the limiting case, Equations (3) and (4) provide a VARMAX(1,1) structure for the dependent variables: The benchmark ice-age model is a special case of the score-driven ice-age model, because if ν → ∞ and Ψ = Γ 1 for Equations (4) and (5), then we obtain Equation (1). This can also be seen for the limiting case for ν → ∞ in Equation (9), using Ψ = Γ 1 . We name Equation (9) the score-driven homoskedastic Gaussian ice-age model.
All score-driven models are estimated using the ML method. For the ML estimation, we assume correct model specifications for all econometric models of the present paper. Under that assumption, the standard errors of the ML estimates are consistently estimated using the outer product of the gradient of the LL function. Another consequence of the correct model specification assumption is that the updating terms of the score-driven filters of this paper, asymptotically and at the true values of parameters, are white noise vectors. For technical details on the statistical inference of score-driven models, we refer to the works of Harvey and Chakravarty (2008), Harvey (2013), Creal et al. (2008Creal et al. ( , 2011Creal et al. ( , 2013, and Blasques et al. (2021). For the multivariate score-driven models (such as the score-driven ice-age models), we also refer to the recent works of Blazsek et al. (2020Blazsek et al. ( , 2021aBlazsek et al. ( , 2021b. Impulse responses-First, we define the vector of the structural-form error terms t (3 × 1). The variance of the reduced-form error term v t ∼ t 3 (0, Σ, ν) is factorized, as follows: Based on that, the following multivariate i.i.d. structural-form error term t is introduced: Furthermore, by substituting Equation (11) into Equation (7), u t as a function of the structural-form error term is: Second, from Equations (4) and (5), the nonlinear MA(∞) representation of y t is: We focus on the impulse responses for the dependent variables y t , because the variables within z t are strongly exogenous to humanity. From y t , we are particularly interested in the dynamic effects of the atmospheric carbon dioxide level on the Antarctic-based land surface temperature, because humanity has influence on the CO 2 emissions. Thus, the new measurement method of impulse responses for the score-driven ice-age model of this paper may have policy implications in relation to carbon dioxide emissions regulation. The contemporaneous and dynamic effects of the structural-form error term t , respectively, are given by the following equations: We compare the IRFs of the score-driven ice-age models and the IRFs of the benchmark ice-age model, which are not reported in the work of Castle and Hendry (2020).
In Equation (15), the dynamic interaction effects are time-dependent due toD t . In the empirical application of the present paper, we replaceD t by the sample average (White 2001), where the sample average is estimated for the last 10 observations of the sample, motivated by the work of Castle and Hendry (2020, p. 111). There are alternative ways for the estimation of dynamic interaction effects for nonlinear models (Lütkepohl 2005), hence the IRF estimation of our paper may be modified in future applications.
Finally, we also report contemporaneous and dynamic effects of the structural-form error term t for the score-driven ice-age model for the multivariate normal distribution: The IRFs of this section are identified using the procedure of Rubio-Ramirez et al. (2010).

Score-Driven Heteroskedastic Ice-Age Model
Model specification-We extend the score-driven ice-age model for the homoskedastic multivariate t-distribution, by considering score-driven conditional heteroskedasticity with constant correlation coefficients for the reduced-form error term v t . The model is specified as follows: where µ t (3 × 1) is the conditional mean of y t given F t−1 , which is defined later in this section, v t is the multivariate i.i.d. reduced-form error term, z t (9 × 1) is the vector of strongly exogenous explanatory variables, and u t (3 × 1) is the vector of scaled score functions. We initialize µ t using y 1 (Harvey 2013). The conditional mean µ t includes the vector of constant parameters γ 0 (3 × 1), and the parameter matrices Γ 1 (3 × 3), Γ 2 (3 × 9), Γ 3 (3 × 9), and Ψ (3 × 3). For the parameters of µ t , we use the same restrictions as for the homoskedastic score-driven ice-age model for the t-distribution. We assume that the maximum modulus of the eigenvalues of Γ 1 , denoted as C µ , is less than one.
is a time-varying diagonal matrix with the score-driven scales of each time series, and R (3 × 3) is the time-invariant correlation matrix. This specification assumes that the correlation coefficients are constant over time, which can be extended according to the results of Creal et al. (2011) to dynamic correlation coefficients, using a score-driven volatility plus correlation model for the t-distribution.
The positive definiteness of R and boundedness of the correlation coefficients for (−1, 1) are ensured using the following specification: is a diagonal matrix in which the elements of the diagonal are the square roots of the elements of the diagonal of the positive definite matrix Q (3 × 3) (Engle 2002). The positive definiteness of Q is ensured using the Cholesky decomposition Q = ΩΩ , in which Ω (3 × 3) is a lower triangular matrix with positive elements in the diagonal. For parameter identification reasons, each element of the diagonal of Ω is restricted to one. Furthermore, D t is specified as follows: We specify the filters λ i,t in Equation (21) as follows: where sgn(·) is the signum function, and α * i for i = 1, 2, 3 measure asymmetric effects in the conditional scale of the dependent variables. We initialize λ i,t for i = 1, 2, 3 using the unconditional mean E(λ i,t ) = ω i /(1 − β i ) for i = 1, 2, 3, respectively (Harvey 2013). In the following, we define the updating terms u t and e i,t . For the covariance stationarity of λ i,t , asymptotically at the true values of parameters, it is required that C i,λ = |β i | < 1 for i = 1, 2, 3.
First, the scaled score function u t is defined as follows. The log conditional density of y t is: The scaled score function u t is defined in the second equality of Equation (24) . Therefore, the scaled score function is bounded by the reduced-form error term: |u t | < |v t |. All elements of u t are bounded functions of v t for ν < ∞, hence all moments of u t are well-defined. The scaled score function u t |(F t−1 ; Θ) has a zero conditional mean and the following conditional covariance matrix: The latter result is an extension of the work of Harvey (2013, p. 206). Second, the updating term e i,t for i = 1, 2, 3 is defined as follows. The conditional distributions of the marginals of y t are y i,t |(F t−1 ; Θ) ∼ t[µ i,t , exp(λ i,t ), ν] for i = 1, 2, 3 (Kibria and Joarder 2006). The log of the conditional density of y i,t |(F t−1 ; Θ) for i = 1, 2, 3 is We define the updating term of the filter λ i,t for i = 1, 2, 3 as follows: Equations (22) and (27) are the Beta-t-EGARCH with leverage effects model of Harvey and Chakravarty (2008) (see also Creal et al. 2013 andHarvey 2013).
Impulse responses-First, we define the structural-form error terms t (3 × 1). The conditional variance of the reduced-form error term v t |(F t−1 ; Θ) ∼ t 3 (0, Σ t , ν) is factorized, as follows: Furthermore, by substituting Equation (29) into Equation (24), u t as a function of the structural-form error term is: Second, from Equations (19) and (20), the nonlinear MA(∞) representation of y t is: The contemporaneous and dynamic effects of the structural-form error term t , respectively, are: whereD t is given by Equation (16). In Equations (32) and (33), the dynamic interaction effects are time-dependent due to D t andD t . We replace D t andD t by their sample averages for the last 10 observations of the sample. The IRFs are identified using the procedure of Rubio-Ramirez et al. (2010). We also refer to the work of Lütkepohl (2005) for alternative estimation methods.

Data
The dependent variables of this study are global ice volume Ice t , atmospheric CO 2,t , and Antarctic land surface temperature Temp t . The data source of global ice volume Ice t is in the work of Lisiecki and Raymo (2005), in which time series of the δ 18 O, obtained from calcium carbonate (CaCO 3 ) shells of foraminifera, are used to approximate temperature. Those authors use benthic records of foraminifera from seafloor sediment, which are collected at 57 globally distributed sites. Those sites are well-distributed in latitude, longitude, and depth in the Atlantic, Pacific, and Indian Oceans. The data source of atmospheric CO 2,t is in the work of Lüthi et al. (2008), in which changes in past atmospheric CO 2 concentrations are determined by measuring the composition of air trapped in ice cores from Antarctica. Within the European Project for Ice Coring in Antarctica (EPICA), two deep ice cores have been drilled at the Kohnen Station and the Concordia Station (Dome C). The drillings were stopped at, or a few meters above, bedrock at a depth of 2774 m and 3270 m, respectively. The data source of Antarctic land surface temperature Temp t is in the work of Jouzel et al. (2007), in which temperature data were obtained within the EPICA at the Concordia Station (Dome C), using deuterium δD ice measurements from the surface down to 3259.7 m.
The exogenous explanatory variables of this study are eccentricity of the Earth's orbit, obliquity of the Earth's rotational axis relative to the ecliptic, and precession of the equinox. The source of those data is in the work of Paillard et al. (1996), in which the AnalySeries software is used, to provide measurements of eccentricity, obliquity, and precession.
In Table 2, the descriptive statistics of the dependent and the explanatory variables are presented. The table shows the definition of each variable, observation period, units of measurement, data sources, and some additional descriptive statistics. In Figures 1 and 2, the evolution of the dependent and explanatory variables, respectively, are presented. According to Figure 1b,c, atmospheric CO 2,t and Antarctic land surface temperature Temp t , respectively, remarkably are in unison. In Figure 1, it can also be noticed that global ice volume Ice t moves in the opposite direction from CO 2,t and Temp t , creating the ice-age and inter-glacial periods periodically. The seasonality of the dependent variables (Figure 1), which is due to the three main interacting orbital changes over time affecting solar radiation ( Figure 2), is clearly observed in these figures.
Additional explanatory variables, which are exogenous to humanity are omitted from the econometric models of this paper. For example, the following variables are omitted: (i) the variations in the Sun's radiation output, (ii) volcanic eruption particles in the atmosphere and ice cover, and (iii) changes in the magnetic poles. In Appendix A, we present details of those exogenous variables, and we also present why those variables are not included in the climate-econometric models.

Estimation Results
In Table 3, the ML parameter estimates for the (i) benchmark ice-age model, (ii) scoredriven homoskedastic ice-age model for the normal distribution, (iii) score-driven homoskedastic ice-age model for the t-distribution, and (iv) score-driven heteroskedastic ice-age model for the t-distribution are reported. According to the table, the parameters for Γ 1 , Γ 2 , Γ 3 , and Ψ for all models are significantly different from zero. For model (iv), significant and asymmetric volatility dynamics are estimated, as α i and β i for i = 1, 2, 3, and α * i for i = 1, 2, are significantly different from zero. In Table 4, the statistical performance metrics and some diagnostic test results for the models of Table 3 are reported. The statistical performances are compared using the LL, Akaike information criterion (AIC), Bayesian information criterion (BIC), and Hannan-Quinn criterion (HQC) metrics (Harvey 2013, p. 56). The statistical performance of the score-driven heteroskedastic ice-age model for the t-distribution is superior to the statistical performances of other specifications of Tables 3 and 4. For all models, the C µ and C i,λ for i = 1, 2, 3 statistics support the covariance stationarity of µ t and λ i,t , respectively. As a diagnostic test of the residuals, we use the Ljung-Box test (Ljung and Box 1978) for all elements of v t , t , u t , and e t . The diagnostic tests for the score functions u t , and e t are motivated by the work of Harvey (2013, p. 55), due to the robustness to extreme observations of the score-function-based Ljung-Box test. The Ljung-Box test results indicate that the ice-age model of Castle and Hendry (2020) is not fully supported, which is noted in the same work (p. 102, footnote 4). From the score-driven ice-age specifications, full support is provided for the most general score-driven heteroskedastic ice-age model for the t-distribution.
In the following, we present the parameter estimates for the benchmark ice-age model (Equation (1)) and the score-driven heteroskedastic ice-age model for the t-distribution (Equation (20)). First, the estimates of Equation (1) are given by ( The estimates correspond to the estimates of Castle and Hendry (2020, p. 104). Second, the estimates of Equation (20) To compare the parameter estimates of Equations (34)-(36) with those of Equations (37)-(39), the dynamic interaction effects for global ice volume, atmospheric CO 2 , and Antarctic land surface temperature are studied using the IRFs.
In Figures 3-6, the IRFs for the (i) benchmark ice-age model, (ii) score-driven homoskedastic ice-age model for the normal distribution, (iii) score-driven homoskedastic ice-age model for the t-distribution, and (iv) score-driven heteroskedastic ice-age model for the t-distribution, respectively, are reported. The IRF figures indicate that the signs of the dynamic interaction effects are coherent with the signs of the same effects of the aforementioned works of Archer et al. By comparing the IRF estimates of the benchmark ice-age model with those of the score-driven ice-age models, for several panels of Figures 3-6, stronger effects are measured for the score-driven ice-age models than for the benchmark ice-age model. The strongest effects are measured for the score-driven heteroskedastic ice-age model for the t-distribution ( Figure 6). We find the following differences: (i) For the benchmark ice-age model, the dynamic effects of a unit Ice t shock (i.e., measured based on the δ 18 O proxy) on CO 2,t are less than −0.25 (i.e., −195 gigatonnes of CO 2 ), while for the score-driven models the same effect is stronger and it is approximately −0.35 (i.e., −273 gigatonnes of CO 2 ) (Figures 3-6, Panel (d)). (ii) For the benchmark ice-age model, the dynamic effects of a unit Temp t shock on CO 2,t are less than 0.25 (i.e., 195 gigatonnes of CO 2 ), while for the score-driven models the same effect is stronger, and it is between 0.30 and 0.35 (i.e., 234 and 273 gigatonnes of CO 2 , respectively) (Figures 3-6, Panel (f)). (iii) For the benchmark ice-age model, the dynamic effects of a unit Ice t shock (i.e., measured based on the δ 18 O proxy) on Temp t are less than −0.35°C, while for the score-driven models the same effect is stronger, reaching an estimate between −0.40 and −0.45°C, respectively (Figures 3-6, Panel (g)). (iv) For the benchmark ice-age model, the dynamic effects of a unit CO 2,t shock (i.e., an increase of 780 gigatonnes of CO 2 in the atmosphere) on Temp t is approximately 0.35°C, while it is above 0.40°C for the score-driven ice-age models (Figures 3-6, Panel (h)).
In Figure 7, we present the scaled score function u t as a function of the structural-form error term t . The figure presents the estimates for the score-driven heteroskedastic ice-age model for the t-distribution. In the three-dimensional graphs of Figure 7, we present the elements of u t from Equation (30) as functions of 1,t and 2,t , where 3,t = 0 for the purpose of illustration. For the D t term of Equation (30), we use the unconditional mean estimate of λ i,t for i = 1, 2, 3, which isÊ(λ i,t ) =ω i /(1 −β i ) for i = 1, 2, 3, respectively.
(a). u 1,t as a function of 1,t and 2,t (b). u 2,t as a function of 1,t and 2,t (c). u 3,t as a function of 1,t and 2,t Figure 7. Robustness of the scaled score function to extreme values. Note: 3,t = 0 is assumed for this figure. Notes: Not available (NA); log-likelihood (LL); Akaike information criterion (AIC); Bayesian information criterion (BIC); Hannan-Quinn criterion (HQC); Ljung-Box (LB). C µ is the maximum modulus of eigenvalues of the matrix Γ 1 , for which C 1 < 1 indicates covariance stationarity of the location filter. C λ,i = |β i | for i = 1, 2, 3, for which C λ,i < 1 for i = 1, 2, 3 indicates covariance stationarity of the scale filter. The lag-order for the LB test is 28 √ T. *** and ** indicate significance at the 1% and 5% levels, respectively.

Forecasting Results
In Table 5, the multi-step-ahead out-of-sample forecasting performances of the benchmark ice-age model and all score-driven ice-age models of the present paper are compared. The following climate variables are predicted: Ice t , CO 2,t , and Temp t . The estimation window is for the period of 798 thousand years ago to 101 thousand years ago (T = 698), for which humanity did not influence the Earth's climate. The multi-step ahead forecasting window is for the last 100 thousand years (T f = 100). We use two loss functions for forecasting performance evaluation: (i) mean square error (MSE), and (ii) mean absolute error (MAE). These loss functions are averaged for different periods of the last 100 thousand years (Table 5). For most of the cases, the MSE and MAE results indicate that for the periods of the last 100 thousand years to the last 30-40 thousand years, the benchmark ice-age model provides the most precise forecasts (Table 5). The results indicate for all variables that for the most recent period of the last 20-30 thousand years the score-driven homoskedastic ice-age model for the t-distribution provides the most precise forecasts (Table 5).
In Figure 8, the multi-step ahead out-of-sample forecasts of Ice t , CO 2,t , and Temp t for all ice-age models of this paper are presented. The figure includes the observed values of Ice t , CO 2,t , and Temp t , the forecasts of these variables, and the forecasts ± one standard deviation estimates of the forecasts. The figure shows the following results for the most recent period of the sample, when humanity impacted the Earth's climate. For the last 10-15 thousand years of the forecasting window, the observed values of global ice volume are below the forecast interval, indicating unexpectedly low levels of global ice volume. For the same period, the observed levels of CO 2 and Antarctic land surface temperature are above the forecast interval, indicating unexpectedly high levels of CO 2,t and Antarctic land surface temperature. These multi-step-ahead forecasting results are robust for the different econometric models (Figure 8), and are consistent with the results in the work of Castle and Hendry (2020).
The reason the observed values of the three climate variables leave the forecasting interval in Figure 8 may be due to model misspecification or the growing influence of humanity on the Earth's climate. To study this, we repeat the multi-step ahead out-ofsample forecasting exercise for the forecasting period of 223 to 124 thousand years ago (i.e., a 100-thousand-year period), using the estimation window of the period of 798 to 224 years ago (Figure 9). For these estimation and forecasting periods, humanity did not influence the Earth's climate. The last part of the forecasting window matches the local maximum and minimum points of the climate variables in the present time (Figure 1). We find that the observed values of the climate variables leave the forecasting interval much less in Figure 9 than in Figure 8. This result is the clearest for the most general score-driven heteroskedastic ice-age model for the t-distribution (i.e., Panels (j), (k), and (l) of Figures 8 and 9).
In Figure 10, the one-step ahead out-of-sample forecasts of Ice t , CO 2,t , and Temp t for all ice-age models of this paper are presented. We use a rolling-window approach for estimation and forecasting. The figure includes the observed values of Ice t , CO 2,t , and Temp t , the forecasts of these variables, and the forecasts ± one standard deviation estimates of the forecasts. For the last 10-15 thousand years, the figure shows a significant decrease in the level of global ice volume and significant increases in CO 2 and Antarctic land surface temperature. For the same period, the observed values of global ice volume are unexpectedly located below the mean forecasts, and the observed values of CO 2 and Antarctic land surface temperature are unexpectedly located above the mean forecasts. These forecasting results are robust for the different models. (a). Benchmark ice-age Ice t (b). Benchmark ice-age CO 2,t (c). Benchmark ice-age Temp t (d). Score-driven homoskedastic Gaussian ice-age Ice t (e). Score-driven homoskedastic Gaussian ice-age CO 2,t (f). Score-driven homoskedastic Gaussian ice-age Temp t (g). Score-driven homoskedastic t ice-age Ice t (h). Score-driven homoskedastic t ice-age CO 2,t (i). Score-driven homoskedastic t ice-age Temp t (j). Score-driven heteroskedastic t ice-age Ice t (k). Score-driven heteroskedastic t ice-age CO 2,t (l). Score-driven heteroskedastic t ice-age Temp t Figure 8. Multi-step ahead out-of-sample forecasts for the period of the last 100 thousand years. Notes: The confidence interval is mean ± one standard deviation.
The true values are indicated by thick lines.
(a). Benchmark ice-age Ice t (b). Benchmark ice-age CO 2,t (c). Benchmark ice-age Temp t (d). Score-driven homoskedastic Gaussian ice-age Ice t (e). Score-driven homoskedastic Gaussian ice-age CO 2,t (f). Score-driven homoskedastic Gaussian ice-age Temp t (g). Score-driven homoskedastic t ice-age Ice t (h). Score-driven homoskedastic t ice-age CO 2,t (i). Score-driven homoskedastic t ice-age Temp t (j). Score-driven heteroskedastic t ice-age Ice t (k). Score-driven heteroskedastic t ice-age CO 2,t (l). Score-driven heteroskedastic t ice-age Temp t Figure 9. Multi-step ahead out-of-sample forecasts for the period of 223 to 124 thousand years ago. Notes: The confidence interval is mean ± one standard deviation.
The true values are indicated by thick lines.
(a). Benchmark ice-age Ice t (b). Benchmark ice-age CO 2,t (c). Benchmark ice-age Temp t (d). Score-driven homoskedastic Gaussian ice-age Ice t (e). Score-driven homoskedastic Gaussian ice-age CO 2,t (f). Score-driven homoskedastic Gaussian ice-age Temp t (g). Score-driven homoskedastic t ice-age Ice t (h). Score-driven homoskedastic t ice-age CO 2,t (i). Score-driven homoskedastic t ice-age Temp t (j). Score-driven heteroskedastic t ice-age Ice t (k). Score-driven heteroskedastic t ice-age CO 2,t (l). Score-driven heteroskedastic t ice-age Temp t Figure 10. One-step ahead out-of-sample forecasts for the last 100 thousand years. Notes: The confidence interval is mean ± one standard deviation. The true values are indicated by thick lines.

Conclusions
We have used data for climate and orbital variables for the period of the last 798 thousand years, to solve the dynamic misspecification of the benchmark ice-age model. We have improved the model specification using robust score-driven models with heteroskedastic errors from the Student's t-distribution. We have compared the statistical and forecasting performance of the benchmark ice-age model and those of our score-driven ice-age models. The statistical performance metrics and diagnostic tests have indicated that the dynamic specification of the benchmark ice-age model is improved and that the evidence of dynamic misspecification is no longer there. We have reported impulse responses for global ice volume, atmospheric CO 2 , and Antarctic land surface temperature.
The forecasting results of the benchmark ice-age model are robust to the dynamic model misspecification, using data for the first 698 thousand years of the sample to forecast global ice volume, atmospheric CO 2 , and Antarctic land surface temperature for the last 100 thousand years of the sample. For the last 10 to 15 thousand years when humanity influenced the Earth's climate, we have found the following results: (i) the forecasts of global ice volume are above the observed global ice volume, (ii) the forecasts of the atmospheric CO 2 level are below the observed CO 2 level, and (iii) the forecasts of Antarctic land surface temperature are below the observed Antarctic land surface temperature. These results may indicate the increasing influence of humanity on the Earth's climate, and provide a motivation to take further proactive actions to significantly reduce the greenhouse gas emissions, while respond to the most important challenge of the 21st century: global warming.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In all ice-age models we include exogenous orbital variables which influence the Earth's climate in the long run, but we omit further exogenous variables that may also influence the climate. The omitted variables are: (i) the variations in the Sun's radiation output, (ii) volcanic eruption particles in the atmosphere and ice cover, and (iii) changes in the magnetic poles. In this Appendix, we present why those variables are not included in the climate-econometric models of the present work: (i) For the variable variations in the Sun's radiation output, according to the literature (e.g., Feulner and Rahmstorf 2010;Jones et al. 2012;Anet et al. 2013;Meehl et al. 2013;Ineson et al. 2015;Maycock et al. 2015), the global warming of the last few decades is too large to be caused by solar activity. There are short-and long-term variations in the Sun's radiation output. The short-term variations are the 11-year cycles (i.e., Schwabe cycles), which are part of the 22-year magnetic cycles (i.e., Hale cycles). For this paper, the long-term variations in the Sun's radiation output, e.g., grand maximums and grand minimums which may last several centuries, and their effects on the Earth's climate are more interesting. An example of the long-term variations is the grand minimum for the period of 1645 to 1715 approximately (i.e., the Maunder minimum), during which the solar magnetism diminished, sunspots were less frequent, and ultraviolet radiation was reduced.
The currently low level of solar activity motivated several works in the literature, to study the effects of a possible grand minimum on global surface temperatures. Feulner and Rahmstorf (2010) find that the temperature offset due to the lower level of solar activity is no more than −0.3°C for the 21st century. The same authors note that the effects of the 11-year solar cycles and the grand minimum are negligible regarding the global warming. Similar results are reported in the work of Anet et al. (2013), indicating −0.2 to −0.3°C offsetting effects for the period of 2081-2100. The work of Jones et al. (2012) uses 9000-year data on the Sun's radiation output and those authors find that the temperature offset due to the lower level of solar activity is no more than −0.06 to −0.1°C for the 21st century. In the work of Meehl et al. (2013), solar irradiance is reduced by 0.25% for the period of 2020 to 2070. The results indicate that after the initial reduction of solar radiation in 2020, global surface air temperature cools by up to several tenths of a degree Centigrade compared to the reference simulation (i.e., it is not significant), and by the end of the grand solar minimum in 2070, the warming nearly catches up to the reference simulation. Meehl et al. (2013) conclude that a future grand solar minimum may slow down but not stop global warming. In the work of Ineson et al. (2015), the effects of a future scenario are investigated, for which the solar activity decreases to Maunder minimum-like conditions by 2050. The results show that the impact of that scenario on winter northern European surface temperatures over the late 21st century would not be significant. Similarly, in the work of Maycock et al. (2015), Maunder minimum-like conditions are considered for the 21st century, and the results show that the impact of such reduced solar activity is a 1.2°C cooling at the stratopause, which is more significant than the results of other works in the literature, but it is much lower than projected temperature increases due to global warming.
With respect to the observation period of solar activity, the short-term solar cycles have been observed since 1610 (i.e., the first telescopic observations). Moreover, we also refer to the proxy records of solar activity in the work of Steinhilber et al. (2008), which creates a time series of solar modulation potential for the last 9300 years, that, to the best of our knowledge, is the longest available observation period of solar activity in the literature.
We omit the variable on solar activity from the econometric models of this paper, because the global warming of the last few decades is too large to be caused by solar activity, and we do not have data on solar activity for the full sample period of the last 798 thousand years.
(ii) The impact of volcanic eruptions on the climate depends on the location of the volcano. Eruptions in the tropics influence the climate more than eruptions at mid or high latitudes which only influence the hemisphere they are within. The gases and dust particles entering the Earth's atmosphere from volcanic eruptions influence the climate. The dust particles cool the planet by shading incoming solar radiation, which can last from months to years. Particularly effective cooling particles emitted during volcanic eruptions are the sulfuric gases, which move into the stratosphere and combining with water form sulfate aerosols that reflect the incoming solar radiation. The sulfate aerosols may stay in the stratosphere for 3-4 years. The sulfate aerosol absorption heats the stratosphere, but a reduction of downward radiation at the tropopause cools the troposphere and the underlying surface (Stenchikov et al. 1998;Kirchner et al. 1999). Volcanic eruptions also influence global warming when greenhouse gases, such as water vapor and CO 2 are released into the atmosphere. The cooling effects of dust and sulfate aerosol particles, and the relatively small volume of heating greenhouse gases emitted by volcanic eruptions (compared to the volume of greenhouse gases in the Earth's atmosphere) do not significantly influence the Earth's climate in the long term. Moreover, data on volcanic eruptions for the period of the last 798 thousand years of this paper are not available. Hence, the omission of the variable volcanic eruption particles in the atmosphere and ice cover in the econometric models.
(iii) According to NASA (2021), changes in the magnetic poles do not influence the Earth's climate, because there is little scientific evidence of any significant links between Earth's drifting magnetic poles and climate. The changes in the magnetic poles can be classified into three categories: First, shifts in magnetic pole locations (with a speed of approximately 16 to 55 kilometers per year). Second, Earth's magnetic north and south poles swap locations (the frequency is variable, but on average it takes place in every 300 thousand years, with the last one taking place about 780 thousand years ago). Third, geomagnetic excursions which are shorter-lived but significant changes in the magnetic field's intensity with duration of a few centuries to a few tens of thousands of years (the last major excursion, named the Laschamps event, took place about 41,500 years ago, when the magnetic field weakened significantly, the poles reversed, and flipped back about 500 years later). The changes in the magnetic poles does not influence the Earth's climate, at least because of two reasons NASA (2021): First, there is insufficient energy in the Earth's upper atmosphere (where electromagnetic currents exist), to influence the Earth's climate. Second, there is no known physical mechanism capable of connecting weather conditions at the Earth's surface with electromagnetic currents in space. Hence, the omission of changes in the magnetic poles.