# System Identification of the Global Climate Temperature by Output Error Method

Stéphane Victor, Rachid Malti, Abir Mayoufi

► To cite this version:

Stéphane Victor, Rachid Malti, Abir Mayoufi. System Identification of the Global Climate Temperature by Output Error Method. 2023 International Conference on Fractional Differentiation and Its Applications (ICFDA), Mar 2023, Ajman, United Arab Emirates. pp.1-6, 10.1109/ICFDA58234.2023.10153210. hal-04317768

HAL Id: hal-04317768

https://hal.science/hal-04317768v1

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# System Identification of the Global Climate Temperature by Output Error Method

Stéphane Victor Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France

stephane.victor@ims-bordeaux.fr Rachid Malti Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France

rachid.malti@ims-bordeaux.fr Abir Mayoufi Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France

abir.mayoufi@ims-bordeaux.fr

## Abstract

Global warming is a fact that has been announced by the Intergovernmental Panel on Climate Change (IPCC) since the eighties. Modeling the global terrestrial temperature through system identification is a very challenging task. The fractional integration operator is well adapted to model diffusion or propagation phenomena with few parameters and thanks to its long memory property. It is also known that fractional models are well suited for modeling thermal systems and especially for solving the heat equation. Also, it is quite natural to use the fractional operator to model the global climate temperature. To this aim, system identification by using fractional order models has been recently proposed for continuous-time multiple-input single-output (MISO) models. When differentiation orders are assumed known, coefficients are estimated by using the output error method for continuous-time fractional models extended to the MISO case. When unknown, a gradient-based algorithm is proposed for differentiation order estimation. Finally, model parameters are estimated on real input/output data and provide a very accurate goodness fit for the global earth temperature modeling.

## Index Terms

system identification, continuous-time, output error method, multiple-input single-output (MISO) system, order optimization, fractional model, global warming

## I. INTRODUCTION

Forecasting the terrestrial global climate temperature is a priority. The Intergovernmental Panel on Climate Change (IPCC) is involved in modeling the global climate change and has warned of a significative increase of the global terrestrial temperature. Such an increasing is very alarming as important climate disasters occur.

The climate modeling used in the IPCC reports are made from general circulatory models, which are very complex. Here, it is proposed to use objective criteria by means of system identification in automatic control to model the global terrestrial temperature. More precisely, a black-box dynamical model is estimated to approximate the climate system behavior through system identification with fractional models.

Developments on system identification with fractional order models (FOM) have significantly increased over the last 20 years (see [10] for the state of the art): model identification with low order of fractional systems [5], subspace method for state-space identification [21], multivariable non commensurate fractional systems [12], FOM with time delays [16], identification in fractional differential equations [7], identification of FOM through block pulse functions [20], system identification based on polynomial modulating function [6], recursive system identification with long memory prediction error method [22], multiple-inputs single-output system identification with optimal instrumental variable mechanism [23], etc.

Concerning the input-output data that characterize the global climate, instrumental measurements, namely from satellites, are available only since the eighties. To go further before, proxies are used to reconstruct data from different geological sources. The current and future evolutions of the global climate temperature is a burning issue as it raises many questions. In this context, a terrestrial climate model was proposed by De Larminat [4] where the dynamic relationship has been carried out between the inputs (solar activity, volcanic activity, human activity) and the averaged measured temperature.

The contribution is to propose system identification with fractional models to model the global climate change from data measurements.

Recalling briefly MISO fractional systems in section II, section III develops system identification methods that will be used to model the MISO climate system. Then, section IV presents the terrestrial climate system and the input/output data selection. Section V presents the results of MISO system identification applied on measurement data with different model structures. Finally, a conclusion is drawn in section VI.
## II. FRACTIONAL CALCULUS FOR MISO SYSTEMS

### A. Fractional calculus

Differentiating a piecewise continuous function $f(t)$, to an arbitrary order $\nu \in \mathbb{R}^+$ according to Grünwald-Letnikov, is defined by:

