Next Article in Journal
Identification in Parametric Models: The Minimum Hellinger Distance Criterion
Next Article in Special Issue
Green Bonds for the Transition to a Low-Carbon Economy
Previous Article in Journal
The Impact of COVID-19 on Airfares—A Machine Learning Counterfactual Analysis
Previous Article in Special Issue
Climate Finance: Mapping Air Pollution and Finance Market in Time Series
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Szabolcs Blazsek
1 and
Alvaro Escribano
2,*
1
School of Business, Universidad Francisco Marroquín, Guatemala 01010, Guatemala
2
Department of Economics, Universidad Carlos III de Madrid, Calle Madrid 126, 28903 Getafe, Spain
*
Author to whom correspondence should be addressed.
Submission received: 19 October 2021 / Revised: 21 January 2022 / Accepted: 11 February 2022 / Published: 16 February 2022
(This article belongs to the Collection Econometric Analysis of Climate Change)

Abstract

:
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 filters of the conditional mean and variance, generalizing the updating mechanism and solving the misspecification 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 find 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 influenced 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.

1. 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 changed, driven by exogenous orbital variables, such as (i) changes in the non-circularity of the Earth’s orbit with a period of 100 thousand years, (ii) changes in the tilt of the Earth’s rotational axis relative to the ecliptic (i.e., the plane of the Earth’s orbit around the Sun) with a period of 41 thousand years, and (iii) circular rotation of the rotational axis itself, which changes the season at which the Earth’s orbit is nearest to the Sun, with a period that is between 19 to 23 thousand years (it is variable due to the changes in the tilt of the rotational axis).
As in the work of Castle and Hendry (2020), we use data for climate and orbital variables for the period of the last 798 thousand years, where the climate variables in the model are global ice volume, atmospheric CO 2 , and Antarctic land surface temperature. The orbital variables (i) to (iii) are exogenous to humanity and are included in the ice-age models. There are additional exogenous variables which may also influence the Earth’s climate and are omitted, such as: (iv) variations in the Sun’s radiation output, (v) volcanic eruption particles in the atmosphere and ice cover, and (vi) changes in the magnetic poles. The new ice-age models of the present paper consider error term specifications that allow for fat-tails and heteroskedasticity, to control for the omitted variables (iv) to (vi).
The main objective of this work is to study the robustness of the results of Castle and Hendry (2020). We improve the model specification of those authors using more general score-driven updates of location and scale in the ice-age model. Score-driven time series models are introduced in the works of Creal et al. (2008) and Harvey and Chakravarty (2008). Those authors name the score-driven models as: generalized autoregressive score (GAS) and dynamic conditional score (DCS) models, respectively. Score-driven models are observation-driven models (Cox 1981), in which the filters are updated using the scaled conditional score functions of the log-likelihood (LL) of the dependent variables.
The use of the score-driven ice-age models is motivated as follows: (i) The updating mechanisms of those models are generalizations of those of the classical time series models such as: ARMA (autoregressive moving average) (Box and Jenkins 1970), GARCH (generalized autoregressive conditional heteroskedasticity) (Engle 2002; Bollerslev 1986), and VARMA (vector ARMA) (Tiao and Tsay 1989). Hence, the score-driven ice-age models of the present paper are extensions of the VAR with explanatory variables (VARX) model of Castle and Hendry (2020). (ii) Score-driven models are robust to outliers and missing observations (Harvey 2013). Hence, the score-driven ice-age models may control for the effects of the omitted exogenous variables. (iii) A score-driven update, asymptotically at the true values of parameters, locally reduces the Kullback–Leibler distance between the true and the estimates values of the score-driven filter in every step, and only score-driven models have this property (Blasques et al. 2015). Hence, score-driven filters use an information-theoretically optimal updating mechanism. We note that the linear updating mechanisms of ARMA and VARMA, and the quadratic updating mechanism of GARCH are optimal from an information-theoretic perspective only if the data generating process (DGP) has a normal distribution. (iv) Under a correct model specification, asymptotically, and at the true values of parameters, the updating terms of the score-driven models are martingale difference sequences (MDS) (Harvey 2013). For the multivariate t-distribution of the present paper, under the same conditions, the same updating terms are white noise processes, due to the well-defined covariance matrix of errors. In the empirical application, we perform model specification tests which support the use of the most general score-driven ice-age model of this paper.
These advantages of the score-driven models motivate their application to climate data. We compare the statistical and 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.
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.

2. Climate Econometrics

2.1. 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 ice-age 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:
y t = μ t + v t = γ 0 + Γ 1 y t 1 + Γ 2 z t + Γ 3 z t 1 + v t
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 t 2 , Ob t 2 , Pr t 2 )
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 ice-age 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 ice-age models of the present paper, because those models are robust to extreme observations (Harvey 2013).
The GUM estimates reported in the work of Castle and Hendry (2020, p. 104) impose restrictions on the matrices Γ 1 , Γ 2 , and Γ 3 . According to those restrictions, the following elements of Γ 1 are not restricted to zero: Γ 1 , 1 , 1 , Γ 1 , 1 , 3 , Γ 1 , 2 , 2 , Γ 1 , 2 , 3 , Γ 1 , 3 , 2 , and Γ 1 , 3 , 3 . Additionally, the following elements of Γ 2 are not restricted to zero: Γ 2 , 1 , 1 , Γ 2 , 1 , 4 , Γ 2 , 1 , 5 , Γ 2 , 2 , 1 , Γ 2 , 2 , 7 , Γ 2 , 3 , 1 , Γ 2 , 3 , 4 , and Γ 2 , 3 , 5 . Furthermore, the following elements of Γ 3 are also not restricted to zero: Γ 3 , 1 , 1 , Γ 3 , 1 , 2 , Γ 3 , 1 , 4 , Γ 3 , 2 , 1 , Γ 3 , 2 , 2 , Γ 3 , 2 , 4 , and Γ 3 , 3 , 4 . The benchmark ice-age model is estimated using the maximum likelihood (ML) method.
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:
y t + j ϵ t = Γ 1 j Ω for j = 0 , ,
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 matrix K ˜ of independent N ( 0 , 1 ) numbers is simulated. (ii) The QR decomposition of K ˜ is performed, and the resulting matrices are denoted as Q ˜ (orthogonal matrix) and R ˜ (upper triangular matrix). (iii) We define Ω ˜ Ω × Q ˜ , for which the estimates of Ω are used. For the IRFs, 3 million simulations of K ˜ 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).

2.2. 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.

