Claudiu Tiberiu Albulescu
email: claudiu.albulescu@upt.ro.

Eugenia Grecu
email: eugenia.grecu@upt.ro.

Claudiu Tiberiu

Tiberiu § Claudiu

Albulescu

Government Interventions and Sovereign Bond Market Volatility during COVID 19: A Quantile Analysis

Keywords: COVID-19, government bond price volatility, government policy responses, international financial markets, containment and closure, economic support, panel quantile regression. JEL classifications: G01, G12, G15, G18, G41, H12, I18

Introduction

A careful mapping of the COVID-19 research shows that much of its efforts and attention has, so far, been focused on the possible impact of government interventions on the equity market. 1 However, the research on sovereign bonds, which account for a large part of global trading volume, has remained relatively sparse. To be specific, only a handful of studies focus on the impact of the pandemic on bond yields, prices, liquidity, or term spreads (see, for example [START_REF] Arellano | Deadly debt crises: COVID-19 in emerging markets[END_REF][START_REF] Gubareva | The impact of Covid-19 on liquidity of emerging market bonds[END_REF][START_REF] He | Treasury inconvenience yields during the COVID-19 crisis[END_REF][START_REF] O'hara | Anatomy of a liquidity crisis: corporate bonds in the COVID-19 crisis[END_REF][START_REF] Sène | Overshooting of sovereign emerging Eurobond yields in the context of COVID-19[END_REF][START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF]Zaremba et al., 2022, respectively). Hence, the primary goal of this study is to improve the understanding of the COVID-19-bond market nexus. Specifically, we scrutinize the effect of the government policy responses to the pandemic and of the COVID-19-induced uncertainty on sovereign bond volatility, showing that this effect is asymmetric, being influenced by the volatility level. Our is the first paper assessing the asymmetric, nonlinear effect of government interventions on sovereign bond market Volatility during COVID 19

During the COVID-19 outbreak, financial markets have experienced extraordinary levels of uncertainty leading to significant price drawdowns, volatility spikes, and liquidity shortages [START_REF] Baker | The unprecedented stock market reaction to COVID-19[END_REF][START_REF] Belaid | Emerging and advanced economies markets behaviour during the COVID-19 crisis era[END_REF][START_REF] Fakhfekh | Volatility dynamics of the Tunisian stock market before and during the COVID-19 outbreak: Evidence from the GARCH family models[END_REF]Fetzer et al., 2020;[START_REF] Lee | The impact of COVID-19 pandemic on hospitality stock returns in China[END_REF][START_REF] Lyócsa | Fear of the coronavirus and the stock markets[END_REF][START_REF] Szczygielski | The only certainty is uncertainty: An analysis of the impact of COVID-19 uncertainty on regional stock markets[END_REF]. Importantly, besides the pandemic itself, which generated a specific form of uncertainty associated with the increased number of new infection cases and deaths [START_REF] Albulescu | Oil-US Stock Market Nexus: Some insights about the New Coronavirus Crisis[END_REF] or with news related to COVID-19 [START_REF] Ftiti | Does non-fundamental news related to COVID-19 matter for stock returns? Evidence from Shanghai stock market[END_REF], global economies have faced unprecedented government policy responses. These interventions may significantly affect financial market volatility; however, the direction of these forces is far from trivial. On the one hand, any government action may induce additional uncertainty (Pastor and Varonesi, 2012), which in turn, leads to an increase in volatility of government bond markets. On the other hand, several other papers consider government interventions as responsible actions that may curb down the adverse effects of crises and uncertainty [START_REF] Amengual | Resolution of policy uncertainty and sudden declines in volatility[END_REF][START_REF] Kizys | From COVID-19 herd immunity to investor herding in international stock markets: The role of government and regulatory restrictions[END_REF]-which can also be the case for sovereign debt.

Furthermore, the interventions may take different forms. Some of them include containment and closure policies that are targeted at curbing the spread of the pandemic; others 1 See, e.g., [START_REF] Albulescu | COVID-19 and the United States financial markets' volatility[END_REF], [START_REF] Alexakis | COVID-19 containment measures and stock market returns: An international spatial econometrics investigation[END_REF], [START_REF] Baig | Deaths, panic, lockdowns and U.S. equity markets: The case of COVID-19 pandemic[END_REF], [START_REF] Duan | COVID-19 Sentiment and the Chinese Stock Market: Evidence from the Official News Media and Sina Weibo[END_REF], [START_REF] Gao | To what extent does COVID-19 drive stock market volatility? A comparison between the U.S. and China[END_REF], [START_REF] Goodell | Did Congress trade ahead? Considering the reaction of U.S. industries to COVID-19[END_REF], [START_REF] James | Association between COVID-19 cases and international equity indices[END_REF], [START_REF] Ozkan | Impact of COVID-19 on stock market efficiency: Evidence from developed countries[END_REF], [START_REF] Seven | World equity markets and COVID-19: Immediate response and recovery prospects[END_REF], [START_REF] Szczygielski | The only certainty is uncertainty: An analysis of the impact of COVID-19 uncertainty on regional stock markets[END_REF], [START_REF] Zaremba | Infected markets: Novel coronavirus, government interventions, and stock return volatility around the globe[END_REF], [START_REF] Zaremba | COVID-19, government policy responses, and stock market liquidity around the world: A note[END_REF], [START_REF] Zhang | Financial markets under the global pandemic of COVID-19[END_REF].

provide economic support to both enterprises and consumers. The impact of these very different actions does not need to be identical [START_REF] Kizys | From COVID-19 herd immunity to investor herding in international stock markets: The role of government and regulatory restrictions[END_REF]. If we consider for instance the containment and closure policies, we expect an immediate negative impact of these measures on the real economy. However, these policies limit COVID-19 propagation and might restore the investors confidence. In this case, the containment policies might reduce the bond price volatility. This is also the case of economic support measures which in the short run generates a positive market sentiment but in the long run they might be associated with fiscal imbalances, increasing thus the market uncertainty. Consequently, we attempt to shed light on this issue, and explore the impact of different government policy responses on government bond volatility.

Towards this end, we examine the behaviour of sovereign bonds in 31 countries during the recent pandemic. Contrary to earlier studies (e.g. [START_REF] Kizys | From COVID-19 herd immunity to investor herding in international stock markets: The role of government and regulatory restrictions[END_REF][START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF], we employ the Canay's (2011) panel-quantiles regression approach with fixed effects to determine whether the relationship is consistent across several parts of the bond's volatility distribution. Otherwise said, we investigate if the impact of interventions on sovereign bond market volatility is different, depending on the volatility level.

We therefore build upon [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF] and we extend their analysis in three ways. First, we posit that the effect of government interventions on bond market volatility is not linear and is influenced by the level of volatility recorded in each market. More precisely, it is well known that countries with more developed financial markets tend to record a reduced volatility level [START_REF] Wang | Financial development and long-run volatility trends[END_REF]. These mature markets are not reacting to news and uncertainty in the same way the emerging financial markets does. Therefore, we expect that the governmental interventions will have a stronger impact on bond price volatility at upper quantiles, that is, for more volatile bond markets. Highlighting the asymmetric effect of governmental interventions on bond market volatility represents main advantage of a panel quantile approach over the classic panel data models. In addition, a quantile approach has other advantages, including its robustness to non-normality, as well as to heteroscedasticity, skewness, and leptokurtosis-all of which are typical financial data features (Canay, 2011).

The estimated conditional quantile functions provide a much more complete image of the covariates' effect on the location, scale, and shape of the distribution of a response variable [START_REF] Rosen | Set identification via quantile restrictions in short panels[END_REF]. Application of this method to study the relationship between COVID-19 and sovereign bond volatility is uncovered by the extant literature. We demonstrate that both the spread of the infections and the policy measures augment the bond market volatility. As novelty, we show that the impact of government interventions increases for upper quantiles, that is, for more volatile markets. The effect is driven principally by containment and closure policies, such as lockdowns or school closing. On the other hand, economic support policies tend to stabilize bond price fluctuations.

Second, we cover the first two waves of the pandemics, while Zaremba et al.'s (2021a) data span only covers the first wave. We investigate the two waves of pandemics (for a description, please refer to [START_REF] Duttilo | Volatility Modeling: An Overview of Equity Markets in the Euro Area during COVID-19 Pandemic[END_REF], given the high level of uncertainty and volatility recorded in 2020. Starting with 2021, financial markets volatility decreased, pointing in the favour of shock accommodation and uncertainty downturn. Moreover, the bond purchases measures (see for example the Federal Reserve quantitative easing program), diminished the market volatility. Third, we check for the "Monday effect" of new infection cases. The new infections are reported at the date "t", for the tests done at date "t-1" [START_REF] Albulescu | COVID-19 and the United States financial markets' volatility[END_REF]. Given that less tests are performed during the week-end, the number of new infections cases is smaller Monday compared with the other days of the week.

Our findings contribute to the literature on the effect of the COVID-19 outbreak on the bond markets volatility (e.g., [START_REF] He | Treasury inconvenience yields during the COVID-19 crisis[END_REF][START_REF] Arellano | Deadly debt crises: COVID-19 in emerging markets[END_REF][START_REF] Sène | Overshooting of sovereign emerging Eurobond yields in the context of COVID-19[END_REF][START_REF] Gubareva | The impact of Covid-19 on liquidity of emerging market bonds[END_REF]) in several ways. 2 In particular, our study is most closely related to [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF], who applied simple panel regressions to delve into the pandemic-bond volatility nexus.

Significantly, our conclusions expand the findings of that study, showing that the interventions impact is influenced by the bond's volatility distribution. Whereas [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF] only found the link between bond volatility and the economic support policies, we also document the essential role of containment and closure interventions-which amplify the price variability.

Consequently, while [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF] find the overall stabilizing effect of the government, we demonstrate their detrimental impact.

The structure of the paper is as follows: Section 2 describes data. Section 3 outlines the methods. Section 4 discusses the empirical findings and Section 5 the robustness checks.

Finally, Section 6 concludes.

Data

As in [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF], the data consists of information on different policy responses from 31 countries that are covered by Datastream: Australia, Austria, Belgium, Canada, China, Czechia, Denmark, Finland, France, Germany, Greece, Hungary, India, Indonesia, Ireland, Italy, Japan, Korea, Mexico, the Netherlands, New Zealand, Norway, Poland, Portugal, Singapore, South Africa, Spain, Sweden, Switzerland, the United Kingdom, and the United States. All of the bond-related data and variables are derived from the Datastream 10-Year Government Bond Total Returns indices. The 10-year maturities are the primary choice in asset pricing literature due to high liquidity and broad international coverage (e.g., [START_REF] Andres | Measuring abnormal credit default swap spreads[END_REF]Baltussen et al., 2021). The sample period encompasses the spread of the pandemic, running from January 1, 2020, through November 3, 2020, covering thus the first two waves of COVID-19. 3 Following the typical approach in international bond pricing studies (e.g., [START_REF] Asness | Value and momentum everywhere[END_REF], we express the market data in U.S. dollars and the risk-free rate is proxied by the U.S. one-month Treasury-bill rate from Kenneth R. French's website. 4To quantify day-to-day changes in volatility, we build on [START_REF] Antonakakis | Dynamic spillovers between commodity and currency markets[END_REF]; [START_REF] Khalifa | Return distributions and volatility forecasting in metal futures markets: Evidence from gold, silver, and copper[END_REF]andZaremba et al. (2020, 2021a-all of whom employ absolute measures of daily returns. Furthermore, to extract the country-specific volatility component, free of the impact of systematic risks, we replace the raw returns with residuals from a factor model. To be precise, in order to capture the multidimensionality of bond returns, we utilize the comprehensive seven-factor model originating from [START_REF] Zaremba | Volatility in international sovereign bond markets: The role of government policy responses to the COVID-19 pandemic[END_REF] risk factors, respectively. The detailed description of factor construction is provided in Table A1 in the Appendix.

We derive look-ahead bias-free absolute daily residuals by performing the following steps.

To begin, for each day t we run the regression (1) using five years of trailing data ending on day t-1. Subsequently, we utilize the coefficient estimates and factor realizations from day t to calculate the expected daily returns. Finally, we compute the residual returns as the difference between the actual return realizations on day t and their expected values implied by the model (1).

Our main explanatory variables are based on the policy response indices from the Oxford COVID-19 Government Response Tracker. 5 The indices aggregate data on different government interventions following the COVID-19 outbreak, such as canceling public gatherings and closing workplaces, social distancing requirements, debt relief, or income support. In our baseline approach, we use three different indices: the Government Response Index (gvt), which incorporates information on all types of policies; as well as the Containment and Health Index (cntm) and the Economic Support Index (eco). The latter two constitute subindices of gvt and reflect different types of policies. Whereas cntm concentrates on containment and closure policies aimed at curbing the pandemic, eco is about economic support to consumers and enterprises during the pandemic.

Besides the primary independent variable, we include a range of additional control variables. These include bond duration (dur), default risk (cred), money market rate (mmr) and convexity (cx), carry (car), momentum (mom), reversal (rev), and "Monday effect" dummies (dummy). The detailed descriptions for all variables are presented in Table A2 in the Appendix.

Table 1 reports the descriptive statistics for the key variables. Though not reported here, all variables are stationary according to [START_REF] Maddala | A comparative study of unit root tests with panel data and a new simple test[END_REF], and [START_REF] Pesaran | A simple panel unit root test in the presence of cross-section dependence[END_REF] unit root tests.

Methods

Quantile regression models are useful to account for unobserved heterogeneity and asymmetry. In addition, when relying on fixed-effects models, researchers can control for unobserved covariates. A combination of these approaches represents the basis of panel quantile fixed effects models that is proposed in the literature (e.g., [START_REF] Koenker | Quantile regression for longitudinal data[END_REF][START_REF] Lamarche | Robust penalized quantile regression estimation for panel data[END_REF][START_REF] Galvao | Quantile regression for dynamic panel data with fixed effects[END_REF][START_REF] Rosen | Set identification via quantile restrictions in short panels[END_REF].

Let us consider the following model:

𝑌 𝑖𝑡 = 𝑋 𝑖𝑡 ′ 𝜃(𝑈 𝑖𝑡 ) + 𝛼 𝑖 , (2)
where: 𝑡 = 1, . . . , 𝑇; 𝑖 = 1, . . . , 𝑛; 𝑌 𝑖𝑡 and 𝑋 𝑖𝑡 represent the observable variables; 𝑈 𝑖𝑡 is an unobservable component; 𝑋 𝑖𝑡 ′ includes a constant term; and 𝜃(𝜏) is the parameter of interest.

It is assumed that the function 𝜏 → 𝑋 ′ 𝜃(𝜏) is increasing in 𝜏 ∈ (0,1). In the case 𝛼 𝑖 is observable, it follows that:

𝑃[𝑌 𝑖𝑡 ≤ 𝑋 𝑖𝑡 ′ 𝜃(𝑈 𝑖𝑡 ) + 𝛼 𝑖 |𝑋 𝑖 , 𝛼 𝑖 ] = 𝜏, (3)
where:

𝑈 𝑖𝑡 ~𝑈[0,1], conditional on 𝑋 𝑖 = (𝑋 𝑖1 ′ , . . . , 𝑋 𝑖𝑇 ′ )′ and 𝛼 𝑖 .
The challenge is the 𝜃(𝜏) identification, which cannot be done by imposing only covariates quantile restrictions [START_REF] Rosen | Set identification via quantile restrictions in short panels[END_REF]

where:

𝔼 𝑛𝑇 (•) ≡ (𝑛𝑇) -1 ∑ ∑ (•) 𝑛 𝑖=1 𝑇 𝑡=1
.

Starting from this framework, we use the first lag of explanatory variables to avoid any endogeneity bias, and we test the following general regression:

𝑅 𝑖𝑡 = 𝛼 0 + 𝛼 1 𝑋 𝑖𝑡-1 + 𝛼 2 𝑍 𝑖𝑡-1 + 𝜇 𝑖 + 𝛾 𝑡 + 𝜀 𝑖𝑡 , (7)
where: Rit is the daily measure of sovereign bond volatility in the country i on day t, i.e., the absolute residuals from the model (1); 𝛼 0 represents a constant term; 𝑋 𝑖𝑡-1 is the vector of COVID-19 variables, represented by new cases of infection and governmental response to the SARS-Cov-2 pandemic; 𝑍 𝑖𝑡-1 is the vector of control variables defined Section 2; and 𝜇 𝑖 are the time-invariant country-specific effects, 𝛾 𝑡 are the time-specific effects and 𝜀 𝑖𝑡 are the error terms. Besides the impact of policy responses to the pandemic, our baseline regression analysis uncovers the role of the pandemic itself: growth in the number of new infections translates into an increase in the bond market volatility. This observation matches similar earlier findings from equity markets (e.g., [START_REF] Zaremba | Infected markets: Novel coronavirus, government interventions, and stock return volatility around the globe[END_REF][START_REF] Baig | Deaths, panic, lockdowns and U.S. equity markets: The case of COVID-19 pandemic[END_REF][START_REF] Albulescu | COVID-19 and the United States financial markets' volatility[END_REF], arguing that COVID-19 induced uncertainty contributes to the instability of stock prices. Interestingly, the COVID-19 figures more strongly influence the sovereign bond prices located at the lower and medium volatility quantiles when compared with high volatility bonds. Consequently, less volatile financial marketstypically found in developed countriesare more sensitive to the changes in COVID-19 figures.

Empirical findings

The overall government response index, as examined in Table 2, encompasses various interventions that may exhibit differing economic impact. Therefore, in the subsequent analysis, we distinguish between containment and closure measures vs. economic support policies. These two categories are measures with cntm and eco, respectively.

Table 3 presents the influence of the containment and closure measures on bond market volatility. This additional analysis unequivocally reveals the underlying source of the government policy responses on the market volatility. Highly significant cntm coefficients indicate that these containment and closure interventions constitute a major contributor to bond price variability. In accordance with previous results, the policy measures generate a more substantial impact when we test the higher quantiles of the distribution, whereas the spread of the pandemic as measured by the inf variable is more powerful in low quantiles. Finally, similar to the previous case, the effect of policy measures (nor the spread of the disease) on bond price volatility is not significant for very volatile markets (i.e., the 0.95 quantile).

Let us now turn to the role of the other category of government interventions: economic support policies (eco). Table 4 demonstrates the results of another set of quantile regressions to capture the role of this category of government actions. Our analysis uncovers a negative impact on bond price volatility for lower and upper quantiles, but not for middle quantiles (Table 4).

For sovereign bonds with smaller and higher volatility, economic support interventions stabilize the markets. The effect is more substantial for the upper quantiles, which is in line with the impact generated by other policy interventions (see, for example, the results reported in Table 3). To sum up our considerations, we find that both containment and closure restrictionsas well as economic support policies -affect the market volatility; however, the direction of the impact is opposite. Whereas the first category tends to boost the market fluctuations, the latter helps to stabilize the market.

Robustness checks

To assure the validity of our findings, we run a number of additional robustness checks.

First, use a different metric to compute the sovereign bond price volatility, relying on absolute raw returns rather than on risk-adjusted returns (residuals). 6 These results are reported in Table 5 and are very similar with those reported in Section 4. Third, in an unreported analysis, we employ a modified measure of the strictness of government policies -namely the Stringency Indexwhich is also sourced from the Oxford COVID-19 Government Response Tracker. None of these extra robustness checks materially affect our findings. Our overall conclusions remain unaffected.

Conclusions

We examine the impact of governments' interventions and COVID-19 numbers on the volatility of sovereign bonds. We apply quantile regressions to a sample of 31 countries to scrutinize the importance of different types of policy responses. We posit that the impact of COVID-19 on sovereign bonds is influenced by the level of market volatility.

Our findings demonstrate that the direction of the effect on government bond returns' volatility depends strongly on the type of interventions. Confinement and closure restrictions increase market uncertainty and, in consequence, drive the return volatility up. In contrast, economic support measures tend to calm the volatility level in trading and enhance market stability. Further, we show that the impact of COVID-19 induced policy measures and related uncertainty is higher in the case of more volatile markets (i.e., at higher quantiles).

The conclusions from this study yield clear, practical implications. Since confinement and closure restrictions amplify volatility, our results imply that governments should be transparent and clear with their plans about this type of interventions in the short and longer terms. The COVID-19 period is overwhelmed with increased uncertainty of which governments interventions may even worsen it. Hence, providing information publicly as soon as possible may calm the adverse effect of closures. In addition, since economic interventions seem to be associated with positive response, it does not mean that transparency about economic steps is not needed, especially if such supportive actions are expected to increase the fiscal deficit.

The findings also imply that investors can exploit this information to better shape their investment decisions. They should be aware that non-economic interventions, which are not directly related to financial markets, may spillover to capital markets and are not limited to the equity markets. Therefore, investors, and particularly those operating in the fixed-income markets should monitor the changes in government policy and make the required adjustments to their portfolios. More precisely, to anticipate the volatility dynamics, they need to analyze the type of interventions and the sovereign bond markets' characteristics.

:

  𝑅 𝑖,𝑡 = 𝛼 𝑖 + 𝛽 𝑖

	𝑀𝐾𝑇 𝑀𝐾𝑇 𝑡 𝐹 + 𝛽 𝑖 𝐷𝑈𝑅 𝐷𝑈𝑅 𝑡 𝐹 + 𝛽 𝑖 𝐶𝑅𝐸𝐷 𝐶𝑅𝐸𝐷 𝑡 𝐹 + 𝛽 𝑖 𝑆𝐼𝑍𝐸 𝑆𝐼𝑍𝐸 𝑡 𝐹 +
	𝛽 𝑖 𝑀𝑂𝑀 𝑀𝑂𝑀 𝑡 𝐹 + 𝛽 𝑖 𝑅𝐸𝑉 𝑅𝐸𝑉 𝑡 𝐹 + 𝛽 𝑖 𝐶𝐴𝑅 𝐶𝐴𝑅 𝑡 𝐹 + 𝜀 𝑖,𝑡 .	(1)

Ri,t indicates the daily payoff on a country government i on day t, αi measures the abnormal return, and εi,t is the error term. The regression coefficients 𝛽 𝑖 𝑀𝐾𝑇 , 𝛽 𝑖 𝐷𝑈𝑅 , 𝛽 𝑖 𝐶𝑅𝐸𝐷 , 𝛽 𝑖 𝑆𝐼𝑍𝐸 , 𝛽 𝑖 𝑀𝑂𝑀 , 𝛽 𝑖 𝑅𝐸𝑉 , and 𝛽 𝑖 𝐶𝐴𝑅 reflect the exposures to the market risk (MKT F ), duration (DUR F ), credit risk (CRED F ), size (SIZE F ), momentum (MOM F ), long-term reversal (REV F ), and carry (CAR F )

Table 1 .

 1 Summary Statistics

		Mean	Std. Dev.	Min	Max
	R1	3.392	4.283	0.001	104.7
	R2	5.607	6.966	0.000	125.0
	gvt	49.32	24.61	0.000	95.54
	cntm	49.28	24.70	0.000	98.96
	stg	47.99	26.62	0.000	100.0
	eco	49.57	35.21	0.000	100.0
	inf	4.578	3.097	0.000	11.49
	dur	8.559	1.124	5.390	10.45
	cred	4.623	3.663	1.000	13.00
	mmr	0.788	1.805	-1.957	7.300
	car	0.693	1.107	-1.269	7.623
	cx	81.47	20.98	34.86	119.9
	size	16.09	0.902	14.13	18.38
	mom	-0.496	0.589	-3.138	3.220
	rev	-0.656	1.547	-16.54	3.680
	Notes: (i) R1 -daily absolute residuals from a seven-factor model, R2 -daily absolute returns in U.S. dollars, gvt
	-government response index, cntm -containment and health index, eco -economic support index, stg -original
	stringency index, inf -new infection cases, dur -duration, cred -credit rating, mmr -money market rate, car -
	yield-based carry, cx -convexity, size -bond market value, mom -momentum, and rev -reversal; (ii) 6,789
	observations. The sovereign bonds price volatility variables (R1 and R2) are adjusted (|ln(1+R)|) and multiplied by
	1,000 before running the regression.

  This way 𝛼 𝑖 is present in the conditional mean of 𝑌 𝑖𝑡 , allowing Canay (2011) to compute the two-step estimator 𝜃 ̂𝜇. First, we obtain a consistent estimator of 𝛼 𝑖 (√𝑇) and 𝜃𝜇 (√𝑛𝑇), with 𝛼 ̂𝑖 ≡ 𝐸 𝑇 [𝑌 𝑖𝑡 -𝑋 𝑖𝑡 ′ 𝜃 ̂𝜇]. Second, we define 𝑌 ̂𝑖 ≡ 𝑌 𝑖𝑡 -𝛼 ̂𝑖 and 𝜃 ̂𝜇 becomes:

	𝑌 𝑖𝑡 = 𝑋 𝑖𝑡 ′ 𝜃𝜇 + 𝛼 𝑖 + 𝑢 𝑖𝑡 ,	𝐸(𝑢 𝑖𝑡 |𝑋 𝑖 , 𝛼 𝑖 ) = 0.	(5)
	𝜃 ̂𝜇 ≡ argmin 𝜃∈Θ	𝔼 𝑛𝑇 [𝜌 𝜏 (𝑌 ̂𝑖𝑡 -𝑋 𝑖𝑡 ′ 𝜃 ̂𝜇],

. If 𝑄 𝑌 (𝜏|𝑋) is the 𝜏-quantile of a random variable Y conditional on X and 𝑒 𝑖𝑡 (𝜏) ≡ 𝑋 𝑖𝑡 ′ [𝜃(𝑈 𝑖𝑡 ) -𝜃(𝜏)], equation (2) can be written as:

𝑌 𝑖𝑡 = 𝑋 𝑖𝑡 ′ 𝜃(𝑈 𝑖𝑡 ) + 𝛼 𝑖 + 𝑒 𝑖𝑡 (𝜏),

(4)

Canay (2011) proves that 𝜃(𝜏) is identified for 𝑇 ≥ 2 under independence restrictions and the existence of moments. When we move from identification to estimation, we get rid of the fixed effects under the assumption that 𝛼 𝑖 is a location shift.

Practically, Canay (2011)

assumes that only 𝜃(𝜏) and 𝑒 𝑖𝑡 (𝜏) depend on 𝜏 and transforms equation (4) as follows:

Table 2

 2 reports the results of the quantile regressions that account for the overall role of the policy responses. The positive and highly significant coefficients on gvt suggest that government interventions amplify bond market volatility. The effect is robust across the majority of quantiles tested. The impact of policy measures increases when we shift from lower

to higher quantiles. In other words, a volatile financial market environment implies a stronger reaction to the COVID-19 induced policy measures. The only exception is the most volatile quantile-where the gvt does not differ significantly from zero.

Table 2 .

 2 Panel Conditional Quantile Regression -Government Response Index

			Lower Quantiles		Middle Quantiles			Upper Quantiles
		0.05	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
	gvt	0.007*** 0.006*** 0.005***	0.007*** 0.010*** 0.012*** 0.014***	0.013*** 0.011** -0.007
		(0.002) (0.002) (0.002)	(0.001) (0.002) (0.002) (0.002)	(0.003) (0.005) (0.012)
	inf	0.201*** 0.188*** 0.192***	0.181*** 0.163*** 0.159*** 0.136***	0.135*** 0.117*** 0.097
		(0.021) (0.017) (0.016)	(0.014) (0.016) (0.018) (0.022)	(0.025) (0.041) (0.100)
	dur	-2.444*** -2.857*** -2.373***	-2.218*** -2.088*** -2.034*** -1.890***	-1.529*** -1.482* -3.128
		(0.475) (0.374) (0.357)	(0.318) (0.361) (0.394) (0.495)	(0.559) (0.899) (2.190)
	cred	-0.379*** -0.382*** -0.329***	-0.327*** -0.326*** -0.326*** -0.321***	-0.328*** -0.343*** -0.443***
		(0.020) (0.015) (0.015)	(0.013) (0.015) (0.016) (0.020)	(0.023) (0.038) (0.092)
	mmr	2.361*** 2.518*** 2.535***	2.558*** 2.617*** 2.689*** 2.753***	2.869*** 2.946*** 3.495***
		(0.045) (0.036) (0.034)	(0.030) (0.034) (0.038) (0.047)	(0.053) (0.086) (0.211)
	car	2.250*** 2.422*** 2.338***	2.435*** 2.484*** 2.598*** 2.696***	2.879*** 3.220*** 4.721***
		(0.071) (0.056) (0.053)	(0.048) (0.054) (0.059) (0.074)	(0.084) (0.135) (0.329)
	cx	0.121*** 0.151*** 0.127***	0.120*** 0.114*** 0.114*** 0.105***	0.085*** 0.083	0.192
		(0.027) (0.021) (0.020)	(0.018) (0.021) (0.023) (0.028)	(0.032) (0.052) (0.127)
	size	-0.445*** -0.506*** -0.540***	-0.621*** -0.717*** -0.844*** -0.936***	-1.073*** -1.274*** -1.568***
		(0.054) (0.042) (0.040)	(0.036) (0.041) (0.044) (0.056)	(0.063) (0.102) (0.248)
	mom -0.534*** -0.689*** -0.616***	-0.673*** -0.773*** -0.859*** -0.930***	-1.102*** -1.562*** -2.664***
		(0.082) (0.065) (0.062)	(0.055) (0.062) (0.068) (0.086)	(0.097) (0.156) (0.381)
	rev	0.565*** 0.482*** 0.470***	0.491*** 0.494*** 0.528*** 0.543***	0.519*** 0.591*** 0.658***
		(0.036) (0.029) (0.027)	(0.024) (0.028) (0.030) (0.038)	(0.043) (0.069) (0.169)
	dummy 0.004	-0.062	-0.087	-0.118	-0.1537 -0.148* -0.156	-0.062	0.0663 -0.110
		(0.108) (0.085) (0.081)	(0.072) (0.082) (0.089) (0.112)	(0.127) (0.205) (0.498)

Notes:

The table reports slope coefficients from panel regressions along with the corresponding standard errors. (i) standard error in parentheses; (ii) *** p<0.01, ** p<0.05, * p<0.1; (iii) 6,788 observations; (iv) gvtgovernment response index, infnew infection cases, durduration, credcredit rating, mmrmoney market rate, caryield-based carry, cxconvexity, sizebond market value, mom -momentum, revreversal, and dummybinary variable that takes value 1 if Monday and 0 for the rest of the weekdays.

Table 3 .

 3 Panel Conditional Quantile Regression -Containment and Health Index

			Lower Quantiles		Middle Quantiles			Upper Quantiles
		0.05	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
	cntm	0.010*** 0.009*** 0.007***	0.009*** 0.013*** 0.015*** 0.017***	0.019*** 0.019*** 0.007
		(0.002) (0.002) (0.002)	(0.002) (0.002) (0.002) (0.002)	(0.003) (0.005) (0.012)
	inf	0.182*** 0.170*** 0.174***	0.166*** 0.143*** 0.135*** 0.116***	0.101*** 0.075*	0.001
		(0.021) (0.017) (0.016)	(0.014) (0.016) (0.017) (0.022)	(0.025) (0.041) (0.101)
	dur	-2.209*** -2.498*** -2.067***	-1.858*** -1.744*** -1.697*** -1.596***	-1.243** -1.115	-2.561
		(0.466) (0.375) (0.364)	(0.325) (0.364) (0.386) (0.480)	(0.562) (0.910) (2.221)
	cred	-0.374*** -0.374*** -0.317***	-0.324*** -0.314*** -0.322*** -0.315***	-0.315*** -0.329*** -0.424***
		(0.019) (0.015) (0.015)	(0.013) (0.015) (0.016) (0.020)	(0.023) (0.038) (0.093)
	mmr	2.293*** 2.470*** 2.475***	2.492*** 2.538*** 2.611*** 2.676***	2.791*** 2.886*** 3.419***
		(0.044) (0.036) (0.035)	(0.031) (0.035) (0.037) (0.046)	(0.054) (0.087) (0.213)
	car	2.170*** 2.314*** 2.249***	2.345*** 2.389*** 2.535*** 2.608***	2.782*** 3.070*** 4.610***
		(0.070) (0.056) (0.054)	(0.049) (0.054) (0.058) (0.072)	(0.084) (0.137) (0.334)
	cx	0.110*** 0.131*** 0.111***	0.099*** 0.095*** 0.095*** 0.089***	0.070** 0.062	0.164
		(0.027) (0.021) (0.021)	(0.019) (0.021) (0.022) (0.028)	(0.032) (0.053) (0.129)
	size	-0.453*** -0.516*** -0.553***	-0.625*** -0.716*** -0.836*** -0.936***	-1.084*** -1.302*** -1.535***
		(0.052) (0.042) (0.041)	(0.036) (0.041) (0.043) (0.054)	(0.063) (0.103) (0.251)
	mom -0.521*** -0.668*** -0.583***	-0.642*** -0.748*** -0.821*** -0.878***	-1.082*** -1.518*** -2.690***
		(0.080) (0.065) (0.063)	(0.056) (0.063) (0.067) (0.083)	(0.097) (0.158) (0.385)
	rev	0.588*** 0.493*** 0.481***	0.496*** 0.508*** 0.539*** 0.554***	0.531*** 0.607*** 0.698***
		(0.035) (0.028) (0.028)	(0.025) (0.028) (0.029) (0.036)	(0.043) (0.070) (0.170)
	dummy -0.012 -0.051	-0.080	-0.131* -0.148* -0.162* -0.161	-0.057	0.044	0.004
		(0.106) (0.085) (0.082)	(0.074) (0.082) (0.087) (0.109)	(0.128) (0.207) (0.505)

Notes:

The table reports slope coefficients from panel regressions along with the corresponding standard errors. (i) standard error in parentheses; (ii) *** p<0.01, ** p<0.05, * p<0.1; (iii) 6,788 observations; (iv) cntm -Containment and Health Index, infnew infection cases, durduration, credcredit rating, mmrmoney market rate, caryield-based carry, cxconvexity, sizebond market value, mom -momentum, revreversal, and dummybinary variable that takes value 1 if Monday and 0 for the rest of the weekdays.

Table 4 .

 4 Panel Conditional Quantile Regression -Economic Support Index

			Lower Quantiles		Middle Quantiles			Upper Quantiles
		0.05	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
	eco		-							-
		-0.002	0.004*** -0.002*	-0.001	-0.001	-0.001	-0.001	-0.004**	0.010*** -0.025***
		(0.001) (0.001) (0.001)	(0.001) (0.001) (0.001) (0.001)	(0.002)	(0.003)	(0.007)
	inf	0.273*** 0.255*** 0.240***	0.236*** 0.239*** 0.238*** 0.224***	0.243*** 0.209*** 0.170**
		(0.016) (0.014) (0.012)	(0.011) (0.013) (0.014) (0.019)	(0.023)	(0.033)	(0.085)
	dur	-3.718*** -4.251*** -3.546*** -3.272*** -3.264*** -3.125*** -3.131*** -2.733*** -2.920*** -3.751
		(0.444) (0.378) (0.337)	(0.316) (0.362) (0.391) (0.507)	(0.614)	(0.884)	(2.282)
	cred -0.366*** -0.370*** -0.314*** -0.314*** -0.325*** -0.331*** -0.332*** -0.348*** -0.345*** -0.424***
		(0.018) (0.015) (0.014)	(0.013) (0.015) (0.016) (0.021)	(0.025)	(0.037)	(0.095)
	mmr	2.343*** 2.546*** 2.572***	2.594*** 2.664*** 2.721*** 2.814***	2.902*** 3.006*** 3.394***
		(0.043) (0.036) (0.032)	(0.030) (0.035) (0.038) (0.049)	(0.059)	(0.086)	(0.221)
	car	2.348*** 2.463*** 2.351***	2.437*** 2.555*** 2.654*** 2.777***	3.006*** 3.217*** 4.479***
		(0.066) (0.056) (0.050)	(0.047) (0.054) (0.058) (0.075)	(0.091)	(0.132)	(0.340)
	cx	0.187*** 0.226*** 0.191***	0.176*** 0.178*** 0.172*** 0.172***	0.148*** 0.161*** 0.215
		(0.025) (0.022) (0.019)	(0.018) (0.021) (0.022) (0.029)	(0.035)	(0.051)	(0.132)
	size	-0.331*** -0.419*** -0.468*** -0.544*** -0.640*** -0.757*** -0.840*** -0.982*** -1.175*** -1.347***
		(0.050) (0.042) (0.038)	(0.035) (0.041) (0.044) (0.057)	(0.069)	(0.100)	(0.258)
	mom -0.542*** -0.697*** -0.630*** -0.660*** -0.770*** -0.850*** -0.928*** -1.016*** -1.392*** -2.258***
		(0.077) (0.066) (0.059)	(0.055) (0.063) (0.068) (0.089)	(0.107)	(0.155)	(0.399)
	rev	0.475*** 0.388*** 0.383***	0.413*** 0.416*** 0.430*** 0.416***	0.377*** 0.427*** 0.442**
		(0.034) (0.029) (0.026)	(0.024) (0.028) (0.030) (0.039)	(0.047)	(0.069)	(0.177)
	dummy 0.001	-0.068	-0.085	-0.107	-0.196** -0.175** -0.185	-0.019	-0.012 -0.100
		(0.101) (0.086) (0.076)	(0.072) (0.082) (0.089) (0.115)	(0.139)	(0.201)	(0.519)

Notes:

The table reports slope coefficients from panel regressions along with the corresponding standard errors. (i) standard error in parentheses; (ii) *** p<0.01, ** p<0.05, * p<0.1; (iii) 6,788 observations; (iv) eco -Economic Support Index, infnew infection cases, durduration, credcredit rating, mmrmoney market rate, caryield-based carry, cxconvexity, sizebond market value, mom -momentum, revreversal, and dummybinary variable that takes value 1 if Monday and 0 for the rest of the weekdays.

Table 5 .

 5 Panel Conditional Quantile Regression -Robustness Analysis Using Daily USD Returns as a Proxy for Bond Price VolatilitySecond, we work with alternative sets of control variables, and we show a similar effect of governmental interventions and COVID-19 related uncertainty (Table6).

	Notes: The table reports slope coefficients from panel regressions along with the corresponding standard errors.
	(i) standard error in parentheses; (ii) *** p<0.01, ** p<0.05, * p<0.1; (iii) 6,788 observations; (iv) gvt -government

response index, infnew infection cases, durduration, credcredit rating, mmrmoney market rate, caryield-based carry, cxconvexity, sizebond market value, mom -moo=mentum, revreversal, and dummybinary variable that takes value 1 if Monday and 0 for the rest of the weekdays.

6

In an unreported analysis we also consider different nested models. The major results remain unaffected.

Table 6 .

 6 Panel Conditional Quantile Regression -Robustness Analysis Using a Different Set of Control Variables

			Lower Quantiles		Middle Quantiles			Upper Quantiles
		0.05	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
	gvt	0.009*** 0.015*** 0.014***	0.012*** 0.012*** 0.016*** 0.018***	0.018*** 0.011** -0.002
		(0.003) (0.002) (0.002)	(0.002) (0.002) (0.002) (0.003)	(0.004) (0.005) (0.012)
	inf	0.095*** 0.060*** 0.042***	0.046*** 0.043*** 0.031*	0.018	0.017	0.026	-0.133
		(0.025) (0.017) (0.016)	(0.014) (0.016) (0.017) (0.022)	(0.028) (0.044) (0.102)
	dur	-0.126** -0.223*** -0.200***	-0.168*** -0.182*** -0.158*** -0.166***	-0.188*** -0.188* -0.084
		(0.055) (0.039) (0.035)	(0.031) (0.036) (0.038) (0.050)	(0.063) (0.097) (0.227)
	cred	0.153*** 0.205*** 0.225***	0.252*** 0.258*** 0.262*** 0.273***	0.261*** 0.284*** 0.218**
		(0.022) (0.015) (0.014)	(0.012) (0.014) (0.015) (0.020)	(0.025) (0.039) (0.090)
	car	0.036	-0.067	-0.052	0.006	0.138*** 0.232*** 0.396***	0.702*** 1.044*** 2.473***
		(0.071) (0.050) (0.045)	(0.040) (0.046) (0.048) (0.064)	(0.081) (0.125) (0.290)
	size	-0.801*** -0.876*** -0.881***	-0.929*** -1.008*** -1.139*** -1.232***	-1.352*** -1.615*** -1.762***
		(0.061) (0.043) (0.039)	(0.034) (0.040) (0.042) (0.056)	(0.070) (0.108) (0.252)
	mom	0.967*** 0.924*** 0.935***	0.851*** 0.734*** 0.599*** 0.470***	0.260** -0.348** -1.264***
		(0.091) (0.064) (0.058)	(0.051) (0.059) (0.062) (0.083)	(0.104) (0.160) (0.371)
	rev	1.242*** 1.138*** 1.140***	1.102*** 1.125*** 1.190*** 1.211***	1.249*** 1.364*** 1.473***
		(0.037) (0.026) (0.023)	(0.020) (0.024) (0.025) (0.033)	(0.042) (0.065) (0.151)
	dummy 0.007	-0.036	-0.092	-0.091	-0.126	-0.141* -0.116	-0.061	0.112	0.080
		(0.124) (0.087) (0.080)	(0.070) (0.080) (0.085) (0.113)	(0.142) (0.218) (0.507)

Notes:

The table reports slope coefficients from panel regressions along with the corresponding standard errors. (i) standard error in parentheses; (ii) *** p<0.01, ** p<0.05, * p<0.1; (iii) 6,788 observations; (iv) gvtgovernment response index, infnew infection cases, durduration, credcredit rating, caryield-based carry, sizebond market value, mom -momentum, revreversal, and dummybinary variable that takes value 1 if Monday and 0 for the rest of the weekdays.

Table A2 .

 A2 Major Variables Used in the Study government policy response index aggregating all types of policies and rescaling them to create a score between 0 and 100 on day t cntm Containment and Health index COVID-19 containment and health index aggregating only containment, closure, and health policies and is rescaled to create a score between 0 and 100 on day t. stg Stringency Index COVID-19 containment and health index aggregating only containment and closure, policies, and is rescaled to create a score between 0 and 100 on day t. eco Economic Support Index COVID-19 economic support index aggregating government only the policy responses targeted and providing economic support and is rescaled to create a score between 0 and 100 on day t. inf New infections The new cases of infection are computed as ln(1+ΔINF'), where INF' is the number of infected cases. Panel C: Control variables dur Duration Average adjusted duration of the bond market index on day t-1. cred Quantified credit rating Numerical sovereign rating of the government bonds in the index on day t-1. The credit risk score for each market is calculated as the average numerical rating from three major rating agencies: Moody's, S&P, and Fitch. To obtain the numerical ratings, we convert all the ratings linearly, so that the top rating (AAA/Aaa) is associated with 1, and the bottom rating (C) is associated with 21. mmr Money market rate Three-month interbank rate that is available in a given country at t-1. car Carry The difference between the yield-to-maturity on 10-year government bonds and the 3-month interbank interest rate cx Convexity Average adjusted convexity of the bond market index on day t-1. size Market value Natural logarithm of the market value of the bond index portfolio expressed in U.S. dollars on day t-1. mom Momentum Change in the yield to maturity level on the government bond index in months t-12 to t-1. rev Reversal Change in the yield to maturity level on the government bond index in months t-60 to t-13. dummy "Monday effect" dummy The variable takes value 1 if the day of the week is Monday and 0 otherwise.

	Symbol	Variable	Description
		Panel A: Dependent variables
	R1	Daily absolute residuals from	R1 represents the residuals from the seven-factor model (1), computed
		a seven-factor model	as |ln(1+R)|.
	R2	Daily absolute returns in U.S.	R2 represents the daily returns computed as |ln(1+R)|.
		dollars
		Panel B: Explanatory variables of interest
	gvt	Government Response index	COVID-19

Notes:

The table presents the variables that are used in the study.

Our focus is on sovereign bond market volatility. A concurrent strand of the literature (e.g.[START_REF] Andrieș | The COVID-19 pandemic and sovereign bond risk[END_REF][START_REF] Augustin | In sickness and in debt: The COVID-19 impact on sovereign credit risk[END_REF][START_REF] Cevik | Contagion of fear: Is the impact of COVID-19 on sovereign risk really indiscriminate?[END_REF] Daehler et al., 2021;[START_REF] Pan | The COVID-19 pandemic and sovereign credit risk[END_REF] investigates the pandemic's effect on sovereign bond risk.

Most of the existing works on this topic focus on the first wave of pandemics (March-May 2020). In our opinion, the second wave of the sanitary crisis (September-November 2020) is equally important to study the impact of COVID-19 on sovereign bond market volatility, given the additional measures imposed by governments to fight against the pandemics. However, the study of the third wave of pandemics (February-March 2021) should be placed in a totally different context given the start of the vaccination campaign.

Http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html (accessed 03 November 2020).

Https://www.bsg.ox.ac.uk/research/research-projects/coronavirus-government-response-tracker (accessed 03 November 2020).

This work was supported by a Grant of the Romanian National Authority for Scientific Research and Innovation, CNCS-UEFISCDI, Project Number PN-III-P1-1.1-TE-2019-0436. We are grateful to David Y. Aharon and Adam Zaremba for their inputs on the earlier versions of this paper.

Appendix Table A1. Construction of the Cross-Sectional Asset Pricing Factors

Symbol

Factor Description MKT F Market risk factor MKT F is the excess return on the market, i.e., the value-weighted return of all the bond indices in the sample at the end of month t minus the risk-free rate, i.e., the one-month T-Bill return.

DUR F Duration factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the highest (lowest) adjusted duration.

CRED F Credit risk factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the highest (lowest) adjusted credit risk score. The credit risk score for each market is calculated as the average numerical rating from three major rating agencies: Moody's, S&P, and Fitch. To obtain the numerical ratings, we convert all the ratings linearly so that the top rating (AAA/Aaa) is associated with 1, and the bottom rating (C) is associated with 21.

SIZE F Size factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the highest (lowest) market value of the relevant bond basket.

MOM F Momentum factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the lowest (highest) change in yield-tomaturity from t-12 to t-1. This corresponds with going long (short) bonds with the highest (lowest) return induced by the change in YTMs.

REV F Reversal factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the highest (lowest) change in the yieldto-maturity (YTM) from t-60 to t-13. This corresponds with going long (short) bonds with the lowest (highest) return induced by the change in YTMs.

CAR F Carry factor

The factor is represented by a long-short zero-investment portfolio that buys (sells) the valueweighted portfolio comprising 30% of bond indices with the highest (lowest) lowest carry. The carry variable is measured as the difference between the yield-to-maturity on 10-year government bonds and the 3-month interbank interest rate.

Notes: The table displays the procedures used to calculate the returns on asset pricing factors used in this study.


domain: phys spi

text:
Pierre Delullier
email: pierre.delullier@onera.fr

Guillaume Druart

Florence De

La Barrière

Laurent Calvez

Matthieu Lancry

FEMTOSECOND LASER DIRECT WRITING OF A GRADIENT INDEX FRESNEL LENS FOR MID-INFRARED IMAGING APPLICATIONS

Keywords: Optical design, Gradient Index, Laser Writing

An important challenge in infrared imaging today consists in addressing the SWaP problem (Size, Weight and Power), for example by simplifying as much as possible the optical system before the sensor. This work presented takes advantage of recent techniques in femtosecond laser direct writing to imprint optical systems. Here we want to simplify an infrared multispectral imaging system, which combines a lens array and a filter array. Each channel has a focal length of 7 𝑚𝑚 and a f-number of 4.

We showed that reaching a phase variation of 𝜋 radians with femtosecond laser writing is possible even for a huge scanning speed (0,01 𝑚𝑚/𝑠 up to 25 𝑚𝑚/𝑠). By overlapping the written lines, we were able to implement a maximum phase variation of 8 radians (with only one written layer), demonstrating the possibility to write quickly a GRIN Fresnel lens with a diameter of some millimeters.

INTRODUCTION

The SWaP problem (Size, Weight and Power) is an important challenge today in infrared imaging. In order to simplify the optical system before the sensor and to keep high image quality, optical designers have to find new degrees of freedom in their designs. Several ways are thus explored like freeform optics [START_REF] Wang | Optimized multielement accommodative intraocular lens with a four-freeform-surface Alvarez lens and a separate aspheric lens[END_REF], GRIN (gradient index) optics [START_REF] Yee | Design of a freeform gradient-index prism for mixed reality head mounted display[END_REF] and metasurfaces [START_REF] Raulot | Méthodes de conception et de fabrication de dispositifs imageurs en optique diffractive à structures sub longueur d'onde[END_REF].

GRIN optics involve creating a spatial variation of the refractive index in the material. This index variation can correct some optical aberrations such as spherical aberration or chromatic aberrations [START_REF] Mait | Dispersion design in gradient index elements using ternary blends[END_REF] or encodes the optical power inside the lens (not just at surfaces).

One of the main limitations of GRIN optics is their manufacturing [START_REF] Gibson | Multispectral IR optics and GRIN[END_REF]. Most of the fabrication techniques can only create classic index variation (axial or radial) but freeform or even aspheric index profiles are very difficult to create. Some solutions exists in visible materials with thin film coating but it is much more complicated for infrared materials. This paper focuses on femtosecond laser direct writing, whose main advantage is to be very flexible. By focusing an ultrafast laser in a substrate, this technique can change locally the density and, as a consequence, the refractive index of the material. It is already well known in visible and near infrared (IR) materials such as silica [START_REF] Beresna | Ultrafast laser direct writing and nanostructuring in transparent materials[END_REF] but few studies concern materials used in mid-IR optics (chalcogenide) [START_REF] Mcmillen | Ultrafast laser fabrication of lowloss waveguides in chalcogenide glass with 0.65 dB/cm loss[END_REF]. This paper presents how to simplify a multispectral camera by merging a lens array with a filter array, with gradient index lenses obtained by femtosecond laser writing inside a dedicated glass substrate.

A GRADIENT INDEX PLATE IN A MULTISPECTRAL CAMERA

Multispectral camera

ONERA has developed a multichannel multispectral camera whose concept is illustrated in Figure 1. It consists in an array of 2×2 lenses with a focal length of 7 mm and a F-number of 4. Each channel has a field of view of 40°×30° and forms a 320×256 pixels image on each quarter of the detector. An array of 2x2 filters is placed after the lens array, so that each optical channel corresponds to a different spectral band with a width of about 300 nm.

Figure 1. Cross-section of the ONERA multispectral camera. From top to bottom, the diaphragm, the lens array, the filter array and the sensor.

The optical system is composed of two objects, a lens array for imaging and a filter array to select a specific spectral band for each channel. These components are integrated in a cold shield of a Detector Dewar Cooler Assembly (DDCA) and are cooled at cryogenic temperature. These additional masses tend to increase the cool-down time of the DDCA. In order to address the SWaP problem for this camera, we want to remove one of the components. The filter array is an interference filter made by a multiple thin layers deposition on a substrate. It is difficult to depose these films on a non-planar surface such as the lens array. So, the solution is to use a unique GRIN optic plate for the focusing and filter functions. In the rest of the paper, we will study only one channel of the multichannel camera.

A GRIN focusing plate needs a radial index variation to create an optical power instead of the spherical surfaces. A classical index profile used in this case is described in Eq. 1:

𝑛(𝑟) = 𝑛 0 + 𝑛 𝑟 𝑟 2 (1)
With 𝑛 0 the index on the optical axis, r the radial position and 𝑛 𝑟 the coefficient of index variation. This index distribution gives to the optical plate a vergence V [START_REF] Boyd | Optical design of multi-material, general rotationally symmetric GRIN lenses[END_REF]:

𝑉 = -2𝑛 𝑟 𝑡 = - Δ𝑛 𝑅 2 𝑡 (2)
With 𝑡 the thickness of the plate, Δ𝑛 the index variation and 𝑅 the radius of the lens. The issue with this index profile is to have a short focal length (ie a huge vergence). The vergence is proportional to the index variation (Eq. 2) so a short focal length as the focal length of the studied camera needs an important index variation. For a 7 𝑚𝑚 focal length in our case, this type of design needs a gradient index with a maximum variation of 𝛥𝑛 = 0,11 for a thickness of 1 mm that is higher than what a femtosecond laser direct writing can create. One of the maximum index variation reached today in mid-IR components with the femtosecond laser writing technique is about 0.05 in As2S3 (a chalcogenide glass) [START_REF] Zoubir | Direct femtosecond laser writing of waveguides in As2S3 thin films[END_REF].

GRIN Fresnel lenses

In order to reduce the index variation and preserve the small thickness of this component, the purpose is to code a Fresnel GRIN lens instead of a classical GRIN lens. This involves cutting a continuous surface into a series of discontinuous rings. The Fresnel versions of a spherical lens and GRIN lens are shown in Figure 2. Given that a GRIN index profile is equivalent to a surface sag of a conventional lens, the modeling of a Fresnel lens can be done in the Zemax software with a "User Define Surface" created for conventional Fresnel lenses. The surface sag used in the software is then converted into an index profile distribution. The surface sag (𝑧 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 ) equation is described in Eq. 3 [START_REF] Grulois | Etude de l'apport des lentilles de Fresnel pour la vision[END_REF].

{ 𝑧 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 (𝑟) = 𝑐𝑟 2 1 + √(1 -(1 + 𝑘)𝑐 2 𝑟 2 mod 𝐻 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 𝐻 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 = 𝜆 0 𝑛 0 -1 (3)
𝑐 is the curvature, 𝑘 the conic constant, mod is the modulo operation, 𝜆 0 is the central wavelength and 𝐻 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 is the etching depth of a first order Fresnel lens. In the case of a GRIN Fresnel lens, the surface sag is replaced by an index profile distribution (𝑛 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 ) and so 𝐻 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 corresponds to the total index variation Δ𝑛 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 :

{ 𝑛 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 (𝑟) = 𝑐𝑟 2 1 + √(1 -(1 + 𝑘)𝑐 2 𝑟 2 mod Δ𝑛 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 Δ𝑛 𝐹𝑟𝑒𝑠𝑛𝑒𝑙 = 𝜆 0 𝑛 0 -1 (4)
By using these Eq. 3 and Eq. 4, the Fresnel GRIN lens is modeled in Zemax and optimized for the camera specifications. Only one channel is modeled and represented in Figure 3 because each channel is the same with a small adaptation of the total index variation as a function of the central wavelength of the filter. 16 levels of discretization could be better, but we would have rings thinner than 1 µ𝑚, which cannot be reached with our laser direct writing technique. We modeled a GRIN Fresnel lens in our material, which is presented in Figure 5. The total gradient index is about 3,7.10 -2 in the material which is reachable with our laser writing technique. With the numerical aperture specification of the camera, the last ring needed is about 30.7 µ𝑚 thick and so we need a spatial resolution of 3.8 µ𝑚 to be able to write this last ring discretized into 8 levels.

FEMTOSECOND LASER DIRECT WRITING OF GRADIENT INDEX

Identification of permanent modifications in glasses

The laser-material interaction must be characterized to be able to reproduce the index profile shown in Figure 5. We studied the chalcogenide glass 75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙). There are different types of permanent modifications due to laser writing. First, the type I modification is a permanent isotropic index change due to a local volume variation [START_REF] Cao | Pulse energy dependence of refractive index change in lithium niobium silicate glass during femtosecond laser direct writing[END_REF]. The second type of effect (type II) corresponds to a permanent anisotropic index change associated to the formation of nanogratings [START_REF] Drevinskas | Giant birefringence and dichroism induced by ultrafast laser pulses in hydrogenated amorphous silicon[END_REF]. The neutral axis of the birefringence created depends on the laser polarization. The last type of permanent change in the material is the heat accumulation which appears with high energy. The phenomenon is the same as type I but there is a melting and quenching of the matter causing an increase of the line size bigger than the laser spot ones.

For a GRIN Fresnel fabrication, we are especially interested in type I modifications. We want the smallest line size, for a better spatial resolution, and a huge index variation. Several parameters, such as laser frequency, pulse duration, pulse energy or scanning speed can influence the index variation type. Figure 6 shows the appearance of different types of modification as functions of writing speed and laser energy. As shown in Figure 6, the faster the writing is the higher the energy needed is to write a permanent change into the substrate. We note that our material has a huge range of energies leading to type I changes even for high writing speed (up to 25 𝑚𝑚/𝑠). It confirms that femtosecond laser direct writing is able to write GRIN optics very quickly. The heat accumulation could also be a way to increase the index variation and to write lines bigger than the laser spot ones enabling faster writing for big surfaces irradiations.

Phase variation with laser energy and writing speed

GRIN Fresnel optics are more easily described with their phase variation than with their index variation because Fresnel optics are defined by a phase variation 𝛥𝜑 = 2𝜋 𝑟𝑎𝑑 for a first order Fresnel lens.

To write a GRIN Fresnel plate, we needed to know the precise phase variation implied by the laser as function of the speed writing and the laser energy. So, several laser-written discs were made at different energies and speeds in the type I range. To write circular lenses or circular objects like monophase discs, a spiral trajectory technique was used. Presented in Figure 7, the writing has a tight step of only 1 µ𝑚 avoiding any groove appearance. These discs were then studied with a Quantitative Phase Measurement (QPM) technique [START_REF] Bélanger | Comparative study of quantitative phase imaging techniques for refractometry of optical waveguides[END_REF]. This technique consists in using the transport of intensity equation (TIE) which relies on intensity variation and phase variation. Using an optical microscope, a set of images was taken: one in-focus and two very slightly (a few microns) positively/negatively defocused [START_REF] Zuo | Transport of intensity equation: a tutorial[END_REF]. A phase image given by QPM of a laser-written simple disc is shown in Figure 7. Figure 7 illustrates the technique capability to write a huge phase variation into a material even at 1 𝑚𝑚/𝑠 scanning speed. Higher writing speeds were tried (up to 2.5 mm/s) with success but a dedicated setup would be more appropriate for accurately writing complex shapes. Therefore, we use only this 1 𝑚𝑚/𝑠 writing speed. For many applications it would be necessary to convert phase variation (𝛥𝜑 in rad) into refractive index changes written by the femtosecond laser. In this view, it is sufficient to measure the thickness of the laser-modified region 𝑑 to have an approximation of the index variation:

Δ𝑛 = 𝜆 2𝜋𝑑 Δ𝜑 (6)
𝜆 is the wavelength used for the quantitative phase measurement. Thus the index variation implemented by the fs-laser in the Figure 7 is about 2.7 × 10 -2 .

Optical losses pre and post irradiation

Another important characterization needed before the writing of a Fresnel lens is the comparison of optical losses pre and post laser irradiation. If our femtosecond laser direct writing technique creates a lot of losses, it would be inappropriate to make an infrared imaging system. So we decided to measure the optical transmission of our material before and after writing a disc which introduces a phase variation of 𝜋 𝑟𝑎𝑑. Two spectrophotometers were used to measure the all IR band: an UV-Vis-NIR spectrophotometer (0.2 -2.5 µ𝑚) and a FTIR spectrophotometer (2 -20 µ𝑚). The discs written in our material had a 1.5 𝑚𝑚 diameter. The Figure 8 represents the spectral transmission of the 75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙) material pre and post irradiation. We remark that the difference before and after irradiation is negligible. The writing of a disc of 𝜋 𝑟𝑎𝑑 doesn't affect the transmission. That confirms the capability of this femtosecond laser direct writing technique to make a GRIN Fresnel plate.

CONCLUSION

In this article, we demonstrated the interest of using gradient index optics and especially gradient index Fresnel lenses to simplify an IR multispectral camera. With a GRIN Fresnel lenses matrix written in a plate, we can reduce the size and the weight of optics integrated in a dewar cooler assembly by fusing GRIN Fresnel lenses inside the array of filters. The femtosecond laser direct writing technique has demonstrated the possibility to make quickly these GRIN Fresnel lenses in a chalcogenide material. Some more trials are necessary to write a first GRIN Fresnel prototype as for instance some experiment concerning movement control. Within future work, we plan to make a first prototype using these GRIN optics and assess its performances experimentally.

Figure 2 .

 2 Figure 2. From a classical lens to a Fresnel lens and from a GRIN lens to a GRIN Fresnel lens.

Figure 3 .

 3 Figure 3. On the top, Zemax 2D layout of a Fresnel lens. Below, the Modulation Transfert Function (MTF) and the spot diagram given for three different field of view (0°, 10° and 20°).

From

  the results given by Zemax modelling, we can deduce the GRIN profile of our GRIN Fresnel plate. With the material 75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙) used in the[7,3 ; 7,5] µ𝑚 bandwidth, we obtain a GRIN lens with a diameter of 3.7 𝑚𝑚 and conic constant of -3.6. This GRIN Fresnel lens is composed of 31 rings with the outer ring has a width of approximately 30.7 µ𝑚. The laser direct writing technique will discretize the profile into several index levels so we need to determine a compromise between the imaging quality of the Fresnel plate and the spatial resolution of the writing technique. For a first order Fresnel lens as here, we can define the diffraction efficiency as a function of the index versus wavelength variation (Eq. 5)[START_REF] Raulot | Méthodes de conception et de fabrication de dispositifs imageurs en optique diffractive à structures sub longueur d'onde[END_REF]. the reference wavelength, 𝑁 is the number of discretization levels. 𝑛(𝑟, 𝜆) is the index at the radial position 𝑟 and the wavelength 𝜆.

Figure 4 .

 4 Figure 4. Graphic representation of the diffraction efficiency as a function of the wavelength for a first order Fresnel lens. The chromatic dispersion is disregarded: 𝑛(𝑟, 𝜆) = 𝑛(𝑟, 𝜆 0 ).

Figure 4

 4 Figure4shows a good diffraction efficiency (up to 90% for the [7.3 ; 7.5] µ𝑚 spectral bandwidth) for a first order Fresnel lens discretized in 8 index levels.

Figure 5 .

 5 Figure 5. Index profile of the GRIN Fresnel lens made in the chalcogenide glass (75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙)): blazed on the top and discretized on the bottom.

Figure 6 .

 6 Figure 6. Threshold identification in the 75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙) material as a function of laser energy and writing speed. These values were found with a pulse duration of 800 fs and a frequency of 100 kHz.

Figure 7 .

 7 Figure 7. Phase map of a laser-written disc. The diameter of the disc is 100 µ𝑚 and was written with a 1 𝑚𝑚/𝑠 scanning speed.

Figure 8 .

 8 Figure 8. Transmission before and after writing a 1,5 𝑚𝑚 diameter disc with a phase variation of 𝜋 𝑟𝑎𝑑 in the 75(𝐺𝑒𝑆 2 ) -15(𝐼𝑛 2 𝑆 3 ) -10(𝐶𝑠𝐶𝑙) material.


Hussein El Khansa

Carmen Gervet

Audrey Brouillet

On the ranking of variable length discords through a hybrid outlier detection approach

Keywords: Matrix profile, Prominent discord discovery, point and patter outlier detection, time series, Climate data

In this paper we are interested in identifying insightful changes in climate observations series, through outlier detection techniques. Discords are outliers that cover a certain length instead of being a single point in the time series. The choice of the length can be critical, leading to works on computing variable length discords. This increases the number of discords, with potential overlapping, subsumption and reduced insightful results. In this work we introduce a hybrid approach to rank variable length discords and extract the most prominent ones, that can yield more impactful results. We propose a ranking function over extracted variable length discords that accounts for contained point anomalies. We investigate the combination of pattern wise anomaly detection, through the Matrix Profile paradigm, with two different point wise anomaly detectors. We experimented with MAD and PROPHET algorithms based on different concepts to extract point anomalies. We tested our approach on climate observations, representing monthly runoff time series between 1902 and 2005 over the West African region. Experimental results indicate that PROPHET combined with the Matrix Profile method, yields more qualitative rankings, through an extraction of higher values of extreme events within the variable length discords.

Introduction

Outlier detection is a field of machine learning to analyze data, covering numerous application domains including fraud detection, insurance [START_REF] Hansson | Insurance fraud detection using unsupervised sequential anomaly detection[END_REF],medical [START_REF] Zhang | Unsupervised deep anomaly detection for medical images using an improved adversarial autoencoder[END_REF],internet of things [START_REF] Sgueglia | A systematic literature review of iot time series anomaly detection solutions[END_REF], cyber security and hydrology [START_REF] Siniosoglou | A unified deep learning anomaly detection and classification approach for smart grid environments[END_REF]. Overall, an outlier can be generally defined as an observation that is significantly dissimilar to other data observation or an observation that does not behave like the expected typical behavior of the other observations. An outlier detection method can be applied to multiple data types such as images, transaction data, sequence data including genomics and time series. There are two main types of anomaly detection, point wise and pattern also called collective anomalies or discords [START_REF] Borges | Anomaly detection in time series[END_REF]. Discords denote the most unusual time series sub sequences, and are detected using similarity measures that compare sub sequences with each other for a given length. Existing approaches can be either approximate (e.g. HOT-SAX [START_REF] Keogh | Hot sax: Efficiently finding the most unusual time series subsequence[END_REF], QUICK MOTIF [START_REF] Li | Quick-motif: An efficient and scalable framework for exact motif discovery[END_REF], Rare Rule Anomaly [START_REF] Senin | Time series anomaly discovery with grammar-based compression[END_REF]) or exact like the matrix profile (MP) paradigm [START_REF] Yeh | Matrix profile i: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[END_REF], used to detect motifs or discord. The choice of the sub sequence length is critical, set as input parameter and often specified with experts' knowledge which clearly influences the outcome. This observation lead more recently to works on computing variable length discords, over a given length interval (e.g. GraphAn [START_REF] Boniol | Graphan: Graph-based subsequence anomaly detection[END_REF], PanMatrix [START_REF] Madrid | Matrix profile xx: Finding and visualizing time series motifs of all lengths using the matrix profile[END_REF]). The increased usage of these approaches has also raised their weaknesses. Indeed, the number of extracted discords is larger, with potential overlapping and subsumption among them, leading to the question of "how can we select the more relevant or actionable ones?". We are not aware of approaches seeking a ranking of variable lengths discords to this date. In this work, we address this issue in the field of climate data analysis [START_REF] Philomène | Improved regional frequency analysis of rainfall data[END_REF]. Climate impact models provide time series based both on historical data and projected data generated from complex simulations of physical processes. The impact models concern many fields, including yields of agricultural crops, biodiversity, runoff of water infiltration. Such models are commonly analyzed using multi-model ensembles means to mitigate uncertainties. With a primary goal of estimating future trends, the use of outlier detection methods has seen limited interest in this field to this date, with some works on hydrological impact models [START_REF] Yu | Time series outlier detection based on sliding window prediction[END_REF][START_REF] Ye | Hydrologic time series anomaly detection based on flink[END_REF][START_REF] Yu | Time series outlier detection based on sliding window prediction[END_REF]. Another particularity of these models is the very large scale of the spatio-temporal data series, covering daily or monthly data over 150 years for each spatial point in a given region. Our objective is to investigate a ranking approach over variable length discords, by integrating two unsupervised anomaly detection frameworks, point anomaly detectors and pattern based anomaly detection using MP. The intuition behind our work, is that a discord will dominate another one, if proportionally to its length, it also contains anomalous points that are greater in numbers and extreme values within the discord sequence. In our work, we compare two types of point outlier detectors, model and predictor based and draw conclusions on their insights to rank variable length discords. We consider respectively the MAD and PROPHET methods. We also define and specify a ranking function over variable length discords. MAD is a general point anomaly detection method, based of how far a point is from the median, while PROPHET is more tailored to time series and detects anomaly points that diverge from the expected values with respect to the projected seasonality and trends derived in the learning phase of PROPHET. Our approach is evaluated on total (surface + subsurface) runoff observations, that reflect the soil water levels. Our main contributions are: 1) a hybrid approach to rank variable length discords, through the combination of exact discord method with point detectors and a ranking function; and 2) the evaluation through an experimental study over large scale climate impact data, showing the differences among point anomalous detectors and their added value as an effective means to rank variable discords. This paper is organized as follows. Section 2 gives a background on the concepts and methods we use, section 3 presents our approach. We conclude in section 4.

Background

In this section we review existing concepts and approaches we make use of, relative to outlier detection in time series, and review the use of point outliers to extract pattern outliers. Point outliers over time series can be specified as global outliers when compared to all values in the time series, and local outliers when compared to its neighboring points. Subsequence or pattern outliers in time series, extract collective consecutive points that behave unusually compared to other subsequences of similar length. In such cases a single point in the sequence might or not be also an anomaly on its own. Methods for point outlier detection can be classified as model outlier detectors or predictor based outlier detectors. Model or statistical outlier detectors derive a probabilistic model which captures the distribution of the time series. If an instance in the time series has a low probability of belonging to the estimated model, it is flagged as an outlier [START_REF] Chandola | Outlier detection: A survey[END_REF]. Examples of model based point outlier detectors are Median Absolute Deviation (MAD) [START_REF] Iglewicz | How to detect and handle outliers[END_REF] and Minimum Covariance Determinant (MCD) [START_REF] Peter | A fast algorithm for the minimum covariance determinant estimator[END_REF]. Predictor based outlier detectors proceed differently. They construct a predictor (like forecast model or a regressor) to learn the normal behavior of the time series, and then exploits the learned information to predict future values. Prediction errors which are predictions that significantly deviate from the true value are labeled as point outliers [START_REF] Li | Anomaly scoring for prediction-based anomaly detection in time series[END_REF]. Many predictor methods have been developed including ARIMA [START_REF] Sanayha | Fault detection for circulating water pump using time series forecasting and outlier detection[END_REF], and SVM [START_REF] Ma | Time-series novelty detection using one-class support vector machines[END_REF] and PROPHET [START_REF] Sean | Forecasting at scale[END_REF]. In this work, we will apply one point outlier from each type to extract point outliers within computed variable length discords, as a means to rank those discords. Note that point outliers have also been used to directly segment data series and thus extract discords. We briefly survey these approaches and their suitability for our goal.

Point outliers detectors and series segmentation

In [START_REF] Ding | An anomaly detection approach based on isolation forest algorithm for streaming data using sliding window[END_REF] the authors proposed IForest-based anomaly detection to detect drift in streaming data, where a fixed sliding window is applied, and in each window the IForest detector is applied to generate the anomaly rate by averaging the instance depth of the forest, if the anomaly rate is smaller than the threshold then their is a drift. Also, in [START_REF] Chen | Unsupervised anomaly detection of industrial robots using sliding-window convolutional variational autoencoder[END_REF] the authors proposed sliding-window convolutional variational autoencoder (SWCVAE) an algorithm to detect point anomalies in robot, by applying sliding windows. Each window is fed to a convolutional variational autoencoder, a trained artificial neural network(ANN), to detect point wise anomaly. The ANN will give each instance in the sliding window an anomaly score, then these scores are averaged to calculate the anomaly score of the discord. If an anomaly score is above the predefined threshold then it is classified as a discord. In [START_REF] Yue | An integrated anomaly detection method for load forecasting data under cyberattacks[END_REF] the author combines both the result of SAX as a pattern outlier detection and combined with point outlier detection second-order-difference (SOD) and Chebyshev Inequality (CI)-based methods, where discord that does not contain point anomaly detected by SOD and CI are discarded. Those approaches contribute both ranking mechanisms of discords through a common concept, that is the segmentation of the time series and the use of a point outlier detector on each subsequence to give it an anomaly score. With respect to our goal of comparing variable length discords over large time series, these approaches offer scoring mechanisms for subsequences, but are applied to fixed size windows, and go through the whole data set. This requires the application of outlier detectors over very large sets of sliding windows, without considering variable length windows. In our approach, we purposely focus on variable length discords, and separate the issue of computing variable length discords, through an exact dedicated approach, with that of exploiting point anomaly detectors for ranking purpose. Part of our goal is to analyse the contribution and added value of point anomaly detectors to rank the discords of variable lengths.

Matrix Profile, MAD and PROPHET

In this subsection we provide the necessary background on the three methods used in our work, namely the STOMP algorithm to compute Matrix Profiles and extract discords, and the point anomaly detectors MAD and PROPHET.

Matrix Profile methods Matrix Profile is a powerful data structure that stores the z-normalized Euclidean distances between each subsequence and its closest neighbor sequence.Two main algorithms have been defined to compute the exact solution of the matrix profile by performing all-pairs-similarity-search on time series, namely STAMP [START_REF] Yeh | Matrix profile i: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[END_REF] and STOMP [START_REF] Zhu | Matrix profile ii: Exploiting a novel algorithm and gpus to break the one hundred million barrier for time series motifs and joins[END_REF]. We adopted STOMP for its runtime and scalability over large time series. STOMP computes the matrix profile by calculating the distance profile of every subsequent in time series T and then selecting the minimum value in each distance profile. It uses a similarity search algorithm to compute the matrix profile faster in time complexity of O(n 2 ).The produced Matrix can be used to detect discords, where the highest K subsequences in the Matrix Profile profile are the top K discords, specifying the highest distance of these subsequences with its closest neighbor. The concepts used in this work are recalled here after, relative to discords and the Matrix Profile data structure.

Definition 1. A time series T is a sequence of real-valued numbers t i : T = [t 1 , t 2 , ..., t n ]
where n is the length of T .

Definition 2 (Matrix Profile).

A matrix profile P m of time series T and given length m is a meta series of the Euclidean distances vector between each subsequence T i,m of given length m where i varies, and its nearest neighbor (closest match) in time series T , together with the corresponding position vector for each closest neighbor associated with min(D i,m ). We denote it P m = [min(D 1,m ), ..., min(D n-m+1,m )], where D i,m (1 ≤ i ≤ n -m + 1) is the distance profile D i,m of time series T for subsequences of length m. Definition 3 (Discord). The discord denoted ∆ j,m is a subsequence T j,m of length m starting at the position j, that has the largest distance to its nearest neighbor. In matrix profile the largest distance corresponds to the maximum distance value in P m . MAD The Median Absolute Deviation outlier detection method is a point outlier detector that uses the Absolute deviation around the median to measure the distance between a data point and the median [START_REF] Iglewicz | How to detect and handle outliers[END_REF]. The median is considered a tendency measurement of the data similar to the mean, but unlike the mean it is more robust to detect relevant point outliers in the data, being immune to noise, and to the sample size [START_REF] Leys | Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median[END_REF]. Given a time series T = x 1 , x 1 , ..., x i , ..., x n such that i is the time unit, MAD first calculates the median of T , this median is denoted by x, and its absolute deviation, by subtracting x from each x in T :

M AD = median(|x i -x|) (1)
For each x i we calculate their modified z-score using the predefined formula [START_REF] Rand R Wilcox | Fundamentals of modern statistical methods: Substantially improving power and accuracy[END_REF]:

M i = 0.6745(x i -x) M AD (2)
In statistics, the z-score for each x i measures the difference between the value and the mean how many standard deviations away a value is from the mean in terms of standard deviation . MAD uses the following formula to calculate a z-score:

Z -Score = x i -µ σ (3)
where: µ is the mean value of the data set, σ is the standard deviation of the data set where if the Z-score is 0 then it means the value of x i is equal to the mean, if z-score is 2 then it is 2 standard deviation (2 * σ) away from the mean.

If the absolute value of the modified z-core of x i is larger than a threshold then is considered as an outlier. [START_REF] Iglewicz | How to detect and handle outliers[END_REF] proposed the value 3.5 as a threshold, also the default one in MAD.

PROPHET It is an open source time series forecasting model, that detects outliers based on the data learned through the forecasting component. Its purpose is to forecast future data, and detect outliers, with an uncertainty estimation. Once PROPHET fits the given time series, it generates a data set containing the following variables with an estimation of the uncertainty through a confidence interval for each value. Then for each value in the Time series, it calculates the error and uncertainty based on the confidence interval of each value. If the error of a value is bigger than its uncertainty it will be flagged as an anomaly.

PROPHET is based on generalized additive model (GAM), which means that the model is the sum of several components summarized below:

y t = g(t) + s(t) + h(t) + ε(t) (4)
Where:

g(t) denotes the growth function, a piece-wise linear function. The number of linear growths (straight lines), is determined by the number of change points in the time series, that can be entered manually or detected automatically by PROPHET.

g(t) = k + a(t) ⊤ δ t + m + a(t) ⊤ γ ( 5
)
where k is the growth rate scalar value, m the offset value, a(t) a binary vector ∈ {1, 0} indicating if there is a growth rate adjustments at time t. δ is a vector that contains the growth rate value. For instance if there is 10 change point then the vector of δ will have 10 values. -γ is vector similar too that contains the offset parameter adjustment s(t) the seasonality, modeled using Fourier Series as a function of time.

h(t) holiday seasons parameter, that we will not consider ε(t) the error estimation using a maximum a posteriori estimation. Not all components are mandatory, and in our case since we will be experimenting on climate impact data. h(t) not considered for instance.

Our approach

Our work is first motivated by the insightful ranking of variable length discords, to extract those that are combined with more eventful changes and potential impact on forecasted data trends. We do so by combining an exact, robust pattern outlier detection approach, with two different point outlier detectors. The exact approach ensures reliable extraction of discords, while the point outlier detectors will contribute to their ranking. We use the Matrix Profile method, STOMP as an exact approach to extract discords, that also offers the qualities of being parameter-free and that does not utilize a similarity or distance threshold. This is important since the point outlier detectors do not have these advantages, and an increase in thresholds tuning would reduce the robustness of the overall results. With respect to the point outlier detectors, we experimented with the generic MAD model approach, and the predictor PROPHET method that also handles trends and seasonality. Our Algorithm 1 is based on four 4 main phases :

1. (Lines 8-12) Extract variable length discords (Matrix Profile method) 2. (Line 13) Detect point outliers (using PROPHET, MAD) over the time series 3. (Lines 14-20) Score the discords based on the anomalous points they contain 4. (Lines 21-25) Filter and sort the discords with a ranking function Extract variable lengths discords. We compute the Matrix Profile for variable window sizes j over T , where 4 < j < length(T )/2. For each window size, we extract the top discord. The output is a list of top discords, one per length. Detect point outliers in the time series. This step seeks point wise outliers (through PROPHET and MAD), over the whole time series T . The outcome is a Boolean list denoting for each data in the time series if it is an anomaly.

Scoring the discords. This step integrates the variable length discords and Boolean anomaly list for scoring purpose. For each discord we count the anomalous points that are included in each discord ∆ i,j . Clearly a simple count is not the only element in the scoring, and the discords' length needs to be taken into account. Note that some discords can overlap, thus share anomalous points. To address this issue we define a scoring function as the ratio of the count over each length.

Filtering and sorting. This procedure filters the discords based on their count, as well as their starting dates, and length. Based on our past work and analysis [START_REF] Khansa | Prominent discord discovery with matrix profile : Application to climate data insights[END_REF] of subsumption among variable length discords, it is important to compare variable length discords with identical starting dates. The filtering stage groups discords sharing a starting date i and for each group selects the discord with highest ratio, meaning it subsumes the others. Then the discords (of variable length) are sorted in descending order of their ratio value. Climate data. In our experimental study we used observed monthly runoff data obtained from Global Runoff Reconstruction dataset (GRUN) [START_REF] Ghiggi | Grun: an observation-based global gridded runoff dataset from 1902 to 2014[END_REF]. In climate and corresponding impact data science, runoff is an impact variable that can be used to quantify flood and drought risks at regional and global scales [START_REF] Arnell | The global-scale impacts of climate change on water resources and flooding under new climate and socio-economic scenarios[END_REF]. It is usually provided in kg/m2/s or in mm/day. The source of the data has a 0.5°× 0.5°spatial grid resolution with a focus on the Sahel region. This data gives indicators on the soil water content over the time period between 1902 and 2005 (i.e. 104 years, 1248 months). This is a standard observation period in historical climate analysis. We spatially average the monthly runoff data over the grid box [5°W-25°E ; 10°N-18°N].

Point anomaly detectors.

A first comparative study between MAD and PROPHET on the runoff data set is shown in Figure 1. Those point anomaly detectors use different methods to extract anomalies, and we can see that some anomalous points are shared and many differ. Over the whole data series, MAD extracts a large set of anomalies (green), essentially in the upper part of the data, thus detecting high water levels. On the other hand, PROPHET extracts fewer anomalies, some shared with MAD for the highest values of anomalous water levels, and most importantly detected the drought periods in terms of anomalous points (purple crosses), in 1900-1910 and 1970-1990. This can be explained by the fact that PROPHET through the learning phase was able to identify a normality within seasons, whereas MAD focuses on extremes in the overall data series. If we zoom on the period 1981-1984 where 1983 is recorded as one of the driest years in the West African region, this aspect is striking. MAD does not extract any anomalous point, whereas PROPHET detects three anomalies displayed in Figure 2. By not accounting for, or detecting, seasonality MAD misses important point outliers in such time series. On the other hand, PROPHET takes into account seasonality through the learning phase, and we can see that it extracts anomalies with respect to months which are learnt to be dry and those that should not be. We note that the anomalous points do not necessarily reflect extreme values overall, but anomalous points with respect to the anticipated values (learnt during the forecasting phase). We also studied the early period 1902-1905, a period of higher runoff values in average. We can see that both MAD and PROPHET detect anomalies during this period, with some common points, showcased in Figures 3 and4. Those figures are interesting as they also how MAD does detect some anomalous points not tagged as such by PROPHET. This reflects the role played by the median function and the z-score defined by MAD. All experiments were run on an Intel(R) Xeon(R) Bronze 3106 CPU processor at 1.70GHz with 8 core with 64 GB of RAM. The main computational cost are the N/2 runs of STOMP to compute the variable length discords, taking around 9 minutes. MAD is in constant time, and PROPHET took 1.5 minutes.

Integrated approach: ranking and outcomes

Having compared the commonalities and differences between two point anomaly detectors, we saw that even though they are generic, they can have a great impact on the outliers' detection. This is related to their models and potentially the nature of the data. We now present the results of our methodology. Both Figure 5 represents respectively for PROPHET and MAD, the top discords ranking with respect to our ratio function. Note that the discords of variable lengths are specified by their starting date, window length and the value of the ratio function. Note that since we are comparing discords of variable lengths, at this stage we only extract the top discord for each given length m, to reduce the combinatorics of our approach. Extracting top k discords per length is part of future work. We recall that the ratio values depend on the length and the number of anomalous points within each discord. In this figure we picked the top 20 discords relative to each detector output ranking. The discords highlighted in yellow, correspond to a filtering, that extracts among those discords the top ones outside a proximity temporal interval of +/-three months. For instance for the year 1982 we have 3 ranked discords starting from July till September with overlapping periods. We kept the one with highest ranking to reduce redundant discord periods. Thus the colored row in the Table 5 1. Point anomalous detectors do matter. While STOMP produces exact matrix profiles to derive variable length discords, the rankings differ between PROPHET and MAD. By detecting the anomalies in the driest period, where droughts occurred (1982)(1983)(1984), PROPHET manages to rank top discords in this period. We can see that the variable length are not negligible ranging in this period from 97 months to 123 months. MAD detects mainly the periods with high water levels (1902)(1903)(1904), yielding high ranking for the corresponding discords. It corresponds to the periods with higher soil water levels and larger inter-annual variability, that comes before the later continuous longterm drying trends observed within Sahel [START_REF] Khansa | Prominent discord discovery with matrix profile : Application to climate data insights[END_REF]. PROPHET extracted and ranked the 1 st , 2 nd and 11 th discords, that represent the intense droughts that occurred in Sahel, while the 3 rd and 19 th discords represent the the period where the high soil water levels and larger inter-annual variability. 2. Rankings are insightful. We can see in the table that the discords vary greatly in length and also that the number of anomalous points extracted in those are not proportional to the length. Thus the integration of exact variable length discords with different anomalous point detectors does bring complementary information and the scoring function is adequate. This information was able to extract the prominent discords over the runoff data that reflected both the driest periods, and periods with highest soil water levels. 3. Climate-related impact insights. Those results are compared to basic statistical analysis of runoff data to demonstrate the thematic insights brought by combining anomalous point detection and discord methods (Figure 6). The first rank of combined PROPHET and discord methods detects an anomalous event/period from August 1983 for 97 months with 12 anomalous points (Table 5, Figure 2). Following discords 2, 3 and 4 exhibit larger time windows but include this first top discord. This is consistent with well known intense droughts that occurred in Sahel between 1982 and 1987 (Figure 6a, 6b), with the most severe drought ever recorded within the African continent in 1983-1984 (e.g. [START_REF] Masih | A review of droughts on the african continent: a geospatial and long-term perspective[END_REF]). Moreover, following ranks of combined PROPHET and discords emphasize top discords starting in 1902 for at least 60 months (Table 5). The first ranked of combined MAD and discord methods also detects "anomalies" in 1904 for 14 months, with following ranked ones including this first top ranked (Table 5). These results coincide with large positive anomalies both in annual maximum and annual mean in 1902 and 1903, and with very large year-to-year variability from 1902 to 1907 (Figure 6a, 6c). These early 1900s result also illustrate the well known soil drying trend resulting from a rainfall deficit observed between 1900 and 2013 within the region ( [START_REF] Ghiggi | Grun: an observation-based global gridded runoff dataset from 1902 to 2014[END_REF]). PROPHET detects both 1980s dry events and 1904 abrupt year-to-year changes, and it is important to note that MAD does not detect any 1980s dry signal in the first top 20 ranked discords (Table 5).

Conclusion and future work

In this paper, we proposed a hybrid approach to rank variable length discords that combine an exact pattern wise outlier approach with anomalous point detectors. Our methodology is generic. We chose the Matrix Profile as the pattern outlier detector because it is exact and parameter free, generated variable length discords, and used point wise anomaly detectors to give each discord a score. This score together with the discord length were used to define a ranking function.

We tested two different point anomaly detectors, a model based detector (MAD) and predictor based one (PROPHET). We run experiments on large scale Monthly runoff data of the Sahel region in West Africa from 1902 till 2005. Our results show that PROPHET is more adequate than MAD, as its top 5 ranked discords did capture the actual discords. The top 5 discords extracted with MAD were only able to reflect the periods with highest water level. They also show that In summary, this work contributed a novel approach to exploit variable length discords and be able to rank them with respect to a complementary data analysis, point anomalous detectors. The ranking function confirmed pattern outliers for the climate impact data at hand. Future work includes running more experimental studies on different impact climate data and other data sets. We also intend to extend our integration with alternative point anomaly detectors such as ARIMA, Random Forest and LSTM, towards a comprehensive approach, given the specificity and impact they have on the extracted point outliers.

Fig. 1 :

 1 Fig. 1: Anomaly points extracted by MAD and PROPHET: similarities and differences

Fig. 2 :

 2 Fig. 2: Anomalous points detected by PROPHET in the subsequence over 1981-1984

Fig. 3 :Fig. 4 :

 34 Fig. 3: Anomalous points from PROPHET in the subsequence 07/1902-01/1905

Fig. 5 :

 5 Fig. 5: Ranking table of top 20 discords

( a )

 a Time series of normalized annual runoff (b) Seasonal cycles (c) Seasonal cycles

Fig. 6 :

 6 Fig. 6: (a) Time series of normalized annual maximum and mean runoff (relative to their respective mean climatologies over 1902-2005 ). (b)-(c) Seasonal cycles of absolute runoff values for each year between 1902 and 2005 (in mm/day).

  represent the top 5 discords resulting from the filtering procedure.

	PROPHET			MAD
	date	window ratio points	date	window ratio points
	1983-08-01	97	0.12 12	1904-06-01	14	0.29	4
	1982-09-01 108 0.12 13	1904-05-01	15	0.27	4
	1982-08-01 121 0.11 13	1902-03-01	43	0.26 11
	1982-07-01 122 0.11 13	1904-04-01	16	0.25	4
	1902-07-01	62	0.1	6	1904-03-01	32	0.25	8
	1902-06-01	63	0.1	6	1904-08-01	24	0.25	6
	1902-05-01	64	0.09	6	1902-02-01	93	0.25 23
	1902-04-01	65	0.09	6	1902-01-01 107 0.23 25
	1902-03-01	66	0.09	6	1903-09-01	35	0.23	8
	1902-02-01

Acknowledgements The authors would like to thank the Occitanie Region, who partially funded this research, and the reviewers for their comments and suggestions.