$$p^{\nu} f(t) = \lim_{\ell \to 0} \frac{1}{\ell^{\nu}} \sum_{h=0}^{\lfloor \frac{\ell}{h} \rfloor} (-1)^{h} \binom{\nu}{h} f(t - h\ell) \quad (1)$$

where $\lfloor \cdot \rfloor$ is the floor operator and $\binom{\nu}{h}$ is the Newton's binomial generalized to fractional orders as follows:

$$\binom{\nu}{h} = \frac{\Gamma(\nu + 1)}{\Gamma(h + 1)\Gamma(\nu - h + 1)} = \frac{\nu(\nu - 1)\dots(\nu - h + 1)}{h!}.$$

As a consequence, the fractional derivative has a global characterization by considering the whole past of the function.

When numerically evaluating the $\nu$-th fractional derivative, $\ell$ in (1) can be replaced by the sampling time $T$ thus leading to:

$$p^{\nu} f(t) = \frac{1}{T^{\nu}} \sum_{h=0}^{\lfloor \frac{T}{h} \rfloor} (-1)^{h} \binom{\nu}{h} f(t - Th) + \mathcal{O}(T). \quad (2)$$

The Laplace transform of the derivative of order $\nu$ of a causal function $f(f(t) = 0, \forall t \le 0)$ is defined by [17]:

$$\mathcal{L}\{p^{\nu}x(t)\} = s^{\nu}X(s), \quad (3)$$

where $s$ is the Laplace variable.

### B. MISO fractional order systems

A MISO fractional order system is governed by (see Fig.1):

$$S: \begin{cases} y_k(t) = G_k(p)u_k(t) \\ y(t) = \sum_{k=1}^{K} y_k(t) \\ y^*(t_h) = y(t_h) + e(t_h) \end{cases} \quad (4)$$

where $K \in \mathbb{N}$ denotes the number of inputs. Let be $u(t) = [u_1(t), \dots, u_K(t)]$ that contains the uncorrelated input signals and $p$ the differential operator ($p = \frac{d}{dt}$). The input-output data $(u(t), y(t))$ are assumed collected at regular samples and large enough to guarantee the convergence of the estimates to the true ones. The quasi-stationary input signals $\{u_k(t), 0 \le t \le t_{final}, k = 1, \dots, K\}$ are applied to the MISO system (4), and the inputs are assumed persistently exciting and give rise to the output signal $\{y(t), 0 \le t \le t_{final}\}$. The output measurement $y(t)$ is corrupted by a discrete-time white noise $e(t_h)$ and $G_k$ is the fractional operator that relates the input signal $u_k(t)$ to its noise-free output $y_k(t)$:

$$G_k(p, \theta) = \frac{B_k(p, \theta)}{A_k(p, \theta)} = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{\beta_{j,k}}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{\alpha_{i,k}}} \quad (5)$$

where the parameter vector $\theta$ is defined in (10) and, for identifiability purpose, the differentiation orders are ordered positive real numbers:

$$\begin{cases} 0 < \alpha_{1,k} < \alpha_{2,k} < \dots < \alpha_{N_k,k} & \forall k = 1, \dots, K. \\ 0 \le \beta_{0,k} < \beta_{1,k} < \dots < \beta_{M_k,k} & \forall k = 1, \dots, K. \end{cases}$$

Assume that the $G_k(p)$ operators are strictly proper with $\alpha_{N_k,k} > \beta_{M_k,k}, \forall k$, that $B_k$ and $A_k$ are coprime polynomials, and that all transfer functions $G_k(p)$ are asymptotically stable [13].

**Definition 2.1 (S-commensurability, [23]):** A SISO system $G_1$, for a given model structure (5), characterized by the number of terms $M_1$ and $N_1$, has a structured-commensurate (or S-commensurate) order

$$\nu = \alpha_{1,1}, \quad (6)$$

if all its differentiation orders are successive integer multiples of $\nu$, namely:

$$G_1(p) = \frac{\sum_{j=0}^{M_1} b_{j,1} p^{j\nu}}{1 + \sum_{i=1}^{N_1} a_{i,1} p^{i\nu}} \quad (7)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. MISO model structure</figcaption></figure>

**Definition 2.2 (Local S-commensurability, [23]):** A MISO system, for a given model structure (5), has local S-commensurate orders $\nu_k$ if each subsystem $G_k$ has its own S-commensurate order (see Definition 2.1)

$$ \nu_k = \alpha_{1,k} \quad \text{for} \quad k = 1, \dots, K. $$

Therefore, each transfer function can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu_k}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu_k}}. \qquad (8) $$

**Definition 2.3 (Global S-commensurability, [23]):** A MISO system, for a given model structure (5), has a global S-commensurate order $\nu$ if all its subsystems have the same S-commensurate order according to Definition 2.2:

$$ \nu = \alpha_{1,1} = \alpha_{1,2} = \dots = \alpha_{1,K}, $$

As a consequence, all transfer functions can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu}}. \qquad (9) $$

### III. OUTPUT ERROR METHOD FOR SYSTEM IDENTIFICATION OF FRACTIONAL MISO MODELS

#### A. Problem formulation

Thanks to $H$ samples of input/output data $\{u_1(t_h), \dots, u_K(t_h), y^*(t_h)\}_{h=1}^H$, between $t_1 = 0$ and $t_H = t_{final}$, the aim is to estimate the MISO model parameters, where the MISO model $S$ is described by equation (4).

The parameter vector, $\theta$, is then defined as

$$ \theta = \begin{bmatrix} \rho \\ \mu \end{bmatrix} \qquad (10) $$

* where $\rho$ has all the coefficients of the MISO transfer function

$$ \rho = [\rho_1^T, \dots, \rho_K^T]^T \qquad (11) $$

with

$$ \rho_k^T = [b_{0,k}, b_{1,k}, \dots, b_{M_k,k}, a_{1,k}, \dots, a_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\rho$ has a total number of coefficients of $\sum_{k=1}^K (N_k + M_k + 1)$

* and where $\mu$ has all the MISO transfer function differentiation orders, that can be defined according to these cases,

Case 1: if the fractional MISO $S$ (4) has a S-commensurate order that is global, then the differentiation order vector is reduced to one parameter

$$ \mu = \nu; \qquad (12) $$

Case 2: if the fractional MISO system $S$ (4) has S-commensurate orders that are local, then the differentiation order vector is extended to $K$ parameters

$$ \mu = [\nu_1, \dots, \nu_K]^T; \qquad (13) $$
Case 3: if the fractional MISO system S (4) is a non commensurate MISO system, then $\mu$ gathers all the differentiation orders

$$ \mu = [\mu_1^T, \dots, \mu_K^T]^T \qquad (14) $$

with

$$ \mu_k^T = [\beta_{0,k}, \beta_{1,k}, \dots, \beta_{M_k,k}, \alpha_{1,k}, \dots, \alpha_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\mu$ has a total number of differentiation orders of $\sum_{k=1}^K (N_k + M_k + 1)$.

Note that fractional order systems should be correctly simulated with negligible simulation errors so that the parameter estimation is consistent.

## B. MISO-oe system identification

The true model structure being assumed known, the parameter vector $\theta$ is estimated using a nonlinear optimization method, based on Levenberg-Marquardt algorithm. The whole developments and convergence properties can be found in [14].

The MISO-oe algorithm is summarized in algorithm 1 where the coefficient vector $\rho$ is estimated along with the differentiation order vector $\mu$.

## Summary of MISO-oe algorithm for all parameter estimation

### Step 1: Initialization

iter = 0

Initialize $\hat{\mu}^0$ and compute $\hat{\rho}^0$ with the srivcf method.

Initialize the parameter vector $\theta^0 = [\rho^0, \mu^0]^T$

Compute the cost function $J(\theta^0)$ from (17)

### Step 2: Iterative Levenberg-Marquardt algorithm

#### do

i. $iter = iter + 1$

ii. Parameter estimation

Compute the parameter vector $\hat{\theta}^{iter}$ as

$$ \hat{\theta}^{iter+1} = \hat{\theta}^{iter} - \left\{ [\mathcal{H} + \xi I]^{-1} \frac{\partial J}{\partial \hat{\theta}} \right\|_{\hat{\theta}^{iter}}, \qquad (15) $$

with

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial J}{\partial \hat{\theta}} = \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \varepsilon(t_h) \quad : \text{gradient} \\[1em] \displaystyle \mathcal{H} \approx \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \frac{\partial \varepsilon(t_h)}{\partial \hat{\theta}} \quad : \text{pseudo-Hessian} \\[1em] \displaystyle \xi \quad : \text{Marquardt parameter} \end{array} \right. \qquad (16) $$

$$ \textbf{while} \sum_{\ell=1}^{\dim \hat{\theta}} \left| \frac{\hat{\theta}_{\ell}^{\text{iter}} - \hat{\theta}_{\ell}^{\text{iter}-1}}{\hat{\theta}_{\ell}^{\text{iter}-1}} \right| > \epsilon $$

where $\hat{\theta}_{\ell}^{\text{iter}}$ corresponds to the $\ell$-th element of the vector $\hat{\theta}^{\text{iter}}$ at iteration iter.

### Step 3: Compute the parametric covariance matrix $P_{\theta}$:

$$ P_{\theta} = \hat{\sigma}^2 \mathcal{H}^{-1} $$

where $\hat{\sigma}^2$ denotes the empirical estimation variance of the noise defined by:

$$ y^*(t) - \sum_{k=1}^{K} \hat{y}_k(t) $$

where $\hat{y}_k(t)$ is computed after convergence and the approximated Hessian $\mathcal{H}$ (see [23, Section 3]).

The error criterion is calculated by:

$$ J(\hat{\theta}) = \frac{1}{2} \| \varepsilon(t, \hat{\theta}) \|_2^2 = \frac{1}{2} E^T E, \qquad (17) $$

where

$$ E = [\varepsilon(t_1, \hat{\theta}), \dots, \varepsilon(t_H, \hat{\theta})]^T, $$
gathers the output error $\epsilon(t, \hat{\theta})$:

$$ \epsilon(t, \hat{\theta}) = y^*(t) - \hat{y}(t, \hat{\theta}) \qquad (18) $$

and the estimated output $\hat{y}(t, \hat{\theta})$ is defined as:

$$ \hat{y}(t) = \sum_{k=1}^{K} \hat{y}_k(t). \qquad (19) $$

In algorithm 1, note that $\epsilon$ is a scalar that is set for a desired precision of the estimates.

## IV. TERRESTRIAL CLIMATE SYSTEM MODELING

Such as proposed by [4], the MISO terrestrial climate model can be illustrated as in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Terrestrial climate model</figcaption></figure>

In the sciences of climatology, a *forcing factor* defines a disturbance in the energetic equilibrium of the earth. The global climate system reacts to these forcing factors through variations of the temperature. In our case, three forcing factors are considered, which correspond to the three input signals: the solar activity, the volcanic activity and the human activity. The output of the climate system is the global annual temperature.

Reconstructions of paleoclimates from proxies are combined with instrumental recordings (such as thermometer and plu- viometer recordings) to enlarge the comprehension of the climate variabilities. These past reconstructions of the climate and environment cover all time scale, variations from one year to another to those produced for million of years.

The input signals can be natural (the variations of radiation intensity (solar activity) and aerosol emissions (volcanic activity)) or anthropogenic through greenhouse gas emissions (carbon dioxide). The climate system reacts to these forcing factors with temperature variations.

### A. Solar activity (Total Solar Irradiance)

From year 1610, the sunspot measurements enabled reconstructing the total solar irradiance, and [18] discovered in 1843 a cycle of 11 years. Starting from 1978, these measurements are now carried out by satellites. Nowadays, the Total Solar Irradiance (TSI) is brought by TIM instrument (Total Irradiance Monitor), as Timv17 ([8]).

### B. Volcanic activity (Aerosol Optical Depth)

Aerosols can change the global climate temperature in several ways. The first effect is observed through the diffusion and absorption of solar radiation, which modifies the irradiative balance of the planet. The Aerosols Optical Depth (AOD) describes how aerosols affect the passage of light through the atmosphere, for a given wavelength (the wavelength is 550 nm). The signal of volcanic activity is considered as the connection of two recordings: satellite recording [19] with calibrated sulphate records in Antarctic ice cores [3].

### C. Human activity (CO2 concentration)

The atmospheric concentration of carbon dioxide, CO₂, is probably the most well-known greenhouse gas. This signal is mostly known to represent the human activity, and is constituted of two records: measurement of the centre of Mauna Loa since 1958 and proxy measurements.

### D. Temperature

Since 1850, the instrumental recordings sufficiently reliable of the temperature have been recorded on a surface covering the whole globe. The recent developments in the observational analyses of the temperature on the surface of earth soils (CRUTEM4) and seas (HadSST3) have been combined as the temperature reconstruction of HadCRUT4: such data have been gathered as a set of data on the evolution of the temperature in the world and in regions since 1850.
## E. Identification period

All used data start from 1850 as it is the recording period that is the most reliable. Moreover, the identification data are available until 2018 as the TSI reconstruction stops at this date. Consequently, [1850, 2018] will be the period of reference such as the IPCC. The Input/Output data are plotted on Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Input/output data set for the identification period [1850, 2018]</figcaption></figure>

## V. APPLICATION

All simulation have been done thanks to the CRONE toolbox<fnref n="1" />, which is an oriented-object toolbox developed for fractional order system modeling and also for system identification with fractional order models. The reader is invited see [9] to have a tutorial of this CRONE toolbox.

Remember that the identification period is set to [1850, 2018], and all means and trends have been removed.

### A. Identification of a MISO first order model

First, it is proposed to identify first order MISO models in order to only estimate the rational model coefficients. For the initialization, the coefficient estimated such first order MISO models in [4] are used. Note that in the models proposed in [4], all the models have the same denominator, which imposes the same dynamics to all forcing inputs. We propose here not to compare the estimated models with the ones proposed in [4], as the author, in the latter, has used different kind of filtering, such as Kalman filtering, to connect data from different kind of proxies.

In this context, let us apply the MISO-oe Alg. 1, for identifying rational models. Again, remember that only the coefficients are estimated here. The estimated models are given in Table I (col. 2) and the outputs are illustrated in Fig. 4. The cost function is computed to determine the fit goodness, specified as the normalized root mean squared error (NRMSE). It can be seen that the Fit cost function is around 45%. Moreover, the volcanic contribution is not well representative by using such rational models. The human activity contribution ($CO_2$) is much more significative compared to the two other input contributions.

### B. Identification of fractional MISO models with a global S-commensurate order

Let us apply the MISO-oe Alg. 1, for identifying a MISO global S-commensurate model. As thermal diffusion has a theoretical differentiation order of $v = 0.5$, it is proposed to identify global S-commensurate MISO-transfer functions where the global S-commensurate order is set to $v = 0.5$ and therefore, only the coefficients are estimated. The estimated models are given in Table I (col. 3) and the outputs are illustrated in Fig. 5. It can be seen that the Fit cost function has significantly increased and is around 56%. Moreover, the volcanic contribution is more representative by using such fractional order models. Note also, that the solar activity contribution catches a global temperature tendency. The human activity contribution ($CO_2$) is also significative compared to the two other input contributions.

<footnote n="1">¹The latest version of the object-oriented CRONE toolbox, usable on Matlab® (version 2020b and further). It can be downloaded at http://archive.ims-bordeaux.fr/CRONE/toolbox/</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. System identification with rational models; Left: global measured (—) and estimated (—) temperature outputs; Right: contributions from each input (TSI, AOD and CO₂ respectively from top to bottom)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. System identification with a global S-commensurate order model; Left: global measured (—) and estimated (—) temperature outputs; Right: contributions from each input (TSI, AOD and CO₂ respectively from top to bottom)</figcaption></figure>

### C. Identification of fractional MISO models with local S-commensurate orders

Applying the MISO-oe Alg. 1, for identifying MISO local S-commensurate models. As each input of the system are completely independent forcings, it is proposed to identify local S-commensurate MISO transfer functions where both the local S-commensurate orders and the coefficients are estimated. The estimated models are given in Table I (col. 4) and the outputs are illustrated in Fig. 6. It can be seen that the Fit cost function has even more increased and is around 67%. The volcanic contribution and the solar activity better catch the global temperature tendency. The human activity contribution (CO₂) is also significative and explains why there is an increase of the global temperature, but it is associated along with both solar and volcanic contributions.

<table><caption>TABLE I
COMPARISON OF THE ESTIMATED MODELS FOR A RATIONAL MODEL, A GLOBAL S-COMMENSURATE MODEL AND LOCAL S-COMMENSURATE MODELS BY USING THE *MISO-oe* ALG. 1</caption><thead><tr><th></th><th>Initial model<br>from [4]</th><th>Rational models</th><th>Global S-commensurate<br>model (ν = 0.5)</th><th>Local S-commensurate<br>model</th></tr></thead><tbody><tr><th>TSI</th><td><sup>4.25</sup>/<sub>86s+1</sub></td><td>G<sub>1</sub>(s) = <sup>0.36</sup>/<sub>7.2s+1</sub></td><td>G<sub>1</sub>(s) = <sup>1.07</sup>/<sub>50.52s<sup>0.5</sup>+1</sub></td><td>G<sub>1</sub>(s) = <sup>0.27</sup>/<sub>7.54s<sup>1.25</sup>+1</sub></td></tr><tr><th>EOA</th><td><sup>3.86</sup>/<sub>86s+1</sub></td><td>G<sub>2</sub>(s) = <sup>3.6</sup>/<sub>86.6s+1</sub></td><td>G<sub>2</sub>(s) = <sup>12.51</sup>/<sub>20.73s<sup>0.5</sup>+1</sub></td><td>G<sub>2</sub>(s) = <sup>2.69</sup>/<sub>80s<sup>1.84</sup>+1</sub></td></tr><tr><th>CO<sub>2</sub></th><td><sup>3.04</sup>/<sub>86s+1</sub></td><td>G<sub>3</sub>(s) = <sup>1.9</sup>/<sub>0.016s+1</sub></td><td>G<sub>3</sub>(s) = <sup>2.34</sup>/<sub>0.01s<sup>0.5</sup>+1</sub></td><td>G<sub>3</sub>(s) = <sup>2.03</sup>/<sub>0.61s<sup>1.82</sup>+1</sub></td></tr><tr><th>% Fit</th><td>32%</td><td>45%</td><td>56%</td><td>67%</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 6. System identification with local S-commensurate order models; Left: global measured (—) and estimated (—) temperature outputs; Right: contributions from each input (TSI, AOD and CO₂ respectively from top to bottom)</figcaption></figure>

## VI. CONCLUSION

This paper has enabled modeling the terrestrial global temperature by using simple MISO fractional models which bring parameter compactness and are well suited to model diffusive phenomena such as thermal systems. Also, the global climate change has been modeled thanks to system identification with fractional models by using the MISO version of the oe algorithm. For future works, an analysis with model order selection and predictions could be carried out. Another perspective would be to identify the climate change on a larger period of time (before 1850) by using proxy reconstructions of the inputs and having other temperatures reconstructions such as Ljungqvist [1], Moberg [15], Loehle [2] or even Mann&Jones [11] ones. Also, another method could be compared by using the optimal instrumental variable method such as the MISO-oosrivcf one, proposed in [23].

## REFERENCES

[1] F. C.Ljungqvist. A new reconstruction of temperature variability in the extra-tropical northern hemisphere during the last two millennia. *Geografiska Annaler: Series A, Physical Geography*, 92(3):339-351, 2010.

[2] C.Loehle. A 2000-year global temperature reconstruction based on non-treering proxies. *Energy & Environment*, 18(7):1049-1058, 2007.

[3] T. Crowley and M. Unterman. Technical details concerning development of a 1200 yr proxy index for global volcanism. *Earth System Science Data*, 5(1):187-197, 2013.

[4] P. de Larminat. Earth climate identification vs. anthropic global warming attribution. *Annual Reviews in Control*, 42:114-125, 2016.

[5] T. Djamah, R. Mansouri, S. Djennoune, and M. Bettayeb. Optimal low order model identification of fractional dynamic systems. *Applied Mathematics and Computation*, 206(2):543-554, 2008.

[6] Z. Gao, X. Lin, and Y. Zheng. System identification with measurement noise compensation based on polynomial modulating function for fractional-order systems with a known time-delay. *ISA transactions*, 79:62-72, 2018.

[7] L. Jing and G. Boling. Parameter identification in fractional differential equations. *Acta Mathematica Scientia*, 33(3):855-864, 2013.

[8] G. Kopp, G. Lawrence, and G. Rottman. The total irradiance monitor (tim): Science results. *Solar Physics*, 230(1):129-139, 2005.

[9] R. Malti and S. Victor. Crone toolbox for system identification using fractional order models. In *17th IFAC Symposium on System Identification (SYSID 2015)*, volume 48, pages 769–774, Beijing, China, 10 2015.

[10] R. Malti, S. Victor, and A. Oustaloup. Advances in system identification using fractional models. *Journal of Computational and Nonlinear Dynamics*, 3:021401.1-021401.7, 2008.

[11] M. Mann and P. Jones. Global surface temperatures over the past two millennia. *Geophysical Research Letters*, 30(15), 2023/01/30 2003.

[12] R. Mansouri, M. Bettayeb, and S. Djennoune. Multivariable fractional system approximation with initial conditions using integral state space representation. *Computers & Mathematics with Applications*, 59(5):1842–1851, 2010.

[13] D. Matignon. Stability properties for generalized fractional differential systems. *ESAIM proceedings - Systèmes Différentiels Fractionnaires - Modèles, Méthodes et Applications*, 5, 1998.

[14] A. Mayoufi, S. Victor, R. Malti, M. Chetoui, and M. Aoun. Output error MISO system identification using fractional models. *Fractional Calculus & Applied Analysis*, 5(24):1601-1618, 2021.

[15] A. Moberg, D. Sonechkin, K. Holmgren, N. M. Datsenko, and W. Karlén. Highly variable northern hemisphere temperatures reconstructed from low- and high-resolution proxy data. *Nature*, 433(7026):613-617, 2005.

[16] A. Narang, S. Shah, and T. Chen. Continuous-time model identification of fractional-order models with time delays. *IET Control Theory & Applications*, 5(7):900-912, 2011.

[17] I. Podlubny. *Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications*. Academic Press, San Diego, 1999.

[18] H. Schwabe and H. Schwabe Herrn. Sonnen—beobachtungen im jahre 1843. *Astronomische Nachrichten*, 21(15):234-235, 1844.

[19] S. Solomon, J. Daniel, R. Neely, J. Vernier, E. Dutton, and L. Thomason. The persistently variable “backgroun” stratospheric aerosol layer and global climate change. *Science*, 333(6044):866-870, 2011.

[20] Y. Tang, H. Liu, W. Wang, Q. Lian, and X. Guan. Parameter identification of fractional order systems using block pulse functions. *Signal Processing*, 107:272–281, 2015.

[21] M. Thomassin and R. Malti. Subspace method for continuous-time fractional system identification. *IFAC Proceedings Volumes - 15th IFAC Symposium on System Identification (SYSID'2009)*, 42(10):880-885, 2009.

[22] S. Victor, J. Duhe, P. Melchior, Y. Abdelmoumen, and F. Roubertie. Long memory recursive prediction error method for identification of continuous-time fractional models. *Nonlinear Dynamics*, 110:635-648, 06 2022.
[23] S. Victor, A. Mayoufi, R. Malti, M. Chetoui, and M. Aoun. System identification of MISO fractional systems: parameter and differentiation order estimation. *Automatica*, 141, 07 2022.