2.2.1. 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:
y t = μ t + v t
μ t = γ 0 + Γ 1 μ t 1 + Γ 2 z t + Γ 3 z t 1 + Ψ u t 1
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:
ln f ( y t | F t 1 ; Θ ) = ln Γ ν + 3 2 ln Γ ν 2 3 2 ln ( π ν )
1 2 ln | Σ | ν + 3 2 ln 1 + v t Σ 1 v t ν
where v t = y t μ t , Θ = ( Θ 1 , , Θ S ) 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):
ln f ( y t | F t 1 ; Θ ) μ t = ν + 3 ν Σ 1 × 1 + v t Σ 1 v t ν 1 v t ν + 3 ν Σ 1 × u t
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:
Var ( u t ) = E ln f ( y t | F t 1 ; Θ ) μ t × ln f ( y t | F t 1 ; Θ ) μ t = ν + 3 ν + 5 × Σ 1
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:
y t = γ 0 + Γ 1 y t 1 + Γ 2 z t + Γ 3 z t 1 + ( Ψ Γ 1 ) v t 1 + v t
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. (2008, 2011, 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. (2020, 2021a, 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:
Var ( v t ) = Σ × ν ν 2 = ν ν 2 1 / 2 × Ω Ω × ν ν 2 1 / 2
Based on that, the following multivariate i.i.d. structural-form error term ϵ t is introduced:
v t = ν ν 2 1 / 2 Ω × ϵ t
where E ( ϵ t ) = 0 , Var ( ϵ t ) = I 3 and ϵ t t 3 [ 0 , I 3 × ( ν 2 ) / ν , ν ] .
Furthermore, by substituting Equation (11) into Equation (7), u t as a function of the structural-form error term is:
u t = [ ( ν 2 ) ν ] 1 / 2 Ω ϵ t ν 2 + ϵ t ϵ t .
Second, from Equations (4) and (5), the nonlinear MA() representation of y t is:
y t = v t + j = 0 Γ 1 j γ 0 + Γ 1 j Γ 2 z t j + Γ 1 j Γ 3 z t 1 j + Γ 1 j Ψ u t 1 j
= ν ν 2 1 / 2 Ω × ϵ t
+ j = 0 ϕ j γ 0 + Γ 1 j Γ 2 z t j + Γ 1 j Γ 3 z t 1 j + Γ 1 j Ψ [ ( ν 2 ) ν ] 1 / 2 Ω ϵ t 1 j ν 2 + ϵ t 1 j ϵ t 1 j
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:
y t ϵ t = ν ν 2 1 / 2 Ω
y t + j ϵ t = Γ 1 j 1 Ψ [ ( ν 2 ) ν ] 1 / 2 Ω D ˜ t for j = 1 , ,
where
D ˜ t = ν 2 + ϵ t ϵ t 2 ϵ 1 , t 2 ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 1 , t ϵ 2 , t ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 1 , t ϵ 3 , t ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 2 , t ϵ 1 , t ( ν 2 + ϵ t ϵ t ) 2 ν 2 + ϵ t ϵ t 2 ϵ 2 , t 2 ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 2 , t ϵ 3 , t ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 3 , t ϵ 1 , t ( ν 2 + ϵ t ϵ t ) 2 2 ϵ 3 , t ϵ 2 , t ( ν 2 + ϵ t ϵ t ) 2 ν 2 + ϵ t ϵ t 2 ϵ 3 , t 2 ( ν 2 + ϵ t ϵ t ) 2
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 to D ˜ t . In the empirical application of the present paper, we replace D ˜ 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:
y t ϵ t = Ω
y t + j ϵ t = Γ 1 j 1 Ψ Ω for j = 1 , ,
The IRFs of this section are identified using the procedure of Rubio-Ramirez et al. (2010).

2.2.2. 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:
y t = μ t + v t
μ t = γ 0 + Γ 1 μ t 1 + Γ 2 z t + Γ 3 z t 1 + Ψ u t 1
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.
The reduced-form error term v t | ( F t 1 ; Θ ) t 3 ( 0 , Σ t , ν ) has a multivariate conditional t-distribution, where degrees of freedom ν > 2 , the scale matrix is Σ t D t R D t , where D t ( 3 × 3 ) 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: R = Δ 1 Q Δ 1 Δ 1 Ω Ω Δ 1 , where Δ ( 3 × 3 ) 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:
D t = exp ( λ 1 , t ) 0 0 0 exp ( λ 2 , t ) 0 0 0 exp ( λ 3 , t )
We specify the filters λ i , t in Equation (21) as follows:
λ i , t = ω i + β i λ i , t + α i e i , t 1 + α i * sgn ( v i , t 1 ) ( e i , t 1 + 1 )
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:
ln f ( y t | F t 1 ; Θ ) = ln Γ ν + 3 2 ln Γ ν 2 3 2 ln ( π ν )
1 2 ln | Σ t | ν + 3 2 ln 1 + v t Σ t 1 v t ν
where v t = y t μ t , F t 1 σ ( y 1 , , y t 1 , z 1 , , z t , μ 1 , λ 1 , 1 , λ 2 , 1 , λ 3 , 1 ) , Θ = ( Θ 1 , , Θ S ) is the vector of constant parameters, which includes the elements of γ 0 , Γ 1 , Γ 2 , Γ 3 , Ψ , Ω , ω 1 , ω 2 , ω 3 , β 1 , β 2 , β 3 , α 1 , α 2 , α 3 , α 1 * , α 2 * , α 3 * , and ν . The partial derivative of ln f ( y t | F t 1 ; Θ ) with respect to μ t is:
ln f ( y t | F t 1 ; Θ ) μ t = ν + 3 ν Σ t 1 × 1 + v t Σ t 1 v t ν 1 v t ν + 3 ν Σ t 1 × u t
The scaled score function u t is defined in the second equality of Equation (24), where v t is multiplied by [ 1 + ( v t Σ t 1 v t ) / ν ] 1 = ν / ( ν + v t Σ 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 ν < , 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:
Var ( u t | F t 1 ; Θ ) = ν + 3 ν + 5 × Σ t 1
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
ln f i ( y i , t | F t 1 ; Θ ) = ln Γ ν + 1 2 ln Γ ν 2 1 2 ln ( π ν ) λ i , t
ν + 1 2 ln 1 + v i , t 2 ν exp ( 2 λ i , t )
where v i , t = y i , t μ i , t . We define the updating term of the filter λ i , t for i = 1 , 2 , 3 as follows:
e i , t ln f i ( v i , t | F t 1 ; Θ ) λ i , t = ( ν + 1 ) v i , t 2 ν exp ( 2 λ i , t ) + v i , t 2 1
Equations (22) and (27) are the Beta-t-EGARCH with leverage effects model of Harvey and Chakravarty (2008) (see also Creal et al. 2013 and Harvey 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:
Var ( v t | F t 1 ; Θ ) = Σ t × ν ν 2 = ν ν 2 1 / 2 D t Δ 1 Ω × Ω Δ 1 D t ν ν 2 1 / 2
Based on that, the following multivariate i.i.d. structural-form error term ϵ t is introduced:
v t = ν ν 2 1 / 2 D t Δ 1 Ω × ϵ t
where E ( ϵ t ) = 0 , Var ( ϵ t ) = I 3 and ϵ t t 3 [ 0 , I 3 × ( ν 2 ) / ν , ν ] . Furthermore, by substituting Equation (29) into Equation (24), u t as a function of the structural-form error term is:
u t = [ ( ν 2 ) ν ] 1 / 2 D t Δ 1 Ω ϵ t ν 2 + ϵ t ϵ t .
Second, from Equations (19) and (20), the nonlinear MA() representation of y t is:
y t = v t + j = 0 Γ 1 j γ 0 + Γ 1 j Γ 2 z t j + Γ 1 j Γ 3 z t 1 j + Γ 1 j Ψ u t 1 j
= ν ν 2 1 / 2 D t Δ 1 Ω × ϵ t
+ j = 0 ϕ j γ 0 + Γ 1 j Γ 2 z t j + Γ 1 j Γ 3 z t 1 j + Γ 1 j Ψ [ ( ν 2 ) ν ] 1 / 2 D t Δ 1 Ω ϵ t 1 j ν 2 + ϵ t 1 j ϵ t 1 j
The contemporaneous and dynamic effects of the structural-form error term ϵ t , respectively, are:
y t ϵ t = ν ν 2 1 / 2 D t Δ 1 Ω
y t + j ϵ t = Γ 1 j 1 Ψ [ ( ν 2 ) ν ] 1 / 2 D t Δ 1 Ω D ˜ t for j = 1 , ,
where D ˜ t is given by Equation (16). In Equations (32) and (33), the dynamic interaction effects are time-dependent due to D t and D ˜ t . We replace D t and D ˜ 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.

3. Empirical Results

3.1. 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 Figure 1 and Figure 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.

3.2. Estimation Results

In Table 3, the ML parameter estimates 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 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 Table 3 and Table 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 (Table 3):
Ice ^ t = 1.3735 + 0.8549 Ice t 1 0.0208 Temp t 1 + 95.8353 Ec t 47.5937 Ec t Ob t 5.2167 Ec t Pr t 93.5393 Ec t 1 0.3706 Ob t 1 + 46.7753 Ec t 1 Ob t 1
CO ^ 2 , t = 1.8718 + 0.8468 CO 2 , t 1 + 0.0136 Temp t 1 + 13.8095 Ec t + 0.2106 Ob t 2 27.1270 Ec t 1 1.1138 Ob t 1 + 5.6423 Ec t 1 Ob t 1
Temp ^ t = 2.6657 + 0.8587 CO 2 , t 1 + 0.8684 Temp t 1 335.9696 Ec t + 254.2055 Ec t Ob t + 26.6287 Ec t Pr t 111.3537 Ec t 1 Ob t 1
The estimates correspond to the estimates of Castle and Hendry (2020, p. 104). Second, the estimates of Equation (20) are given by (Table 3):
Ice ^ t = 1.1817 + 0.8824 Ice t 1 0.0172 Temp t 1 + 84.6136 Ec t 42.7119 Ec t Ob t 4.8681 Ec t Pr t 83.4944 Ec t 1 0.3307 Ob t 1 + 42.4409 Ec t 1 Ob t 1 + 0.9651 u 1 , t 1 0.0289 u 3 , t 1
CO ^ 2 , t = 1.4068 + 0.8528 CO 2 , t 1 + 0.0122 Temp t 1 + 11.9581 Ec t + 0.1270 Ob t 2 26.3809 Ec t 1 0.7289 Ob t 1 + 6.1482 Ec t 1 Ob t 1 + 1.3943 u 2 , t 1 + 0.0166 u 3 , t 1
Temp ^ t = 0.6955 + 0.1382 CO 2 , t 1 + 0.9377 Temp t 1 272.0190 Ec t + 194.7069 Ec t Ob t + 17.5525 Ec t Pr t 78.4289 Ec t 1 Ob t 1 + 4.7059 u 2 , t 1 + 0.9860 u 3 , t 1
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 Figure 3, Figure 4, Figure 5 and Figure 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. (2004), Jouzel et al. (2007), Lüthi et al. (2008), Qin and Buehler (2012), Bronselaer et al. (2018), Wadham et al. (2019), and Castle and Hendry (2020). The IRF estimates are persistent and are consistent with the estimates of the long-run solutions of Equation (1) reported in the work of Castle and Hendry (2020).
By comparing the IRF estimates of the benchmark ice-age model with those of the score-driven ice-age models, for several panels of Figure 3, Figure 4, Figure 5 and Figure 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 ) (Figure 3, Figure 4, Figure 5 and Figure 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) (Figure 3, Figure 4, Figure 5 and Figure 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 (Figure 3, Figure 4, Figure 5 and Figure 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 (Figure 3, Figure 4, Figure 5 and Figure 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 E ^ ( λ i , t ) = ω ^ i / ( 1 β ^ i ) for i = 1 , 2 , 3 , respectively. The figure indicates that extreme values of ϵ t are discounted by u t . This supports the outlier-robustness of the score-driven ice-age models of our paper.

3.3. 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-of-sample 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 Figure 8 and Figure 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.

4. 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.

Author Contributions

Conceptualization, S.B. and A.E.; formal analysis, S.B. and A.E.; writing—original draft preparation, S.B. and A.E.; writing—review and editing, S.B. and A.E. All authors have read and agreed to the published version of the manuscript.

Funding

Blazsek acknowledges funding from Universidad Francisco Marroquín. Escribano acknowledges funding from Ministerio de Economía, Industria y Competitividad (ECO2016-00105-001 and MDM 2014-0431), Comunidad de Madrid (MadEco-CM S2015/HUM-3444), and Agencia Estatal de Investigación (2019/00419/001).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

We greatly appreciate that Jennifer L. Castle and David F. Hendry provided us with the dataset of their work (i.e., Castle and Hendry 2020, chp. 6). Data and computer codes are available from the authors of the present paper upon request.

Acknowledgments

All views expressed in this paper are strictly those of the authors and do not reflect the official views or position of any institution with which they may be affiliated or of anyone at those institutions. This paper was presented at the Annual Conference of the Hungarian Economic Society, 21–22 December 2021, Budapest. The authors are thankful for the comments and suggestions of Matthew Copley, Zsolt Darvas, David F. Hendry, and Miklós Vörös. All remaining errors are our own.

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.

References

  1. Anet, Julien G., E. V. Rozanov, Stefan Muthers, T. Peter, Stefan Brönnimann, F. Arfeuille, J. Beer, A. I. Shapiro, C. C. Raible, F. Steinhilber, and et al. 2013. Impact of a potential 21st century “grand solar minimum” on surface temperatures and stratospheric ozone. Geophysical Research Letters 40: 4420–25. [Google Scholar] [CrossRef] [Green Version]
  2. Archer, David, Pamela Martin, Bruce Buffett, Victor Brovkin, Stefan Rahmstorf, and Andrey Ganopolski. 2004. The importance of ocean temperature to global biogeochemistry. Earth and Planetary Science Letters 222: 333–48. [Google Scholar] [CrossRef]
  3. Blasques, Francisco, Janneke van Brummelen, Siem Jan Koopman, and Andre Lucas. 2021. Maximum likelihood estimation for score-driven models. Journal of Econometrics. [Google Scholar] [CrossRef]
  4. Blasques, Francisco, Siem Jan Koopman, and Andre Lucas. 2015. Information-theoretic optimality of observation-driven time series models for continuous responses. Biometrika 102: 325–43. [Google Scholar] [CrossRef]
  5. Blazsek, Szabolcs, Alvaro Escribano, and Adrian Licht. 2020. Identification of seasonal effects in impulse responses using score-driven multivariate location models. Journal of Econometric Methods 10: 53–66. [Google Scholar] [CrossRef]
  6. Blazsek, Szabolcs, Alvaro Escribano, and Adrian Licht. 2021a. Co-integration with score-driven models: An application to US real GDP growth, US inflation rate, and effective federal funds rate. Macroeconomic Dynamics, 1–21. [Google Scholar] [CrossRef]
  7. Blazsek, Szabolcs, Alvaro Escribano, and Adrian Licht. 2021b. Multivariate Markov-switching score-driven models: An application to the global crude oil market. Studies in Nonlinear Dynamics & Econometrics. [Google Scholar] [CrossRef]
  8. Bollerslev, Tim. 1986. Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31: 307–27. [Google Scholar] [CrossRef] [Green Version]
  9. Box, George E. P., and Gwilym M. Jenkins. 1970. Time Series Analysis, Forecasting and Control. San Francisco: Holden-Day. [Google Scholar]
  10. Bronselaer, Ben, Michael Winton, Stephen M. Griffies, William J. Hurlin, Keith B. Rodgers, Olga V. Sergienko, Roland J. Stouffer, and Joellen L. Russell. 2018. Change in future climate due to Antarctic meltwater. Nature 564: 53–58. [Google Scholar] [CrossRef]
  11. Castle, Jennifer, and David F. Hendry. 2020. Climate econometrics: An overview. Foundations and Trends in Econometrics 10: 145–322. [Google Scholar] [CrossRef]
  12. Cox, David R. 1981. Statistical analysis of time series: Some recent developments. Scandinavian Journal of Statistics 8: 93–115. [Google Scholar]
  13. Creal, Drew, Siem Jan Koopman, and Andre Lucas. 2008. A General Framework for Observation Driven Time-Varying Parameter Models. Tinbergen Institute Discussion Paper 08-108/4. Available online: https://www.tinbergen.nl/discussion-paper/2649/08-108-4-a-general-framework-for-observation-driven-timevarying-parameter-models (accessed on 25 December 2021).
  14. Creal, Drew, Siem Jan Koopman, and Andre Lucas. 2011. A dynamic multivariate heavy-tailed model for time-varying volatilities and correlations. Journal of Business & Economic Statistics 29: 552–63. [Google Scholar] [CrossRef] [Green Version]
  15. Creal, Drew, Siem Jan Koopman, and Andre Lucas. 2013. Generalized autoregressive score models with applications. Journal of Applied Econometrics 28: 777–95. [Google Scholar] [CrossRef] [Green Version]
  16. Engle, Robert. 2002. Dynamic conditional correlation: A simple class of multivariate generalized autoregressive conditional heteroskedasticity models. Journal of Business & Economic Statistics 20: 339–51. [Google Scholar] [CrossRef]
  17. Feulner, Georg, and Stefan Rahmstorf. 2010. On the effect of a new grand minimum of solar activity on the future climate on Earth. Geophysical Research Letters 37. [Google Scholar] [CrossRef] [Green Version]
  18. Harvey, Andrew C. 2013. Dynamic Models for Volatility and Heavy Tails: With Applications to Financial and Economic Time Series. Econometric Society Monographs. Cambridge: Cambridge University Press. [Google Scholar]
  19. Harvey, Andrew C., and Tirthankar Chakravarty. 2008. Beta-t-(E)GARCH. Cambridge Working Papers in Economics 0840. Cambridge: Faculty of Economics, University of Cambridge, Available online: http://www.econ.cam.ac.uk/research/repec/cam/pdf/cwpe0840.pdf (accessed on 25 December 2021).
  20. Ineson, Sarah, Amanda C. Maycock, Lesley J. Gray, Adam A. Scaife, Nick J. Dunstone, Jerald W. Harder, Jeff R. Knight, Mike Lockwood, James C. Manners, and Richard A. Wood. 2015. Regional climate impacts of a possible future grand solar minimum. Nature Communications 6: 443. [Google Scholar] [CrossRef]
  21. Intergovernmental Panel on Climate Change. 2021. Sixth Assessment Report. Available online: https://www.ipcc.ch/report/ar6/wg1/#SPM (accessed on 25 December 2021).
  22. Jones, Gareth S., Mike Lockwood, and Peter A. Stott. 2012. What influence will future solar activity changes over the 21st century have on projected global near-surface temperature changes? Journal of Geographical Research 117. [Google Scholar] [CrossRef]
  23. Jouzel, J., V. Masson-Delmotte, O. Cattani, G. Dreyfus, S. Falourd, and G. E. Hoffmann. 2007. Orbital and millennial Antarctic climate variability over the past 800,000 years. Science 317: 793–97. [Google Scholar] [CrossRef] [Green Version]
  24. Kibria, B. M. Golam, and Anwar H. Joarder. 2006. A short review of multivariate t-distribution. Journal of Statistical Research 40: 59–72. [Google Scholar]
  25. Kirchner, Ingo, Georgiy L. Stenchikov, Hans-F. Graf, Alan Robock, and Juan Carlos Antuna. 1999. Climate model simulation of winter warming and summer cooling following the 1991 Mount Pinatubo volcanic eruption. Journal of Geophysical Research 104: 19039–55. [Google Scholar] [CrossRef]
  26. Lisiecki, Lorraine E., and Maureen E. Raymo. 2005. A pliocene-pleistocene stack of 57 globally distributed Benthic δ18O records. Paleoceanography 20. [Google Scholar] [CrossRef] [Green Version]
  27. Ljung, Greta M., and George E. P. Box. 1978. On a measure of lack of fit in time-series models. Biometrika 65: 297–303. [Google Scholar] [CrossRef]
  28. Lüthi, Dieter, Marthine Le Floch, Bernhard Bereiter, Thomas Blunier, Jean-Marc Barnola, Urs Steigenhaler, Dominique Raynaud, Jean Jouzel, Hubertus Fischer, Kenji Kawamura, and et al. 2008. High-resolution carbon dioxide concentration record 650,000–800,000 years before present. Nature 453. [Google Scholar] [CrossRef]
  29. Lütkepohl, Helmut. 2005. New Introduction to Multivariate Time Series Analysis. Berlin and Heidelberg: Springer. [Google Scholar]
  30. Maycock, A. C., S. Ineson, L. J. Gray, A. A. Scaife, J. A. Anstey, M. Lockwood, N. Butchart, S. C. Hardiman, D. M. Mitchell, and S. M. Osprey. 2015. Possible impacts of a future grand solar minimum on climate: Stratospheric and global circulation changes. JGR Atmospheres 120: 9043–58. [Google Scholar] [CrossRef] [Green Version]
  31. Meehl, Gerald A., Julie M. Arblaster, and Daniel R. Marsh. 2013. Could a future “Grand Solar Minimum” like the Maunder Minimum stop global warming? Geographical Research Letters 40: 1789–93. [Google Scholar] [CrossRef]
  32. NASA. 2021. Flip flop: Why Variations in Earth’s Magnetic Field Aren’t Causing Today’S Climate Change. Available online: https://climate.nasa.gov/ask-nasa-climate/3104/flip-flop-why-variations-in-earths-magnetic-field-arent-causing-todays-climate-change/ (accessed on 25 December 2021).
  33. Paillard, Didier, Laurent D. Labeyrie, and Pascal Yiou. 1996. Macintosh program performs time-series analysis. Eos Transactions AGU 77: 379–79. [Google Scholar] [CrossRef]
  34. Qin, Zhao, and Markus J. Buehler. 2012. Carbon dioxide enhances fragility of ice crystals. Journal of Physics D: Applied Physics 45. [Google Scholar] [CrossRef]
  35. Rubio-Ramirez, Juan F., Daniel Waggoner, and Tao Zha. 2010. Structural vector autoregressions: Theory for identification and algorithms for inference. Review of Economic Studies 77: 665–96. [Google Scholar] [CrossRef] [Green Version]
  36. Ruddiman, William. 2005. Plows, Plagues and Petroleum: How Humans Took Control of the Climate. Princeton: Princeton University Press. [Google Scholar]
  37. Steinhilber, F., J. A. Abreu, and J. Beer. 2008. Solar modulation during the Holocene. Astrophysics and Space Sciences Transactions 4: 1–6. [Google Scholar] [CrossRef] [Green Version]
  38. Stenchikov, Georgiy L., Ingo Kirchner, Alan Robock, Hans-F. Graf, Juan Carlos Antuna, R. G. Grainger, Alyn Lambert, and Larry Thomason. 1998. Radiative Forcing from the 1991 Mount Pinatubo volcanic eruption. Journal of Geophysical Research 103: 13837–57. [Google Scholar] [CrossRef] [Green Version]
  39. Tiao, George C., and Ruey S. Tsay. 1989. Model specification in multivariate time series. Journal of the Royal Statistical Society 51: 157–213. [Google Scholar] [CrossRef]
  40. Wadham, J. L., J. R. Hawkings, L. Tarasov, L. J. Gregoire, R. G. M. Spencer, M. Gutjahr, A. Ridgwell, and K. E. Kohfeld. 2019. Ice sheets matter for the global carbon cycle. Nature Communications 10: 3567. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. White, Halbert. 2001. Asymptotic Theory for Econometricians. revised edition. San Diego: Academic Press. [Google Scholar]
Figure 1. Evolution of Ice t , CO 2 , t , and Temp t from 798 thousand years ago to 1 thousand years ago.
Figure 1. Evolution of Ice t , CO 2 , t , and Temp t from 798 thousand years ago to 1 thousand years ago.
Econometrics 10 00009 g001
Figure 2. Evolution of Ec t , Ob t , and Pr t from 798 thousand years ago to 1 thousand years ago.
Figure 2. Evolution of Ec t , Ob t , and Pr t from 798 thousand years ago to 1 thousand years ago.
Econometrics 10 00009 g002
Figure 3. IRFs up to 20 lags for the benchmark ice-age model. Notes: The confidence interval is mean ± one standard deviation that is estimated for 2302 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Figure 3. IRFs up to 20 lags for the benchmark ice-age model. Notes: The confidence interval is mean ± one standard deviation that is estimated for 2302 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Econometrics 10 00009 g003
Figure 4. IRFs up to 20 lags for the score-driven ice-age model for the normal distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 1771 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Figure 4. IRFs up to 20 lags for the score-driven ice-age model for the normal distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 1771 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Econometrics 10 00009 g004
Figure 5. IRFs up to 20 lags for the score-driven homoskedastic ice-age model for the t-distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 2289 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Figure 5. IRFs up to 20 lags for the score-driven homoskedastic ice-age model for the t-distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 2289 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Econometrics 10 00009 g005
Figure 6. IRFs up to 20 lags for the score-driven heteroskedastic ice-age model for the t-distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 1698 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Figure 6. IRFs up to 20 lags for the score-driven heteroskedastic ice-age model for the t-distribution. Notes: The confidence interval is mean ± one standard deviation that is estimated for 1698 out of the 3 million simulations under the restrictions of Table 1. Ice uses the δ 18 O proxy; 1 CO 2 is 780 gigatonnes; 1 Temp is 1 °C.
Econometrics 10 00009 g006
Figure 7. Robustness of the scaled score function to extreme values. Note: ϵ 3 , t = 0 is assumed for this figure.
Figure 7. Robustness of the scaled score function to extreme values. Note: ϵ 3 , t = 0 is assumed for this figure.
Econometrics 10 00009 g007
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.
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.
Econometrics 10 00009 g008
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.
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.
Econometrics 10 00009 g009
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.
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.
Econometrics 10 00009 g010
Table 1. Sign restrictions on contemporaneous impact responses.
Table 1. Sign restrictions on contemporaneous impact responses.
Ice t Shock CO 2 , t Shock Temp t Shock
Ice t +
CO 2 , t ++
Temp t ++
Table 2. Descriptive statistics.
Table 2. Descriptive statistics.
(a) Dependent Variables Ice t CO 2 , t Temp t
VariableIce volumeAtmospheric CO 2 Antarctic-based land surface temperature
Start date798 thousand years ago798 thousand years ago798 thousand years ago
End date1 thousand years ago1 thousand years ago1 thousand years ago
Data frequency1 thousand years1 thousand years1 thousand years
MeasurementBased on the δ 18 O proxy1 unit = 780 gigatonnes of CO 2 1 unit = 1 Celsius degree
Data sourceLisiecki and Raymo (2005)Lüthi et al. (2008)Jouzel et al. (2007)
Sample size798798798
Minimum 3.1000 1.7269 10.2530
Maximum 5.0800 2.9500 3.7662
Mean 4.1707 2.2382 5.2892
Standard deviation 0.4467 0.2546 2.9009
(b) Explanatory Variables Ec t Ob t Pr t
VariableEccentricity of the Earth’s orbitObliquityPrecession of the equinox
Start date798 thousand years ago798 thousand years ago798 thousand years ago
End date1 thousand years ago1 thousand years ago1 thousand years ago
Data frequency1 thousand years1 thousand years1 thousand years
MeasurementPeriodicity deriving from thePeriodicity deriving from thePeriodicity deriving from the
changing non-circularity of the Earth’s orbitchanges in the tilt of the Earth’s rotational axisprecession of the equinox
(zero denotes circularity).relative to the ecliptic (1 unit = 10 degrees).(1 unit = 1 degree).
Data sourcePaillard et al. (1996)Paillard et al. (1996)Paillard et al. (1996)
Sample size798798798
Minimum 0.0042 2.2076 0.0008
Maximum 0.0500 2.4455 0.3593
Mean 0.0271 2.3342 0.1802
Standard deviation 0.0119 0.0591 0.1039
Table 3. In-sample ML parameter estimates.
Table 3. In-sample ML parameter estimates.
BenchmarkScore-Driven HomoskedasticScore-Driven HomoskedasticScore-Driven Heteroskedastic
Ice-Age ModelGaussian Ice-Age Modelt Ice-Age Modelt Ice-Age Model
γ 0 , 1 1.3735 *** ( 0.3009 ) γ 0 , 1 1.2697 *** ( 0.2884 ) γ 0 , 1 1.2745 *** ( 0.2798 ) γ 0 , 1 1.1817 *** ( 0.2812 ) ω 1 1.1773 *** ( 0.2845 )
γ 0 , 2 1.8718 *** ( 0.3023 ) γ 0 , 2 1.6589 *** ( 0.3552 ) γ 0 , 2 1.8197 *** ( 0.3579 ) γ 0 , 2 1.4068 *** ( 0.3259 ) ω 2 0.7409 *** ( 0.1416 )
γ 0 , 3 2.6657 *** ( 0.6905 ) γ 0 , 3 1.4256 + ( 0.9049 ) γ 0 , 3 1.4366 + ( 0.9252 ) γ 0 , 3 0.6955 ( 0.6592 ) ω 3 0.0455 *** ( 0.0139 )
Γ 1 , 1 , 1 0.8549 *** ( 0.0144 ) Γ 1 , 1 , 1 0.8733 *** ( 0.0155 ) Γ 1 , 1 , 1 0.8738 *** ( 0.0150 ) Γ 1 , 1 , 1 0.8824 *** ( 0.0162 ) β 1 0.5266 *** ( 0.1144 )
Γ 1 , 1 , 3 0.0208 *** ( 0.0020 ) Γ 1 , 1 , 3 0.0175 *** ( 0.0024 ) Γ 1 , 1 , 3 0.0175 *** ( 0.0024 ) Γ 1 , 1 , 3 0.0172 *** ( 0.0025 ) β 2 0.7628 *** ( 0.0446 )
Γ 1 , 2 , 2 0.8468 *** ( 0.0176 ) Γ 1 , 2 , 2 0.8410 *** ( 0.0287 ) Γ 1 , 2 , 2 0.8291 *** ( 0.0286 ) Γ 1 , 2 , 2 0.8528 *** ( 0.0231 ) β 3 0.8796 *** ( 0.0286 )
Γ 1 , 2 , 3 0.0136 *** ( 0.0016 ) Γ 1 , 2 , 3 0.0125 *** ( 0.0026 ) Γ 1 , 2 , 3 0.0137 *** ( 0.0026 ) Γ 1 , 2 , 3 0.0122 *** ( 0.0022 ) α 1 0.0500 *** ( 0.0192 )
Γ 1 , 3 , 2 0.8587 *** ( 0.2556 ) Γ 1 , 3 , 2 0.3714 ( 0.3333 ) Γ 1 , 3 , 2 0.3771 ( 0.3396 ) Γ 1 , 3 , 2 0.1382 ( 0.2458 ) α 2 0.0930 *** ( 0.0154 )
Γ 1 , 3 , 3 0.8684 *** ( 0.0231 ) Γ 1 , 3 , 3 0.9008 *** ( 0.0297 ) Γ 1 , 3 , 3 0.8995 *** ( 0.0304 ) Γ 1 , 3 , 3 0.9377 *** ( 0.0228 ) α 3 0.1070 *** ( 0.0172 )
Γ 2 , 1 , 1 95.8353 *** ( 30.3012 ) Γ 2 , 1 , 1 88.7529 *** ( 29.2107 ) Γ 2 , 1 , 1 88.1550 *** ( 28.6783 ) Γ 2 , 1 , 1 84.6136 *** ( 28.7107 ) α 1 * 0.0513 *** ( 0.0163 )
Γ 2 , 1 , 4 47.5937 *** ( 12.4804 ) Γ 2 , 1 , 4 45.5349 *** ( 12.0258 ) Γ 2 , 1 , 4 45.9952 *** ( 11.8462 ) Γ 2 , 1 , 4 42.7119 *** ( 12.0538 ) α 2 * 0.0199 + ( 0.0124 )
Γ 2 , 1 , 5 5.2167 *** ( 1.0663 ) Γ 2 , 1 , 5 5.4569 *** ( 1.0321 ) Γ 2 , 1 , 5 5.6132 *** ( 1.0064 ) Γ 2 , 1 , 5 4.8681 *** ( 0.9172 ) α 3 * 0.0094 ( 0.0123 )
Γ 2 , 2 , 1 13.8095 *** ( 3.8487 ) Γ 2 , 2 , 1 14.5128 *** ( 4.8910 ) Γ 2 , 2 , 1 14.8536 *** ( 4.8133 ) Γ 2 , 2 , 1 11.9581 *** ( 3.5543 )
Γ 2 , 2 , 7 0.2106 *** ( 0.0446 ) Γ 2 , 2 , 7 0.1688 *** ( 0.0518 ) Γ 2 , 2 , 7 0.1876 *** ( 0.0515 ) Γ 2 , 2 , 7 0.1270 *** ( 0.0477 )
Γ 2 , 3 , 1 335.9696 *** ( 38.6135 ) Γ 2 , 3 , 1 311.6342 *** ( 43.7423 ) Γ 2 , 3 , 1 326.4885 *** ( 44.2197 ) Γ 2 , 3 , 1 272.0190 *** ( 32.9400 )
Γ 2 , 3 , 4 254.2055 *** ( 28.4900 ) Γ 2 , 3 , 4 240.1468 *** ( 30.8981 ) Γ 2 , 3 , 4 258.2711 *** ( 30.4240 ) Γ 2 , 3 , 4 194.7069 *** ( 22.0173 )
Γ 2 , 3 , 5 26.6287 *** ( 8.3118 ) Γ 2 , 3 , 5 25.0557 *** ( 7.8012 ) Γ 2 , 3 , 5 28.2720 *** ( 7.6689 ) Γ 2 , 3 , 5 17.5525 *** ( 5.9466 )
Γ 3 , 1 , 1 93.5393 *** ( 31.3704 ) Γ 3 , 1 , 1 83.6310 *** ( 30.7764 ) Γ 3 , 1 , 1 81.9865 *** ( 30.4458 ) Γ 3 , 1 , 1 83.4944 *** ( 30.2589 )
Γ 3 , 1 , 2 0.3706 *** ( 0.1233 ) Γ 3 , 1 , 2 0.3519 *** ( 0.1164 ) Γ 3 , 1 , 2 0.3544 *** ( 0.1131 ) Γ 3 , 1 , 2 0.3307 *** ( 0.1135 )
Γ 3 , 1 , 4 46.7753 *** ( 12.9406 ) Γ 3 , 1 , 4 43.5364 *** ( 12.6666 ) Γ 3 , 1 , 4 43.5483 *** ( 12.5889 ) Γ 3 , 1 , 4 42.4409 *** ( 12.6434 )
Γ 3 , 2 , 1 27.1270 *** ( 6.7356 ) Γ 3 , 2 , 1 29.1279 *** ( 7.5563 ) Γ 3 , 2 , 1 30.6263 *** ( 7.4143 ) Γ 3 , 2 , 1 26.3809 *** ( 6.1789 )
Γ 3 , 2 , 2 1.1138 *** ( 0.2191 ) Γ 3 , 2 , 2 0.9235 *** ( 0.2529 ) Γ 3 , 2 , 2 1.0217 *** ( 0.2528 ) Γ 3 , 2 , 2 0.7289 *** ( 0.2339 )
Γ 3 , 2 , 4 5.6423 ** ( 2.2875 ) Γ 3 , 2 , 4 6.2517 ** ( 2.4809 ) Γ 3 , 2 , 4 6.7580 *** ( 2.4337 ) Γ 3 , 2 , 4 6.1482 *** ( 2.1916 )
Γ 3 , 3 , 4 111.3537 *** ( 24.0465 ) Γ 3 , 3 , 4 107.2597 *** ( 26.1398 ) Γ 3 , 3 , 4 119.1342 *** ( 26.2451 ) Γ 3 , 3 , 4 78.4289 *** ( 16.6672 )
Ψ 1 , 1 , 1 NA Ψ 1 , 1 , 1 0.7982 *** ( 0.0343 ) Ψ 1 , 1 , 1 0.9020 *** ( 0.0438 ) Ψ 1 , 1 , 1 0.9651 *** ( 0.0551 )
Ψ 1 , 1 , 3 NA Ψ 1 , 1 , 3 0.0295 *** ( 0.0041 ) Ψ 1 , 1 , 3 0.0325 *** ( 0.0050 ) Ψ 1 , 1 , 3 0.0289 *** ( 0.0047 )
Ψ 1 , 2 , 2 NA Ψ 1 , 2 , 2 1.1351 *** ( 0.0300 ) Ψ 1 , 2 , 2 1.3418 *** ( 0.0545 ) Ψ 1 , 2 , 2 1.3943 *** ( 0.0599 )
Ψ 1 , 2 , 3 NA Ψ 1 , 2 , 3 0.0197 *** ( 0.0025 ) Ψ 1 , 2 , 3 0.0221 *** ( 0.0030 ) Ψ 1 , 2 , 3 0.0166 *** ( 0.0025 )
Ψ 1 , 3 , 2 NA Ψ 1 , 3 , 2 4.1436 *** ( 0.5376 ) Ψ 1 , 3 , 2 4.9368 *** ( 0.7038 ) Ψ 1 , 3 , 2 4.7059 *** ( 0.6182 )
Ψ 1 , 3 , 3 NA Ψ 1 , 3 , 3 0.8723 *** ( 0.0347 ) Ψ 1 , 3 , 3 1.0084 *** ( 0.0480 ) Ψ 1 , 3 , 3 0.9860 *** ( 0.0538 )
Ω 1 , 1 0.0891 *** ( 0.0019 ) Ω 1 , 1 0.0887 *** ( 0.0019 ) Ω 1 , 1 0.0861 *** ( 0.0021 ) Ω 1 , 1 NA
Ω 2 , 1 0.0102 *** ( 0.0018 ) Ω 2 , 1 0.0085 *** ( 0.0018 ) Ω 2 , 1 0.0096 *** ( 0.0018 ) Ω 2 , 1 0.1698 *** ( 0.0410 )
Ω 2 , 2 0.0524 *** ( 0.0011 ) Ω 2 , 2 0.0492 *** ( 0.0010 ) Ω 2 , 2 0.0475 *** ( 0.0010 ) Ω 2 , 2 NA
Ω 3 , 1 0.1755 *** ( 0.0267 ) Ω 3 , 1 0.1600 *** ( 0.0262 ) Ω 3 , 1 0.1672 *** ( 0.0274 ) Ω 3 , 1 0.2343 *** ( 0.0449 )
Ω 3 , 2 0.3808 *** ( 0.0249 ) Ω 3 , 2 0.3440 *** ( 0.0239 ) Ω 3 , 2 0.3476 *** ( 0.0247 ) Ω 3 , 2 0.5107 *** ( 0.0424 )
Ω 3 , 3 0.6711 *** ( 0.0145 ) Ω 3 , 3 0.6684 *** ( 0.0146 ) Ω 3 , 3 0.6458 *** ( 0.0152 ) Ω 3 , 3 NA
ν NA ν NA ν 38.2860 *** ( 6.1464 ) ν 27.7970 *** ( 5.1586 )
Notes: Not available (NA). Standard errors are reported in parentheses. ***, **, *, and + indicate significance at the 1%, 5%, 10%, and 15% levels, respectively.
Table 4. Model performance and diagnostics for in-sample estimates.
Table 4. Model performance and diagnostics for in-sample estimates.
BenchmarkScore-Driven HomoskedasticScore-Driven HomoskedasticScore-Driven Heteroskedastic
Ice-Age ModelGaussian Ice-Age Modelt Ice-Age Modelt Ice-Age Model
LL 1.5094 1.5801 1.6026 1.7644
AIC 2.9435 3.0699 3.1126 3.4134
BIC 2.7675 2.8587 2.8955 3.1435
HQC 2.8759 2.9887 3.0291 3.3097
C μ 0.9663 0.9453 0.9444 0.9543
C λ , 1 NANANA 0.5266
C λ , 2 NANANA 0.7628
C λ , 3 NANANA 0.8796
LB v 1 , t (p-value) 20.5200 ( 0.8448 ) 17.8671 ( 0.9294 ) 19.1502 ( 0.8934 ) 19.3022 ( 0.8886 )
LB v 2 , t (p-value) 118.9560 *** ( 0.0000 ) 28.4522 ( 0.4407 ) 29.2620 ( 0.3993 ) 30.8129 ( 0.3255 )
LB v 3 , t (p-value) 41.8520 ** ( 0.0448 ) 30.2398 ( 0.3518 ) 33.0351 ( 0.2345 ) 37.1410 ( 0.1158 )
LB ϵ 1 , t (p-value) 20.5200 ( 0.8448 ) 17.867 ( 0.9294 ) 19.1502 ( 0.8934 ) 16.9958 ( 0.9487 )
LB ϵ 2 , t (p-value) 98.1577 *** ( 0.0000 ) 28.6594 ( 0.4300 ) 27.1411 ( 0.5106 ) 21.1773 ( 0.8179 )
LB ϵ 3 , t (p-value) 47.0696 ** ( 0.0135 ) 44.3814 ** ( 0.0255 ) 45.6695 ** ( 0.0189 ) 33.1683 ( 0.2296 )
LB u 1 , t (p-value)NANA 19.9979 ( 0.8645 ) 21.7319 ( 0.7935 )
LB u 2 , t (p-value)NANA 27.8256 ( 0.4737 ) 27.7605 ( 0.4772 )
LB u 3 , t (p-value)NANA 33.5007 ( 0.2178 ) 35.0719 ( 0.1678 )
LB e 1 , t (p-value)NANANA 34.1775 ( 0.1951 )
LB e 2 , t (p-value)NANANA 25.8774 ( 0.5798 )
LB e 3 , t (p-value)NANANA 22.7776 ( 0.7441 )
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.
Table 5. Multi-step ahead forecasts for the last 100 thousand years.
Table 5. Multi-step ahead forecasts for the last 100 thousand years.
Score-Driven Score-Driven
HomoskedasticScore-DrivenScore-Driven HomoskedasticScore-DrivenScore-Driven
BenchmarkGaussianHomoskedasticHeteroskedasticBenchmarkGaussianHomoskedasticHeteroskedastic
Ice-Age ModelIce-Age Modelt Ice-Age Modelt Ice-Age ModelIce-Age ModelIce-Age Modelt Ice-Age Modelt Ice-Age Model
Ice t MSEMSEMSEMSEMAEMAEMAEMAE
last 100,000 years 0 . 0917 0.0982 0.0969 0.1005 0.2388 0.2520 0.2555 0 . 2380
last 90,000 years 0 . 1003 0.1075 0.1058 0.1098 0.2535 0.2680 0.2713 0 . 2515
last 80,000 years 0 . 1081 0.1169 0.1148 0.1205 0 . 2661 0.2830 0.2867 0.2664
last 70,000 years 0 . 1186 0.1284 0.1258 0.1352 0 . 2818 0.3008 0.3044 0.2881
last 60,000 years 0 . 1259 0.1370 0.1317 0.1526 0 . 2848 0.3063 0.3066 0.3097
last 50,000 years 0 . 1416 0.1549 0.1472 0.1765 0 . 3005 0.3265 0.3238 0.3373
last 40,000 years 0 . 1712 0.1868 0.1765 0.2150 0 . 3444 0.3738 0.3689 0.3914
last 30,000 years 0.2148 0.2292 0 . 2134 0.2691 0 . 3946 0.4187 0.4080 0.4471
last 20,000 years 0.3049 0.3164 0 . 2911 0.3767 0.5037 0.5143 0 . 4948 0.5575
last 10,000 years 0.4889 0.5027 0 . 4527 0.6162 0.6935 0.7032 0 . 6667 0.7810
CO 2 , t MSEMSEMSEMSEMAEMAEMAEMAE
last 100,000 years 0 . 0399 0.0432 0.0424 0.0466 0 . 1471 0.1503 0.1506 0.1556
last 90,000 years 0 . 0440 0.0479 0.0470 0.0517 0 . 1580 0.1642 0.1641 0.1709
last 80,000 years 0 . 0460 0.0509 0.0494 0.0550 0 . 1595 0.1675 0.1664 0.1755
last 70,000 years 0 . 0513 0.0568 0.0552 0.0623 0 . 1719 0.1810 0.1797 0.1920
last 60,000 years 0 . 0590 0.0653 0.0634 0.0702 0 . 1900 0.1995 0.1986 0.2064
last 50,000 years 0 . 0692 0.0769 0.0746 0.0831 0 . 2128 0.2254 0.2240 0.2352
last 40,000 years 0 . 0842 0.0935 0.0902 0.1015 0 . 2462 0.2622 0.2598 0.2754
last 30,000 years 0 . 1104 0.1214 0.1165 0.1313 0 . 3089 0.3263 0.3212 0.3404
last 20,000 years 0.1269 0.1338 0 . 1224 0.1467 0.3261 0.3353 0 . 3212 0.3528
last 10,000 years 0.1891 0.1978 0 . 1733 0.2221 0.4280 0.4379 0 . 4088 0.4664
Temp t MSEMSEMSEMSEMAEMAEMAEMAE
last 100,000 years 4 . 1809 4.5136 4.5168 4.8663 1 . 6976 1.7704 1.7906 1.8453
last 90,000 years 4 . 4536 4.8600 4.8533 5.2546 1 . 7714 1.8574 1.8754 1.9406
last 80,000 years 4 . 5747 5.1051 5.0628 5.4916 1 . 7939 1.9107 1.9207 1.9848
last 70,000 years 5 . 0599 5.6824 5.6177 6.0398 1 . 9174 2.0581 2.0621 2.1106
last 60,000 years 5 . 7960 6.5246 6.4482 6.8760 2 . 1167 2.2824 2.2886 2.3143
last 50,000 years 6 . 4533 7.3840 7.2948 7.8500 2 . 2670 2.4787 2.4871 2.5283
last 40,000 years 7 . 3939 8.3779 8.1927 9.0115 2 . 4913 2.7003 2.6894 2.7802
last 30,000 years 8 . 8750 9.7268 9.3292 10.5446 2 . 8029 2.9516 2.8995 3.0520
last 20,000 years 10.0692 10.5303 9 . 5685 11.9260 2.9566 3.0268 2 . 8845 3.2107
last 10,000 years 14.0302 14.6377 12 . 8303 16.8965 3.7069 3.7892 3 . 5384 4.0858
Notes: Mean square error (MSE); mean absolute error (MAE). The lowest loss function values are indicated by bold numbers.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Blazsek, S.; Escribano, A. Robust Estimation and Forecasting of Climate Change Using Score-Driven Ice-Age Models. Econometrics 2022, 10, 9. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics10010009

AMA Style

Blazsek S, Escribano A. Robust Estimation and Forecasting of Climate Change Using Score-Driven Ice-Age Models. Econometrics. 2022; 10(1):9. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics10010009

Chicago/Turabian Style

Blazsek, Szabolcs, and Alvaro Escribano. 2022. "Robust Estimation and Forecasting of Climate Change Using Score-Driven Ice-Age Models" Econometrics 10, no. 1: 9. https://0-doi-org.brum.beds.ac.uk/10.3390/econometrics10010009

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop